{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# GQRS Algorithm\n", "### An implementation for non-linealy sampled signal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook we're going to presenet a ri-formulation of the gqrs algorithm by George B. Moody (george@mit.edu) and persent in the [WFDB](https://www.physionet.org/physiotools/wag/gqrs-1.htm) library.\n", "\n", "The classial gQRS algorithm achieve above-the-average performances in the detection of the QRS complex.\n", "In this notebook we will give an alternative formulation that can work with a sub-class of the non-linarly sampled signal while achiveing good results (when compared to the original ones).\n", "\n", "In partiular, the non-linearly sampled signal where generated from the ECG signal in the [MIT-BIH arithmia database](https://physionet.org/physiobank/database/mitdb/). For each signal, the sub-sampling is based on the idea to keep the error between the real(original) curve and the linearly interpolated signal below a given threshold. \n", "\n", "This give us a strong prior: each pair of point in the sub-sampled signal can be linearly interpolated while keeping the error below the selected threshold. We will not execute a linear interpolation during this work but we will use this fact to obtain an efficient re-formulation.\n", "\n", "It's important to notice the total absence of the original smothing filter present in the gQRS algorithm. This is because the sub-sampling can be thought as a strog pice-wise smothing.\n", "\n", "The original python implementation from which this work departed was originally posted on github from: [MIT Laboratory for Computational Physiology](https://github.com/MIT-LCP/wfdb-python) " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2020-01-06T16:22:58.586110Z", "start_time": "2020-01-06T16:22:58.392597Z" } }, "outputs": [], "source": [ "#Let's import some usefull standard library\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import collections\n", "import copy\n", "import sys\n", "import os\n", "import time\n", "import pprint\n", "pp = pprint.PrettyPrinter(indent=4)\n", "\n", "#In order to write the results found to a standard ANNOTATION file that can be checked with the WFDB bxb tool\n", "scripts = '../data_parsing'\n", "if scripts not in sys.path:\n", " sys.path.insert(0,scripts)\n", "import MITAnnotation as MITA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## File opening and showing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each sub-sampled file is saved inside the data directory, inside the directory of the corresponding threshold.\n", "First let's select a threshold, after that we will be able to select the interested file(s).\n", "If we want to select all thresholds (and, after, all files) we can insert -1.\n", "If we want to select a custom location we can insert -2 and in the file selection, insert the path.\n", "\n", "The file format is a standard CSV file with 2 columns: time and value" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "ExecuteTime": { "end_time": "2020-01-06T16:32:00.879432Z", "start_time": "2020-01-06T16:31:58.072371Z" }, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Directories:\n", "Select using the number (-1 to use all the directories):\n", "0 --> ../data/dataOut/lvlCrossingBits5Hist0\n", "-1\n", "Using:\n", "['../data/dataOut/lvlCrossingBits5Hist0']\n" ] } ], "source": [ "# We find the folder in 'dataOut' containing the data to use as input for the gQRS algorithm\n", "# change \"dataRaw\" with the output folder used for the ADC\n", "Folder = \"../data/dataOut\"\n", "\n", "dirOut = []\n", "for f in os.listdir(Folder):\n", " path = os.path.join(Folder,f)\n", " if 'lvlCrossingBits' in f and os.path.isdir(path):\n", " dirOut.append(path)\n", "\n", "print(\"Directories:\")\n", "print(\"Select using the number (-1 to use all the directories):\")\n", "for i,name in zip(range(len(dirOut)),dirOut):\n", " print(str(i)+\" --> \"+str(name))\n", " \n", "selectedFolder = int(input())\n", "if selectedFolder == -1:\n", " pass\n", "elif selectedFolder >= len(dirOut):\n", " print(\"Invalid, used the firs as default\")\n", " dirOut = [dirOut[0]]\n", "else:\n", " dirOut = [dirOut[selectedFolder]]\n", "print(\"Using:\")\n", "print(dirOut)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2020-01-06T16:23:57.018341Z", "start_time": "2020-01-06T16:23:50.918984Z" }, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ '221',\n", " '212',\n", " '122',\n", " '228',\n", " '109',\n", " '230',\n", " '207',\n", " '119',\n", " '209',\n", " '112',\n", " '115',\n", " '234',\n", " '106',\n", " '205',\n", " '124',\n", " '202',\n", " '219',\n", " '123',\n", " '201',\n", " '101',\n", " '222',\n", " '217',\n", " '231',\n", " '117',\n", " '200',\n", " '111',\n", " '220',\n", " '118',\n", " '203',\n", " '210',\n", " '100',\n", " '114',\n", " '215',\n", " '214',\n", " '103',\n", " '121',\n", " '233',\n", " '107',\n", " '113',\n", " '213',\n", " '116',\n", " '208',\n", " '105',\n", " '108',\n", " '232',\n", " '223',\n", " -1]\n", "Select the wanted file (-1 to use all thresholds): \n", "-1\n", "Opening: ../data/dataOut/lvlCrossingBits5Hist0\n", "(\"Slected file(s): ['221', '212', '122', '228', '109', '230', '207', '119', \"\n", " \"'209', '112', '115', '234', '106', '205', '124', '202', '219', '123', '201', \"\n", " \"'101', '222', '217', '231', '117', '200', '111', '220', '118', '203', '210', \"\n", " \"'100', '114', '215', '214', '103', '121', '233', '107', '113', '213', '116', \"\n", " \"'208', '105', '108', '232', '223']\")\n" ] } ], "source": [ "fileList = [os.path.splitext(name)[0] for name in os.listdir(dirOut[0])]\n", "fileList.append(-1)\n", "pp.pprint(fileList)\n", "print(\"Select the wanted file (-1 to use all thresholds): \")\n", "select = int(input())\n", "if select != -1 and str(select) not in fileList:\n", " select = fileList[0]\n", " print('not present, used: ',select)\n", "select = str(select)\n", "\n", "\n", "\n", "data = []\n", "if select != '-1':\n", " for dire in dirOut:\n", " fileThisThreshold = os.path.join(dire,select+'.csv')\n", " data.append([pd.read_csv(fileThisThreshold,header = None)])\n", " selectedFileNames = [select]\n", "elif select == '-1':\n", " for dire in dirOut:\n", " print(\"Opening: \"+str(dire))\n", " data.append([pd.read_csv(os.path.join(dire,os.listdir(dire)[j]),header = None) \n", " for j in range(len(os.listdir(dire)))])\n", " selectedFileNames = [fileList[j]for j in range(len(os.listdir(dirOut[0])))]\n", "\n", "pp.pprint(\"Slected file(s): \"+str(selectedFileNames))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Let's see one file:\n", "We ordered the \"data\" structure as it follow:\n", "\n", "[Threshold][File][DataFrame]\n", "\n", "Select a starting and a stop time (in seconds)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2020-01-06T16:24:13.222615Z", "start_time": "2020-01-06T16:24:09.429261Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Start: \n", "5\n", "Stop: \n", "15\n" ] } ], "source": [ "print('Start: ')\n", "start = int(float(input())*360)\n", "print('Stop: ')\n", "stop = int(float(input())*360)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we define a function that return us the indexes to use given the time vector:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2020-01-06T16:24:14.617418Z", "start_time": "2020-01-06T16:24:14.606408Z" } }, "outputs": [], "source": [ "def getIdxs(time, start, stop):\n", " startIdx = 0\n", " stopIdx = len(t)-1\n", " for i in range(len(time)-1):\n", " if time[i]<=start and time[i+1]> start:\n", " startIdx = i\n", " if time[i]<=stop and time[i+1]> stop:\n", " stopIdx = i\n", " \n", " return startIdx,stopIdx" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2020-01-06T16:24:16.064468Z", "start_time": "2020-01-06T16:24:15.171604Z" }, "scrolled": false }, "outputs": [ { "ename": "NameError", "evalue": "name 'thresholdComplete' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_title\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ECG at threshold = '\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthresholdComplete\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mth\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m35\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 7\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_xlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Time\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m35\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_ylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"ADC value\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m35\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'thresholdComplete' is not defined" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABrcAAAmxCAYAAAD7GqMSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3U+oZvddx/HPtxmjUKsFcwXJRBtwSh1EiF5ioQsrrTDpIrMRyUARpXQ2RhcWIaJEiSvbRUGIf4JIVbAxdqGDjGQhFUFMyQ3V0iQEhvgnlwid1pJNsTHwczE39fbmTu+T6ZM2H+b1goHnnPN9zv2uns2bc2bWWgEAAAAAAIAGb/l2LwAAAAAAAACbErcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGifGrZn5k5n5wsx8/jrXZ2Z+b2auzMznZubHt78mAAAAAAAAbPbk1ieSnPsG1+9Jcubg38Ukf/DNrwUAAAAAAACvdWLcWmv9Y5L//gYj55P82brmiSRvn5kf2NaCAAAAAAAA8Kpt/J9btyd54dDx/sE5AAAAAAAA2KpTW7jHHHNuHTs4czHXXl2Yt771rT/xrne9awt/HgAAAAAAgCZPPfXUF9daOzfy3W3Erf0kdxw6Pp3kxeMG11qPJHkkSXZ3d9fe3t4W/jwAAAAAAABNZuY/bvS723gt4aUkPz/XvDvJS2ut/9rCfQEAAAAAAODrnPjk1sx8Msl7k9w2M/tJfivJdyTJWusPk1xO8oEkV5J8JckvvlHLAgAAAAAAcHM7MW6ttS6ccH0l+aWtbQQAAAAAAADXsY3XEgIAAAAAAMC3hLgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABAjY3i1sycm5nnZubKzDxwzPUfnJlPz8xnZ+ZzM/OB7a8KAAAAAADAze7EuDUztyR5OMk9Sc4muTAzZ4+M/WaSx9ZadyW5L8nvb3tRAAAAAAAA2OTJrbuTXFlrPb/WejnJo0nOH5lZSb7n4PP3JnlxeysCAAAAAADANZvErduTvHDoeP/g3GG/neSDM7Of5HKSXz7uRjNzcWb2Zmbv6tWrN7AuAAAAAAAAN7NN4tYcc24dOb6Q5BNrrdNJPpDkz2fmNfdeaz2y1tpda+3u7Oy8/m0BAAAAAAC4qW0St/aT3HHo+HRe+9rBDyV5LEnWWv+c5LuS3LaNBQEAAAAAAOBVm8StJ5OcmZk7Z+bWJPcluXRk5j+TvC9JZuZHci1uee8gAAAAAAAAW3Vi3FprvZLk/iSPJ3k2yWNrradn5qGZufdg7CNJPjwz/5rkk0l+Ya119NWFAAAAAAAA8E05tcnQWutykstHzj146PMzSd6z3dUAAAAAAADg623yWkIAAAAAAAB4UxC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1Ngobs3MuZl5bmauzMwD15n5uZl5Zmaenpm/2O6aAAAAAAAAkJw6aWBmbknycJKfSbKf5MmZubTWeubQzJkkv57kPWutL8/M979RCwMAAAAAAHDz2uTJrbuTXFlrPb/WejnJo0nOH5n5cJKH11pfTpK11he2uyYAAAAAAABsFrduT/LCoeP9g3OHvTPJO2fmn2bmiZk5d9yNZubizOzNzN7Vq1dvbGMAAAAAAABuWpvErTnm3DpyfCrJmSTvTXIhyR/PzNtf86W1Hllr7a61dnd2dl7vrgAAAAAAANzkNolb+0nuOHR8OsmLx8z8zVrrf9da/5bkuVyLXQAAAAAAALA1m8StJ5OcmZk7Z+bWJPcluXRk5q+T/HSSzMxtufaawue3uSgAAAAAAACcGLfWWq8kuT/J40meTfLYWuvpmXloZu49GHs8yZdm5pkkn07ya2utL71RSwMAAAAAAHBzmrWO/vdZ3xq7u7trb2/v2/K3AQAAAAAA+PaZmafWWrs38t1NXksIAAAAAAAAbwriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBobxa2ZOTczz83MlZl54BvM/ezMrJnZ3d6KAAAAAAAAcM2JcWtmbknycJJ7kpxNcmFmzh4z97Ykv5LkM9teEgAAAAAAAJLNnty6O8mVtdbza62Xkzya5Pwxc7+T5KNJ/meL+wEAAAAAAMDXbBK3bk/ywqHj/YNzXzMzdyW5Y631t1vcDQAAAAAAAL7OJnFrjjm3vnZx5i1JPp7kIyfeaObizOzNzN7Vq1c33xIAAAAAAACyWdzaT3LHoePTSV48dPy2JD+a5B9m5t+TvDvJpZnZPXqjtdYja63dtdbuzs7OjW8NAAAAAADATWmTuPVkkjMzc+fM3JrkviSXXr241npprXXbWusda613JHkiyb1rrb03ZGMAAAAAAABuWifGrbXWK0nuT/J4kmeTPLbWenpmHpqZe9/oBQEAAAAAAOBVpzYZWmtdTnL5yLkHrzP73m9+LQAAAAAAAHitTV5LCAAAAAAAAG8K4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1NopbM3NuZp6bmSsz88Ax1391Zp6Zmc/NzN/PzA9tf1UAAAAAAABudifGrZm5JcnDSe5JcjbJhZk5e2Tss0l211o/luRTST667UUBAAAAAABgkye37k5yZa31/Frr5SSPJjl/eGCt9em11lcODp9Icnq7awIAAAAAAMBmcev2JC8cOt4/OHc9H0ryd8ddmJmLM7M3M3tXr17dfEsAAAAAAADIZnFrjjm3jh2c+WCS3SQfO+76WuuRtdbuWmt3Z2dn8y0BAAAAAAAgyakNZvaT3HHo+HSSF48Ozcz7k/xGkp9aa311O+sBAAAAAADA/9vkya0nk5yZmTtn5tYk9yW5dHhgZu5K8kdJ7l1rfWH7awIAAAAAAMAGcWut9UqS+5M8nuTZJI+ttZ6emYdm5t6DsY8l+e4kfzUz/zIzl65zOwAAAAAAALhhm7yWMGuty0kuHzn34KHP79/yXgAAAAAAAPAam7yWEAAAAAAAAN4UxC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1xC0AAAAAAABqiFsAAAAAAADUELcAAAAAAACoIW4BAAAAAABQQ9wCAAAAAACghrgFAAAAAABADXELAAAAAACAGuIWAAAAAAAANcQtAAAAAAAAaohbAAAAAAAA1BC3AAAAAAAAqCFuAQAAAAAAUEPcAgAAAAAAoIa4BQAAAAAAQA1xCwAAAAAAgBriFgAAAAAAADXELQAAAAAAAGqIWwAAAAAAANQQtwAAAAAAAKghbgEAAAAAAFBD3AIAAAAAAKCGuAUAAAAAAEANcQsAAAAAAIAa4hYAAAAAAAA1NopbM3NuZp6bmSsz88Ax179zZv7y4PpnZuYd214UAAAAAAAAToxbM3NLkoeT3JPkbJILM3P2yNiHknx5rfXDST6e5He3vSgAAAAAAABs8uTW3UmurLWeX2u9nOTRJOePzJxP8qcHnz+V5H0zM9tbEwAAAAAAAJJTG8zcnuSFQ8f7SX7yejNrrVdm5qUk35fki4eHZuZikosHh1+dmc/fyNIAN7HbcuS3FYAT+e0EeH38bgK8fn47gf9j7/5CLL/LO45/HrOmgo0K7hYkuzaBbqqpFGKHkOJFI7YlyUVyE2QD4h+CuWmUtiJEKirxqpYiCFG7rSFV0DT1QhdJyUWNWMRIJkiDGwks0ZohQraa5kZsTPv0YiZhnJ3snN3MzM7Teb1g4fx+53vOea6+zOx7fr/Dufvd833hInFrsyuw+jzWpLuPJzmeJFW13N1LC3w+AGvsnQDnzt4JcG7smwDnzt4JcO6qavl8X7vIbQlXkhxZd3w4yVMvtaaqDiR5bZKfn+9QAAAAAAAAsJlF4tbDSY5W1eVVdXGSY0lObFhzIsl71h7fnOSb3X3GlVsAAAAAAADwcmx5W8K179C6PckDSS5Kcnd3n6yqO5Msd/eJJF9I8qWqOpXVK7aOLfDZx1/G3AD7lb0T4NzZOwHOjX0T4NzZOwHO3XnvneUCKwAAAAAAAKZY5LaEAAAAAAAAsCeIWwAAAAAAAIyx43Grqq6rqser6lRV3bHJ879RVf+09vz3quqynZ4JYK9bYO/8y6p6rKoerap/rarfvhBzAuwVW+2b69bdXFVdVUu7OR/AXrTI3llV71z7ufNkVX15t2cE2GsW+H39jVX1YFV9f+139hsuxJwAe0VV3V1VT1fVD17i+aqqz6ztq49W1VsXed8djVtVdVGSu5Jcn+TKJLdU1ZUblt2a5Jnu/p0kn07y1zs5E8Bet+De+f0kS939+0m+muRTuzslwN6x4L6ZqrokyQeTfG93JwTYexbZO6vqaJKPJHlbd/9ekj/f9UEB9pAFf+78aJL7uvuqJMeSfHZ3pwTYc+5Jct1Znr8+ydG1f7cl+dwib7rTV25dneRUdz/R3c8luTfJTRvW3JTkH9cefzXJO6qqdngugL1sy72zux/s7l+sHT6U5PAuzwiwlyzyM2eSfDKrfwzwy90cDmCPWmTvfH+Su7r7mSTp7qd3eUaAvWaRvbOTvGbt8WuTPLWL8wHsOd397SQ/P8uSm5J8sVc9lOR1VfWGrd53p+PWpUmeXHe8snZu0zXd/XySZ5O8fofnAtjLFtk717s1yb/s6EQAe9uW+2ZVXZXkSHd/YzcHA9jDFvmZ84okV1TVd6rqoao621/cAuwHi+ydn0jyrqpaSXJ/kg/szmgAY53r/4UmSQ7s2DirNrsCq89jDcB+svC+WFXvSrKU5I92dCKAve2s+2ZVvSKrt79+724NBDDAIj9zHsjq7WGuzeqdAv6tqt7S3f+1w7MB7FWL7J23JLmnu/+2qv4wyZfW9s7/3fnxAEY6r0a001durSQ5su74cM68FPfFNVV1IKuX657tEjWA/+8W2TtTVX+c5K+S3Njd/71LswHsRVvtm5ckeUuSb1XVj5Nck+REVS3t2oQAe8+iv69/vbt/1d0/SvJ4VmMXwH61yN55a5L7kqS7v5vkVUkO7sp0ADMt9H+hG+103Ho4ydGquryqLs7qlyie2LDmRJL3rD2+Ock3u9uVW8B+tuXeuXZ7rb/Latjy3QfAfnfWfbO7n+3ug919WXdfltXvKryxu5cvzLgAe8Iiv69/Lcnbk6SqDmb1NoVP7OqUAHvLInvnT5K8I0mq6s1ZjVund3VKgFlOJHl3rbomybPd/dOtXrSjtyXs7uer6vYkDyS5KMnd3X2yqu5MstzdJ5J8IauX557K6hVbx3ZyJoC9bsG982+S/GaSf66qJPlJd994wYYGuIAW3DcBWGfBvfOBJH9aVY8l+Z8kH+7un124qQEurAX3zg8l+fuq+ous3lbrvf6QH9jPquorWb3N9cG17yP8eJJXJkl3fz6r3094Q5JTSX6R5H0Lva+9FQAAAAAAgCl2+raEAAAAAAAAsG3ELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAACi7cWrAAAgAElEQVQAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMbYMm5V1d1V9XRV/eAlnq+q+kxVnaqqR6vqrds/JgAAAAAAACx25dY9Sa47y/PXJzm69u+2JJ97+WMBAAAAAADAmbaMW9397SQ/P8uSm5J8sVc9lOR1VfWG7RoQAAAAAAAAXrAd37l1aZIn1x2vrJ0DAAAAAACAbXVgG96jNjnXmy6sui2rty7Mq1/96j9405vetA0fDwAAAAAAwCSPPPLIf3b3ofN57XbErZUkR9YdH07y1GYLu/t4kuNJsrS01MvLy9vw8QAAAAAAAExSVf9xvq/djtsSnkjy7lp1TZJnu/un2/C+AAAAAAAA8Gu2vHKrqr6S5NokB6tqJcnHk7wySbr780nuT3JDklNJfpHkfTs1LAAAAAAAAPvblnGru2/Z4vlO8mfbNhEAAAAAAAC8hO24LSEAAAAAAADsCnELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMRaKW1V1XVU9XlWnquqOTZ5/Y1U9WFXfr6pHq+qG7R8VAAAAAACA/W7LuFVVFyW5K8n1Sa5McktVXblh2UeT3NfdVyU5luSz2z0oAAAAAAAALHLl1tVJTnX3E939XJJ7k9y0YU0nec3a49cmeWr7RgQAAAAAAIBVi8StS5M8ue54Ze3cep9I8q6qWklyf5IPbPZGVXVbVS1X1fLp06fPY1wAAAAAAAD2s0XiVm1yrjcc35Lknu4+nOSGJF+qqjPeu7uPd/dSdy8dOnTo3KcFAAAAAABgX1skbq0kObLu+HDOvO3grUnuS5Lu/m6SVyU5uB0DAgAAAAAAwAsWiVsPJzlaVZdX1cVJjiU5sWHNT5K8I0mq6s1ZjVvuOwgAAAAAAMC22jJudffzSW5P8kCSHya5r7tPVtWdVXXj2rIPJXl/Vf17kq8keW93b7x1IQAAAAAAALwsBxZZ1N33J7l/w7mPrXv8WJK3be9oAAAAAAAA8OsWuS0hAAAAAAAA7AniFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBgLxa2quq6qHq+qU1V1x0useWdVPVZVJ6vqy9s7JgAAAAAAACQHtlpQVRcluSvJnyRZSfJwVZ3o7sfWrTma5CNJ3tbdz1TVb+3UwAAAAAAAAOxfi1y5dXWSU939RHc/l+TeJDdtWPP+JHd19zNJ0t1Pb++YAAAAAAAAsFjcujTJk+uOV9bOrXdFkiuq6jtV9VBVXbfZG1XVbVW1XFXLp0+fPr+JAQAAAAAA2LcWiVu1ybnecHwgydEk1ya5Jck/VNXrznhR9/HuXurupUOHDp3rrAAAAAAAAOxzi8StlSRH1h0fTvLUJmu+3t2/6u4fJXk8q7ELAAAAAAAAts0icevhJEer6vKqujjJsSQnNqz5WpK3J0lVHczqbQqf2M5BAQAAAAAAYMu41d3PJ7k9yQNJfpjkvu4+WVV3VtWNa8seSPKzqnosyYNJPtzdP9upoQEAAAAAANifqnvj12ftjqWlpV5eXr4gnw0AAAAAAMCFU1WPdPfS+bx2kdsSAgAAAAAAwJ4gbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIyxUNyqquuq6vGqOlVVd5xl3c1V1VW1tH0jAgAAAAAAwKot41ZVXZTkriTXJ7kyyS1VdeUm6y5J8sEk39vuIQEAAAAAACBZ7Mqtq5Oc6u4nuvu5JPcmuWmTdZ9M8qkkv9zG+QAAAAAAAOBFi8StS5M8ue54Ze3ci6rqqiRHuvsb2zgbAAAAAAAA/JpF4lZtcq5ffLLqFUk+neRDW75R1W1VtVxVy6dPn158SgAAAAAAAMhicWslyZF1x4eTPLXu+JIkb0nyrar6cZJrkpyoqqWNb9Tdx7t7qbuXDh06dP5TAwAAAAAAsC8tErceTnK0qi6vqouTHEty4oUnu/vZ7j7Y3Zd192VJHkpyY3cv78jEAAAAAAAA7Ftbxq3ufj7J7UkeSPLDJPd198mqurOqbtzpAQEAAAAAAOAFBxZZ1N33J7l/w7mPvcTaa1/+WAAAAAAAAHCmRW5LCAAAAAAAAHuCuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAD8H3t3F6LpXd5x/Hc1axR8BXcLkk004FpNRdAOweKBirYkOcieWElAfCG4J42lVYSIopIeqRRBiC9bFF9AY/RAl7KSA02xiJFsSBtMJLDE1gwREjXmRDSmvXowkzqZTHaeXZ+Znavz+cDCc9/3f565jv7M7Hfu+wFgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDEWiltVdUVV3VdVp6vqhi2uv7eq7q2qu6vqu1X14uWPCgAAAAAAwH63bdyqqguS3JTkyiSXJbm2qi7btOyuJCvd/aok30zy8WUPCgAAAAAAAIvcuXV5ktPdfX93P5bk5iRHNy7o7tu6+zfrh7cnObzcMQEAAAAAAGCxuHVRkgc2HK+un3s61yX5zlYXqupYVZ2qqlMPP/zw4lMCAAAAAABAFotbtcW53nJh1duSrCT5xFbXu/t4d69098qhQ4cWnxIAAAAAAACSHFhgzWqSizccH07y4OZFVfXmJB9M8vru/t1yxgMAAAAAAIA/WOTOrTuSHKmqS6vqwiTXJDmxcUFVvTrJ55Jc3d0PLX9MAAAAAAAAWCBudffjSa5PcmuSnyS5pbvvqaobq+rq9WWfSPKcJN+oqn+vqhNP83YAAAAAAABwzhZ5LGG6+2SSk5vOfXjD6zcveS4AAAAAAAB4ikUeSwgAAAAAAAB7grgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwxkJxq6quqKr7qup0Vd2wxfVnVtXX16//qKpesuxBAQAAAAAAYNu4VVUXJLkpyZVJLktybVVdtmnZdUke6e6XJvlkko8te1AAAAAAAABY5M6ty5Oc7u77u/uxJDcnObppzdEkX1p//c0kb6qqWt6YAAAAAAAAsFjcuijJAxuOV9fPbbmmux9P8miSFy5jQAAAAAAAAHjCgQXWbHUHVp/DmlTVsSTH1g9/V1U/XuD7A/AHB5P84nwPATCMvRPg7Ng3Ac6evRPg7P3ZuX7hInFrNcnFG44PJ3nwadasVtWBJM9P8qvNb9Tdx5McT5KqOtXdK+cyNMB+Ze8EOHv2ToCzY98EOHv2ToCzV1WnzvVrF3ks4R1JjlTVpVV1YZJrkpzYtOZEknesv35Lku9191Pu3AIAAAAAAIA/xrZ3bnX341V1fZJbk1yQ5AvdfU9V3ZjkVHefSPL5JF+pqtNZu2Prmp0cGgAAAAAAgP1pkccSprtPJjm56dyHN7z+bZK/Ocvvffws1wNg7wQ4F/ZOgLNj3wQ4e/ZOgLN3zntneXogAAAAAAAAUyzymVsAAAAAAACwJ+x43KqqK6rqvqo6XVU3bHH9mVX19fXrP6qql+z0TAB73QJ753ur6t6quruqvltVLz4fcwLsFdvtmxvWvaWquqpWdnM+gL1okb2zqt66/nPnPVX11d2eEWCvWeD39Uuq6raqumv9d/arzsecAHtFVX2hqh6qqh8/zfWqqk+t76t3V9VrFnnfHY1bVXVBkpuSXJnksiTXVtVlm5Zdl+SR7n5pkk8m+dhOzgSw1y24d96VZKW7X5Xkm0k+vrtTAuwdC+6bqarnJvm7JD/a3QkB9p5F9s6qOpLkA0le191/nuTvd31QgD1kwZ87P5Tklu5+dZJrknx6d6cE2HO+mOSKM1y/MsmR9X/HknxmkTfd6Tu3Lk9yurvv7+7Hktyc5OimNUeTfGn99TeTvKmqaofnAtjLtt07u/u27v7N+uHtSQ7v8owAe8kiP3MmyT9m7Y8BfrubwwHsUYvsne9OclN3P5Ik3f3QLs8IsNcssnd2kuetv35+kgd3cT6APae7v5/kV2dYcjTJl3vN7UleUFUv2u59dzpuXZTkgQ3Hq+vntlzT3Y8neTTJC3d4LoC9bJG9c6PrknxnRycC2Nu23Ter6tVJLu7uf9nNwQD2sEV+5nxZkpdV1Q+q6vaqOtNf3ALsB4vsnR9N8raqWk1yMsl7dmc0gLHO9v9CkyQHdmycNVvdgdXnsAZgP1l4X6yqtyVZSfL6HZ0IYG87475ZVX+Stcdfv3O3BgIYYJGfOQ9k7fEwb8jakwL+rape2d2/3uHZAPaqRfbOa5N8sbv/qar+MslX1vfO/9n58QBGOqdGtNN3bq0muXjD8eE89Vbc/1tTVQeydrvumW5RA/j/bpG9M1X15iQfTHJ1d/9ul2YD2Iu22zefm+SVSf61qv4zyWuTnKiqlV2bEGDvWfT39W939++7+6dJ7sta7ALYrxbZO69LckuSdPcPkzwrycFdmQ5gpoX+L3SznY5bdyQ5UlWXVtWFWfsQxROb1pxI8o71129J8r3uducWsJ9tu3euP17rc1kLWz77ANjvzrhvdvej3X2wu1/S3S/J2mcVXt3dp87PuAB7wiK/rwhzi4MAACAASURBVH8ryRuTpKoOZu0xhffv6pQAe8sie+fPkrwpSarqFVmLWw/v6pQAs5xI8vZa89okj3b3z7f7oh19LGF3P15V1ye5NckFSb7Q3fdU1Y1JTnX3iSSfz9rtuaezdsfWNTs5E8Bet+De+Ykkz0nyjapKkp9199XnbWiA82jBfROADRbcO29N8tdVdW+S/07y/u7+5fmbGuD8WnDvfF+Sf66qf8jaY7Xe6Q/5gf2sqr6WtcdcH1z/PMKPJHlGknT3Z7P2+YRXJTmd5DdJ3rXQ+9pbAQAAAAAAmGKnH0sIAAAAAAAASyNuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMMa2cauqvlBVD1XVj5/melXVp6rqdFXdXVWvWf6YAAAAAAAAsNidW19McsUZrl+Z5Mj6v2NJPvPHjwUAAAAAAABPtW3c6u7vJ/nVGZYcTfLlXnN7khdU1YuWNSAAAAAAAAA8YRmfuXVRkgc2HK+unwMAAAAAAIClOrCE96gtzvWWC6uOZe3RhXn2s5/9Fy9/+cuX8O0BAAAAAACY5M477/xFdx86l69dRtxaTXLxhuPDSR7camF3H09yPElWVlb61KlTS/j2AAAAAAAATFJV/3WuX7uMxxKeSPL2WvPaJI9298+X8L4AAAAAAADwJNveuVVVX0vyhiQHq2o1yUeSPCNJuvuzSU4muSrJ6SS/SfKunRoWAAAAAACA/W3buNXd125zvZP87dImAgAAAAAAgKexjMcSAgAAAAAAwK4QtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhjobhVVVdU1X1Vdbqqbtji+iVVdVtV3VVVd1fVVcsfFQAAAAAAgP1u27hVVRckuSnJlUkuS3JtVV22admHktzS3a9Ock2STy97UAAAAAAAAFjkzq3Lk5zu7vu7+7EkNyc5umlNJ3ne+uvnJ3lweSMCAAAAAADAmkXi1kVJHthwvLp+bqOPJnlbVa0mOZnkPVu9UVUdq6pTVXXq4YcfPodxAQAAAAAA2M8WiVu1xbnedHxtki929+EkVyX5SlU95b27+3h3r3T3yqFDh85+WgAAAAAAAPa1ReLWapKLNxwfzlMfO3hdkluSpLt/mORZSQ4uY0AAAAAAAAB4wiJx644kR6rq0qq6MMk1SU5sWvOzJG9Kkqp6RdbilucOAgAAAAAAsFTbxq3ufjzJ9UluTfKTJLd09z1VdWNVXb2+7H1J3l1V/5Hka0ne2d2bH10IAAAAAAAAf5QDiyzq7pNJTm469+ENr+9N8rrljgYAAAAAAABPtshjCQEAAAAAAGBPELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGWChuVdUVVXVfVZ2uqhueZs1bq+reqrqnqr663DEBAAAAAAAgObDdgqq6IMlNSf4qyWqSO6rqRHffu2HNkSQfSPK67n6kqv50pwYGAAAAAABg/1rkzq3Lk5zu7vu7+7EkNyc5umnNu5Pc1N2PJEl3P7TcMQEAAAAAAGCxuHVRkgc2HK+un9voZUleVlU/qKrbq+qKZQ0IAAAAAAAAT9j2sYRJaotzvcX7HEnyhiSHk/xbVb2yu3/9pDeqOpbkWJJccsklZz0sAAAAAAAA+9sid26tJrl4w/HhJA9usebb3f377v5pkvuyFruepLuPd/dKd68cOnToXGcGAAAAAABgn1okbt2R5EhVXVpVFya5JsmJTWu+leSNSVJVB7P2mML7lzkoAAAAAAAAbBu3uvvxJNcnuTXJT5Lc0t33VNWNVXX1+rJbk/yyqu5NcluS93f3L3dqaAAAAAAAAPan6t788Vm7Y2VlpU+dOnVevjcAAAAAAADnT1Xd2d0r5/K1izyWEAAAAAAAAPYEcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAPhf9u4vxLO7vOP45zFrFDQqmL2QbDQBY3UrQuwQUrzQopQkF9mbtCQQ/EMwN43SKoUUJUq8UimCEP+kGKJCjTEX7SIrudCIpTQhK2mDiQSW1JolQtZqcxM0Tfv0YiZxnEx2zm5mJvN0Xy9Y+J1zvvOb5+rLzL7nnB8AAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMsShuVdVlVfVIVR2rqhtPsu6qquqqWtm+EQEAAAAAAGDVlnGrqs5KckuSy5McTHJNVR3cZN05ST6a5L7tHhIAAAAAAACSZXduXZLkWHc/2t1PJ7kjyaFN1n0myeeS/GYb5wMAAAAAAIDnLIlb5yV5bN3x8bVzz6mqi5Oc393fPdkbVdX1VXW0qo6eOHHilIcFAAAAAADgzLYkbtUm5/q5i1UvS/KFJB/f6o26+9buXunulf379y+fEgAAAAAAALIsbh1Pcv664wNJHl93fE6Styf5YVX9LMmlSQ5X1cp2DQkAAAAAAADJsrh1f5KLqurCqjo7ydVJDj97sbuf7O5zu/uC7r4gyb1JruzuozsyMQAAAAAAAGesLeNWdz+T5IYkdyf5aZI7u/uhqrq5qq7c6QEBAAAAAADgWfuWLOruI0mObDh30wusfc+LHwsAAAAAAACeb8ljCQEAAAAAAGBPELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGWBS3quqyqnqkqo5V1Y2bXP9YVT1cVQ9W1fer6k3bPyoAAAAAAABnui3jVlWdleSWJJcnOZjkmqo6uGHZA0lWuvsdSe5K8rntHhQAAAAAAACW3Ll1SZJj3f1odz+d5I4kh9Yv6O57uvuptcN7kxzY3jEBAAAAAABgWdw6L8lj646Pr517Idcl+d5mF6rq+qo6WlVHT5w4sXxKAAAAAAAAyLK4VZuc600XVl2bZCXJ5ze73t23dvdKd6/s379/+ZQAAAAAAACQZN+CNceTnL/u+ECSxzcuqqr3JflEknd392+3ZzwAAAAAAAD4nSV3bt2f5KKqurCqzk5ydZLD6xdU1cVJvprkyu5+YvvHBAAAAAAAgAVxq7ufSXJDkruT/DTJnd39UFXdXFVXri37fJJXJ/lOVf1rVR1+gbcDAAAAAACA07bksYTp7iNJjmw4d9O61+/b5rkAAAAAAADgeZY8lhAAAAAAAAD2BHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgjEVxq6ouq6pHqupYVd24yfVXVNW3167fV1UXbPegAAAAAAAAsGXcqqqzktyS5PIkB5NcU1UHNyy7Lsmvu/vNSb6Q5LPbPSgAAAAAAAAsuXPrkiTHuvvR7n46yR1JDm1YcyjJ19de35XkvVVV2zcmAAAAAAAALItb5yV5bN3x8bVzm67p7meSPJnk9dsxIAAAAAAAADxr34I1m92B1aexJlV1fZLr1w5/W1U/WfD9Afidc5P88qUeAmAYeyfAqbFvApw6eyfAqfuD0/3CJXHreJLz1x0fSPL4C6w5XlX7krw2ya82vlF335rk1iSpqqPdvXI6QwOcqeydAKfO3glwauybAKfO3glw6qrq6Ol+7ZLHEt6f5KKqurCqzk5ydZLDG9YcTvKBtddXJflBdz/vzi0AAAAAAAB4Mba8c6u7n6mqG5LcneSsJLd190NVdXOSo919OMnXknyzqo5l9Y6tq3dyaAAAAAAAAM5MSx5LmO4+kuTIhnM3rXv9myR/dorf+9ZTXA+AvRPgdNg7AU6NfRPg1Nk7AU7dae+d5emBAAAAAAAATLHkM7cAAAAAAABgT9jxuFVVl1XVI1V1rKpu3OT6K6rq22vX76uqC3Z6JoC9bsHe+bGqeriqHqyq71fVm16KOQH2iq32zXXrrqqqrqqV3ZwPYC9asndW1Z+v/dz5UFX9/W7PCLDXLPh9/Y1VdU9VPbD2O/sVL8WcAHtFVd1WVU9U1U9e4HpV1RfX9tUHq+qdS953R+NWVZ2V5JYklyc5mOSaqjq4Ydl1SX7d3W9O8oUkn93JmQD2uoV75wNJVrr7HUnuSvK53Z0SYO9YuG+mqs5J8tEk9+3uhAB7z5K9s6ouSvI3Sd7V3X+Y5C93fVCAPWThz52fTHJnd1+c5OokX9rdKQH2nNuTXHaS65cnuWjt3/VJvrzkTXf6zq1Lkhzr7ke7++kkdyQ5tGHNoSRfX3t9V5L3VlXt8FwAe9mWe2d339PdT60d3pvkwC7PCLCXLPmZM0k+k9U/BvjNbg4HsEct2Ts/nOSW7v51knT3E7s8I8Bes2Tv7CSvWXv9DRrFRgAAIABJREFU2iSP7+J8AHtOd/8oya9OsuRQkm/0qnuTvK6q3rDV++503DovyWPrjo+vndt0TXc/k+TJJK/f4bkA9rIle+d61yX53o5OBLC3bblvVtXFSc7v7u/u5mAAe9iSnznfkuQtVfXPVXVvVZ3sL24BzgRL9s5PJ7m2qo4nOZLkI7szGsBYp/p/oUmSfTs2zqrN7sDq01gDcCZZvC9W1bVJVpK8e0cnAtjbTrpvVtXLsvr46w/u1kAAAyz5mXNfVh8P856sPingn6rq7d39Xzs8G8BetWTvvCbJ7d39t1X1x0m+ubZ3/u/Ojwcw0mk1op2+c+t4kvPXHR/I82/FfW5NVe3L6u26J7tFDeD/uyV7Z6rqfUk+keTK7v7tLs0GsBdttW+ek+TtSX5YVT9LcmmSw1W1smsTAuw9S39f/8fu/u/u/vckj2Q1dgGcqZbsndcluTNJuvtfkrwyybm7Mh3ATIv+L3SjnY5b9ye5qKourKqzs/ohioc3rDmc5ANrr69K8oPuducWcCbbcu9ce7zWV7Matnz2AXCmO+m+2d1Pdve53X1Bd1+Q1c8qvLK7j7404wLsCUt+X/+HJH+SJFV1blYfU/jork4JsLcs2Tt/nuS9SVJVb8tq3Dqxq1MCzHI4yftr1aVJnuzuX2z1RTv6WMLufqaqbkhyd5KzktzW3Q9V1c1Jjnb34SRfy+rtuceyesfW1Ts5E8Bet3Dv/HySVyf5TlUlyc+7+8qXbGiAl9DCfROAdRbunXcn+dOqejjJ/yT56+7+z5duaoCX1sK98+NJ/q6q/iqrj9X6oD/kB85kVfWtrD7m+ty1zyP8VJKXJ0l3fyWrn094RZJjSZ5K8qFF72tvBQAAAAAAYIqdfiwhAAAAAAAAbBtxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMbYMm5V1W1V9URV/eQFrldVfbGqjlXVg1X1zu0fEwAAAAAAAJbduXV7kstOcv3yJBet/bs+yZdf/FgAAAAAAADwfFvGre7+UZJfnWTJoSTf6FX3JnldVb1huwYEAAAAAACAZ+3bhvc4L8lj646Pr537xcaFVXV9Vu/uyqte9ao/eutb37oN3x4AAAAAAIBJfvzjH/+yu/efztduR9yqTc71Zgu7+9YktybJyspKHz16dBu+PQAAAAAAAJNU1X+c7tcu+cytrRxPcv664wNJHt+G9wUAAAAAAIDfsx1x63CS99eqS5M82d3PeyQhAAAAAAAAvFhbPpawqr6V5D1Jzq2q40k+leTlSdLdX0lyJMkVSY4leSrJh3ZqWAAAAAAAAM5sW8at7r5mi+ud5C+2bSIAAAAAAAB4AdvxWEIAAAAAAADYFeIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGIviVlVdVlWPVNWxqrpxk+tvrKp7quqBqnqwqq7Y/lEBAAAAAAA4020Zt6rqrCS3JLk8ycEk11TVwQ3LPpnkzu6+OMnVSb603YMCAAAAAADAkju3LklyrLsf7e6nk9yR5NCGNZ3kNWuvX5vk8e0bEQAAAAAAAFYtiVvnJXls3fHxtXPrfTrJtVV1PMmRJB/Z7I2q6vqqOlpVR0+cOHEa4wIAAAAAAHAmWxK3apNzveH4miS3d/eBJFck+WZVPe+9u/vW7l7p7pX9+/ef+rQAAAAAAACc0ZbEreNJzl93fCDPf+zgdUnuTJLu/pckr0xy7nYMCAAAAAAAAM9aErfuT3JRVV1YVWcnuTrJ4Q1rfp7kvUlSVW/Latzy3EEAAAAAAAC21ZZxq7ufSXJDkruT/DTJnd39UFXdXFVXri37eJIPV9W/JflWkg9298ZHFwIAAAAAAMCLsm/Jou4+kuTIhnM3rXv9cJJ3be9oAAAAAAAA8PuWPJYQAAAAAAAA9gRxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwCA/2vvjkI0u8s7jv8edxt7YWrB3ULJJibgpnQbhJQhWLxQ0ZbEi92bVBKQagnmKpVWESIWW+KVShEKaWuKohVsjLnQRSK50EhFmpAtFjGRwBLFLCkk1Zgb0Zj26cVM4jiZ7JxZZ2bn6Xw+sPCec/7zznP1Z2a/c84LAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjLEoblXV9VX1WFWdrarbX2bNO6rq0ap6pKo+v7NjAgAAAAAAQHJ4qwVVdSjJnUn+OMm5JA9X1enufnTdmuNJPpjkjd39TFX9zm4NDAAAAAAAwMG15M6t65Kc7e7Hu/u5JHcnObVhzXuS3NndzyRJdz+1s2MCAAAAAADAsrh1WZIn1h2fWzu33tVJrq6qb1XVg1V1/U4NCAAAAAAAAC/Y8rGESWqTc73J+xxP8uYkx5J8s6qu6e6f/MobVd2a5NYkueKKK7Y9LAAAAAAAAAfbkju3ziW5fN3xsSRPbrLmy939i+7+fpLHshq7fkV339XdK929cvTo0QudGQAAAAAAgANqSdx6OMnxqrqqqi5JclOS0xvWfCnJW5Kkqo5k9TGFj+/koAAAAAAAALBl3Oru55PcluT+JN9Lck93P1JVd1TVybVl9yf5UVU9muSBJB/o7h/t1tAAAAAAAAAcTNW98eOz9sbKykqfOXPmonxvAAAAAAAALp6q+o/uXrmQr13yWEIAAAAAAADYF8QtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxlgUt6rq+qp6rKrOVtXt51l3Y1V1Va3s3IgAAAAAAACwasu4VVWHktyZ5IYkJ5LcXFUnNll3aZL3Jnlop4cEAAAAAACAZNmdW9clOdvdj3f3c0nuTnJqk3UfSfKxJD/bwfkAAAAAAADgRUvi1mVJnlh3fG7t3Iuq6tokl3f3V873RlV1a1WdqaozTz/99LaHBQAAAAAA4GBbErdqk3P94sWqVyT5RJL3b/VG3X1Xd69098rRo0eXTwkAAAAAAABZFrfOJbl83fGxJE+uO740yTVJvlFVP0jyhiSnq2plp4YEAAAAAACAZFncejjJ8aq6qqouSXJTktMvXOzuZ7v7SHdf2d1XJnkwycnuPrMrEwMAAAAAAHBgbRm3uvv5JLcluT/J95Lc092PVNUdVXVytwcEAAAAAACAFxxesqi770ty34ZzH36ZtW/+9ccCAAAAAACAl1ryWEIAAAAAAADYF8QtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMYQtwAAAAAAABhD3AIAAAAAAGAMcQsAAAAAAIAxxC0AAAAAAADGELcAAAAAAAAYQ9wCAAAAAABgDHELAAAAAACAMRbFraq6vqoeq6qzVXX7JtffV1WPVtV3quprVfXanR8VAAAAAACAg27LuFVVh5LcmeSGJCeS3FxVJzYs+3aSle5+fZJ7k3xspwcFAAAAAACAJXduXZfkbHc/3t3PJbk7yan1C7r7ge7+6drhg0mO7eyYAAAAAAAAsCxuXZbkiXXH59bOvZxbknx1swtVdWtVnamqM08//fTyKQEAAAAAACDL4lZtcq43XVj1ziQrST6+2fXuvqu7V7p75ejRo8unBAAAAAAAgCSHF6w5l+TydcfHkjy5cVFVvS3Jh5K8qbt/vjPjAQAAAAAAwC8tuXPr4STHq+qqqrokyU1JTq9fUFXXJvlkkpPd/dTOjwkAAAAAAAAL4lZ3P5/ktiT3J/leknu6+5GquqOqTq4t+3iSVyX5YlX9Z1Wdfpm3AwAAAAAAgAu25LGE6e77kty34dyH171+2w7PBQAAAAAAAC+x5LGEAAAAAAAAsC+IWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGOIWwAAAAAAAIwhbgEAAAAAADCGuAUAAAAAAMAY4hYAAAAAAABjiFsAAAAAAACMIW4BAAAAAAAwhrgFAAAAAADAGOIWAAAAAAAAY4hbAAAAAAAAjCFuAQAAAAAAMIa4BQAAAAAAwBjiFgAAAAAAAGMsiltVdX1VPVZVZ6vq9k2uv7KqvrB2/aGqunKnBwUAAAAAAIAt41ZVHUpyZ5IbkpxIcnNVndiw7JYkz3T365J8IslHd3pQAAAAAAAAWHLn1nVJznb34939XJK7k5zasOZUks+uvb43yVurqnZuTAAAAAAAAFgWty5L8sS643Nr5zZd093PJ3k2yWt2YkAAAAAAAAB4weEFaza7A6svYE2q6tYkt64d/ryqvrvg+wPwS0eS/PfFHgJgGHsnwPbYNwG2z94JsH2/d6FfuCRunUty+brjY0mefJk156rqcJJXJ/nxxjfq7ruS3JUkVXWmu1cuZGiAg8reCbB99k6A7bFvAmyfvRNg+6rqzIV+7ZLHEj6c5HhVXVVVlyS5KcnpDWtOJ3nX2usbk3y9u19y5xYAAAAAAAD8Ora8c6u7n6+q25Lcn+RQkk939yNVdUeSM919Osmnknyuqs5m9Y6tm3ZzaAAAAABJKbdsAAAElUlEQVQAAA6mJY8lTHffl+S+Dec+vO71z5L86Ta/913bXA+AvRPgQtg7AbbHvgmwffZOgO274L2zPD0QAAAAAACAKZZ85hYAAAAAAADsC7set6rq+qp6rKrOVtXtm1x/ZVV9Ye36Q1V15W7PBLDfLdg731dVj1bVd6rqa1X12osxJ8B+sdW+uW7djVXVVbWyl/MB7EdL9s6qesfaz52PVNXn93pGgP1mwe/rV1TVA1X17bXf2d9+MeYE2C+q6tNV9VRVffdlrldV/f3avvqdqvrDJe+7q3Grqg4luTPJDUlOJLm5qk5sWHZLkme6+3VJPpHko7s5E8B+t3Dv/HaSle5+fZJ7k3xsb6cE2D8W7pupqkuTvDfJQ3s7IcD+s2TvrKrjST6Y5I3d/QdJ/nLPBwXYRxb+3PnXSe7p7muT3JTkH/Z2SoB95zNJrj/P9RuSHF/7d2uSf1zyprt959Z1Sc529+Pd/VySu5Oc2rDmVJLPrr2+N8lbq6p2eS6A/WzLvbO7H+jun64dPpjk2B7PCLCfLPmZM0k+ktU/BvjZXg4HsE8t2Tvfk+TO7n4mSbr7qT2eEWC/WbJ3dpLfWnv96iRP7uF8APtOd/9bkh+fZ8mpJP/Sqx5M8ttV9btbve9ux63Lkjyx7vjc2rlN13T380meTfKaXZ4LYD9bsneud0uSr+7qRAD725b7ZlVdm+Ty7v7KXg4GsI8t+Znz6iRXV9W3qurBqjrfX9wCHARL9s6/TfLOqjqX5L4kf7E3owGMtd3/C02SHN61cVZtdgdWX8AagINk8b5YVe9MspLkTbs6EcD+dt59s6pekdXHX797rwYCGGDJz5yHs/p4mDdn9UkB36yqa7r7J7s8G8B+tWTvvDnJZ7r776rqj5J8bm3v/N/dHw9gpAtqRLt959a5JJevOz6Wl96K++Kaqjqc1dt1z3eLGsD/d0v2zlTV25J8KMnJ7v75Hs0GsB9ttW9emuSaJN+oqh8keUOS01W1smcTAuw/S39f/3J3/6K7v5/ksazGLoCDasneeUuSe5Kku/89yW8mObIn0wHMtOj/Qjfa7bj1cJLjVXVVVV2S1Q9RPL1hzekk71p7fWOSr3e3O7eAg2zLvXPt8VqfzGrY8tkHwEF33n2zu5/t7iPdfWV3X5nVzyo82d1nLs64APvCkt/Xv5TkLUlSVUey+pjCx/d0SoD9Zcne+cMkb02Sqvr9rMatp/d0SoBZTif5s1r1hiTPdvd/bfVFu/pYwu5+vqpuS3J/kkNJPt3dj1TVHUnOdPfpJJ/K6u25Z7N6x9ZNuzkTwH63cO/8eJJXJfliVSXJD7v75EUbGuAiWrhvArDOwr3z/iR/UlWPJvmfJB/o7h9dvKkBLq6Fe+f7k/xzVf1VVh+r9W5/yA8cZFX1r1l9zPWRtc8j/Jskv5Ek3f1PWf18wrcnOZvkp0n+fNH72lsBAAAAAACYYrcfSwgAAAAAAAA7RtwCAAAAAABgDHELAAAAAACAMcQtAAAAAAAAxhC3AAAAAAAAGEPcAgAAAAAAYAxxCwAAAAAAgDHELQAAAAAAAMb4P7Yh8E/b2YoFAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(4, 1, sharey='row',figsize=(30,15*3))\n", "i = 0\n", "th = [0,1,2,3]\n", "\n", "for i in range(4):\n", " ax[i].set_title('ECG at threshold = ' + str(thresholdComplete[th[i]]), fontsize=35)\n", " ax[i].set_xlabel(\"Time\", fontsize=35)\n", " ax[i].set_ylabel(\"ADC value\", fontsize=35)\n", " ax[i].tick_params(labelsize = 35)\n", "\n", " t = data[th[i]][0][0].values\n", " v = data[th[i]][0][1].values\n", " startI,stopI = getIdxs(t,start,stop)\n", " #print(selectedFileNames[i])\n", " ax[i].plot(t[startI:stopI]/360,v[startI:stopI])\n", " ax[i].grid()\n", "\n", " #ax[i].legend(\"bablebba\", fontsize=35)\n", "#fig.savefig('/home/zanoli/Desktop/figSaved/eventSampled.pdf',format='pdf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## gQRS Method\n", "Let's now see the gQRS method implementation:" ] }, { "cell_type": "markdown", "metadata": { "heading_collapsed": true }, "source": [ "### Basics of the standard gQRS detector:" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "The standard gQRS is composed by 4 main phases. during these phases a set of parameter is constantly re-adapting to find the (estimated and local) optimal values. The 4 phases can be summarized as:\n", "\n", "1. Smoothing and amplification\n", "2. \"Matched\" filter\n", "3. Peak detection and saving\n", "4. Peak lookback\n", "\n", "Each of this phases are important to understand in order to be able to re-formulate this algorithm as we like.\n", "\n", "Before going into details, is usefull to discuss about some meaningfull parameter that the algorithm use to detect the QRS complex (Some learned, some fixed). All of them are described in number of samples instaead of a physical value, this is not a problem: the inital values, given by the user, are physical values that get transformed in sample numbers by the algorithm using the information about the ADC used :\n", "\n", "* rrmean = Average RR interval\n", "* rrdev = Standard deviation of the RR interval\n", "* rrmin = Minimum RR interval\n", "* rrmax = Maximum RR interval\n", "* rrinc = Increment step of the estimated RR interval (for learning)\n", "* dt = One fourth of the standard duration of the QRS complex, other 3 variable like this are present (1, 1/2, 1/3 of the duration of the QRS complex)\n", "* pthr = Threshold for a cuspid to be considered peak\n", "* qthr = Threshold for a peak to be considered a candidate for QRS\n", "* pthmin = Lower possible value for pthr\n", "* qthmin = Lower possible value for qthr\n", "* tamean = Average amplitude of a T-Wave\n" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "#### Smoothing and amplification\n", "During this phase the signal pass throught a trapezoidal (R-C) filter and get amplified by the fixed delta time that the algorithm suppose to be a standard QRS duration. This filter mainly try to remove the noise present" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "#### \"Matched\" filter\n", "Even if it's not immidate to get from the code, the \"Matched\" filter is a FIR filter that operate upon the integrated smothed signal, it's shape is customly designed to obtain a strong output if a QRS complex is present and almost zero otherwise (except for strong T-Waves). We wont get any further in this section mainly because this filter is one of the key to the good behaviour of this algorithm and was one of the main blocks that needed to be re-formulated. We will see it in detail in the next section" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "#### Peack detection and saving\n", "The peak detetion is a standard check that control the left and right (previous and next) point of the filtered signal and check if it's higher than both and higher than pthr.\n", "\n", "Every point that pass this condition is saved in a circular structure that will be analyzed duing the lookback phase: in this algorithm the detection of the QRS complex is not immidatate but happen ater a variable (but bounded) time delay that we will see in the next phase" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "#### Peak lookback\n", "Whenever a peak is detected the algorithm look inside the cicular peak buffer analyzing all the peaks that fall between the last QRS complex found untill the time of the present peak minus the standard time between QRS and T wave (this exclude the present detected peak from the analysis and the previous QRS complex if we detected a t-wave).\n", "\n", "If the peak analyzed is higher than qthr, is the dominant peak in his dierct neighbourhood and happened after rrmin from the last QRS complex then the peak is considered to be the next QRS complex, in this case the algorithm also check the next peak in the buffer to search for a possible T-Wave. All the parameters (estimated intervals and thresholds) get updated in a classical way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Customization of the gQRS\n", "Let's now discuss on the problematique of this algorithm and how they where solved in order to adapt it to our enviroment.\n", "\n", "As discussed in the indroduction of this notebook, the idea of this study is to find a method to apply state-of-the-art algorithm to non linarly sampled signals. In particular, the prior that the subsampled signal (from now on, the signal) can be expressed as a sequence of pice-wise linear segment between each point without exceding a fixed error, will be extremly helpful." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Smoothing filter\n", "Because of the method used for the linear sampling, no smoothing was required. The only operation performed was just the amplification by $4\\delta$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The problem with non linear sampling\n", "The original gQRS algorithm relay on a subsequent summation of the output of the matched filter:\n", "\n", "$y[t] = \\sum_{i=4\\delta}^{t} \\sum_{j=-4}^{4}h[j]x[i-j\\delta]$\n", "\n", "where $\\delta = \\frac{QRS duration}{4}$ is the delay of the FIR filter and $h[j]$ is the $j^{th}$ tap of the filter. The Impulsive response of the filter is the following:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T13:43:20.094571Z", "start_time": "2019-08-18T13:43:19.980741Z" } }, "outputs": [], "source": [ "taps = [1,-2,-4,10,0,-10,4,2,-1]\n", "plt.plot(taps)\n", "f = np.fft.fftshift(np.fft.fft(taps)) \n", "fig = plt.gcf()\n", "fig.set_size_inches(10.5, 6.5)\n", "plt.grid()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The filter composed of 9 taps: each of them dealyed $\\delta$ from the previous one, covering a time-span of $8\\delta$ around the point under examination" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This require to know all the points of the original signal: to know the value of the filter at $t$ we mustknow the output of it at $t-1$. Is possible to show that, if we simply apply the filter to each of point of the sub-sampled signal (lineraly interpolating the 8 point around the interesting one) the results, in average, diverge.\n", "\n", "One solution could be to simply interpolate linearly all the points between 2 sample and apply the algorithm to the signal obtained in this way. Instaed we want to have the fewer possible operations. We will search, hence, to find a method to use only the needed points (so, 8 points for each new sample)\n", "\n", "In order to solve this we can re-write our filter as:\n", "\n", "$y[t] = \\sum_{j=-4}^{4} h[j] \\sum_{i=4\\delta}^{t} x[i-j\\delta]$\n", "\n", "Now we define $s_{j}[t] = \\sum_{i=4\\delta}^{t} x[i-j\\delta]$\n", "\n", "If we assume to start from $-8\\delta$ with a signal filled with 0, then we can wrtie:\n", "\n", "$s_{j}[t] = \\sum_{i=0}^{t+j\\delta} x[i]$\n", "\n", "Wich is the numerical integration of the signal from 0 to $t+j\\delta$\n", "\n", "Hence: $y[t] = \\sum_{j=-4}^{4} h[j] s_{j}[t]$\n", "\n", "From here we can see that the result of the filtering is a FIR filter applied to the integrated signal.\n", "\n", "The idea from here is to accumulate and save the integral of the signal at each new point that arrive and apply the filter to this saved signal in order to avoid the need to accumulate the output of the filter. This will require to \"interpolate\" only 8 points from the integrated signal. Of course we will not be able to perform a linear interpolation of it and neither to find an interpolator formula as we will show shortly. Let's write the generic numerical integration of a digital signal:\n", "\n", "$i[t] = \\sum_{i=0}^{t-1} x[i]$\n", "\n", "(Note: The sum need to reach $t-1$ otherwise $i[t]+i[t:t+k] \\neq i[t+k]$)\n", "\n", "we can now use the prior that our signal is pice-wise linear and find a formula that give us the generic integral value between 2 given points on a linear signal:\n", "\n", "$i[a:b] = \\sum_{i=a}^{b-1} m*i+q = \\frac{m}{2}(b(b+1)-a(a+1))+q(b-a)$\n", "\n", "Now we have a vector of $i[t]$ where the $t_{s}$ are not equaly separeted. Now we only need to find a method to obtain the 8 points (spaced by $\\delta$) around any point of our signal\n", "\n", "Is it possible, given $i[t_{1}]$ and $i[t_{2}]$ to find $i[t_{k}]$ with $t_{1} \n", " # [-4dt,-3dt,-2dt,-1dt,0,+1dt,+2dt,+3dt,+4dt]\n", " #print(pts)\n", " \n", " #FILTER PART (qrst shape)\n", " dv1 = int(pts[-4] - pts[3])#dt1\n", " dv = dv1 << 1\n", " dv -= int(pts[-3]- pts[2])#dt2\n", " dv = dv << 1\n", " dv += dv1\n", " dv -= int(pts[-2]- pts[1])#dt3\n", " dv = dv << 1\n", " dv += int(pts[-1]- pts[0])#dt4\n", " \n", " \n", " self.v1 = dv\n", " \n", " v0 = int(self.v1 / self.c.v1norm)\n", " self.OUTFILTER.append(v0**2)\n", " self.TIMES.append(self.data.t[self.i])\n", " self.qfv_put(self.j, v0 * v0) \n", " \n", " def add_peak(self, peak_time, peak_amp, type):\n", " p = self.current_peak.next_peak\n", " p.time = peak_time\n", " p.amp = peak_amp\n", " p.type = type\n", " self.current_peak = p\n", " p.next_peak.amp = 0 \n", " def peaktype(self, p):\n", " # peaktype() returns 1 if p is the most prominent peak in its neighborhood, 2\n", " # otherwise. The neighborhood consists of all other peaks within rrmin.\n", " # Normally, \"most prominent\" is equivalent to \"largest in amplitude\", but this\n", " # is not always true. For example, consider three consecutive peaks a, b, c\n", " # such that a and b share a neighborhood, b and c share a neighborhood, but a\n", " # and c do not; and suppose that amp(a) > amp(b) > amp(c). In this case, if\n", " # there are no other peaks, a is the most prominent peak in the (a, b)\n", " # neighborhood. Since b is thus identified as a non-prominent peak, c becomes\n", " # the most prominent peak in the (b, c) neighborhood. This is necessary to\n", " # permit detection of low-amplitude beats that closely precede or follow beats\n", " # with large secondary peaks (as, for example, in R-on-T PVCs).\n", " if p.type:\n", " return p.type\n", " else:\n", " a = p.amp\n", " t0 = p.time - self.c.rrmin\n", " t1 = p.time + self.c.rrmin\n", "\n", " if t0 < 0:\n", " t0 = 0\n", "\n", " pp = p.prev_peak\n", " while t0 < pp.time and pp.time < pp.next_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and self.peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.prev_peak\n", "\n", " pp = p.next_peak\n", " while pp.time < t1 and pp.time > pp.prev_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and self.peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.next_peak\n", "\n", " p.type = 1\n", " return p.type\n", " def find_missing(self, r, p):\n", " if r is None or p is None:\n", " return None\n", "\n", " minrrerr = p.time - r.time\n", "\n", " s = None\n", " q = r.next_peak\n", " while q.time < p.time:\n", " if self.peaktype(q) == 1:\n", " rrtmp = q.time - r.time\n", " rrerr = rrtmp - self.c.rrmean\n", " if rrerr < 0:\n", " rrerr = -rrerr\n", " if rrerr < minrrerr:\n", " minrrerr = rrerr\n", " s = q\n", " # end:\n", " q = q.next_peak\n", "\n", " return s\n", " \n", " class Conf(object):\n", " \"\"\"\n", " Initial signal configuration object for this qrs detector\n", " \"\"\"\n", " #adc gain is the lvl/mV, can be found in the .hea file, default is 200\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " def __init__(self, fs = 360, adc_gain = 200, hr=75,\n", " RRdelta=0.2, RRmin=0.28, RRmax=2.4,\n", " QS=0.07, QT=0.35,\n", " RTmin=0.25, RTmax=0.33,\n", " QRSa=750, QRSamin=150,\n", " thresh=1.0):\n", " self.fs = fs\n", "\n", " self.sps = int(1 * fs + 0.5)\n", " self.spm = int(60* fs + 0.5)\n", "\n", " self.hr = hr\n", " self.RR = 60.0 / self.hr\n", " self.RRdelta = RRdelta\n", " self.RRmin = RRmin\n", " self.RRmax = RRmax\n", " self.QS = QS\n", " self.QT = QT\n", " self.RTmin = RTmin\n", " self.RTmax = RTmax\n", " self.QRSa = QRSa\n", " self.QRSamin = QRSamin\n", " self.thresh = thresh\n", "\n", " self._NORMAL = 1 # normal beat\n", " self._ARFCT = 16 # isolated QRS-like artifact\n", " self._NOTE = 22 # comment annotation\n", " self._TWAVE = 27 # T-wave peak\n", " self._NPEAKS = 10 # number of peaks buffered (per signal)\n", " self._BUFLN = 32768 # must be a power of 2, see qf()\n", "\n", " self.rrmean = int(self.RR * self.sps)\n", " self.rrdev = int(self.RRdelta * self.sps)\n", " self.rrmin = int(self.RRmin * self.sps)\n", " self.rrmax = int(self.RRmax * self.sps)\n", "\n", " self.rrinc = int(self.rrmean / 40)\n", " if self.rrinc < 1:\n", " self.rrinc = 1\n", "\n", " self.dt = int(self.QS * self.sps / 4)\n", " if self.dt < 1:\n", " raise Exception('Sampling rate is too low. Unable to use signal.')\n", "\n", " self.rtmin = int(self.RTmin * self.sps)\n", " self.rtmean = int(0.75 * self.QT * self.sps)\n", " self.rtmax = int(self.RTmax * self.sps)\n", " \n", " #dv: min peak to peak val for QRS complex in uVolts\n", " #QRSamin: min peak to peak val for QRS complex in samples\n", " dv = adc_gain * self.QRSamin * 0.001\n", " self.pthr = int((self.thresh * dv * dv) / 6)\n", " self.qthr = self.pthr << 1\n", " self.pthmin = self.pthr >> 2\n", " self.qthmin = int((self.pthmin << 2) / 3)\n", " self.tamean = self.qthr # initial value for mean T-wave amplitude\n", "\n", " # Filter constants and thresholds.\n", " self.dt2 = 2 * self.dt\n", " self.dt3 = 3 * self.dt\n", " self.dt4 = 4 * self.dt\n", "\n", " self.smdt = self.dt\n", " self.v1norm = self.smdt * self.dt * 64\n", "\n", " self.smt = 0\n", " self.smt0 = 0 + self.smdt \n", " class Peak(object):\n", " def __init__(self, peak_time, peak_amp, peak_type):\n", " self.time = peak_time\n", " self.amp = peak_amp\n", " self.type = peak_type\n", " self.next_peak = None\n", " self.prev_peak = None\n", " def __str__(self):\n", " return \"Peak object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nAmp: \"+str(self.amp)+\\\n", " \"\\nTyppe: \"+str(self.type)\n", " class Annotation(object):\n", " def __init__(self, ann_time, ann_type, ann_subtype, ann_num):\n", " self.time = ann_time\n", " self.type = ann_type\n", " self.subtype = ann_subtype\n", " self.num = ann_num\n", " def __str__(self):\n", " return \"Annotation object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nNum: \"+str(self.num)+\\\n", " \"\\nTyppe: \"+str(self.type) \n", " class DataInteg:\n", " def __init__(self,length):\n", " self.i = collections.deque( maxlen=length)\n", " self.m = collections.deque( maxlen=length)\n", " self.q = collections.deque( maxlen=length)\n", " self.t = collections.deque( maxlen=length)\n", " self.allInteg = []\n", " self.lastVal = None\n", " self.idx = 0\n", " self.len = length\n", " \n", " def put(self,t,v):\n", " presentVal = [t,v]\n", " if self.lastVal == None:\n", " self.lastVal = presentVal\n", " return 0\n", " #print(self.lastVal)\n", " #print(presentVal)\n", " integ,m,q = self.computeIMQ(self.lastVal,presentVal)\n", " self.lastVal = presentVal\n", " #print(integ)\n", " \n", " \n", " if len(self.i)>0:\n", " newInteg = self.i[len(self.i)-1]+integ\n", " else:\n", " newInteg = integ\n", " #print(t,' --> M: ',m,' Q: ',q, 'Integ: ',integ,' total integ: ',newInteg)\n", " self.allInteg.append(newInteg)\n", " self.i.append(newInteg)\n", " self.m.append(m)\n", " self.q.append(q)\n", " self.t.append(t)\n", " return len(self.i)\n", " \n", " \n", " def computeIMQ(self,pt1,pt2):\n", " m = (pt2[1]-pt1[1])/(pt2[0]-pt1[0])\n", " q = pt1[1]-pt1[0]*m\n", " t0 = pt1[0]\n", " t1 = pt2[0]\n", " #integral version not considering the fist point (for accumulation )\n", " integ = m/2*(t1*(t1+1)-t0*(t0+1)) + q*(t1-t0)\n", " #integ = q*(t1-t0+1)+m/2*(t1**2-t0**2+t1+t0)\n", " #print(t0,' --> M: ',m,' Q: ',q)\n", " return integ,m,q\n", " \n", " def reSampler(self, index, dt):\n", " t = self.t[index]\n", " newTimes = [t-4*dt]+[t-3*dt]+[t-2*dt]+[t-1*dt]+[t]+[t+1*dt]+[t+2*dt]+[t+3*dt]+[t+4*dt]\n", " \n", " newSamples = []\n", " firstIdx = None\n", " lastIdx = None\n", " \n", " for t in newTimes:\n", " idx = 0\n", " #can be optimized for t<0\n", " if t <= self.t[-2]:\n", " present = False\n", " for i in range(len(self.t)-1):\n", " if self.t[i] ==t:\n", " idx = i\n", " present = True\n", " break\n", " elif self.t[i+1] ==t:\n", " idx = i+1\n", " present = True\n", " break\n", " elif self.t[i] <=t and self.t[i+1]>t:\n", " idx = i\n", " break\n", " if present:\n", " new = self.i[idx]\n", " else:\n", " lastInteg = self.i[idx]\n", " t0 = self.t[idx]\n", " m = self.m[idx+1]\n", " q = self.q[idx+1]\n", " new = lastInteg+m/2*(t*(t+1)-t0*(t0+1)) + q*(t-t0)\n", " if firstIdx == None:\n", " firstIdx = idx\n", " if t == newTimes[-1]:\n", " lastIdx = idx\n", " else:\n", " if firstIdx == None:\n", " firstIdx = -1\n", " if t == newTimes[-1]:\n", " lastIdx = -1\n", " new = self.i[-1]\n", " t0 = self.t[idx]\n", " m = self.m[idx+1]\n", " q = self.q[idx+1]\n", " #print(t,' ---> ',t0,' -- ',self.t[idx+1],('M: ',m,'Q: ',q))\n", " newSamples.append(int(new))\n", "# print(firstIdx,lastIdx)\n", "# print(\"Saved times: \",list(self.t)[firstIdx:lastIdx+2])\n", "# print(\"Saved samples: \",list(self.i)[firstIdx:lastIdx+2])\n", "# print(\"Found samples: \",newSamples)\n", "# print(\"(for times t = \",newTimes,\")\")\n", " return newSamples\n", " \n", " \n", " #ENTRY POINT, NEED A CONF OBJ FOR INSIDE THE GQRS CLASS\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " #adc_zero in this case is 0 because we already removed the average\n", " \n", " def __init__(self, adc_zero = 1024):\n", " self.c = gQRS.Conf()\n", " \n", " self.annotations = []\n", " self.sample_valid = False\n", " self.adc_zero = adc_zero\n", "\n", " self.qfv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.smv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.OUTFILTER = []\n", " self.TIMES = []\n", " self.v1 = 0\n", " \n", " #self.data = self.Data(self.c.dt4*4+6)\n", " self.data = self.DataInteg(self.c.dt4*4+6)\n", " \n", "\n", " t0 = 0\n", " self.t = -self.c.dt4\n", " self.i = 0\n", " self.j = 0\n", "\n", " self.annot = gQRS.Annotation(-self.c.rrmin, \"NOTE\", 0, 0)\n", " # Cicular buffer of Peaks\n", " first_peak = gQRS.Peak(0, 0, 0)\n", " tmp = first_peak\n", " for _ in range(1, self.c._NPEAKS):\n", " tmp.next_peak = gQRS.Peak(0, 0, 0)\n", " tmp.next_peak.prev_peak = tmp\n", " tmp = tmp.next_peak\n", " tmp.next_peak = first_peak\n", " first_peak.prev_peak = tmp\n", " self.current_peak = first_peak\n", " \n", "\n", " self.state = \"RUNNING\"\n", "# self.gqrs(t0, self.tf)\n", "# return self.annotations\n", "\n", " #State variable\n", " self.q0 = None\n", " self.q1 = 0\n", " self.q2 = 0\n", " self.rr = None\n", " self.rrd = None\n", " self.rt = None\n", " self.rtd = None\n", " self.rtdmin = None\n", "\n", " self.p = None # (Peak)\n", " self.q = None # (Peak)\n", " self.r = None # (Peak)\n", " self.tw = None # (Peak)\n", "\n", " self.last_peak = 0\n", " self.last_qrs = 0\n", " self.lastdV = 0\n", " def isPeak(self,t,val,ver=False):\n", " #print('-'*40)\n", " #print(t)\n", " # VALUE TO RETURN\n", " # [integral, square, time, reDo, ok] is the classical return\n", " # for gqrs we wil return \n", " # [type, 0, time, False, ok]\n", " # type = dominant\n", " ok = False\n", " noneFlag = False\n", " rVal = 0\n", " noneVal = 0\n", " tVal = 0\n", " toReturn = []\n", " \n", " \n", " val = (val-1024)*25\n", " dataLen = self.data.put(t,val)\n", " if dataLen >= self.data.len//2:\n", " self.i = dataLen//2\n", " #self.qf()\n", " self.qfInteg()\n", " q0 = self.qfv_at(self.j)\n", " q1 = self.qfv_at(self.j - 1)\n", " q2 = self.qfv_at(self.j - 2)\n", " tPres = self.data.t[self.i-1]\n", " if ver:\n", " print(q0, q1, q2)\n", " pass\n", " \n", " if q1 > self.c.pthr and q2 < q1 and q1 >= q0 :\n", " #print(\"T: \", tPres, \"Filter Values: \", q0, q1, q2,\"\\tThreshold: \",self.c.qthr)\n", " tPeak = self.data.t[self.i-1]\n", " self.add_peak(tPeak, q1, 0)\n", " self.last_peak = tPeak\n", " #print(self.current_peak)\n", " p = self.current_peak.next_peak\n", " if ver:\n", " print(\"\\t\\tPEAK\")\n", " print(\"TIME: \",tPeak)\n", " print(\"THRESHOLD: \",self.c.pthr)\n", " print(\"VALUE: \",q1)\n", " print(\"TIME NEXT: \",p.time ,' Tpres-rtMax: ',tPres - self.c.rtmax)\n", " l = 0\n", " while p.time < tPres - self.c.rtmax:\n", " l+=1\n", " if ver:\n", " #print('*-*-'*10)\n", " #print(p.time,' < ' ,tPres, '-', self.c.rtmax)\n", " pass\n", " \n", " \n", "\n", " \n", " #peaks after the last annotation to the current one (dominant)\n", " if p.time >= self.annot.time + self.c.rrmin and self.peaktype(p) == 1:\n", " \n", " if p.amp > self.c.qthr:\n", " ok = True\n", " if ver:\n", " print(\"p.amp:\",p.amp,\" > self.c.qthr:\",self.c.qthr) \n", " self.rr = p.time - self.annot.time\n", " q = self.find_missing(self.r, p)\n", " if self.rr > self.c.rrmean + 2 * self.c.rrdev and \\\n", " self.rr > 2 * (self.c.rrmean - self.c.rrdev) and \\\n", " q is not None:\n", " p = q\n", " self.rr = p.time - self.annot.time\n", " self.annot.subtype = 1\n", " self.rrd = self.rr - self.c.rrmean\n", " if self.rrd < 0:\n", " self.rrd = -self.rrd\n", " self.c.rrdev += (self.rrd - self.c.rrdev) >> 3\n", " if self.rrd > self.c.rrinc:\n", " self.rrd = self.c.rrinc\n", " if self.rr > self.c.rrmean:\n", " self.c.rrmean += self.rrd\n", " else:\n", " self.c.rrmean -= self.rrd\n", " if p.amp > self.c.qthr * 4:\n", " self.c.qthr += 1\n", " #print(\"QTHR INCREASED: \",self.c.qthr)\n", " elif p.amp < self.c.qthr:\n", " self.c.qthr -= 1\n", " #print(\"QTHR DECREASED: \",self.c.qthr)\n", " if self.c.qthr > self.c.pthr * 20:\n", " self.c.qthr = self.c.pthr * 20\n", " #print(\"QTHR FIXED: \",self.c.qthr)\n", " self.last_qrs = p.time\n", "\n", "##################################################################################################################\n", " #ANNOTATION, HERE THE VALUE TO RETURN\n", " self.annot.time = p.time - self.c.dt2\n", " #self.annot.time = p.time - self.c.dt2\n", " self.annot.type = \"NORMAL\"\n", " qsize = int(p.amp * 10.0 / self.c.qthr)\n", " if qsize > 127:\n", " qsize = 127\n", " self.annot.num = qsize\n", " self.putann(self.annot)\n", " toReturn.append([\"NORMAL\", p.type, self.annot.time, False, ok])\n", " #self.annot.time += self.c.dt2\n", " self.annot.time += self.c.dt2\n", " # [type, 0, time, False, ok]\n", "##################################################################################################################\n", " \n", " #Look for this beat's T-wave\n", " self.tw = None\n", " self.rtdmin = self.c.rtmean\n", " q = p.next_peak\n", " while q.time > self.annot.time:\n", " self.rt = q.time - self.annot.time - self.c.dt2\n", " if self.rt < self.c.rtmin:\n", " # end:\n", " q = q.next_peak\n", " continue\n", " if self.rt > self.c.rtmax:\n", " break\n", " self.rtd = self.rt - self.c.rtmean\n", " if self.rtd < 0:\n", " self.rtd = -self.rtd\n", " if self.rtd < self.rtdmin:\n", " self.rtdmin = self.rtd\n", " self.tw = q\n", " # end:\n", " q = q.next_peak\n", " if self.tw is not None:\n", " tmp_time = self.tw.time - self.c.dt2\n", " tann = gQRS.Annotation(tmp_time, \"TWAVE\",\n", " 1 if tmp_time > self.annot.time + self.c.rtmean else 0,\n", " self.rtdmin)\n", " # if self.state == \"RUNNING\":\n", " # self.putann(tann)\n", " self.rt = tann.time - self.annot.time\n", " self.c.rtmean += (self.rt - self.c.rtmean) >> 4\n", " if self.c.rtmean > self.c.rtmax:\n", " self.c.rtmean = self.c.rtmax\n", " elif self.c.rtmean < self.c.rtmin:\n", " self.c.rtmean = self.c.rrmin\n", " self.tw.type = 2 # mark T-wave as secondary\n", " \n", " self.r = p\n", " self.q = None\n", " self.annot.subtype = 0\n", " elif tPres - self.last_qrs > self.c.rrmax and self.c.qthr > self.c.qthmin:\n", " self.c.qthr -= (self.c.qthr >> 4)\n", " #print(\"QTHR DECREASED: \",self.c.qthr)\n", " # end:\n", " p = p.next_peak\n", " \n", " elif tPres - self.last_peak > self.c.rrmax and self.c.pthr > self.c.pthmin:\n", " self.c.pthr -= (self.c.pthr >> 4)\n", " #print(\"PTHR DECREASED: \",self.c.pthr)\n", " self.j += 1 \n", " else:\n", " self.j = 0\n", " self.i = 0\n", " #print(t)\n", " \n", " if len(toReturn)==0:\n", " toReturn.append([0,noneVal,0,noneFlag,ok])\n", " \n", " return toReturn" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:18:54.173810Z", "start_time": "2019-08-18T13:43:27.471921Z" }, "scrolled": true }, "outputs": [], "source": [ "p = gQRS(adc_zero = 1024)\n", "foundAllTh = []\n", "l = 1\n", "padT = 0\n", "padV = 0\n", "tailLen = 100\n", "tai=[]\n", "poi = []\n", "for hFiles in data:\n", " print('---------',l,\"/\",len(data))\n", " l+=1\n", " foundAll = []\n", " n=1\n", "\n", " verbose = False\n", "\n", " for h in hFiles:\n", " startT = time.time()\n", "\n", " print(n,\"/\",len(hFiles), end = \" \")\n", " n+=1\n", "\n", " p.__init__()\n", "\n", " times,points = h[0].values,h[1].values\n", " found = []\n", " #for tm,pt in tqdm(zip(times,points),total=len(times)):\n", " for tm,pt in zip(times,points):\n", " #print('-----------------------------------')\n", " \n", " a = p.isPeak(tm,pt,ver=verbose)\n", " if verbose:\n", " print(\"*\"*50)\n", " #print(a)\n", " if type(a[0]) != list:\n", " if a[3]:\n", " found[-1] = a\n", " elif a[4] == True:\n", " found.append(a)\n", " else:\n", " for elem in a:\n", " if elem[4] == True:\n", " found.append(elem)\n", " #TAIL\n", " for i in range(1,tailLen):\n", " tm = times[-1]+10*i\n", " pt = points[-1]\n", " tai.append(tm)\n", " poi.append(pt)\n", " a = p.isPeak(tm,pt,ver=verbose)\n", " if verbose:\n", " print(\"*\"*50)\n", " #print(a)\n", " if type(a[0]) != list:\n", " if a[3]:\n", " found[-1] = a\n", " elif a[4] == True:\n", " found.append(a)\n", " else:\n", " for elem in a:\n", " if elem[4] == True:\n", " found.append(elem)\n", " lastT = tai[-1]\n", " for i in range(1,60):\n", " tm = lastT+i*20\n", " pt = points[-1]+50*(-1)**i\n", " tai.append(tm)\n", " poi.append(pt)\n", " a = p.isPeak(tm,pt,ver=verbose)\n", " if verbose:\n", " print(\"*\"*50)\n", " #print(a)\n", " if type(a[0]) != list:\n", " if a[3]:\n", " found[-1] = a\n", " elif a[4] == True:\n", " found.append(a)\n", " else:\n", " for elem in a:\n", " if elem[4] == True:\n", " found.append(elem)\n", " foundAll.append(found)\n", " end = time.time()\n", " print(\"---> Elapsed: \",end - startT)\n", " foundAllTh.append(foundAll)\n", "foundSampleOneComplete = foundAllTh[0][0]\n", "print(\"FINISHED\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:18:56.158808Z", "start_time": "2019-08-18T14:18:56.156311Z" } }, "outputs": [], "source": [ "len(foundSampleOneComplete)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:18:58.277747Z", "start_time": "2019-08-18T14:18:58.123114Z" } }, "outputs": [], "source": [ "t = data[0][0][0].values\n", "v = data[0][0][1].values\n", "\n", "tPeak = np.array(np.asarray(foundSampleOneComplete)[:,2], dtype=np.int64)\n", "startIdxMarks = (np.where(tPeak>= start) [0]).min()\n", "stopIdxMarks = (np.where(tPeak <= stop) [0]).max()\n", "foundSample = foundSampleOneComplete[startIdxMarks:stopIdxMarks+1]\n", "\n", "fig = plt.gcf()\n", "fig.set_size_inches(18.5, 10.5)\n", "\n", "plt.plot(t[startI:stopI]/360,v[startI:stopI])\n", "for i in range(len(foundSample)):\n", " plt.axvline(x=foundSample[i][2]/360,color = 'red',linewidth = .5)\n", " plt.axvline(x=foundSample[i][2]/360+0.075,color = 'green',linewidth = .5)\n", " plt.axvline(x=foundSample[i][2]/360-0.075,color = 'green',linewidth = .5)\n", "plt.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:19:00.694014Z", "start_time": "2019-08-18T14:19:00.266205Z" } }, "outputs": [], "source": [ "toSaveAllTh = []\n", "l = 1\n", "for allFileOneTh in foundAllTh:\n", " print(l,\"/\",len(foundAllTh))\n", " l+=1\n", " toSave = []\n", " for sample in allFileOneTh:\n", " sampleToSave = []\n", " for elem in sample:\n", " sampleToSave.append(MITA.MITAnnotation(time = elem[2]))\n", " toSave.append(sampleToSave)\n", " toSaveAllTh.append(toSave)\n", "toSaveAllTh[0][0][0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:19:04.312749Z", "start_time": "2019-08-18T14:19:04.310770Z" } }, "outputs": [], "source": [ "cwd = os.path.dirname(os.getcwd())\n", "sourceDir = cwd+'/data/dataRaw/'\n", "outDir = cwd+'/data/results/'\n", "tempdir = cwd+'/temp/'" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:19:09.361474Z", "start_time": "2019-08-18T14:19:08.158953Z" } }, "outputs": [], "source": [ "# All the Files, independent number of thresholds\n", "l = 0\n", "for th in toSaveAllTh:\n", " print(\"----------------\",l+1,\"/\",len(toSaveAllTh))\n", " print(thresholdComplete[l])\n", " ext = 'out'+str(thresholdComplete[l])+'bits'+str(hyst[0])+'HystLvlCrossing_gQRS'\n", " evalFile = 'val'+str(thresholdComplete[l])+'bits'+str(hyst[0])+'HystLvlCrossing_gQRS.txt'\n", " for sample,name in zip(th,selectedFileNames):\n", " pathOut = sourceDir+name+'.'+ext\n", " #print(\"saving \",pathOut)\n", " MITA.save_annotations(sample,pathOut)\n", " print(\"---------------- EVALUATING ----------------\")\n", " command = cwd+'/scripts/eval_mitdb.sh -f ' + evalFile +' -o ' + ext\n", " command += ' -s '+sourceDir+' -d '+outDir\n", " r = os.system(command)\n", " print('Executed: ',command,'\\nResult: ',r)\n", " l+=1\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Re-sample" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "data\n", "thresholdComplete = [4,5,6,7]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:42:41.648555Z", "start_time": "2019-08-18T14:42:41.643126Z" } }, "outputs": [], "source": [ "#RESAMPLER\n", "def getMidPoints(pt1,pt2):\n", " m = (pt2[1]-pt1[1])/(pt2[0]-pt1[0])\n", " q = pt1[1]-pt1[0]*m\n", " t1 = pt1[0]\n", " t2 = pt2[0]\n", " points = []\n", " for t in np.arange(t1,t2):\n", " points.append(m*t+q)\n", " return points" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:54:09.020609Z", "start_time": "2019-08-18T14:51:42.681397Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "th\n", "th\n", "th\n", "th\n" ] } ], "source": [ "dataResam = []\n", "for th in data:\n", " dataAllFileResamp = []\n", " print(\"th\")\n", " for file in th:\n", " dataFileResamp = []\n", " p = file.values\n", " for i in range(1,len(p)):\n", " p1 = p[i-1]\n", " p2 = p[i]\n", " res = getMidPoints(p1,p2)\n", " dataFileResamp.extend(res)\n", " if(len(dataFileResamp)<649999):\n", " p1 = p[-1]\n", " p2 = [649999,p1[1]]\n", " res = getMidPoints(p1,p2)\n", " dataFileResamp.extend(res)\n", " dataAllFileResamp.append(np.array([list(range(len(dataFileResamp))),dataFileResamp]))\n", " dataResam.append(np.array(dataAllFileResamp))\n", "dataResam = np.array(dataResam)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T14:54:12.846150Z", "start_time": "2019-08-18T14:54:12.843025Z" } }, "outputs": [ { "data": { "text/plain": [ "(4, 48, 2, 649999)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dataResam.shape" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T15:00:21.769117Z", "start_time": "2019-08-18T15:00:21.682245Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvXd8HOW1//8+29Rly5ZckHvHxgZsY0wwLYUeCBj4UnIhDUICpNwkNxByUy+/QJKbe28SSgghtAAhFUMI4CQGA8YYG3DFRbjKRZYtq0tbn98fO7tarXZVV9Jq5rxfL700enZ292hmdz5zynMeMcagKIqiOA/XYBugKIqiDA4qAIqiKA5FBUBRFMWhqAAoiqI4FBUARVEUh6ICoCiK4lBUABRFURyKCoCiKIpDUQFQFEVxKJ7BNqAzSktLzaRJkwbbDEVRlCHFunXrjhhjyrraL6sFYNKkSaxdu3awzVAURRlSiMie7uynISBFURSHogKgKIriUFQAFEVRHIoKgKIoikNRAVAURXEoKgCKoigORQVAURTFoagAKI7gpc2HOFzfOthmKEpWoQKg2J5QOMLnH1/HVb96c7BNUZSsQgVAsT3+UASA3UebB9kSRckuVAAU2xMTAEVR2qMCoNiegAqAoqSkSwEQkYdF5LCIbEoY+6GIbBCR90TkZRE5zhoXEfm5iFRYj89PeM4NIrLD+rmhf/4dRemIPxQebBMUJSvpjgfwCHB+0thPjDHzjDEnAc8D37HGLwCmWz83AfcDiMgI4LvAqcAi4LsiUtJn6xWlG6gHoCip6VIAjDErgZqksfqEPwsAY21fCjxmoqwGhovIWOA8YLkxpsYYcwxYTkdRGRACoQjPvrcfY0zXOyu2IDEHoGKgKG30ej0AEbkLuB6oA86xhsuBfQm7VVpj6cYHnJ8t384Dr35Aca6Xc2aNGgwTlAEmUQCONQcYXZw7iNYoSvbQ6ySwMeZOY8x44HfArdawpNq1k/EOiMhNIrJWRNZWV1f31ry07KhqALQyxEkk5gCONPoH0RJFyS4yUQX0JLDU2q4Exic8Ng440Ml4B4wxDxpjFhpjFpaVdbmiWY9p8IcAcLtSaZJiRxLFvqYpMIiWKEp20SsBEJHpCX9eAmy1tpcB11vVQIuBOmPMQeAl4FwRKbGSv+daYwNOQ2tUAJosIVDsT6BdCCg4iJYoSnbRZQ5ARJ4CzgZKRaSSaDXPhSIyE4gAe4Cbrd1fAC4EKoBm4NMAxpgaEfkh8La13w+MMe0SywNFLBzQqALgGBI9gKCG/hQlTpcCYIy5JsXwb9Lsa4Bb0jz2MPBwj6zrB1wSDf2oB+AcEj2AYFgFQFFiOG4mcCz0rwLgHBKTwCoAitKG4wQgFg5o9OvsUKeQ6AEEwjr/Q1FiOE4AWgLRC796AM7BryEgRUmJ8wQgGBWAYEQvBE4hoElgRUmJ4wSg1RKAkIYCHIM/FI7nftQDUJQ2et0KYigSDEcIWhf+kHoAjiEQipDjcROOGM0BKEoCjvIAYnf/QFwIFPvjD0XI8brwukU9AEVJwFEeQEuCAIT0QuAYoh5A9F5HBUBR2nCUALQG2r78oYh6AE7BH4rg87iIGBUARUnEWQKgE4IcSSwHEIlAIKTCrygxHJUDSCwH1Cog5+APhfG5NQegKMk4SgASv/xBDQE5hrYksEurvxQlAYcJQPSin+t1aRLYQfhDEcsDcGkIqBPW7anhlLv+waG61sE2RRkgHCYA0Yt+vs+jISAHEQhFyPG68XpcGgLqhM0H6qlu8PPS5kODbYoyQDhKAGI5gDyvW1tBOIiYB+DTHECnHG2Mrpa2fEvVIFuiDBTOEoC4B+BWD8BB+EPheA4goL2A0nKsOSoAq3cepb5VV05zAo4SgHgIKMejOQAHEQhFyHG7cLtE5390wtGmQPwYvbqterDNUQYAZwqA161VQA4isQoorOc9LTWNAeaNG8bIAp+GgRyCswTAqgCJhoDUA3AKASsH4HZpDqAzGvxBhud5+fCsUazYdliPlQNwlAD4rQ90nuYAHEU0B+DG6xb1ADqhORAmP8fDx2aPpqE1xJpdNYNtktLPOEoAYouB5Pu0CsgpGGMSPACX5gA6odkfJt/rZsn0UnI8Lg0DOQBnCYDOA3AcoYghYiDH48LjEp0J3AnNgRD5Pjf5Pg9nTC9l+ZYqjNHviZ1xqAC4CUWMfrgdQKzs02cJQFiFPy0twTB5vmh/yLNmlLG/toWDOivY1jhKAGKrQeV53YC2hHYCsQXhczwuPG7R6q80xFbLy/dFvxulhTkA1LXofAA74ywBsGLBHnf039YwkP1p8wDceFxaBpqO5kC0VXpMAIpyvQA0tIYGzSal/3GUAATDEbxuwWOtEB7WEJDt8VtrQOR4tAy0M5oD0Qt9vhUCKs6L/q5XD8DWOE8APC5cMQFQD8D2JOYAtAw0PWk9AL8KgJ1xnAD43C71ABxEYg7A7XJp2C8NLZYA5MUFIOYBaAjIzjhKAAIhg9fd5gFoSaD9iQuA161loJ0Q8wAKrBBQTAAatCmcrXGWAIQj8XJAAL0W2J9YDiCa/BciBiIaBupAk5UDiHkAOR43Po9Lk8A2x1ECEAxFk8Bu9QAcQyDuAbQJv5b/dqQlKQcQ224JhgfLJGUAcJYAhCN43S7cYuUA9EJge2IhoMTyXz3vHUlOAkN0vkxMGBR74igBiIeA3CoATiHmAeQmeADaB6ojLUlloGAJgHoAtsZRAhDzAFzqATiGNg/A3Vb9pZVAHUjpAfjUA7A7jhKA+ExgLQN1DIk5ALcVAlIPoCNNsTJQb1IISD0AW9OlAIjIwyJyWEQ2JYz9RES2isgGEfmLiAxPeOwOEakQkW0icl7C+PnWWIWI3J75f6VrgmGD1y1tZaB6J2h72lUBudTzS0dLIESut61EGiwPQAXA1nTHA3gEOD9pbDlwgjFmHrAduANARGYDVwNzrOfcJyJuEXED9wIXALOBa6x9B5RYCEgvBM7Bn6oKSIW/A82BcLv4P0CuJoFtT5cCYIxZCdQkjb1sjIkVCK8GxlnblwJPG2P8xphdQAWwyPqpMMbsNMYEgKetfQeUWBLYrSEgxxBoVwWkZaDpiAqAu91YntdNq3oAtiYTOYDPAH+3tsuBfQmPVVpj6cYHlFgrCLd6AI7BHwrjdgketwuPK1YGqjmAZOpaggzL87Yb0xyA/emTAIjInUAI+F1sKMVuppPxVK95k4isFZG11dXVfTGvA4GQNQ9ABcAxxBL/QFsZqIaAOlCfSgC0Csj29FoAROQG4GLgOtO2tFYlMD5ht3HAgU7GO2CMedAYs9AYs7CsrKy35qUkGDZ4PaITwRyEPxQhx2sJgE4ES0tdS5Di3BQCoB6AremVAIjI+cA3gUuMMc0JDy0DrhaRHBGZDEwH1gBvA9NFZLKI+Igmipf1zfSeE7Q8AI0FO4fUHoCGgJKpbw3G1wCIked1EwwbPV42xtPVDiLyFHA2UCoilcB3iVb95ADLJXo3vdoYc7MxZrOIPANsIRoausUYE7Ze51bgJcANPGyM2dwP/0+ntCWBoxcEbQpmfxI9AA39pae+JZQyBwDQGgzjdTtqypBj6FIAjDHXpBj+TSf73wXclWL8BeCFHlmXYeJJYFEPwCm08wDU80tJIBShJRjuEALKtaqCWoLh+AIxir1wjKyHwhEiBk0COwx/KIzPE72QxaqAdB5Ae+qtnv/D8tN4AAENAdkVxwhArPJDBcBZ+EMRcjzJHoBe0BKps9b9TfYA8hM8AMWeOEYAAlYiK3E9AJ0IZn8CoWjeB9CZwGmILfyeLgcQWzBesR+OEYBYJUOOJ7EVhN4J2p1AOMEDiIWA1PNrR7216ldyFVCuVz0Au+M4AUgMAemdoP2JTf4DDQGloy6dB+BrqwJS7IljBCDWEyZRACIaArI9scov0DLQdNSnyQHEQkAtmgS2LY4RgLgHkNAMTkMB9icxB+DVKqCUxJPAaXIAGgKyL44RgEAo+qX3JSSBdSKY/UkUALeGgFJS3xrE53HFY/4x8rQKyPY4RgBiHoAvIQmsHoD9ic3+BvDqeU9JqlnAkCAAWgVkWxwnAF5326pHGgu2P/5QxxyAhoDaU98SpDi3Y1OAXEs4NQdgXxwjAIlJYF0RzDkE2k0E0zLQVNS3dmwFDdHj5XO7NARkY5wjAIkegPYCcgTGGALhhDLQuAegd7SJ1LUEOySAY+R6XVoGamMcIwCxVhCJi4NrEtjehCMGY2ibCazN4FKSajGYGLoojL1xkAC0JYG1DNQZBBLOOZCwJKSe90SONgUYnk4AdFlIW+M4AfC6BRHBJXohsDuJC8IDWLqvIaAEWoNhGlpDjCrOTfl4ns9Ds3oAtsUxAuBPSAJD9G5Qm8HZm7gAWB6AiOB1i3p+CVQ3+AEoK8pJ+Xie5gBsjWMEIJgUDnC51AOwO/4kAYBoKagKQBuHG1oBGJVOAHRdYFvjHAFI5QHohcDWBBI6wMbwulw6DyCB6oYAAKWF6TwATQLbGecIQKwKyNMWD1YBsDeBJNGHaDsIbQXRRqPfagWdZsnHXK9bQ0A2xjECkLggDEQnueiFwN7Ew34JAuBxuTQElECjtRxkQY475eNaBWRvnCMAsbtBV1tbAC0GsTfJSWCITgYLawgoTpMV3ilM0QoCostCahWQfXGMAATDETwuifcBcovoimA2J5UAuF1CUM97nEZ/CK9byPGk9gByNQlsaxwlAO1iweoB2B5/uKMAeN2iuZ8EGltDFOSkvvuHaAgoEIroMbMpDhIA0+FOUD0Ae5M8EQysMlANAcVp8oco8HUuAKDLQtoVxwhAIMkD8LgEvQ7Ym5gAtCsD1eR/Oxr9IYrSxP9BF4WxO84RgFAEn1UBBOoBOIGUZaDqAbSj0d95CCg3vi6wCoAdcYwABMMRvMkzQvVCYGuSZ39DrPxXz3uMpi4EIF89AFvjLAFIuhOMaC8gW5PcDRSioT8NAbXR6A9R1EUSGNQDsCuOEYBAyHRMBuqdoK1JVwaqnl8b0RBQ6hJQSBAA9QBsiWMEIFUISEvb7I0/RRWQloG2p8kf7jwHoCEgW+MYAUhOAntUAGxP6jJQF0E970B0RbymQPdCQK0aArIljhGA5ByASzQEZHeipb9ts78BvFr9Fac5GMYYupwIBuoB2BXHCoDHLbomsM0JhNqfc9AcQCJNVifQ7lQBaT8ge+IYAQgkzQRWD8D+BEKRdpPAIDYRTM87tLWCLuxGDkBnAtsTxwhAMBxJagusZaB2JxCKtBN9iHkAGgICqG+JtoIuztMyUKfSpQCIyMMiclhENiWMXSkim0UkIiILk/a/Q0QqRGSbiJyXMH6+NVYhIrdn9t/ommg4oP1MYA0F2JtAuKMAeLT8N059a9QDGJaXejEYiHpMHpdoDsCmdMcDeAQ4P2lsE3A5sDJxUERmA1cDc6zn3CcibhFxA/cCFwCzgWusfQeM1N1A9UJgZ/yhcIc2xx4tA41TF/MA0qwGFkMXhbEv6X0/C2PMShGZlDT2PoCIJO9+KfC0McYP7BKRCmCR9ViFMWan9bynrX239MX4npA8D8DjchHWEJCtiZb+JoeAXPHlQZ1OLATUmQcA0TyA5gDsSaZzAOXAvoS/K62xdOMDRvLFwKUegO3xp8gBRCeCaQ4AEjyALgRAVwWzL5kWgA4uAWA6Ge/4AiI3ichaEVlbXV2dMcOS1wPQiWD2J5UAaO6njfrWID6PK97xMx15XrcmgW1KpgWgEhif8Pc44EAn4x0wxjxojFlojFlYVlaWMcNik4JiuEQFwO5oGWjn1DYFuwz/QLQltOYA7EmmBWAZcLWI5IjIZGA6sAZ4G5guIpNFxEc0Ubwsw++dlnDEEI6YDgvCaFdIe5NKANx63uMcafRTVpjT5X55Xs0B2JUuk8Ai8hRwNlAqIpXAd4Ea4BdAGfA3EXnPGHOeMWaziDxDNLkbAm4xxoSt17kVeAlwAw8bYzb3xz+Uilhf+HZVQG5dE9juaBlo5xxp9FNa1A0B8Lk53BAcAIuUgaY7VUDXpHnoL2n2vwu4K8X4C8ALPbIuQ8QEIPFu0C2aDLQ7/lC4QxWQx+XCmGgjtMQeQU6kusHPtFFFXe6nOQD74oiZwLGyP50H4CyiIaCO8wAAgg4Xf2MMRxoDlBb5utw3z6cCYFccIQDp1oZVAbA3qVpBeKy7fqef+/qWEIFwpNs5AE0C2xNHCEBbDiBpPQCdCGZr0vUCAhw/Gay60Q9AWTdzACoA9sQRApBqbVj1AOxP6olg0b+dfu6rGywB6IYHkOt10xqMaPt0G+IIAYh5ALomsHOIRAyhiElZBgo4viPoEcsD6FYVkDVRLLbEpmIfnCEAodRJ4Fg1iGI/Unl90JYDcLr4xwWgWzmA6DHUMJD9cIQABMLRD643qQwU0DyATUm1IDyAR0NAQGIn0C4rwcn3RfdpDoT61SZl4HGGAMQ9gIT1ANxaDWJn/KGo6Ock9bnxxJPAzg5n1LUEKczxxAWxM3RVMPviDAFIkQPQckB7Eyv9zengAeh5h2gZaHf6AEHiqmDOFk074ggBiE1iSex66BKNBduZmACkywE4vQy0riVIUTfCP5AgAAkeQDAcYc2umn6xTRk4HCEAMdc139cmALELgSaB7Yk/rQBoDgCiraC77QH4OiaBH1y5k6t+9SZv7TzaL/YpA4MjBCC2mEUsmQUJ5YAOvxDYlXgIKLkMVFtBANHVwLpaCCZGboqF4Q/UtgCw5WB95o1TBgyHCEC0eiEvIQTktu4EI1oFZEu6KgN1vAfQ0n0PIHbjlJgEHmmVj1bV+zNvnDJgOEIAYh/cvBQhIPUA7EkgXRmoJfxOXxWsvjXU5WLwMWI3TonLQsaO776a5swbpwwYjhCA5kAYj0va3Q3GWgGHHX4hsCuxMtAOHoA7JvzODQGFwhEa/b2oAkrwABpao/MIKg43Zt5AZcBwjADkpakH14lg9qQtB5D6vDvZ82tojYZEi/O6VwWUayWBW9sJQPQ1dh1pwuh3aMjiCAFoDYbbhX8gwQNw8J2gnemqCsjJIaB66+69ux6Az+3C65b4RR/aPIBAOBLPtyhDD0cIQHMg3K4EFBKTgYNhkdLfpKsCapsI5twT39YGonsCICIMy/PFnwe0E4PWoHOP5VDHMQKQmxQCapsIph9eOxKLV+d4dSJYMnEB6KYHADAsz0NdSyD+d6IA+LVFxJDFEQLQ5A9RmNM+3qnlgPamrjl1mMOt552Dda0AjCnO7fZzhucnewBt29omeujiDAEIhChMmvauzeDsTV1LkHyfu0MSONYS3MlJ4P3HWhCBMcN6IAB5XmqbEwTAH4qvJqZN4oYujhCAxtYQBUkeQLwdtIMvBHamNs1EJ10QBvbXtjCqKKdDgrwzhuV54x5AJGJo9Ifiq4lpDmDo4gwB8Ico0hCQo6hLIwBt8wCce973H2uhfHhej55TVpTD4QY/oXCEpkAIY9pWE4vNuVCGHo4QgCZ/Rw/ApQJga+qa0whAvAzUuXetB+paKC/J79FzZowuIhCKsPtoczwBrB7A0Mf2AhCJGJoC4Q4CoBOC7E1tS4Dh+eoBJBOJGA7WtvbYA5g1tgiArYfq4wJQWuQDNAcwlLG9ADRZjeCSQ0BunQlsa+paggzP83UYd7rwVzf6CYQjlJf0TACmjSrE7RK2HmygtjlaDjq6KJpE1iqgoYvtBaDRHxWADklg7QVka2qbgwxL4QE4vQy08li0jXP58O5XAEG0pcaU0gK2HmrgmFUNNNaqIlIPYOhiewFosgSgQxmow+8E7UxrMIw/FOk0B+BUAdhfGxOAnuUAAGaNLWbrofq4BxArI23VJPCQxfYCEItXFuakrgd36oXAzsTKFVMJgKX7jhX+/TEPoIchIIBZY4qoPNbCvmPRFtAxAfBrEnjIYnsBaPJH704Kc1LPCNVWEPYjNmEpVRJYRHC7xLG9gPbXNjMsz9thZnx3mDUmmghevbMGn8dFSb6VBFYPYMhiewFoywEkeQDaFdK2xEIUqZLAEBX/bPMADtW1xpdZ7E96MwcgxqyxxQCs23OMkQW+eKM9LQMduvT8NmCIEROAomQPQBcGsS3VjdFlCmNlisl4XEIkiwSgNRjm8vveoKrBz8fnjeULZ09jpnW3nWn217YwYURBr5573LBcinI9NLSGmFs+DBEhx+PSiWBDGNt7AE1pPACnlwPameqGqACMKkpd6ZJtHsATq/dwoK6Vi+aO5eUtVZz3vyv53KNvs27PsQ77+kNhfvvGrvjnuicYY9h/rIVxvYj/QzR8Fgv7LJo8Aoi229YcwNDFMR5A2olgGgKyHYcb/HjdwvA07Y49Lsma5H9Da5B7V1RwxvRSfn7NyRxrCvDom7t5ZNVult6/ikWTR/DFs6dy1owyRIQn39rL95/bQpM/xK0fnt6j9zrc4KcpEGbSyJ5XAMWIrbOwcFJUAHK9bi0DHcLY3gNo9IfwuqXjwiAu7QppVw7X+yktzIm3+0jG7XJlzXn/9Wu7ONYc5BvnzQSgpMDHVz46g1W3f5j/vHg2+2qa+dRv3+biX7zOsvUHeHDlTgAee3NP/GLcXWLr904f3fvw0v9efRLnzh7NCcdF8wE5XpdOBBvCdCkAIvKwiBwWkU0JYyNEZLmI7LB+l1jjIiI/F5EKEdkgIvMTnnODtf8OEbmhf/6djsQ6gYq0vxjoylD2pbrRzyirUVkq3K7smAB4pNHPQ6/t5MK5Y5g3bni7x/J9Hj67ZDKvfuMcfnzFPFqCYb701LscrGvls0smc7jBz/MbDvTo/WICMG1UYa9tXjxlJA9evxCPVUad61EPYCjTHQ/gEeD8pLHbgX8aY6YD/7T+BrgAmG793ATcD1HBAL4LnAosAr4bE43+JtViMNBWBurklaHsyuH6VsrSxP8h6v1lgwdw74oK/KEIXzt3Ztp9fB4XVy0cz/KvnsX9183nWxfO4s4Lj2dqWQGPrtrdo/erONxIUY6nU3HsKRoCGtp0KQDGmJVATdLwpcCj1vajwCcSxh8zUVYDw0VkLHAesNwYU2OMOQYsp6Oo9AuNaQRAJ4LZl+oGf3yxklS4XUJkkHtA7a9t4Xer93LlgnFMLev6jtztEi6YO5abzpyKyyVcf9ok1lfW8d6+2m6/Z8XhRqaOKuzgDfeFaBWQetFDld7mAEYbYw4CWL9HWePlwL6E/SqtsXTj/U46AYjPCB3ktsChcIQ7/ryBZ9/bD8DGyjqOWmWMSs8JhiPUNAc6vcv1ZEEV0H0rKjAYbvtIzxK5MS6fX06Bz81j3fQCQuEIG/fXMduK3WcK9QCGNpmuAkp1a2E6Ge/4AiI3EQ0fMWHChD4b1OQPMTy/Yz24iOB1D/6F4Of/3MFTa/ax7L0D/Pmd/by6vZpJI/P57BlTuPzk8g7VS0rnbKisxRgYPyJ9pctgzwQ+UNvCM2v3ceXC8b2elFWU62XpgnE8vWYfd150PCMLOw/rbNxfR6M/xIemjuzV+6Uj1+uipkk9gKFKbz2AKiu0g/X7sDVeCYxP2G8ccKCT8Q4YYx40xiw0xiwsKyvrpXltNPg7rgccY7DrwVdVHOEXKyo4Z2YZgXCEV7dXA7D7aDP/+ddN3PzEOp1k00Meem0Xw/K8XHDCmLT7uF0yqOW/97/yAQBfPHtqn17n+tMmEghHeGL13i73fXPnUSCaxM0kOR63toIYwvRWAJYBsUqeG4BnE8avt6qBFgN1VojoJeBcESmxkr/nWmP9TpM/RKEvtQB4Xa5BuxAcafTz5d+/x+TSAn557XzuvPB4Lpw7hsvnl3P6tJF85aPTeW3HEb72zHrNU3STPUebeHHzIT65eEKnnpPHPXjzAA7WtfD7t/dxxYJxjOvhqlzJTBtVxLmzR/PLFTtYsys5TdeeNyqOMGtMEaVdeAo9JcerE8GGMl3GF0TkKeBsoFREKolW89wNPCMinwX2Aldau78AXAhUAM3ApwGMMTUi8kPgbWu/HxhjOv/EZoBAKEJ9S8flIGO43TIorSAiEcPXnllPXUuQxz6ziIIcD586fTKfOn1yu/3yvG5+9PetlOT7+MGlczKavLMjv3l9F16XixtOm9Tpfm4ZPM/vV6/uJGIMXzx7WkZe7ydXnMhl97/B5x9fy7O3LGFCiklezYEQb+86xg0fmpiR90wk1+tWL3UI06UAGGOuSfPQR1Lsa4Bb0rzOw8DDPbKujzy6ajctwTBnzihN+fhglQP++rWdvLq9mh9+4gSOH5s+Kff5s6ZS0xTgVyt3UlaUw5d6mTB0AjVNAZ5Zu49PnHwco4o7X+xksKqAqupbeXLNXi6fX95pjqInDMv38psbTuET977BZx99mz9/8UMU5bafAf3WrhoC4QhnTO97SDWZ6DwA9QCGKradCXy4vpX/++cOPjxrFGfPHJVyH49LBrwK6J29x/jJS9u44IQxfPLUrpPct18wi8vnl/Oz5dt5ek3XsV6n8sTqPbQGI3zujCld7usZpNDfA69+QDhiuPWczAr55NIC7r9uPruONPGlp97tEN5aub2aHI8r3r8nk+R6XVoFlGECoQh3/W0LB+v6vzusbQXg7he3EghF+M7Fs9Pu4xngKqC6liBfeupdRhfncvfSed0K6YgI9yydx9kzy/jWXzayfEsVkYhh3Z5jGIeuZ7x+Xy2twbC1QHmQ1mCYR1ft5pyZZczoRpsD9yD0Aqqqb+XJt/Zy2cnlKcM0feVD00r5/qVzWLGtmh+98H67x17bcYRTp4wk1+tO8+zek+NxE4qYQS+nthOvbq/m16/t4tt/2dT1zn3EljWGO6sb+fM7+/nC2VOZVJq+9a1ngKtBHlz5AftrW/jjzR9KuVpVOrxuF/ddN59rHlzNrU++w6LJI3htxxH+5/+dyGUnj8uojU+v2cuTa/Zy77XzMxKmMMbw0uZDnDp5JCUFqdszpyMYjnCgtoW9Nc3srWlm1phiWoNhrnvoLRZNHsG7e48xfVQRF5wwhqNNAW48s+u7f4gK/0DGrVdsPcwdf96IMXDLOZmJ/afiulMnsqOqkYde38VDr+/ilnOmcu2pE6k43MjVp4zv+gV6Qa43eg/pD0Xi7SGGOqG0sJmDAAAYQklEQVRwhOVbqph9XDFV9X7yvG6mjy7ktqfepa45yMUnjuXCuWPjCfVwxPDse/s5fmwxY4fl0hqM0BoM0xIM0xoMR/8OhfFb27FxfyjC/Akl3P9KBXk+N5db3+U3Ko4AsMNq3dGf2FIAppQV8tSNi5k3blin+3ncrgG7E2wNhnlqzT4+evxoFkzseReMfJ+Hhz91Clc88Cav7Yh+QJ5YvZePzzsuI1+8YDjCfz2/hUff3APA1/6wnqdvXJy2oVp3MMZw94tb+dWrO7lq4Th+fMWJHR6vbQ7GL/B7a5rZl7B9oLaFxNPjdQuTSwsoyvHEq162HKxny8F6TpsyktO6WeI4UB5AXUuQHz6/hT+uq2TG6EJ+9W8LmNzJDUkm+PZFx7OhspZ39tZy74oPeGVbtLT4zBmZj/8Dca+iNRge8nNWIhHD8xsP8r/Lt7PzSBOTSwvYV9PMrLFFjCzIYeWOaqaUFvCdZzfzvWWbOX1aKRfPG8vzGw7Gv5O9IdfrImLghY2H4mMTRuTz3G1LMvFvdcrQPmOdcFo3Jrx4XEJwgFzX59YfoKYpwKc/NKnXrzGyMIenb1rM2t3HaA2G+dof1nPPi1u586L0Ya7ucKwpwC1PvsOqD45y4xmTmTaqkG/+aSO/WrmTL/ShVv1//7GDX726kxEFPp7fcJC55cPYd6yFvUfbLvYNSX3tSwt9jB+Rz4KJJVx2cjRZOmFEPiMKfPzbb95ie1Uj377oeIpzvTy34QDXLJrAgyt38qPL53a7SmogqoD+tbWKO/68kSONAW45Zypf+sh0cjyZD8Ek43G7eOyzp7JyezWvVxzhybf2MqY4l+l9aADXGTEPoHUIt4MwxrB8SxU/W76drYcamDm6iPkThvPO3lpcApv21wPw46XzuOqU8Ww71MCy9ft5bv1Bvvmnjfg8Lr5/yRwg2l04z+sm1+siN/bb4ybH604ad/PajmoeXbWbuy6bi9ctbKisw+0SHn5jNzefOaVHUYLeYlsB6A4DdSdojOGRVbuZMbqwW8LUGaOLc7lo3lgA1lfW8uvXdjF33HAuOfG4Xr3etkMN3PjYWg7VtfLfV57I0gXjMMawcscRfvLSVk4aP7xXNt//ygf83z93cOWCcVy9aAJL71/Ffz67GZ/HxfiSPCaMyOeUSSXxC/yEkfmML8nv9C7yvusW8ORbe7n21Ank+zxcZYU1Lpw7tke29ed5r2sO8oPnt/CndyqZObqIh64/hbldeKKZpjDHw4Vzx3L+nDGUFuZw3LDcfishTvQAhhrGGF6vOMJPX97O+n21TC4t4P+uPomPzzuOTQfqeGTVbk6ZNII7/ryR718yJ/55mzmmiG+MmcXXz53Jpv31DMvz9iqvc+lJ5Vx6UltHnGmjiuLjA4WjBcDjdhEcAAFYt+cYmw/Uc9dlJ2T0i/jti2bz/sF6vvnHDUwfVdhpSWkqlm+p4itPv0t+joenP7+Y+ROioalY4vn9g/Xc9tS7vPClJV2WViby2zd2cc+LW7nkxOO4e+k83C7hH/9+FoVWJ8rehpUWTCzpVfgsmf6aCPaPLVV86y8bOdoU4LYPT+PWD08bkLv+dLhcwr9/bEa/vkfs/xtqk8He3l3DT17axppdNZQPz+OepXNZOn9cPJw6b9xwfnbVSUQihiXTSlPmw0RkwMU909gja9NLvANUBvrbVbspzvVw2cmZVXafx8W9182nOM/D5x9fF18MvSuMMfzyXzu46fG1TB1VyHO3Lolf/GMU5nh44JMLoitPPflut0NlsRWrzpszmv++6sR42+1powoZMyy3TzmFTOF2RXM/H1Q38st/7aA50HF5xUZ/iIbWYLder7Y5wL///j0+99haRhT4ePaW0/nauTMH9eI/ULSFgIaGB7Cxso4bHl7DlQ+8ya4jTXz/kjn86+tn8f9OmZAyl+ZyScbmbGQjjhYAn8dFIBThufUH+NHf3++XssqDdS28uOkQVy+Khi0yzaiiXO67bgEH61r48tPvdbiz3Xygrl3FS2swzG1PvctPX97OpScexzOfP40xw1Lf3c8YXcTdS+eyxrpb6oq/vFvJnX/dyNkzy/j5NSfHW25nG7FuoL9/ex8/fXk7l927il1HmgC475UKnlt/gBsfXcsV97/Z4QbhQG1Lu26tK7Ye5mP/s5Jl6w/wpY9MZ9mtSzihfGjfFfaEmMhlewhoe1UDNz++jo//8nXWV9Zy+wWzWPmNc7jhQ5McIdTpcHQIyOdx0egP8cd1lby6vZphed6MTdHfcqCe7z23mUjEEDGGf1uc+Wn4MRZMLOF7l8zhzr9s4rvLNjG3fBjTRxexfl8t339uC+fMLKMlGGZqWSFjh+Xy/IaDfPP8Wdx81pQuQ1KXnlTOuj3HeHDlTg7VtfL+wXpGF+cyriSP8uF5lFu/RxbmcPufNnLq5BE88MkFWf2lclue35FGP0U5Hg43tHLJL15nYml+POEX45m1lWyvamBYnpczZ5Sx9P5VzBxdxJULx/Gnd/bz/sF6Zo0p4refOsVRF/4YBTnR89wSyE4BONro54fPb+HZ9Qco8Hn4yken85klkynO7f8E61DA2QLgjnoAVfWtAPzkpW0cP7aYc9LMHO4uz60/wDf+uD4+Rf7CuWP63Y28dtEENuyr69AZcnRxDiu2VVOcGy2djBgoH57Xo+qeOy86nvWVdSxbf4CTJwynwR/iH+8f5kiKdQu++tEZ/TLhKJN4rdzPkcYAU0YVcu+1J/PF373Dhsq6dvvNHF3E957bTCAUiVZnvL4LgG1VDfzX39omW33rwuMdefEH4kn7Rn/HMFo28K2/bGTF1mpuOnMKN585tcdzUeyOswXACgEdbQpwxYJxbD5Qzy2/e4d7ls7j472oqglHDD99eRv3v/IBCyeWsGjyCN4/WM/dS+f1g/XtERG+f+kcjjYFOH5sEXOOG8ax5gAfP/E4nl9/gI8cP5qq+lb+vxfe73FL4ByPm4euX8ir26u57OTyeFy/NRhmf20L+4+1sL+2hYgx/dJuINP43BI9741+xhTnMq4knz/cfBqX3buKLQejHsCpk0fwzQtmcfl9qzhlUgnrK+sozvPyhy+cxppdNSyZVsr+2hZ+/OI2Fk4akNVNs5IiSwAaWrNPAFZsO8xLm6v4j/NnZsyztxuOF4D61iB1LUEmlxbwjfNm8sXfvcNtT73L2t01fOui47sdyqhrCfLlp9/llW3VXHvqBL738Tn4PAMbA8/1unnohoUdxq9eFO05VFaUw5M3Lu7Va5cV5XDFgvazjnO9bqaWFXZrScNswudxEbRCQHOsFbJyPG6ev20Jqz44yp/eqeTy+eXMn1DCX285nemjCtlztJlRxTmUFuYwa0z0OVPKCvulwdpQIrbWRrZ5AK3BMN9btpkpZQV8bkn3Zog7EUcLQI7HxZHGaOXMiAIfo4tzefqmxdzz96089Pou3qus495rT+6yb/uOqgZuenwdlceaueuyE7ju1P6L9yt9x+uOCkBdS7DdanEul7BkeilLprd1jz1p/HCAjC+laBfyvG5cEl13I5t4cOVO9hxt5onPnjrgN2JDCUcfGV9ClUqRdSfjdbv49sWzeeCT89l5uJGLf/E6K7YdTvcSvLz5EJfdt4qG1hBP3rhYL/5DgKgAGFqDkZTrRSvdR0QoyPFkVQhoX00z966o4KJ5Y9uJudIRRwtATkKyMvlCcP4JY1l22xLGFOfy6d++zU9f2kY4YqhvDXKgNtqm9fHVe7jp8XVMKSvgudtO55RJ2R//Vmh3R6gC0HeKcjxZFQL6/nNbcLuEb190/GCbkvU4+tPf3gPoWBY2ubSAv95yOt99djO/XFHBO3uPUdscZF9NM7+78VR+/PetnDG9lF9fvzDrK1+UNhLPuwpA34l6AN2bNNff/GtrFf94v4o7LpjF2GF5g21O1uPoT3/inWBRmoXjc71u7rliHgsmlfCff92EPxQh1+viigfeJBCKcPsFs/TiP8TwutvmPhSmOe9K9xme76WuJTsE4InVexlXksenk5ZXVVLj6E9/dwQgxlULx3PiuOFsPlDHiAIfn37kbT48axRzjnNm/fdQxqshoIxSku9jz9HmwTYDiM74PXlCiSZ+u4mjP/09DQXMHFPEzDHRjn1/vPlDTC3r397uSv+QeN6Heg/7bGBkoY939tYOthm0BKLzUq5c0D+L39gRR3/6E+8SCnrYpycTXSmVwUGTwJllRIGPY80BjDH91na6O3xQ3Ygx0caDSvdwtJ+UeCHIhi6VysCQ2KQuXSM8pfuU5PuiFXItg1sJ9EF1dAnF6aNVALqLowVAk7fOJDEENBCrLtmd2B33qg96vyxiJth6qAG3S5jYi8VZnIqjBWCs3v05klgS2K1eX0Y4Y3oZ40ry+O2q3YNqx6oPjnLS+OFZ3Yk223C0AOidgjOJlYGWD9c68Uzgdgk3nDaJNbtq2HKgvusn9AN1zUE2VtayZJrO/O0JjhaAssIcQD0BpxGwFjAfV6ICkCmuWjiePK+bRwfJC3hz5xEiBm390EMcXQIhIvzpCx/SC4HDiM35+PCsvq37oLQxLN/LBSeM4eUth/hRZO6AF1W8XnGEAp873rxP6R6O9gAgWs45ugcLnitDnwUTR/C3Ly3hs0t0tmgmWTx1JMeag/FqnIHCGMPK7UdYPGVk1i5Dmq3o0VIcyZzjhg1qzbodWWQ1Q3xrV82Avu/OI03srWnmbPXoeowKgKIoGWHiyHxGFeWwZoAFYMXWaLv2c2Y6e3Ge3qACoChKRhARFk0ewZpdNRhjBux9V2w7zIzRhV0u3KR0RAVAUZSMsXjKSA7Vt7K9amDyAI3+EGt21XDOTA3/9AYVAEVRMsb5J4zB7RL++t7+AXm/13ccIRg2nKPx/16hAqAoSsYoLczhzOmlPPvufiKR/g8DvbLtMEU5Hm3O2EtUABRFySifOLmcA3WtrN51tF/fxxjDim2HOWNGqZZ/9pI+HTUR+bKIbBKRzSLyFWtshIgsF5Ed1u8Sa1xE5OciUiEiG0Rkfib+AUVRsotzZ4+hMMfDX9/t3zDQloP1VNX7Nf7fB3otACJyAnAjsAg4EbhYRKYDtwP/NMZMB/5p/Q1wATDd+rkJuL8PdiuKkqXk+dycf8IYnllbyY/+/n6/hYJi5Z9naflnr+mLB3A8sNoY02yMCQGvApcBlwKPWvs8CnzC2r4UeMxEWQ0MF5GxfXh/RVGylMtOLgfgV6/u5PWK/mkT/XrFEU4oL2ZUkc7k7y19EYBNwJkiMlJE8oELgfHAaGPMQQDrd8w/Kwf2JTy/0hpTFMVmLJ4yMr79xOo9GX99YwxbDzUwt1x7//SFXguAMeZ94B5gOfAisB7obEmgVPPuO/iGInKTiKwVkbXV1dW9NU9RlEHE7RI2fu9cPrtkMv94v4p9NZldNL660U9tc5AZuvpXn+hTEtgY8xtjzHxjzJlADbADqIqFdqzfh63dK4l6CDHGAQdSvOaDxpiFxpiFZWUa21OUoUpRrpdPnz6JXK+bzz++jobWYMZee/uh6ESzmaOLMvaaTqSvVUCjrN8TgMuBp4BlwA3WLjcAz1rby4DrrWqgxUBdLFSkKIo9GVeSz/2fXMD2qgZufmJdfC2GvrKtqgGAGWNUAPpCX4tn/yQiW4DngFuMMceAu4GPicgO4GPW3wAvADuBCuDXwBf7+N6KogwBzppRxj1L5/FGxVG+/of1GakK2lHVwIgCH6XWok5K7+jTgjDGmDNSjB0FPpJi3AC39OX9FEUZmixdMI6qhlZ+/OI2RhfncOdFs/v0etuqGjT+nwEcvSKYoigDxxfOmkpVXSu/fm0Xo4tz+dwZU3r1Ov5QmB1VjSydr0WEfUUFQFGUAUFE+M7H51Dd6Oe//vY+o4tz+fiJx/X4dX784jYa/SE+cvzofrDSWWgDDUVRBgy3S/jZVSexYGIJ33l2E63BcI+e/8q2w/zm9V1cf9pEzpyhVYJ9RQVAUZQBJdfr5msfm8Gx5iAvbOx+IWB1g5+v/2E9M0cX8a0Lj+9HC52DCoCiKAPOaVNHMqWsgMe7OUs4EjF87Q/raWgN8YtrTybX6+5nC52BCoCiKAOOiPDJUyfy7t5aNu2v63L/h9/Yxcrt1Xz74tnM0MlfGUMFQFGUQWHpgnHkel08/PqulBPEjjUFANi0v457XtzKubNH88lTJwy0mbZGq4AURRkUhuV5ufTEcn6/dh/PbzzII586hdW7amgNhtl2qIFXt1dzzaIJvFFxhJEFOdyzdB4iqVqKKb1FBUBRlEHj306byO/X7iMQinDtQ28BIBIVh8VTRvDUmr2UD8/jl9eeTEmBb5CttR8qAIqiDBonlA/jtf84h+I8L7/45w5On1bKnPJi8n0e8r1udh5pZEppIS6X3vn3ByoAiqIMKuNH5APw7Ys7toeYNkoTvv2JJoEVRVEcigqAoiiKQ1EBUBRFcSgqAIqiKA5FBUBRFMWhqAAoiqI4FBUARVEUh6ICoCiK4lAkulRvdiIi1UD3+sWmphQ4kiFz+hO1M/MMFVvVzswzVGztTzsnGmO6XDEnqwWgr4jIWmPMwsG2oyvUzswzVGxVOzPPULE1G+zUEJCiKIpDUQFQFEVxKHYXgAcH24BuonZmnqFiq9qZeYaKrYNup61zAIqiKEp67O4BKIqiKGkYUgIgIjNF5L2En3oR+YqI/EREtorIBhH5i4gMT3jOHSJSISLbROS8hPHzrbEKEbl9gOz8oWXjeyLysogcZ+0vIvJzy5YNIjI/4bVuEJEd1s8NmbSzM1sTHv+6iBgRKR1MWzs5pt8Tkf0J4xcmPCdrzr312G3W+24WkR8Ppp2d2Soiv08Y2y0i7w2mrZ3YeZKIrLbG1orIImv/bPuMnigib4rIRhF5TkSKE54zKOc+jjFmSP4AbuAQMBE4F/BY4/cA91jbs4H1QA4wGfjAep7b2p4C+Kx9Zg+AncUJ418CHrC2LwT+DgiwGHjLGh8B7LR+l1jbJQNxTK2/xwMvEZ2LUZottiYd0+8BX0+xT7ad+3OAfwA51mOjssXOVOc+Yfy/ge9ki61Jx/Rl4IKEz+UrWfoZfRs4yxr/DPDDbDmeQ8oDSOIjwAfGmD3GmJeNMSFrfDUwztq+FHjaGOM3xuwCKoBF1k+FMWanMSYAPG3t29921ieMFwCxBMylwGMmympguIiMBc4Dlhtjaowxx4DlwPn9ZGc7W62//wf4jwQ7s8XWZDtTkVXnHvgCcLcxxg9gjDmcRXYm2wpE76SBq4CnssjWRDsNELubHgYcSLAzmz6jM4GV1vhyYGmCnYN6PIeyAFxN2wczkc8QVX+AcmBfwmOV1li68f6gnZ0icpeI7AOuA76TRXa2s1VELgH2G2PWJ+2TDbYmn/tbLVf/YREpyVI7ZwBniMhbIvKqiJySRXYm2xrjDKDKGLPD+jsbbE208yvAT6zv00+BO7LUzk3AJdb2lUQ966ywc0gKgIj4iB7QPySN3wmEgN/FhlI83XQynlFS2WmMudMYM96y8dZssBPa2yoi+cCdtAlUu13T2DRYx/R+YCpwEnCQaMiCTuwZLDs9RMMOi4FvAM9Yd9hZde6THrqG9qKQbcf0C8BXre/TV4HfZKmdnwFuEZF1QBEQyAY7YYgKAHAB8I4xpio2YCV0LgauM1aAjahyjk943jiibmK68X63M4EnaXMFB9tOaG/rVKIxyfUistt633dEZEwW2NrumBpjqowxYWNMBPg1UfeZbLPTet8/W2GJNUCEaC+YwbYzla2IiAe4HPh9wn6DbWuynTcAf7a2/0CWnntjzFZjzLnGmAVEBfWDLLFzaCaBicbEPp3w9/nAFqAsab85tE+y7CSaYPFY25NpS7LMGQA7pyds3wb80dq+iPZJqzWmLWm1i+idY4m1PWIgjmnSY7tpSwIPqq0pjunYhO2vEo2pZuO5vxn4gbU9g6iLL4NtZ7pzb32nXk0ay7Zj+j5wtrX9EWBdln5GYwl/F/AY8JlsOJ7GmKEnAEA+cBQYljBWYX2h3rN+Hkh47E6iirsNq2LAGr8Q2G49ducA2fknovHADcBzQLk1LsC9li0bgYUJz/mM9f9VJH9J+9PWpMd30yYAg2ZrmmP6uGXHBmAZ7QUhm869D3jCOv/vAB8ebDs7O/fAI8DNKfbPpmO6BFhnXSDfAhZk6Wf0y9ax2Q7cjTUBd7DPvTFGZwIriqI4laGaA1AURVH6iAqAoiiKQ1EBUBRFcSgqAIqiKA5FBUBRFMWhqAAoiqI4FBUARVEUh6ICoCiK4lD+f2wEwdB0Q8bWAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(dataResam[3][0][0][20*360:22*360],dataResam[3][0][1][20*360:22*360])" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T15:00:42.808786Z", "start_time": "2019-08-18T15:00:42.717815Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# ORIGINAL METHOD\n", "def time_to_sample_number(seconds, frequency):\n", " return seconds * frequency + 0.5\n", "def relu(num):\n", " if num<0:\n", " return 0\n", " return num\n", "\n", "class GQRSorig(object):\n", " \"\"\"\n", " GQRS detection class\n", " \"\"\"\n", " class Conf(object):\n", " \"\"\"\n", " Initial signal configuration object for this qrs detector\n", " \"\"\"\n", " #adc gain is the lvl/mV, can be found in the .hea file, default is 200\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " def __init__(self, fs = 360, adc_gain = 200, hr=75,\n", " RRdelta=0.2, RRmin=0.28, RRmax=2.4,\n", " QS=0.07, QT=0.35,\n", " RTmin=0.25, RTmax=0.33,\n", " QRSa=750, QRSamin=150,\n", " thresh=1.0):\n", " self.fs = fs\n", "\n", " self.sps = int(time_to_sample_number(1, fs))\n", " self.spm = int(time_to_sample_number(60, fs))\n", "\n", " self.hr = hr\n", " self.RR = 60.0 / self.hr\n", " self.RRdelta = RRdelta\n", " self.RRmin = RRmin\n", " self.RRmax = RRmax\n", " self.QS = QS\n", " self.QT = QT\n", " self.RTmin = RTmin\n", " self.RTmax = RTmax\n", " self.QRSa = QRSa\n", " self.QRSamin = QRSamin\n", " self.thresh = thresh\n", "\n", " self._NORMAL = 1 # normal beat\n", " self._ARFCT = 16 # isolated QRS-like artifact\n", " self._NOTE = 22 # comment annotation\n", " self._TWAVE = 27 # T-wave peak\n", " self._NPEAKS = 10 # number of peaks buffered (per signal)\n", " self._BUFLN = 32768 # must be a power of 2, see qf()\n", "\n", " self.rrmean = int(self.RR * self.sps)\n", " self.rrdev = int(self.RRdelta * self.sps)\n", " self.rrmin = int(self.RRmin * self.sps)\n", " self.rrmax = int(self.RRmax * self.sps)\n", "\n", " self.rrinc = int(self.rrmean / 40)\n", " if self.rrinc < 1:\n", " self.rrinc = 1\n", "\n", " self.dt = int(self.QS * self.sps / 4)\n", " if self.dt < 1:\n", " raise Exception('Sampling rate is too low. Unable to use signal.')\n", "\n", " self.rtmin = int(self.RTmin * self.sps)\n", " self.rtmean = int(0.75 * self.QT * self.sps)\n", " self.rtmax = int(self.RTmax * self.sps)\n", " \n", " #dv: min peak to peak val for QRS complex in uVolts\n", " #QRSamin: min peak to peak val for QRS complex in samples\n", " dv = adc_gain * self.QRSamin * 0.001\n", " self.pthr = int((self.thresh * dv * dv) / 6)\n", " self.qthr = self.pthr << 1\n", " self.pthmin = self.pthr >> 2\n", " self.qthmin = int((self.pthmin << 2) / 3)\n", " self.tamean = self.qthr # initial value for mean T-wave amplitude\n", "\n", " # Filter constants and thresholds.\n", " self.dt2 = 2 * self.dt\n", " self.dt3 = 3 * self.dt\n", " self.dt4 = 4 * self.dt\n", "\n", " self.smdt = self.dt\n", " self.v1norm = self.smdt * self.dt * 64\n", "\n", " self.smt = 0\n", " self.smt0 = 0 + self.smdt\n", "\n", " class Peak(object):\n", " def __init__(self, peak_time, peak_amp, peak_type):\n", " self.time = peak_time\n", " self.amp = peak_amp\n", " self.type = peak_type\n", " self.next_peak = None\n", " self.prev_peak = None\n", " def __str__(self):\n", " return \"Peak object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nAmp: \"+str(self.amp)+\\\n", " \"\\nTyppe: \"+str(self.type)\n", "\n", " class Annotation(object):\n", " def __init__(self, ann_time, ann_type, ann_subtype, ann_num):\n", " self.time = ann_time\n", " self.type = ann_type\n", " self.subtype = ann_subtype\n", " self.num = ann_num\n", " def __str__(self):\n", " return \"Annotation object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nNum: \"+str(self.num)+\\\n", " \"\\nTyppe: \"+str(self.type)\n", "\n", "\n", " def putann(self, annotation):\n", " self.annotations.append(copy.deepcopy(annotation))\n", " \n", " def pad(self,lenFinal,signal):\n", " s = list(copy.deepcopy(signal))\n", " delta = lenFinal-len(s)\n", " if delta>0:\n", " pad = list(np.zeros(delta//2))\n", " if delta%2 == 0:\n", " s = pad+s+pad\n", " else:\n", " s = pad+s+pad+[0.0]\n", " else:\n", " delta = -delta\n", " rem = delta//2\n", " s = s[rem:rem+lenFinal]\n", " return s\n", " \n", "\n", " #ENTRY POINT, NEED A CONF OBJ FOR INSIDE THE GQRS CLASS\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " def detect(self, x, conf, adc_zero = 1024):\n", " \"\"\"\n", " Run detection. x is digital signal\n", " \"\"\"\n", " self.c = conf\n", " self.annotations = []\n", " self.sample_valid = False\n", " self.outTimes = [] # self.outTimes[i] = [0,0,T] for compatibility\n", "\n", " if len(x) < 1:\n", " return []\n", "\n", " self.x = x\n", " self.adc_zero = adc_zero\n", "\n", " self.qfv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.smv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.v1 = 0\n", " \n", "\n", " t0 = 0\n", " self.tf = len(x) - 1\n", " self.t = 0 - self.c.dt4\n", "\n", " self.annot = GQRSorig.Annotation(-self.c.rrmin, \"NOTE\", 0, 0)\n", " # Cicular buffer of Peaks\n", " first_peak = GQRSorig.Peak(0, 0, 0)\n", " tmp = first_peak\n", " for _ in range(1, self.c._NPEAKS):\n", " tmp.next_peak = GQRSorig.Peak(0, 0, 0)\n", " tmp.next_peak.prev_peak = tmp\n", " tmp = tmp.next_peak\n", " tmp.next_peak = first_peak\n", " first_peak.prev_peak = tmp\n", " self.current_peak = first_peak\n", "\n", " if self.c.spm > self.c._BUFLN:\n", " if self.tf - t0 > self.c._BUFLN:\n", " tf_learn = t0 + self.c._BUFLN - self.c.dt4\n", " else:\n", " tf_learn = self.tf - self.c.dt4\n", " else:\n", " if self.tf - t0 > self.c.spm:\n", " tf_learn = t0 + self.c.spm - self.c.dt4\n", " else:\n", " tf_learn = self.tf - self.c.dt4\n", "\n", "# self.countdown = -1\n", "# self.state = \"LEARNING\"\n", "# self.gqrs(t0, tf_learn)\n", "\n", " self.rewind_gqrs()\n", "\n", " self.state = \"RUNNING\"\n", " self.t = t0 - self.c.dt4\n", " self.gqrs(t0, self.tf)\n", " return self.outTimes\n", "\n", " def rewind_gqrs(self):\n", " self.countdown = -1\n", " self.at(self.t)\n", " self.annot.time = -self.c.rrmin\n", " self.annot.type = \"NORMAL\"\n", " self.annot.subtype = 0\n", " self.annot.num = 0\n", " p = self.current_peak\n", " for _ in range(self.c._NPEAKS):\n", " p.time = 0\n", " p.type = 0\n", " p.amp = 0\n", " p = p.next_peak\n", "\n", " def at(self, t):\n", " if t < 0:\n", " self.sample_valid = True\n", " return self.x[0]\n", " if t > len(self.x) - 1:\n", " self.sample_valid = False\n", " return self.x[-1]\n", " self.sample_valid = True\n", " return self.x[t]\n", "\n", " #Because _BUFLN is a power of 2 (2**N), _BUFLN-1 is (in binary)\n", " #all 1 from 0 to N-1, hence, self.smv[t & (self.c._BUFLN - 1)] \n", " #is just a fancy pants method for self.smv[t%_BUFLN] \n", " def smv_at(self, t):\n", " return self.smv[t & (self.c._BUFLN - 1)]\n", "\n", " def smv_put(self, t, v):\n", " self.smv[t & (self.c._BUFLN - 1)] = v\n", "\n", " def qfv_at(self, t):\n", " return self.qfv[t & (self.c._BUFLN - 1)]\n", "\n", " def qfv_put(self, t, v):\n", " self.qfv[t & (self.c._BUFLN - 1)] = v\n", " \n", "\n", "\n", " def sm(self, at_t):\n", " # implements a trapezoidal low pass (smoothing) filter\n", " # (with a gain of 4*smdt) applied to input signal sig\n", " # before the QRS matched filter qf().\n", " # Before attempting to 'rewind' by more than BUFLN-smdt\n", " # samples, reset smt and smt0.\n", "\n", " # Calculate samp values from self.smt to at_t.\n", " smt = self.c.smt\n", " smdt = int(self.c.smdt)\n", "\n", " v = 0\n", " while at_t > smt:\n", " smt += 1\n", " # from dt+1 onwards\n", " if smt > int(self.c.smt0):\n", " tmp = int(self.smv_at(smt - 1) + \\\n", " self.at(smt + smdt) + self.at(smt + smdt - 1) - \\\n", " self.at(smt - smdt) - self.at(smt - smdt - 1))\n", "# print('Pos: ',self.at(smt + smdt), self.at(smt + smdt - 1), \\\n", "# 'Neg: ',self.at(smt - smdt), self.at(smt - smdt - 1), \\\n", "# 'Pos-Neg: ',self.at(smt + smdt) + self.at(smt + smdt - 1) - \\\n", "# self.at(smt - smdt) - self.at(smt - smdt - 1))\n", " \n", " self.smv_put(smt, tmp)\n", " self.SIG_SMOOTH.append(tmp)\n", " # from 1 to dt. 0 is never calculated.\n", " else:\n", " v = int(self.at(smt))\n", " for j in range(1, smdt):\n", " smtpj = self.at(smt + j)\n", " smtlj = self.at(smt - j)\n", " v += int(smtpj + smtlj)\n", " self.smv_put(smt, (v << 1) + self.at(smt + j+1) + self.at(smt - j-1) - \\\n", " self.adc_zero * (smdt << 2))\n", "\n", " self.SIG_SMOOTH.append((v << 1) + self.at(smt + j+1) + self.at(smt - j-1) - \\\n", " self.adc_zero * (smdt << 2))\n", " self.c.smt = smt\n", "\n", " return self.smv_at(at_t)\n", "\n", " def qf(self):\n", " # evaluate the QRS detector filter for the next sample\n", "\n", " # do this first, to ensure that all of the other smoothed values needed below are in the buffer\n", " \n", " \n", " dv2 = self.sm(self.t + self.c.dt4)\n", " dv2 -= self.smv_at(self.t - self.c.dt4)\n", " dv1 = int(self.smv_at(self.t + self.c.dt) - self.smv_at(self.t - self.c.dt))\n", " dv = dv1 << 1\n", " dv -= int(self.smv_at(self.t + self.c.dt2) - self.smv_at(self.t - self.c.dt2))\n", " dv = dv << 1\n", " dv += dv1\n", " dv -= int(self.smv_at(self.t + self.c.dt3) - self.smv_at(self.t - self.c.dt3))\n", " dv = dv << 1\n", " dv += dv2\n", " self.v1 += dv\n", " v0 = int(self.v1 / self.c.v1norm)\n", " self.qfv_put(self.t, v0 * v0)\n", " \n", " self.SIG_QRS.append(v0 ** 2)\n", "# #DEBUG \n", "# temp = []\n", "# temp.append(self.smv_at(self.t - self.c.dt4))\n", "# temp.append(self.smv_at(self.t - self.c.dt3))\n", "# temp.append(self.smv_at(self.t - self.c.dt2))\n", "# temp.append(self.smv_at(self.t - self.c.dt))\n", "# temp.append(self.smv_at(self.t))\n", "# temp.append(self.smv_at(self.t + self.c.dt))\n", "# temp.append(self.smv_at(self.t + self.c.dt2))\n", "# temp.append(self.smv_at(self.t + self.c.dt3))\n", "# temp.append(self.smv_at(self.t + self.c.dt4))\n", "# print('T: ',self.t,' >>>>',temp,' ----> ', dv,' (', self.v1,')')\n", "# #DEBUG\n", "\n", "\n", " def gqrs(self, from_sample, to_sample):\n", " q0 = None\n", " q1 = 0\n", " q2 = 0\n", " rr = None\n", " rrd = None\n", " rt = None\n", " rtd = None\n", " rtdmin = None\n", "\n", " p = None # (Peak)\n", " q = None # (Peak)\n", " r = None # (Peak)\n", " tw = None # (Peak)\n", "\n", " last_peak = from_sample\n", " last_qrs = from_sample\n", "\n", " self.SIG_SMOOTH = []\n", " self.SIG_QRS = []\n", "\n", " def add_peak(peak_time, peak_amp, type):\n", " p = self.current_peak.next_peak\n", " p.time = peak_time\n", " p.amp = peak_amp\n", " p.type = type\n", " self.current_peak = p\n", " p.next_peak.amp = 0\n", "\n", " def peaktype(p):\n", " # peaktype() returns 1 if p is the most prominent peak in its neighborhood, 2\n", " # otherwise. The neighborhood consists of all other peaks within rrmin.\n", " # Normally, \"most prominent\" is equivalent to \"largest in amplitude\", but this\n", " # is not always true. For example, consider three consecutive peaks a, b, c\n", " # such that a and b share a neighborhood, b and c share a neighborhood, but a\n", " # and c do not; and suppose that amp(a) > amp(b) > amp(c). In this case, if\n", " # there are no other peaks, a is the most prominent peak in the (a, b)\n", " # neighborhood. Since b is thus identified as a non-prominent peak, c becomes\n", " # the most prominent peak in the (b, c) neighborhood. This is necessary to\n", " # permit detection of low-amplitude beats that closely precede or follow beats\n", " # with large secondary peaks (as, for example, in R-on-T PVCs).\n", " if p.type:\n", " return p.type\n", " else:\n", " a = p.amp\n", " t0 = p.time - self.c.rrmin\n", " t1 = p.time + self.c.rrmin\n", "\n", " if t0 < 0:\n", " t0 = 0\n", "\n", " pp = p.prev_peak\n", " while t0 < pp.time and pp.time < pp.next_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.prev_peak\n", "\n", " pp = p.next_peak\n", " while pp.time < t1 and pp.time > pp.prev_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.next_peak\n", "\n", " p.type = 1\n", " return p.type\n", "\n", " def find_missing(r, p):\n", " if r is None or p is None:\n", " return None\n", "\n", " minrrerr = p.time - r.time\n", "\n", " s = None\n", " q = r.next_peak\n", " while q.time < p.time:\n", " if peaktype(q) == 1:\n", " rrtmp = q.time - r.time\n", " rrerr = rrtmp - self.c.rrmean\n", " if rrerr < 0:\n", " rrerr = -rrerr\n", " if rrerr < minrrerr:\n", " minrrerr = rrerr\n", " s = q\n", " # end:\n", " q = q.next_peak\n", "\n", " return s\n", "\n", " r = None\n", " next_minute = 0\n", " minutes = 0\n", " while self.t <= to_sample + self.c.sps:\n", "# print('-'*40)\n", "# print(self.t)\n", " if self.countdown < 0:\n", " if self.sample_valid:\n", " self.qf()\n", " else:\n", " self.countdown = int(time_to_sample_number(1, self.c.fs))\n", " self.state = \"CLEANUP\"\n", " else:\n", " self.countdown -= 1\n", " if self.countdown < 0:\n", " break\n", " #print(self.t,self.countdown,self.sample_valid)\n", "\n", " q0 = self.qfv_at(self.t)\n", " q1 = self.qfv_at(self.t - 1)\n", " q2 = self.qfv_at(self.t - 2)\n", " #print(q0,q1,q2)\n", " # state == RUNNING only\n", " \n", " if q1 > self.c.pthr and q2 < q1 and q1 >= q0 and self.t > self.c.dt4:\n", " \n", " \n", " add_peak(self.t - 1, q1, 0)\n", " last_peak = self.t - 1\n", " p = self.current_peak.next_peak\n", "# if debug:\n", "# print(\"\\t\\tPEAK\")\n", "# print(\"TIME: \",self.t-1)\n", "# print(\"THRESHOLD: \",self.c.pthr)\n", "# print(\"VALUE: \",q1)\n", "# print(\"TIME NEXT: \",p.time ,' LIMIT TIME: ',self.t - self.c.rtmax)\n", " \n", " while p.time < self.t - self.c.rtmax:\n", "# if debug:\n", "# print(p.time,'\\t' ,self.annot.time , self.c.rrmin,'\\t',p.amp,' ---> ',self.c.qthr)\n", " if p.time >= self.annot.time + self.c.rrmin and peaktype(p) == 1:\n", " if p.amp > self.c.qthr:\n", "# if debug:\n", "# print(\"p.amp:\",p.amp,\" > self.c.qthr:\",self.c.qthr) \n", " rr = p.time - relu(self.annot.time)\n", " q = find_missing(r, p)\n", " if rr > self.c.rrmean + 2 * self.c.rrdev and \\\n", " rr > 2 * (self.c.rrmean - self.c.rrdev) and \\\n", " q is not None:\n", " p = q\n", " rr = p.time - self.annot.time\n", " self.annot.subtype = 1\n", " rrd = rr - self.c.rrmean\n", " if rrd < 0:\n", " rrd = -rrd\n", " self.c.rrdev += (rrd - self.c.rrdev) >> 3\n", " if rrd > self.c.rrinc:\n", " rrd = self.c.rrinc\n", " if rr > self.c.rrmean:\n", " self.c.rrmean += rrd\n", " else:\n", " self.c.rrmean -= rrd\n", " if p.amp > self.c.qthr * 4:\n", " self.c.qthr += 1\n", " elif p.amp < self.c.qthr:\n", " self.c.qthr -= 1\n", " if self.c.qthr > self.c.pthr * 20:\n", " self.c.qthr = self.c.pthr * 20\n", " last_qrs = p.time\n", "\n", " if self.state == \"RUNNING\":\n", " self.annot.time = p.time - self.c.dt2\n", " self.annot.type = \"NORMAL\"\n", " qsize = int(p.amp * 10.0 / self.c.qthr)\n", " if qsize > 127:\n", " qsize = 127\n", " self.annot.num = qsize\n", " self.putann(self.annot)\n", "###############################################################################################################\n", " ##OUTPUT HERE##\n", " self.outTimes.append([0,0,self.annot.time])\n", "###############################################################################################################\n", " self.annot.time += self.c.dt2\n", "\n", " #Look for this beat's T-wave\n", " tw = None\n", " rtdmin = self.c.rtmean\n", " q = p.next_peak\n", " while q.time > self.annot.time:\n", " rt = q.time - self.annot.time - self.c.dt2\n", " if rt < self.c.rtmin:\n", " # end:\n", " q = q.next_peak\n", " continue\n", " if rt > self.c.rtmax:\n", " break\n", " rtd = rt - self.c.rtmean\n", " if rtd < 0:\n", " rtd = -rtd\n", " if rtd < rtdmin:\n", " rtdmin = rtd\n", " tw = q\n", " # end:\n", " q = q.next_peak\n", " if tw is not None:\n", " tmp_time = tw.time - self.c.dt2\n", " tann = GQRSorig.Annotation(tmp_time, \"TWAVE\",\n", " 1 if tmp_time > self.annot.time + self.c.rtmean else 0,\n", " rtdmin)\n", " # if self.state == \"RUNNING\":\n", " # self.putann(tann)\n", " rt = tann.time - self.annot.time\n", " self.c.rtmean += (rt - self.c.rtmean) >> 4\n", " if self.c.rtmean > self.c.rtmax:\n", " self.c.rtmean = self.c.rtmax\n", " elif self.c.rtmean < self.c.rtmin:\n", " self.c.rtmean = self.c.rrmin\n", " tw.type = 2 # mark T-wave as secondary\n", " r = p\n", " q = None\n", " self.annot.subtype = 0\n", " elif self.t - last_qrs > self.c.rrmax and self.c.qthr > self.c.qthmin:\n", " self.c.qthr -= (self.c.qthr >> 4)\n", " # end:\n", " p = p.next_peak\n", " \n", " elif self.t - last_peak > self.c.rrmax and self.c.pthr > self.c.pthmin:\n", " self.c.pthr -= (self.c.pthr >> 4)\n", "\n", " self.t += 1\n", " if self.t >= next_minute:\n", " next_minute += self.c.spm\n", " minutes += 1\n", " if minutes >= 60:\n", " minutes = 0\n", "\n", " if self.state == \"LEARNING\":\n", " return\n", "\n", " # Mark the last beat or two.\n", " p = self.current_peak.next_peak\n", " i = 0\n", " while p.time <= p.next_peak.time:\n", " if p.time >= self.annot.time + self.c.rrmin and p.time < self.tf and peaktype(p) == 1:\n", " self.annot.type = \"NORMAL\"\n", " self.annot.time = p.time\n", " self.putann(self.annot)\n", "###############################################################################################################\n", " ##OUTPUT HERE##\n", " self.outTimes.append([0,0,self.annot.time])\n", "###############################################################################################################\n", " # end:\n", " p = p.next_peak\n", " if i>= self.c._NPEAKS:\n", " break\n", " i+=1" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T15:32:58.028800Z", "start_time": "2019-08-18T15:11:06.465153Z" }, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--------- 1 / 4\n", "1 / 48 ---> 6.890849590301514\n", "2 / 48 ---> 6.863418102264404\n", "3 / 48 ---> 6.817744493484497\n", "4 / 48 ---> 6.838284969329834\n", "5 / 48 ---> 6.903501033782959\n", "6 / 48 ---> 6.895578861236572\n", "7 / 48 ---> 6.855767488479614\n", "8 / 48 ---> 6.994442462921143\n", "9 / 48 ---> 6.927514314651489\n", "10 / 48 ---> 6.795306205749512\n", "11 / 48 ---> 6.865689516067505\n", "12 / 48 ---> 6.827620267868042\n", "13 / 48 ---> 6.76933479309082\n", "14 / 48 ---> 6.846712589263916\n", "15 / 48 ---> 6.794588327407837\n", "16 / 48 ---> 6.822516918182373\n", "17 / 48 ---> 6.825845003128052\n", "18 / 48 ---> 6.697320461273193\n", "19 / 48 ---> 6.845311880111694\n", "20 / 48 ---> 6.808375597000122\n", "21 / 48 ---> 6.8282880783081055\n", "22 / 48 ---> 7.032347917556763\n", "23 / 48 ---> 6.869121074676514\n", "24 / 48 ---> 6.766818046569824\n", "25 / 48 ---> 6.8810715675354\n", "26 / 48 ---> 6.9044458866119385\n", "27 / 48 ---> 6.894068956375122\n", "28 / 48 ---> 6.85528826713562\n", "29 / 48 ---> 6.864752292633057\n", "30 / 48 ---> 6.868131875991821\n", "31 / 48 ---> 6.816994667053223\n", "32 / 48 ---> 6.791010856628418\n", "33 / 48 ---> 6.8068931102752686\n", "34 / 48 ---> 6.835047960281372\n", "35 / 48 ---> 6.9153571128845215\n", "36 / 48 ---> 7.119139671325684\n", "37 / 48 ---> 6.811431646347046\n", "38 / 48 ---> 6.676790714263916\n", "39 / 48 ---> 6.809314489364624\n", "40 / 48 ---> 6.793010711669922\n", "41 / 48 ---> 6.808601140975952\n", "42 / 48 ---> 6.869395017623901\n", "43 / 48 ---> 6.811517715454102\n", "44 / 48 ---> 7.025799512863159\n", "45 / 48 ---> 6.8990068435668945\n", "46 / 48 ---> 6.9022862911224365\n", "47 / 48 ---> 6.7199695110321045\n", "48 / 48 ---> 6.832911252975464\n", "--------- 2 / 4\n", "1 / 48 ---> 6.830714464187622\n", "2 / 48 ---> 6.855673789978027\n", "3 / 48 ---> 6.830801010131836\n", "4 / 48 ---> 7.065879821777344\n", "5 / 48 ---> 6.799373626708984\n", "6 / 48 ---> 6.733932256698608\n", "7 / 48 ---> 6.69924259185791\n", "8 / 48 ---> 6.745116233825684\n", "9 / 48 ---> 6.793877363204956\n", "10 / 48 ---> 6.7591211795806885\n", "11 / 48 ---> 6.731123208999634\n", "12 / 48 ---> 6.7189621925354\n", "13 / 48 ---> 7.100652456283569\n", "14 / 48 ---> 6.823315620422363\n", "15 / 48 ---> 6.780345916748047\n", "16 / 48 ---> 6.823434352874756\n", "17 / 48 ---> 6.8489532470703125\n", "18 / 48 ---> 6.760432243347168\n", "19 / 48 ---> 6.797024488449097\n", "20 / 48 ---> 6.749071836471558\n", "21 / 48 ---> 6.8140037059783936\n", "22 / 48 ---> 6.812416076660156\n", "23 / 48 ---> 7.0913848876953125\n", "24 / 48 ---> 6.772224187850952\n", "25 / 48 ---> 6.7818074226379395\n", "26 / 48 ---> 6.8762853145599365\n", "27 / 48 ---> 6.811401844024658\n", "28 / 48 ---> 6.854294538497925\n", "29 / 48 ---> 6.853086471557617\n", "30 / 48 ---> 6.861169338226318\n", "31 / 48 ---> 6.8596110343933105\n", "32 / 48 ---> 6.8102147579193115\n", "33 / 48 ---> 6.816570520401001\n", "34 / 48 ---> 6.806317567825317\n", "35 / 48 ---> 6.902137279510498\n", "36 / 48 ---> 6.899842739105225\n", "37 / 48 ---> 6.812600374221802\n", "38 / 48 ---> 6.811528921127319\n", "39 / 48 ---> 6.927834987640381\n", "40 / 48 ---> 6.9459686279296875\n", "41 / 48 ---> 6.959559679031372\n", "42 / 48 ---> 6.865916013717651\n", "43 / 48 ---> 7.055027723312378\n", "44 / 48 ---> 6.833699941635132\n", "45 / 48 ---> 7.105437278747559\n", "46 / 48 ---> 7.131114959716797\n", "47 / 48 ---> 7.099620342254639\n", "48 / 48 ---> 7.163788557052612\n", "--------- 3 / 4\n", "1 / 48 ---> 7.0224714279174805\n", "2 / 48 ---> 6.863124847412109\n", "3 / 48 ---> 6.8707451820373535\n", "4 / 48 ---> 7.027323484420776\n", "5 / 48 ---> 6.8426642417907715\n", "6 / 48 ---> 6.85205078125\n", "7 / 48 ---> 6.948765754699707\n", "8 / 48 ---> 6.8398683071136475\n", "9 / 48 ---> 6.867389917373657\n", "10 / 48 ---> 6.796664476394653\n", "11 / 48 ---> 6.775071620941162\n", "12 / 48 ---> 6.838229417800903\n", "13 / 48 ---> 6.7921836376190186\n", "14 / 48 ---> 6.774097204208374\n", "15 / 48 ---> 6.750752925872803\n", "16 / 48 ---> 6.828451871871948\n", "17 / 48 ---> 6.8614091873168945\n", "18 / 48 ---> 7.185412645339966\n", "19 / 48 ---> 6.818621635437012\n", "20 / 48 ---> 6.81509256362915\n", "21 / 48 ---> 6.817326545715332\n", "22 / 48 ---> 6.876378297805786\n", "23 / 48 ---> 6.819307565689087\n", "24 / 48 ---> 6.750303268432617\n", "25 / 48 ---> 6.798572540283203\n", "26 / 48 ---> 6.916138410568237\n", "27 / 48 ---> 6.786556243896484\n", "28 / 48 ---> 6.716574192047119\n", "29 / 48 ---> 6.733654022216797\n", "30 / 48 ---> 6.768553972244263\n", "31 / 48 ---> 6.7146806716918945\n", "32 / 48 ---> 6.811197996139526\n", "33 / 48 ---> 6.684150695800781\n", "34 / 48 ---> 6.681123971939087\n", "35 / 48 ---> 6.797106742858887\n", "36 / 48 ---> 6.759219169616699\n", "37 / 48 ---> 6.704828500747681\n", "38 / 48 ---> 6.67307448387146\n", "39 / 48 ---> 6.775623083114624\n", "40 / 48 ---> 6.753901958465576\n", "41 / 48 ---> 6.685752630233765\n", "42 / 48 ---> 6.781758069992065\n", "43 / 48 ---> 6.747058629989624\n", "44 / 48 ---> 6.779052257537842\n", "45 / 48 ---> 7.092763900756836\n", "46 / 48 ---> 6.678934097290039\n", "47 / 48 ---> 6.67268705368042\n", "48 / 48 ---> 6.742933750152588\n", "--------- 4 / 4\n", "1 / 48 ---> 6.711636543273926\n", "2 / 48 ---> 6.736102819442749\n", "3 / 48 ---> 6.716845512390137\n", "4 / 48 ---> 6.749647378921509\n", "5 / 48 ---> 6.7956156730651855\n", "6 / 48 ---> 6.857494592666626\n", "7 / 48 ---> 6.8462724685668945\n", "8 / 48 ---> 6.845669507980347\n", "9 / 48 ---> 6.849987745285034\n", "10 / 48 ---> 6.815703392028809\n", "11 / 48 ---> 6.806895971298218\n", "12 / 48 ---> 6.854324579238892\n", "13 / 48 ---> 6.8294899463653564\n", "14 / 48 ---> 6.815604209899902\n", "15 / 48 ---> 6.774214029312134\n", "16 / 48 ---> 6.702061891555786\n", "17 / 48 ---> 6.738706588745117\n", "18 / 48 ---> 6.688871383666992\n", "19 / 48 ---> 6.6972455978393555\n", "20 / 48 ---> 6.71265721321106\n", "21 / 48 ---> 6.698553085327148\n", "22 / 48 ---> 6.755293130874634\n", "23 / 48 ---> 6.804404258728027\n", "24 / 48 ---> 6.799153566360474\n", "25 / 48 ---> 6.805914402008057\n", "26 / 48 ---> 6.802731275558472\n", "27 / 48 ---> 7.223186016082764\n", "28 / 48 ---> 6.8150670528411865\n", "29 / 48 ---> 6.840357542037964\n", "30 / 48 ---> 6.888453245162964\n", "31 / 48 ---> 6.892460107803345\n", "32 / 48 ---> 6.814990997314453\n", "33 / 48 ---> 6.696547269821167\n", "34 / 48 ---> 6.731197834014893\n", "35 / 48 ---> 6.8276731967926025\n", "36 / 48 ---> 6.728882074356079\n", "37 / 48 ---> 6.727901935577393\n", "38 / 48 ---> 6.704946994781494\n", "39 / 48 ---> 6.782262086868286\n", "40 / 48 ---> 6.747831583023071\n", "41 / 48 ---> 6.785724639892578\n", "42 / 48 ---> 6.9242260456085205\n", "43 / 48 ---> 6.796584367752075\n", "44 / 48 ---> 6.832125902175903\n", "45 / 48 ---> 6.824783802032471\n", "46 / 48 ---> 6.782539367675781\n", "47 / 48 ---> 6.789820194244385\n", "48 / 48 ---> 6.801103830337524\n" ] }, { "ename": "NameError", "evalue": "name 'methodsOut' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 30\u001b[0m \u001b[0mfoundSampleOneComplete\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfoundAllTh\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 31\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfoundSampleOneComplete\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint64\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m \u001b[0mpp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmethodsOut\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdefaultSample\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNameError\u001b[0m: name 'methodsOut' is not defined" ] } ], "source": [ "foundAllTh = []\n", "l = 1\n", "hAllTh = dataResam\n", "for hFiles in hAllTh:\n", " print('---------',l,\"/\",len(thresholdComplete))\n", " l+=1\n", " foundAll = []\n", " n=1\n", "\n", " verbose = False\n", "\n", " for h in hFiles:\n", " startT = time.time()\n", " c = GQRSorig.Conf()\n", " detector = GQRSorig()\n", " \n", " print(n,\"/\",len(hFiles),end = \" \")\n", " n+=1\n", " start = 0\n", " end = h[1][-1]\n", " times,points = h[0],h[1]\n", "\n", " found = detector.detect(points,c,1024)\n", "\n", " stopT = time.time()\n", " print(\" ---> \",stopT-startT)\n", " foundAll.append(found)\n", " foundAllTh.append(foundAll)\n", " \n", "foundSampleOneComplete = foundAllTh[0][0]\n", "t = np.array(np.asarray(foundSampleOneComplete)[:,2], dtype=np.int64)\n", "pp.pprint(methodsOut[0][0][defaultSample])" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T16:37:41.711628Z", "start_time": "2019-08-18T16:37:40.786175Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 / 4\n", "2 / 4\n", "3 / 4\n", "4 / 4\n" ] }, { "data": { "text/plain": [ "207 1 0 0 0 None" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "toSaveAllTh = []\n", "l = 1\n", "for allFileOneTh in foundAllTh:\n", " print(l,\"/\",len(foundAllTh))\n", " l+=1\n", " toSave = []\n", " for sample in allFileOneTh:\n", " sampleToSave = []\n", " for elem in sample:\n", " sampleToSave.append(MITA.MITAnnotation(time = elem[2]))\n", " toSave.append(sampleToSave)\n", " toSaveAllTh.append(toSave)\n", "toSaveAllTh[0][0][0]" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "ExecuteTime": { "end_time": "2019-08-18T16:39:21.244373Z", "start_time": "2019-08-18T16:39:19.561083Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---------------- 1 / 4\n", "4\n", "---------------- EVALUATING ----------------\n", "Executed: /home/zanoli/Desktop/thesis/silvio-zanoli/scripts/eval_mitdb.sh -f val4bits1HystLvlCrossingOrig_gQRS.txt -o out4bits1HystLvlCrossingOrig_gQRS -s /home/zanoli/Desktop/thesis/silvio-zanoli/data/dataRaw/ -d /home/zanoli/Desktop/thesis/silvio-zanoli/data/results/ \n", "Result: 0\n", "---------------- 2 / 4\n", "5\n", "---------------- EVALUATING ----------------\n", "Executed: /home/zanoli/Desktop/thesis/silvio-zanoli/scripts/eval_mitdb.sh -f val5bits1HystLvlCrossingOrig_gQRS.txt -o out5bits1HystLvlCrossingOrig_gQRS -s /home/zanoli/Desktop/thesis/silvio-zanoli/data/dataRaw/ -d /home/zanoli/Desktop/thesis/silvio-zanoli/data/results/ \n", "Result: 0\n", "---------------- 3 / 4\n", "6\n", "---------------- EVALUATING ----------------\n", "Executed: /home/zanoli/Desktop/thesis/silvio-zanoli/scripts/eval_mitdb.sh -f val6bits1HystLvlCrossingOrig_gQRS.txt -o out6bits1HystLvlCrossingOrig_gQRS -s /home/zanoli/Desktop/thesis/silvio-zanoli/data/dataRaw/ -d /home/zanoli/Desktop/thesis/silvio-zanoli/data/results/ \n", "Result: 0\n", "---------------- 4 / 4\n", "7\n", "---------------- EVALUATING ----------------\n", "Executed: /home/zanoli/Desktop/thesis/silvio-zanoli/scripts/eval_mitdb.sh -f val7bits1HystLvlCrossingOrig_gQRS.txt -o out7bits1HystLvlCrossingOrig_gQRS -s /home/zanoli/Desktop/thesis/silvio-zanoli/data/dataRaw/ -d /home/zanoli/Desktop/thesis/silvio-zanoli/data/results/ \n", "Result: 0\n" ] } ], "source": [ "# All the Files, independent number of thresholds\n", "cwd = os.path.dirname(os.getcwd())\n", "sourceDir = cwd+'/data/dataRaw/'\n", "outDir = cwd+'/data/results/'\n", "tempdir = cwd+'/temp/'\n", "l = 0\n", "for th in toSaveAllTh:\n", " print(\"----------------\",l+1,\"/\",len(toSaveAllTh))\n", " print(thresholdComplete[l])\n", " ext = 'out'+str(thresholdComplete[l])+'bits'+str(hyst[0])+'HystLvlCrossingOrig_gQRS'\n", " evalFile = 'val'+str(thresholdComplete[l])+'bits'+str(hyst[0])+'HystLvlCrossingOrig_gQRS.txt'\n", " for sample,name in zip(th,selectedFileNames):\n", " pathOut = sourceDir+name+'.'+ext\n", " #print(\"saving \",pathOut)\n", " MITA.save_annotations(sample,pathOut)\n", " print(\"---------------- EVALUATING ----------------\")\n", " command = cwd+'/scripts/eval_mitdb.sh -f ' + evalFile +' -o ' + ext\n", " command += ' -s '+sourceDir+' -d '+outDir\n", " r = os.system(command)\n", " print('Executed: ',command,'\\nResult: ',r)\n", " l+=1\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.7.1" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }