{
 "metadata": {
  "name": "",
  "signature": "sha256:82dc553075965f93218083b1b984bc8eb49c721f0b6215409d0977596e81483a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Setup"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compatibility with Python 3\n",
      "from __future__ import print_function  # print('me') instead of print 'me'\n",
      "from __future__ import division  # 1/2 == 0.5, not 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There is a four-dimensional (X, Y, Z, t) BOLD image in the same directory as this notebook, called `ds107_sub012_t1r2.nii`.   Import the `nibabel` module, and load the image with nibabel to create an image object.  Turn off memory mapping with the load command."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load image object using nibabel, turning off memory mapping.\n",
      "fname = 'ds107_sub012_t1r2.nii'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the usual way get the array data from this image and show the image shape."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get image array data from image object\n",
      "import nibabel as nib\n",
      "img = nib.load(fname, mmap=False)\n",
      "data = img.get_data()\n",
      "#data = np.asarray(img.get_data())\n",
      "type(data)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "numpy.ndarray"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "V = data[...,9]\n",
      "V.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "(64, 64, 35)"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Select the 10th volume (time index 9) from 4D image data array, by slicing over the last dimension.  What shape is it?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the 10th volume and show shape\n",
      "Vstd = V.std()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Get the standard deviation across all voxels in this 3D volume:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Standard deviation across all voxels for 10th volume\n",
      "stdlist = [V.std() for V in data.transpose(3,0,1,2)]\n",
      "plt.plot(stdlist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "[<matplotlib.lines.Line2D at 0x33a6510>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXlYU2f2x0/CIuDCUhcQtFAQNSgGyqIt1qhY2lKxuFVb\nqa22tupMO+24lP7qVKqjVtuZqdNlxhmdOnasS9Vq64qSuNQFKCAKLogwgIAKGkQEAsn9/XHmNpdw\nsxJyIzmf58kD2e492b733O973vOKGIYBgiAIomsjFjoAgiAIovMhsScIgnAASOwJgiAcABJ7giAI\nB4DEniAIwgEgsScIgnAADIp9U1OTW2xs7DmpVJonkUgKU1NTVwMAnD9/fsSoUaPOhIeH5yclJe2r\nr6/vyX1eWVnZwB49etz/7LPPft+ZwRMEQRCmYVDs3dzcmuRy+di8vDxpfn5+uFwuH3vq1Km4119/\n/Z9r165dkp+fH56cnLxn3bp1i7nPe++99/6UmJi4v3NDJwiCIEzFqI3j4eHxAABApVK5qtVqJ29v\n77tFRUWDRo8efRIAID4+/uiuXbumsI//4YcfXnjssceuSySSws4LmyAIgjAHo2Kv0WjEUqk0r1+/\nfjfHjh0rDwsLKwgLCyvYu3fvJACAnTt3TisvLx8AAHD//v0ea9euXbJ8+fLlnRw3QRAEYQbOxh4g\nFos1eXl50rq6Os+EhITDCoVCtmnTpjlvv/32+hUrVixLSkra5+rqqgIAWL58+fJ33333zx4eHg8Y\nhhHp26ZIJKIeDQRBEBZgSFuNPdHky8cff7xs3bp1i7i3XblyJTQ2NvYswzAwevToE4GBgSWBgYEl\nXl5ed318fGq//PLLBbrbwd3aFx999JHQIbTDHmNiGPuMi2IyDVNjKi5mmJEjOzcWFnt8nxjGeFzP\nPsswP/1km1hY/qedZuk2ezGY2dfU1PR2dnZu9fLyUjY2Nrqnp6dP+Oijj9Ju377dp0+fPrc1Go14\n5cqVH7711lt/AwA4ceLEU+xz09LSPurZs2f9ggULvrLoKEQQhGDcvQtQViZ0FPZNQwNA9+5CR2E6\nBj37qqoqv3HjxmVIpdK82NjYcxMnTvxx/Pjxx7Zu3frS4MGDrwwdOvRSQEBAxauvvvqNjeIlCMIG\nqFQASqXQUXQeNTUABQUd28bDJvYWnQ509AJ2aOPI5XKhQ2iHPcbEMPYZV1eM6b33GOb6devEwmJq\nTMePMwwAwzQ3W3f/uly5wjCPP25aTNbkyy8ZZvhwhtFo9D/G2Hs1ZAjDFBRYNy5jQAdsHBEjQD97\nkUjECLFfgniY6NsX4IMPAH73O9vv+9gxgPh4gFu3APr06bz9LF0KsHEjZtq2ZNkygJUrAc6cARg5\n0rJtDBwIcPIkwKOPWjc2Q4hEIosHaKldAkHYIffuAdy+DXD0qDD7b2nBv3fvdt4+WlsB/v1vtIts\nnftVVQGEhgJs2GD5Nh42G4fEniDskOJigAEDMHNUqWy/f1bsO9O3P3QIICgIoFs3FE5bUl0NsGQJ\nwO7dlr9GEnuCIDrMtWsAUVEAgwYBnDtn+/3bQuz/9S+A114D8PICqKvrvP3wUVUFEB4O8PTTADt2\nmP98tRrfIzc368fWWZDYE4QdUlwMEBKCvnl6uu3339k2zr17+LpefBHF3taVP1VVAH5+eEC9etX8\n5zc0AHh4AIgsm94kCCT2BGGHXLsGEBwMMGGCML59Z2f2x48DxMQA9OoF4OlpW7FXq3E8pF8/HHy+\nfdv8bTxsFg4AiT1B2CVsZv/kkwD5+bb3tNlxgs4S4fR0PJAB2D6zr6nBfbq4oNhbUglEYk8QhFVg\nM3s3N7QbKipsu39LbJymJoDt2017rJBiX12N7ykAQO/elNkTBCEQjY0oQAMG4HU/P/SYbUlLC4CT\nk3kifOECwFtvGS+jrKjA1yeV4nVbD9BWVQH4+uL/ZOMQBCEYJSU4UcfJCa/37w9QWWnbGFpaMOs1\nR+yrq/HxxsTz6FGA8eMBxP9TH1t79uzgLIDlYv/gAYk9QRAd5No19OtZ/PyEEfs+fcyzcW7exL+X\nLxt+HNfCARDWxunZE19rY6N526DMniCIDlNcjH49S//+wtg4ffqYJ8Ks2F+5Yvhxp08DjBmjvW5r\nsefaOCIRnsGYO0hLYk8QRIfRzeyFsnH69jXfxvHzM5zZq9X4Wrj9ZGzt2XMzewDLrBy2zv5hgsSe\nIOyMggKAIUO014UaoLXExpHJDGf21dUAPj4Arq7a24TM7AEosycIQgAYBiAvDyAiQnubEJm9SqW1\ncUxtUlZdjfaMocy+rAy7RXIRcoAWwPLMnsSeIByYmhrLpt+zlJYC9OjRtq2wUJl9z55YEWTq4OXN\nmwBPPIGllc3N/I8pL9eWlLIIOUALQGJPEIQFbN4MsHCh5c/PzdXWn7P06gWg0QDU13csNnNoacEZ\npt7epls51dUAAQHox1+7xv8YfWJvK8++vh7PVHr00N5GYk8QBC9nzmAvdj5KSgBOnAC4f9+ybefm\ntrVwALBixNbZPSv2pmbdTU148fLC8QZ9Vg6fjWPLzJ7N6rkNzMizJwiCl8mTAXJy+O8rKUEhyciw\nbNt8Yg9ge9/eXLG/eRMbi4lEAIMH6x+k5cvs3dww225q6njcxqiowPeSC2X2/6OpqcktNjb2nFQq\nzZNIJIWpqamrAQDOnz8/YtSoUWfCw8Pzk5KS9tXX1/cEAEhPT58QFRWVHR4enh8VFZUtl8vHdvaL\nIAhbce8eZof6suzSUoAZMwAOHLBs+/rEXqjM3lQbp7oaxR7AcGbPJ/YAtsvuL1/GgxEXEvv/4ebm\n1iSXy8fm5eVJ8/Pzw+Vy+dhTp07Fvf766/9cu3btkvz8/PDk5OQ969atWwwA0KdPn9s//fTT8/n5\n+eGbN2+enZKSsqXzXwZB2AY2Y+XLshkGxX7BAoCDB81fau/WLRSRwMD29z0MmT1bzhgYiHYNH3w2\nDoD+/TQ24niFtbh0CWDo0La3WdIMrUuKPQCAh4fHAwAAlUrlqlarnby9ve8WFRUNGj169EkAgPj4\n+KO7du2aAgAglUrzfH19qwEAJBJJYWNjo3tLS4tLZ70AgrAlbMbKl2XfugXg7g4QHY19XwoLzdt2\nXh4OzvItiGHvnj03s+/bF98LXZqbcVvs47joG6R95RWAI0fMix0A3/vnngMoKmp7O5/YW9LmuMuK\nvUajEUul0rx+/frdHDt2rDwsLKwgLCysYO/evZMAAHbu3DmtvLy83cnZrl27pjz++OO/uLi4tFg7\ncKJzWb8e4OuvhY7C/rhyRX+vmpISXFNVJEKh+ekn87atW1/PxdaZvUqFE5+8vc3z7AH0iz3rl4t5\nVEdfrf2dO9o2DKaydSvW+zc2Anz5Zdv7+MT+kUdw32q16ft4GMXe2ZQHicViTV5enrSurs4zISHh\nsEKhkG3atGnO22+/vX7FihXLkpKS9rm6urZZFrmgoCDs/fffX5Oenj6Bb5vLly//9X+ZTAYymawD\nL4OwNmVlwix0be9cuQIwdqxhsQcAmDIFYPFigKVLTd/25csAI0fy3ydkZm/Kfm/e1HrhrHi2tgI4\ncxRGn4UDYNjGMdfLX7MG15V97DGAxx8HWLUKWxvU1eH4g24MTk64/9paPFCZgq26XioUClAoFFbZ\nlkliz+Lp6VmXmJi4Pzs7O2rRokWfHj58OAEA4OrVq6H79+9PZB9XUVERMHny5N1btmxJCQoKKuHb\nFlfsCftDpbKs9WtX58oVgLffBvjii/b3lZZqxX7MGIAbN9r3uTFEcTHAyy/z3yekZ3/pkvHHV1cD\nPPUU/u/khC0RamratiXQNzgLYFjszanBv3MHD7pxcRj/yJEo/K++igfTIUP4zyzY8ktTxd5Wmb1u\nIpyWlmbxtozaODU1Nb2VSqUXAEBjY6N7enr6hIiIiNzbt2/3AUCLZ+XKlR/Onz//awAApVLplZiY\nuP+TTz5ZOmrUqDMWR0YISnMzib0uGg2Kt0zGn+1yM3snJ4CpU/lXblKpAP797/a363a75BIQgDaI\nNQcrDcGtxjF3gBYALR1dK8eY2POJurli//PPKPAu/xslfPNNgL/9Df/ns3BYzK3IeRhtHKNiX1VV\n5Tdu3LgMqVSaFxsbe27ixIk/jh8//tjWrVtfGjx48JWhQ4deCggIqHj11Ve/AQD44osvflNcXByc\nlpb2UURERG5ERERuTU1N705/JYRVsYfMvrUV4J13bCdwxigrw4w1KAgzyBadkaiSkraVNDNmAGzb\n1n47J04AzJvX9nU1NmJmqU8Me/bEfZeWdvRVmAY3sze39BIAM2Rdr91SG8ccsT95EmD0aO31557D\nOE6fNi72fOMMfDAM2jgPW9dLozbO8OHDL+Tk5ETq3v7OO+98/s4773yue/uHH3648sMPP1xprQAJ\nYbCHzL6qCgeKp0/HhbeFhq3RdnJCMauubivO3MweAPvEKJVYGSKRaG8/fBjf31u3tNlwSQkKIbs6\nFR9hYbitxx6z7uvioyOllwD8g7QVFQATJ/I/39OTf51dc8X+xAn07FmcnABSUwFWrMDX8+qr/M8L\nCgK4ft20fTQ14eC1oc/KHqEZtAQvzc2YaZpbK25NbtzAv99/b/pzGAb99O++A/jvf60bz5Ur2kFI\nXQ9drUabgtunXSwGmDatvZVz6BBmhdxadEMWDktYGLY/tgXmiP3t21iB1KuX9jY+sa+sbD97lcUa\nmX1DA66DGxvb9vbZs/F9k8v1Z/ahoaY3sHsYe9kDkNgTelCp8Advy0UldKmsxB/h99+bbuVcuQKQ\nlgawaxf+6M0pp9NHaSnA3r0Ax4+3FXuub19ZiYN8bm5tn/vii2jlsAfNGzfwsRMmtD0Y2aPYs6WX\nxmycU6fwLIY7P0AIsT93DmDECJzrwKVbN6yKamrSP1huqMWDLg+jXw9AYk/ogW1RK6SVc+MGQHw8\n+tWZmaY95/hx9Gm//x4F59y5jsexdCnAypV4psMup6eb2RcXt7VwWGJi8MB5/jxeP3IEX1NQkPmZ\nvURiO7FXqTCz79ULW0QYOtieOtXWJwdoL/YtLXjQ4LZu5sJXWtrSggdrU0svT5xoHwfL3LkA//mP\nduBWFxJ7wmGxF7Hv3x+rWky1chQKrSBPnGj+xCZdGAa3+f33KCZs+2HdiVVZWVjTrYtIhNk9a+Uc\nOgSQkIB2j7mZvUSC4wa2GLBmbRxnZxQ2Q+2VT57EUkcuutU4N2+i0OvzuQMC0AbjwvbRNzWzZw+k\nfLi54diPPnx98Tt/547+x9y7h58TiT3RpVCpcHKMkGJfWQng7w8waRKKpDEYBjN7tiz5+ecBfvyx\n7WMKC/kHAvVRWIg/bK4XD9A+sz97Vv+EKNbK+etfcXDWUrHv1Qs/kxLemSvWhRV7AMNWzv37eLYR\nE9P2dt3M3pCFA4AHgvv32y6U0tionQxljLt3AS5e1J/ZG0MkMu7bf/cdQEoKib3DUFWFMyO7Os3N\nKLRCZ/b+/uizlpbyDxYzDB4IlEqsgXdy0topMTGYUXLLFd9/Hyt8dKmt5Y9BLscZs7pwPXuGwR73\n+sReKkWRlssx+/T3x8ob1sZRq1H4TamyYStyOhuu2BsapD13Dl+f7liFbuml7lKAuohE+L5wD8SN\njfi+tbYan8197BieXejGYQ7GxL60FA/qpaUk9g5Berp51SGdzV//CvDVV5gZWpPmZjy1Fjqz798f\ny/KcnNoLzrVrAOPHA8yahTNPWQuHHSh0cmrbo0atRivm5Mm222lowNfKVv9wMST2bGZfXo7b5utW\nCYDxZGcD7N6tzYC5mX1FBQ7u6g4s8mGNQdpNm7BCxdByg6aKvW5dOwub2bMHaGOZPQCWseqKvbs7\nfv7GsnvWHusIxnz70lI8mOzYQWLvEPz8M3p39oBGA/DBB5hdRUVZd/BOpRJe7NnMHgCFQLdtbloa\nNg6rqMDT+NRUrV/PkpQEsGcP/n/+PE5MunABJ8WwFBZipcauXW2fq9GgLWRM7M+eBRg1ir9bpT4e\neQQPqPX1plk4LPoGaRcs4J+ty8f+/fidiY/H94JvDMCYjaNSoW1y8GB7vx5AK4YNDfjXFLHX9e1N\nFXuG0dpjHWHwYOOZ/axZmPCR2DsAp0/jF0/I+nOW4mLMCDdvRtG3ZrshoW2c+no8fff0xOsDBrQf\nwLt2DSA5GbOt7dvxsePHt33Mc89hN8nycsz8n3kGYPjwtlU6BQU4QLdjR9vnXriABwf2gMOld28U\no6IiwxaOPkQirZWTk6O//luX4cMB8vPb3sYwAD/80D5+fWRnA+zbB/Dss/j+9evX/mzVWGafnIxj\nKcHB2jESLiJRW9/emI0DoG0JwWKq2F+6hGdxuouSmEtoqPHMfsECPDiS2HdxlEr8wJ2dDZ8C2wq2\n/zkALnL988+40pE1ENrGYQdn2WyZT+yvX9f63AMGoPjr+t5ubjixacsWFHuZDG0HrpVTWIg9VC5d\nais2J060P1NgEYsB1q7FlggnTpgv9gBo5ZSVYUngtGmmPWf4cDzAcM9MiopQnI8e1VZR6ePWLTwz\nHTQI4MMP8T07dAjgt7/F1sAAePBobdUv9i0teMaTk4PP4S7ezYUr9qbaOPoye302UnExllVOm2be\nmRUfoaH4fvCd6TQ24tlNeDiexZHYd3HOnkW7RF/TJluTl4eTSACwaiE1FeAPf7DOtoW2cdiyS5aB\nA9sKwf37KFrcKfr6fuyzZwN88w0K/Jgx7cW+oADtoEmT2ma42dntZ2NymT8fBTsvD78X5jJwIFoq\nNTX82TEfbm5o5eTlaW87fhzPWMLC8MBjiOxsjJX7Xj3+OFoTb72FsbS0YELDPkbXxrl4EV83e9al\nD275ZWdk9oWFeJCdPh3gk08Mb9sUevTA18q3ylZZGR6MxGLM7qOjO74/W0NibwY//4wzBT097cO3\n52b2AABvvIGZvTWy+862cTQawxUWbGbPouvZX7+OVTd87Wp1YbNuPz8UoCefRBuntRVvLyhAoZw+\nva0Vkp3NXzvPIhIBbNyIC2Toy24N8eijAP/8Jw4um9NnJToa6/pZ2DOW55/Hwei6Omy0xlcbn5XF\nf2AaNgwz8bt321o4AO0ze0NlplysldnrS6727QN46SWAd9817XtgCpMn4/Z0s/vSUu0A/EsvmX4m\nZk+Q2JvB6dMoFL162Udmf/58W7F3cwP4/e9xsYaOwrVxGAZ/WNZcEHr/foDXXtN/P3dwFqC9EFy/\nbvqgpkiENsULL+B1Hx8U2pwcPEO4fRsPHGPHYtZ65w4OLJaU4EHAEN7emBFbwsCB+D6npJj3PK7Y\ns3MLxoxBsd+3D89Qvv0WBy11YTN7PtzdUWCtKfY3bxqfPcuim9k3NRnO7NnXbU0++wzLcP/v/9re\nzhX7hxUSexNRq3HK/siR9pHZ376NQqU72WfePDyVN2XBCX2o1SgivXqhUF67hrNY2VPllhbTJjkZ\ngl1kQh+6mSCf2JvT/fG3vwX44x+115OT0ccvLNR2suzWDStL5HI8awoLw/4wnUVYGFpKxg4ounDF\nni25DQ5GP1+jQWtr9Wo8oHJhGHyePgvCw4Nf7HVtnLNnDdtbLP7+OJ5gbPYsCzuxih2PMGTjtLZi\n8sUumGItXF2xKuvHH7Fa6Zdf8HYSeweithYzZx8f+8js2axe16fu3h1XUvr0U8u33dyMwicS4Q9w\n+XKAxESADRvwILN8OZY0dqTJWHOz4eXudDN7tg6ePb02V+wB2r5X8+bhwOi5c23bD8fHo3/9yy+W\n+fDmEBlp3GPnQyLBg6FSiZOJ2LkFIhFmu//+N7aKOHCgrR3Bvn/6eua7u6PQGsrsa2tx36YcoCZP\nxgZyV68at3AAtBOr2PkOhsQ+NxdfR+9OWCmjTx/c/vTp+H24e5fE3qGorcXaaAD7yOy5g7O6PPkk\nZuOWwi42DYBf/O3bAdatw8qTWbNwUk7PnvyTkEyluRn7wesrYdUdoGV/9KwHbE5tOh8BAehzr1rV\nVrgmTNCKvSG/XkicnHBAec8ePPDOm6e9LzAQP7vHHsPva3a29r4zZzCr1zeQbYqNk5mJ2zBljKF/\nf3yP//Qn44OzLNwzOENi3xkWDhcXF3xfn3kGLTESeweipkabRdhDZv/zz/rF3tlZO/hoCWxmD4Bi\nP3EitixITcUqln/8A+vCO9KjpbkZPVl972NFRfv6du4grSWZvS4LFuABhyv2w4ahX79/v/2KPQAK\n7ptv4mCiPtFjB2xZtm7FhdD14eGhzey59hV3acJz58wrM503D99LUzJ7gLa+vaHSy84We5Y33sDv\nu+4qZA8jJPYmopvZCyn2W7bgxJrkZP77nZw6brGwYj9rFs5UBcAf4q1bKCJBQR0XewB+K0epxFNn\n3SXs2KyP7SXT0R/fuHFoT3HtGpEIT93r68330m3JxIk4wG2oT1NiIg7YMgzab3I5jr3ow1Bmz3r2\nhgZ4+Xj6afwcTRV7UzJ7tRrbKlvbr+dDJsMD4J07pp+d2CtGlyUkEK7Y9+pl2G/uTE6cAHjvPYCM\nDMy4+OhoZs+1cV5+ue19bImhtcS+urr97NHz53GwUdcqYGvtb9zAsZOOrhYkFvO3QJ4wAX3mzhyc\n7SgymfHa/CefRI/+m2/w4DVxYtvVpHRhB2jZXvYsPXrgWVhLC9qHERGmx+nkhONH+sYJdPH31zZ6\n01d6WVKCt3HnWHQWYjHA669jia21yjuFwmD4TU1NbrGxseekUmmeRCIpTE1NXQ0AcP78+RGjRo06\nEx4enp+UlLSvvr6+J/uc1atXpw4aNKhoyJAhl48cOfK0OcFoNChi9gjXxhEisy8sxBr/l18G+Ne/\nUAz1Yc3MXh9BQR1b/NpQZq9PUNisz5yyS0uYOdP0PjP2jLMzWjdLlmCnz9mzDT9e3wCtSITf+WvX\nsFrG3DOqadNMt358fLQ95fVl9mVl7avQOpM332xbyfWwYlDs3dzcmuRy+di8vDxpfn5+uFwuH3vq\n1Km4119//Z9r165dkp+fH56cnLxn3bp1iwEACgsLJdu3b3+xsLBQcujQoWcWLFjwlUajMfl4eP26\n/U5W0M3sbT1Au2EDdkwsKUEbxRDW9Oz1ERjY8czeyQkze110J4uxhITgpKdVqzp30W1XV/5Vpx5G\nhg1DG06t5m/oxkWfjQOAZ5FyOX8FmDXhlnnqE/vyctPPFKwVk6GFTx4WjAqxh4fHAwAAlUrlqlar\nnby9ve8WFRUNGj169EkAgPj4+KO7du2aAgCwd+/eSTNnzvzOxcWlJTAwsDQkJORaZmZmjKHtc6mu\nto+eM3yY6tmnpZm+Sr05/PwzDq45m2C8OTt3LLPn2jj6sIaNExDAn9nn5vJn9i+8gJaETAYwZ47l\n+3Y0FizAM0NjFTTcAVpdsffyQrE3x8KxBHsU+66CUenQaDTiyMjInOLi4uD58+d/HRYWVhAWFlaw\nd+/eSZMmTdq7c+fOaeXl5QMAACorK/uPHDnyLPvcgICAihs3bvD0DARYzmnRKJPJQCaT/Sr2DNO5\n2YMl6Fbj8GX2Z85gKZy/v+HM8+efcULI736nf01MLvfv44/V1H4cTk6dn9mzg7WmPJaPpiY8FdcV\ne5UKOw8OG9b+OSIRdrXU7WxJGMeUXvnu7vi91if2CkXH5m+YAp/Yu7mhxdvUhP+XlXX+QcdeUCgU\noFAorLIto2IvFos1eXl50rq6Os+EhITDCoVCtmnTpjlvv/32+hUrVixLSkra5+rqqrfLiUgk4q2k\n5oo9C3tKz06TtieMZfYMg5URI0YYn736z3/iZJjvvsMFLYx5oJmZ/KsB6cMWNo6zMx7Uysqwg6Il\n+wgMbF+rX1CAB0p7+/wdAa6No3tm5+2NCY8QmT07ZlBXh7+B8nKc1OcIsIkwSxpbGmcBJvvpnp6e\ndYmJifuzs7OjBg8efOXw4cMJ2dnZUTNmzNgWHBxcDADg7+9/g83yAQAqKioC/P39TZ56w2Z59mjl\nGMvs9+7F29LSjIu9QoHtBqKisBe9MU6d4l8gQh8dHaA1xcYB6JiVw4q9bmZvbrUHYT2M2Tiurqb3\n3bcUtqafYbRiz+6frbUnG8cyDIp9TU1Nb6VS6QUA0NjY6J6enj4hIiIi9/bt230A0OJZuXLlh/Pn\nz/8aACApKWnftm3bZqhUKteSkpKgoqKiQTExMZmmBsNm9txe3faCscx+61ZsQhYWZljsS0vxSzx0\nKPZFuXzZ+L5//hnL6EzF1Mz+8GH+GaymWjMdqchhxV53gFbf4CzR+RgaoPXyQmvNFNuxI7i4YPZe\nX99W7IODtatIkdhbhkGxr6qq8hs3blyGVCrNi42NPTdx4sQfx48ff2zr1q0vDR48+MrQoUMvBQQE\nVLz66qvfAABIJJLC6dOn75BIJIXPPvvswa+++mqBPhuHD/aH35HMXqk0TUDNQaPBU0sfH7zeowce\nkLjZM9sTPSgIGz+xy7HpwrajFYkAhgzRvzJOfT3Wm7e2YuOpJ54wPV5TMvvmZlypiF1aT/c+U8S+\nIxU5zc040aa+vu2CGzk5JPZCwWb2unX2APjdt9UZF2vlcMVeKsVE4N49/G7rm2NC6MegZz98+PAL\nOTk5kbq3v/POO5+/8847n/M954MPPlj1wQcfWNRk1xqZ/d69eNm92/Jt6KJUYi8YthJGLEbBr6/H\njKe5GStw2O6JISEo4pHt3jmt2ANo17zUaLQTNlpasOpGLkfbSKVCb9ychk+mZPbl5ZjVl5a2b0tg\njo3z44+mx8WluRkzuH798HN/9FF87bm52kW5CdtiKLOfN8/w+gPWhE/sIyJw7gOb1dtbAcfDgF3N\nCauuxmyvI5n9gwfWX3CDa+GwcK2cq1cxy2Wz4aFD9Vs5XLHv1Qsv3EHKw4fxzKCqCg8g+/Zhjb05\nmFJ6ydovfDaMOZn9f/9rXmy6+/Dz0x7kz5/HA4ihWZ5E52FI7B95xHbtAgxl9mThWI7dtEvQaFCk\no6M7ltl3htjX1LQXe+4g7cWLbUsF9Yl9aSmKHHdh5CFD0HZiv8DffIPTs9m2BJY04zKl9JIVeT4b\nxlSxf+Q2RDv8AAAgAElEQVQR7WxHc2H34eurHaQ9cwbX9ySEwdAArS3hE/uQEPxdX7hAYm8pdpPZ\n19aigHp62p/Y19a2t1G4mX1BgWlif+IENm/inoIOHqwdY6itxUWjOzpbz9TM3t+fP7M31cbRXcHI\nHPgyexJ7YTFUemlLvL0xiWCtPgBMYIYPx15GJPaWYTdizy5IzDZjshS2Q11Li/Vi47NxdDN7bodE\nfWLPJ2bcQdpt2wCee874Qs7GMDWzl8k6ZuOwYq+vJ70h2H3074+rGQGQ2AuNPWX2VVV4wOE2H4uI\nwDJkEnvLsBuxr67GU3q2GZOlsM+trbVOXAD8Ng43s9e1cUJD0R7RPeDwiRmb2TMMTrYy1qzKFEwZ\noC0pwV4pHbFx3NzwLKWpyfwY2X3MmoVzDX75Bd/P0FDzt0VYB0OevS3x9sYqMd2JdVKp4ZW2CMPY\nndhbI7MHsI6VU1iIQsZn47CZ/YMHOMAaEqK9z80Nr589q72tvh4zWN3yNTaz37kTs5gJEzoet1iM\nBw/uknS6lJaipVRR0d7yMdXGAWi7sIU5sGIfHIwLdk+ejJ0RH/Y2sg8z9i727G+HxN4y7OanZe3M\n3hpiP3kywJdfGq7GuXQJs1HdBmVz5wL87W/a61lZmJnoiujAgbj9xYux74g1xE4kMlxr39yMZyuP\nPYYHMd1ae3P63Vjq23P3kZqKsZKFIyzsb4+vzt6W6BP7YcMwySKxtwy7E3trZPYuLh0X+7o6XOf0\niy+wFFKfZ5+fz9+0a/ZsXPD55k28rs+PFouxt8yIEcZb0JqDIbEvK8NGZk5OWD6p69urVLYV++7d\ncd3XhQvN3w5hPVxc8PvY0CC82FdVtRd7ti8OW6lGmIddin1HM/uBA7ULU1tKdjZAbCzGdPgwfzXO\nf/8LsGIFf/WMtzcuAbdxI143NPi4bBnA57xT1CzHkG/PXTyZbxZsc7PpNg53yTpT0WjaV3wMHaqd\noUwIh4cHJjFCi311NX8zPJqDYTl2J/asb2gpDx7gbMyOZvZZWTiT8+230bfny+y//RaXeps0iX8b\n8+fjmcGPP6J/r0/sp061/mIZhsovdcVeN7PvbBuH3T7NgrQ/2DbHQoq9jw9+d6nzqXWxG7Hnll52\nNLO3lthHR2PrgthYtD24hITggtXr1unfRmQkwGef4ZJm3t6mL7psDQyVX3LFnq+ZWWfbOJb2wCc6\nH3d3tDCFrrNnYyGsh92IfXU19kmxxgCtNcXe1RWzcl2LIS4Oe9Ib+1HMnInPLyjoWDzmYsjGKSmx\nro1DYt91sBcbB4DE3trYhdhrNPgF8/KyzgBtYGDHxP7mTVwdypqLWts6UzI0QGuK2FNm75iwmb2Q\nYu/igoP2JPbWxS7EvqEBP1gnJ/vI7LOycGGRh9lT1pfZNzbiJLARI/D6o4/iWRW3zbA5No4ldfYk\n9vaLPYg9AH6vSOyti100Qquv146yWyOzf/TR9tU4DQ1YPXP7Nl5cXPQPrLIWzsOMvsz+9GnsMcK+\n3y4uWL10/bp2FSKycRwXe7BxAEjsOwO7EPt797BfPEDHBmjZpcwCAlCA1GoUPYbBxT8aG3FcoE8f\ngIwM/g56DAOwaxdW0TzM6Mvsjx1rv2B3aCi2aeaKPdk4jom7OyZfJPZdD7uwcbiZfUdKL5uaUERc\nXPDLwvbHyclBD/7KFYCTJ3FhkylTsEWBLidO4BjCmDGWxWAv6Cu95BP7QYO0S74BmF+NY26dPYm9\n/eLhgQkPiX3Xwy7E/t69tjaOpZn9gwf4fADM3lnf/j//AXjppbYe/IsvAuzY0X4bX34JsGDBw+3X\nA/CXXiqV2N5Bt96fzexZyMZxXFiBFbL0EoDEvjOwG7FnbZyOZPZ8Yq9WY+vgl19u+9ixY7EdArfG\nvLIS+8m/8opl+7cn+DJ7hQKFXldoQ0O1bYYByMZxZFiBFTqzl0rbNhckOo5BsW9qanKLjY09J5VK\n8yQSSWFqaupqAIDMzMyYmJiYzIiIiNzo6OisrKysaPbxM2fO/C48PDxfIpEUrlmz5n1TgtAdoLVW\nZn/rFq7l6ueH3SW5uLgAJCe3tXL+9S/M+LvClGy+zD4jo72FA9A+s7dkUpU5Pe25i1IQ9gX7+xFa\n7N95p+OL+BBtMSj2bm5uTXK5fGxeXp40Pz8/XC6Xjz116lTc0qVLP1mxYsWy3NzciI8//vgPS5Ys\nWQsAsG3bthkAAPn5+eG//PLL43//+9/fLCsrG2gsiM7I7Pv2BTh0CLtJ6svUX3wRYOtWrVBt397+\nDOBhhW+AtqyMP1vy90fBrq/H6+bYOK6ueOE7QN+9i2MlulBmb7/YS2ZPWB+jNo6Hh8cDAACVSuWq\nVqudvL297/r6+lbX1dV5AgAolUovf3//GwAAfn5+VQ0NDd3VarVTQ0NDd1dXV1WvXr3uGdsHN7N3\nc0MxMNSLXR9csQ8KAjh+HDMEfd0Ux47FmuKsLFxApLYWq3a6Anyll2o1/49YLMaDAGvlmCvGfFaO\nWo3r5/btixO3Cgu195HY2y/2ktkT1sdo6aVGoxFHRkbmFBcXB8+fP//rsLCwgjVr1rwfFxd3atGi\nRZ9qNBrx6dOnnwAASEhIOLxly5YUPz+/qgcPHnj85S9/+Z2Xlxevo7t8+fJf/y8qksGwYTIAQOFx\nc8Psvnt3814MV+wXLcKLoYFWsRiblX31Fc6WnTq16yyewZfZt7biQYAP1rePjDTPxgHQir2/v/a2\nAwdQ6K9dwxW4pk0DyMzEz5StmiLsD8rs7QuFQgEKhcIq2zIq9mKxWJOXlyetq6vzTEhIOKxQKGQr\nV678cP369W8nJyfv2blz57S5c+duTE9Pn/Dtt9/OamxsdK+qqvK7c+eOz+jRo0+OHz/+WFBQULvF\n77hiv3Ch1sYB0Fo5HRF7U6tpXnsNs9revQG++ca8/dkzfJl9a2v7RVZYWN+eYcxfvIIvs//iC4Df\n/AYPnm+8gZO5fvMbHBehzN5+IbG3L2QyGchksl+vp6WlWbwtk/NYT0/PusTExP3Z2dlRmZmZMcnJ\nyXsAAKZOnfp9ZmZmDADA6dOnn0hOTt7j5OSk7tOnz+0nn3zy5+zs7Chj2+baOACWD9Jyxd5UevcG\nSErCg0tXsXAA+DN7tdq42LNCb84Zjm6t/ZUrAHl5mM0D4IF3/XosgWUYEnt7hv39CF16SVgfgz/p\nmpqa3kql0gsAoLGx0T09PX2CVCrNCwkJuXb8+PExAAAZGRnjQkNDrwIADBky5HJGRsY4AICGhobu\nZ8+eHTl06NBLfNvmDsJy6+wBLB+ktUTsAQD+8AeAv/6161g4APyll8ZsHFbszRVi3cz+b3/DZRm5\n2+nVC2N68IDE3p6hzL7rYtDGqaqq8ps9e/ZmjUYj1mg04pSUlC3x8fFHN2zYMG/hwoVfNjc3d3N3\nd2/csGHDPACAN9988+9z587dOHz48AsajUY8Z86cTcOGDbvIt+3KSm1Xyfr6tjaOLTN7ALRxulpN\nL1/ppSk2jjmVOCxcsddosJw1Pb3947y98QyAxN5+IbHvuhgU++HDh1/IycmJ1L09Kioq+9y5c7G6\nt3fr1q3522+/nWXKjrliL3Rm3xXhy+wN2TjsSlyVlR3L7H/5BdcIZfvscPH2BrhzB8Xe09O8fRC2\ngapxui6CGRc3bmj/59bZA9g+s++K6Mvs9dk4IhFm9xcumC/23DbHP/wA8MIL/I/z8aHM3t6hzL7r\nYhdirztAS5l9xzF3gBZAK/aW2DhsH6I9e3BmMh9k49g/lNl3XexC7HVtHMrsO465pZcAKPYXL5ov\nxGPGoEe/ZAlm+PrWAuDaOCT29gll9l0XwcS+shL/qtXta+rNEfuWFpwBC0Biz8XcSVUAlov9oEEA\np04B/PgjwOTJ+quauJk99caxT1ixN5QUEA8ngn2kbGZ//z4O6HEnQZlq45SW4oLeubko9CT2Wswd\noAVA0f7vf9vOhDWVxx7DdQMMQZ69/ePhgVn9w97im2iP4DaO7uAsgOHMvrUV4Kmn8DFDh+IiJJ6e\nuI4qib0WcwdoAVDsASwXYnd3wz3Iycaxf3r2xHknRNdDsMy+shJnU+oOzgIYzuw3bUKb4PZtHEh0\ncQH4/nvM8knsteibVGUos+/ZE9tBd5YQ0wCt/ePkBPDhh0JHQXQGgom9uztmebqDswAo2Gy7XS71\n9QAffQTw009tPf6gIICSEhJ7LnyZvTEbBwB9+86aKk82DkEIh2A2jr8/Wjm6s2cB9Gf2n30G8PTT\n2DqXS2AgZfa6WDJAC4Bi39mZPXW9JAjbI5jY9++PYq8vs+fz7NPTsUulLpTZt8eS0ksA9O07U+zJ\nsycIYRBM7IcOBSgo4M/sPTzaZ/YMg48PC2u/Lcrs22PJpCoAXL1r/vzOiYlsHIIQDsHEPjoaF7Pg\ny+zd3dtn9pWV6CX36dN+W0FBJPa66Mvsjdk4Awd2XqtntocO2TgEYXsEE/uYGP1iz2fjFBbyZ/UA\nKFDl5VizT2KP6PPshZws4+KCB/KaGhJ7grA1gol9SAiu/3rtGv8Ara7YFxQASCT82+rWDRchUatp\n0QUWSyZV2QIfH7TuSOwJwrYIJvZiMVo5x461z+x79sQsnYs+v54lKAizepr5h+iWXjIMir0xG6ez\n8fbGv9QugSBsi6BrM8XEAFRUtM/sPT0x6+diTOwDA8nC4aKb2avVeIAV+mDIij1l9gRhWwQXe4D2\nmb2u2DOMYc8eQJvZE4huZm8PFg4A2jgAZLcRhK0RVOzZVrj6xJ5h8DpbidO7t/5tUWbfFt0BWlMq\ncWyBtzd+lkKfYRCEoyGo2Pv5YSWNl1fb27t1Q8uhqQmvG7NwAHCJwx49OifOhxHd0kt7yey9vcnC\nIQghMCj2TU1NbrGxseekUmmeRCIpTE1NXQ0AkJmZGRMTE5MZERGRGx0dnZWVlfXrchX5+fnho0aN\nOjNs2LCL4eHh+c3NzQZ/2sePA4wY0f52rpVjitg/9RTAd98ZfowjwZfZ24PY+/iQ2BOEEBj8+bu5\nuTXJ5fKxHh4eD1pbW53j4uJOnTp1Km7ZsmUrVqxYsSwhIeHwwYMHn12yZMlauVw+trW11TklJWXL\nt99+O2v48OEX7t696+3i4tJiaB+Bgfy3s2Lv64s19PoexyIWaxcwJ9oP0NqTjUNiTxC2x6iN4+Hh\n8QAAQKVSuarVaidvb++7vr6+1XV1dZ4AAEql0svf3/8GAMCRI0eeDg8Pzx8+fPgFAABvb++7YrFY\nY0lg3Mz+zh2ARx6xZCuOi70O0JLYE4QwGP35azQacWRkZE5xcXHw/Pnzvw4LCytYs2bN+3FxcacW\nLVr0qUajEZ85c2YUAEBRUdEgkUjEPPPMM4du377dZ8aMGdsWL168jm+7y5cv//V/mUwGMpmszf1c\nsa+t1VZxEKZhr5k92TgEYToKhQIUCoV1NsYwjEkXpVLpGRsbe1Yul8vGjx9/dPfu3ckMw8COHTum\nxcfHpzMMA+vWrVsUFBR0vba21ufBgwfuo0aNOn3s2LFxutvC3Rpm8mSG2bkT/3/iCYY5edLoUwgO\nW7YwzEsvaa9fv84wjz4qWDi/UlDAMHFxQkdBEA8n/9NOk3WbezG5GsfT07MuMTFxf3Z2dlRmZmZM\ncnLyHgCAqVOnfp+ZmRkDADBgwIDyp5566oSPj88dd3f3xueee+5ATk5OpCUHId3Mnmwc8+CbVGUP\nNo5EAiCXCx0FQTgeBsW+pqamt1Kp9AIAaGxsdE9PT58glUrzQkJCrh0/fnwMAEBGRsa40NDQqwAA\nTz/99JELFy4Mb2xsdG9tbXU+fvz4mLCwsAJLAiMbp2Poevb2YuMA2MdBhyAcDYM/u6qqKr/Zs2dv\n1mg0Yo1GI05JSdkSHx9/dMOGDfMWLlz4ZXNzczd3d/fGDRs2zAPAAdn33nvvT9HR0VkikYhJTEzc\n/+yzzx60JDBW7DUa7IFOYm8e9lp6SRCEMIgYdpqqLXcqEjHG9vvnP2OP+rQ0nHh1755tYusq7NsH\n8I9/APz4I14/fx7glVfwL0EQDycikQgYhrFo/rmgM2gNwWb2VHZpGfbaLoEgCGGwe7GnwVnL4Cu9\nJBuHIByXh0Lsya83H3udVEUQhDDYvdiTjWMZ9jqpiiAIYbB7safM3jIosycIgovdiz1l9pZBpZcE\nQXCxe7GnAVrL0O1nTzYOQTg2div2bm64mtGNG2TjWIJuZk82DkE4NnYr9gCY3ZeUUGZvCVR6SRAE\nF7sX++vXKbO3BHvujUMQhO2xe7Gvq6PM3hLsteslQRDCYPdiD0BibwmU2RMEwcXuxV4k0oo+YTrk\n2RMEwcXuxd7bmzJSS6BJVQRBcLF7safBWcugrpcEQXCxe7Env94y+CZVUWZPEI6L3Ys9ZfaWQZOq\nCILgYtdi7+sL4O8vdBQPJ9T1kiAILnad6734IsCUKUJH8XBCA7QEQXAxmNk3NTW5xcbGnpNKpXkS\niaQwNTV1NQBAZmZmTExMTGZERERudHR0VlZWVjT3eWVlZQN79Ohx/7PPPvt9R4JzcgLo1q0jW3Bc\nqPSSIAguBn/+bm5uTXK5fKyHh8eD1tZW57i4uFOnTp2KW7Zs2YoVK1YsS0hIOHzw4MFnlyxZslYu\nl49ln/fee+/9KTExcX/nh0/ogyZVEQTBxWiu5+Hh8QAAQKVSuarVaidvb++7vr6+1XV1dZ4AAEql\n0svf3/8G+/gffvjhhccee+x69+7dGzovbMIYYjEAwwBoNPg/2TgE4dgY/flrNBpxZGRkTnFxcfD8\n+fO/DgsLK1izZs37cXFxpxYtWvSpRqMRnz59+gkAgPv37/dYu3btkqNHj8avW7dusaHtLl++/Nf/\nZTIZyGSyDr4UgotIpC2/FIsxs3d1FToqgiDMQaFQgEKhsMq2jIq9WCzW5OXlSevq6jwTEhIOKxQK\n2cqVKz9cv37928nJyXt27tw5be7cuRvT09MnLF++fPm77777Zw8PjwcMw4gMbZcr9kTnwJZfurjg\nX3d3oSMiCMIcdBPhtLQ0i7clYhjG5AevWLFimbu7e+PHH3/8h3v37vUCAGAYRuTl5aWsq6vzfOqp\np06Ul5cPAEB7RywWa1asWLFswYIFX7XZqUjEmLNfwjJ69ACorsa/S5YA9O6NfwmCeDgRiURgLJHW\nh8HMvqamprezs3Orl5eXsrGx0T09PX3CH/7wh49DQkKuHT9+fMyYMWOOZ2RkjAsNDb0KAHDixImn\n2OempaV91LNnz3pdoSdsB3eQlqpxCMKxMfjzr6qq8ps9e/ZmjUYj1mg04pSUlC3x8fFHN2zYMG/h\nwoVfNjc3d3N3d2/csGHDPFsFTJgOt/ySqnEIwrExy8ax2k7JxrEJffsCXLgA0K8fwMKFABIJ/iUI\n4uGkIzaOXbdLIDoGZfYEQbCQ2HdhyLMnCIKFxL4Lw+18SZOqCMKxIbHvwnB72pONQxCODYl9F4Yy\ne4IgWEjsuzC6A7Qk9gThuJDYd2F0B2jJxiEIx4XEvgvDzezJxiEIx4bEvgtDpZcEQbCQ2HdhaFIV\nQRAsJPZdGG5mTzYOQTg2JPZdGG7pJdk4BOHYkNh3YcjGIQiChcS+C0M2DkEQLCT2XRiaVEUQBAuJ\nfReGJlURBMFCYt+FoUlVBEGwkNh3YWhSFUEQLCT2XRjdzJ5sHIJwXIyKfVNTk1tsbOw5qVSaJ5FI\nClNTU1cDAGRmZsbExMRkRkRE5EZHR2dlZWVFAwCkp6dPiIqKyg4PD8+PiorKlsvlYzv7RRD8UGZP\nEMSvMAxj9NLQ0ODBMAy0tLQ4x8bGnj158mScTCaTHzp0KIFhGDhw4MCzMplMzjAM5ObmSquqqnwZ\nhoGLFy+G+fv7V+huD3dLdDavvsowGzfi/wMGMExpqbDxEATRMf6nnSbptu7FpFzPw8PjAQCASqVy\nVavVTt7e3nd9fX2r6+rqPAEAlEqll7+//w0AAKlUmsc+TyKRFDY2Nrq3tLS4uLi4tFj9SEUYhAZo\nCYJgMennr9FoxJGRkTnFxcXB8+fP/zosLKxgzZo178fFxZ1atGjRpxqNRnzmzJlRus/btWvXlMcf\nf/wXPqFfvnz5r//LZDKQyWQdeBkEH2TjEMTDjUKhAIVCYZ2NmXMaoFQqPWNjY8/K5XLZ+PHjj+7e\nvTuZYRjYsWPHtPj4+HTuYy9evBgWHBx87fr160G62wGycWzCwoUM89e/4v8+Pgxz+7aw8RAE0TGg\nAzaOWdU4np6edYmJifuzs7OjMjMzY5KTk/cAAEydOvX7zMzMGPZxFRUVAZMnT969ZcuWlKCgoBLr\nHJYIc6F2CQRBsBgV+5qamt5KpdILAKCxsdE9PT19glQqzQsJCbl2/PjxMQAAGRkZ40JDQ68CoH+f\nmJi4/5NPPlk6atSoM50bPmEIapdAEASL0Z9/VVWV3+zZszdrNBqxRqMRp6SkbImPjz+6YcOGeQsX\nLvyyubm5m7u7e+OGDRvmAQB88cUXvykuLg5OS0v7KC0t7SMALMfs3bt3TWe/GKIt1C6BIAgWEdpA\nNt6pSMQIsV9H44MPALp3B/i//wNwcQF48AD/EgTxcCISiYBhGJElz6UZtF0Y1sZhGMrsCcLRIbHv\nwrA2jkYDIBIBiOnTJgiHhX7+XRg2s6dKHIIgSOy7MGxmT5U4BEGQ2Hdh2Mye/HqCIEjsuzBsZk82\nDkEQJPZdGG5mT2JPEI4NiX0XhpvZk41DEI4NiX0XxsUFQKWizJ4gCBL7Lk2vXgD19TRASxAEiX2X\nxssLQKmkAVqCIEjsuzSs2JONQxAEiX0XxssLoK6ObByCIEjsuzRk4xAEwUJi34Xx9ESxb2khsScI\nR4fEvgvj5obdLu/fJxuHIBwdEvsujpcXQG0tZfYE4eiQ2HdxvLwAampI7AnC0TEo9k1NTW6xsbHn\npFJpnkQiKUxNTV0NAJCZmRkTExOTGRERkRsdHZ2VlZUVzT5n9erVqYMGDSoaMmTI5SNHjjzd2S+A\nMAwr9mTjEISDwzCMwUtDQ4MHwzDQ0tLiHBsbe/bkyZNxMplMfujQoQSGYeDAgQPPymQyOcMwUFBQ\nIBkxYkSeSqVyKSkpCQwODr6mVqvFutvE3RK2ICGBYd57j2HGjxc6EoIgOsr/tNOobvNdjNo4Hh4e\nDwAAVCqVq1qtdvL29r7r6+tbXVdX5wkAoFQqvfz9/W8AAOzdu3fSzJkzv3NxcWkJDAwsDQkJuZaZ\nmRnTiccqwgjk2RMEAQBgVAI0Go04MjIyp7i4OHj+/Plfh4WFFaxZs+b9uLi4U4sWLfpUo9GIz5w5\nMwoAoLKysv/IkSPPss8NCAiouHHjhn9nvgDCMF5eADdukI1DEI6OUbEXi8WavLw8aV1dnWdCQsJh\nhUIhW7ly5Yfr169/Ozk5ec/OnTunzZkzZ1N6evoEvueLRCKG7/bly5f/+r9MJgOZTGbhSyAM4eUF\ncP48QN++QkdCEIS5KBQKUCgUVtmWySf3np6edYmJifuzs7OjMjMzY44ePRoPADB16tTvX3/99X8C\nAPj7+98oLy8fwD6noqIigLV4dOGKPdF5eHriAG3//kJHQhCEuegmwmlpaRZvy6BnX1NT01upVHoB\nADQ2Nrqnp6dPkEqleSEhIdeOHz8+BgAgIyNjXGho6FUAgKSkpH3btm2boVKpXEtKSoKKiooGxcTE\nZFocHdFhWM+ebByCcGwMZvZVVVV+s2fP3qzRaMQajUackpKyJT4+/uiGDRvmLVy48Mvm5uZu7u7u\njRs2bJgHACCRSAqnT5++QyKRFDo7O7d+9dVXC/TZOIRt8PICuHuXBmgJwtERYTWPjXcqEjFC7NcR\nOXgQ4LnnAGbNAtiyRehoCILoCCKRCBiGEVnyXJpB28Xx8sK/lNkThGNDYt/FIbEnCAKAxL7Lw4o9\nDdAShGNDYt/FocyeIAgAEvsuj5sbgIsLiT1BODok9l0ckQize7JxCMKxIbF3ALy8KLMnCEeHxN4B\nILEnCILE3gEgG4cgCBJ7B4Aye4IgSOwdABJ7giBIAhyAxx4DeOQRoaMgCEJIqBEaQRDEQwI1QiMI\ngiAMQmJPEAThAJDYEwRBOAAk9gRBEA4AiT1BEIQDQGJPEAThABgU+6amJrfY2NhzUqk0TyKRFKam\npq4GAHjxxRe3R0RE5EZEROQGBQWVRERE5LKPnzlz5nfh4eH5EomkcM2aNe/b4kVYA4VCIXQI7bDH\nmADsMy6KyTQoJtOx17gsxaDYu7m5Ncnl8rF5eXnS/Pz8cLlcPvbUqVNx27dvfzE3NzciNzc3YsqU\nKbumTJmyCwBg27ZtMwAA8vPzw3/55ZfH//73v79ZVlY20BYvpKPY4wdrjzEB2GdcFJNpUEymY69x\nWYrRGbQeHh4PAABUKpWrWq128vHxucPexzCMaMeOHdPlcvlYAAA/P7+qhoaG7mq12qmhoaG7q6ur\nqlevXvc6L3yCIAjCFIx69hqNRiyVSvP69et3c+zYsXKJRFLI3nfy5MnR/fr1uxkcHFwMAJCQkHC4\nV69e9/z8/KoCAwNLFy9evM7Ly0vZmS+AIAiCMAGGYUy6KJVKz9jY2LNyuVzG3vbWW299/ac//eld\n9vqWLVtmTZ48eVdra6vTrVu3+gwePPjy9evXg3S3BQAMXehCF7rQxfyLqZqtezG5EZqnp2ddYmLi\n/uzs7CiZTKZobW113rNnT3JOTk4k+5jTp08/kZycvMfJyUndp0+f208++eTP2dnZUUFBQSXcbVna\n24EgCIKwDIM2Tk1NTW+lUukFANDY2Oienp4+ga28OXr0aPzQoUMv9e/fv5J9/JAhQy5nZGSMAwBo\naGjofvbs2ZFDhw691JkvgCAIgjCOwcy+qqrKb/bs2Zs1Go1Yo9GIU1JStowfP/4YAMD27dtfnDlz\n5nfcx7/55pt/nzt37sbhw4df0Gg04jlz5mwaNmzYxc58AQRBEIQJWOr/WHI5ePDgM4MHD74cEhJS\ntHqaDYEAAAfJSURBVGbNmqW23Df3UlZWNkAmk8klEklBWFjYxc8///xthmGgtrbWJz4+Pn3QoEFX\nJ0yYcOTu3bteto6ttbXVSSqV5j7//PM/2kNMd+/e9ZoyZcr3Q4YMuTR06NDCs2fPxgod06pVq1Il\nEknBsGHDLsycOXNrU1NTN1vH9Nprr23q27fvzWHDhl1gbzMUw6pVq1JDQkKKBg8efPnw4cNP2zKu\nRYsWrRsyZMil8PDw88nJybuVSqWnLePii4m9fPrpp78XiUSa2tpaH3uIaf369b8dMmTIpbCwsItL\nliz5ROiYzp07FxMdHZ0plUpzo6KisjIzM6MtjanTfgy6l9bWVqfg4OBrJSUlgSqVymXEiBF5hYWF\nQ221f+6lqqrKNzc3V8owDNTX1/cIDQ29UlhYOHTx4sVrP/nkkyUMw8CaNWuWLl26dI2tY/vss8/e\ne+mll/4zceLEfQzDgNAxvfLKK5s3btw4h2EYaGlpcVYqlZ5CxlRSUhIYFBR0vampqRvDMDB9+vTt\n33zzzWxbx3TixInROTk5Edwfpr4YCgoKJCNGjMhTqVQuJSUlgcHBwdfUarXYVnEdOXJkAru/pUuX\nrrF1XHwxMQwmXQkJCYcCAwNLWLEXMqaMjIyx8fHx6SqVyoVhGLh161YfoWMaM2aM4tChQwkMw8CB\nAweelclkcktj6rQfg+7l9OnToxISEg6x11evXv3+6tWr37fV/g1dJk2a9EN6enr84MGDL1dXV/dj\nGDwgDB48+LIt4ygvLw8YP3780YyMjLFsZi9kTEql0jMoKOi67u1CxlRbW+sTGhp65c6dO94tLS3O\nzz///I9HjhyZIERMJSUlgdwfpr4YVq1alco9k01ISDh05syZkbaKi3vZvXt38ssvv/ytrePii2nq\n1Kk7z58/H84VeyFjmjZt2o5jx46N032ckDHNmDHju+3bt09nGAa2bt06syOfnc1649y4ccN/wIAB\n5ez1gICAihs3bvjbav/6KC0tDczNzY2IjY09d/PmzX79+vW7CQDQr1+/mzdv3uxny1jefffdP69b\nt26xWCzWsLcJGVNJSUlQnz59br/22mv/ioyMzHnjjTf+0dDQ0F3ImHx8fO78/ve//2zgwIFl/fv3\nr/Ty8lJOmDAhXejPDkD/Z1VZWdk/ICCggn2ckN/9TZs2zXnuuecOCB3X3r17JwUEBFSEh4fnc28X\nMqaioqJBJ06ceGrkyJFnZTKZIjs7O0romNasWfM++31fvHjxutWrV6daGpPNxF4kEjG22pep3L9/\nv8eUKVN2ff755+/07NmznnufSCRibBnzTz/99Hzfvn1vRURE5DJ6SlNtHVNra6tzTk5O5IIFC77K\nycmJ7N69e4NuvyNbx1RcXBz8l7/85XelpaWBlZWV/e/fv9/j22+/nSVkTHwYi0GI+P74xz/+n6ur\nq+qll17aqu8xtojrwYMHHqtWrfogLS3tI/Y2fd95W8UEgN/3u3fvep89e3bkunXrFk+fPn2H0DHN\nnTt34/r1698uKysb+Oc///ndOXPmbLI0JpuJvb+//43y8vIB7PXy8vIB3COTrWlpaXGZMmXKrpSU\nlC0vvPDCDwCYjVVXV/sCYCVS3759b9kqntOnTz+xb9++pKCgoJKZM2d+l5GRMS4lJWWLkDEFBARU\nBAQEVERHR2cBAEydOvX7nJycSF9f32qhYsrOzo564oknTj/yyCO1zs7OrZMnT9595syZUULGxKLv\ns9L97ldUVAT4+/vfsGVs33zzzasHDhx47j//+c/L7G1CxVVcXBxcWloaOGLEiPNBQUElFRUVAY8/\n/vgvN2/e7CfkexUQEFAxefLk3QAA0dHRWWKxWFNTU9NbyJgyMzNjkpOT9wDg7y8zMzMGwLLPzmZi\nHxUVlV1UVDSotLQ0UKVSuW7fvv3FpKSkfbbaPxeGYURz587dKJFICn/3u9/9hb09KSlp3+bNm2cD\nAGzevHk2exCwBatWrfqgvLx8QElJSdC2bdtmjBs3LmPLli0pQsbk6+tbPWDAgPKrV6+GAuDcirCw\nsIKJEyf+KFRMQ4YMuXz27NmRjY2N7gzDiI4ePRovkUgKhYyJRd9nlZSUtG/btm0zVCqVa0lJSVBR\nUdGgmJiYTFvFdejQoWfWrVu3eO/evZPc3NyauPEKEdfw4cMv3Lx5s19JSUlQSUlJUEBAQEVOTk5k\nv379bgr5Xr3wwgs/sPOErl69GqpSqVx79+5dI2RMISEh144fPz4GACAjI2NcaGjoVQALP7vOGiTi\nuxw4cODZ0NDQK8HBwddWrVqVast9cy8nT56ME4lEmhEjRuRJpdJcqVSae/DgwWdqa2t9xo8ff1TI\n0kuGYUChUIxhq3GEjikvL29EVFRUFrdsT+iYPvnkkyVs6eUrr7yyWaVSudg6phkzZnzn5+dX6eLi\nogoICCjftGnTa4Zi+OMf//hBcHDwtcGDB19mqytsEdfGjRvnhISEFA0cOPC/7Hd9/vz5X9kyLjYm\nV1fXZva94t4fFBR0nVt6KVRMKpXKZdasWVuGDRt2ITIy8hduaxhbxsT9TmVlZUXFxMScGzFiRN7I\nkSPP5OTkRFgak4hh7M5KJwiCIKwMrVRFEAThAJDYEwRBOAAk9gRBEA4AiT1BEIQDQGJPEAThAJDY\nEwRBOAD/DxjGNB3JeuQQAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x32ef950>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Loop over all volumes in the 4D image and store the standard deviation for each volume in a list.  Plot these to see if there are any volumes with particularly unusual standard deviation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get standard deviation for each volume; then plot the values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "SPM uses a measure called \"global signal\" to give a very rough estimate of the average signal value within the brain.  The idea is that you need a threshold to select voxels with signal high enough to be inside the brain, and not in the air around the brain, which has very low (near 0) signal.\n",
      "\n",
      "The algorithm goes like this.\n",
      "\n",
      "* Get a single 3D volume V;\n",
      "* Calculate the mean signal of the voxels in V; call that M;\n",
      "* Make a threshold T where T = M / 8.\n",
      "* Select all voxel values in V that are greater than T; call these values W;\n",
      "* Return the mean of all values in W.\n",
      "\n",
      "See [SPM global image signal](http://imaging.mrc-cbu.cam.ac.uk/imaging/PrinciplesStatistics#Global_image_signal).\n",
      "\n",
      "In the SPM code, the algorithm is implemented in a MATLAB function called `spm_global`.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I used the MATLAB script `get_global_signals.m` (in this directory) to run the `spm_global` MATLAB function on the volumes of `ds107_sub012_t1r2.nii`.  The script saved the SPM global values to a text file `global_signals.txt` in this directory.   The first four lines of the `global_signals.txt` file look like this:\n",
      "\n",
      "```\n",
      "376.53\n",
      "375.75\n",
      "375.26\n",
      "376.01\n",
      "```\n",
      "\n",
      "Read these global values calculated by SPM into a list, and plot the values."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Read global signal values calculated by SPM, and plot\n",
      "def spm_mean(fname):\n",
      "    return [V[V > V.mean()/8].mean() for V in nib.load(fname).get_data().T]\n",
      "\n",
      "spm_globals = spm_mean('ds107_sub012_t1r2.nii')\n",
      "plt.plot(spm_globals)\n",
      "print(spm_globals[:5])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[376.53382753187225, 375.75022022330285, 375.26113600893456, 376.00757342351704, 376.83132797616491]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXlcU2f2/0/CIqssYlm1Im6ExYAg2lIbRUst1bp0xKXU\nqe3MVNvOdHEpzsxPmG/H3U6X6eZ3amu1Vey4tloVhajUBRAQFRVErKKIoIQqBELI8/vjfG9zCUlI\nAiQhnPfrdV+Q5C7nZvnccz/PeZ5HwBgDgiAIwrYRWjoAgiAIovshsScIgugFkNgTBEH0AkjsCYIg\negEk9gRBEL0AEnuCIIheQIdi39TU5BQXF3dGLBYXiUSiktTU1FUAALm5uaNHjx6dGxUVVRgbG5uX\nl5cXCwCQmZk5KSYmJj8yMrI4JiYmPzs7e3x3nwRBEATRAYyxDpeGhgYXxhi0tLTYx8XFnT5x4kS8\nRCLJPnjwYCJjDA4cODBZIpFkM8agsLBQXFVV5ccYgwsXLoQFBgZWGnIMWmihhRZaum+xN+SC4OLi\n0ggAoFAoHFtbW+28vLzq/Pz87tTX13sAAMhkMs/AwMBbAABisbiI204kEpXI5XLnlpYWBwcHh5bu\nuFgRBEEQHWOQ2KtUKmF0dHRBeXl5yMKFCz8LCwu7uHr16nfj4+NzFi9evF6lUglPnTo1VnO7nTt3\nzhw1atRZEnqCIAgLY8xtgEwm84iLizudnZ0tSUhIOLJr167pjDHYsWPH7yZOnJjJX/fChQthISEh\nV69duxasuR8AYLTQQgsttBi/mGrjGL3BP/7xj7+vW7dusbu7+6/ccyqVStC3b9967vHNmzeDhg0b\nduXkyZNjtR4UA7YqVqxYYekQ2mGNMTFmnXFRTIZBMRmONcbVGbHvsBqntrbWRyaTeQIAyOVy58zM\nzElisbhoyJAhV48dO/YkAEBWVtaEYcOGlQKgf5+UlLR/zZo1y8aOHXuqo/0TBEEQ3U+Hnn1VVZX/\n/PnzN6tUKqFKpRKmpKRsmThx4pGNGzf+8bXXXvukubm5j7Ozs3zjxo1/BAD497///Xp5eXlIenr6\nivT09BUAWI7p4+NT290nQxAEQejA1FuCzixghTZOdna2pUNohzXGxJh1xkUxGQbFZDjWGBd0wsYR\n4PbmRSAQMEsclyAIoicjEAiAMSYwZVsaLoEgCKIXQGJPEATRCyCxJwiC6AWQ2BMEQfQCSOwJgiB6\nAST2BEEQJtDaCtCTigpJ7AmCIEzg978H2LrV0lEYDok9QRCECTQ2Ari4WDoKwyGxJwiCMAG5nMSe\nIAjC5qHMniAIohdAYk8QBNELaGwEcHa2dBSGQ2JPEARhApTZEwRB9AJI7AmCIHoBNiX2TU1NTnFx\ncWfEYnGRSCQqSU1NXQUAkJubO3r06NG5UVFRhbGxsXl5eXmx3DarVq1KHTp0aNmIESMuHz58+Knu\nPgGCIAhLIJf3LM++w9lNGhoaXBhj0NLSYh8XF3f6xIkT8RKJJPvgwYOJjDE4cODAZIlEks0Yg4sX\nL4pGjhxZpFAoHCoqKgaFhIRcbW1tFWruE6xwpiqCIAhDUSgYc3Aw/3GhOyccd3FxaQQAUCgUjq2t\nrXZeXl51fn5+d+rr6z0AcILxwMDAWwAAe/fufW7OnDnbHBwcWgYNGnR9yJAhV3Nzc0d347WKIAjC\n7PQ0CwfAgAnHVSqVMDo6uqC8vDxk4cKFn4WFhV1cvXr1u/Hx8TmLFy9er1KphKdOnRoLAHD79u2A\nMWPGnOa2DQoKqrx161agtv2mpaX99r9EIgGJRNLpkyEIgjAH5iq7lEqlIJVKu2RfHYq9UChUFRUV\nievr6z0SExMPSaVSyXvvvfe3jz766M/Tp0/f/f333/9uwYIFmzIzMydp214gEGgdF44v9gRBED0J\nc2X2molwenq6yfsyuBrHw8OjPikpaX9+fn5Mbm7u6OnTp+8GAHj++ef/y1k1gYGBt27evDmA26ay\nsjKIs3gIgiCsCZUK4IsvAE6cMH7bnmjj6BX72tpaH5lM5gkAIJfLnTMzMyeJxeKiIUOGXD127NiT\nAABZWVkThg0bVgoAMHXq1H3bt2+frVAoHCsqKoLLysqGjh49Orf7T4MgCMJw7t4FmDgRYOlSgP/+\n1/jte6LY67Vxqqqq/OfPn79ZpVIJVSqVMCUlZcvEiROPbNy48Y+vvfbaJ83NzX2cnZ3lGzdu/CMA\ngEgkKpk1a9YOkUhUYm9vr/z0008X6bJxCILovTAG8PnnKLrDhgHMmWPe43/7LUC/fgBffgmwZYvx\n2/e0ES8BAATMAlOtCAQCZonjEgRhHRQXAzz9NMDcuSi8VVXmPf6rrwJERAA89hhOQnLunHHbHzgA\n8MknAPv3d0t4OhEIBMAYE5iyLfWgJQiiHYwB1NV13/537waYPRvgvfcA7t83//R+ly8DDB8OMGgQ\nwPXrxh+/pw2CBkBiTxBWjVIJsHkzNiaakytXAEZ3Yw+ZPXsApk0DcHICsLcHaGjovmNp48oVgBEj\nADw98bFMZtz2PdGzJ7EnCCtFqQSYNw9thtu3zXvsX38FuHoVs+6u5vp1gFu3AB5/HB/36wdw717X\nH0cX9fUADx4ABAYCCATq7N4YSOwJgugyFi5E0R0xAqCmxrzHVijwb2Fh1+97zx6AKVMA7OzwsbnF\n/soVbBQW/J/zTWJPEITFUCgAtm8H+O47zEDNLfbNzfi3oMDwberqAN59t+P19u5FC4ejX7/uuYPQ\nBWfhcJDYEwRhMfLyAIYOBfDyAujfH6C21rzH5zJ7Y8T+yhWANWsASkr0r3fxIkBsrPqxt7d5M3uu\ncZbDFLHviaWXJPYEYYUcOwbw5JP4f//+lsnsBw82Tuy5C9KXX+pep6UFG0P791c/Zwkbp7NiT9U4\nBEF0CdYg9iNHAlRWYruBIdTUAMTHYycl7s5Ak7t3AXx81H49gGXEnmwcgiAsTksLwKlTAOPG4WNL\niL1CAeDqChAZCVBUZNg2NTUAY8YAhIcD7NunfZ07dwD8/No+5+2t27M/eLBrO1y1tmKV0dCh6udM\nqbUnsScIotMUFAAEB6MIAlgus3d0BIiONtzKqanBrH3+fGxc1kZVFYC/f9vn9GX2H34IcPy44XHz\n2b+//R3GL7/g++nqqn7OlFp7EnvCZtizx/xdwW2BW7c6X64olaotHADLZfZ9+hgn9rW1GOuwYfg+\naENbZq9P7JubsS7eGBgDWL4c4NlnAU6fbvva1asAQ4a0fU4gAHj0UeOsHBJ7wmb4+Wf0jQnj+P57\nrEjpDGfOqDscAVg2sw8IQJ/dEGpqMFYfH93VQ8Zm9qaI/cqVAIcOAUyaBHDjRtvXbt4EGDiw/TbG\n+vYk9r2Aq1cBnnjC0lF0P83N5q19thVqarBRszNcvYrZMYePj2XEvk8fgL59jWug5cReV7zGevam\niP2hQwDr1wNERaG486msBBgwoP02AQHGtQ1Q6WUv4No141vueyJNTd07EJat0lmxZwy/YyEh6ue8\nvVHwlMrOx2conI3j7o5DCxgC59l7eOBYNy0t7depqjLexjF23BpO0AcO1J7ZaxN7Pz+8EBkKlV72\nAu7cMf7L1xNpaqLM3hRqa9GvNnXgsjt3MGPs21f9nJ0ddq4yZ3kiZ+MYI/acZy8U6u4odedOexvH\nywsvZq2t2uMwJrNXqfD9Dww0Tux9fQGqqw0/Dtk4vYDqaoCHD82bZVkCEnvTqKnB74aptkt5edus\nnsPcvWiNzeybm3HhLlK64tWW2dvZ4XG0JVHGin1tLe7L2Vm32AcFtd+u14t9U1OTU1xc3BmxWFwk\nEolKUlNTVwEAJCcnZ0RFRRVGRUUVBgcHV0RFRRVy68+ZM2dbZGRksUgkKlm9erUBI2X0LLhbPWN9\nxJ4Gib1p1NSg0Jhq5egTe3P69sZm9pyFww0upq2RljHtnj2A7vFxjBV7vic/YEBbz56xrrVxeprY\n652W0MnJqSk7O3u8i4tLo1KptI+Pj8/JycmJz8jISObWWbx48XpPT08ZAMD27dtnAwAUFxdHyuVy\nZ5FIVDJ37tzvBg4ceEPXMXoa3NW/vh6/oLZKczN59qZQW4sdkSorAUaNMn57baWBAOYXey6z79MH\nH3MNtrrgxJ5Dm9jX1wM4OLStcefgfHt+ZyfuuMaKPZe5e3nhXVZ9PbYjcPvx8Gi/nSmZvc159i4u\nLo0AAAqFwrG1tdXO29v7t+svY0ywY8eOWXPmzNkGAODv71/V0NDg2traatfQ0ODq6Oio6Nu3r4Ft\n+T0D7upv6759UxM2snGjHxId09qK3wtumAFTsKbMnhN3Q7J7zq/n0Cb22vx6Dl2NtMY20PLFXiBA\nK4fL7rmsX6BlUj9O7A3pRdvSgu0SDg6Gx2UN6M3sAQBUKpUwOjq6oLy8PGThwoWfiUSi38a0O3Hi\nxBO+vr7VISEh5QAAiYmJh7Zs2ZLi7+9f1djY6PLBBx+8yWX9mqSlpf32v0QiAYlE0umTMQfV1Zgx\n9AaxB8DsXtttN9Gee/ewN+ajj3a92Ju7/JKzcQDUYs/P3DXhyi45tIm9Nr+eQ1eDbmcyewC1bx8e\nrtvCAQBwc0MBf/CgbeO4NsyZ1UulUpBKpV2yrw7FXigUqoqKisT19fUeiYmJh6RSqUQikUgBALZt\n2zZn7ty533Hrbt269QW5XO5cVVXlf//+fe8nnnjiREJCwtHg4OAKzf3yxb4nUV2Ngyj1Bs8eAH1U\nEnvD4LLboCCAzEzT9nH1qu7MvrS0c/EZA2fjAKDYd1Rrr03sKzR+9R1l9pqevVKJ1TUPHmDGrS0j\n16SyEiAhQf2Y30irq3GWg8vuDRF7c/n1molwenq6yfsyuBrHw8OjPikpaX9+fn4MAIBSqbTfvXv3\n9OTk5AxunZMnTz42ffr03XZ2dq39+/evefzxx3/m1rcFlErMdENCekdmr6+zC9EezrcOCjIts5fJ\nUGQfeaT9a5ZqoAVA8evIxjHEs9fVOAug3cZpbsY5ap2csALOELRl9pyNoy+zBzDct++JjbMAHYh9\nbW2tj0wm8wQAkMvlzpmZmZO4ypsjR45MDA0NvRQQEPDb7JgjRoy4nJWVNQEAoKGhwfX06dNjQkND\nL3XnCZiTmhoUwH79bF/sm5sxC6NGWsPhsltTxZ6zcLRlsJZqoAXoOs/eWBuHi8HDw/Dfm6bYDxjQ\nNrPXJ/aGVuTYpNhXVVX5T5gwIUssFhfFxcWdmTJlyg8JCQlHAQAyMjKSuYZZjj/96U9fKBQKx4iI\niPOjR4/OXbBgwabw8PAL3XkC5oTLTDw9e4eNExBAmb0xcIIXGIiiY8yQuQC6LRwAy2b2hoi9IZ69\nsQ20XCOxob83xvB9DwxUP6dp4/TmzF6vZx8REXG+oKAgWttrX3311Uuaz/Xp06d569atL3RVcNZG\ndTV+ITw8cKhUW4bE3ng4wXN1Revh/n3jynP1if0jj6BYGupddxZjq3EMEfubN/E7pQ1f3/ZZNRcD\nv2xSH3V1uL6bm/o5vtjrGheHw8/PtsWeetAaAT+zt3Ubp6kJszAS+/ZERwNc0HK/yvetg4J0D/Or\ni9xc3Lc2/PzQOy8uNm6fpmKsjWOIZ3/5ctsZoviEhQGcP9/2bshYsde0cADw8e3bAGVlhjXQGmLj\nyOU9r8YegMTeKLjM3tZtHMbIs9dFfT2OV3/gQPvX+Nmtsb69SgWQk6N7RFWBAGDGDIBdu4yP2RQ6\na+O4umK/g8ZGfCyTYSOrLrH19cXj8Xu8doXY9+kD8Le/4ZDRzs54LrqwdRuHxN4I+DaOLWf2CgX+\n8HR1Ye/NFBZiZxptpZX8RsqgoPbjsnD88gtO38cf+OvSJUwi+H6zJjNmAOzcaXrsxmBM6eXDh/g6\nN7MWAF6cfHzUPvylS5jV67OgRo4EOHdO/Zjv2Rvye9Mm9gAAf/87vpaXp3/7jmyc5ma8SyCx7wX0\nFhunqQk9Zyq9bM/ZswBz5uAMSHJ529f4VsbTTwNs2KD9e7JjB05Qwhef48c7nichLg4/jytXOncO\nhmBMZv+f/wBMmwZgr9ECyO8IdukSQGio/mPqEntDM/tr13TfOTg6AgwerH97bTYO31Zavx4/g1u3\nSOxtnt5i45DY6+bsWYDx43H8m59/bvsa38qYORNg8mSAefPaD3f8/fcAMTFtp308caJjsRcKAaZP\nB9i9u/1rXI1+V8FvoNVXZ69Q4EVt2bL2r/FHvuxusVcqAbZtw6kITUVzyIS8PByrh+vU9e23OG7R\nP/9JYm/zcJm9tdo4164BvPde5/fT1IQ/Mm9vy3r2KhXAli3GlTBWVeFk3Z9+avqY8vo4exYHOJs4\nsa2Vw1j7WvMNG9DeWLtW/dz169izdM0atdgzhpn9uHEdH1+X2CcnA2zcaNIpacXQBtpt29Ce0Tbo\nG7+RtrvFft8+tMBMGXyOw9UV704ePMBjv/QSfjb/+7/YMC6X42c2aBCJvc3D9+x//dX4Ouru4N13\n8ct48SJmnP/4h+GTTeiC67lo6cz+5k2AF18EyMjoeF2ODRsARo/Gi8TMmV0bz6+/ovcbGorzm/LF\n/sEDtAqcnNTPOThgNvj++wBFRfjc99+jYI8bh8J/+zb+VSq1j3apyRNPoHDyP5d79wCOHgXIyjLs\nPFpbcUL58ePx7mPTpvazShli4zQ14Xyv2rJ6AOPFfvhwfH8bGtQxGNqp6uOPAf78Z/3rGIKvL35G\nS5fi57FjB8AHHwBs3oz2nYsLfu6vv975Y5kbEnsDkcvVwxrb22PLvrYu3IZcAB48APjhB5wr09D5\nPbVRVwfwySeY6UVGYtVBbCxAfr7p+wRQ2zicXdUdGbIhlJXhj+/dd9v743xUKnzf791D4dqwAW2R\nc+fQW+8qiooAIiLw84+LQ2G6eBFf0yw95Bg4EOOZNw8vyl99BfC73+E+nnoKxWTFCoAJEwyrn+/T\nBy8UR46on9u7F+M5dsywz2rNGmy0fPVVgN//HuDzzwFWrWq7jiGZ/d/+hu8HfywaPpzYy+V4UdPV\nh4DDwQHvEs6fx8eGdqrKy8Nxg7ri4r5wITaCl5Tg3eGoUTgf8IcfAsydi+v4+/fM8aKsSux//RXg\nlVcsHYV2du3CH6SdHT7W1kgrl6MX21GGtXw5/thefx2/UIaQl4cZBx+pFOCxx7AM8MEDgD/8AbPa\n3FzD9qkLTuzt7LCDSn09Zo5cGZ25KCsDmDoVa88/+ED7OpcuYUeZxx8HeOMN/MEHBaGYvvEG/ki7\nCs7CAUBhSk3FBQDtI76Fw+eFFwBSUvBzmTEDM2oAgKQkgLfewvfZ0O8BAF4kDh9WP/7+e8xq+/XT\nXv+vyY8/4vuZnIzLf/8L8NFHaAMC4AWjtVXd4KpN7KVStHA+/1z3RWrkSExELlzAxlHNBlxd23BW\nTkc2jlKJF6nJkwHWreuaIYffeQd/T5mZ6g5gy5cDxMfjyJk9GsaY2Rc8bHuyshhzddX6ksWRSBj7\n/nv147AwxoqL267z1luMeXoytmiR7v20tDDm68tYWRlj//kPY/PnG3b8uXMZ8/Bg7P599XMLFzK2\nbl3b9bZtY2z6dMP2qYtjxxiLj8f/g4MZO3yYMXt7xr7+unP7NZY338TzKylhzM+PMZWq7evnzzPm\n749xbdnCWFwcY1evql+XyRjz8mKsslL9XFERY7/8YlwcVVWMffYZYyNHMvbll+rnm5oYe/RRxrZu\nZWzoUMbef9+4/TY3M5aba9w2jDF26RJjAwbg+3H/PmPu7oz9+itjf/gDYx9+qH9bmYwxNzfG5PK2\nz69axVhSEu5TLmesTx/1a5WV+D5zyOWMDR7M2A8/6D+WSsXY1KmMRUYy9vzzhp3bhg2MvfEG/r9x\nI2Mvv8zYuXOMhYe33/fLLzM2bhxj164Ztu/OoPndsxT/p50m6a5VZfbFxejXaZt42JKUleFt3dSp\n6uc0by2lUvSW9+zBzECXnSOVYiY6ZAjesl6+3PHxW1oAfvoJs4sNG9TPHzmCDYV8ujKzB0Df/s03\nMavZs6dz+zWWsjKshggNRduspKTt6ytXYiY2fz5mz6dPt7UKPDzw+U8+wceMoW3BZeOGwBjA889j\npvenP6lv5QEw63zvPTzGwoWYpRuDoyPabsYyfDhm05cvA3zzDdoo7u4AEgl+vxjD74u235FUijX+\n/LYFAIC338b37/bt9rNSadbZr12LGXhHlS8CAcAXX6jbOQyBm3wcQH9mv2YNQEEBNpgGBxu2785g\njiEquh1TrxKdWUBHZv/SS4wBMFZX1wWXwC5k2TLGFi9u+1xSUtvMZvJkxr75BjOAAQMw+9LGggWM\nrV+P/9fU4J1AR1nD0aOMxcYydv06Y97ejFVX4//9+zPW2tp2XZWKMR+fttmssezbx9izz+L/kybh\nncgvvzDWty9jDQ2m79dYhg5l7OJF/P+VV9pnrSNHMpafr38f5eWM9euH7/XJk4wNHIh3SPfuqddp\nbsbzralpv/3Bg4yNGMGYUql9/62teLdgbl55hbGoKDyfwkJ87tYt/H689BJjdnaMffJJ++1ef52x\nNWu073PECHy/q6vxO8TR2sqYQIB/KyrwGNevGx5rbi5uZwgZGYzNnIn/r1+Pd8t1dXj3wnHjBsZw\n65bhMdgKYCuZPefVWUsNe3MzwL/+hSVtmm0J/AoBuRy7uk+ZghnAM89o707f3IzZcfL/zeDr44N+\n7d277ddtaQHIzsb/9+7Fu4pHHwVYsAD3//nnmNEJNT5BgaDz2T0/sx88GOB//gcbGkeNatsw2J0o\nldgDlcvUExKw3YCjtRUb5XSNtcIxeDDArFmYCX72Gfr4zzyDVTIchYXoY2vOC8EYtq2kp6vbajQR\nCjHLNTfz5uF4MoWFAGIxPhcQgO0Gd+9iH4C0tPYjZWq7G+Rwd8eiA37jLACeo4sLvrZ2LcBrr+F3\n0VBiY7Fc0RDc3NpX4/Tt2/aO//JlPGddg6oR2rEasVcqsbHt0Uc7V6HSlXC37zk5eOvMh2/jHDuG\nP3hPT3w8ebJ2sT9yBEAkatvLT5eV85//YINwejrWED/3HD6/di0K/gcfYEOdNuLisIemqfDF/rPP\nsOEXAHtJmsvKuX4dqx440ZkwAd9npRIfV1RgpY62yas1+fvfsUpn3z6snX7lFayM4ay2nByA2bMB\ntm9vaxX99BMKzvPPd+mpdQkSCZaX8ocoAECbZs8e/A7Mm9e2LLKyEsWfuzho4uamrjHnyi45uI5V\nZ87gxbK7cHVtL/ZCoTo2AN2TshP6sRqxLy3FThEBAdYh9ioVdnTZuhUFWhN+Nc5PP6HAcyQk4I9C\n8w7l4MH2Pufw4e3FvqEBs+kDB9RjoXCVAAIBwKJFWDmRkqI99uhozPhMhetUxR2P47nnsGT09m3t\n23UlpaXo13M88gjeXXBlpYbUbXP4+2PlU3IyVqxIJFhZxA1XkJODF7LlywGWLFFvt2MHwB//2P7u\nyZrx81NXvaSlYQnq6tV49/nWW3iHqOt8dGX23Gv37+P7HhnZffG7uqpLmvltB3zfnsTeNCz2Ndac\nqKC4GLPjvn2tQ+yvXUNB18ycOPg2jqbYu7lhCSC/5yQA1tUnJrZ9bsSI9mOdfPQRdp6ZPBkztZ07\n2zcQ+fvrLmVzd9dfl94R/Myez6OP4oUmPBzg//0/0/dvCGVlWN/MZ+JEtZVTUqL9IqyLtDS8SwFA\nsXv5Zbx7YgzF/vHH8dxOncIMWKXCz7U7s9juxsMDvz9ff43vlb29/hJPfmavTezPnMHG0O7sPaot\nswfAuzhuyGgSe9PQK/ZNTU1OcXFxZ8RicZFIJCpJTU1dBQCQnJycERUVVRgVFVUYHBxcwU1VCABQ\nXFwcOXbs2FPh4eEXIiMji5ubm/to27dm5nnunFrsrcGz53uh2uBsnPJy/IForvvee+irc0O2VlTg\n+ppZkaaNU1eHPS7/53/wsbe37jHOdWFn17mKJq4HrTbS07Fu+quv2nZt72o0M3sAtK24IQaMFXuB\noG1GO38+1qefPYsiFxSEwjJtGj5fWIjvvTkqPbqTwEAU/L/9DdspdH2uAOrMXpuN4+6ONlpUVLeG\nq1Pso6LUmqFvkhdCN3q7OTg5OTVlZ2ePd3FxaVQqlfbx8fE5OTk58RkZGcncOosXL17v6ekpA8BJ\nyFNSUrZs3br1hYiIiPN1dXVeDg4OLdr2XVjYtqHo3Dksbbtxwzoy+6Ii/V9sLy+0NH7+GTvHaGbe\nQUGYKS5fjt7qoUMoVpq30Jpiv2YNdrzRzGqNobNiryuz5wgIwM/qk0+6djwWPmVl+L7ySUhA66qi\nAsX+1VdN339AAPZE/ctfsKSVY/Zs9PgfPmx7t9aT8fPDO5mO4DJ7XTbOsWNYhtud6BL76Gi03VQq\nvOsmsTeeDm0cFxeXRgAAhULh2Nraauft7f3bqByMMcGOHTtmcXPRHj58+KnIyMjiiIiI8wAAXl5e\ndUKhUGsHbs3M3tpsnMJC/WL/zDNoA3zxhe5emsuW4cXgP//BHo+aFg4AZo63b6PtcusWNhx21iKx\nt1c3ZOrj00/bj4kC0LHYA2Cj7fffd8+AcAoFdpnXbBR3cMChBr77Di+Qhnr2unjlFYCTJ9uK/YQJ\neDHZtKlnWzim4OamP7O/ebP7M3tt1TgAWAl29iz+Rry8DGuYJ9rSYQdmlUoljI6OLigvLw9ZuHDh\nZyKR6Ld6hRMnTjzh6+tbHRISUg4AUFZWNlQgELCnn376YE1NTf/Zs2dvX7JkyTpt+83MTIO0NPx/\n5EgJNDRI4NFHrcvG0ffFdnVtn3lq4uaGGf24cXhOn3/efh17eywP/PJLLLV8+WX9E1gYgiGZ/YMH\nWELn7Y3ZLJ+mJvR79eHri2L49demZXs//4xd6devb//a++9jJqdt/PG5c7E6pm9fdfWTqUyejLba\nhAnq5+ztcf/ffNPxkMO2hrs7iqmuzB6g+8Xe0RGzd4WirdhHRKC1d/Fi7/LrpVIpSKXSLtlXh2Iv\nFApVRUVF4vr6eo/ExMRDUqlUIpFIpAAA27ZtmzN37tzvuHVbWloccnJy4vPz82OcnZ3lCQkJR0eN\nGnV2woQtzZJCAAAgAElEQVQJ7UaLaWhIgyVLUDQPHMArt0CAIlNRYfoJtbbil6QzjUh37uA+9E1O\nbChDh2IN9+bNWFGijfnzscTT19e43p26METsr17FTHnDBqxS4dtQ+jx7PgsWAPz1r6aJ/dWr2BC9\nbFnbMWWuX8cLQG6u9l6Ljz2GAqDp55uCvT3adZrHefVVjElT8GwdfmavTeyDgzt/ge0IgUBt5fDj\ncHLCz3zPnt4l9hKJBCQSyW+P0zU7gxiBwdU4Hh4e9UlJSfvz8/NjANCf37179/Tk5OTfBqAdMGDA\nzXHjxh339va+7+zsLH/mmWcOFBQUaG1eDA1VT5589iwOIAbQeRtn506spe4MnF/fVV2kR41CYdPF\nsmXYcWrjRrxF7SyGiH1ZGWbmMhlWo/AxxMYBwB+fMfOs8mlsRAtp8+a2z6emoo+ua1YhoRB9+67q\nyKTtM46MbN/BqjfADXimq86+u7N6Dm1iD4C/o507e5fYdyV6xb62ttZHJpN5AgDI5XLnzMzMSVzl\nzZEjRyaGhoZeCggI+K3qOjEx8dD58+cj5HK5s1KptD927NiTYWFhF7Xt+7HHsEoAoO1ogp21ce7f\nV48dbiodWTjWjqFiP3w41l7zx9sBMFzsfX2xt6YpQyDL5djxZ+NGdecmxrBto6ORT9PS2g/JS3Qe\nLrPXZuM8/TSO/2MOdIl9dDQOmUyNs6ahV+yrqqr8J0yYkCUWi4vi4uLOTJky5YeEhISjAAAZGRnJ\nXMMsh6enp+ztt99+PzY2Ni8qKqpw1KhRZydPnvyTtn1PnYo9GgGwowwn9tzEIKbS0IAWQVOT6fuQ\nSlGIeiqGNNBydewvvIACy784GCr2ffqoO9sYi1yOXnmfPuqLfkUF2m/+/vq3tbdvn3kSnUdfZh8b\nq3uYha6Ga6TVJvYAlNmbil7PPiIi4rwuG+arr77SapbMmzfv23nz5n2r7TU+Tz6JFRXnzuEtPVfP\n3Fkbp6EBM83SUu09/ZRK7ByiqwfqzZtY4mXuER67EkMy+9JS9Nz79sUM/do1tQ/O70HbEX5+2Mah\nbeIOfTQ2orC/9BL2Uh4/Hn16U0aBJLoGfZm9OdGV2Y8ciUkIib1pWKwHraMjliKmp6sbZwE6b+Nw\nE2xc1GIeVVbiLeBrr2H3+evX26/zzTc4cJazs+kxWBpDbRxO3EWitu+XoQ20ACj2VVXGxyiX43s8\ndSo20KtUKPajRxu/L6Jr4Gf2lhZ7bQ3Frq5499dRpRihHYuO+jF1KpbfcY2zAF1j43h5tRd7uRzn\n/ly0CD39mTOxUZQPY1hK2NkGXkujS+xra/HORibD7J2bWi0srO0AYIbaOADqzN5YuMx+yBC8wBcW\nkthbGn119uZEV2YP0DOnA7QWLCr2kyej/8qfEb4rbJzY2PZi/5e/oLAsXYqPp0/HCw2fnBwsR+zp\ngqNL7BctwtEyy8rwveDupsLC2r5fxoi9v79pYs9l9gDYX2HvXrwI878LhHnhMntrtXGIzmFRsffy\nwtmGuDk5AfCDbmoyrAeoNhoaUKz54vXLL1iy9cUXaoGbOBGzSf543//6Fw4D0NNnpbG31y72cjkO\nBnblSts6dUtm9gAo9p9+iv0a6BbdcnAi29RkvZk9YToWH7x1yRIcdpZDIGg/DZoxNDbioGQ3bqgr\ncj75BKej69tXvZ6zM8CkSdjhCQDLP8+cwSFtezp2dtovltyEIJ980nbsHW7kTe4CYUoDrbHwM/sn\nnsBssqffUfV07O3xc5fJLCuyXDWOpe8wbA2Li702OuPbNzRgL7+QEBSwhw9xnJPXX2+/7vTpOG7N\n3bs4Hk1qas9umOXQZeO0tOBQAKdPt83s3dzUFTkA5m2gBcAsMikJhxkmLIubGw4/bmkbR1dPXsJ0\nOhwuwRJ0xrdvaEB7ICwMx6W5dw/HptE2VO3MmTh++YgReDexa1fn4rYWdIm9UonlpocOtR9kjLNy\nhg41v40DgA3jusbnJ8wHJ/aWtnHu3iWx72qsMrPvTPllYyN+WZ58Egceu3JFPTa8Jk5OAP/+Nw7b\n8NNPtvPF0if2ffsCFBS07zTGb6Q1h9jzM3sAfO91zfNKmA93d+vI7Ovr8fvQk2YJs3asMpfqrI3j\n6oq19K+9Ztg2/DlhbQFdDbQtLfiatsmfRSLsSQtgnNh7e5t2y62Z2RPWgZsbtutYOrO/f992ki9r\nwSqvm11h4/Rm9DXQ6rJKRCKcX5QbXtbQH7tQiH5/dbVxMWpm9oR1YA2ZvTW0G9giViv2nbVxejPc\nra/mAGVKJfYj0MawYVh/z2XoxpSfGtJIq1K1tXvkcrooWyNubvjZWNrGocy+67EpG4cxsgc4ON+e\n73lyNo42PDzwfbt+3XALh0Obb9/aiiNaVlSgzfPjj/iZcjNbNTZSZm+NcJOUkI1je1htZq9P7LOy\n2naG4mhqQjGjqg7tjbT6bBwAzO6Lizsv9r/8AiCR4LSF/fphhc/+/dg7U6nEi3JTE4m9NeLmhn8p\ns7c9rFIW+/bFYYq1UViIE3I7OeG45mVlWB++ezdZOHy0NdLqs3EA1GJv7I9Mc8iEP/wBxT49ve2d\nBdcd38kJM0eqtLA+rCWzt7SVZItY5c9Nl43T0oKDlH30EY6lsm8fZoc//ojZIleJQ2hvpNVn4wCg\n2J8/b1pmz3n29+5hT+TU1PZizn2u1DhrvVhLZm/pGGwRvWLf1NTkFBcXd0YsFheJRKKS1NTUVQAA\nycnJGVFRUYVRUVGFwcHBFdzsVRw3btwY6Obm9nDDhg3vmBKULhvn/fcBAgKwY1BcHA6N+957+KV4\n8IAqcfiY08aJjQU4eBAvJnv34jAU2j4H7nOldhXrhcvsLV2NY+kYbBG9No6Tk1NTdnb2eBcXl0al\nUmkfHx+fk5OTE5+RkZHMrbN48eL1np6eMv52b7/99vtJSUn7TQ1KVzXOrl04GbVmpYi3N0BdHdk4\nfHSJfUc2zi+/GD9GTWwszhm7bRsOOPfCC9rX48TeyYkye2uFE1pL2zgAJPZdTYeevYuLSyMAgEKh\ncGxtbbXz9vb+bRI6xphgx44ds7Kzs38bt3LPnj3TBg8efM3V1bXB1KD69UM7gI9CAXDhgva5Yb28\nUOzJxlGjTew7snFCQvBCamxmD4C2zWuvoXe/bZv2dTixd3enzN5asYbM3tFRPSgb0XV06NmrVCqh\nWCwu8vX1rR4/fny2SCT6bTDcEydOPOHr61sdEhJSDgDw8OFDt7Vr1y5NS0tL60xQ2kr5Ll7E8W24\nzIOPtze23pPYq9HVQKtP7J2cAB591LQf2cSJ6MlLJG1HF+XDt3Eos7dOrCGzB8DfsaVjsDU6zOyF\nQqGqqKhIXF9f75GYmHhIKpVKJBKJFABg27Ztc+bOnfsdt25aWlraW2+99S8XF5dGxpjebjn864FE\nIgGJRPLb4/79Ubz54sSflFwTLrNnjDJGDm0NtB3ZOABo5Zgi9gIBwFdf4WegC07sqYHWerEWv9zV\n1fIxWANSqRSkUmmX7Mvg0ksPD4/6pKSk/fn5+TESiUSqVCrtd+/ePZ0/IXlubu7onTt3zly6dOla\nmUzmKRQKVc7OzvJFixZ9qrk/fcm/vT1m6zU1WNYHgGLPn76QD5fZ9+lDmT2HKTYOAIq9sUMfcISH\n63+dGmitH2sovQTAiw6JfftEOD093eR96f3p19bW+tjb2ys9PT1lcrncOTMzc9KKFSvSAQCOHDky\nMTQ09FJAQMBtbv3jx4+P4wW1wt3d/YE2oTcEzsrhxP7sWZyARBtcZu/uTmLPoSn2KhUuHY0sOWxY\n56aF1Adl9tYPZfa2i16xr6qq8p8/f/5mlUolVKlUwpSUlC0JCQlHAQAyMjKS58yZo6MprvP4+akz\nzOZmHGt95Ejt63KZvVBIGSOHpti3tmJW39GYN/PnAzz3XPfE1LcvjqhI4+JYL+7u+DuydC90Evuu\nR+9HGhERcZ5v0/D56quvXtK3LXcHYCr8Rtrz57HLvS6B8PLCMV1cXCiz59BsoDXEwgFAQdbVwNpZ\nuJJaaqC1Xjw8AKK1/uLNC4l912OVPWgBcNhcTuz1+fUAVI2jDc0G2o4qccwB38ahzN466dMHIC/P\n0lGQ2HcHViv2/My+uFi3hQPQts6eRATRtHEMqcTpbqj0kjAUEvuup0eI/eXLAKGhutf18sLMnnrQ\nqtEUe0NtnO6EGmgJQ6FqnK7HKke9BGgr9pcu4aTguuCGSyAbR422zN5axL6xUV1lRRDamD4dwNPT\n0lHYFlYv9jIZTn6hb55YsnHao9lAa002DmX2REdMmmTpCGwPq7ZxqqsBrlwBGD5cf8mghwdeEB48\noMyeQ7OB1tpsHLooE4R5sVqx9/LC2/2iIv0WDgDWBfftC3DrFok9hzXaOO7uePf18CFl9gRhbqxW\n7AUCLL88dqxjsQdA376yksSewxqrcYRC/HxqaiizJwhzY7ViD4BWjlRqmNh7eWFPWxIRxBqrcQDw\nDuzOHcrsCcLcWLXY+/ridHeGZvYAlNlzaGugtRaxr64msScIc2PVYu/nh7f+Q4d2vK6XF/4lsUe0\n9aC1tI0DgGIvk9EdGEGYG6sX+8GDDetcwWX2JCKINds4AJTZE4S5sXqxN8TCAcDMvk+fjofw7S1Y\nYzUOgFrs6aJMEObFCn7+unn+eYD4eMPW9fYmC4ePNVbjAFBmTxCWwqrF3tcXF0Pw8qJskY+pQxx3\nN5TZE4RlsGobxxi8vCiz52PtNg5l9gRhXjoU+6amJqe4uLgzYrG4SCQSlaSmpq4CAEhOTs6Iiooq\njIqKKgwODq6IiooqBADIzMycFBMTkx8ZGVkcExOTn52dPb67TwKAbBxNrLkax9GR2lYIwtx0mOs5\nOTk1ZWdnj3dxcWlUKpX28fHxOTk5OfEZGRnJ3DqLFy9e7+npKQMA6N+/f82PP/74rJ+f352LFy+G\nJSYmHqqsrNQzjFnXMHYswBdfdPdReg7WXI1DWT1BmB+Dfv4uLi6NAAAKhcKxtbXVztvb+z73GmNM\nsGPHjllcBi8Wi4u410QiUYlcLnduaWlxcHBwaOnq4Pk4Ouqfzaq3Yc02Dok9QZgfgzx7lUolFIvF\nRb6+vtXjx4/PFolEJdxrJ06ceMLX17c6JCSkXHO7nTt3zhw1atTZ7hZ6oj3WXI1DjbMEYX4MyvWE\nQqGqqKhIXF9f75GYmHhIKpVKJBKJFABg27Ztc+bOnfud5jYXL14Me/fdd1dnZmZqHZk6LS3tt/8l\nEglIJBJT4id0YM3VOJTZE4RhSKVSkEqlXbIvo37+Hh4e9UlJSfvz8/NjJBKJVKlU2u/evXt6QUFB\nm/noKysrg2bMmLFry5YtKcHBwRXa9sUXe6LrscYJxwEAwsMB3nzT0lEQRM9AMxFOT083eV8d2ji1\ntbU+MpnMEwBALpc7Z2ZmTuIqb44cOTIxNDT0UkBAwG1ufZlM5pmUlLR/zZo1y8aOHXvK5MiITmGt\nNo6nJ8Arr1g6CoLofXQo9lVVVf4TJkzIEovFRXFxcWemTJnyQ0JCwlEAgIyMjOQ5c+Zs46//73//\n+/Xy8vKQ9PT0FVxpZm1trU93nQChHWutxiEIwjIIGGPmP6hAwCxx3N7E8uXY7+Cvf8XHa9fipCHr\n1lk2LoIgTEcgEABjTM8krbqxmR60RFu0NdBag41DEIRlILG3Uay1gZYgCMtAYm+jWGunKoIgLAOJ\nvY2irYGWbByC6L2Q2NsolNkTBMGHxN5G0TbhOGX2BNF7IbG3UTQbaKnOniB6NyT2NgrZOARB8CGx\nt1GsdbgEgiAsA4m9jULDJRAEwYfE3kbR1kBLYk8QvRcSexvFWuegJQjCMpDY2yhk4xAEwYfE3kah\nahyCIPiQ2NsoVI1DEAQfEnsbxVrnoCUIwjLoFfumpianuLi4M2KxuEgkEpWkpqauAgBITk7O4Gah\nCg4OruCmKQQAWLVqVerQoUPLRowYcfnw4cNPdfcJENqhIY4JguCj9+fv5OTUlJ2dPd7FxaVRqVTa\nx8fH5+Tk5MRnZGQkc+ssXrx4vaenpwwAoKSkRJSRkZFcUlIiunXrVuDEiROPlJaWDhMKharuPhGi\nLWTjEATBp0Mbx8XFpREAQKFQOLa2ttp5e3vf515jjAl27Ngxi5uHdu/evc/NmTNnm4ODQ8ugQYOu\nDxky5Gpubu7o7guf0AVV4xAEwafDn79KpRJGR0cXlJeXhyxcuPAzkUhUwr124sSJJ3x9fatDQkLK\nAQBu374dMGbMmNPc60FBQZW3bt0K1LbftLS03/6XSCQgkUg6cRqEJlSNQxA9H6lUClKptEv21eHP\nXygUqoqKisT19fUeiYmJh6RSqUQikUgBALZt2zZn7ty53+nbXiAQaJ1ZnC/2RNdDQxwTRM9HMxFO\nT083eV8GV+N4eHjUJyUl7c/Pz48BAFAqlfa7d++enpycnMGtExgYeOvmzZsDuMeVlZVBgYGBt0yO\njjAZGuKYIAg+esW+trbWRyaTeQIAyOVy58zMzElc5c2RI0cmhoaGXgoICLjNrT916tR927dvn61Q\nKBwrKiqCy8rKho4ePTq3e0+B0AbZOARB8NH786+qqvKfP3/+ZpVKJVSpVMKUlJQtCQkJRwEAMjIy\nkrmGWQ6RSFQya9asHSKRqMTe3l756aefLtJl4xDdC1XjEATBR8CY+bVYIBAwSxy3N5GbC/DaawB5\nefh4wACAn38GGDjQsnERBGE6AoEAGGMCU7alHrQ2Cg1xTBAEHxJ7G4VsHIIg+JDY2yhUjUMQBB8S\nexuFqnEIguBDYm+jkI1DEAQfEnsbhYY4JgiCD4m9jcLP7FUqXOzsLBsTQRCWg8TeRuE30La2YlYv\nMKk6lyAIW4DE3kbhZ/Zk4RAEQWJvo/DFnipxCIIgsbdR+A20VIlDEASJvY1CNg5BEHxI7G0UfgMt\n2TgEQZDY2yianj3ZOATRuyGxt1HIxiEIgg+JvY2i2UBLYk8QvRu9Yt/U1OQUFxd3RiwWF4lEopLU\n1NRV3Gsff/zxG6GhoZfCw8MvLFu2bA23/pw5c7ZFRkYWi0SiktWrV7/b3SdAaEcoBGAMe86SjUMQ\nhN58z8nJqSk7O3u8i4tLo1KptI+Pj8/JycmJb2lpcdi3b9/U4uLiSAcHh5aampr+AADbt2+fDQBQ\nXFwcKZfLnUUiUcncuXO/Gzhw4A1znAzRFs7KIRuHIIgObRwXF5dGAACFQuHY2tpq5+XlVff555+/\nmpqausrBwaEFAKB///41AAD+/v5VDQ0Nrq2trXYNDQ2ujo6Oir59+/7avadA6IITe7JxCILoUAJU\nKpUwOjq6oLy8PGThwoWfhYWFXSwtLR12/PjxccuXL1/p5OTUtH79+sUxMTH5iYmJh7Zs2ZLi7+9f\n1djY6PLBBx+86enpKdO237S0tN/+l0gkIJFIuuykCIQv9mTjEETPQyqVglQq7ZJ9dSj2QqFQVVRU\nJK6vr/dITEw8JJVKJUql0r6urs7r9OnTY/Ly8mJnzZq149q1a4O3bt36glwud66qqvK/f/++9xNP\nPHEiISHhaHBwcIXmfvliT3QPXCMt2TgE0TPRTITT09NN3pfB1TgeHh71SUlJ+/Pz82OCgoIqZ8yY\nsQsAIDY2Nk8oFKpqa2t9Tp48+dj06dN329nZtfbv37/m8ccf/zk/Pz/G5OiITkE2DkEQHHrFvra2\n1kcmk3kCAMjlcufMzMxJUVFRhdOmTduTlZU1AQCgtLR0WEtLi4OPj0/tiBEjLnPPNzQ0uJ4+fXpM\naGjope4/DUIbXC9asnEIgtCb71VVVfnPnz9/s0qlEqpUKmFKSsqWhISEo+PGjTu+YMGCTREREecd\nHR0Vmzdvng8A8Kc//emLl19++cuIiIjzKpVKuGDBgk3h4eEXzHMqhCZUjUMQBIeAMWb+gwoEzBLH\n7W0EBADk5QHk5wNs2gSwd6+lIyIIojMIBAJgjJk0DRH1oLVhuAZa8uwJgiCxt2HIxiEIgoPE3oah\nBlqCIDhI7G0YKr0kCIKDxN6GIRuHIAgOEnsbht9ASzYOQfRuSOxtGLJxCILgILG3YbgGWrJxCIIg\nsbdhaNRLgiA4SOxtGLJxCILgILG3YWiIY4IgOEjsbRiycQiC4CCxt2HIxiEIgoPE3obhqnEaGgBc\nXCwdDUEQloTE3obhMvt79wD69bN0NARBWBISexuGa6C9fx/A29vS0RAEYUn0in1TU5NTXFzcGbFY\nXCQSiUpSU1NXca99/PHHb4SGhl4KDw+/sGzZsjXc88XFxZFjx449FR4efiEyMrK4ubm5T3eeAKEb\nLrMnsScIQm+znZOTU1N2dvZ4FxeXRqVSaR8fH5+Tk5MT39LS4rBv376pxcXFkQ4ODi01NTX9AQCU\nSqV9SkrKlq1bt74QERFxvq6uzsvBwaHFPKdCaEJiTxAER4c2jouLSyMAgEKhcGxtbbXz8vKq+/zz\nz19NTU1dxQl5//79awAADh8+/FRkZGRxRETEeQAALy+vOqFQqOrOEyB0wzXQktgTBNFhQZ5KpRJG\nR0cXlJeXhyxcuPCzsLCwi6WlpcOOHz8+bvny5SudnJya1q9fvzgmJia/rKxsqEAgYE8//fTBmpqa\n/rNnz96+ZMmSddr2m5aW9tv/EokEJBJJl50UgdjZYYeq+noAT09LR0MQhLFIpVKQSqVdsq8OxV4o\nFKqKiorE9fX1HomJiYekUqlEqVTa19XVeZ0+fXpMXl5e7KxZs3Zcu3ZtcEtLi0NOTk58fn5+jLOz\nszwhIeHoqFGjzk6YMCFLc798sSe6B3t7zOrd3KjOniB6IpqJcHp6usn7Mrgax8PDoz4pKWl/fn5+\nTFBQUOWMGTN2AQDExsbmCYVCVW1trc+AAQNujhs37ri3t/d9Z2dn+TPPPHOgoKAg2uToiE5hZwdQ\nU0MWDkEQHYh9bW2tj0wm8wQAkMvlzpmZmZOioqIKp02bticrK2sCAEBpaekwhULh6OPjU/vUU08d\nPn/+fIRcLndWKpX2x44dezIsLOyiOU6EaA+JPUEQHHpv7quqqvznz5+/WaVSCVUqlTAlJWVLQkLC\n0XHjxh1fsGDBpoiIiPOOjo6Kb7755kUAbJB9++2334+Njc0TCAQsKSlp/+TJk38yz6kQmtjZAdy9\nS2JPEASAgDFm/oMKBMwSx+1tLFoEUFwMEBQEsH27paMhCKKzCAQCYIwJTNmWetDaMPb2ZOMQBIGQ\n2Nsw5NkTBMFBYm/D2NkB1NWR2BMEQWJv09jZ4V8Se4IgSOxtGBJ7giA4SOxtGK7XLIk9QRAk9jYM\nZfYEQXCQ2NswJPYEQXCQ2NswnNh7eVk2DoIgLA+JvQ1jZwfg6grQh+YKI4heD4m9DWNvTxYOQRAI\nib0NY2dHYk8QBEJib8PY2QH062fpKAiCsAZI7G0YyuwJguAgsbdhgoMBYmMtHQVBENaAXrFvampy\niouLOyMWi4tEIlFJamrqKu61jz/++I3Q0NBL4eHhF5YtW7aGv92NGzcGurm5PdywYcM73RV4V9NV\nk/p2JZ2N6emnAZYu7ZpY+Njie9UdUEyGYY0xAVhvXKaiV+ydnJyasrOzxxcVFYmLi4sjs7Ozx+fk\n5MRnZ2eP37dv39Ti4uLICxcuhC9evHg9f7u33377/aSkpP3dG3rXYo0frDXGBGCdcVFMhkExGY61\nxmUqeqclBABwcXFpBABQKBSOra2tdl5eXnX/+Mc//l9qauoqBweHFgCA/v3713Dr79mzZ9rgwYOv\nubq6NnRf2ARBEIQxdOjZq1QqoVgsLvL19a0eP358dlhY2MXS0tJhx48fHzdmzJjTEolEmp+fHwMA\n8PDhQ7e1a9cuTUtLS+v2yAmCIAjDYYwZtMhkMo+4uLjT2dnZkvDw8PN//vOfP2SMQW5ubmxwcPA1\nxhi8884763fs2PE7xhisWLEibf369e9o2xcAMFpooYUWWoxfDNVszaVDG4fDw8OjPikpaX9+fn5M\nUFBQ5YwZM3YBAMTGxuYJhUJVbW2tT25u7uidO3fOXLp06VqZTOYpFApVzs7O8kWLFn3K35epE+YS\nBEEQpqFX7Gtra33s7e2Vnp6eMrlc7pyZmTlpxYoV6e7u7g+ysrImPPnkk8dKS0uHKRQKRx8fn9rj\nx4+P47ZNT09f4e7u/kBT6AmCIAjzo1fsq6qq/OfPn79ZpVIJVSqVMCUlZUtCQsLRcePGHV+wYMGm\niIiI846OjopvvvnmRXMFTBAEQZiAqf6PKctPP/309PDhwy8PGTKkbPXq1cvMeWz+cuPGjQESiSRb\nJBJdDAsLu/Dhhx/+mTEG9+7d8544cWLm0KFDSydNmnS4rq7O09yxKZVKO7FYXPjss8/+YA0x1dXV\nec6cOfO/I0aMuBQaGlpy+vTpOEvHtHLlylSRSHQxPDz8/Jw5c75ramrqY+6YXnrppU2PPPJIdXh4\n+HnuOX0xrFy5MnXIkCFlw4cPv3zo0KGnzBnX4sWL140YMeJSZGTkuenTp++SyWQe5oxLW0zcsn79\n+ncEAoHq3r173tYQ00cfffTGiBEjLoWFhV1YunTpGkvHdObMmdGxsbG5YrG4MCYmJi83NzfW1Ji6\n7ceguSiVSruQkJCrFRUVgxQKhcPIkSOLSkpKQs11fP5SVVXlV1hYKGaMwYMHD9yGDRt2paSkJHTJ\nkiVr16xZs5QxBqtXr162bNmy1eaObcOGDW/PnTv32ylTpuxjjIGlY3rxxRc3f/nllwsYY9DS0mIv\nk8k8LBlTRUXFoODg4GtNTU19GGMwa9asjK+//nq+uWM6fvz4EwUFBVH8H6auGC5evCgaOXJkkUKh\ncKioqBgUEhJytbW1VWiuuA4fPjyJO96yZctWmzsubTExhklXYmLiwUGDBlVwYm/JmLKyssZPnDgx\nU6FQODDG4O7du/0tHdOTTz4pPXjwYCJjDA4cODBZIpFkmxpTt/0YNJeTJ0+OTUxMPMg9XrVq1bur\nVkzPn4kAAAWXSURBVK1611zH17c899xzezIzMycOHz788p07d3wZwwvC8OHDL5szjps3bwYlJCQc\nycrKGs9l9paMSSaTeXCVVvzFkjHdu3fPe9iwYVfu37/v1dLSYv/ss8/+cPjw4UmWiKmiomIQ/4ep\nK4aVK1em8u9kExMTD546dWqMueLiL7t27Zo+b968reaOS1tMzz///Pfnzp2L5Iu9JWP63e9+t+Po\n0aMTNNezZEyzZ8/elpGRMYsxBt99992cznx2Zhsb59atW4EDBgy4yT0OCgqqvHXrVqC5jq+L69ev\nDyosLIyKi4s7U11d7evr61sNAODr61tdXV3ta85Y3nrrrX+tW7duiVAoVHHPWTKmioqK4P79+9e8\n9NJLX0VHRxf84Q9/+N+GhgZXS8bk7e19/5133tkwcODAGwEBAbc9PT1lkyZNyrT0Zweg+7O6fft2\nQFBQUCW3niW/+5s2bVrwzDPPHLB0XHv37n0uKCioMjIyspj/vCVjKisrG6qt/5AlY1q9evW73Pd9\nyZIl61atWpVqakxmE3uBQMDMdSxDefjwodvMmTN3fvjhh39xd3d/wH9NIBAwc8b8448/PvvII4/c\njYqKKmQ6SlPNHZNSqbQvKCiIXrRo0acFBQXRrq6uDatXr37XkjGVl5eHfPDBB29ev3590O3btwMe\nPnzotnXr1hcsGZM2OorBEvH985///Kujo6Ni7ty53+laxxxxNTY2uqxcuXJ5enr6Cu45Xd95c8UE\ngN/3uro6r9OnT49Zt27dklmzZu2wdEwvv/zylx999NGfb9y4MfBf//rXWwsWLNhkakxmE/vAwMBb\nN2/eHMA9vnnz5gD+lcnctLS0OMycOXNnSkrKlmnTpu0BwGzszp07fgBYifTII4/cNVc8J0+efGzf\nvn1Tg4ODK+bMmbMtKytrQkpKyhZLxhQUFFQZFBRUGRsbmwcA8Pzzz/+3oKAg2s/P746lYsrPz495\n7LHHTvbr1++evb29csaMGbtOnTo11pIxcej6rDS/+5WVlUGBgYG3zBnb119//fsDBw488+23387j\nnrNUXOXl5SHXr18fNHLkyHPBwcEVlZWVQaNGjTpbXV3ta8n3Slf/IUvGlJubO3r69Om7AfD3l5ub\nOxrAtM/ObGIfExOTX1ZWNvT69euDFAqFY0ZGRvLUqVP3mev4fBhjgpdffvlLkUhU8uabb37APT91\n6tR9mzdvng8AsHnz5vncRcAcrFy5cvnNmzcHVFRUBG/fvn32hAkTsrZs2ZJiyZj8/PzuDBgw4GZp\naekwAIAjR45MDAsLuzhlypQfLBXTiBEjLp8+fXqMXC53ZowJjhw5MlEkEpVYMiYOXZ/V1KlT923f\nvn22QqFwrKioCC4rKxs6evToXHPFdfDgwafXrVu3ZO/evc85OTk18eO1RFwRERHnq6urfSsqKoIr\nKiqCg4KCKgsKCqJ9fX2rLfleTZs2bU9WVtYEAAB+/yFLxjRkyJCrx44dexIAICsra8KwYcNKAUz8\n7LqrkUjbcuDAgcnDhg27EhIScnXlypWp5jw2fzlx4kS8QCBQjRw5skgsFheKxeLCn3766el79+55\nJyQkHLFk6SVjDKRS6ZNcNY6lYyoqKhoZExOTxy/bs3RMa9asWcqVXr744oubFQqFg7ljmj179jZ/\nf//bDg4OiqCgoJubNm16SV8M//znP5eHhIRcHT58+GWuusIccX355ZcLhgwZUjZw4MBfuO/6woUL\nPzVnXFxMjo6Ozdx7xX89ODj4Gr/00lIxKRQKhxdeeGFLeHj4+ejo6LPZ2dkSS8TE/07l5eXFjB49\n+szIkSOLxowZc6qgoCDK1JgEjFmdlU4QBEF0MTRTFUEQRC+AxJ4gCKIXQGJPEATRCyCxJwiC6AWQ\n2BMEQfQCSOwJgiB6Af8fEEWldHghXnwAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3327c90>"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now implement the algorithm above to recalculate the SPM global signal for the first volume (volume index zero).  Hint: you will likely need to index using a boolean (mask) array.  Remember, the steps are:\n",
      "\n",
      "* Get a single 3D volume V;\n",
      "* Calculate the mean signal of the voxels in V; call that M;\n",
      "* Make a threshold T where T = M / 8.\n",
      "* Select all voxel values in V that are greater than T; call these values W;\n",
      "* Get the mean of all values in W.\n",
      "\n",
      "You should get the same value as SPM - the first value you read from `global_signals.txt`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Apply algorithm for SPM global calculation to first volume"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make a function called `spm_global` that accepts a 3D array as input, and returns the global signal using the SPM algorithm.  Call that function on the first volume to show that it is working (as in `print(spm_global(data[:, :, :, 0]))`)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make a `spm_global` function that accepts a 3D array as input,\n",
      "# and returns the global mean for the volume according to the SPM algorithm"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make a function called `get_spm_globals` that accepts an image filename as an argument.  The function will load the image, get the array data for the image, use your new `spm_global` function calculate the global value for each volume, and return these values as a list.  Finally, show this is working by plotting the values for the `ds107_sub012_t1r2.nii` image with something like: \n",
      "\n",
      "```\n",
      "all_globals = get_spm_globals('ds107_sub012_t1r2.nii')\n",
      "plt.plot(all_globals)\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Write a function `get_spm_globals` that returns the global values \n",
      "# for each volume"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}