Kaggle-NDSB-2-MXNet-R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transforming How We Diagnose Heart Disease\n",
    "\n",
    "\n",
    "Cardiovascular disease, also commonly referred to as heart disease, takes on many forms - heart attacks, stroke, heart failure and arrhythmia - to name a few.  According to a report published on behalf of the American Heart Association, 610,000 Americans die each year from heart disease.\n",
    "\n",
    "Booz | Allen | Hamilton and Kaggle offered the [2015 Data Science Bowl](https://www.kaggle.com/c/second-annual-data-science-bowl) as an opportunity for data scientist to help automate the diagnosis of heart disease through a competition.  Part of the initial overview of the challenge included the following image identifying some of the metrics used for determining heart disease.\n",
    "\n",
    "![alt text](img/fig2.jpg \"Only One Unit of Time That Matters\")\n",
    "\n",
    "End-systolic and end-diastolic volumes yield the ejection fraction (EF) noted above in the image and represented below in the formula where $V_S$ denotes the volumes at systole and $V_D$ denotes the volumes at diastole.  Volumes and EF can be used as an indication of the presence of heart disease.\n",
    "\n",
    "$$100\\cdot\\frac{V_D-V_S}{V_D}$$  \n",
    "\n",
    "In addition, Magnetic Resonance Imaging (MRI) correctly assess the heart's squeezing ability.  However, evaluation of MRIs is a manual, time consuming process - deep neural networks can automate this process.  Variations in anatomy, function, image quality, and acquisition make automated quantification of left ventricle size a challenging problem. You will encounter this variation in the competition dataset, which aims to provide a diverse representation of cases. \n",
    "\n",
    "<img src=\"img/fig8.png\" alt=\"Manual determination of ejection fraction\" style=\"width: 700px;\"/>\n",
    "<p style=\"text-align: center;\">*C.M.S Nambakhsh et al., Medical Image Analysis 17(2013) 1010-1024*</p>\n",
    "\n",
    "\n",
    "## The Dataset\n",
    "The National Heart, Lung, and Blood Institute (NHLBI) provided the MRI images for the Data Science Bowl. The [dataset](https://www.kaggle.com/c/second-annual-data-science-bowl/data) consists of hundreds of cardiac MRI images in [DICOM](https://en.wikipedia.org/wiki/DICOM) format. These 2D cine images contain approximately 30 images across the cardiac cycle. Each slice is acquired on a separate breath hold. This is important since the registration from slice to slice is expected to be imperfect.\n",
    "\n",
    "\n",
    "<img src=\"img/fig3.gif\" alt=\"MRI imaging of heartbeat\" style=\"width: 500px;\"/>\n",
    "\n",
    "The dataset contains patients from young to old, images from numerous hospitals, and hearts from normal to abnormal cardiac function. A computational method which is robust to these variations could both validate and automate the cardiologists' manual measurement of ejection fraction.  Each case has an associated directory of DICOM files. The exact number of images will differ from case to case, either varying in the number of slices, the views which are captured, or the number of frames in the time sequences.\n",
    "\n",
    "The main view for assessing ventricle size is the short axis stack (PSAX), which contains images taken in a plane perpendicular to the long axis (PLAX) of the left ventricle: \n",
    "\n",
    "<img src=\"img/fig6.png\" alt=\"PSAX diagram\" style=\"width: 250px;\"/>\n",
    "<p style=\"text-align: center;\">*image credit: fpnotebook.com*</p>\n",
    "\n",
    "These have the prefix \"sax_\" in the dataset. Most cases also have alternative views."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation\n",
    "Performance is evaluated using the Continuous Ranked Probability Score (CRPS). For each MRI, a cumulative probability distribution is predicted for both the systolic and diastolic volumes (two separate distributions per case). The CRPS is computed as follows:\n",
    "\n",
    "$$CRPS = \\frac{1}{600\\cdot N}\\sum_{m=1}^{N}\\sum_{n=0}^{599} (P(y\\le n)-H(n-V_{m}))^2$$\n",
    "\n",
    "where $P$ is the predicted distribution, $N$ is the number of rows in the test set (equal to twice the number of cases), $V$ is the actual volume (in mL) and $H(x)$ is the Heaviside step function ($H(x\\lt0)=0$ and $H(x\\ge 0)=1$).  While it is not simple to visualize the CRPS, the shaded area on the figure below may be a helpful guide for understanding the error term between the predicted distribution and actual volume:\n",
    "<img src=\"img/fig4.png\" alt=\"CRPS error between predicted and actual volume\" style=\"width: 500px;\"/>\n",
    "Note that the entry will not score if any of the predicted values has $P(y \\le k) > P(y \\le k+1)$ for any $k$.  That is the CDF is non-decreasing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# End-to-End Deep Learning for NDSB-II\n",
    "\n",
    "In this example, we will show how to use GPU accelerated MXNet library to build an end-to-end deep learning system to help diagnose heart disease.  Keep in mind, this is a very simple model without any network structure optimizations or hyper parameter tuning.  However, it is possible to build fantastic networks based on the example solution.\n",
    "\n",
    "### About MXNet\n",
    "[MXNet](https://github.com/dmlc/mxnet) is a deep learning framework designed for both efficiency and flexibility by DMLC group. MXNet will fully utilize all the resources to solve the problem under limited resource constraint, with a flexible programming interface. You can use it for all purposes of data science and deep learning tasks with R, Julia, Python and more. To run on multiple GPU with huge network, or questions about saving network parameters etc, please refer [MXNet docs](https://mxnet.readthedocs.org/en/latest/)\n",
    "\n",
    "## General Overview of Model\n",
    "### Input Data\n",
    "The dataset itself contains 500 training studies with on average 10 unique SAX observation at various points along the the PLAX of the left ventricle.  That’s roughly 5000 observations in total.  Each of these SAX observation usually contain 30 DICOM images (“frames”) in a time sequence which captures an entire heartbeat (expansion and contraction).   The idea here is to pack each of the 30 image frames of an observation into a 64x64x30 tensor.  While this might sound complicated, just think of stacking playing cards to form a deck.  Each card has the same rectangular dimensions (say, 120x90) and then cards are stacked one on top of the other to form a deck.  It’s the same thing here.  We’re going to take each of the 30 square image frames of 64x64 pixels and stack them one-by-one to form an input “deck” to the deep neural network.  The only difference with the card deck analogy is that the image frames are organized in time so that the video clip of the heartbeat is not all garbled.  This tutorial is based on this simple idea: we first accumulate all suitable observations having 30 frames, then feed to the deep neural network to learn the target directly.\n",
    "\n",
    "The label data set contains only 500 labels. That is, each of the roughly 10 SAX observation for a particular study get the same diastole and systole volume.  Again, each SAX observation of a study is just a different view of the same heart (i.e. a single heartbeat as measured from different cross sections perpendicular to the long axis of the heart).  Therefore, in the data preprocessing step, each label must be duplicated for each unique observation.\n",
    "\n",
    "Additionally, the same 64x64x30 deck of images is used to predict both systole and diastole volumes.  Therefore, we’re going to build a network that ingests each observation tensor to predict systole volume.  Then a separate network that ingests that same image deck to predict diastole volume.\n",
    "\n",
    "Another idea used in this tutorial is taking a frame-by-frame difference to measure change per frame (in time).  By using MXNet symbolic interface, we can dynamically difference the input inside of the network. It helps a little in the final result.\n",
    "\n",
    "### Network Objective\n",
    "For the network, we use the well documented LeNet style convolution network with batch normalization and dropout. This is a basic network with a generic configuration.  In this challenge, we are asked to predict a CDF value of 600 data-point. Therefore the problem is formulated as a regression problem. We ask the neural-net that given a stack of images 64x64x30 to output 600x1 vector - one predicted value for each of the 600 points in the CDF.  Note here that the label is just a single floating point value like 83.3 (mL).  This label is transformed into a step function having 600 discrete values where all the y-values of the CDF with x less than 83 are 0 and all y-values with x greater than 83 are 1.\n",
    "\n",
    "\n",
    "### Preprocessing\n",
    "We first run a preprocessing step, to pack the data into a csv file (```train-64x64-data.csv```).  Each line of this csv file corresponds to a 64 x 64 x 30 tensor, which gives 30 frames of images. We can also use other inputs besides csv. The CSV is used here since this format is common for all languages and is easy to parse.\n",
    "\n",
    "The input dataset is quite big (5293 observations of size 64x64x30). While this data set can likely fit into memory of a big machine, we want to be safe, so we will use the ```CSVIter``` from ```mxnet``` to load data from disk on-the-fly during training, without loading all the data into memory at once.\n",
    "\n",
    "The labels for the training data are stored in ```train-systole.csv``` and ```train-diastole.csv``` where each line is a step function of 0s and 1s as described above. Line ```i``` in the training CSV data file is label by the associated line ```i``` in the label CSV files.  Again, when training the systole network we will use labels from ```train-systole.csv``` and likewise when training the diastole network we will use labels from ```train-diastole.csv```.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we load the necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: mxnet\n",
      "Loading required package: data.table\n",
      "Loading required package: ggplot2\n"
     ]
    }
   ],
   "source": [
    "# import libraries\n",
    "require(mxnet)\n",
    "require(data.table)\n",
    "require(ggplot2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define a function that will architect (not trained yet) neural networks.  This network is a classic network that has it's origins from Yann LeCun's paper many years ago and hence has the nastalgic name \"LeNet\".  \n",
    "\n",
    "Notice here that we create a variable using ```mx.symbol.Variable``` which *represents* the input data.  Just like regular variables we can manipulate the ```source``` variable.  For example, the first thing we do to the image data input is to normalize the pixel values. Next we utilize ```mx.symbol.SliceChannel``` so that we can actually access each individual frame of the input.  Recall that the input is 64x64x30 so that each indexed frame is of size 64x64. The purpose of using the frames is so that differences between successive frames can be calculated. That is, instead of training the network parameters using the images directly, the network is trained on the differences between each successive frame.  This is pretty cool since we can actually manipulate the data within the network definition rather than having to preprocess the whole dataset when we want to experiment with a new idea.  Once each frame delta has been computed the individual frames are put back into a 64x64x29 tensor input using ```mxnet:::mx.varg.symbol.Concat```.  Since we've taken successive deltas we actually only have 29 total deltas now!  \n",
    "\n",
    "From here we just chain various layers together, such as ```mx.symbol.Convolution```, to form the network described in LeCun's paper.  There is no magic here.  Feel free to add layers in between or remove layers etc.  That's part of the fun of working with DL; various recipies have been described in papers over the years but you're free to construct anything you'd like.  Play with the network structure to get a feel for how training performance changes. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create LeNet style network\n",
    "get.lenet <- function() {\n",
    "    \n",
    "    # create data variable (i.e. symbol)\n",
    "    source <- mx.symbol.Variable(\"data\")\n",
    "    \n",
    "    # normalize values of input data \n",
    "    source <- (source-128) / 128\n",
    "    \n",
    "    # SliceChannel is a symbol that can be indexed just like an array\n",
    "    frames <- mx.symbol.SliceChannel(source, num.outputs = 30);\n",
    "    \n",
    "    # init list of differences\n",
    "    diffs <- list()\n",
    "  \n",
    "    # compute differences for each \"frame\"\n",
    "    for (i in 1:29) {\n",
    "        diffs <- c(diffs, frames[[i + 1]] - frames[[i]])\n",
    "    }\n",
    "  \n",
    "    # set the property for number of arugments in diffs variable\n",
    "    diffs$num.args = 29\n",
    "    \n",
    "    # concatonate frame variables into single variable\n",
    "    source <- mxnet:::mx.varg.symbol.Concat(diffs)\n",
    "    \n",
    "    # convolution layer with 5x5 kernel dimention\n",
    "    net <- mx.symbol.Convolution(source, kernel = c(5, 5), num.filter = 40)\n",
    "    \n",
    "    # normalization layer: simply subtract mean divide by std\n",
    "    net <- mx.symbol.BatchNorm(net, fix.gamma = TRUE)\n",
    "    \n",
    "    # activation layer using Rectified Linear Unit (relu) activation function\n",
    "    net <- mx.symbol.Activation(net, act.type = \"relu\")\n",
    "    \n",
    "    # Max pooling layer with a 2x2 kernel and no overlap since stride is 2x2\n",
    "    net <- mx.symbol.Pooling(net, pool.type = \"max\", kernel = c(2, 2), stride = c(2, 2))\n",
    "    \n",
    "    # continue building the network ...\n",
    "    net <- mx.symbol.Convolution(net, kernel = c(3, 3), num.filter = 40)\n",
    "    net <- mx.symbol.BatchNorm(net, fix.gamma = TRUE)\n",
    "    net <- mx.symbol.Activation(net, act.type = \"relu\")\n",
    "    net <- mx.symbol.Pooling(net, pool.type = \"max\", kernel = c(2, 2), stride = c(2, 2))\n",
    "    \n",
    "    # flatten the features to single variable\n",
    "    flatten <- mx.symbol.Flatten(net)\n",
    "    \n",
    "    # add a drop out layer where 50% of data gets dropped out at training time\n",
    "    flatten <- mx.symbol.Dropout(flatten)\n",
    "    \n",
    "    # add the final fully connected layer for output of dimension 600\n",
    "    fc1 <- mx.symbol.FullyConnected(data = flatten, num.hidden = 600)\n",
    "    \n",
    "    # Name the final layer as softmax so it auto matches the naming of data iterator\n",
    "    # Otherwise we can also change the provide_data in the data iter\n",
    "    return(mx.symbol.LogisticRegressionOutput(data = fc1, name = 'softmax'))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to know more about a particular layer, pull up the help listing for that layer. For example,?mx.symbol.Activation provides the help information for the activation layer where additional functions are listed for use such as 'sigmoid', 'tanh' and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "?mx.symbol.Activation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next the training batch size is specified.  The size can be smaller or larger arbitrarily but since we're using GPUs for the training we have to be careful not to use a really large batch size since the GPU has limited memory.  The batch size is also a factor of how large the input images are.  Since each input is a 64x64x30 tensor, 16 is a good size. It is not the end of the world if batch size is too large.  If the batch size is too large, the training command will just fail and complain about something like \"device out of memory\" which is likely your queue to reduce the batch size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "batch_size <- 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we load the training data using an iterator so that we use host memory efficiently.  This is nice since it does not consume huge amounts of data but the trade off is that it does increase the training time since we have to keep loading data over and over. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# CSVIter is used here, since the data can't fit into memory\n",
    "data_train <- mx.io.CSVIter(\n",
    "   data.csv  = \"./train-64x64-data.csv\", data.shape  = c(64, 64, 30),\n",
    "  label.csv  = \"./train-systole.csv\"   , label.shape = 600,\n",
    "  batch.size = batch_size\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets now initialize a LeNet style network using our network function.  This network has structure but none of the parameters are determined yet.  Think of it as raw materials that will be shaped into meaningful values/parameters during the training phase.  Notice that when this function is called all of the parameters are initialized to random values which does actually consume memory on the system.  If the network has 10 million parameters in it then we have to initialize 10 million values in memory.  Also, all the network parameters must actually fit into the GPU memory also (with the input data too!) so if the network structure is too large then it is possible to run out of memory on the GPU.  So start conservatively and grow the network as you work through the performance requirements. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "network <- get.lenet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can train the network parameters we must specify a *loss function* also called a *cost function*.  This function tells the training algoritm how closely the current set of parameters reproduce the target output.  In general, for a given cost function the training algorithm will try to choose parameters for which the cost function is minimized.  That is, a lower cost is better than higher cost. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Custom evaluation metric using CRPS.\n",
    "costfun <- function(label, pred) {\n",
    "  pred  <- as.array(pred )\n",
    "  label <- as.array(label)\n",
    "  return(sum((label - pred) ^ 2) / length(label))\n",
    "}\n",
    "\n",
    "# create custom mxnet metric for training\n",
    "mx.metric.CRPS <- mx.metric.custom(\"CRPS\", costfun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you are ready to train your architected network with training data to build a model.  Below we offer two options - use our pretrained model provided in the first code block below that was created with 65 training rounds (epochs) or train your own model in the second code block below using a smaller number of training rounds (five training rounds will likely take 10+ minutes to execute).  Once you feel comfortable, feel free to play around with the learning rate, momentum etc. and see how the training performance changes.\n",
    "\n",
    "NOTE: No changes are needed in the first of two code blocks that follow this explanation and you only need to execute one of the two code blocks.\n",
    "\n",
    "Notice that one of the arguments provided is ```ctx = mx.gpu(0)``` which effectively tells MXNet to use GPU with device id equal to 0.  If you want to perform the training with the CPU instead then set the context to ```ctx=mx.cpu()```.  Furthermore, an array of GPU devices can be passed into the context for multi-gpu training if there are more than one GPU devices present on the machine.  Keep in mind here that using multiple GPU devices for training incurs some overhead - especially with large datasets.  It's not unusual that training time actually increases when using multiple GPU due to the additional data transfer required. So definitely start small and grow the training strategy and always profile execution times etc. to keep an eye on learning performance. \n",
    "\n",
    "If you want to save some time, simply load the pretrained model using the cell below and skip the cell that follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# load the pretrained result\n",
    "systole_model <- mx.model.load(\"stytole_model_BN\",65)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute the cell above if you would like to use the pretrained results OR execute the code below if you would like to train the model yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Training the systole model (this takes a long time).\n",
    "\n",
    "# set the random seed so the results are reproducible\n",
    "mx.set.seed(0)\n",
    "\n",
    "# start the clock\n",
    "start_time <- proc.time()\n",
    "\n",
    "# train the model parameters with the data\n",
    "systole_model <- mx.model.FeedForward.create(\n",
    "  X             = data_train,\n",
    "  ctx           = mx.gpu(0) ,\n",
    "  symbol        = network   ,\n",
    "  num.round     = 5         ,\n",
    "  learning.rate = 0.01      ,\n",
    "  wd            = 0.0001    ,\n",
    "  momentum      = 0.9       ,\n",
    "  eval.metric   = mx.metric.CRPS\n",
    ")\n",
    "\n",
    "# get the time delta\n",
    "proc.time() - start_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you've trained a network and want to save the results then export the trained model so it can be reloaded without having to perform training phase next time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Save the trained neural network model to disk.\n",
    "# This creates two files: \n",
    "#    1. systole_model-0065.params\n",
    "#    2. systol_model-symbol.json\n",
    "#\n",
    "# WARNING: the 65 here must match num.rounds parameter provided to mx.model.FeedForward.create\n",
    "#          So if the model was saved at 65, then it must be loaded with 65.\n",
    "#\n",
    "mx.model.save(systole_model,\"systole_model_BN\",5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to get a better feel for the data, lets extract all the data from the CSV iterator.  Doing this reads in the entire dataset in to host memory.  This can be prohibative with a large dataset.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# this reads in all the data at once (not a great idea, but we can do it)\n",
    "labels <- mx.io.extract(data_train,\"label\")\n",
    "data   <- mx.io.extract(data_train,\"data\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the data has been loaded we can check the dimensions.  As advertized we have ~ 5000 training samples each of size 64x64x30."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>600</li>\n",
       "\t<li>5293</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 600\n",
       "\\item 5293\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 600\n",
       "2. 5293\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]  600 5293"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>64</li>\n",
       "\t<li>64</li>\n",
       "\t<li>30</li>\n",
       "\t<li>5293</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 64\n",
       "\\item 64\n",
       "\\item 30\n",
       "\\item 5293\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 64\n",
       "2. 64\n",
       "3. 30\n",
       "4. 5293\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]   64   64   30 5293"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# \n",
    "dim(labels)\n",
    "dim(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the trained model we can pass input samples into the model to predict their systolic volume.  To select various samples we index into the data variable ```data[,,,1:20]``` to extract some samples. Notice here we're selecting input data along the 4th dimension.  That is we're selecting ```N``` input samples with dimension 64x64x30 and by leaving the first three index empty we are effectively saying \"give me everything\" along that dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "   0.00    9.00   24.00   37.86   51.00  255.00 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "summary(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create index for selecting samples\n",
    "sample_index = c(1:20);\n",
    "\n",
    "# apply the model to just a few data to generate some predictions\n",
    "model_predictions <- predict(systole_model, data[,,,sample_index])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets have a look at what was predicted by the model for these input.  Try selecting a different index for the plot to view other predictions and their label.  Don't hesitate to predict more input in the cell above to investigate even further."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "0.0143224532680549"
      ],
      "text/latex": [
       "0.0143224532680549"
      ],
      "text/markdown": [
       "0.0143224532680549"
      ],
      "text/plain": [
       "[1] 0.01432245"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAIAAAByhViMAAAACXBIWXMAABJ0AAASdAHeZh94\nAAAgAElEQVR4nOzde3xThf3/8XeS3qEtLZQ7yKByn7pRBWWA3Lyg3OQyZaIoXvZVBN104k9U\nHKig23TqUAGVqeAAHVSUqaCoKFQt3gFR6bgjBVpKr9Am+f2RLqQhvaVJTk76ej548Eg+5/bJ\n6Wny7rnF4nQ6BQAAAPOzGt0AAAAAAoNgBwAAECEIdgAAABGCYAcAABAhCHYAAAARgmAHAAAQ\nIQh2AAAAEYJgBwAAECGijG4gpJxOZ1lZWbAXYbFYgroIICDYVmEWbKswBdc3PoRgW7VYLHFx\ncdUNbVzB7vjx4//4xz9atGgRvEXY7Xar1cp7UGCF7LelUamoqIiKalzvAMHGhhoMTqfT4XDY\nbDajG4kobKvB4HA4JFmtwT0WmpubGxMT86c//am6ERrd23qrVq2mTp0avPkXFhbGxcVFR0cH\nbxGNjdPprKiosFqtvLMHVl5eXmpqqtFdRBS73e5wOKKiovi8DCC73V5YWNisWTOjG4koFRUV\nTqeTj6rAKi4uttlsNexLC4innnqqtLS0hhE4xw4AACBCEOwAAAAiBMEOAAAgQhDsAAAAIgTB\nDgAAIEIQ7AAAACIEwQ4AACBCEOwAAAAiBMEOAAAgQhDsAAAAIgTBDgAAIEIQ7AAAACIEwQ4A\nACBCEOwAAAAiBMEOAAAgQhDsAAAAIgTBDgAAIEIQ7AAAACIEwQ4AACBCEOwAAAAiBMEOAAAg\nQhDsAAAAIgTBDgAAIEIQ7AAAACIEwQ4AACBCEOwAAAAiBMEOAAAgQhgQ7GbNmjVixIg2bdpY\nLJZZs2bVPPKePXvGjRuXlJSUnJw8fvz4ffv21X0oAABAo2JAsHviiSeOHTt2xRVX1DpmUVHR\n4MGDd+zY8dJLL/3zn//8/vvvhwwZUlJSUpehAAAAjU1U6Bd5/Phxq9UqacGCBTWPuWjRol27\ndu3YsSM9PV1S7969u3Xrtnjx4unTp9c6FKgLp1PHjys5WcXFSkiQxXJqkN2uAwfUpo2iqv6W\n2O367DOVlCgtTb16qaJCeXmy29W2rawefyiVl6uwUCdPqnXrysrhw/r+e1VUqGtXpaQoIcG7\nmZMnVVHho+5WWKj//ldpaUpMlNWqmJhTvf38s3bu1PHjattWJ07oyy9VVqa8PEVFKTZW6emy\nWlVaKkkVFSovV26ucnPjnU6lpys9Xbt2acsWHTmioiI1b66OHVVaquPHKyt2u06ckN0um01l\nZSotVWmp7HZZLLJaZbfL6ZTTKYtFUVGVjx0OWSyVRUk2W+XKcVUsFlVUVFYcjsrx3ZNYrYqK\nktUqp7Py1bkGORxVFufqxzW+a6hrWa6pHI7Kx66nVmvlDN0ju4ru+bsac83TVXGN6Vqc1Vo5\nieslu+d/2hJtki0q6lTdNR/XZmOzVTbs7s3ds8Wi6Ogqs7VaVV5euWjXq3aN6VpRrgldXXmu\nWHfn7lXnepl2u6xW2WynXo57BPdQ9w/R1YB7/q7V5X7gGuSelavotZJd/bifej52v0DXA9e2\n4V6HNpvKy0+tH9ePwGq1WSxJribdv6HuB+75uBbkXquukV0zca8910xcq8j1s7DbK+fmei2u\nOXj+yrsmcW9snhuV+yflfpmuumuzcW8DrjXm/ul4rn+b7dSm6JrcvYV4Tu7eOD23QNc8XS/B\nvVzP9ePeql1F90L/90sUdfprcfNcuvun417Pnr/aXpM4HJUvwWsSefy6uSuei3a/HPdmdvqs\nPJfr+ZI9N073y3RvY+6nXnzWPV+d+4fi+Ro95+9+43K1arEoPj6ha1fnNdfommuqvN4QMyDY\nWev8ctesWdOvXz9XbpOUnp7er1+/zMxMV3SreShQUqKlS7VunYqL1auXbrtNb7+tL75QWZks\nFu3fr6++0pEjpz4YYmKUmKiSEp04caroYrUqNlZJSbJadejQqaHuz3svrt9592P3R7vXOO4x\n3R+T+t97vfutMDq68r3bXaxuQfUX7/eUQGhxOjhMwbJ/v2XDBi1frrfeMizbGRDs6m7r1q1j\nxozxrPTq1WvNmjV1GYoIsHevXnhBP/xg6djRdtFFzrPPVmpqnSZ88kk9+aRyck6FnrVr9dhj\ntUx18qSOHvU9yOGo3E11et0nr79WvSKd5ziu/ysqTtU9R3Y4dOJETT03INUBAILi7be1cKF+\n/3tjlh7Wfwbl5+enpKR4VlJTU/Py8uoy1OXIkSMZHiZOnOhwOE4Ek91uP3nyZFAX4ZPjvPPs\nkyZ5F0eOdAwb5lW0T53q7NPnRGFhleL/+3/OPn1O5uR4FiueeMLZp8/JzZs9i+Wvvurs06c8\nM7NK8f33nX36VCxa5Fk8+d13zj59KubOrbL43Fxnnz726dO9W73gAsdvf+tZWbu2/ItfjBs4\ne8iyZZo3zzpkiK15c0VH68X4W35qljH6ouLWrZ0pKTrrLOfdd1e8dd4DOxIzJl2Qc+GFjthY\nzZihnTt1s/PZbGX01aeem8QYrc5Wxii94Vk8X5uzlXGjFnkWOysnWxmzNNez2ETF2cp4Srd5\nbasf6MKVmuBVXK7ffqSBXsUnNT1bGYkq9Czeq4eylZGunzyLN2hxtjL66xPP4kityVbGFfq3\nZ/FcfZ6tjP/TM57FjtqTrYzZmu1ZjNWJbGU8o//z6mq9hq3SWK/iUv1uky7wKv5Nf8hWRory\nPYt3a362Mrrre8/iFC3JVsYgfehZvFT/yVbGRK3wLP5KX2Yr4zY95VlsqwPZypirKtdXRaki\nWxmLdKNXV+/o4jUa6VV8Sdd8pvO8io/prmxlpOmwZ/FO/SVbGb31nWfxar2SrYyhes+zOFzr\nspUxScs8i2fpm2xl3KHHPYutdChbGfM006uBbGW8qOu8im/psrUa4VV8QddnK8OiKuH9Ed2T\nrYzW+tmzeLueyFbG2fras3il/pWtjIv0rmdxsDZkK2OyXvYs9tS2bGXcqb94FpvraLYyHtNd\nXl1lqd/LmuxVzNTod3WRV3GhbspWRrTKPYtzdF+2MtqryoVu0/R0tjL6aItncbxey1bGCK31\nLA7QxmxlXKcXPYtd9UO2MmZqnmcxWQXZynhcd3h19bF+86qu8iq+rnHva4hXcYFuyVZGvKr8\nPfeAHsxWRift8iz+Xo3lreY8fXb6W80Z2n36W02cyny+1bynoae/1SzTpNPfah7XHae/1czU\nvIa81fxaX5z+VtNO+7OVMUf3eRZreKt5Q6O8ii/pmk/V16u4alWwwoaztj/ow3qPXcPZ7fac\nnBz309LSUlcxeEt0OBwOhyOoi/ApZtcuZ1yc13Kj9+61FBd7FaMOHlROjr2iovLcFkmSLTdX\nOTn2sjKnx8jWo0eVk+MoLnZ4FC3Hjiknx3H8eJXZFhYqJ8eRl+dZtJaWKifHefiwZ9FSXq6c\nHHXp4tVVzK5dTqv1/ff17bcWu10//GBdtizqY/ve5qqyD62iQikVB5uX5axfZymRRdKxY5Zv\nv43qqKMXKCd7s/1Hj79VUpTfWTle78tJOt5ZOV5vdgkq6awcr3eQGJ3srByvBGCVo7NydqmT\nquqkXbHy3rfWXvs6aK9XsbV+7qwcq6rs62uhI52VE6OTnsVmOtZZOQmqcj1QUxV1Vk6SjnsW\n41XaWTmpqvJXTbTKOyunpXI9ixY5Oytnv9p5dXWGdhco2avYTvt/of96FVvpUGfl2FTlxxfi\n/g+p1en9lynO7/6b62hn5Xj9BF39N1Hx6f0nq8CzGKeyzsrx2lajVNFZOa10yKuBzsrJk/ee\n547a4xXgJLXVgc7K8Sq2VG5n5USpwrOYqrzOyolTmWcxWQWdldNURZ7FJirurJxmOuZZjNWJ\nzsppoSOeRZvsnZWzXT28GviF/lsi7/NAO2iv14qS1EYHOyvH63Wl6XBn5XilPVf/Xr+qAen/\nR515ev+O0/ZodJD3W029+uetpl5vNcfUzKtY31/VuvSfqMKG9/+zWquqM7S79LSTWHz2X1rq\nDH0SqOQ0jqR77723hhFatmx50003eVZuvPHG1q1b12WoT8eOHVu8eLG//dbJ8ePHT548GdRF\nRKQTJ5wjRjgl/vGPf/zjH/9M/+/uu4P1cfnkk0/Onz+/hhHC+lBsr169tm7d6lnZtm1bz549\n6zK0Ufj8c61cqcLC2scMew88oLVrax8NAIAwd8YZuucew5Ye1sHu8ssvz8rKch9L3blzZ1ZW\n1qhRo+oytFFYsEATJ2r/fqP78Mfx49q8WRs3atcuvfmmnnvO6IbQYKffUCB4c3BfU1zDHDxv\nM+F5Jwg/+qlucp+Lq7mr0wfVtzHPBXmuB58vtro1cPokNffjc4mn170m8fnYa/KaX2PNap68\nukX7XDmnP/Wsn77QGl716Sunhtflc4k+H1c3js+fXa0biU772Xk9qG5Wp//zmsTn4+q2vdNH\nqG7yGsbxuaDqWq1hUK2vtLqfpvuf1aqmTTVhgjZuVLL3uS2hY8A5dq+99pr78ffff+96On78\n+MqGoqLuvffeBx98UNKNN9741FNPjR49eu7cuZJmzZrVqVOnG264wTVmzUMbhQkT1L270tKM\n7qPeHn1Us2f7uMg09CwWRUfLbq9yLarr5iZWq+LiFBurEydUWKiKiirXwMbFyW5XebksFiUk\nKDZW5eWKiVFsrJxOuc5gtFp18qSKiysveo2NVbNmSkuT1arDh3X8eOWN61y3K1PVe1m5Zlhe\nXlm3/O+eSU2bqmVLWSzKy1NRkZxOxcerWTOlpCglpfI2cvHxatlSTZsqNVUWiw4e1NGjKilR\nUlJl27GxSk5WSkqJzZbwwQc6dkzR0UpP1+TJSk7WTz/p8GEVFal7d7VvryZNlJKi4mJZrZW3\ng0lObtB7VlGRmjb1f/JwZrfbHQ5HVFSU5fQPfPjLbrcXFhY2a+Z9YhYaoqKiwul0Rrtun4gA\nKS4uttlscXHe5/uGmAHBbsKEU1f0vP7666+//rok5/+u8rDb7e7zDRMTEzds2HDHHXdMnjxZ\n0rBhw5544okmTZrUZWijMGKERnhfSRfO8vL0zjuaM0fbt4d0ue47w8XEaPBgjRmjAwfUurWG\nDFH37pXjbNmib7+Vw6G+fdWrl+/55Ofr3XdVUqLevXXuuZJUWlqZAhvi4EF98YWSkpSRoeho\nlZWdyj3FxYqNPXW71MBGory8stTUhJneV22qT5+ALcKnSE11ABAODAh2zhqv1PUa2qlTp1Wr\nVlU3cs1DEQ7Ky/Xoo/rHP5Sb6/t2bkHSooWefVYXXCBJbdpI0rFjquFv/j59ag80KSn67W+r\nVOIDcYvfNm102WWnnnrmHq+/U4hEAICaRfjtTmAsh0OXX6533619zLqIilJ8/KlrRXr0UGqq\ntm2TxaIOHTRxoiRt367oaF11lYYO9b7rN0dyAAARj2CHIFq61P9Ul5am0aNlsejIEWdqqrNV\nK8u4cZb0dH3+uQoL1a2benjfYAsAgMaOYGdmu3bp6FH17t3Qk7yC4PhxPf64XnnFz8nbtdPO\nnZUvy+lURYXdarXabDZJQ4cGrksAACJLWN/uBLW4/35lZGjPHqP7qKK8XDNmqFkzzZ6tn36q\nfXyfHnwwDMMqAADhjj12ZnbZZWrbVlW/MNdw992nJ59s0BzOOMP7GgUAAFAX7LEzs9/+VvPm\nqUULo/s45fhx/e1v9Rg/LU1vvqmuXU9VOnbUihVc/gkAgD/YY4dA2rWr8p66ddGmjVauVP/+\nuuQSvfee9uxRq1YaOlQJ3t8wDgAA6oRgh0Cq17dgrFqlvn0lyWbTRRcFqSMAABoRDsUikNq0\n0fDhdRrzggt03nlB7gYAgEaGYGdmubnKyanHsc8gKyvTjBnasKH2Mfv21b/+5eP7swEAQEMQ\n7MzsjjvUpYv27jW6j0rXXKMnn6z8wvvqXHGFNm7Upk3q0CFUbQEA0Ghwjp2Z9e2r8nLv7xM1\nyAcfaOXKWsYZNkzLl5/6SnsAABBY7LEzs+nTtWKFWrUyug9J+uKLagdZLGrXTn//u9atI9UB\nABBEfMwiAMrLlZVV7dDNmyuvfgUAAEFFsENDVVRo6FBt3Oh76CWXcPUrAAAhwqFYNNTTT1eb\n6oYN08svc/UrAAAhwh47NNTSpT6KiYn6/HN16xbybgAAaMTYY2dmV14pi0W7dxvYwsaNys72\nUU9OJtUBABBqBDsz+8Uv1KePYmIMbGHuXN/17t1D2wcAACDYmdsjjyg7W23aGNjC9u2+6+vX\n6/33Q9sKAACNHsEO/tu+XQcOVDt02rQQtgIAAAh2aIhbbpHdXu3Q7dtVUhLCbgAAaPQIdqg3\nh0PPPaf+/fXhhzWNFh1t7Ol/AAA0OtzuBPV2221asKD20S67jC8QAwAgpNhjZ2bTpqlLF/38\ncyiX+cUXdUp1nTvrmWeC3w0AAPDAHhXUz6ef1jS0UycNGKBzz9UNNyg+PlQ9AQAASQQ7c3v6\n6dAvs7qjqwMHatIkTZ3K4VcAAAzDhzDqZ+hQxcWprKxKsWNHvfcekQ4AAINxjh3qp3Nn/fnP\nVSrR0XrhBVIdAADG49MY9XbXXTr7bL38snbuVG6unE49/LD27tW118piMbo5AAAaMYId/HHR\nRbJaNXx45dOcHL3/vj79lCthAQAwEodizeyuu5SRoSNHQrzYigqVlGjqVO/6s89q06YQ9wIA\nAE4h2JnZzp3askXl5SFbYE6Oxo5VYqKSkrRnj48RPvooZL0AAABvHIo1s3//O5RLy8vThRdq\n796axrHylwIAAMbhcxh1NWNGLalO0pAhIWkFAAD4QrBDnXz7rV59tZZxxoxRRkZIugEAAL4Q\n7FAns2fLbq9lnKuuCkkrAACgGgQ71Mn27bWMEBOjCy4ISSsAAKAaBDsze/RRTZyogoIQLCo1\ntZYRHnpI7duHoBEAAFAtgp2ZbdqklSu9v7c1OH73u5qGLl+uO+8MQRcAAKAmBDsze/ZZ7dyp\nFi1CsKjf/97HHYndLr00BC0AAIBaEOzMrHVrde4smy0Ei7JYtHix/v53H4P69lViYghaAAAA\ntSDYoR6mT9cNN1SpJCbq+ecN6gYAAFTFN0+gfhYu1ODBeuMNHT+url11111q187ongAAgCSC\nHerLYtGkSZo0yeg+AADAaTgUa2ZLlmjmTBUXG90HAAAICwQ7M3v9dc2fr9JSo/sAAABhgWBn\nZo8/ruxsNWtmdB8AACAscI6dmaWnG90BAAAII+yxAwAAiBAEOwAAgAhBsAMAAIgQBDszW7lS\n8+errMzoPgAAQFgg2JnZiy9q5kxudwIAAFwIdmZ2zz1asUJNmhjdBwAACAvc7sTMBgwwugMA\nABBG2GMHAAAQIQh2AAAAEYJDsagTh0Ovv65vv1VCgoYNU0aG0Q0BAIDTEOzM7KOPdOiQxoxR\ndHRQl1NUpGHD9OmnlU/vuUezZmnOnKAuEwAA1BuHYs3soYc0cWIIbnfypz+dSnUuc+dq3bpg\nLxYAANQPwc7Mpk7VvHmKjQ32claurGsRAAAYiEOxZjZxYmiWc/y4j2JBQWgWDgAA6oo9dqhd\nr14+ir/8Zcj7AAAANSLYoXaPPupdOeMM3XabEa0AAIDqEexQu2HDtHKlOneWpKgoDRumd95R\ncrLRbQEAgKo4x87MfvhBhYX61a9kDXpAHz9e48eruFgxMcG+uQoAAPATe+zM7JZblJERgtud\nuDVpQqoDACB8scfOzC67TF26KIofIgAAkAh25nbHHUZ3AAAAwgiHYgEAACIEwQ4AACBCEOwA\nAAAiBMHOzAoLlZ8vpzPYy/nuO02apB49NHiwFi+WwxHsBQIAAH9w8YSZjRqlDz5Qaani4oK3\nkI0bNXy4TpyQpO+/1wcf6JNP9OKLwVsgAADwE3vszGzQIE2YIJstqAuZOrUy1bktWaL33w/q\nMgEAgD/YY2dms2cHewkHD+rHH33UP/pIQ4YEe+EAAKB+2GOHmlT3XWXB/w4zAABQb3w+oyat\nWql3bx/1oUND3goAAKgNwQ61WLJECQlVKjNmqH9/g7oBAADVa1zn2DmdTqfTWV5eHrxFOByO\nioqK4M0/9M46S19/bfnrX61ffWVp1co5aZJz7FhHMFehbw6Hw8F9VgIq2L8LjVaEvQMYzvW7\nz7YaDKzVwLLb7SF4X3XWdo+zxhXsLBaL+/+gLiXYi3CxDR5s+fjjipISRUcHdUGdOumppzxD\nVShenSfXdhyatdqosEoDiw01eFirgcW2GgyW/zG2jcYV7CRZLJaoqCC+aqvVarPZgrqIU7p3\nV2lpVHS0QrM4gzidzoqKCteKNbqXiBLs34VGyG63OxwOm81m+Dt7JLFYLFarlW01sCoqKpxO\nJ2s1sEITAGp9e+GHamaLFhndAQAACCNcPAEAABAhCHYAAAARgmAHAAAQITjHDr599pnee08V\nFerbV8OHi3PBAQAIfwQ7M5swQV98oR9/DPg3fM2apYceOvV04EC9/bbi4wO7EAAAEGAcijWz\nwkLl5wd8rhs2VEl1kj76SPffH/DlAACAACPYmdnbbysvL+C76/71Lx/FRYv05psqLQ3sogAA\nQCAR7OAtL89HsaBAI0eqZ09t2RLyhgAAQN0Q7OCtR49qB+3apQkTVFwcwm4AAECdEezgbdo0\npaZWO/S//9X774ewGwAAUGcEO1Rht2vKFN9HY91yc0PVDQAAqA+CnZnddpsuuiiws3ziCf3n\nP7WM07VrYJcJAAACg/vYmdnnn+uzzwI7y8zMWkYYPFj9+wd2mQAAIDDYY2dmWVlyOAI7y5IS\nH0XX105YLBo3Tq+8EvD7qwAAgMDgIxpV/PrXvotPP63Dh/Xaa2rbNuQ9AQCAuiHYoYr771dK\nindxyxZNm6b+/fXjj0b0BAAA6oZghyrat9fGjRo7Vs2aeR9y3bFD11wjp9OgzgAAQG0IdvDW\nq5f+/W8tWeLj/L2sLG3bZkRPAACgDgh2ZnbPPbrqqiDNu7qb1R05EqQFAgCAhiLYmdn69Xrt\ntSDNu317H0WLRWeeGaQFAgCAhuI+dma2erVOngzGjB0O/e1vPuo33cRVsQAAhC+CnZm1axek\nGa9apfXrvYtNm/pOewAAIExwKBY+fPGFj2JRkY4dC3krAACgzthjhypKSvT001q1yscgi0UJ\nCSFvCAAA1BnBDqccPaq+fbVzp++hgwapWbPQNgQAAOqDQ7Fm9uSTmjUrgPP74x+rTXXt2unF\nFwO4KAAAEHgEOzNbskR/+UsA5/ef//gonnGG/vpXbd+uTp0CuCgAABB4HIo1s2eeUXFxAOfn\n894pffroD38I4EIAAECwsMfOzPr21ZAhAZxfv34+iuefH8AlAACAICLY4ZQnnlDTplUq55yj\nadMM6gYAANQTwQ6nREXpH//QxRerQwd166Y779QHHyguzui2AABA3XCOHSQpP1/XX6/Vqyuf\nnneeXnmFr4UFAMBk2GNnZq+9phdeCMicbrzxVKqT9NlnGjdOZWUBmTcAAAgRgp2ZzZmj225r\n+Gz27tXrr3sXv/3Wx9fFAgCAcMahWDN77DEVFTV8Nnv3+q7v3t3weQMAgNAh2JnZRRcFZDZt\n2/qun3FGQGYPAABChEOxUKdOGjPGu9ijh4YNM6IbAADgL4IdJGnhwiq7/84+WytWcKMTAABM\nhkOxkKS0NL3zjnbs0Pbt6thR55wjK5kfAACzIdiZ2TvvqKREY8cGan7duqlbt0DNDAAAhBq7\nZczsj3/Utdca3QQAAAgX7LEzszvu0PHjRjcBAADCBcHOzKZONboDAAAQRjgUCwAAECEIdgAA\nABGCYAcAABAhCHZmtnWrvvrK6CYAAEC44OIJMxs/Xj//rPx8o/sAAABhgWBnZhMncrsTAADg\nRrAzswcfNLoDAAAQRjjHDgAAIEIQ7AAAACIEwQ4AACBCEOzM7PhxHTtmdBMAACBcEOzMLCND\n6elGNwEAAMIFV8Wa2QUXqLDQ6CYAAEC4INiZ2ZIlwZhrUZHi42WzBWPeAAAgiDgUi1MWLlTb\ntkpMVEKCxo3Tvn1GNwQAAOqDPXaotHixbr658vHJk/r3v/XTT8rKUny8oW0BAIA6Y48dJMnh\n0P/7f97Fb77RsmVGdAMAAPxCsIMkLVumw4d91LdtC3krAADAXwQ7M/v1r9W1a8NnM3++Jk/2\nPSg1teGzBwAAIcI5dmbWoYNKSxs4j927dd99vgfFxWncuAbOHgAAhA7BzswyMxs+j6wslZf7\nqMfE6Ikn1L17w5cAAABChGDX2DmdvusvvqhJk0LbCgAAaBjOsWvsfvMbxcR4F5OTNXKkEd0A\nAIAGINg1duvWyeHwLj79tBITjegGAAA0AMGuUfvxR02bpoqKKsVf/UpXX21QQwAAoAEIdmY2\ncKDOP78hM1izRiUl3sUvv/R9TzsAABDmCHZmVlCgY8caMoPCQt/148cbMlcAAGAMroo1s6+/\nbuAMfvlLH8VmzdSxYwNnDAAADMAeu0Zt9GhdeKF38bHHFB1tQDMAAKCBCHaNms2mf/9bt9yi\n5GRJOvNMvfSSbrjB6LYAAIBfOBTb2KWk6B//0NNPq6RETZoY3Q0AAGgA9thBkiwWUh0AAKZH\nsDOzSZN01VVGNwEAAMIFh2LNbONGozsAAABhhGBnZt99Z3QHAAAgjBDszMx1LSsAAIAkQ86x\n27Nnz7hx45KSkpKTk8ePH79v377qxrT4EhcXV8MI37ETCwAANFah3mNXVFQ0ePDg+Pj4l156\nSdKsWbOGDBny1VdfJSQknD7y5s2bPZ/++OOP11xzzbhx4zyLU6ZMufnmm91Pu3TpEpzGAQAA\nwl2og92iRYt27dq1Y8eO9PR0Sb179+7WrdvixYunT59++sj9+vXzfLpq1SpJN910k2exXbt2\nXqMBAAA0TqE+FLtmzZp+/fq5Up2k9PT0fv36ZWZm1jpheXn5kiVLunbtOmjQIK9BTqcz8I2a\nwt136557jG4CAACEi1AHu61bt/bu3duz0qtXr23bttU64RtvvJGbm3vjjTd61T2xXBcAACAA\nSURBVBcsWBAbG9ukSZOhQ4d+8sknXkMdDke+h+PHjzew//DyyitatszoJgAAQLgI9aHY/Pz8\nlJQUz0pqampeXl6tEy5atCgmJubaa6/1LE6dOnXEiBFt27bds2fPY489Nnjw4Pfee2/AgAHu\nEQ4dOtS2bVv30/bt2993331BjXcnT5602+1WaygSs/Xdd2WxOCIsrfricDhcF8cY3UhEqaio\niLQ/dYzmdDqdTmdofv0bD6fTefLkSbbVwHI4HJLYVgOrvLzcYrGcPHkyqEtx/exqYI7bneze\nvXvdunXjx49PS0vzrC9evNj1oF+/fpdddlnv3r3vv//+DRs2uEeIi4ubMGGC+2lSUpLVam0S\n5C/Pio2NjYoKyYrt3j0USzGa0+m02+0Wi8VmsxndS0QpLy8P9u9CY2O3251Op81m44+QALLb\n7Q6Hg201sFzbaog+qhqNkpISm80WGxsb1KXUGsdD/UNNSUnJz8/3rOTl5aWmptY81fPPP+9w\nOE4/DuupSZMmI0eOfP75570Wt2LFCvfTgoKC1157Laj5wGKxWK1WIkgAufeCsFYDi6wcDA6H\ng2AXcGyrAed6X2WtBpbVag2Hj6pQ74bt1avX1q1bPSvbtm3r2bNnDZPY7fYXXnihc+fOQ4cO\nDXJ3AAAAJhbqYHf55ZdnZWXl5OS4nu7cuTMrK2vUqFE1TPKf//xn//79N9xwQ81/BBcVFb3x\nxht9+/YNZLsAAADmEepgd+ONN3bo0GH06NGZmZmZmZljxozp1KnTDTfc4B4hKirqgQce8Jxk\n0aJFUVFR1113ndesHn/88dtvvz0zM/Pjjz9+9dVXBw4cePDgwblz54biZYSJuXP16KNGNwEA\nAMJFqINdYmLihg0b0tPTJ0+ePHny5DPPPPP999/3PCvWbrfb7Xb304MHD7711lsjR45s3bq1\n16y6dOny6aefTpkyZdCgQbfddluHDh0++eSTCy64IESvJBw8/bSefdboJgAAQLgw4IqYTp06\nub5Dwievuw23adOmoqLC55ijRo2q+Rhu5FuxQkafpAkAAMIHlzqb2cCBRncAAADCCDcnBAAA\niBAEOwAAgAhBsGvU9u3TlVcqJUWxserTR++8Y3RDAACgAQh2ZvbSS/rXv/yeurRUF12k5ct1\n7JhOntQXX+iSS7RxYwD7AwAAIUWwM7M//lH33ef31AsXavt27+IddzSoIwAAYCCuijWzv/5V\nMTF+T/3FFz6KX30lh0O1fcUwAAAIRwQ7M7vmmoZMnZDgu0iqAwDApPgMb7x83t159OiQ9wEA\nAAKEYNd4XXqpbr+9SqV7d/397wZ1AwAAGoxDsY3a449r5Ei9845KStS7t6ZMUWys0T0BAAB/\nEezM7MMPFR2tCy5oyDyGDNGQIYFqCAAAGIlgZ2Zjx6plS33/vdF9AACAsECwM7MHHvB9aSsA\nAGiUCHZmNmOG0R0AAIAwwlWxAAAAEYJgBwAAECEIdgAAABGCYGdmX32lrVuNbgIAAIQLLp4w\ns0GD1L492Q4AALgQ7MxsyhSlpBjdBAAACBcEOzPji10BAIAHzrEDAACIEAQ7AACACEGwAwAA\niBAEOzM7dkzHjxvdBAAACBcEOzNr316/+Y3RTQAAgHDBVbFmNniw2rQxugkAABAuCHZmtmaN\n0R0AAIAwwqFYAACACEGwAwAAiBAEOwAAgAhBsAMAAIgQBDsza9tWgwYZ3QQAAAgXBDsza9ZM\niYlGNwEAAMIFtzsxs23bjO4AAACEEfbYAQAARAiCHQAAQIQg2AEAAEQIgh0AAECEINiZWb9+\nmjzZ6CYAAEC44KpYM/vvfxUfb3QTAAAgXBDszOzQIaM7AAAAYYRDsQAAABGCYAcAABAhCHYA\nAAARgmAHAAAQIQh2ZnbxxZo2zegmAABAuOCqWDNbv17HjhndBAAACBcEOzM7ckRR/AQBAEAl\nYoGZpaQY3QEAAAgjnGMHAAAQIQh2AAAAEYJgBwAAECEIdqbldOrmmzV/vtF9AACAcEGwM7OF\nC5WZaXQTAAAgXHBVrGlZLMrOVpMmRvcBAADCBcHOzPr0MboDAAAQRjgUCwAAECEIdgAAABGC\nYAcAABAhCHam5XBo/nwtXWp0HwAAIFwQ7EzL4dDMmXr2WaP7AAAA4YKrYk3LZtO6dWrWzOg+\nAABAuCDYmZbFomHDjG4CAACEEQ7FNkarV+ussxQTo7g4paaqe3ddf7327TO6LQAA0DDssWt0\n1q7V2LGnnp44ofx87dihtWv1zTdq2dK4zgAAQMOwx67RmTHDd/3QId17b2hbAQAAAUWwMy27\nXQsXas2aek1UVKSffqp26MaNDW0KAAAYiGBnWhUVuvlm/fWv9ZooJkZR1R9+t9ka2hQAADAQ\n59iZVlSUnntObdrUa6KYGI0Zo9de8z2Uq2wBADA1gp1p2Wy66SY/pnv6aWVna9cu73q3bpoz\np+FtAQAAwxDsGp1WrbR1q154QVu3ascOVVSoWTOdd57uuENNmhjdHAAAaACCXWOUkKBp04xu\nAgAABBoXTwAAAEQIgp1p2e1av15bthjdBwAACBcEO9M6cULDh+tPfzK6DwAAEC44x860oqN1\n993q0sXoPgAAQLgg2JlWdLTmzTO6CQAAEEYaY7BzOp3Bnn+wFxEQH3ygTz5RbKwuuUS9exvd\nTfXcK9MUa9VcWKWB5VqfrNXAYq0GD2s1sJz/Y2wbjSvYudZ4RUVF8BbhcDjsdrvFYgneIhqu\nokKTJtlWr648w/JPf9J99znuu89ubFc1czgcDofD6C4iSrB/Fxotuz2sf5VMx/W7z7YaDKzV\nwHJ9SAV7rdYaHBtXsHPlLas1iJeMWCwWq9Ua1EU03KOPWtypzmXOHOt55+nSS8P0rzeHw2Gx\nWMI8LptRmG+opuP605G1GlhOp9P1vmp0IxHFFUFYq4EVmgBQ60dh4wp2kiwWiy2Y33Xv+rkG\ndRGVHA7t2qX4+Pp+XaykV1/1UVy2zHr55QHoK+CcTqcr2IVirTYmrNKAs9vtrmDHHyGBxbYa\ncK4/QlirgeVKdYavVdK6aRUXq0sXXXutH5Pm5dW1CAAATKTR7bGLHNHRuukm9ejhx6Q9eujw\nYe9iz54BaAoAABiIYGdacXF67jn/Jn3oIQ0YUKWSlqY77wxAUwAAwEAcim2MfvMbvfmmunaV\nJItFAwZo/Xq1bWt0WwAAoGHYY9dIXXaZLrtMhw8rJkbJyUZ3AwAAAoFg16ilpRndAQAACBwO\nxZqW06n8fBUVGd0HAAAIFwQ70yooUGqqJkwwug8AABAuCHamFRWlYcP0q18Z3QcAAAgXnGNn\nWk2bat06o5sAAABhhD12AAAAEYJgBwAAECEIdgAAABGCYAcAABAhCHamlZen1FRdeaXRfQAA\ngHBBsDMtq1UpKWra1Og+AABAuOB2J6bVrJl27jS6CQAAEEYIdo1OcbFefFEHD+qMMzR6tFq1\nMrohAAAQIAS7RuTkSc2cqSeflN1eWfnjH/XSSxo71tC2AABAgHCOXSNy++16/PFTqU5SUZGm\nTlVurnE9AQCAwCHYNRY//KBnnvFRz8/XW2+FvBsAABAEBDvTys9XRob+8Ic6jv7ttzXNCQAA\nRADOsTOt8nJt2aJ27eo4emJitYN69QpMRwAAwFgEO9Nq2VJOZ91H799fbdvqwAHv+uDBGj48\nkH0BAACjcCi2sWjSRMuWKSXlVMVi0ciReu01WdkKAACICOyxa0QGDdKOHXr5ZeXkqH17TZ2q\ntDSjewIAAIFDsGtc0tLqfrkFAAAwGQ7CAQAARAiCnWnl5Wn4cM2ebXQfAAAgXHAo1rROnND6\n9TXdxQQAADQyBDvTatVKeXmKjja6DwAAEC4IdqZltVa5eQkAAGj0OMcOAAAgQhDsAAAAIgTB\nDgAAIEIQ7Ezr2DHdfLOee87oPgAAQLjwM9i9/fbbF198cfPmzW02m+U0gW0RvhUXa+FCrVtn\ndB8AACBc+HNV7PLly6+88kpJMTExrVq1svId8oZo2VLZ2VwYCwAA3PwJdnPnzo2Li3vhhRcm\nTJgQFcUNUwwSHa0+fYxuAgAAhBF/YtkPP/xw7bXXXnXVVQHvBgAAAH7z5yhq8+bNU1NTA94K\nAAAAGsKfYDdx4sT169c7nc6AdwMAAAC/+RPsXOfYXXvttfv27Qt4Q6ir48c1f74yM43uAwAA\nhAt/zrHr3r17RUXFJ5988vLLL7do0SI2NtZrBAJfKBw7ppkzNXGiRo82uhUAABAW/Al2+/fv\ndz8+cuRI4JpBfbRooRUr1KGD0X0AAIBw4U+w4+y6sJCQoAkTjG4CAACEEe4tDAAAECEIdgAA\nABHC/2D36quvDh8+vEWLFtHR0WlpaRdddNHy5csD2BkAAADqxZ9g53A4rrzyykmTJq1fv764\nuLh58+ZFRUXr1q278sorf/e733EGXogUF2vlSn36qdF9AACAcOFPsHv22WeXL1+ekZHx4Ycf\nFhUV/fzzz0VFRR999FFGRsayZcuee+65gHcJH3JzNXGi/v53o/sAAADhwp9g98ILL3Ts2HHD\nhg0DBw602WySbDbbgAEDNmzY0LFjx+effz7QTcKXFi303HOaOtXoPgAAQLjwJ9ht27btiiuu\naNq0qVe9adOmY8eO3bZtWyAaQ20SE3XTTRo61Og+AABAuOCqWAAAgAjhT7Dr0aPH6tWrS0pK\nvOrFxcWrV6/u0aNHIBoDAABA/fgT7K677rpdu3YNGzZs06ZNDodDksPh+OSTT4YMGbJ79+7r\nr78+0E0CAACgdv58pdj//d//ffDBB6+//nr//v3j4+OTk5MLCgpKS0slTZw48fe//32gm4Qv\nJSXatEmtW6t3b6NbAQAAYcGfPXY2m23lypUvv/zy0KFD4+PjDx8+HB8fP3To0KVLly5fvtxq\n5by9kDhwQMOHa/58o/sAAADhwp89dpIsFsvVV1999dVXB7Yb1ENKiu6+W7/6ldF9AACAcOFn\nsIPxmjfXvHlGNwEAAMIIh00BAAAiRF332LVv317Shx9+2KVLF9fjGuzbt6+hfQEAAKCe6hrs\n9u/fL6m8vNz9GAAAAGGlrsHO6XT6fAzDlJdr714lJiotzehWAABAWOAcO9PKyVGXLrrrLqP7\nAAAA4YJgZ1qJiZowQeeeW5dxS0p0//3KyFCPHrriCn3zTbCbAwAABvAn2FksltmzZ/scNHfu\nXIvF0qCOUEdt22rFCt16a81jFRbqgQfUvr3mzNGWLfr+e61apfPO01dfhaZLAAAQOtzHLpLl\n5ysjQzk53vUTJ3TLLdq0yYieAABA0AT4UGxZWVlMTExg5wm/3XOPj1TnkpWl8vLQdgMAAIIs\nYHvsnE7nvn371q5d26FDh0DNEw303nvVDrJaxTFzAAAiTD322Fn+R9KDDz5oqcpqtXbs2PHL\nL7+84oorgtYt6sfhqHbQhRcqiuPwAABElnp8trdr1871YP/+/YmJiUlJSZ5DbTZbixYtRowY\nMWvWrEA2iOps26ZevXT99Xr++epGGTDA96HYlBQ9+2wQWwMAAIaoR7Bzf1GYxWL5wx/+UN2F\nsQiRJk00bJh69qxhlPnz9e67OnjwVCUxUddfr5kz1bp10BsEAAAh5s/ROL55IiyccYbWrat5\nlFat9PXXmjNHGzcqOlqXXKK771aTJqHpDwAAhBqnWUW4tDQ9+aTRTQAAgJDw53Ynn3322ezZ\ns3Nzc73qhw4dmj17dnZ2diAaQ4BVVOjJJ9Wnj1q10vDh+uADoxsCAACB5k+we+SRR5YtW9ay\nZUuvesuWLV955ZVHHnkkEI0hwKZO1YwZ+uIL5eZq/XoNHqx//9vongAAQED5E+w2b948ZMiQ\n0+sWi2XIkCGbN29ucFcIsI0b9dJL3sXf/14VFUZ0AwAAgsOfYHfkyJEWLVr4HJSWlnb48OGG\ntYS62bZNqam6/fa6jJuV5aN4+LB++inATQEAAAP5E+ySk5P379/vc9D+/fubNm3asJZQNzab\nUlKUkFCXcav7mrfY2EB2BAAAjOVPsDv33HNXrVp16NAhr/qhQ4dWr16dkZFR8+R79uwZN25c\nUlJScnLy+PHj3bfH88lymu+++86/WUWabt20c6cefrgu4w4f7qPYo4c6dQpsTwAAwEj+BLtb\nb721oKBgyJAh7733nt1ul2S32997770hQ4YUFBRMmzathmmLiooGDx68Y8eOl1566Z///Of3\n338/ZMiQkpKSGiaZMmXKZg9dunTxe1aNVs+emj+/SiUxUS+/zNfFAgAQUfy5j93IkSPvuuuu\nxx57bNiwYXFxccnJyQUFBWVlZZLuuuuu0aNH1zDtokWLdu3atWPHjvT0dEm9e/fu1q3b4sWL\np0+fXt0k7dq169evX0Bm1Zj96U86/3y9/LL27dNZZ2n6dLVta3RPAAAgoPzZYyfp0Ucffeut\nt0aMGNG0adOjR482bdp0xIgRa9euffTRR2uecM2aNf369XNFMUnp6en9+vXLzMyseSqf33Xh\n36waswEDtHCh1q7VvHmkOgAAIpD/3zwxYsSIESNG1HeqrVu3jhkzxrPSq1evNWvW1DDJggUL\nHn300ejo6H79+v35z3/u379/3WdVUlLyyiuvuJ9arVYLRx8BAECECvVXiuXn56ekpHhWUlNT\n8/Lyqht/6tSpI0aMaNu27Z49ex577LHBgwe/9957AwYMqOOsCgoKbr75ZvfT9u3bz5o16+jR\no4F5Mb7Y7fYTJ06EID7adu5setNNJ0eOLK3bHU9MzbXLllAeWOXl5UH9XWiE2FCDwel0VlRU\nsK0GFttqMNjtdovFUlxcHOyl1DxCuH9X7OLFi10P+vXrd9lll/Xu3fv+++/fsGFDHSdPTk5+\n7rnn3E9de+yaN28e+Eb/p7CwMC4uLjo6OniLqLR3r77+Our88xOC+XLCgett3Wq12mw2o3uJ\nKHl5eampqUZ3EVHsdrvD4YiKiuLzMoDsdnthYWGzZs2MbiSiVFRUOJ3OUHxUNSbFxcU2my0u\nLi6oS6n1o7Cuwa59+/aSPvzwwy5durge16CG246kpKTk5+d7Vur+6dKkSZORI0c+//zzdZ9V\nQkLCTTfd5H5aUFDw2muv1WVZJnDOOfJ16iEAAGi06hrsXHckLi8vdz/2T69evbZu3epZ2bZt\nW8+ePY2dFQAAQASo61WxTqfT6XR2797d/bgGNczn8ssvz8rKysnJcT3duXNnVlbWqFGj6tJD\nUVHRG2+80bdv34bPCgAAIPL4ebsTv914440dOnQYPXp0ZmZmZmbmmDFjOnXqdMMNN7hHiIqK\neuCBB1yPH3/88dtvvz0zM/Pjjz9+9dVXBw4cePDgwblz59ZxVgAAAI1KqINdYmLihg0b0tPT\nJ0+ePHny5DPPPPP9999v0qSJewS73e6+4qNLly6ffvrplClTBg0adNttt3Xo0OGTTz654IIL\n6jgrAACARsWAq2I7deq0atWq6oZ6HskdNWpUzYdWa55VhMvJ0cyZuvhiTZ1qdCsAACAs1O+q\n2Dqq4apYBEx+vlauVMuWRvcBAADCRf2uikUYOecc5eUpNtboPgAAQLio31WxLmVlZaNGjera\ntevy5csPHTpUXl5+6NCh5cuXd+3addSoUWVlZUHtGJVsNqWkKCHB6D4AAEC48OfiiTlz5nz+\n+eebN2+eOHFiy5Yto6KiWrZsOXHixM2bN3/++edz5swJeJcAAAColT/BbunSpePHjz/96yJS\nU1PHjx+/dOnSQDQGAACA+vEn2B04cKC6L5iLjo4+cOBAw1oCAACAP/wJdu3atcvMzDz9XLqy\nsrLVq1fX6/pZ+G/XLt18s5YvN7oPAAAQLvwJdlOmTNm5c+fw4cOzsrIcDockh8ORlZU1bNiw\nnJycKVOmBLhH+JSbq4ULtWmT0X0AAIBw4c8NimfOnLlly5Y33njj/PPPj4+PT05OLigoKC0t\nlTR69Oi777470E3Cl169lJ3NfewAAICbP3vsYmJiVq9evXTp0mHDhsXHxx8+fDg+Pn7YsGFL\nly5dtWpVTExMwLuED02aqE8fdehgdB8AACBc+PmVYhaLZdKkSZMmTQpsNwAAAPCbP3vsAAAA\nEIb8DHbl5eWPP/74ueeem5iYaLFYXMVvvvlm2rRp27dvD1x7AAAAqCt/DsWWlpZefPHFGzdu\njIqKivX4rtIuXbq8+OKLcXFxf/nLXwLXIapx4IBeflkZGRo61OhWAABAWPBnj928efM2btx4\nxx13HD169M4773TXmzRpMnDgwA8//DBw7aF6u3dr5ky99ZbRfQAAgHDhT7B79dVX+/fv/7e/\n/S0pKclrUI8ePX744YdANIbadO2qFSt07bVG9wEAAMKFP4did+/ePX78eJ+DEhISSkpKGtYS\n6qZ5c02YYHQTAAAgjPizxy4hIaGwsNDnoL1797Zo0aJhLQEAAMAf/gS7jIyMNWvWnP5dsceO\nHVuzZs35558fiMYAAABQP/4Eu+nTp+/evXvcuHH//e9/3cXt27ePGjUqPz9/xowZgWsPAAAA\ndeXPOXYjR4685557Hnnkkc6dO8fFxUlKTU3Nz8+XNHfu3EGDBgW4R/h05Ig2bFDXrjr7bKNb\nAQAAYcHPGxQ//PDD77777ujRo5OSkqKioqxW64gRI95999177703sP2hWt9/r4kT9corRvcB\nAADChZ/fFStp+PDhw4cPD2ArqJ9OnTRvns47z+g+AABAuPAn2M2ePTs9Pf3qq68OeDeoh/bt\ndffdRjcBAADCiD+HYh966KFvv/024K0AAACgIfwJdm3bti0vLw94KwAAAGgIf4LdmDFj3nrr\nrYqKioB3AwAAAL/5E+z+/Oc/JyUljRkzZuvWrQFvCHVVVKQtW7R/f13G3blTixfrmWf0zTfB\nbgsAABjGn4snevXqVV5enp2d/dZbbyUmJiYlJXmNsG/fvkD0hhpt2aILL9Tdd2vevJpHnDVL\nf/mLTpyofHrNNXrxRVn9vNENAAAIX/4Eu/0ee4kKCwur+95YBFfHjrr7btV2O+i1a/XQQ1Uq\nL72kjAzddlsQWwMAAIbwZ7+NszYB7xI+/OIXmjdPl15a81hLlvgoPv98MBoCAAAG44BchDt8\n2EfxyJGQ9wEAAIKPYBfhzjzTR7Fbt5D3AQAAgs//YPfRRx9dd91155xzzi9+8Ytzzjnnuuuu\n++ijjwLYGQLizjvVtKl38b77jGgFAAAEmZ/n2N16662DBg1asmTJ119/vWvXrq+//nrJkiWD\nBg269dZbOccuREpLlZOjvLyax+raVW++qe7dK5+2a6eVK3XhhcFuDgAAGMCfYPfUU08tWLCg\ne/fuK1euzM3NLS8vz83NXblyZffu3RcsWPDUU08FvEv4sHmzunTRX/9a64iDBmn7duXmav9+\n7dun8eND0BwAADCAP8HumWeeadu27aZNm8aPH5+WlhYVFZWWljZ+/PhNmza1bdv22WefDXiX\n8KFlS02YoN696zh6Wpratg1qQwAAwGD+BLucnJxx48alpKR41VNSUsaNG5eTkxOIxlCb3r21\nYoWuusroPgAAQLjwJ9i1aNEiOjra56Do6OjmzZs3rCUAAAD4w59gN3bs2MzMzLKyMq96WVlZ\nZmbm2LFjA9EYAAAA6sefYPfQQw81b958+PDhWVlZDodDksPhyMrKGj58eIsWLR5++OFANwkA\nAIDa+fNdsb169XJdCXv++efHx8cnJycXFBSUlpZKatmyZc+ePT1H3rdvX2A6BQAAQI38CXb7\n9+93Py4tLXVFOpfc3NwANIW6WLdOF12kWbM0Z47RrQAAgLDgT7DjFsRhISlJffqoXTuj+wAA\nAOHCn2CHsNC3r7KzjW4CAACEEf+/KxYAAABhhWAHAAAQIQh2AAAAEYJgBwAAECEIdqb10Ufq\n0kV//7vRfQAAgHBBsAMAAIgQ3O7EtAYO1M6dRjcBAADCCHvsAAAAIgTBDgAAIEIQ7AAAACIE\nwQ4AACBCEOxM6+OPlZGhxYuN7gMAAIQLgp1pFRRoyxYdPGh0HwAAIFxwuxPTuuwyOZ1GNwEA\nAMIIe+wAAAAiBMEOAAAgQhDsAAAAIgTBDgAAIEIQ7EwrO1sTJ2rVKqP7AAAA4YJgZ1oHDmjl\nSm3fbnQfAAAgXHC7E9MaNkw7dyo11eg+AABAuCDYmVZCgjp3NroJAAAQRjgUCwAAECEIdgAA\nABGCYAcAABAhCHamtXWrZs7Uhx8a3QcAAAgXBDvT+uEHzZ+vrCyj+wAAAOGCq2JNa/BgZWer\nXTuj+wAAAOGCYGdazZqpTx+jmwAAAGGEQ7EAAAARgmAHAAAQIQh2kcZu14IFyshQq1YaMkTr\n1xvdEAAACBXOsTOtHTu0erUGD9Z553mWb7lFCxdWPs7N1YYNWrZMV11lQIMAACDE2GNnWt99\nd/p97D7//FSqc7vlFp04Ebq+AACAURrdHjun02m324M6f4fDEdRFVMrIsPzrX85f/lIey9q0\nyXJ6WD92TFu3Os4+2xn0loIp2D+4RohVGnBOp1OSw+EwupGIYrfb2VYDzrWtslYDy+FwWCwW\nw9dq4wp2IXjbdQW7ULyzt22rK66QJI9lRUf73gUbHe1wOEwf7Fw/PgQQESQYWKuB5frdZ60G\nA2s1sEITAGr9KGxcwc5isVgslqioIL5qq9Vqs9mCuogaDB+uuDiVlVUpdu6snj1tVtMedXf9\nsW61Wq3mfQ1hKdi/C42Qa9+SzWazWCxG9xI5XL/+bKuB5dpWWauBFZoAUOvbS2P8oQb7PdcV\nH4O6iOqceabmz9eMGacqTZrolVdks0XCxwwflgHHKg0si8XidDoNfAeISK6VySoNBtZqYFn+\nx9g2GmOwi2zTp+u88/Tii9q7V716acYMtW9vdE8AACAkCHamtXevsrJ09tnq2tVrSL9+6tfP\nkJ4AAICROGnJtDZv1sSJevNNo/sAAADhgmBnWmedpXnzNGCA0X0AAIBw40qhJwAAIABJREFU\nwaFY0+reXd27G90EAAAII+yxAwAAiBAEOwAAgAhBsAMAAIgQBDvTysvTli3KzTW6DwAAEC4I\ndqb1zjvKyNCrr54+xOHQokU6+2w1b66zztKCBTL6K4kBAEAocFWsaaWn66ab9Mtfnj7k4Yd1\n332Vj/PydOutOnBAc+eGtDsAABB6BDvTOvdcnXvu6eWDB/XnP3sXH3pIN9ygTp1C0BYAADAM\nh2IjzZdfqrzcRz07O+StAACA0CLYRZSKCv3nP74H3XefvvsutN0AAIDQIthFlEmT9PTTvgd9\n/70uvVRHj4a2IQAAEEIEO9M6eVL5+Tpxwl14802tXFnTFPv2adGioPcFAACMQrAzrRUrlJrq\nmdQ2bKh9oh9+CGJHAADAWAQ70zrjDE2YoPR0d8HprH2ili2D2BEAADAWwc60BgzQihW65BLP\nQs3i43X11cFtCgAAGIhgFznGjNHYsVUqNtupx4mJWrRIvXuHuCkAABA63KA4clgsmjJFBQX6\n8UfFxKh/f82apSNH9PnnSknR8OFq3droFgEAQDAR7CLHNdfo5ZdPPU1NVbt2OvNMnX++cT0B\nAIAQ4lBshHjllSqpTtLnn+veew3qBgAAGIFgZ1r//KcsFj3zjOtZZqaPUVavDmlHAADAWAQ7\n02reXH36KC3N9aykxMcoxcUh7QgAABiLYGdal1+u7GyNH+96ds45Pkb59a9D2hEAADAWwS5C\n3HWXOnSoUomO1v79SktT9+568EGVlhrUGQAACBWCXYRo1kzvv6+RI5WUpNhYpaervFzffacj\nR7Rjh2bP1uTJRrcIAACCjGAXOdLT9cYbKijQ8eM6csR76Ouva/16I9oCAAChQrCLQDk5OnbM\nRz07O+StAACAECLYmdbKlerSRf/61+lDEhJ8T9GkSXA7AgAAxiLYmdbJk8rP14kTpw/p2NHH\nRbJxcbrkklD0BQAAjEKwM63f/U55ebr2Wp8DX3xRSUlVKg8/rDPPDEVfAADAKHxXbGQ65xzt\n3KnFi/Xtt2rfXr/9Lfe0AwAg8hHsIlaLFpo50+gmAABACHEoFgAAIEIQ7AAAACIEwc601q7V\n8OF6+22j+wAAAOGCYGda+/dr/Xrt3290HwAAIFwQ7EzrxhvldGrqVKP7AAAA4YJgBwAAECEI\ndgAAABGCYAcAABAhCHYAAAARgmBnWuvWaeJEffCB0X0AAIBwQbAzrZ07tXKldu0yug8AABAu\n+K5Y05o0SRddpLQ0o/sAAADhgmBnWklJSkoyugkAABBGOBQLAAAQIQh2AAAAEYJgBwAAECEI\ndqa1ebNmztSXXxrdBwAACBcEO9P64gvNn6/vvjO6DwAAEC64Kta0Ro9Wt27q1cvoPgAAQLgg\n2JlW+/Zq396zsH27srNltapvX6WnG9UWAAAwDMEuQvzxj/+/vTuPq7LO+z/+OQuy74SiuOOY\n4i4ajktDimWp6KjZWJa2MelM6eTdNGmU1V15O7emlpraNFqZlcvDrTRRZlyS0m61FDETFZHE\nn4LEDmf5/XHseEIEa84518Lr+eiPrs/1Pdf5cPyCb7/XgixaJDU1IiJNmsjzz8usWUr3BAAA\nvItr7PTggw9k3ryrqU5Eqqvl+efls88U7QkAAHgdwU4Pliypo7h0qdf7AAAAiiLYadbRo7Js\nmXz/vYgUFNSx/8IFb3cEAACURbDTrF27JDVVvvpKRNq1q2N/nUUAAKBjBDvNGjZMPv5YBgwQ\nkWeeqb3Tz0+eflqBpgAAgIIIdprVoYOMGyetWonI4MGycqVERl7d07y5fPihJCQo2R0AAPA+\nHneiEw8+KOPHy6lTYjRKu3bSpInSDQEAAK8j2OmHr6907qx0EwAAQDmcigUAANAJgp1mnTwp\nn3wi584p3QcAAFALgp1mbdki994rX3yhdB8AAEAtCHaaNWCAvP66dO+udB8AAEAtuHlCs/r0\nkT59lG4CAACoCCt2AAAAOkGwAwAA0AmCHQAAgE4Q7DSroEC+/lqKipTuAwAAqAXBTrPee08S\nEmTnTqX7AAAAakGw06zu3eXxx6VdO6X7AAAAasHjTjQrOVmSk5VuAgAAqAgrdgAAADpBsAMA\nANAJgh0AAIBOKBDscnNzx4wZExISEhoaOnbs2Ly8vBuNnDt37uDBg6Oionx8fFq2bDl16tTL\nly+7DjBc5+jRo57/CtShslKKiqSmRuk+AACAWng72JWWliYlJZ04cWLVqlUrV67Mzs6+4447\nysvL6xz8wgsvtGjRYtGiRbt27UpLS1u/fv2AAQNqDZ40adJ+F+3bt/fK16ECCxZIRIRs2aJ0\nHwAAQC28fVfs8uXLz5w5c+LEibi4OBHp0qVLx44dV6xY8eSTT14/+PTp002bNnX8/8CBA+Pj\n4/v37//xxx9PmjTJOaZFixaJiYle6V1l2rSRIUMkOlrpPgAAgFp4e8Vu8+bNiYmJjlQnInFx\ncYmJiRs3bqxzsDPVOfTq1UtEzp07V2uY3W73QKeqN3687Ngh/fsr3QcAAFALbwe7Y8eOdenS\nxbUSHx+flZV1M6/dvn27Y7xrcfHixb6+voGBgYMHD963b58bWwUAANAWb5+KLSoqCg8Pd61E\nREQUFhY2+MLCwsJp06Z17do1JSXFWXzkkUfuvvvu5s2b5+bmzp07NykpaefOnQMHDnR9u9TU\nVOdmSEhIYmKi1Wp1x5dSN7vdbrPZPPoWjY3jI22k67KeZLfbmajuZbVaHZ+qwWBQuhf9cH6q\nSjeiK44fqnyq7mWz2QwGg+KfqjZ+80R5eXlKSkpxcfGnn35qMpmc9RUrVjj+JzEx8Z577unS\npUtaWlpGRoZzQGVl5SeffOLcjI2N7dOnT0lJiedarampsVgsRiPPkXEnx3cLf1m6l8Vi8ej3\nQiNkt9vtdjvf/u5lt9tramqYq+5ls9lEhLnqXhaLxWAwVFdXe/RdHH929fB2sAsPDy8qKnKt\nFBYWRkRE1POSioqKESNGHDlyJD09vVOnTjcaFhgYOGLEiHfeece1GB0dferUKddDZWZmhoWF\n/dr2G1ZSUuLn5+fj4+O5t2hs7Ha7Iyu7Znr85woLCz36vdAIWa1Wm81mNpv5R4gbWa3WkpIS\n5qp7WSwWu93OX1XuVVZWZjKZ/Pz8PPouDcZxb6f1+Pj4Y8eOuVaysrI6d+58o/GVlZUpKSmZ\nmZlbtmzp27fvL307k8nUzkVsbOyvaVqdXntNDAbZtEnpPgAAgFp4O9gNHz48MzMzJyfHsXnq\n1KnMzMyRI0fWObiqqmrUqFF79uzZuHHjoEGD6j9yaWnppk2bbrvtNjd3rFrNmknv3sK/YgEA\nwE+8fSr2scceW7RoUUpKyiuvvCIis2bNatOmzaOPPnqtIbN55syZs2fPFpExY8Zs3749LS0t\nKCgoMzPTMSA2Ntax8DZ//vyzZ88mJSVFRkaeO3du7ty5P/zww+rVq738FSlm8mSZPFnpJgAA\ngIp4e8UuODg4IyMjLi5u4sSJEydO7NChw65duwIDA50DrFar846SrVu3ishLL73Uz8XSpUsd\ne9u3b//ll19OmjTp9ttv//Of/9yyZct9+/b99re/9fJXBAAAoBIK3BXbpk2bDRs23Giv61Mt\n6n/CxciRI290DhcAAKAR4lZnAAAAnSDYAQAA6ATBTrPefFPatxeXpzEDAIBGjmCnWZWVUlQk\nNTV798pdd0mzZtKzp/zv/0pNjdKNAQAAhWjjV4qhDjNmyIwZW7fK8J9+NW5BgRw+LAcOyJo1\nijYGAAAUwoqdhtlskppau/jRR7JjhxLdAAAApRHsNOz77+X8+Trq+/d7vRUAAKACBDtNysuT\nMWMkPr7uvU2aeLcbAACgDlxjpz0VFXL33fLttzcccOedXuwGAACoBit22rN6tXz7rUyU93ZI\nch85UGvvyy9Lz56K9AUAABRGsNOe48dFRNrLqSGSHimXnfUJEyQjQ2bNUqwxAACgLE7Fak9E\nhIjIHPnrAnmqVIIcRbNZli+XgAAlGwMAAMpixU57xoyRgACpEP8iCa8RH9ciAABozAh22tOx\noyxd+rMY16ePLFmiXEMAAEAdOBWrSRMnyu23y+bNcvGi9OolI0aIkYgOAECjR7DTqlatZOpU\npZsAAABqwjqPZm3aJKmpV2+RBQAAINhp2MGDsmyZ5OUp3QcAAFALgp1mTZ8up07JgAFK9wEA\nANSCa+w0KzxcwsOdW2VlkpMjAQHSurWY+VMFAKBRYsVOD15/XaKjpVs3iYuTrl3liy+UbggA\nACiBYKd5y5fL3/4m5eVXN7OzZcQIOXdO0Z4AAIASCHaa9/rrtSuFhbJsmRKtAAAARRHsNGv7\ndnn2Wfv3p86cqWNnTo632wEAAIoj2GnW7t0yZ44h9+wtt9Sxs1kzr/cDAACURrDTrMmTZccO\n6dHjiSdq7/Hzk4cfVqIlAACgKB6MoVlxcRIXJyLPPSdHj8ratVfLgYHy1lsSH69kawAAQBEE\nO83z8ZFPPpGDB+XoUQkIkAEDpHlzpXsCAABKINjpREKCJCQo3QQAAFAU19gBAADoBMFOs776\nSpYtk/x8pfsAAABqQbDTrA0bJDVVvvtO6T4AAIBacI2dZo0eLW3bym9+o3QfAABALQh2mtW3\nr/Ttq3QTAABARTgVCwAAoBMEOwAAAJ0g2AEAAOgEwU6zTp6U9HS5ckXpPgAAgFoQ7DRr+XJJ\nTpZvv1W6DwAAoBbcFatZw4ZJZKS0bat0HwAAQC0IdpqVlCRJSUo3AQAAVIRTsQAAADpBsAMA\nANAJgh0AAIBOEOw0Ky9Pvv5aSkuV7gMAAKgFwU6z5s+XhAT55hul+wAAAGpBsNOsxER5/HFp\n2lTpPgAAgFrwuBPNGjdOxo1TugkAAKAirNgBAADoBMFO28rL5cwZqa5Wug8AAKACBDutunhR\nxo+XwEBp21aCguTpp6WyUumeAACAorjGTpOsVnlwbPlXe6rMEmwRc02NzJsnVVXy5ptKdwYA\nAJTDip0mZWTIsD3PFUpEb/naWVyyRAoKFGwKAAAojGCnSSdPynfym3QZUiyhzqLNJt9/r2BT\nAABAYZyK1aTISFksUxbLlFr1qChF2gEAAKrAip0m3XWXNG9eu9i/v3TsqEQ3AABAHQh2mhQS\nIu+/L5GR1ypms/TuzY2xAAA0agQ7rUpKkocfvrZpscjChTJ9unINAQAApRHstOrCBZk/v3Zx\n6VI5dkyJbgAAgAoQ7LSq6qlnLloiesjhWvVDhxRpBwAAKI9gp1XGoIAiCbeKqVY9MFCRdgAA\ngPIIdlrVdMmLg1qc+la6uhbDwuT225XqCAAAKIxgp1VNmsjq1RIcfK3i7y/vvisREcr1BAAA\nFMUDijVs0CDJzpZ//lNyc6V5c3ngAWnXTumeAACAcgh22ta8uTz3nNJNAAAAdeBULAAAgE4Q\n7DQrLU3at5esLKX7AAAAasGpWI2xWGTFCsnMlPsOliddKDJXWWo/7wQAADRWBDstqamR22+X\n/ftFRFbK30X+PuBJ2bVLfHyU7gwAAKgAp2K1ZO7cq6nOae9emTdPoW4AAIDKEOy0ZNOmOoob\nN3q9DwAAoEoEOy2prLzZIgAAaIQIdpqxZIl8910d9d69vd4KAABQJYKdNixdKlOmSEXFtco0\neWOHJCdE5MyerVxbAABATbgrVgOsVpk1q3axkxwfIulhc0ubN1eiJwAAoD6s2GnAxYty+XLt\n4nSZHyGFA/8YP22a2O1KtAUAAFSGYKcBwcFivO4PqlwCiiS8ssa0YIEsWKBEWwAAQGUIdhoQ\nFCQpKfUNWLzYW60AAAAVI9hpw5IlYrrx7w7Lz/diKwAAQK0IdtrQtKn4+d1wb7t2XmwFAACo\nVaO7K9Zut1ssFs8d32azWa1Wg8HgxmPa7VJeLr16mfbsuXbYB2VVf9k3W17Il+Z/+YvVYtH5\nDRQ2m83OTSJu5envhUbIMUWtVqvSjeiK1Wq12WzMVfdyzFU+Vfey2Wzi+U+1wb8KG1ewc3wc\nns4HdrvdXW9RWChpaaaPPjIWF0urVvYmTaS6+uquQbL7EXlnuf+Tqc82vf9+W2PIPAQ7t+Mj\n9QQ+VU/gU/UEPlX3sv9E2TYaV7AzGAwGg8HHx8dzb2E0Gs1ms1vewmKRceNkz56rm7m5BhHp\n3FnOnpWgIDnQ8/lOdz+RcV+7oFtMIje+/k77HAtLRqPRVM9lhvjlPP290Ag51pbMZrN71+wb\nOavVajQamavuZbFY7HY7n6p7VVdXm0wmT3+qDf54aVzBTlvWr7+W6pwuXpTiYseNFK1FWivQ\nFgAAUCtunlCvb7+to3jpkpw/7/VWAACAFhDs1Cs4uI6iwVB3HQAAgGCnXikp4u9fu5icLOHh\nSnQDAABUj2CnXh07ysKF4ut7rRIXJ//4x08bn34qc+bU8UtkAQBAY8XNE6r26KMycKCsWyc/\n/CA9e8r997vkvE8+kX/+U1JSJDJSyRYBAIBqEOzUrmNHee65unbMmCH33y8tW3q7IQAAoFYE\nO82Kj5f4eKWbAAAAKsI1dgAAADpBsAMAANAJgh0AAIBOEOw0KyNDli2T4mKl+wAAAGpBsNOM\n0tKfb7/7rqSmysWLynQDAADUh2CndhUVkpYmt9wiwcESGytvvCFWq4iITJ4sb78t0dEK9wcA\nAFSDx52oXWqqvPfe1f8/f16mT5fCQnnpJZGkJElKUrQ1AACgLqzYqdr//d+1VOf06qucgAUA\nAHUg2Knat9/WUbRa5ehRr7cCAABUj2CnasHBdddDQ73bBwAA0AKCnaolJUlUVO1ihw7So4fI\n0aOSni4VFUr0BQAA1Ihgp2rh4fLeexIScq3StKl89JGYTCKvvSbJyfLDD8p1BwAA1IW7YtXu\nrrskO1s++EDOnJGOHeXBB386DztihLRsKWFhCvcHAABUg2CnATExMmPGddX77pP77lOgGwAA\noFacigUAANAJgh0AAIBOEOwAAAB0gmCnWRcvSk6OWCxK9wEAANSCYKdZTz4p7dvL+fNK9wEA\nANSCu2I1a/BgCQ2VoCCl+wAAAGpBsNOsxx6Txx5TugkAAKAinIoFAADQCYIdAACAThDsAAAA\ndIJgpxmXLsmhQ5Kf/9N2WZkUFYnNpmRPAABATQh2GlBSIg89JNHR0quXtGghw4bJ+fMikyZJ\nRIRL0AMAAI0dd8VqwLRpsmrVtc1t22TCBPnXHV0NV66In59yfQEAAHVhxU7t8vLkH/+oXdy9\nWzIGpsmOHRIVpURTAABAjQh2apeTU3f9+++92wcAAFA9gp3aRUfXXW/a1Lt9AAAA1SPYqd2t\nt0pSUu1iXJwMHapENwAAQMUIdhrwz39K9+7XNlu3lg8/FH9/5RoCAACqxF2xqmaxSHa2lJRI\nRoYcOCC5udK0qQweLAEBIhMnytatcvw4J2UBAIADwU51vvtOsrKkVSspKZHHH5fvvhMRCQyU\nF1+UGTNcxgUFSXi4GAwKtQkAAFSHYKcily7JxImybdvVTZNJrNar/19WJv/1XxIZKZMn/zR6\nyRLvdwgAANSMa+xU5IknrqU6kWupzum117zZDgAA0BiCnVrk58vatQ2MOX1a7HavdAMAADSI\nYKcWeXkNj4mJ4Zo6AABwQwQ7tYiNbTi0TZnilVYAAIA2EezUonlzGTeuvgFPPCHPPOOy/Ze/\nSEKCFBZ6uC8AAKAZBDsVef31uutDh8qJE7J4sRhd/7gKCiQnR2w2r7QGAAA0gMedqEirVhIS\nIj/+WLs+bJj85jfXjf7gA680BQAANIMVOxUxmX7+CGIREWnWTB54QIluAACA1hDs1OVvf5M/\n/1l8fK5uduki69dLVJSiPQEAAI0g2KmL2SwLF8rFi7Jnjxw/LkeOSL9+SvcEAAA0gmvsVOTH\nH2XdOsnJkQ4dZOxYCQhQuiEAAKApBDu1yMyU3/9efvjh6ubzz8uWLdK1641fMGuWfPmlrFsn\nISFeaRAAAKgdp2JVobxc7rvvWqoTkdxcufdesVhu/JrDhyU9XWpqPN8dAADQBlbsVGHPHjl7\ntnYxO1uGDpW2bWXIEBk//ucPsRORNWukpkbCwrzVIwAAUDtW7FShqKjuekaG/OMfMmGCDB16\n3dpcUJCEh/O7YwEAgBPBThU6dWpgwM6dMmeOV1oBAACaRbBThe7d5Q9/aGDMhg1eaQUAAGgW\nwU4tli2TadPEz++GAyorvdgNAADQIIKdWgQFyfz5UlIiHTrUPaB3759vL1woqalSXu751gAA\ngDYQ7NTFbJYnn6yjHhIi//3fPy999pksWybV1V7pCwAAaADBTnWmTpXp069tGgzSu7ccOCAt\nW/583MKFcvCgBAd7tzsAAKBePMdOdQwGmTdPpk6VL78Uf38ZOFCiouoad6NTtgAAoLEi2KlU\n+/bSvr3STQAAAE0h2KmC1SrZ2XLlinTtyq9+BQAAvxLX2Clv/37p1k26dJEBAyQmRl59VemG\nAACANhHsFHb+vIwcKVlZVzfLy2XmTFmy5CZe+dFHMmcOd8UCAAAngp3CVqyQS5dqF2s/2aRO\n77wjzz4rVVUeaAoAAGgSwU5hp0/XUTx/Xg4dauiVs2bJxx+Lv78HmgIAAJrEzRMKi4mpu/7Y\nY3LwYL2vHDTIA+0AAAANY8VOYZMmSZMmddS//lry8rzeDQAA0DKCncI6dpRx4+reVVHh3VYA\nAIDGEeyU165dHcXISGnb1uutAAAALeMaO2WUlcn778s338ju3XL0aB0DFi4Uc/1/OP/+t1y8\nKKNHNzQOAAA0FmQCb6uulueekzfeEKu17gExMbJ0qYwc2dCBXn5Zdu6U0lKCHQAAcCATeNvT\nT8ubb9Y3IDLyJlKdiEyfLvfeK76+buoLAABoHsHOq77/voFUJyLl5Td3rHvu+Y/bAQAAusLN\nE15V5+V0tfTu7fk+AACAHikQ7HJzc8eMGRMSEhIaGjp27Ni8eh/XVv/gX3QoxVkssmNHA2N8\nfOT1173SDQAA0B1vB7vS0tKkpKQTJ06sWrVq5cqV2dnZd9xxR/kNzj7WP/gXHUpxly7J4MGy\neHEDw3r3rvvpJwAAAA3y9jV2y5cvP3PmzIkTJ+Li4kSkS5cuHTt2XLFixZNPPvlLB/+iQ3nH\n+fMyY4Zferr5xx/FYhGDQez2q//dpD59bvrNvv5aiorkjjvEyPl0AAAg4v0Vu82bNycmJjqi\nmIjExcUlJiZu3LjxVwz+RYfygkuXpHdvWbPG59IlQ3W12GxitYrN9gtSXViYzJhx0+83Y4Yk\nJ0t19a9qFgAA6JC3g92xY8e6dOniWomPj8/KyvoVg3/RobzghRekoOBXvtZgkNtuk61bpVWr\nm37NH/4gf/0rD7EDAABO3o4FRUVF4eHhrpWIiIjCwsJfMfhmDlVQUNCpUyfnZkxMzLRp0zx0\nHd6//+3nCMo+UlMgTb+Q3w6XLa4DsuXWCvHvKYdci7sMgweFf/P/juWGhNjlp2edNHn0UdNn\nn1UeOmSPjnaO9ElLM7/zTtXWrbYePUREHnhARKS6WveLdna73WazGQwGIyed3cpqtar2mlSN\nstlsdrvdaDQaDAale9EPm81msViYq+7lmKs1NTVKN6Ir1dXVRqPRZrN59F3sDZ0H1Pl6j8lk\naudyM0JkZKSj6Jn3uvqj3C6GHGn3g8TUGnBWWldJ7ecJh3RqbgwsCQ//WWQxNGsm7dqZmjSx\nu7RqiIqSdu2M/v4Gz/SvWo5JbDQaCXbuZTAYPPS90GgZDAabzWYymQh2bmQwGJirnmC32/lU\n3cvx95Tin6q3g114eHhRUZFrpbCwMCIi4lcMvplDRUVFHTx40LlZXFy8du1aX8/8toYhQ+Sb\nb0RELGJOkIPXD7hTtrtumkzyzDPSbfZ7Bp/r4t68eSLSpFbx2Wfl2WdrFxsBu91usVjU8N2i\nM2VlZR76Xmi0rFarzWYzm80EOzeyWq1VVVXMVfeyWCx2u93Hx0fpRnTFYrGYTCZPz9UGf7x4\newkkPj7+2LFjrpWsrKzOnTv/isG/6FBe8OKL4nLWtwF33y0//iivvip8WwEAAHfxdrAbPnx4\nZmZmTk6OY/PUqVOZmZkjb/C7Uesf/IsO5QXBwXLggMycWRUfbw8JET8/8fOTwEAJDZXoaOna\nVYYOlQkT5PXX5fBh2bpVAgKU6hQAAOiTocGr8NyrpKSkW7duQUFBr7zyiojMmjWroqLiyJEj\ngYGBjgFms3nmzJmzZ89ucHCDh7qe41TsI4884tEv0M/Pj/VtN+JUrIfUcxUEfh1OxXqC1Wot\nKSkJCwtTuhFd4VSsJ5SVlZlMJj8/P4++y6JFiyoqKp555pkbDfD2il1wcHBGRkZcXNzEiRMn\nTpzYoUOHXbt2uUYxq9VqtVpvZnCDhwIAAGhUFLgrtk2bNhs2bLjR3loriPUPrn8vAABAo8Lz\nIwAAAHSCYAcAAKATBDsAAACdINgBAADoBMEOAABAJwh2AAAAOkGwAwAA0AmCHQAAgE4Q7AAA\nAHSCYAcAAKATBDsAAACdINgBAADoBMEOAABAJwh2AAAAOkGwAwAA0AmCHQAAgE4Q7AAAAHSC\nYAcAAKATBDsAAACdINgBAADoBMEOAABAJwh2AAAAOkGwAwAA0AmCHQAAgE4Q7AAAAHSCYAcA\nAKATBDsAAACdMCvdgLcVFBQsW7bMc8e3Wq1Go9FgMHjuLRohu90uInyq7mWxWMzmRvcTwKOY\nqJ5gt9ttNpvJZFK6EV1hrnqCzWYTEaPRs0tmxcXFTZo0qWdAo/uxXl1dnZ+f76GD19TUzJs3\nr02bNuPHj/fQWwBucfz48U2bNg0ZMqR3795K9wLUZ/PmzVlZWampqWFhYUr3AtRn/vz5ISEh\njzzyiEffxWg0xsTE1DPA4IjtcIvS0tLg4OAhQ4bs2LFD6V6A+nz44YcTJkx44403nnrqKaV7\nAepz//33r169+tSpU+3atVO6F6A+ISEhrVq1Onr0qLJtcI0dAADUht65AAAKWUlEQVSAThDs\nAAAAdIJgBwAAoBNcY+dOdrv9ypUrPj4+QUFBSvcC1Ke6urqsrCwgIMDX11fpXoD6lJWVVVdX\nh4aGevpmQ+A/dOXKFaPRGBISomwbBDsAAACd4B9AAAAAOkGwAwAA0AmCndvk5uaOGTMmJCQk\nNDR07NixeXl5SneERmru3LmDBw+Oiory8fFp2bLl1KlTL1++7Dqg/rnKTIaX2e32QYMGGQyG\nWbNmudaZqFCVVatWJSQkBAQEBAcH9+vX79ixY85dqpqrXGPnHqWlpd27d/f393/llVdEZNas\nWdXV1YcPHw4ICFC6NTQ6AQEBY8eOHTZsWGxsbHZ2dlpaWlhY2Ndff+2YjfXPVWYyvO/NN998\n7bXX8vPzZ86c6Zh4wkSFyjz33HNLliz561//2q9fv6qqqszMzLvvvjshIUFUOFftcId58+YZ\njcaTJ086Nk+ePGk0GhcsWKBsV2icLly44Lq5b98+EXn33Xcdm/XPVWYyvOz06dNBQUHr168X\nkZkzZzrrTFSoxxdffGEymb744os696ptrrJi5x533HFHVVWV429Qh/79+/v5+e3cuVPBrgAR\nqays9Pf3f+mll55//nlpaK4yk+FlycnJwcHB69evNxgMrit2TFSox4MPPnjmzJndu3fXuVdt\nc5Vr7Nzj2LFjXbp0ca3Ex8dnZWUp1Q/gtH37dhGJj493bNY/V5nJ8KYVK1YcOHDgrbfeun4X\nExXqsXfv3q5du86cOTM6OtrHx6dHjx4bNmxw7lXbXCXYuUdRUVF4eLhrJSIiorCwUKl+AIfC\nwsJp06Z17do1JSXFUal/rjKT4TXnz5+fMWPG//zP/8TExFy/l4kK9cjPz3/vvfeOHTu2evXq\nTz/9tGXLlr///e+3bt3q2Ku2uWr23KEBKKu8vDwlJaW4uPjTTz81mUxKtwP8zB//+Mfu3bs/\n9thjSjcCNMBms7Vq1WrdunWOH6RJSUndunV79dVX77nnHqVbqwPBzj3Cw8OLiopcK4WFhRER\nEUr1A1RUVIwYMeLIkSPp6emdOnVy1uufq8xkeMfatWu3bdu2d+/e4uJiZ7GqqurKlSshISFG\no5GJCvWIjIxMTk52/vPYbDYPHTp0+fLljk21zVVOxbpHfHy86yNtRCQrK6tz585K9YNGrrKy\nMiUlJTMzc8uWLX379nXdVf9cZSbDO44ePWqxWBITE8N/IiJ///vfw8PDHZcfMVGhHs5rlG+0\nV1VzlWDnHsOHD8/MzMzJyXFsnjp1KjMzc+TIkcp2hcapqqpq1KhRe/bs2bhx46BBg2rtrX+u\nMpPhHZMmTcr4ORF54IEHMjIy2rZtK0xUqMno0aM///xzq9Xq2LRYLNu3b+/Tp49jU21zlced\nuEdJSUm3bt2CgoKcTyCsqKg4cuRIYGCg0q2h0Rk+fPjWrVvT0tKGDRvmLMbGxsbGxkpDc5WZ\nDKXUetwJExXqUVlZ2atXr/bt2//pT3+y2+0LFy78/PPP09PTf/e734kK56rnHpHX2Jw+fXrU\nqFHBwcHBwcGjR48+e/as0h2hkarzO9310a/1z1VmMhRRa5bamahQk/z8/AkTJoSFhfn6+vbr\n1y89Pd11r6rmKit2AAAAOsE1dgAAADpBsAMAANAJgh0AAIBOEOwAAAB0gmAHAACgEwQ7AAAA\nnSDYAcDNMhgMjuc8A4A6EewAAAB0gmAHAACgEwQ7AAAAnSDYAdCzzMxMg8EwevTo63d16tTJ\n19e3sLDQsWmz2RYvXpyQkBAYGBgYGNinT5+lS5fabLb6j1/nVXexsbEGg+H6YVeuXPnTn/7U\nokWLJk2a9OzZc8eOHSJy5cqVKVOmNG3a1NfXt2fPnp9//vn177J169Y777wzIiLC19e3Y8eO\nL774YkVFxc1/CAAaD35XLACdu/XWW0+fPp2fnx8ZGeksfvXVV7fddtuYMWPWrl3rqIwfP/7j\njz9u3br1qFGjRGTDhg25ubkTJkz44IMPnK8yGAwtWrTIy8urpyIisbGx58+fd/3pajAYoqOj\nW7VqZbVahw4dWlBQ4Djsv/71rylTphiNxuTk5IsXLzqKBw8e7Natm/O1aWlpL7/8cqtWrYYP\nHx4SErJnz559+/YNHDhw586dPj4+7v2sAGieHQB07dVXXxWRRYsWuRanTJkiIps2bXJsOhJV\n3759S0tLHZWSkpKEhAQRWbNmjfNVItKiRQvX41xfsdvtLVq0qPXT1fHzdsKECTabzVFZsGCB\niAQFBV1ffPTRR50vdKzq3XXXXeXl5c7i9OnTRWTevHm//MMAoHOs2AHQuby8vNatW/fq1evA\ngQOOSnV1dUxMjNlsPn/+vNlsFpEhQ4bs3Llz7969/fv3d75w9+7dt99+e3JysvP06H+yYici\nubm5LVu2dFQKCgqaNWsmIufOnXOezL1w4UJMTEz37t0PHz7sqIwaNWrjxo05OTlt27Z1Hu3H\nH3+MjIzs0aOH8ysCAAez0g0AgGfFxsYOHjx4x44dWVlZnTt3FpHNmzcXFhZOnz7dkepE5NCh\nQ/7+/v369XN9Yf/+/X19fQ8dOuSWNqKiopypTkSio6MdRddL9Jo2bSoiRUVFzsr+/fv9/PxW\nrlxZ62i+vr7Z2dluaQyAnhDsAOjfpEmTduzYsXLlyjlz5oiIIyc99NBDzgHFxcW33HKL0fiz\n+8lMJlNYWNjly5fd0oOvr6/rpmMNr86i1Wp1VgoLCy0Wy+zZs93SAwDd465YAPo3evTokJCQ\n999/32q1Xrx48bPPPuvevXv37t2dA0JDQ4uLi2vdA2u1WouLi0NDQ+s5ssFQxwUtVVVV7uo8\nNDQ0PDz8RhfTuOtdAOgGwQ6A/vn7+9977735+fnp6emrV6+2WCyuy3Ui0rNnz4qKiszMTNfi\n/v37Kysre/bsWc+Rw8PDL1++bLFYnJVz585dunTJXZ3fdtttRUVFx48fd9cBAegbwQ5AozBp\n0iQRWbVq1apVq8xm8/3333/93hkzZpSXlzsqZWVlTz/9tIg8/PDD9Rw2ISGhqqpq3bp1jk27\n3f7ss8+6se2nnnpKRFJTU69cueJaLysr++abb9z4RgD0gbtiATQWHTp0OHv2bE1NzYgRIzZt\n2uS6y263jxs3bt26dW3atHE+x+7s2bPjx49fs2aNc9j198Bu27Zt2LBhfn5+EydOjIiI+Pzz\nzwMDA0+ePFlQUFDrrtjrb569yeJLL730wgsvREREDB8+vGXLlhUVFVlZWbt37548efKbb77p\nps8GgE6wYgegsXjooYdqamrk57dNOBgMho8++mjRokURERFvv/3222+/HRUV9dZbb61evbr+\nY951111r1qyJi4tbuXLlqlWrBg0atG3bNufNtm6Rlpa2c+fOQYMGbd++fc6cOcuXL79w4UJq\naurUqVPd+C4A9IEVOwAAAJ1gxQ4AAEAnCHYAAAA6QbADAADQCYIdAACAThDsAAAAdIJgBwAA\noBMEOwAAAJ0g2AEAAOgEwQ4AAEAnCHYAAAA6QbADAADQCYIdAACAThDsAAAAdOL/Ax9Bp85h\nPXy6AAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# have a look at what the label and associated prediction looks like\n",
    "ix = 14;\n",
    "x  <- seq(1, 600, 1)\n",
    "y1 <- model_predictions[,ix]\n",
    "y2 <- labels[,ix]\n",
    "df <- data.frame(volume=x,prediction=y1,label=y2)\n",
    "\n",
    "# basic graphical object\n",
    "ggplot(df, aes(volume)) +\n",
    "  \n",
    "  # the predicted CDF using the trained model\n",
    "  geom_point(aes(y=prediction), color=\"blue\") +  \n",
    "\n",
    "  # the label for this set of images (i.e. systole volume)\n",
    "  geom_line(aes(y=label), color=\"red\", linetype=\"dotted\") + \n",
    "\n",
    "  # make the plot theme simple \n",
    "  theme_bw()\n",
    "\n",
    "# evaluate the cost function for this prediction\n",
    "costfun(labels[,ix],model_predictions[,ix])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the prediction(blue curve) is a logistic like curve (CDF) while the label is a step function.  Certainly we can use our cost function to get a sense for how \"close\" the set of predictions are to their labels.  However, notice with this formulation of the problem it's not very easy to look at the prediction and call out what the actual volume predicted is.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.0155942614742691"
      ],
      "text/latex": [
       "0.0155942614742691"
      ],
      "text/markdown": [
       "0.0155942614742691"
      ],
      "text/plain": [
       "[1] 0.01559426"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "costfun(labels[,sample_index],model_predictions[,sample_index])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## OK Now it is your Turn!\n",
    "\n",
    "Lets use what we've learned from creating the ```systole_network``` and create a network to predict diastole volume.\n",
    "\n",
    "Step 1. Create a training data set using the diastole labels \n",
    "\n",
    "Step 2. Create network \n",
    "\n",
    "Step 3. Train the network\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# step 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# step 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# step 3\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discovery Requires Experimentation ...\n",
    "\n",
    "There are many ways to explore and even improve this model:\n",
    "\n",
    "1. Maybe try removing batch normalization layers ([ref](https://arxiv.org/abs/1502.03167)) [hint: modify get.lenet()]\n",
    "2. Try increasing or decreasing the number of features in the convolution layer [hint: modify ```num.filter``` in get.lenet()]\n",
    "3. How does batch size effect training?\n",
    "4. Have a go at modifying the learning rate and momentum of the training phase [hint: learning.rate, wd, and momentum]\n",
    "5. Notice that the CRPS function uses the residual squared.  What else might we try (hint: try abs in costfun)\n",
    "6. Try using different activation function (i.e. other than relu).  How does this effect performance?\n",
    "7. Maybe try using different pooling functions (i.e. other than max)[hint: use help to see other functions]\n",
    "8. Maybe try removing dropout layer or modifying the percentage dropout (i.e. default is 50%)\n",
    "9. Maybe try double differencing the data or not differencing the data at all\n",
    "10. What are other ways we could formulate the network output (i.e. other than discrete CDF)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NDSB Competition Results\n",
    "\n",
    "There were nearly [200 participants](https://www.kaggle.com/c/second-annual-data-science-bowl/leaderboard) in the NDSB-II with a variety of different approaches to solving this challange. The CRPS scores for the leaderboard top ten ranged from 0.009485 in first place to 0.012611 in tenth place.  To get a better feel for how your results stackup against the competion, we can create an empirical CDF from the leaderboard CRPS scores\n",
    "<img src=\"img/fig10.png\" alt=\"CRPS ECDF\" style=\"width: 500px;\"/>\n",
    "Therefore, achieving a CRPS score of 0.03 is at about the 80th percentile w.r.t. the overall competition results.  \n",
    "\n",
    "So how much better is a CRPS of 0.009485 than 0.012611?  To get a loose sense for how CRPS effects ejection fraction calculations, the competition provided this infographic to performers\n",
    "\n",
    "<img src=\"img/fig11.png\" alt=\"CRPS ECDF\" style=\"width: 700px;\"/>\n",
    "\n",
    "Analysis of the competion results seem more favorable than this graphic might imply and noted that top performer models agreed well with previous studies of human performance.  For a medical perspective on the competion results check out this [kaggle blog](https://www.kaggle.com/c/second-annual-data-science-bowl/forums/t/19839/a-medical-perspective-on-the-quality-of-the-left-ventricular-volume-and) by Dr. Andrew Arai of the NIH.  Additional analysis of the clinical applicability provided by Jonathan Mulholland of Booze Allen Hamilton [here](http://www.datasciencebowl.com/leading-and-winning-team-submissions-analysis/). These analysis conclude, based on the top four models submitted to the competion, that \"*The models keep the diagnosis categories pretty tightly grouped together. While the models are not right 100% the time, there is a very low probability of a severely abnormal EF being incorrectly categorized in the mild to hyperdynamic range. The normal to mild diagnoses are very likely to stay within their domain ... This is a pretty good sign pointing towards suitability for clinical applications*\"."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}