{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", " \n", "## [mlcourse.ai](mlcourse.ai) – Open Machine Learning Course \n", "###
AEGo, ODS Slack nickname: AEGo\n", " \n", "##
Individual data analysis project" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Research plan**\n", " - Dataset and features description\n", " - Exploratory data analysis\n", " - Visual analysis of the features\n", " - Patterns, insights, pecularities of data\n", " - Data preprocessing\n", " - Feature engineering and description\n", " - Cross-validation, hyperparameter tuning\n", " - Validation and learning curves\n", " - Prediction for hold-out and test samples\n", " - Model evaluation with metrics description\n", " - Conclusions" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import json\n", "from sklearn.model_selection import StratifiedKFold, cross_val_score, GridSearchCV\n", "from sklearn.metrics import mean_absolute_error\n", "from sklearn.metrics import classification_report, confusion_matrix \n", "from sklearn.metrics import accuracy_score, precision_score, balanced_accuracy_score\n", "from sklearn.metrics import label_ranking_average_precision_score, label_ranking_loss, coverage_error\n", "from scipy.sparse import csr_matrix, hstack\n", "from scipy.stats import probplot\n", "import pickle\n", "import matplotlib.pyplot as plt\n", "get_ipython().run_line_magic('matplotlib', 'inline')\n", "import seaborn as sns\n", "import gc\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "import time\n", "import random\n", "import itertools\n", "from scipy.signal import resample\n", "from sklearn.utils import shuffle\n", "from sklearn.preprocessing import OneHotEncoder,StandardScaler\n", "from sklearn.linear_model import LogisticRegression\n", "from sklearn.model_selection import train_test_split\n", "\n", "from keras.models import Model\n", "from keras.layers import Input, Dense, Conv1D, MaxPooling1D, Softmax, Add, Flatten, Activation\n", "from keras import backend as K\n", "from keras.optimizers import Adam\n", "from keras.callbacks import LearningRateScheduler, ModelCheckpoint\n", "from keras.wrappers.scikit_learn import KerasClassifier\n", "\n", "from time import gmtime, strftime\n", "from scipy.signal import butter, lfilter\n", "import pywt\n", "from scipy.signal import medfilt\n", "\n", "import math\n", "import os\n", "\n", "from numpy.random import seed\n", "from tensorflow import set_random_seed\n", "\n", "sns.set(style=\"darkgrid\")\n", "\n", "import tensorflow as tf\n", "\n", "K.clear_session()\n", "\n", "session_conf = tf.ConfigProto(intra_op_parallelism_threads=1,\n", " inter_op_parallelism_threads=1)\n", "\n", "# The below tf.set_random_seed() will make random number generation\n", "# in the TensorFlow backend have a well-defined initial state.\n", "# For further details, see:\n", "# https://www.tensorflow.org/api_docs/python/tf/set_random_seed\n", "\n", "tf.set_random_seed(1234)\n", "\n", "sess = tf.Session(graph=tf.get_default_graph(), config=session_conf)\n", "K.set_session(sess)\n", "\n", "seed(17)\n", "set_random_seed(17)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 1. Dataset and features description" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[ECG Heartbeat Categorization Dataset](https://www.kaggle.com/shayanfazeli/heartbeat)\n", "\n", "This dataset is composed of two collections of heartbeat signals derived from two famous datasets in heartbeat classification, the MIT-BIH Arrhythmia Dataset and The PTB Diagnostic ECG Database. \n", "\n", "The signals correspond to electrocardiogram (ECG) shapes of heartbeats for the normal case and the cases affected by different arrhythmias and myocardial infarction. These signals are preprocessed and segmented, with each segment corresponding to a heartbeat.\n", "Content\n", "Arrhythmia Dataset\n", "\n", "| Physionet's MIT-BIH Arrhythmia Dataset | |\n", "|-|-|\n", "| Number of Samples | 109446 |\n", "| Number of Categories | 5 |\n", "| Sampling Frequency | 125Hz |\n", "| Classes |['N': 0, 'S': 1, 'V': 2, 'F': 3, 'Q': 4] |\n", "\n", "| Class | Description |\n", "|-|-|\n", "| N | Normal |\n", "| S | Supraventricular premature beat |\n", "| V | Premature ventricular contraction |\n", "| F | Fusion of ventricular and normal beat |\n", "| Q | Unclassifiable beat |\n", "\n", " \n", "| Physionet's PTB Diagnostic Database | |\n", "|-|-|\n", "| Number of Samples | 14552 |\n", "| Number of Categories | 2 |\n", "| Sampling Frequency | 125Hz |\n", "\n", "Remark: All the samples are cropped, downsampled and padded with zeroes if necessary to the fixed dimension of **188**.\n", "\n", "Electrocardiogram (**ECG**) is a graphical representation of the electric activity of the heart and has been commonly used for cardiovascular disease diagnosis.\n", "\n", "Data Files\n", "\n", "This dataset consists of a series of CSV files. Each of these CSV files contain a matrix, with each row representing an example in that portion of the dataset. The final element of each row denotes the class to which that example belongs." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "PATH = \"../../data\"\n", "signal_frequency = 125" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "train_df = pd.read_csv(os.path.join(PATH, 'mitbih_train.csv'), header=None)\n", "test_df = pd.read_csv(os.path.join(PATH, 'mitbih_test.csv'), header=None)\n", "full_df = pd.concat([train_df, test_df], axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 2. Exploratory data analysis" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "train_df.shape, test_df.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "full_df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are 109446 samples" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "full_df.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### There are 87554 train samples" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "train_df[187].value_counts(normalize=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "train_df.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### There are 21892 test samples" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "test_df.info()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "test_df[187].value_counts(normalize=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "M = full_df.values\n", "X = M[:, :-1]\n", "y = M[:, -1].astype(int)\n", "\n", "M_train = train_df.values\n", "XX_train = M_train[:, :-1]\n", "yy_train = M_train[:, -1].astype(int)\n", "\n", "M_test = test_df.values\n", "XX_test = M_test[:, :-1]\n", "yy_test = M_test[:, -1].astype(int)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Delete unused variable" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "del M, M_train, M_test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 3. Visual analysis of the features" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_count(data, x, ax, title='Data', xlabel='Class', tick_labesl=['N', 'S', 'V', 'F', 'Q']):\n", " ax1 = sns.countplot(x=x, data=data, ax=ax)\n", " ncount = len(data)\n", " for p in ax1.patches:\n", " x=p.get_bbox().get_points()[:,0]\n", " y=p.get_bbox().get_points()[1,1]\n", " ax1.annotate('{:.1f}%'.format(100.*y/ncount), (x.mean(), y), \n", " ha='center', va='bottom') # set the alignment of the text\n", " ax.set_xticklabels(['N', 'S', 'V', 'F', 'Q']);\n", " ax.set_title(title)\n", " ax.set_xlabel(xlabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Distribution of all data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 5))\n", "plot_count(full_df, 187, ax, title='All Data')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Distribution of train/test data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_count(train_df, 187, ax[0], title='Train')\n", "plot_count(test_df, 187, ax[1], title='Test')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C0 = np.argwhere(y == 0).flatten()\n", "C1 = np.argwhere(y == 1).flatten()\n", "C2 = np.argwhere(y == 2).flatten()\n", "C3 = np.argwhere(y == 3).flatten()\n", "C4 = np.argwhere(y == 4).flatten()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.arange(0, 187)*8/1000\n", "\n", "plt.figure(figsize=(20,8))\n", "plt.plot(x, X[C0, :][0], label=\"Cat. N\")\n", "plt.plot(x, X[C1, :][0], label=\"Cat. S\")\n", "plt.plot(x, X[C2, :][0], label=\"Cat. V\")\n", "plt.plot(x, X[C3, :][0], label=\"Cat. F\")\n", "plt.plot(x, X[C4, :][0], label=\"Cat. Q\")\n", "plt.legend()\n", "plt.title(\"1-beat ECG for every category\", fontsize=20)\n", "plt.ylabel(\"Amplitude\", fontsize=15)\n", "plt.grid(True)\n", "plt.xlabel(\"Time (ms)\", fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 4. Patterns, insights, pecularities of data " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Electrocardiogram (ECG) is a graphical representation of the electric activity of the heart and has been commonly used for cardiovascular disease diagnosis." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Some usefull functions (plot data, filter data, detect R-peaks, HR (heartrate))**:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import signal\n", "from scipy.signal import butter, lfilter\n", "\n", "def butter_bandpass(lowcut, highcut, fs, order=5):\n", " \n", " nyq = 0.5 * fs\n", " low = lowcut / nyq\n", " high = highcut / nyq\n", " b, a = butter(order, [low, high], btype='band')\n", " return b, a\n", "\n", "\n", "def butter_bandpass_filter(data, lowcut, highcut, fs, order=5): \n", " b, a = butter_bandpass(lowcut, highcut, fs, order=order)\n", " y = lfilter(b, a, data)\n", " return y\n", "\n", "\n", "def create_fir_filter(cutoff_hz, fs, n=0, window=\"hamming\"): \n", " N = int(fs / 2)\n", " if n == 0:\n", " n = N\n", " return signal.firwin(n, cutoff=cutoff_hz / N, window=window)\n", "\n", "\n", "def fir_filter(fir, data):\n", " return lfilter(fir, 1.0, data)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# code from https://github.com/c-labpl/qrs_detector\n", " \n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from time import gmtime, strftime\n", "from scipy.signal import butter, lfilter\n", "\n", "\n", "LOG_DIR = \"logs/\"\n", "PLOT_DIR = \"plots/\"\n", "\n", "\n", "class QRSDetectorOffline(object):\n", " \"\"\"\n", " Python Offline ECG QRS Detector based on the Pan-Tomkins algorithm.\n", " \n", " Michał Sznajder (Jagiellonian University) - technical contact (msznajder@gmail.com)\n", " Marta Łukowska (Jagiellonian University)\n", "\n", "\n", " The module is offline Python implementation of QRS complex detection in the ECG signal based\n", " on the Pan-Tomkins algorithm: Pan J, Tompkins W.J., A real-time QRS detection algorithm,\n", " IEEE Transactions on Biomedical Engineering, Vol. BME-32, No. 3, March 1985, pp. 230-236.\n", "\n", " The QRS complex corresponds to the depolarization of the right and left ventricles of the human heart. It is the most visually obvious part of the ECG signal. QRS complex detection is essential for time-domain ECG signal analyses, namely heart rate variability. It makes it possible to compute inter-beat interval (RR interval) values that correspond to the time between two consecutive R peaks. Thus, a QRS complex detector is an ECG-based heart contraction detector.\n", "\n", " Offline version detects QRS complexes in a pre-recorded ECG signal dataset (e.g. stored in .csv format).\n", "\n", " This implementation of a QRS Complex Detector is by no means a certified medical tool and should not be used in health monitoring. It was created and used for experimental purposes in psychophysiology and psychology.\n", "\n", " You can find more information in module documentation:\n", " https://github.com/c-labpl/qrs_detector\n", "\n", " If you use these modules in a research project, please consider citing it:\n", " https://zenodo.org/record/583770\n", "\n", " If you use these modules in any other project, please refer to MIT open-source license.\n", "\n", "\n", " MIT License\n", "\n", " Copyright (c) 2017 Michał Sznajder, Marta Łukowska\n", "\n", " Permission is hereby granted, free of charge, to any person obtaining a copy\n", " of this software and associated documentation files (the \"Software\"), to deal\n", " in the Software without restriction, including without limitation the rights\n", " to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n", " copies of the Software, and to permit persons to whom the Software is\n", " furnished to do so, subject to the following conditions:\n", "\n", " The above copyright notice and this permission notice shall be included in all\n", " copies or substantial portions of the Software.\n", "\n", " THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n", " IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n", " FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n", " AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n", " LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n", " OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE\n", " SOFTWARE.\n", " \"\"\"\n", "\n", " def __init__(self, ecg_data_path, ecg_data_raw=None, signal_frequency=125,\n", " filter_lowcut=0.1, filter_highcut=15.0, filter_order=1,\n", " verbose=True, log_data=False, plot_data=False, show_plot=False,\n", " save_plot=False):\n", " \"\"\"\n", " QRSDetectorOffline class initialisation method.\n", " :param string ecg_data_path: path to the ECG dataset\n", " :param bool verbose: flag for printing the results\n", " :param bool log_data: flag for logging the results\n", " :param bool plot_data: flag for plotting the results to a file\n", " :param bool show_plot: flag for showing generated results plot - will not show anything if plot is not generated\n", " \"\"\"\n", " # Configuration parameters.\n", " self.ecg_data_path = ecg_data_path\n", "\n", " self.signal_frequency = signal_frequency # Set ECG device frequency in samples per second here.\n", "\n", " self.filter_lowcut = filter_lowcut\n", " self.filter_highcut = filter_highcut\n", " self.filter_order = filter_order\n", "\n", " self.integration_window = 8 #15 # Change proportionally when adjusting frequency (in samples).\n", "\n", " self.findpeaks_limit = 0.35\n", " self.findpeaks_spacing = 25 #50 # Change proportionally when adjusting frequency (in samples).\n", "\n", " self.refractory_period = 60 #120 # Change proportionally when adjusting frequency (in samples).\n", " self.qrs_peak_filtering_factor = 0.125\n", " self.noise_peak_filtering_factor = 0.125\n", " self.qrs_noise_diff_weight = 0.25\n", "\n", " # Loaded ECG data.\n", " self.ecg_data_raw = None\n", "\n", " # Measured and calculated values.\n", " self.filtered_ecg_measurements = None\n", " self.differentiated_ecg_measurements = None\n", " self.squared_ecg_measurements = None\n", " self.integrated_ecg_measurements = None\n", " self.detected_peaks_indices = None\n", " self.detected_peaks_values = None\n", "\n", " self.qrs_peak_value = 0.0\n", " self.noise_peak_value = 0.0\n", " self.threshold_value = 0.0\n", "\n", " # Detection results.\n", " self.qrs_peaks_indices = np.array([], dtype=int)\n", " self.noise_peaks_indices = np.array([], dtype=int)\n", "\n", " # Final ECG data and QRS detection results array - samples with detected QRS are marked with 1 value.\n", " self.ecg_data_detected = None\n", "\n", " # Run whole detector flow.\n", " if ecg_data_raw is not None:\n", " self.ecg_data_raw = ecg_data_raw \n", " else: \n", " self.load_ecg_data()\n", " self.detect_peaks()\n", " self.detect_qrs()\n", "\n", " if verbose:\n", " self.print_detection_data()\n", "\n", " if log_data:\n", " self.log_path = \"{:s}QRS_offline_detector_log_{:s}.csv\".format(LOG_DIR,\n", " strftime(\"%Y_%m_%d_%H_%M_%S\", gmtime()))\n", " self.log_detection_data()\n", "\n", " if plot_data:\n", " self.plot_path = \"{:s}QRS_offline_detector_plot_{:s}.png\".format(PLOT_DIR,\n", " strftime(\"%Y_%m_%d_%H_%M_%S\", gmtime()))\n", " self.plot_detection_data(show_plot=show_plot, save_plot=save_plot)\n", "\n", " \"\"\"Loading ECG measurements data methods.\"\"\"\n", "\n", " def load_ecg_data(self):\n", " \"\"\"\n", " Method loading ECG data set from a file.\n", " \"\"\"\n", " self.ecg_data_raw = np.loadtxt(self.ecg_data_path, skiprows=1, delimiter=',')\n", "\n", " \"\"\"ECG measurements data processing methods.\"\"\"\n", "\n", " def detect_peaks(self):\n", " \"\"\"\n", " Method responsible for extracting peaks from loaded ECG measurements data through measurements processing.\n", " \"\"\"\n", " # Extract measurements from loaded ECG data.\n", " ecg_measurements = self.ecg_data_raw#[:, 1]\n", "\n", " # Measurements filtering - 0-15 Hz band pass filter.\n", " self.filtered_ecg_measurements = self.bandpass_filter(ecg_measurements.flatten(), \n", " lowcut=self.filter_lowcut,\n", " highcut=self.filter_highcut, \n", " signal_freq=self.signal_frequency,\n", " filter_order=self.filter_order)\n", "# self.filtered_ecg_measurements[:5] = self.filtered_ecg_measurements[5]\n", "\n", " # Derivative - provides QRS slope information.\n", " self.differentiated_ecg_measurements = np.ediff1d(self.filtered_ecg_measurements)\n", "\n", " # Squaring - intensifies values received in derivative.\n", " self.squared_ecg_measurements = self.differentiated_ecg_measurements ** 2\n", "\n", " # Moving-window integration.\n", " self.integrated_ecg_measurements = np.convolve(\n", " self.squared_ecg_measurements, \n", " np.ones(self.integration_window))\n", "\n", " # Fiducial mark - peak detection on integrated measurements.\n", " self.detected_peaks_indices = self.findpeaks(data=self.integrated_ecg_measurements,\n", " limit=self.findpeaks_limit,\n", " spacing=self.findpeaks_spacing)\n", "\n", " self.detected_peaks_values = self.integrated_ecg_measurements[self.detected_peaks_indices]\n", "\n", " \"\"\"QRS detection methods.\"\"\"\n", "\n", " def detect_qrs(self):\n", " \"\"\"\n", " Method responsible for classifying detected ECG measurements peaks either as noise or as QRS complex (heart beat).\n", " \"\"\"\n", " if self.detected_peaks_indices is None:\n", " return;\n", " for detected_peak_index, detected_peaks_value in zip(self.detected_peaks_indices, self.detected_peaks_values):\n", "\n", " try:\n", " last_qrs_index = self.qrs_peaks_indices[-1]\n", " except IndexError:\n", " last_qrs_index = 0\n", "\n", " # After a valid QRS complex detection, there is a 200 ms refractory period before next one can be detected.\n", " if detected_peak_index - last_qrs_index > self.refractory_period or not self.qrs_peaks_indices.size:\n", " # Peak must be classified either as a noise peak or a QRS peak.\n", " # To be classified as a QRS peak it must exceed dynamically set threshold value.\n", " if detected_peaks_value > self.threshold_value:\n", " self.qrs_peaks_indices = np.append(self.qrs_peaks_indices, detected_peak_index)\n", "\n", " # Adjust QRS peak value used later for setting QRS-noise threshold.\n", " self.qrs_peak_value = self.qrs_peak_filtering_factor * detected_peaks_value + \\\n", " (1 - self.qrs_peak_filtering_factor) * self.qrs_peak_value\n", " else:\n", " self.noise_peaks_indices = np.append(self.noise_peaks_indices, detected_peak_index)\n", "\n", " # Adjust noise peak value used later for setting QRS-noise threshold.\n", " self.noise_peak_value = self.noise_peak_filtering_factor * detected_peaks_value + \\\n", " (1 - self.noise_peak_filtering_factor) * self.noise_peak_value\n", "\n", " # Adjust QRS-noise threshold value based on previously detected QRS or noise peaks value.\n", " self.threshold_value = self.noise_peak_value + \\\n", " self.qrs_noise_diff_weight * (self.qrs_peak_value - self.noise_peak_value)\n", "\n", " # Create array containing both input ECG measurements data and QRS detection indication column.\n", " # We mark QRS detection with '1' flag in 'qrs_detected' log column ('0' otherwise).\n", " measurement_qrs_detection_flag = np.zeros([len(self.ecg_data_raw), 1])\n", " measurement_qrs_detection_flag[self.qrs_peaks_indices] = 1\n", " self.ecg_data_detected = np.append(self.ecg_data_raw, measurement_qrs_detection_flag, 1)\n", "\n", " \"\"\"Results reporting methods.\"\"\"\n", "\n", " def print_detection_data(self):\n", " \"\"\"\n", " Method responsible for printing the results.\n", " \"\"\"\n", " print(\"qrs peaks indices\")\n", " print(self.qrs_peaks_indices)\n", " print(\"noise peaks indices\")\n", " print(self.noise_peaks_indices)\n", "\n", " def log_detection_data(self):\n", " \"\"\"\n", " Method responsible for logging measured ECG and detection results to a file.\n", " \"\"\"\n", " with open(self.log_path, \"wb\") as fin:\n", " fin.write(b\"timestamp,ecg_measurement,qrs_detected\\n\")\n", " np.savetxt(fin, self.ecg_data_detected, delimiter=\",\")\n", "\n", " def plot_detection_data(self, show_plot=False, save_plot=False):\n", " \"\"\"\n", " Method responsible for plotting detection results.\n", " :param bool show_plot: flag for plotting the results and showing plot\n", " \"\"\"\n", " def plot_data(axis, data, title='', fontsize=10):\n", " axis.set_title(title, fontsize=fontsize)\n", " axis.grid(which='both', axis='both', linestyle='--')\n", " axis.plot(data, color=\"salmon\", zorder=1)\n", "\n", " def plot_points(axis, values, indices):\n", " axis.scatter(x=indices, y=values[indices], c=\"black\", s=50, zorder=2)\n", "\n", " plt.close('all')\n", " fig, axarr = plt.subplots(6, sharex=True, figsize=(15, 18))\n", "\n", " plot_data(axis=axarr[0], data=self.ecg_data_raw, title='Raw ECG measurements')\n", " plot_data(axis=axarr[1], data=self.filtered_ecg_measurements, title='Filtered ECG measurements')\n", " plot_data(axis=axarr[2], data=self.differentiated_ecg_measurements, title='Differentiated ECG measurements')\n", " plot_data(axis=axarr[3], data=self.squared_ecg_measurements, title='Squared ECG measurements')\n", " plot_data(axis=axarr[4], data=self.integrated_ecg_measurements, title='Integrated ECG measurements with QRS peaks marked (black)')\n", " plot_points(axis=axarr[4], values=self.integrated_ecg_measurements, indices=self.qrs_peaks_indices)\n", " plot_data(axis=axarr[5], data=self.ecg_data_detected[:, 1], title='Raw ECG measurements with QRS peaks marked (black)')\n", " plot_points(axis=axarr[5], values=self.ecg_data_detected[:, 1], indices=self.qrs_peaks_indices)\n", "\n", " plt.tight_layout()\n", " if save_plot:\n", " fig.savefig(self.plot_path)\n", "\n", " if show_plot:\n", " plt.show()\n", "\n", " plt.close()\n", "\n", " \"\"\"Tools methods.\"\"\"\n", "\n", " def bandpass_filter(self, data, lowcut, highcut, signal_freq, filter_order):\n", " \"\"\"\n", " Method responsible for creating and applying Butterworth filter.\n", " :param deque data: raw data\n", " :param float lowcut: filter lowcut frequency value\n", " :param float highcut: filter highcut frequency value\n", " :param int signal_freq: signal frequency in samples per second (Hz)\n", " :param int filter_order: filter order\n", " :return array: filtered data\n", " \"\"\"\n", " nyquist_freq = 0.5 * signal_freq\n", " low = lowcut / nyquist_freq\n", " high = highcut / nyquist_freq\n", " b, a = butter(filter_order, [low, high], btype=\"band\")\n", " y = lfilter(b, a, data)\n", " return y\n", "\n", " def findpeaks(self, data, spacing=1, limit=None):\n", " \"\"\"\n", " Janko Slavic peak detection algorithm and implementation.\n", " https://github.com/jankoslavic/py-tools/tree/master/findpeaks\n", " Finds peaks in `data` which are of `spacing` width and >=`limit`.\n", " :param ndarray data: data\n", " :param float spacing: minimum spacing to the next peak (should be 1 or more)\n", " :param float limit: peaks should have value greater or equal\n", " :return array: detected peaks indexes array\n", " \"\"\"\n", " len = data.size\n", " x = np.zeros(len + 2 * spacing)\n", " x[:spacing] = data[0] - 1.e-6\n", " x[-spacing:] = data[-1] - 1.e-6\n", " x[spacing:spacing + len] = data\n", " peak_candidate = np.zeros(len)\n", " peak_candidate[:] = True\n", " for s in range(spacing):\n", " start = spacing - s - 1\n", " h_b = x[start: start + len] # before\n", " start = spacing\n", " h_c = x[start: start + len] # central\n", " start = spacing + s + 1\n", " h_a = x[start: start + len] # after\n", " peak_candidate = np.logical_and(peak_candidate, np.logical_and(h_c > h_b, h_c > h_a))\n", "\n", " ind = np.argwhere(peak_candidate)\n", " ind = ind.reshape(ind.size)\n", " if limit is not None:\n", " ind = ind[data[ind] > limit]\n", " return ind" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class FilterHelper():\n", " def __init__(self, \n", " signal_frequency=125,\n", " filter_lowcut=0.01, \n", " filter_highcut=50.0, \n", " filter_order=1,\n", " cutoff_hz=15.0,\n", " verbose=False, \n", " log_data=False, \n", " plot_data=False, \n", " show_plot=False,\n", " save_plot=False):\n", " self.signal_frequency = signal_frequency\n", " self.filter_lowcut = filter_lowcut\n", " self.filter_highcut = filter_highcut\n", " self.filter_order = filter_order\n", " self.cutoff_hz=cutoff_hz\n", " self.verbose=verbose \n", " self.log_data=log_data\n", " self.plot_data=plot_data\n", " self.show_plot=show_plot\n", " self.save_plot=save_plot\n", " \n", " def rr_peaks_detect(self, x):\n", " qrs_detector = QRSDetectorOffline(ecg_data_raw=x.reshape(-1,1), \n", " filter_lowcut=self.filter_lowcut,\n", " filter_highcut=self.filter_highcut,\n", " filter_order=self.filter_order,\n", " signal_frequency=self.signal_frequency,\n", " ecg_data_path=None, \n", " verbose=self.verbose,\n", " log_data=self.log_data, \n", " plot_data=self.plot_data, \n", " show_plot=self.show_plot)\n", " return qrs_detector\n", "\n", " def qrs_filter(self, x):\n", " qrs_detector = QRSDetectorOffline(ecg_data_raw=x.reshape(-1,1), \n", " filter_lowcut=self.filter_lowcut,\n", " filter_highcut=self.filter_highcut,\n", " filter_order=self.filter_order,\n", " signal_frequency=self.signal_frequency,\n", " ecg_data_path=None, \n", " verbose=self.verbose,\n", " log_data=self.log_data, \n", " plot_data=self.plot_data, \n", " show_plot=self.show_plot)\n", " return qrs_detector.filtered_ecg_measurements.flatten()\n", "\n", " def fir_filter(self, x):\n", " # The cutoff frequency of the filter.\n", " flt = create_fir_filter(self.cutoff_hz, self.signal_frequency)\n", " return fir_filter(flt, x)\n", "\n", " def bandpass_filter(self, x):\n", " return butter_bandpass_filter(\n", " x, \n", " self.filter_lowcut, \n", " self.filter_highcut, \n", " self.signal_frequency, \n", " order=self.filter_order)\n", " \n", " def compose_filters(self, x):\n", " return self.fir_filter(self.bandpass_filter(x))\n", " \n", "def calculate_heart_rate(qrs_detector, signal_frequency=125):\n", " sec_per_minute = 60\n", " hr = 0\n", " i1 = qrs_detector.qrs_peaks_indices[0] if len(qrs_detector.qrs_peaks_indices) > 0 else 0\n", " i2 = qrs_detector.qrs_peaks_indices[1] if len(qrs_detector.qrs_peaks_indices) > 1 else 0\n", " if len(qrs_detector.qrs_peaks_indices) > 1:\n", " hr = signal_frequency / (i2 - i1) * sec_per_minute\n", " return hr, i1, i2\n", "\n", "def plot_ecg(data, filter_helper=None, show_filtered=True, show_hr=True, title='ECG', ax=None):\n", " if ax is not None:\n", " plt.axes(ax)\n", " else:\n", " plt.figure(figsize=(20,6))\n", " \n", " if show_filtered:\n", " flt = filter_function(data)\n", " \n", " plt.plot(data, color='b', label='Original ECG')\n", " \n", " if show_hr:\n", " qrsdetector = filter_helper.rr_peaks_detect(data)\n", " hr, i1, i2 = calculate_heart_rate(qrsdetector)\n", " print('Heart Rate: {}, R1: {}, R2: {}'.format(hr, i1, i2))\n", "\n", " plt.scatter(x=qrsdetector.qrs_peaks_indices, \n", " y=qrsdetector.ecg_data_detected[:, 1][qrsdetector.qrs_peaks_indices], \n", " c=\"black\", \n", " s=50, \n", " zorder=2)\n", "\n", " if show_filtered:\n", " plt.plot(flt, color='r', label='Filtered')\n", " plt.title(title, fontsize=16)\n", " plt.grid(True)\n", " plt.legend()\n", " if ax is None:\n", " plt.show();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Examples of ECG for each class" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_ecg(X[np.argwhere(y==0)[0][0], :], ax=axes[0],\n", " title='Normal ECG', show_filtered=False, show_hr=False)\n", "plot_ecg(X[np.argwhere(y==1)[0][0], :], ax=axes[1],\n", " title='Class: Supraventricular premature beat', \n", " show_filtered=False, show_hr=False)\n", "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_ecg(X[np.argwhere(y==2)[0][0], :], ax=axes[0],\n", " title='Class: Premature ventricular contraction', \n", " show_filtered=False, show_hr=False)\n", "plot_ecg(X[np.argwhere(y==3)[0][0], :], ax=axes[1],\n", " title='Class: Fusion of ventricular and normal beat', \n", " show_filtered=False, show_hr=False)\n", "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_ecg(X[np.argwhere(y==4)[0][0], :], ax=axes[0],\n", " title='Class: Unclassifiable beat', \n", " show_filtered=False, show_hr=False)\n", "fig.delaxes(axes[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 5. Data preprocessing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The baseline of the signal is substracted. Additionally, some noise removal can be done.\n", "\n", "Two median filters are applied for this purpose.\n", "Some ideas from [this site](https://github.com/mondejar/ecg-classification)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def baseline_ecg(ecg):\n", " # Remove Baseline\n", " baseline = medfilt(ecg, 71) \n", " baseline = medfilt(baseline, 215) \n", " return ecg - baseline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Example of ECG with substracted baseline:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "orig_ecg = X[C4, :][0]\n", "\n", "flt_bl = baseline_ecg(orig_ecg)\n", " \n", "plt.figure(figsize=(20,8))\n", "plt.plot(x, orig_ecg, label=\"Original\")\n", "plt.plot(x, flt_bl, label=\"Baselined\")\n", "plt.legend()\n", "plt.title(\"'Baselined' 1-beat ECG\", fontsize=20)\n", "plt.ylabel(\"Amplitude\", fontsize=15)\n", "plt.grid(True)\n", "plt.xlabel(\"Time (ms)\", fontsize=15)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "XX_train_bl = np.apply_along_axis(baseline_ecg, 1, XX_train)\n", "XX_test_bl = np.apply_along_axis(baseline_ecg, 1, XX_test)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "XX_train_bl.shape, XX_train.shape, XX_test.shape, XX_test_bl.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Data augmentation**\n", "\n", "To train properly the model, we sould have to augment all data to the same level. Nevertheless, for a first try, we will just augment the smallest class to the same level as class 1.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def stretch(x):\n", " l = int(187 * (1 + (random.random()-0.5)/3))\n", " y = resample(x, l)\n", " if l < 187:\n", " y_ = np.zeros(shape=(187, ))\n", " y_[:l] = y\n", " else:\n", " y_ = y[:187]\n", " return y_\n", "\n", "def amplify(x):\n", " alpha = (random.random()-0.5)\n", " factor = -alpha*x + (1+alpha)\n", " return x*factor\n", "\n", "def augment(x):\n", " result = np.zeros(shape= (4, 187))\n", " for i in range(3):\n", " if random.random() < 0.33:\n", " new_y = stretch(x)\n", " elif random.random() < 0.66:\n", " new_y = amplify(x)\n", " else:\n", " new_y = stretch(x)\n", " new_y = amplify(new_y)\n", " result[i, :] = new_y\n", " return result" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,10))\n", "plt.plot(X[0, :], label='original')\n", "plt.plot(amplify(X[0, :]), label='amplify')\n", "plt.plot(stretch(X[0, :]), label='stretch')\n", "plt.title(\"'Amplified/Stretched' ECG\", fontsize=20)\n", "plt.ylabel(\"Amplitude\", fontsize=15)\n", "plt.xlabel(\"Time (ms)\", fontsize=15)\n", "plt.legend()\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "result = np.apply_along_axis(augment, axis=1, arr=X[C3]).reshape(-1, 187)\n", "classes = np.ones(shape=(result.shape[0],), dtype=int)*3\n", "X = np.vstack([X, result])\n", "y = np.hstack([y, classes])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "subC0 = np.random.choice(C0, 800)\n", "subC1 = np.random.choice(C1, 800)\n", "subC2 = np.random.choice(C2, 800)\n", "subC3 = np.random.choice(C3, 800)\n", "subC4 = np.random.choice(C4, 800)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "use800_test = False" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if use800_test:\n", " X_test = np.vstack([X[subC0], X[subC1], X[subC2], X[subC3], X[subC4]])\n", " y_test = np.hstack([y[subC0], y[subC1], y[subC2], y[subC3], y[subC4]])\n", "\n", " X_train = np.delete(X, [subC0, subC1, subC2, subC3, subC4], axis=0)\n", " y_train = np.delete(y, [subC0, subC1, subC2, subC3, subC4], axis=0)\n", "else:\n", " X_test_orig = XX_test\n", " X_test = XX_test_bl #XX_test\n", " y_test = yy_test\n", "\n", " X_train_orig = XX_train\n", " X_train = XX_train_bl #XX_train\n", " y_train = yy_train\n", "\n", "X_train, y_train, X_train_orig = shuffle(X_train, y_train, X_train_orig, random_state=17)\n", "X_test, y_test, X_test_orig = shuffle(X_test, y_test, X_test_orig, random_state=17)\n", "\n", "\n", "del X\n", "del y\n", "del XX_train, XX_test, XX_train_bl, XX_test_bl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### X_train_orig, X_test_orig - original(raw) data\n", "\n", "#### X_train, X_test - 'baselined' data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ptrain_df = pd.DataFrame(y_train, columns=['class'])\n", "ptest_df = pd.DataFrame(y_test, columns=['class'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Class counts of train/test data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(ptrain_df['class'].value_counts())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(ptest_df['class'].value_counts())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_count(ptrain_df, 'class', ax[0], title='Train')\n", "plot_count(ptest_df, 'class', ax[1], title='Test')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# delete unneccessary variables\n", "del ptrain_df\n", "del ptest_df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 6. Feature engineering and description " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Filter ECG signal using Butterworth filter**\n", "\n", "The high-pass and low-pass filters together are known as a bandpass filter, literally allowing only a certain frequency band to pass through. We will use bandpass filter with cutoff=50." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "filter_helper = FilterHelper(\n", " signal_frequency=125,\n", " filter_lowcut=0.01, \n", " filter_highcut=50.0, \n", " filter_order=1,\n", " cutoff_hz=30.,\n", " verbose=False, \n", " log_data=False, \n", " plot_data=False, \n", " show_plot=False,\n", " save_plot=False)\n", "\n", "filter_function = filter_helper.bandpass_filter\n", "\n", "\n", "def calc_hr(data):\n", " qrsdetector = filter_helper.rr_peaks_detect(data)\n", " hr, i1, i2 = calculate_heart_rate(qrsdetector)\n", " return hr, i1, i2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A typical ECG-based heartbeat mainly consists of three waves including P-wave, QRS complex ([wiki](https://en.wikipedia.org/wiki/QRS_complex)), and T-wave. The QRS complex is the most prominent feature and it can be used to obtain additional useful clinical information from ECG signals, such as RR interval, QT interval, and PR interval, etc. Thus, QRS detection is critical for ECG-based health evaluation.\n", "\n", "The QRS complex is the most noticeable feature in the electrocardiogram (ECG) signal, therefore, its detection is critical for ECG signal analysis.\n", "\n", "See for more detail info - [Electrocardiography](https://en.wikipedia.org/wiki/Electrocardiography)\n", "\n", "We use **QRSDetectorOffline** for detecting R-peaks (QRS)." ] }, { "attachments": { "%D0%B8%D0%B7%D0%BE%D0%B1%D1%80%D0%B0%D0%B6%D0%B5%D0%BD%D0%B8%D0%B5.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAFJCAIAAABb9M4+AAAgAElEQVR4nO3dd1gTWdsH4DOTEENTmjSxIGLB3ntDxYKKFYRFURRpIgiIBQWkCEpvUi0g9oa9d1RAYNd3d3XfLW5Tt7j2XQto+P442fnyJhBCgMxJ8jzX79qLzJwJMzl472TKGYSUtkyNjBbNmJEfEXGrsPDJ5ctv7typqapqrlRXVv557dpXR48eSUpas2TJQBsbiqKaZbV1tPgO4/snr3a+mLv65/MJz0syqqvymzFPr6d9d2rLmW2BMX5zJwztoaHBbZbV5nBp2wndo2Nmnz3n/8PD2Ocv0gS1+fXl/Yfsp3+lVFRuLNqzzNNrbCdLo2ZZB4SQQTv9ce4jvXctiS5bl/NnYuG7zH2CXLHs+js967eEiJtr3Ld9NsxxkJaeVnP99vZd+0x1DfTbenjzwS+zrj7bVV5TWCkgIdnXXiSe+CE4/ew8n2ibweNpDqfpG0vTnCE9Jq2cn5AVfPVMwm83Mt/eyfl4MeXZ7o1VEe6F9sPdWmsbNv23QDVQdiNGnExPf3/3bjPq1mAeFBevXrxYV1tb7tXubW2xPXLpi1uZzaub9Dy5nJwQ6NTe1EDu1TYz19uydd5vvydK0U16PgnybpasdVowhObQ8q0DRVEDpvcJvbhqz8dsSd2kp/Bdpv+B5V2Gdpb7E9BoxZ/o6Lv54JesoyZj0i/85uS/Va+tmXzbq8nXdZ+24dTWR+V5tVJSsu195LI9lmY2cn+wUNJqQI8eJQUFijROLH9eu+a7YAFNN+4frbmx3p645R8q8xTJnGj+LsuOD3TU1eY3arW1tFvFxs199z5LbubEcv9B1MRJPRq1DgihLkM7b64Ibaxxkgk5tdLEyrhRv5qiqFHTF6Wdf8K6X3Jk+513jn6xrTQbt2M7eYjL+aQ/pTMnmjs5H4OdMzRb6TTqt0BJK5qmw7y8FLw3V19u7NzZ3tRUxjV3tBvy1410tpgTzfentwzrYyXjag8c1PH7H2KbiznRZOcs5GtqyNTpHNo5bo4ce3P1peBt5gTPMTJ+AroGbYPTz7JuVhOTcPz7TjYDZdlenobmpqVFsjMnmiMx31ma95Txg4WSVlqamifT01k3TjS/Xb48pFcv6atNUVSM31zWjRPNP2XZC6ePaPADd3Qa9P5DdktIh1N+d0NbY13p68DX5YdeWtVczIlmWY5rg/vm5p17JJ/+mXWqmiXb77wbZuckfXu1+K3z19ySTzqcK2mv+nYZ1eCfFpS00uTzr27fzrpuknl569bIfv2krHnaWhfWdZPMh8o8z/njpay22+KRHz/ltpx0OA++iZbiHV+HH122riWkw1m530OKdxZdemVefso6Us2YXXc/jp65uL7t1dBolRtysynS4VxLf9O9o0x7kVB1FE3Tx1NTWXetvvx57Vp3S8s61zzUYwbrrtWX9xW5s2z717naU6f1qq7JaWnpcO5WbKjz+yzNoUMvtsg+nWjcUhfU+QnoG7dLPfeIdZ6aPTvLqvuMnFrnJoe7FzRdOpzT8Y8NWst6hAfqf2r9smWsiyY9nx88qMkXP/A/cZjN+4pc1lGTkr9upFu1Fz9a376DwV/PUhQjHU5u3iLJTneKmdXS0uEMXzBY7FfTHE5o/g3WYWqhbLvyl6Fpe7FNthvs3FzS4SSuONm0f/RqWTadO78tK2Odswaz2d9fdLV1tPgPz25lnbMGc3X7GrGLB0+d9lekdDjT7HuLroPlgI57aprtjIT0bH+RomfaRvS3T1sUzDpJLZrVGedEt1eTr3sm4bfmxa48r3b8gLnNjYGq16XcXNYhkyVvy8qsO3RgVjvabw7rkMkYp8lDmNWe6dBP8dIJavP/+20Mh/v/h88i76xVjHQ4y0V2LVsbGOfefM26Ry2d/mNnMpvsOSuq2aUrz6s9HP1tYy/PUusa1qcP64rJnvyICLza+q21n91s5psiWi73DkdStHDnrvzuBlawE9TmL1w0HK9DH7ueipRunyB394dtbTsJb/BwWhnHukQKSNTeKrxH34qndSH5r5bADnbuGlf7tmxhnTDZ809ZmaGeHkLI33US64Q1KhOH2SCEho+wYks6QW1++d0NuNNXn1yhYOz2CXKd4+YghLg8XtbVZ6xLpJh07TcSITR1qGsLSVeeV5vqf67+f9xQItVaR6d573VVQLzmz0cIle0NY92vRmVH5FKE0LYsVxaxE9Tmd+tu2tpIt6g6S/HYZfyyhaKowbZzWDdIYVm8PgshlOB7vOWwu51Vrast/02KalQzx41jHa/G5mhysrFBaxbvCZMvjy4lURT18Mc4drHzD5g4wnmI4qXDad+r3dKNeawbpLAknfqJpukrqS9bDrvyvNqx/WaxDYky1JZVq1jHq7F5ev36nAkDWcdLjowZ1pVd6QS1+UeOertnurCFnZ3PuC1HHrBukCIzsO+kFpWuPK92xdytbEOiDEXazWEyJi7AiXW55Ej4agfWsbv/IGrj1SC2sHPPXLizrJp1gBSZxc4tch5WNFt9jrENiTLUF4cOsS6XHDmavIp1ueTIvrzlrGP3z9ttqT9sZgu79WfCWddHwQkLPNDS2BVuqGQbEmWoh2fOsC6XHLmcF8q6XHLk9EEWriWWTO7TJLawi76ZzLo+Ck5c6NmWxu5w9LdsQ6IM9ejiRdblkiMlBRGsyyVHLhUHsy6doDZ/x+tUtrDbUpbJuj4KTkLElZbG7kTcz2xDogwF2CkWu9WsSwfYKTiJgB0h1bzYSb6/Dp8/edjw89nZ5GBX10pyJw/rfj47UHozhJCeNm/GmN5X8kOaC7vff09aHTLJ2tqAw6E4HKprN8N166f8+WeyWLM6V0Zfn+fg0Pv6jRDRlvf+E+7hMbJHj7a6uhoIIQ6HsrY28PEZ8/0Pm2XELu5e2DiPESY9jHi6GgghikMZWRtM8BmT8n0M00b6HxUh2ElfSdKww2sl+3TArtHVEtgxL1/dvn1j584enTohhHbHbCYKO+bly9vbru9c26OTMUJod8yy+ppVV+U/u5lxLivQ0kwPIbQ3dnnTsbt9Z52WFqd9B93jx31fvsp49Trj2DHvdu10dHS4dys2SGInOuXlq4yLFwM7d9ZDCB04IDzvcfiwF0LIyIh/8ODyp3+lvv+Q/fMvWzIyFmBJy8pDG8Qu4LAnQkjTiO930CP3r+TC99vSf45zy3CiOBTFoaLK10sugteN8D07vJKs/GrAjpRqUexwyvfsQQh1NDUlEzucsj0bEEIdTdtIb1ZdlX97dyhCqJOpXhOxe/Y8TV+fp63Nefw4QVSxXx/Ft+JzjIz4L16mS8EOp6w8FCFkaamHX1pY6CCELl8OEmu2e7c7QmjEiE4NYtfaQgchFHpZfJA7791LEEJWIzoBdoCdspYCsHt1+zZCiJKYThR2L29vwyvZIHYvb2UihGiaaiJ2sbGzEEKbImdIErYxbBpCKCFhboPYvfk7EyHE4VD4JZdLIYT+/mebWLMP1Tk3S9Y8fZrSIHY0l0II7fonQ2z67g9Z4SUhOXWdxsXrBtgBdqSXArArLSpCCHXv2JFk7O4UhSKEunds2yB2l3KDEUKDbdo3EbshQ9ojhCqrNta3v8bsiEnB7urVYITQkCHt8cvpM3ohhJKS5st9gqLvjJ4Ioc+S5smIF2AH2ClNtSh2b0pLbxUU9LS0RAgdSUoiE7vXpVklBet6WpoghI4k+kjB7nlJxtltqyyMdGiaku8chSh2fD4HIfTqdYYkRi9epiOEtLQ4UrB79TrjwoVVFhY6HA7FnKP444/k8eO7IISsrQ1C1tgVF/uIfUduELvsPxK7jbdCCBlZG9iHTAos9s58vBWwA+xUoVr6bCyPy7UdNOgC2WdjeVzadpD1hewg6c0QQjSF3GYM+8/hSPl+tSh2+A1rPtbxwJ33H7LxXLHG4itDo8VLhn19P1Js8bsVG0I3TB09xhJ/qzU31/b3t334Y6ws2OHEVIQ6hE6xHmOJv9XqmmtP9h+f+rDu+y7wygB2gB3ppYCvsS2RZv8aK0uz2eP7IoQ+Pyj/r5bcs3v+Ik0Su19+jUcN7dnNndsXIfSfLyOkfF39+Cnvi3vh0TEOeno8Doe6dm21jNgx2fMxJ+5e2PyYmTw9HsWhNl6r46ZavG6AHWBHegF2sje7uy8MITRjTO9mwW748I4Iodt31kkidev2WoTQyJHSjtlVfR6GEHJw6C398BzOvf+EI4R69zZpLHZM4u6FIYTMepsAdi2EHYUohFBpzkfRibezqxFCNMUB7JqhALtGNZs0tBtCqGzPhqZjt3XrHITQ6pBJkjbhs7EpKY5SsBPU5tvZdUMIVVRuYKbk5y/y9h799K9UsZZv32UhhGgaNYjdsvyFtt6jc/9KFpte8DYTIUTRdaCG1w2wayJ2pnqdEUKntv4iOvFIzHcIISvTPoBdMxRg16hmV7eHIITshnZvOnbPX6Thb5fffR8jCtOjxwlaWhwTE83XbzKlY3fjZghCaPLk7swUF5fBCKHY2FliLY8c8UII9etv1iB2w10GIYQcYx3Epgcc8UQIWfQzA+xaCLtl9uEIIZ9ZsZITAx1TAbtmKMCusc2G9eqIELqxa20TscNfV/H1w4cPe+E7KE6eXNG+g66+Pu+Le+GiLfHKSO4D4u/Ct26vxS9/fRRvaamHEAqPsP/l1/jqmtwnTxLz8xfx+RwOhyq5tbZB7DJ+3WJgqYcQmh1hn/HLlqLq7G1P4pflL+TwORSHiri1BrBrIexuZv7T1XwgTXFCF+ZfTnl2NuHJijlbEEKDuk68nV0N2DVDAXaNbXYqwx8hNLq/VdOxE9TmP36c4O9v26Gj8JmqXboYrFs/hbn6t0Hszp71RwiNHfv/D/H561nqxrBpffua8Xg0QojLpXr0aBsQ0IizsXnPkmeFTbXoa0bzaIQQzaVMehhNDoCzsS2LXXlebcm2d14zo9sbdqcQRVMca7P+QU5pd3JqGlwQsJOpYNQTRUbKqCfZ2a4IoXPnAmQ54dDEwKgnZGIndwA7mQqwIwS7H3+KQwj5+Y0D7FQsgB0pBdgRgp2gNn/BggEcDnX9+mrJm1sBO+UNYEdKyYedjaXli5KSmqqqp9ev0xT17ObNmqqqFyUlXdu3Jwq7azvWaPI4WjxObpgb69JJYvf8RdrMmb0oCvXo0bayauPLVxn203tSFNLR4TYKL1TX4TxZsAsvCeFqcbjaHI/ti5QaOxPLnrk33xRWCrKuvaAoOufGq8JKQe7NN0YWXUnDrpNxz+sZb8rzaq+kvqAQfS39VXle7fWMNxaG3QC7li35sFu7ZAm+13VvbNywXr0OJiTUVFUdTkz0dXIiCjtjPc37xdFlezbMn9i/hfwSHTugsdh5e4+OjJr5/EXa3YoNCxYMkHtPDcmLnZYJP+m/UTGVoYMd+7eQbiY9jBSAnf2Sdf6JxYWVAp/YA5a9hq+IP1JYKfBLOGbr6Ecadm5T1sf7FJfn1cYsP9Cr4/AtXkfK82q3eh+bP9YPsGvZkg+70qIidweHmqqqhfb2J9LS3GbMqKmqWjJz5sn0dI/Zs7kcjgaHsysy6s9r14bY2DBLDe3Z889r106kpRnq6hrr6V3MyWlp7Ez0tbJCF769m4tfnkhbaajLN9bTvJgTVF2VfzDeq5UGZ7WbHU2h9xV5XA69zn2qib7W4QTvlS62PC69LdRVbKn3FXmtNDjr3KfSNJW9YeEq14kIIVMDbfmwMzHR/OXXeNEpxcU+enq81q01jhzx+vgpT0ODDt0w1dRU69gx74AAWx6PzslxrfmYS9MoMHACn8+5WbKGwe70aT9DQ76JieaVK0GC2vxnz9O6dzd89z773fvsHj3aPnueJomdtqnW0lzXoups/HL16RWahnwtE37olVX7BLkBRzw5rTj2ayZRNNrzMYfWoB1Cp2ibaq065jU5YDzNo91zXMSW2vMxh9OK4xA6heJQS/NcpwZNQAjpmmm3NHYRu++OdlhWWCkYbu8WmHp6xIwlhZWCkTOXBqWdGTvbk+Zo0Bze8sjCrKvPO9oMZZbq1HOYZ1QRX9dQS89kTfYVxWBXEHp35ohl5Xm104a6Jfudnj5sSXle7YzhS1NWnpk1ypOmNTg0L8K98HLq8x4dhjJL2XQYtmlpkS7fsI2mybZA8e/FgJ1MJR921ZWVA7p1+1BR0b9r15qqqv5du36oqLAwMnp95w5u8M3x4zwut6aqys/Z+ey2bTVVVeezs/2cnWuqqizNzG4VFNwuLMTLtih2ZXs2jBvYxVCXX7TZo7oq39JMr6Rg3a3C9f27mldX5Rvq8q/tWHMs2RchVF2VjxA6kbbyyyORNE0VJ/t9eSSSppDkUgih46l+Xx6J1OFz8Uu59+woCr17ny06pUPHNteurb5+fXWHjm2wYqdP+91/EMnhUCdPrrj/IBLfBYEQuno1+MbNkLZtNRnsOnfWu1O6rrRsff8B5vjdAgMnbN06Z8uWOWL3aTDYxVSGdhtvpWnIX7Fv2T5BrkFnvU2layPL1rUfYL5PkKtpyA8vCQk64YMQ2ifIRQitPr0i/v4mikMFn/SNv78J31AhthRCKPiUb/z9TRo6XPxSAXt2BRWfLLoNLLj70cJ6QGGlwMJ6QMHdj7pGFvm3/8EN4ou/o7m8wkrBRGf/1ZkXCisFIVmXJjr765t1DttVGl5QhhdUAHZluYKu5gNLcz9Zmw8oz6u1Nh9QmvvJSMeiZNtb3OBozPdcmleeV+tk658ecKE8rzZz1SUnW39Tvc471pXuXF+GFwTsGl1yn6DwcXTcERGxaPr0mqqqz6ZO3RUZNWHw4Jqqqj2xcZo8Hn7zmqqqh2fO2A0dWlNVNW3UqB9On66pqqIoCs+lKaqlscO5tmONLl+juir/39+MsGIIobd3c1+XZqF/sXtfkSf2Q51Lic5FTcDOxETzh4exYvy9/5D9/kM2RQlR+/gpT+wH/N8P1TkfqnOYl3hZ4Xr+e1vY7TvrOnRsY2GhI3YHrtgJivCSEJ6uxj5BLvr3HbBiCKGi6mx8oxh+uedjjtgPdS71P3MVgl1hpcB2/oplEbtGTF9cWCkYNtV1eWRh98ETCysF3rH7uTwtvIaFlYLk0z/bDJ1cWCnoPco+6dRP6N9PjaJoxWBXnlc7b+yKsMW77IctLs+rnTx4YYR74aCuE8vzamM89vM4wlUtz6s9ueXnId0nl+fVjrSZfiLuJ+rfD5pCNGAnT8mN3cWcHHNDw4Ko6Jqqqh0RER1MTJKCgmqqqoz19EqLisqKirgczuNLl2qqqpynTDmYkLBg8mS8YGdz87t798rNXKOwWzZ75H8OR35xaBPGrrO5fvnejczc1loatwrXn0z3R//LltgPYkuJzUUI/XktTT7s/PzGrQ+d8ux52oNvopycBjB7dlevBnfq1Eb0YJzYDwih6zdCbpasMTPTZiZaWemLDQXq7T06PX1BcvJ8H58xdWI3bvmI+K83bf0yAmNnaKW/uWoDAxOvtUZk2bqQM37of9kS+0FsKbG5CKH85ykKwG5N9hVdw3aeUUWFlYJlEbv0TDq6BKUUVgq09Ewidt/dtLuC4nDTL/5eWCkYOtllRfyRIXbOhZUCA3OryD1VLbRK9WG3LfCKoU67TUuLyvNqwxbvMm7TcZVjSnlebRtNk4LQuwUbKmiKez7x9/K8WrtBLlu8jkwa6FyeV2uub7V7YxUcs5O/5Mbu3d27WjzeL+fP11RV/Xj2LELoq6NHa6qqdkVGcWg6OTg4xM2NpumaqqrK/ft5XC4D3JnMTGM9PZqmV7m6tjR2uWFuOnyuFo+zY9OS6qr8M5kBxnqaNE2tcp1YXZW/PWKxBodev3Qqkoqd2FJic+P85/K4tHzYvXqdMWtWH5pGNjbG5XdD8TG71q019PV5x4/7Ssdu5crxfD4H77LhiefPB5iYaHI4VFDwREFt/l/PUm1sjD9U57x7n92li0Gdx+w8ti/S0OFytTleBYv3CXLXnl+pZcKnONTUoAn7BLmeu9xoHu2wcQqSip3YUmJznePn0DxaAdjtKq/m8rRTzz0urBSknPkFIRR35EFhpWB5ZCFFc1yCU6e6raFoTmGlIHrfFzSXh40LzjivpWdC0ZzJrkEKw+5OTg2Po30m/nF5Xu2prb8ghA5FfVOeVxvhXkhTnEDH1IV2a/FIJ3vC73FpHjYuzf98G00TmuK4TAwC7OQpdb7O7tqONU8uJ5/OCNDT5jX93eTATmGnX8Ui43V24SUh2X8mrjm3spU+r7lOzsJ1di0RwE6mUmfs1rlP5XKo1loaBVFLATvJOIROobkUr7WGT5E7YAfYKX2pM3aKT3Nh18TAHRSAnToWYAfYAXaAnVoUYAfYAXaAnVoUYKdE2H0S5F24sEpsZGPAjvAAdqQUYKdE2MXFzUYIcTjUb78lAXbKEsCOlALslAg7B4feuNfw3a+AnVIEsCOlADslwm7kyE64106c8AXslCWAHSkF2CkRdt27G+Je27NnKWCnLAHsSCnATomwMzLi417LznYF7JQlgB0pBdgpEXYcjnDci/j4uYCdsgSwI6UAO2XB7t37bKbXNoZNA+yUJYAdKQXYKQt2v/2WxPRaQIAtYKcsAexIKcBOWbC7/yCS6bWlS4cDdsoSwI6UAuyUBbuSW2uZXps/vz9gpywB7EgpwE5ZsDt1agXTa5MndwfslCWAHSkF2CkLdrt3uzO9Nnx4R8BOWQLYkVKAnbJgl5bmxPRaz57GgJ2yBLAjpQA7ZcEuPMKe6TULCx3ATlkC2JFSgJ2yYLdy5Xim13R1NQA7ZQlgR0oBdsqCnevCIUyvMU+GBezID2BHSgF2yoLdNHsb0Y77UJ0D2ClFADtSCrBTFuyGDu0g2nGiz4EF7EgOYEdKAXbKgl2XLgaiHffLr/GAnVIEsCOlADtlwU5fnyfacV/fjwTslCKAHSkF2CkLdhT1Px1XWrYesFOKAHakFGCnFNj9/c82sY67dCkQsFOKAHakFGCnFNg9epwg1nHHjnkDdkoRwI6UAuyUArv/fBkh1nGFhe6AnVIEsCOlADulwO76jRCxjsvMdAbslCKAHSkF2CkFdidO+Ip1XGzsLMBOKQLYkVKAnVJgt3PnYtxfzAPG1odOAeyUIoAdKQXYKQV2SUnzcX8NGdIe/+DnNw6wU4oAdqQUYKcU2G3YOBX3l6Njf/yD2+JhgJ1SBLAjpQA7pcDO13cM7q+QNXb4hzlz+gJ2ShHAjpQC7JQCO2fngbi/MjIW4B8mTuwK2ClFADtSCrBTCuzs7Lrh/jp5UvjYnSFD2gN2ShHAjpQC7JQCu0GDLHB/fXEvHP/QvbshYKcUAexIKcBOKbCztNTD/fXseRr+wcxMG7BTigB2pBRgpxTY6ehwEUIUhQT/Dn+ipcUB7JQigB0pBdiRj90nQR7uLD09HgMfQvI/hgKwA+zUsQA78rF7+SoDd5aVlb6gNr9dOx388t37bMCO/AB2pBRgRz52P/0chzsLn4G1sTHGL58+TQHsyA9gR0oBduRj9/kXYbizpk7rIajNHzZM+OSdhz/GAnbkB7AjpQA78rG7ejUYd9ZnroMFItfc3ftPOGBHfgA7UgqwIx+7I0e8cGfhm//nzeuHX5bcWgvYkR/AjpQC7MjHLj9/Ee6s8Ah7QW2+u/tw/PLcuQDAjvwAdqQUYEc+dvHxc3FnpaY6Cmrz/f1t8ctDhzwBO/ID2JFSgB352K1bPwV3Fn7uBDPc044dboAd+QHsSCnAjnzsPD1H4c46eXKFoDZ/69Y5+GVamhNgR34AO1IKsCMfu/nzhQN23ixZI6jN37bNGb+MjnEA7MgPYEdKAXbkYzdhgjXurK/vRwpq83fvdscv16y1A+zID2BHSgF25GPXr78Z7qwnTxIFtfnFxT74pY/PGMCO/AB2pBRgRz52FhbCm2HfvssS1OZfvhyEX7ouHALYkR/AjpQC7MjHjs/nIIRoWjjMSVl5KO47B4fegB35AexIKcCOcOw+fhKO72RoyMdT7j+IxFPGj+8C2JEfwI6UAuwIx44ZmpgZh/3XR/F4ysCB5oAd+QHsSCnAjnDsfngYi3tqxIhOeMqLl+l4SpcuBoAd+QHsSCnAjnDsKio34J6aPqMXnlLzMRdPadtWE7AjP4AdKQXYEY7dxYuBuKcWuQ1lJnI4FEKIx6MBO/ID2JFSgB3h2B08uBz31KpVE5iJeno8PPGTIA+wIzyAHSkF2BGOXU6OK+6pyKiZzMQOHdvgiX//sw2wIzyAHSkF2BGOXWzsLNxTGRkLmIm9e5vgib//ngTYER7AjpQC7AjHbnXIJNxTe/YsZSaOHNkJT/zu+xjAjvAAdqQUYEc4dsuWjcA9dfasPzNx6rQeeGLV52GAHeEB7EgpwI5w7ObM6Yt76k7pOmaik9MAPPH6jRDAjvAAdqQUYEc4dmPHWuGe+ua/0cxED4+ReOLp036AHeEB7EgpwI5w7Hr1Ep6L+OOPZGZiYOAEPHH/fg/AjvAAdqQUYEc4dmZm2rinPlTnMBPDI+zxxLy8hYAd4QHsSCnAjnDseDwaIcTlUqITExPn4e5LSpoP2BEewI6UAuxIxq66RngbrKmpluj03NyFeHrEphmAHeEB7EgpwI5k7J4+TcHd1LOnsej0ffs88PSg4ImAHeEB7EgpwI5k7P77bTTuptFjLEWnnzq1Ak9fvnwkYEd4ADtSCrAjGbvSsvW4m2bP7iM6/fr11Xi6s/NAwI7wAHakFGBHMnbnzgXgblq6dLjo9MqqjXi6/fSegB3hAexIKcCOZOyYY3PBqyfJ8vUWsCMwgB0pBdiRjF1mpjPups2bZ4lOf/IkEU/v29cMsCM8gB0pBdiRjF1UtAPupqwsF9Hpr99k4umWlnqAHeEB7EgpwI5k7Jjbwg4cWC46/ZNA+HxFAwMeYEd4ADtSCrAjGbvFS4bhbrpwYZXYLA2NOu6sAOwIDGBHSgF2JPra8b4AACAASURBVGM3c2Yv3E3ld0PFZhka8vGsj58a/RgKwA6wU8cC7EjGjhmR+PsfNovN6txZD8969ToDsCM5gB0pBdiRjF2PHm1xN/31LFVsVr/+ZnjW48cJgB3JAexIKcCOZOyMjITfVWs+5orNqnNQT8COwAB2pBRgRzJ2+GHYrfgcyVnTZwgP592t2ADYkRzAjpQC7IjF7v2HbNxH7drpSM51cRmM5169GgzYkRzAjpQC7IjF7rffknAf9eljKjnXy2sUnnvihC9gR3IAO1IKsCMWu/sPInEf2dpaS86t83mygB2BAexIKcCOWOxu3V6L+2jevH6ScyOjZuK52dmugB3JAexIKcCOWOykj9CZkuKI58bHzwXsSA5gR0oBdsRit3u3O+6jtesmS87dvt0Nz90YNg2wIzmAHSkF2BGLXVqaE+6jrVvnSM49eHA5nhsQYAvYkRzAjpQC7IjFLmLTDNxHdT4c9uxZfzxXbBBjwI60AHakFGBHLHb+/ra4jw4f9pKce7NkDZ47f35/wI7kAHakFGBHLHauC4fgPrp8OUhy7hf3wvHcyZO7A3YkB7AjpQA7YrGbZm+D+6jq8zDJuT88jMVzhw/vCNiRHMCOlALsiMVu2LAOuI9+/ClOcu6ffybjuWLPzwbsSAtgR0oBdsRi16WLAe6jFy/TJee+fZeF51pY1HHnLGBHTgA7UgqwIxY7fX0e7qNPgrrHIsZzdXU1ADuSA9iRUoAdsdjRNEIIaWvXMb4TDp/PQQjRNALsSA5gR0oBdmRi9/c/23AHdejYpr42pqZauE11jfjQnoAdOQHsSCnAjkzsHj1OwB00cKB5fW26djPEbZ6/SAPsiA1gR0oBdmRi9+VXm3AHTZrUrb42gwZZ4Da//BoP2BEbwI6UAuzIxO76jRDcQU5OA+prY2trjdt8fT8SsCM2gB0pBdiRid3x4764g3x8xtTXZtasPrhNadl6wI7YAHakFGBHJna7di3GHRS6YWp9bRYuEt5PdulSIGBHbAA7UgqwIxO7pKT5uIMSE+fV18bXdwxuc+yYN2BHbAA7UgqwIxO7DRun4g7ascOtvjZr103GbQoL3QE7YgPYkVKAHZnYybLXtnnzLNwmM9MZsCM2gB0pBdiRiZ2z80DcQdev17tIevoC3CYubjZgR2wAO1IKsCMTu8mTu+MOuvef8PraMCcx1odOAeyIDWBHSgF2ZGInywXDR4544TZ+fuMAO2ID2JFSgB2Z2Fla6uEOevN3Zn1tLlxYhdu4LR4G2BEbwI6UAuzIxE5XVwN3kJQ2t++sw23mzOkL2BEbwI6UAuzIxA73jp4eT0ob5v7ZiRO7AnbEBrAjpQA7ArF79ToD946Vlb6UZj/9HIebDRnSHrAjNoAdKQXYEYidjIr99SwVN+ve3RCwIzaAHSkF2BGI3edfhOHemTJF2mMSP1Tn4GZmZtqAHbEB7EgpwI5A7K5eDca985nrYOktKQohhLS06h26HbBjPYAdKQXYEYjd0aPeuHcavIBOR4eLWwJ2xAawI6UAOwKxy89fhHsnLNxeest27XRwy3fvswE7MgPYkVKAHYHYJSTMxb2TnDxfesuePY1xyz//TAbsxJJ9/aVf/NG5vrG281fYzl9hv2Td8sjCLUe/AezUtAA7ArFbHzoF905BwRLpLYcN64BbPvwxFrBjsibrss3QyfX9zRtaWDsHJuXefAPYqVcBdgRi5+09GvfOiRO+0lva2XXDLaWMF6BW2CWe+KHrQFtZ/vJ52vruYdsLKj4BdupSgB2B2Dk69se9c+NmiPSW8+b1wy1Lbq0F7Lxj91M0h/nbpjka/cfPdQlO9U8s9ks4tixil+38FbpGFqJ//71H2Wdc+hOwU4sC7AjEbuLErrh3vvxqk/SW7u7Dcctz5wLUHDsHj3Dmr5rm8hxXxufceCXZrKDiU0jWJYtuA5nGuobtth77L2Cn+gXYEYhd/wHmuHceP06Q3tLfX/iV7dAhT3XGzn7JOuZP2nrg+NSzj6S3L6j4tGhdFkXReBENvm7Uvs8BOxUvwI5A7Dp0bIN755+3WdJbyvKoCpXHzmlVIvP3PH6+767yahkXjN73BV/XEC+oodW6hfbvADtSCrAjEDttbQ5CiKYbvlR469Y5uB9TUx3VE7uA5BPMH7PdZ4GNXTzp1E+tjdvjxXUN22VefgrYqWwBdqRh9/FTHu4aQ0N+g423bXPGjaNjHNQQu5Qzv3B5WvgTGDjRUb5Tq8mnf9bWN8VvYtlr+I7S94CdahZgRxp2z1+k4a6xtjZosHFR0VLcOGSNnbpht+tuTUeboXjzO9gM2X7nndxvFbXvc+Y07oxlYYCdahZgRxp2D3+MxV0zfHjHBhsfP+6LG3t7j1Y37FxD0vG2c3j8pFM/NfHdfGIPMP8oInbfBexUsAA70rCrrNqIu8Z+es8GG1+5EoQbNzg+iophl3ruMc3h4W33itnbLO853N4Nv6Fhuy7b77wF7FStADvSsLt4MRB3zSK3oQ02LisPxY0dHHqrFXb9xwtvH+42aEJzqZR9/SVz8G6yaxBgp2oF2JGG3aFDnrhr/P1tG2x8/0Ekbjx+fBf1wS4o7Qzeaoqi44u/ay6VCisFIdsuMv80IvdUAXYqVYAdadjl5i7EXbMpckaDjX99FI8bDxxoribYFdz9aGQhvMNkjnd0M0qHM3aO8Gm8nfuMBOxUqgA70rDbskV46Vx6+oIGG794KTxI36VLw6duVQM7z+g9eJN1Dds1+2UihZWCbVefMZezrIg/AtipTgF2pGEXssYOd82ePUsbbFzzMRc3bttWUx2w23W3Rt+sM97kZRG7ml06nM9WC6/+0TVs1/QzFYAdKQXYkYadh8dI3DWnT/vJIheHQyGEeDxaHbBbslH4RF19s8677ta0EHY7yz4YmFvhXzR/5VbATkUKsCMNu7lz++KuuXVbplGb9PWFV2B8EuSpNnY7yz7oGJjhjfWO3d9C0uEw50C4PK06R08B7JSvADvSsBs/vgvumgffRMkiV6dOwlED3vydqdrYeUYV4S017tij4O7HFsWusFJg1W8M/nVzfGIAO1UowI407Hr1MsFd88cfMj1WondvYfvff09SbexMLXvhLfWM3tPS0hVWCkK3l+Bf18SdO8COlALsSMPOzEwbd837DzI9MGzkyE64/Xffx6gwdsEZ5/FmarZuu7PsgwKwK2ymnTvAjpQC7EjDjsejEUIcDiWjXNPsbXBXVlZtVGHsug2agDfTaVWiYqQrrBSsy72GfymXp51z8zVgp9wF2BGFHXMpiYmJrJeSODkNwItcv9HAAyuUF7vo/ffwNtJcXhNPFzQ2zM6d3MgCdqQUYEcUdk//SsX9YmNjLKNcjb1URRmxGzN7Od7GKYtWK1K6QpEbyOT++gzYkVKAHVHYfftdDO6X0WMsZZQrMFD4/W7/fg+VxC7r2gtmsLnEkw8VjF1hpcDMqg/+7csjCwE7JS7AjijsSsvW436ZNauPjHKFR9jjRfLyFqokds5ByXgDe4+aoXjpCisFXjF78QoYtusix2DIgB0pBdgRhd358wG4XxYvGSajXImJ8/AiSUnzVQ+7gopPeqad8AYGZ5xnBbtd5dXMxcxBaWcAO2UtwI4o7Pbv98D9EhQ8UUa5mFFSIjY1PEqK0mHHXHGiZ9pJvkdMNEtcgoXHUq0HjgfslLUAO6Kwy8pywf0i+wN09u1rtI9KhF2/sbPw1jkHJbMlXWGlIPfmG45GK2HX7L8H2CllNRd2v164UBAVvWLBgmmjRg3s3r2PtXUfa+sB3bpNGTHCY/bs5ODg8j173ldUAHbSGYqOccD9sm2bs4xynTq1Ai+yfPlIFcMu9ewjvGkURW+7+oxF7AorBVMWrcYrM3KGO2CnlNVE7J5cvhy23NO6fXtZfheXw5kweHBhdMzfpaWAXZ0JCp6IP6t9+2Q9tXr9uvAf4YIFA1QMu1nLN+FNG27vxq50hZWC5NNCUCiKzrj4B2CnfCU3dr9euODj6EjTtBy/tI2W1np399+vXgXsxOLuPhx/ROfPB8goF/OAnmn2NqqE3a7yas3WbfGmhe0qZR27wkrBwAnz8frMWr4JsFO+kg+7gwkJWjye6PuY6Ou7zZhREBV9u7Dwh9On/7x27c9r174/deru3r2HExPXLlkyvHdvmqJEF9Hl8wuiogE70cyaJbykq6w8VEa5mEvzRo2W9dI8pcAuIOk43i5Ty16sM4cTtvMOXiWelp7sj6kF7EipxmL3trx8ycyZou8wZ/z40qIiWZb9u7R0T2zcsF69RBe3HTToyeXLgB3OqNGW+GP59jtZ7+p/8iQRL9K3r5kqYddjqHDEZrfQHNaZY9Khx2C8VrKPkwzYkVKNwu71nTtTRoxglh3Uo8dXR4/KsV9WvmfPYBsb5n3aGxvfLy4G7AS1+TY2xvgzefo0RUa5Xr/JxItYWuqpDHaJJ37AG0VzNOS+A78lwjxL28SyJ2CnZCU7dm9KS0f168csGL7c893du/J9Ca2pqvpQUZGxdi2XI7wNSIvHu7N7N2BnYqKJP5DqmlwZ5fokyMOLGBjwVAa7qW5r8EaNnePFOnCi2VVeraUnHEBwTfYVwE6ZSkbsqisrnadMYZbaEREhN3OiqTpwwFhPD7+nFo8n+36iqmInxwMlBLX5Gho0QojLlXVUKMKx21H6XkOrNf6riNr7OevAiWX+yq143XqPsgfslKlkxC5mxQpmkUMJCc0iHc5P585ZtWuH39nM0PDXCxfUFrv3H7Lx52Burt0o7IyM+HjBj59kegwF4dh5b96HN6eDzRDWaZNM1tXnzMAEcUceAHZKU7JgV7JrF9M+ITCwGaXDeXjmTNs2wgcpjOzbV5Zrj1USu99/T8IfQu/eJo3CzspKHy/48lWGCmDHnASQb5QRBcTW0U/45zpzKWCnNNUgdi9v3epoaoobu06b1uzS4Xxx6BBz/C7C00s9sXvwTRT+BMaN69Io7PoPMBf25uMEZccuorAcb0ujLu9QcBJPPhT++6GotPNPADvlqAaxW+kivFvT3NDweUlJC2FXU1W1O2Yzs1YNXsuiktjdvrMOb/7cuX0bhd24cY17IBnJ2A2dLPx7c/AIZx01KRk8aQFez6luawA75Sjp2H1x6BDT8mp+fstJh7No+nT8u3paWko/1auS2J05448338ND1rtccWbOFF66WH5XpkuRicUu7fwThPCV5w3vMbGbqL2f48+8wZHiATtSSjp2Y/r3x82WzJzZ0tLVVFX9deOGoa4u/o3xq1apG3ZFRUvxtoessWsUdp+5Cg9yXbkSpNTYMVecDJ36GeucNRjmGUCzPCMBOyUoKdgVp6TgNq00NB5fuqQA7Gqqqg7Ex+NfyuNyf7tyRa2wS08XfjOKi5vdKOy8vUfjBY8f91Ve7LKuPqc5wnsQI/dUsW5Zg1mfdwOvLUejVda1F4Ad6VUfdu8rKrr+O5ZJS5yBlRLbQYPw7/V1clIr7DZFzsAbnpsr6wDrOCFrhHdWFRUtVV7sZnlG4q2wGTqZdchkTNeBtnidZywLA+xIr/qwy90oHEvDzNCw6SMyNSpVBw4wq/egntvIVBI7f3/hv5zDh70ahV1klPBu5exsVyXFLrfkby5fB29F6PYS1hWTMczQABSHm3HpT8CO6KoTu9d37hj9e+Gb3AOTNCUL7YUPkZkzfrz6YLdw0RC81ZcuBTYKu5QUR7xgfPxcJcVujnc03gSrfmNYJ6xRYQYssHX0azp2Z+IfF8c+PJ/0B2DX/FUndsz9Ej0tLT803/DCsufHs2eZ8aBu7tqlJtjZT++JN7mickOjsNu+3Q0vuDFsmjJil3r2EXNPgow3nJKTTUWVzL+mOkdsbxC7K6kvAh1T+3UeS6H/Hx2SpjWGdLPb5L77dnY1YNc8JYnd40uXNP69vvdMZqbipcNZs3gxXod+XbtKgquS2A0f3hFv8g8PYxuF3cGDwmdIBwTYKiN2w6a64vVn62GJTQzzDO/OvUdIPhVICnZ3cmr85yXStIaUf6GG2uahC/PLcgWAXVNLEjuP2bPxrIlDhrAlXU1V1bObN/W0tfGabA8PVwfsunYzxNv7/EVao7A7e1Z4gd7SpcOVDjvmsBdC1NZj37IulxzZduUvLk8Lb4PkLW71YVcU9nmHtj1E/zHSFMfKtM8Aa9uOxjaie3kIoZE20y8mPwXsmlRi2N0qKGBmfXHoEIvY1VRVbQ8Px2uip6397OZNlcfO0LBx9/MzuVkivDxt/vz+yoVdzo1XrY2FJ/2nLFrNOltyx229cBAHmqOx5eg3DWIX7bGPpjjMv7XB3Sal+Z+/lfWBaVCy7W2Mx/5Oxj2ZNoba5keivwXs5C9R7N6UlnbvKPwmJeWyD4XlQ0VF/65d8fr4u7ioPHb4eR5aWpxGSSeozf/invD/CpMnd1cu7JinOugatpN+HwLh2XW3xqrfGLwtbTt0zy35uz7synIFS6ZuYP4BttE0ifcprm+XrSxXsMYli9nL09RoXRT2OWAnZ4lit3iG8DovozZtxPak2IrogCunMzJUGLt/3mbhzWzfQbex2P3wMBYvO3x4RyXCjnnyNEJURGE562A1MekXfmulY4C3Z+BEx113aySxu5X13m6QC/MnPcBq/KWUvxo8/7An/J4uX3iIQ4PDL9xQCdjJUwx2693dmYksnpeQjNc84f/8NXm88j17VBW7x48T8Gb2H2DeWOz+/DMZL9uzp7GyYOcaks78vbkEpbBOVbNkXd51ZqOGTnbBo7Yw2J2Of9TNYhDTYN4Y3zs5NTJeXHI89se2uh3wglwOP2f1DcCu0fXo4sX/njgxe9w4ZkqUry/rwInm79LS3lZWeN14XG7uxo3vKypUD7uvvhY+I3XixK6Nxe7tO+FeoYWFDvnYpZ59NMTOmfl7GzxpgeQZTOWNKOLtuvQN21WaGHHldnZ1iPO2VlwdZlbA/CQZmWNyOv6RiZ7weUwUolcvyCjN/QTYNaKYUdFx+bu4VFdWsg6cWH4+d66zuTmzklwOx9LMoI+1mdKlR9e2hw551gnQjZsheOscHWU6ySAWvKyurkaDLHp4jDTtZWzR14yVGHQ0Ev17GzjRcVd5NetCNW9cglJEt1FHy1iDw///v16al7TiVGOlwzmb8ET0lEUbTZMuZn27mPUVO6sL1UDRNJ0YFMS6a/Xl8aVLI/v2ZftDaobi8egP1TmSDJ044YsbeHmNkgM7Pp+DEKJpJL3Zrl2L2dz4/625PptVaZ9ONMHp55jHaIhWn04j6zujKmOupr0cYWOv+M5Skepsbu7j6PjL+fOsi9Zg9m/ZOqZ/f2ZMYyWt23fWSTJUULAEz10fOkUO7ExNhdd51Skpk0VuQ9nbbmFxNFqNmL44vvg71klq0eTceOXgEa5rKHy+Su9OI6ScdW1sNi8/2NHYRvrnDCVeD06cYJ2wxuZ9RcXhpIBvT8YpV/xdhMOfxcbOkmSIub81IUGm+1vFwlyQ/Oy5tAuS23cQDhcYVb4+9YfNik/IsU07yz6wLpEisyn4iCy3fMmR6xlvimMfFsc+zAu5qXg6lK8a9ZBscqKMJygu5gThz3zSpG6SDG0Mm4bn7tjhJgd2gwZZ4MV/+TW+vja//SZ8oI+OqRbrZ2PVJzDqCSkF2Cksr0uF50y5XEryy6avr/CS1GPHvOXAztbWGi/+9f3I+towt9AOWzAAsAPs1K4AO0XGppsx/tglD9s5Ow/Es65fb+BB2nVm1qw+ePHSsvX1tVm5cjxusyhtPmAH2KldAXaKjKuj8PyA5GG7yZO741n3/hMuB3ayjIXXp4/wkZjR5esBO8BO7QqwU2QSo4V3g0geths8WHg/vJSDblLS4Lfg128ycQMtLc72lymAHWCndgXYKTJnD63Cd/tLHraztBRe3f36TaYc2K1dNxkvXljoXmeDCxdW4QYzZ/Yi4XYx9QlgR0oBdorMpeLVI0d2wp+82GE7XV3hCI5ySCeozd+8eRZePDPTuc4GzNne+Pi5gB1gp44F2CkYu7Bw4eXvYoft8MTWrRu436u+NPgYxrFjhfcX376zDrAD7NSxADsFY3flSh1X2716nYEnWlnpy4cdcx9YnTdgfKjO4XAohBBNo/cfsgE7wE4dC7BTMHZv32VJHrb7+ZctuDuGDGkvH3ZHjnjhd/DzGyc5t7RsPZ47arSlgO1RT1jXB7BT0wLsFIydoDZf8rDd51+E4SlTpsg01LBkmPMPbouHSc5NSJiL527YOBWwA+zUtAA7xWMnedju6tVgPMXFZbB82N2+sw6/w5w5fSXnOjj0xnPPnQsA7AA7NS3ATvHYSR62O3rUG0+p80uoLPnyq3rH/vwkyGvdWniq98XLdMAOsFPTAuwUj53kYTvmKddh4fbyYffTz3H4HSSP+j34JgrP6tXLBE8B7AA7dSzATvHYCWrzR4zohD//O6XrBLX58fHCY2opKY7yYffXM+HDa7p3NxSblZ+/CM9asWIcYAfYqW8Bdqxgx1ziiy+LW7d+Cn5ZULBEPuw+VOfgdzAz0xabxQzYeeDAcsAOsFPfAuxYwe7yZeFhOzu7boLafC+vUfjlyZMr5MNOUJtPUQjV9dhZZsDOXx/FA3aAnfoWYMcKdmKH7Rwd++PuuFmyRm7sdHS4+E1EJzIDdrZr9/8PHgPsADt1LMCOFezEDttNmCAcevOrrzfJjV27dsIn9b17n81MPHTIE090XTgEsAPs1LoAO7awEz1s16+/Gf75yZNEubGzsRGODPr0aQozkRmwMzd3IWAH2Kl1AXZsYSd62I45rPb2XZbc2A0bJnxi/MMfY5mJzICdosO1A3aAnToWYMcWdv+8zcKnFLhcisejkQxPfZUeO7tuuE+ZsY6Z8QW0tTmfBHmAHWCn1gXYsYWd6GE7XIaG/KZgN29eP/w+t26vxVOYG2anz+gl2hKwA+zUsQA7FrFjDtvh6tpN/HrgRmXxkmH4fc6fDxB7//j4uYAdYKfuBdixiB1z2A7X8OEdm4Idcy7i8GEvPEV0wE7ADrBT9wLsWMSOOWyHS+zLZmMTumEqfp+dOxcLJAbsBOwAO3UvwI5F7MQO2y1yG9oU7LZsmYPfJy3NSSAxYCdgB9ipewF27GInetguIMC2Kdht2+aM3yc6xkFQm5+YOA+/DN0wFbAD7KAAO5axEz1stylyRlOw273bHb/PmrV2AokBOwE7wE7dC7BjFzvRw3bp6Quagl1xsQ9+Hx+fMZIDdgJ2gJ26F2DHLnaih+2KipY2BTtmJ9F14ZBv/huNf+7Z01iyJWAH2KljAXasYxcbK3y+dfnd0KZgV1Yeit/HwaG35ICdgB1gp+4F2LGO3fsP2ampjkeOeDVFOkFt/v0HkbhPx4/v4rZYeIHx/v0egB1gB4UQYEcAds2VXx/F4z4dONCcGVngl1/jATvADgohwE6FsHvxMh33qbY2B/9gbi4+RDtgB9ipbwF2KoNdzcdcsc4VHbATsAPs1L0AO5XBTlCbj+8PYyonxxWwYz2AHSkF2KkSdnp6PNHOrW+Qd8AOsFPHAuxUCbsOHdswPaulxfn4Ka/OZoAdYKeOBdipEna9e5swPStlDBXADrBTxwLsVAm7kSM7MT27descwI6EAHakFGCnSthNndaD6VlmcHbADrCDQgiwUy3snJwG4G6l6f95eixgB9hBAXYqhZ2Hx0jcrZIDdgJ2gJ26F2CnStgFBk7A3bo+dApgR0gAO1IKsFMl7HbuXIy79WbJGsCOkAB2pBRgp0rYfRLkHT3qfflykPRmgB1gp44F2KkSdjIGsAPs1LEAO8AOsAPs1KIAO8AOsAPs1KIAO8AOsAPs1KIAO8AOsAPs1KIAO8AOsAPs1KJkxO5DRcUf1679cPr0V0ePVh04wHoKYryqDoSzm3uHN90vjv75fMLLW5nNjt2bvzMfP0749ruYr77edO8/4c2byDtr4+6FsZKAPeuj999T3mw+9NXWY9+mnnucW/I3YKdkJR27DxUVv5w/f2f37pPp6cUpKeRk66oFxcl+5ORcVuC9w5ue3cxoInavXmd8fT/y0qXAkydXtFxWHloWfNKXlbgleq1KPqkaCcm6FHvo65ybrwE75aj6sKuurPzh9Onz2dmsu6YU2DG5VbheCnlSsHv1OqP8bmiLGgfYtVAidt+VQh5gR0rVid3zkpKr+fmsi6aM2BUn+x1P9bt3eNP7ijwZsfv4Ke/BN1GnT/spRjrAriUSmHo67siDgrsfATtySxK7n8+dI+1Lq3Jhh3MlP+R1aVaD2L17n32zZI3CmAPsWjSh20u233kH2BFaYth9d+oU65CpBnbFyX7nswPFzl2IYffP26wrV4IULB1g16JZk30l79bfgB2JJYrdj2fPsq6YKmFXnOx3ITtIdP9OFLt377NZkQ6wU4B3ovt3gB0pxWD39Pr1E2lprCumYtgVJ/td3R7CHL9jsPskyLt9Zx0r0gF2CsjGnbcLKj4BdmQVxu59RcUFUk+8Kjt2xcl+Xx2NEsPu2+9i2JIOsFNMth77FrAjqzB2Xx45wrpfKozdibSVL0oyGez+/mfbmTP+gJ1qJyjtDD54B9iRUo8uXvy7tJT8069KjV1xsl/Zng0Mdp9/EcaidICdwhK173PAjqB6dPHiV0ePso6XymN3PNXv5e1tl4pXv32XpchL6gA7FhOYenr7nbeAHSn164ULxN4moUrYFSf73S+OvlS8+vsfNrMrHWCnyCQc/x6wI6W+OnaMdbkkI7meGhxOP2vrSC9vUezq3CK+BmdQj/YxK2azrptYruSvvlS8WvGXECsMO+/dSwbM7q1tokXRCFFI21Rr0Ny+PkXuom2k/zWyblOzJ3R7CWBHSpUUFrJOW33YMS8PbNkat3KlmaEhQmiVq6sYdqKg7I1dHuE106g1HyEUtGgS68CJ5czBANa/w7YQMRFEoQAAC/FJREFUds7xc2gNmq/Pcwid4nfAY+XB5TPXTebp8Wge/VnyvDr37HD3se5RyyblVELYJcCOiDqXlcU6bQ1ih5MYFIQQMtTVlYIdztZV8xBCRrp81nUTy+7sZaxL1xLYrdi3jMPn0Bq0V8Fi0emeO90oDsXV5vgd8FBT7JJPxqw9AdgRUceSk1mnTUbsDmzZykyXjt3+LZ4IIYqqYxa7yYz/jHXpWgK7EQsHI4SGOg+UnDXYsT9CaNSSoWqLXVjgQcCOiNofH886bTJil7BqFULI1MCgQeyifGYhhDqb67Oum1iSohxZl64lsDO00kcIuabMl5zlkjgXIWTc1VBtsQv13wfYEVF7t2xhnbYGsTuYkLA1IMDCyAghtG7JEinY7YvzDPecqafNoyhE4DmKLRFzWZeuJbCjuRRCaOXB5ZKz/A54IIRoDVptsVu3cjdgR0SRjJ1o0TTd09Iy0scHN5ByNpZCaPxA6/S1LqzTVgd24aqJHf7kg477SM5adcwLz1Vb7Nb6AXZkFMnYSWlQ357d0J4dEEKpIYRehaeq2OE9uxX763hbz51uSL337AA7UkrFsEsKdkQIDerRnnXX1Ao7I2sDhJBz/BzJWc7xc5B6H7MD7EgpFcOuONmvt5UpQighyJF12tQHu5GLhyCEBszuLTkLn40du2w4YAfYsVyqh12M32yEUB9rM9ZpUx/sVuxfxmnFoWi0NNdVdLrXrsW0Bq2hw115aDlgB9ixXKqHXXGyn3U7A4RQnP9c1nVTE+zw11V8/fD0NXb4DopZG6fy9XkcPmdRmiPTDLAD7FgrlcRuo8d0hFD3jm1Z1019sAs+6etVsLivfU++AQ/3i5YJf9Dcvr57loq2AewAO9aKTOwajDKOeqLy2DGx9R6NEJqzabrkLHUb9QSwI6gAO8Cu2eOxfRFCqM80G8AOsCOoADvAriViNbITRSOn2FkBhz0BO8COiALsALuWiN8Bj/YDzRFCNI8G7AA7Iooc7JDUMxJyY4fqOoMB2DWLaOjfm8BkjyzYIdU6awHYkVKEY4fHOAHslAg7XTNtwA6wI7GIwu5YcjKHpufZ2lIU5e3oOGPMGIRQaz6/OCVl47JlmjyeDp8f5eNzNCmJQsjCSKdos4foOE5W7QyKNntMGtKNphBNUytdJhCHXcQc1qVrXuyCTvhQHGrw/P6IQhN9x/Zz6IUQ4rXWCD7pOzvcnqvF0dDhzo+ZGXTch+JQrS10FkQtMTC3YiAwbNfFe/O+HkPsKIqmaM4klwAVxA5GPSGkijZvZp050T07hNCGZcsy1q7V4HBEd/eM2rTZGhAQHxDQ0dQUT49fNb842W/aqF7hnjOLk/02ec2cNqoXw0pW6EKapkjDLi5sNuvSnTy5wneve3Nhh/87K2zqkm3O+PAc+nd3T9OI75Iw1yVxrl6n1nj6Z0nzPotz7z3KfrZn5Krkk7O9onuPsmdQWByaS9Ec1cNutc/OlsbucPS3bAGiTLU9IoJ15sSww4MnMy+ZubgoplmSb3GyX16YG74zbEA389wwt+Jkv6CFk7gcGjcmDbuIkJmsS3f4sOeyPNfmxS7ohI/oS2Yu02dMs3nhTu5hO9p16bsq+aRFt4HuYdtXJZ+csnA1zdHAbVUPu5XLMlsau8INlS3shEpUQlAQ68yJYSf2smhzbHFKSts2bZJEVhWJEDaqn2XI4skj+1jilzp8bvyq+QmB82kK7Yp0RyRh5714HOvYZWY6O252aF7sxF6u2Lcs+KSvVlu+6NjFeO4Er7ErE49Z9RszbfFaqz6jMAdcvs6CVUnOgckURS+P3I1UCzun2WtaGrutPsfY4UO5KnDRItaZk4LdounTaZouTkkJX+6pw+dTFDVjzJjdsbFIhLDk1U40TSX+O8zJSpcJNIWWOoyYNa4vRSFEEnYjBlqxjt26dVNtvUe3HHajFg+lOFTwSd85m6Zr6HARhfo59GLm9p3Wy21d9mer0yia4xKUgjmY5BJAUfRoB4/+4+YgikKqhV3/3hNaGju/uVtY8kOpavLw4awz19isX7ZMV5t3LNmXdbwalb1xXhRFFRUtZRe7GTP7dBtj1SzYyRGDDvoTnVawDpDC4hW9l6bpK6kvWxS7sf1msQ2JMpSutvahxETW/WpUpowYgRBKDHRi3a9GJcLbASEUGjqNXezMzVu30uWtOualeOk8dyxCCHXuPYx1gxQWh2UbEELxvsUtJ93trGpdbQN2GVGaivHzY90v2XM4IUGrVSuE0IyxfVn3q1EZM6ArQmjw4E4sSpeYOA93+uyN9orHbrTbMIQQomivmL2sM6SYdO41FCE0ZehnLYddiv9Z1uxQuhrUsyfrhMkevwUL8Gq30uDs3ezBOmEyZvumxTRN4TXPy1vIFnbjxlnjdejQr52CpVt1zEvHSAv/9lHT3VhnSAFZujGfoiiEEE9D80LS0xbCblz/2Qo3Q5krZfVq1hWTcbfO1NCQWW3XacNYV0zGzBjbl1ltOzsbVqTLznYV7fQ6HxnRcpnkO5b51a20dX3jDrCOUUun1zA7ZpOXz9jUEtIdivovTdMtaYPKVfdOnY4mJbFuWYNZaG8vutpcDpW70Y11yBpM+loX/H94ppKSWHha9qBB7UXXwdjKKLDYWzHSrdi7rJUuT/S39x83i3WMWjQuwami26vZSudM/JNmxw5OTchTy+fPZ90y6UkJWcPlcMRWu5eV2ZFEH9Y5k5KD8d5W7Y3FVrtdO71DhzwVKZ2Pz1gkUSNchygGu+7/fn0WrTne0ayT1EJZseWQgWkHse0d139280qX4Hu8Cf/i1bhomt66ahXrotWX3Zs3mxkZ1bnm8ycMZF20+nI0aYXd8J51rratbbfiYh/FSJeYOE/ifxPCmhtZx3jCzZtxy0bW+as1dfWXRuxgHaZmj39ScY9B4+vc5I2LdzaXdCe3/KqvK/4/UShZi9+qVfratay7Jpn9cXHdO3WSsuYec0ax7ppkjiX7Ok8bKmW1580bcPKkb0tLl5Xlovu/XyFFi8PjfJY4r+Wkm77arr5fjRDSN7bwjFapM7MByScHT5hX3/ZqcFvlrL7RdOmupb/p1mGAlA8WquHS5PPjAwNZ1000u6KirCwsGlxz0k5WHE3wnTOh4T/H6dP7FBd7t+Q+3Xwp0uHi8DjzIme0hHSTvOv47ixWbYza4TtkVSArE472Gz1d+vZq8nWb6N2V1Je9rUY0+MFCNVw0TXs7OR1NTmadueKUlM1+Kw10dGRc8xF9rAi5GKUgetmAHh1lXO0+fS0KC91bQjpv7zH1fXuVrFGLhgaeaLbzFf6Hl/ea1F3GX63B13HwCGOdqibGI7Kwfde+DW8tQhoarcKW7JJPukNR/+1k1kPWToWSpXpbW2euX88ic7tjYqaOrPtYj5TS0+UHLrRj8U6yowm+gYvstLVaNWq1NbV4AQETjh5ttl28zEzn3r3bNfbTM7Y2ckmc23TpZm2Ypmus3djf3m3guGURhaybJd8OnZ1zAI+v1ajtnTTI6VziH7IzdyfnY5BzGr9Voz9YKJlq3KBBicHBCt7Ly1wfOtvWlsflyr3a7U30VzpP2B+3XJHM7YtbHuw22bytntyrbWraJijIbv9+D7mNO3HCJy5uzujRVnKvA0LIcnDHeZEz5NjL8z+8fHrIZGPrus8jyVQUp88o+4Uh6az7JWN8YvdPcQ1qY2gq3+ZqttJZMm39yS2/Smfu5rZ3m5YWwQ6dIqqtgcHciRM3enpuCw3dGxd3uFnvpT2WnLw7NjZ9zdq1S5bMsbXt3K7R+yP1lQaXHtyz09LZo6J8Z20PX7K3ue07uNW7MGZpUrDTSpcJowd05XJl/sYotTgcetSoLitXTkhOdiwqWnrkiJcU3Y4c8dqzZ2lS0vzAwIl2k22MjBq3cyGlNPX4NhO6Tg4Y/1niPJ/dS/wPL5fUbeUhD+/CJU6xs2y9R1uPtOLy5f//k1gZmHToN9bBfvEa19XpXpv3rkw8xrprOH5bjnhGFjkHJU9asNK63yiaboZOp2l6cPcJK+ZuzQy8fDr+8Y2Mf25n11xI/qtgQ0XYkl1Th7oq3a2v/wd6SXsD949WxAAAAABJRU5ErkJggg==" } }, "cell_type": "markdown", "metadata": {}, "source": [ "![%D0%B8%D0%B7%D0%BE%D0%B1%D1%80%D0%B0%D0%B6%D0%B5%D0%BD%D0%B8%D0%B5.png](attachment:%D0%B8%D0%B7%D0%BE%D0%B1%D1%80%D0%B0%D0%B6%D0%B5%D0%BD%D0%B8%D0%B5.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examples of **original/filtered** ECGs with **R-peaks** (QRS)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_ecg(X_train[np.argwhere(y_train==0)[0][0], :], filter_helper, ax=axes[0],\n", " title='Normal ECG', \n", " show_filtered=True, show_hr=True)\n", "plot_ecg(X_train[np.argwhere(y_train==1)[0][0], :], filter_helper, ax=axes[1],\n", " title='Class: Supraventricular premature beat', \n", " show_filtered=True, show_hr=True)\n", "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_ecg(X_train[np.argwhere(y_train==2)[0][0], :], filter_helper, ax=axes[0],\n", " title='Class: Premature ventricular contraction', \n", " show_filtered=True, show_hr=True)\n", "plot_ecg(X_train[np.argwhere(y_train==3)[0][0], :], filter_helper, ax=axes[1],\n", " title='Class: Fusion of ventricular and normal beat', \n", " show_filtered=True, show_hr=True)\n", "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", "plot_ecg(X_train[np.argwhere(y_train==4)[0][0], :], filter_helper, ax=axes[0],\n", " title='Class: Unclassifiable beat', \n", " show_filtered=True, show_hr=True)\n", "fig.delaxes(axes[1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_hr = np.apply_along_axis(calc_hr, 1, X_train)\n", "X_test_hr = np.apply_along_axis(calc_hr, 1, X_test)\n", "\n", "scaler_hr = StandardScaler()\n", "X_train_hr = scaler_hr.fit_transform(X_train_hr)\n", "X_test_hr = scaler_hr.transform(X_test_hr)\n", "\n", "X_train_hr.shape, X_test_hr.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_filtered = np.apply_along_axis(filter_function, 1, X_train)\n", "X_test_filtered = np.apply_along_axis(filter_function, 1, X_test)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**X_train_filtered, X_test_filtered - filtered data**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"X_train\", X_train.shape)\n", "print(\"X_train_filtered\", X_train_filtered.shape)\n", "print(\"y_train\", y_train.shape)\n", "print(\"X_test\", X_test.shape)\n", "print(\"X_test_filtered\", X_test_filtered.shape)\n", "print(\"y_test\", y_test.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Wavelets**\n", "\n", "The wavelet transforms have the capability to allow information extraction from both frequency and time domains, which make them suitable for ECG description. The signal is decomposed using wave_decomposition function using family db1 and 3 levels." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def compute_wavelet_descriptor(beat, family='db1', level=3):\n", " \"\"\" Compute the wavelet for a ecg \"\"\"\n", " wave_family = pywt.Wavelet(family)\n", " coeffs = pywt.wavedec(beat, wave_family, level=level)\n", " return coeffs[0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_wv_only = np.apply_along_axis(compute_wavelet_descriptor, 1, X_train_orig)\n", "X_test_wv_only = np.apply_along_axis(compute_wavelet_descriptor, 1, X_test_orig)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scaler = StandardScaler()\n", "X_train_wv_only_scaled = scaler.fit_transform(X_train_wv_only)\n", "X_test_wv_only_scaled = scaler.transform(X_test_wv_only)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Heart rate feature" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_orig_hr = np.hstack([X_train_orig, X_train_hr])\n", "X_test_orig_hr = np.hstack([X_test_orig, X_test_hr])\n", "X_train_orig_hr.shape, X_test_orig_hr.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_wv_scaled = np.hstack([X_train, X_train_wv_only_scaled, X_train_hr])\n", "X_test_wv_scaled = np.hstack([X_test, X_test_wv_only_scaled, X_test_hr])\n", "X_train_wv_scaled.shape, X_test_wv_scaled.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_wv_filtered_scaled = np.hstack([X_train_filtered, X_train_wv_only_scaled, X_train_hr])\n", "X_test_wv_filtered_scaled = np.hstack([X_test_filtered, X_test_wv_only_scaled, X_test_hr])\n", "X_train_wv_filtered_scaled.shape, X_test_wv_filtered_scaled.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**All our DataFrames:**\n", "\n", "| DatFrame | Type of data |\n", "|-|-|\n", "| X_train_orig, X_test_orig | Original(raw) data |\n", "| X_train, X_test | Baselined data |\n", "| X_train_orig_hr, X_test_orig_hr | Original data + HR |\n", "| X_train_filtered, X_test_filtered | Baselined & Filtered data |\n", "| X_train_wv_scaled, X_test_wv_scaled | Baselined data with wavelet coefficients + HR |\n", "| X_train_wv_filtered_scaled, X_test_wv_filtered_scaled | Baselined & Filtered data with wavelet coefficients + HR |" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_test_for_pred = X_test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 7. Cross-validation, hyperparameter tuning" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Some usefull model creation/training functions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def grid_clf(estimator, param_grid, Xtrain, ytrain, cv, scoring='accuracy'):\n", " \"\"\"\n", " CV using GridSearchCV for given estimator.\n", " \"\"\" \n", " grid_nn = GridSearchCV(estimator=estimator, scoring=scoring, param_grid=param_grid, cv=cv)\n", " grid_result_nn = grid_nn.fit(Xtrain, ytrain)\n", " \n", " # summarize results\n", " print(\"Best: %f using %s\" % (grid_result_nn.best_score_, grid_result_nn.best_params_))\n", " means = grid_result_nn.cv_results_['mean_test_score']\n", " stds = grid_result_nn.cv_results_['std_test_score']\n", " params = grid_result_nn.cv_results_['params']\n", " for mean, stdev, param in zip(means, stds, params):\n", " print(\"%f (%f) with: %r\" % (mean, stdev, param))\n", " \n", " return grid_result_nn\n", " \n", "def grid_logit_model(\n", " Xtrain, \n", " ytrain, \n", " scoring='accuracy',\n", " cv=StratifiedKFold(n_splits=3),\n", " param_grid=None):\n", " \"\"\"\n", " CV using GridSearchCV for LogisticRegression model.\n", " \"\"\" \n", " clf = LogisticRegression(\n", " multi_class='ovr', \n", " solver='saga',\n", " random_state=17, \n", " n_jobs=-1)\n", "\n", " if param_grid is None:\n", " Cs = [1, 0.01]\n", " param_grid = dict(\n", " C=Cs,\n", " multi_class=['ovr']\n", " )\n", " \n", " return grid_clf(clf, param_grid, Xtrain, ytrain, cv, scoring=scoring)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_confusion_matrix(cm, classes,\n", " normalize=False,\n", " title='Confusion matrix',\n", " cmap=plt.cm.Blues):\n", " \"\"\"\n", " This function prints and plots the confusion matrix.\n", " Normalization can be applied by setting `normalize=True`.\n", " \"\"\"\n", " if normalize:\n", " cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n", " print(\"Normalized confusion matrix\")\n", " else:\n", " print('Confusion matrix, without normalization')\n", "\n", " plt.imshow(cm, interpolation='nearest', cmap=cmap)\n", " plt.title(title)\n", " plt.colorbar()\n", " tick_marks = np.arange(len(classes))\n", " plt.xticks(tick_marks, classes, rotation=45)\n", " plt.yticks(tick_marks, classes)\n", "\n", " fmt = '.2f' if normalize else 'd'\n", " thresh = cm.max() / 2.\n", " for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n", " plt.text(j, i, format(cm[i, j], fmt),\n", " horizontalalignment=\"center\",\n", " color=\"white\" if cm[i, j] > thresh else \"black\")\n", "\n", " plt.tight_layout()\n", " plt.ylabel('True label')\n", " plt.xlabel('Predicted label')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def show_total_report(ypred, ytrain, ytest, additional_title='', \n", " show_classification_report = True, show_confusion_plot=True):\n", " \"\"\"\n", " Show base classification metrics.\n", " \"\"\" \n", "\n", " ohe = OneHotEncoder()\n", " ytrain_ = ohe.fit_transform(ytrain.reshape(-1,1))\n", " ytest_ = ohe.transform(ytest.reshape(-1,1))\n", " ypred_ = ohe.transform(ypred.reshape(-1,1))\n", " \n", " print(\"ranking-based average precision : {:.3f}\".format(\n", " label_ranking_average_precision_score(ytest_.todense(), ypred_.todense())))\n", " print(\"Ranking loss : {:.3f}\".format(label_ranking_loss(ytest_, ypred_.todense())))\n", " print(\"Coverage_error : {:.3f}\".format(coverage_error(ytest_.todense(), ypred_.todense())))\n", "\n", " if show_classification_report:\n", " print(classification_report(ytest_.toarray().argmax(axis=1), ypred_.argmax(axis=1)))\n", "\n", " if show_confusion_plot:\n", " # Compute confusion matrix\n", " cnf_matrix = confusion_matrix(ytest_.argmax(axis=1), ypred_.argmax(axis=1))\n", " np.set_printoptions(precision=2)\n", "\n", " # Plot non-normalized confusion matrix\n", " plt.figure(figsize=(10, 10))\n", " plot_confusion_matrix(cnf_matrix, classes=['N', 'S', 'V', 'F', 'Q'],\n", " title=additional_title + ' Confusion matrix, without normalization')\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
**Train LogisticRegression model**\n", "\n", "Our base model is LogisticRegression. For the multiclass we use LogisticRegression with a one-vs-one scheme." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def show_clf_results(\n", " clf, \n", " Xtest, \n", " ytest, \n", " ytrain,\n", " additional_title='',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True):\n", " \n", " ypred = clf.predict(Xtest)\n", " \n", " if show_metrics:\n", " print('accuracy: {}'.format(accuracy_score(y_pred=ypred, y_true=ytest)))\n", " print('precision: {}'.format(precision_score(y_pred=ypred, y_true=ytest, average='macro'))) \n", " \n", " if show_classification_report or show_confusion_plot:\n", " show_total_report(ypred, ytrain, ytest, additional_title=additional_title, \n", " show_classification_report = show_classification_report, \n", " show_confusion_plot=show_confusion_plot)\n", " ypred\n", " \n", "def train_logit(\n", " Xtrain, \n", " ytrain, \n", " Xtest, \n", " ytest, \n", " C=1,\n", " multi_class='ovr',\n", " additional_title='[Logit]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True):\n", " \n", " logit = LogisticRegression(\n", " C=C,\n", " multi_class=multi_class, \n", " random_state=17, \n", " n_jobs=-1)\n", "\n", " logit.fit(Xtrain, ytrain) \n", " \n", " return show_clf_results(\n", " logit, \n", " Xtest, \n", " ytest, \n", " ytrain, \n", " additional_title,\n", " show_metrics, \n", " show_classification_report, \n", " show_confusion_plot\n", " ), logit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Let's try LogisticRegression with original data - it is our **baseline**." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_orig,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_logit, logit = train_logit(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid, \n", " C=1,\n", " additional_title='[Logit]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For base LogisticRegression with **original** data we have:\n", "\n", "Accuracy: **0.905**\n", "\n", "Precision: **0.821**\n", "\n", "Ranking-based average precision : **0.924**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**HR features improve model quality:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_orig_hr,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_logit, logit = train_logit(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid, \n", " C=1,\n", " additional_title='[Logit HR]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### LogisticRegression with **Baseline+Filtered** data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_filtered,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_logit, logit = train_logit(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid, \n", " C=1,\n", " additional_title='[Logit Filtered]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "LogisticRegression with **Baseline+Filtered** data we have:\n", "\n", "Accuracy: **0.911**!\n", "\n", "Precision: **0.830**!\n", "\n", "Ranking-based average precision : **0.929**!\n", "\n", "So using Baseline & Filter functions improves model quality." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### LogisticRegression with **Baselined+Filtered+Wavelet+HR** data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_wv_filtered_scaled,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_logit, logit = train_logit(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid, \n", " C=1,\n", " additional_title='[Logit BL+Filtered_WV+HR]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "LogisticRegression with **Baselined+Filtered+Wavelet+HR** data we have:\n", "\n", "Accuracy: **0.933!**\n", "\n", "Precision: **0.882!**\n", "\n", "Ranking-based average precision : **0.947!**\n", "\n", "And Wavelets also add some improvement." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " #### LogisticRegression with **Baselined + Wavelet + HR** data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_wv_scaled,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_logit_flt, logit_flt = train_logit(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid, \n", " C=1,\n", " additional_title='[Scaled BL+Wavelet+HR Logit]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For **BL + Wavelet + HR** data with LogisticRegression results are:\n", "\n", "Accuracy: **0.933**\n", "\n", "Precision: **0.882**\n", "\n", "Ranking-based average precision : **0.947**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "#### The number of samples in both collections is large enough for training a deep neural network!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class KerasModelCreator():\n", " def __init__(\n", " self, \n", " feature,\n", " depth,\n", " filters=32, #32\n", " pool_size=5, # 5\n", " random_state=17):\n", " \n", " self.feature=feature\n", " self.depth=depth\n", " self.filters=filters\n", " self.pool_size=pool_size\n", " self.random_state=random_state\n", " \n", " def __call__(\n", " self, \n", " optimizer='adam', \n", " init='glorot_uniform', \n", " loss='categorical_crossentropy', #'sparse_categorical_crossentropy'\n", " metrics=['accuracy'],\n", " **sk_params):\n", " \n", " seed(self.random_state)\n", " set_random_seed(self.random_state)\n", "\n", " filters=self.filters\n", " pool_size=self.pool_size\n", "\n", " inp = Input(shape=(self.feature, self.depth))\n", " C = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init)(inp)\n", "\n", " C11 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(C)\n", " A11 = Activation(\"relu\")(C11)\n", " C12 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(A11)\n", " S11 = Add()([C12, C])\n", " A12 = Activation(\"relu\")(S11)\n", " M11 = MaxPooling1D(pool_size=pool_size, strides=2)(A12)\n", "\n", "\n", " C21 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(M11)\n", " A21 = Activation(\"relu\")(C21)\n", " C22 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(A21)\n", " S21 = Add()([C22, M11])\n", " A22 = Activation(\"relu\")(S11)\n", " M21 = MaxPooling1D(pool_size=pool_size, strides=2)(A22)\n", "\n", "\n", " C31 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(M21)\n", " A31 = Activation(\"relu\")(C31)\n", " C32 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(A31)\n", " S31 = Add()([C32, M21])\n", " A32 = Activation(\"relu\")(S31)\n", " M31 = MaxPooling1D(pool_size=pool_size, strides=2)(A32)\n", "\n", "\n", " C41 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(M31)\n", " A41 = Activation(\"relu\")(C41)\n", " C42 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(A41)\n", " S41 = Add()([C42, M31])\n", " A42 = Activation(\"relu\")(S41)\n", " M41 = MaxPooling1D(pool_size=pool_size, strides=2)(A42)\n", "\n", "\n", " C51 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(M41)\n", " A51 = Activation(\"relu\")(C51)\n", " C52 = Conv1D(filters=filters, kernel_size=5, strides=1, kernel_initializer=init, padding='same')(A51)\n", " S51 = Add()([C52, M41])\n", " A52 = Activation(\"relu\")(S51)\n", " M51 = MaxPooling1D(pool_size=pool_size, strides=2)(A52)\n", "\n", " F1 = Flatten()(M51)\n", "\n", " D1 = Dense(filters, kernel_initializer=init)(F1)\n", " A6 = Activation(\"relu\")(D1)\n", " D2 = Dense(filters, kernel_initializer=init)(A6)\n", " D3 = Dense(5, kernel_initializer=init)(D2)\n", " A7 = Softmax()(D3)\n", "\n", " model = Model(inputs=inp, outputs=A7)\n", " model.compile(loss=loss, optimizer=optimizer, metrics=metrics, **sk_params)\n", " return model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def grid_nn_model( \n", " Xtrain, \n", " ytrain, \n", " scoring='accuracy',\n", " batch_size=100,\n", " epochs=15,\n", " filters=32,\n", " pool_size=5,\n", " cv=StratifiedKFold(n_splits=3),\n", " param_grid=None):\n", " \"\"\"\n", " CV using GridSearchCV for NN model.\n", " \"\"\" \n", " n_obs, feature, depth = Xtrain.shape\n", " model_creator = KerasModelCreator(\n", " feature=feature, \n", " depth=depth,\n", " filters=filters,\n", " pool_size=pool_size,\n", " random_state=17)\n", "\n", " clf = KerasClassifier(build_fn=model_creator, \n", " epochs=2, \n", " batch_size=batch_size, \n", " verbose=1)\n", " if param_grid is None:\n", " optimizers = ['rmsprop', 'adam']\n", " init = ['glorot_uniform', 'normal', 'uniform']\n", " metrics=[['accuracy']]#'['categorical_accuracy', 'accuracy']\n", " epochs = [1, 2]\n", " batches = [50, 100]\n", " param_grid = dict(\n", " optimizer=optimizers, \n", " epochs=epochs, \n", " batch_size=batches, \n", " init=init,\n", "# metrics=metrics\n", " )\n", " \n", " return grid_clf(clf, param_grid, Xtrain, ytrain, cv, scoring=scoring)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def exp_decay(epoch):\n", " initial_lrate = 0.001\n", " k = 0.75\n", " t = n_obs//(10000 * batch_size) # every epoch we do n_obs/batch_size iteration\n", " lrate = initial_lrate * math.exp(-k*t)\n", " return lrate\n", "\n", "lrate = LearningRateScheduler(exp_decay)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use sklearn wrapper from keras framework - **KerasClassifier**." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def train_nn_model(\n", " Xtrain, \n", " ytrain, \n", " Xtest, \n", " ytest,\n", " batch_size=100,\n", " epochs=15,\n", " filters=32,\n", " pool_size=5,\n", " additional_title='[NN]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True\n", " ):\n", " \n", " seed(17)\n", " set_random_seed(17)\n", " \n", " n_obs, feature, depth = Xtrain.shape\n", " model_creator = KerasModelCreator(\n", " feature=feature, \n", " depth=depth,\n", " filters=filters,\n", " pool_size=pool_size,\n", " random_state=17)\n", " \n", " model_nn1 = KerasClassifier(\n", " build_fn=model_creator, \n", " epochs=epochs, \n", " batch_size=batch_size, \n", " verbose=0)\n", " \n", " history = model_nn1.fit(Xtrain,\n", " ytrain, \n", " epochs=epochs,\n", " batch_size=batch_size, \n", " verbose=2, \n", " # validation_data=(X_test, y_test_nn), \n", " # callbacks=[lrate]\n", " )\n", " \n", " return show_clf_results(\n", " model_nn1, \n", " Xtest, \n", " ytest, \n", " ytrain, \n", " additional_title,\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False\n", " ), model_nn1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# NN:\n", "ohe = OneHotEncoder()\n", "y_train_nn = ohe.fit_transform(y_train.reshape(-1,1))\n", "y_test_nn = ohe.transform(y_test.reshape(-1,1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "EPOCHS = 2\n", "BATCH_SIZE = 100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Our first **Deep Neural Net** with original data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_nn1 = np.expand_dims(X_train_orig, 2)\n", "X_test_nn1 = np.expand_dims(X_test_orig, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_nn1,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_nn1, model_nn1 = train_nn_model(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid, \n", " additional_title='[NN]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False,\n", " batch_size=BATCH_SIZE,\n", " epochs=EPOCHS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For **original** data with DNN model we have:\n", "\n", "Accuracy: **0.974**\n", "\n", "Precision: **0.957**\n", "\n", "Ranking-based average precision : **0.979**\n", "\n", "
Base DNN model is better than LogisticRegression." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### NN with **Baselined+Filtered** data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_nn2 = np.expand_dims(X_train_filtered, 2)\n", "X_test_nn2 = np.expand_dims(X_test_filtered, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_nn2,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_nn2, model_nn2 = train_nn_model(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid, \n", " additional_title='[NN Filtered data]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False,\n", " batch_size=BATCH_SIZE,\n", " epochs=EPOCHS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For DNN model with **Baselined+Filtered** data we have:\n", "\n", "Accuracy: **0.972**\n", "\n", "Precision: **0.957**\n", "\n", "Ranking-based average precision : **0.977**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Neural Network with **Baselined+Filtered+Wavelet** data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_nn3 = np.expand_dims(X_train_wv_filtered_scaled, 2)\n", "X_test_nn3 = np.expand_dims(X_test_wv_filtered_scaled, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_nn3,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_nn3, model_nn3 = train_nn_model(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid,\n", " additional_title='[NN BL+Filtered+Wavelet]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False,\n", " batch_size=BATCH_SIZE,\n", " epochs=EPOCHS) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For DNN model with **Baselined+Filtered+Wavelet** data we have:\n", "\n", "Accuracy: **0.972**\n", "\n", "Precision: **0.966**\n", "\n", "Ranking-based average precision : **0.978**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Neural Network with **Filtered+Wavelet** data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_nn4 = np.expand_dims(X_train_wv_scaled, 2)\n", "X_test_nn4 = np.expand_dims(X_test_wv_scaled, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_nn4,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_nn4, model_nn4 = train_nn_model(\n", " X_train_part, \n", " y_train_part, \n", " X_valid, \n", " y_valid,\n", " additional_title='[NN Filtered+Wavelet]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False,\n", " batch_size=BATCH_SIZE,\n", " epochs=EPOCHS) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For DNN model with **Filtered+Wavelet** data we have:\n", "\n", "Accuracy: **0.975**\n", "\n", "Precision: **0.960**\n", "\n", "Ranking-based average precision : **0.980**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**LogisticRegression cross-validation**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "skf = StratifiedKFold(n_splits=3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_with_err(x, mu, std, legend=False, grid=False, title=None, xlabel=None, **kwargs):\n", " lines = plt.plot(x, mu, '-', **kwargs)\n", " plt.fill_between(x, mu - std, mu + std, edgecolor='none',\n", " facecolor=lines[0].get_color(), alpha=0.2)\n", " if legend:\n", " plt.legend()\n", " if grid:\n", " plt.grid(True);\n", " if title is not None:\n", " plt.title(title)\n", " if xlabel is not None:\n", " plt.xlabel(xlabel)\n", " \n", "def plot_learning_curve(extimator, X, y, cv):\n", " train_sizes = np.linspace(0.05, 1, 20)\n", " N_train, val_train, val_test = learning_curve(estimator,\n", " X, y, \n", " train_sizes=train_sizes, \n", " cv=cv,\n", " scoring='accuracy')\n", " plot_with_err(N_train, val_train.mean(1), val_train.std(1), label='training scores')\n", " plot_with_err(N_train, val_test.mean(1), val_test.std(1), label='validation scores')\n", " plt.xlabel('Training Set Size'); \n", " plt.ylabel('Accuracy')\n", " plt.legend()\n", " plt.grid(True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### LogisticRegression for our best data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "param_grid = dict(\n", " C=np.arange(1, 3, 0.25),\n", " multi_class=['ovr'],\n", " )\n", "\n", "\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_wv_filtered_scaled,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "grid_logit = grid_logit_model(\n", " X_train_part, \n", " y_train_part, \n", " cv=skf,\n", " param_grid=param_grid)\n", "\n", "print(grid_logit.best_params_)\n", "\n", "y_pred_logit_gr = show_clf_results(\n", " grid_logit, \n", " X_valid, \n", " y_valid, \n", " y_train, \n", " additional_title='[Logit CV]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So optimal C value is **2.5**." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = 2.5\n", "grid_logit.best_params_" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# with open('grid_logit.pkl', 'wb') as f:\n", "# pickle.dump(grid_logit, file=f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Neural Net cross-validation**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_nn3,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "param_nn_grid = dict(\n", " epochs=[1, 2, 10, 30, 50], \n", " batch_size=[50], \n", ")\n", "\n", "grid_nn = grid_nn_model(\n", " X_train_part, \n", " y_train_part, \n", " cv=StratifiedKFold(n_splits=3),\n", " param_grid=param_nn_grid)\n", "\n", "grid_nn.best_params_\n", "\n", "y_pred_nn_gr = show_clf_results(\n", " grid_nn, \n", " X_valid, \n", " y_valid, \n", " y_train, \n", " additional_title='[NN CV]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**So optimal epochs value is 30**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "grid_nn.best_params_" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# with open('grid_nn.pkl', 'wb') as f:\n", "# pickle.dump(grid_nn, file=f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 8. Validation and learning curves" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### LogisticRegression validation/learning curves" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_with_err(grid_logit.param_grid['C'], \n", " grid_logit.cv_results_['mean_train_score'], \n", " grid_logit.cv_results_['std_train_score'], \n", " legend=True,\n", " grid=True,\n", " xlabel='C',\n", " label='Training scores')\n", "plot_with_err(grid_logit.param_grid['C'], \n", " grid_logit.cv_results_['mean_test_score'], \n", " grid_logit.cv_results_['std_test_score'], \n", " legend=True, \n", " grid=True,\n", " xlabel='C',\n", " title='Logit',\n", " label='Test scores')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### NN validation/learning curves" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_with_err(grid_nn.param_grid['epochs'], \n", " grid_nn.cv_results_['mean_train_score'], \n", " grid_nn.cv_results_['std_train_score'], \n", " legend=True,\n", " grid=True,\n", " xlabel='epochs',\n", " label='Training scores')\n", "plot_with_err(grid_nn.param_grid['epochs'], \n", " grid_nn.cv_results_['mean_test_score'], \n", " grid_nn.cv_results_['std_test_score'], \n", " legend=True, \n", " grid=True,\n", " xlabel='epochs',\n", " title='NN',\n", " label='Test scores')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 9. Prediction for hold-out and test samples " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Build model using full train data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "# grid_logit_full = grid_logit.fit(\n", "# X_train_wv_filtered_scaled,\n", "# y_train)\n", "# our train/valid for logit\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_wv_filtered_scaled,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_logit_full, grid_logit_full = train_logit(\n", " X_train_wv_filtered_scaled, \n", " y_train, \n", " X_test_wv_filtered_scaled, \n", " y_test, \n", " C=grid_logit.best_params_['C'],\n", " additional_title='[Logit BL+Filtered_WV+HR]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**LogisticRegresion hold-out results:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_pred_logit_valid_full = show_clf_results(\n", " grid_logit_full, \n", " X_valid, \n", " y_valid, \n", " y_train, \n", " additional_title='[Logit hold-out]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**LogisticRegresion test results:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_pred_logit_test_full = show_clf_results(\n", " grid_logit_full, \n", " X_test_wv_filtered_scaled,\n", " y_test,\n", " y_train,\n", " additional_title='[Logit Test]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "#### Build NN model using full train data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "X_train_part, X_valid, y_train_part, y_valid = \\\n", " train_test_split(\n", " X_train_nn3,\n", " y_train,\n", " test_size=0.3,\n", " random_state=17)\n", "\n", "y_pred_nn3_full, grid_nn_full = train_nn_model(\n", " X_train_nn3, \n", " y_train, \n", " X_test_nn3,\n", " y_test,\n", " additional_title='[NN BL+Filtered+Wavelet]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=False,\n", " batch_size=grid_nn.best_params_['batch_size'], \n", " epochs=grid_nn.best_params_['epochs'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**NN hold-out results:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_pred_nn_valid_full = show_clf_results(\n", " grid_nn_full, \n", " X_valid, \n", " y_valid, \n", " y_train, \n", " additional_title='[NN hold-out]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**NN test results:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "grid_nn_test_full = show_clf_results(\n", " grid_nn_full, \n", " X_test_nn3,\n", " y_test,\n", " y_train,\n", " additional_title='[NN Test]',\n", " show_metrics=True, \n", " show_classification_report=True, \n", " show_confusion_plot=True\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 10. Model evaluation with metrics description" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use accuracy metrics in our model. In medical diagnosis, test sensitivity is the ability of a test to correctly identify those with the disease - [wiki](https://en.wikipedia.org/wiki/Sensitivity_and_specificity#Medical_examples), so we will also use **precision** metric (we calculate other classification metrics, such as f1-score, ...).\n", "\n", "For LogisticRegression we get accuracy=0.934, precision=0.908 for test data. But NN is much better - accuracy=0.985, precision=0.943." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part 11. Conclusions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Deep Neural Network model's accuracy is 0.968! It is perfect result. \n", "\n", "But we can improve our model, some ideas for improvement:\n", "* Use additional ECG features, such as P, Q, R, S, T waves\n", "* Our data is inbalanced, so try data oversampling (see data augmentation above)\n", "* Using original ECG will be better solution (in our current DB we have only 'one R-peaks')\n", "* Use FFT analysis\n", "* Tune hyperparameters (Neural Network tuning need powewfull hardware)\n", "* Try another models (SVM, ...)\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }