{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Figure 4 Accuracy of the SFM \n",
      "\n",
      "In the notebook, we examne the statistics of the accuracy of the Sparse Fascicle Model (SFM) across the voxels of the white-matter.\n",
      "\n",
      "We start by importing the modules we will use: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nibabel as ni\n",
      "\n",
      "import osmosis.model.analysis as oza\n",
      "import osmosis.model.sparse_deconvolution as ssd\n",
      "import osmosis.model.dti as dti\n",
      "import osmosis.viz.mpl as mpl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We set the data-path to point to the installation-specific location of the data-files"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import osmosis as oz\n",
      "import osmosis.io as oio\n",
      "oio.data_path = os.path.join(oz.__path__[0], 'data')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We generate strings pointing to all of the data-files for this subject:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subject = 'SUB1'\n",
      "data_1k_1, data_1k_2 = oio.get_dwi_data(1000, subject)\n",
      "data_2k_1, data_2k_2 = oio.get_dwi_data(2000, subject)\n",
      "data_4k_1, data_4k_2 = oio.get_dwi_data(4000, subject)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Wewill only look at white-matter voxels. This loads the white-matter mask:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "wm_mask = np.zeros(ni.load(data_1k_1[0]).shape[:3])\n",
      "wm_nifti = ni.load(oio.data_path + '/%s/%s_wm_mask.nii.gz'%(subject, subject)).get_data()\n",
      "wm_mask[np.where(wm_nifti==1)] = 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following are the parameters defining the performance of the Elastic Net solver , based on the calculations in AppendixA\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# This is the best according to rRMSE across bvals: \n",
      "l1_ratio = 0.8\n",
      "alpha = 0.0005 \n",
      "solver_params = dict(l1_ratio=l1_ratio, alpha=alpha, fit_intercept=False, positive=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Initialize the class instances of the SFM objects, using AD/RD derived from the corpus callosum. Parameters are not saved (hence `temp`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ad_rd = oio.get_ad_rd(subject, 1000)\n",
      "SD_1k_1 = ssd.SparseDeconvolutionModel(*data_1k_1, mask=wm_mask,  params_file ='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'], solver_params=solver_params)\n",
      "SD_1k_2 = ssd.SparseDeconvolutionModel(*data_1k_2, mask=wm_mask,  params_file ='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'], solver_params=solver_params)\n",
      "ad_rd = oio.get_ad_rd(subject, 2000)\n",
      "SD_2k_1 = ssd.SparseDeconvolutionModel(*data_2k_1, mask=wm_mask,  params_file ='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'], solver_params=solver_params)\n",
      "SD_2k_2 = ssd.SparseDeconvolutionModel(*data_2k_2, mask=wm_mask,  params_file ='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'], solver_params=solver_params)\n",
      "ad_rd = oio.get_ad_rd(subject, 4000)\n",
      "SD_4k_1 = ssd.SparseDeconvolutionModel(*data_4k_1, mask=wm_mask,  params_file ='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'], solver_params=solver_params)\n",
      "SD_4k_2 = ssd.SparseDeconvolutionModel(*data_4k_2, mask=wm_mask,  params_file ='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'], solver_params=solver_params)\n",
      "\n",
      "#TM_4k_1 = dti.TensorModel(*data_4k_1, mask=wm_mask)\n",
      "#TM_4k_2 = dti.TensorModel(*data_4k_2, mask=wm_mask)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvecs\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "CallculaterRMSEfor each b value , from a cross-prediction: fit the model to one data-set in each b value and predict the other data set. Calculations average across both directions (1=>2"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rrmse_1k = oza.cross_predict(SD_1k_1, SD_1k_2)\n",
      "rrmse_2k = oza.cross_predict(SD_2k_1, SD_2k_2)\n",
      "rrmse_4k = oza.cross_predict(SD_4k_1, SD_4k_2)\n",
      "#rrmse_tensor_4k = oza.cross_predict(TM_4k_1, TM_4k_2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r",
        "SparseDeconvolutionModel.model_params [****************100%******************]  64463 of 64464 complete "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = mpl.probability_hist(rrmse_1k[np.isfinite(rrmse_1k)], label='b=%s'%str(1000), color=[0.8, 0.8, 0.8])\n",
      "fig = mpl.probability_hist(rrmse_2k[np.isfinite(rrmse_2k)], fig=fig, label='b=%s'%str(2000), color=[0.59, 0.59, 0.59])\n",
      "fig = mpl.probability_hist(rrmse_4k[np.isfinite(rrmse_4k)], fig=fig, label='b=%s'%str(4000), color=[0.32, 0.32, 0.32])\n",
      "# Add one of the tensor curves from Figure 2 and put it in the background as reference:    \n",
      "# fig = mpl.probability_hist(rrmse_tensor_4k[np.isfinite(rrmse_4k)], fig=fig, color='gray', label='Tensor model at b=%s'%str(4000))\n",
      "#fig.set_size_inches([10, 8])\n",
      "fig.axes[0].plot([1,1], fig.axes[0].get_ylim(), '--k')\n",
      "fig.axes[0].plot([1/np.sqrt(2),1/np.sqrt(2)], fig.axes[0].get_ylim(), '--k')\n",
      "fig.axes[0].set_xlim([0.6,1.4])\n",
      "plt.legend()\n",
      "\n",
      "fig.savefig('figures/Figure4_histogram.svg')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD9CAYAAABOd5eOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl002W+P/B39jZt0iTdl3SlhZalVOCyY3FEkTsgKnpF\nrguOjqPHuaLj9hs9ss24XL3XI3pHEQdFcRkV5wKjcD2IBQShQrVgCy0tXemaNM3SJm2afH9/dPI1\nIWmTpmnWz+ucntLv95tvPk3aTx+e5/M8D4dhGAaEEEJCGjfQARBCCBk/SuaEEBIGKJkTQkgYoGRO\nCCFhgJI5IYSEAUrmhBASBkZN5vfeey+Sk5Mxffp09tgTTzyBwsJCFBcX4+abb4ZWq53wIAkhhIxu\n1GS+fv16HDx40OHYddddh6qqKlRWVqKgoAAvvPDChAZICCHEvVGT+eLFiyGXyx2OLVu2DFzu8MPm\nzp2L1tbWiYuOEEKIR8bVZ75z506sWLHCV7EQQgjxEt/bB/75z3+GUCjEHXfc4fI8h8PxOihCCIlk\n3qyy4lXL/L333sNXX32FDz/80G1Awf6xceNGl/8Oto9gji3U4rz66qsDHkO4vJYUp+8/vDXmZH7w\n4EG8/PLL2Lt3L6Kiorx+4mC0efPmQIdA/ODIkSOBDoEQnxs1ma9duxYLFixATU0NlEoldu7cid//\n/vcwGAxYtmwZSkpK8NBDD/krVkIIISMYtc/8448/djp27733TlgwgVBaWhroEDxCcUaeUHktKc7g\nwGHG00kz2o05nHH1/wRCKMZMxo7eZxLMvP359LqahRASvhQKBTQaTaDDCGtyuRw9PT0+ux8lczsb\nN24MdAjED+h9dk+j0dD/XiaYr8u3qZuFEOKEfn8n3kivsbevPa2aSAghYYCSOSGEhAFK5oSQkJGd\nnY1vvvkm0GEEJUrmhJCQweFwxjVwaDabsWbNGuTk5IDL5bqcDfzUU08hISEBCQkJePrppx3ONTY2\nYunSpYiJiUFhYaHTH5aPPvoIWVlZiI2NxU033eTXiiBK5nY2bdoU6BCIH9D7HNmWLFmC3bt3IyUl\nxekPw/bt27F3716cPXsWZ8+exf79+7F9+3b2/Nq1azFr1iz09PTgz3/+M9asWQOVSgUAqKqqwu9+\n9zt8+OGH6OzshFgs9u8MeWaCTOCtJ0woxkzGjt5n94L1NcrOzmZeeOEFpqioiJHL5cz69esZk8nk\n1b0yMjKYI0eOOBybP38+s2PHDvbrnTt3MvPmzWMYhmFqamoYkUjEGAwG9vySJUuYt956i2EYhvl/\n/+//MevWrWPP1dfXM0Kh0OF6eyO9xt6+9tQyJ4SEDIZh8NFHH+Hrr79GfX09amtr8ac//QktLS2Q\nyWSQy+UuPz755BOP7l9dXY3i4mL26xkzZqCqqgrAcMs7NzcXMTEx7Pni4mKH8/aPzc3NhUgkQm1t\nrS++dbdo0hAhZMzOnDnjk/vMmjVrTNdzOBw8/PDDSE9PBwA888wz+P3vf4+tW7eit7d33PEYDAbE\nxcWxX0ulUhgMBpfnAEAikaC9vR0A0NfX53ReKpVCr9ePOy5PUDInhIzZWJOwLymVSvbfmZmZaGtr\n89m9Y2NjodPp2K+1Wi1iY2NdnrOdl0gk7PkrN7i3Pz/RqJuFEBJSmpubHf6dnp6OlpYWxMbGQiKR\nuPxwtQKsK1OnTsVPP/3Efl1ZWYlp06ax5y5dusS21G3np06dyp6vrKxkz9XX12NwcBAFBQXj+n49\n5lVPuwcm8NYTZuPGjYEOgfgBvc/uBevvb1ZWFjN9+nSmtbWVUavVzMKFC5lnnnlmTPcwmUyM0Whk\nMjIymK+//poxGo3subfeeospLCxkLl++zLS2tjJFRUXM9u3b2fPz5s1jHn/8ccZoNDJ79uxhZDIZ\no1KpGIZhmKqqKkYqlTLHjh1jDAYDs3btWmbt2rUjxjHSa+zta0/JnBDiJFh/f7Ozs5kXX3yRKSoq\nYmQyGXPPPfc4JGNPZGVlMRwOh+FyueznpqYm9vyTTz7JKBQKRqFQME899ZTDYxsbG5nS0lImOjqa\nmTJlCvPNN984nP/oo4+YzMxMJiYmhlm9ejWj0WhGjMPXyZwW2iKEOKHf34lHC20RQghxQsmcEELC\nACVzQggJA5TM7dCaHZGB3mcSjmgA1E4oxkzGjt5n9+g1mng0AEoIIcQJJXNCCAkDlMwJISQMUDIn\nhIQM2jZuZJTM7WzcuDHQIRA/oPc5dI1327iTJ09i2bJliI+PR1JSEm677TZ0dHQ4XEPbxoUBKlmL\nDPQ+R67e3l787ne/Q1NTE5qamiCRSLB+/Xr2fNhuG7d+/XomKSmJmTZtGntMrVYz1157LZOfn88s\nW7ZsxIVk3NyaEBLEgvX315fbxjEMw5w5c4aRSCTs12G7bdz69etx8OBBh2Mvvvgili1bhtraWvzq\nV7/Ciy++OGF/aEINwzA4c+YMPv30U+zfvx/fffcd1eoS4kOMj7eNO3r0KLteORDG28YtXrwYjY2N\nDsf27duHI0eOAADuvvtulJaWUkIHMDQ0hG+//RY6nQ5LliyB2WzG0aNHUVBQgKSkpECHR4hPvfnm\nmz65z4MPPjim6325bdzZs2exdetW7Nu3jz0WUdvGdXZ2Ijk5GQCQnJyMzs5OnwcVir7++msIBALc\neOON4POHX9a8vDw0NDRQMidhZ6xJ2Jd8sW1cXV0dVqxYgW3btmHhwoXs8YjdNs7dyPKmTZvYj7Ky\nsvE8lV94OzDW19eH9vZ2LF26lE3kAJCTk4OGhgYfRUd8hQZAQ9t4t41ramrCsmXL8Nxzz2HdunUO\n9w7EtnFlZWUOudJr7jrVGxoaHAZAJ0+ezLS3tzMMwzBtbW3M5MmTfdqJH0jexlxZWckcOnTI6bjV\namV27do16m4jxP9C8WfT34L1NRrvtnGtra1Mbm4u88orr7g8H8rbxo25Zb5q1Srs2rULALBr1y6s\nXr3a+78kYaK+vh6TJk1yOs7hcJCdnY1Lly4FICpCwg+Hw8G6detw3XXXIS8vD/n5+Xj22Wc9fvw7\n77yDhoYGbNq0iW21S6VS9vwDDzyAlStXYvr06ZgxYwZWrlyJ3/72t+z5Tz75BKdPn4ZCocAzzzyD\nPXv2ID4+HgBQVFSEt956C+vWrUNycjKMRiP+8pe/+O6bd2PUVRPXrl2LI0eOQKVSITk5GVu2bMGN\nN96I2267Dc3NzcjOzsann34KmUzmfOMQXHXNm5gNBgM+++wz3HXXXeDxeE7nW1paUF5ejltuucVX\nYZJxCsWfTX+j12ji+XrVxFEHQO37mewdOnRozE8Ururq6pCTk+MykQNAWloatFotDAYDO5BCCCG+\nRjNAx6m+vh55eXkjnufxeMjMzHQq8SSEEF+iZG5nrGt2GAwG6HQ6tuZ1JDk5OZTMgwitzULCESVz\nO2MtC9JoNIiPjweXO/rLqFQq0dHRgcHBwXFER3yFShNJOKJkPg5ardZpxpcrQqEQKSkpDvWxhBDi\nS5TMx6G3t9ejZA5QVwshZGJRMh8HT1vmwPCi+s3NzbBYLBMcFSEkElEyH4exJPOYmBjExcWxi/IQ\nQogvUTK3M5aBMavVCoPB4DB7zJ3s7GxaqyUI0ABo6KJt40ZGydzO5s2bPb5Wr9dDLBY7LKzlTn5+\nPurq6vy2JCZxbSzvMwku4902zt6WLVvA5XJx+PBhh+O0bVyEGamLRa1Wo76+HmfPnnXqH5dKpZg+\nfTqOHz/urzAJIS7U19fj888/R1pamsPxUN42jpK5l66sZNHr9Xj11Vfx2GOP4c0338S2bdvw/fff\nOz2upKQEPT09VNlCiJfKy8sxdepUKBQK3HvvvRgYGBjzPR5++GG89NJLEAgEDsd37dqFxx9/HGlp\naUhLS8Pjjz+O9957DwBQW1uLH3/8EZs3b4ZIJMLNN9+MGTNmYM+ePQCADz/8EKtWrcKiRYsQExOD\nrVu34osvvkBfX9+4v2dPUDL3kn3LvLq6Go8++ijkcjnefvttvPLKK7jzzjvZHZns8Xg8LFmyBN99\n9x1VthAyRowPto377LPPEBUVhRtuuMHp/mG7bRwZmVarZXc82bNnD+644w5cc8017Pm5c+finXfe\nQW9vr9OqkhkZGeDxeNBqtVAoFH6NmxBf8NUqoLZWrafGu22cXq/HM888M+JigRG1bVw4G8uaHVqt\nFjKZDBaLBbW1tdiwYYPD+aioKMyZMwffffcdfv3rXzs93vYmUzL3P1qbZfzGmoR9aTzbxm3atAl3\n3nknMjMz2WP2y81G7LZx4cbTkjWLxQKDwQCJRIJLly4hMTHR5RtWWlrqsqsFGP6LTlUtgUGliaFt\nPNvGHT58GNu2bUNqaipSU1PR0tKC2267DS+//DKAwGwb5zNe7U/kgQm8dcBpNBpm9+7dDMMwzN//\n/ndmx44dLq8bGhpi7rvvPqa5udnpXEVFBXP8+PEJjZMQbwXr7+94t41Tq9VMZ2cn09nZyXR0dDBK\npZL5/PPPmb6+PoZhImzbOOJYyVJdXc3+Zb4Sj8fD4sWLXfbPUcuckLEb77ZxCoUCSUlJSEpKQnJy\nMng8HuRyOcRiMYAw3jZuXDcO422nKisrodPpsGDBAtxzzz144403RpzW39PTgz/84Q/YsmWLQ19f\nZ2cnjh07hjVr1vgrbEI8Fs6/v8HC19vGUcvcCxqNBnK5HI2NjVAoFKOuz6JQKHDbbbdh+/btsFqt\n7HFqmRNCfImSuR1PB8a6u7uRmJiIqqoqFBUVub3+uuuuw9DQEL799lv2WHR0NIaGhmA2m70Nl3iJ\nBkBJOKJkbseTNTssFgt6e3uhUChQXV3NjnSPhsfj4YEHHsDu3bvZ2WocDgexsbHUOg8AWpuFhCNK\n5mPU09MDiUQCHo+H8+fPe9QyB4Y3p8jOzkZFRQV7jLpaSLCSy+Xsolb0MTEfcrncp+8ZTRoaI5VK\nhcTERLS3t0MsFo/pDVm0aBGOHTuG+fPnA/Dv7DBCxqKnpyfQIZAxopb5GKlUKiQkJKCurg6TJk0a\n02PnzZuHc+fOsQvvUMucEOIrlMzHaDzJPCYmBtOnT8epU6cAUDInhPgOJXM77tbssFqtbDK/ePEi\n8vPzx/wcixcvxrFjxwBQMg8UWpuFhCNK5nbclaxptVqIxWJwuVw0NTUhJydnzM9x1VVXob6+HhqN\nhpJ5gFBpIglHlMzHwNYqb25uRkpKCqKjo8d8D5FIhNmzZ+PkyZOIjo7G4OAg1ZoTQsaNkvkYdHd3\ne91fbm/OnDn44YcfwOFwqHVOCPEJr5P5Cy+8gKlTp2L69Om44447vNq6KdTYyhIvXrw4rmQ+c+ZM\n1NbWor+/n5I5IcQnvErmjY2N2LFjByoqKnDu3DlYLBaHbZnCVU9PD+Lj41FXV+fV4KdNdHQ0Jk+e\njMrKSkrmhBCf8CqZS6VSCAQC9Pf3Y2hoCP39/ew2TqFstIExW982h8NBV1eXw04l3pg9ezZOnz6N\nuLg4p91JyMSiAVASjryaAapQKPCHP/wBmZmZiI6OxvXXX49rr73W6Tr7X5rS0lKUlpZ6G6dfbN68\necRfdNsGzpcuXUJWVhb4/PFNnp09ezY+/fRTrFixYkzbXpHxG+19JsTfysrKUFZWNu77eLWeeX19\nPVauXIljx44hLi4Ot956K9asWYN169b9cuMQXA95tJjr6+tx8eJF6HQ6GI1G3H333W7vNzQ0hLa2\nNlgsFlitVjAMA6vVCqvVipSUFGzZsgVr167FxYsXcccdd/j62yEjCMWfTRI5/Lqe+enTp7FgwQLE\nx8eDz+fj5ptvxokTJ7y5Vciwtcx//PFHlJSUePQYtVoNk8kEqVQKuVyOhIQEpKSkICEhAa2trZg9\nezbOnz8Pg8EAi8Uywd8BISSceZXMp0yZgpMnT8JoNIJhGBw6dMjj1QNDlU6nA5/PR0tLC6ZMmeLR\nY3p6epCSkoL4+HgoFArIZDJIpVLEx8eDx+OhoKAAZ8+edbnrNyGEjIVXyby4uBh33XUXZs+ejRkz\nZgCAwz554Uir1aK9vR1FRUUQCoVurzcajTCbzZBIJE7nOBwOUlJSIBKJcPnyZcTGxqK3t3ciwiaE\nRAivR/GefPJJPPnkk76MJeBGW7NDq9VCrVaPqYtFoVCAw+G4PC+TydDa2ork5GSYzWaqaPEjWpuF\nhCOaAWpnpAqHoaEhGI1G/Pzzzx4lc4Zh2Jr0kXA4HCQnJyM5ORl6vZ5a5n5ElSwkHFEy94BOp4PF\nYkF0dDRSUlLcXq/X6yEQCNyu3ZKQkIDExER0dnZSMieEjAslcw/odDpotVrMnDnTo+ttXSzucLlc\nFBUVoaWlhZI5IWRcKJl7QKvVoqurC8XFxW6vtVgs0Gq1HiVzYHidlu7ubgwMDETE+jaEkIlBydwD\nWq0WGo0G2dnZbq/t7e1FbGwsBAKBR/cWi8VISkqC2Wym1jkhxGuUzO2MNDCmVqsxMDAw6oCm/bWe\ntspt8vPzodfraRNdP6EBUBKOKJnb2bx5s8vjbW1tSEpKApc7+ss1ODiI/v5+yGSyMT1vQUEBDAYD\n2tvbx/Q44p2R3mdCQhklczcsFgtUKhUyMjLcXtvT0wO5XO426V8pNzcXvb29UKlUtGYIIcQrlMzd\nMBgMMJvNbpf4ZRjGqy4WAMjKyoJarUZfXx8MBoO3oRJCIhglczd0Oh1MJhNSU1NHva6/vx9WqxWx\nsbFjfg6RSISUlBT09PSgq6vL21AJIRGMkrkb/f396OvrQ1pa2qjX2WZ8jjR935309HQIhUI0NTV5\n9XhCSGSjZG7H1ZodfX190Gq1bpN5b28v5HK518+dmpoKLpeLnp4emEwmr+9D3KO1WUg4omRux1XJ\nmlqtBpfLdbn6oY3ZbIbFYkFUVJTXz52WlgaTyQSdTkf15hOMShNJOKJk7kZ7ezsSExNHvaa/vx9i\nsdjrLhZguGWuUqkgEAhoGzlCyJhRMneju7sbycnJo17T19eHmJiYcT1PWloa2trakJmZifb2dgwN\nDY3rfoSQyELJ3A2VSuW2LNEXyVwqlcJqtSI+Ph4Gg4F2HiKEjAklcze0Wi0yMzNHPM8wDPr6+iAW\ni8f1PBwOB6mpqeBwONDpdDS1nxAyJpTM7Vw5MGY2m9HX1zdqMh8cHASXy/VoKzl30tLSoFKpIJfL\n0dbWRrNBJwgNgJJwRMnczpVrdvT19aG/v3/UCUO+6GKxSU1NRVtbG+RyOYaGhmg26AShtVlIOKJk\nPorLly9DJBKNWnLoy2SelpaG9vZ2SKVSMAxD+4ISQjxGyXwUbW1tkEqlo17ji/5yG1vLXCqV0ibP\nhJAxoWQ+iu7ubsTFxY14nmEYGI1Gn3az2FrmfX19GBoaot2HCCEeoWQ+Cttg5EiMRiOEQiF4PJ5P\nni8mJgZCoRAMw0Cn0yEuLo5mgxJCPELJ3M6Va3bYFs8aiS/7y23S0tLQ09MDi8WC6Oho6mqZALQ2\nCwlHlMztXFmyptPpRp3Kb5vG70tpaWno6Ohg++r7+vpgsVh8+hyRjkoTSTiiZD4KvV6PlJSUEc9P\nVMu8ubkZUqkUBoMBsbGxNBuUEOIWJfNR9PX1jVhjbrFYMDAwgOjoaJ8+Z0lJCcrLyyGRSKDVaqnf\nnBDiEUrmIxgYGMDg4OCILfP+/n5ER0ePeb9Pd7KysiAUCqHX66HT6SCTyaDVamk2KCFkVF5not7e\nXqxZswaFhYUoKirCyZMnfRlXwHV0dCAqKgoCgcDl+YnoLweG12hZsGAB6urqoNPpIBQKIRKJaDYo\nIWRUXifzRx55BCtWrMD58+dx9uxZFBYW+jKugLAfGGtraxu1P3wi+sttFixYgMrKSraSRSaTQaPR\nTMhzRSIaACXhyKtkrtVqcezYMdx7770AAD6fP+rkmlBhv2ZHZ2fnqJszT2QyVyqViI6ORltbGywW\nC2QyGXp7e6mrxUdobRYSjvjePKihoQGJiYlYv349KisrMWvWLLz22mtO3Q72LaDS0lKUlpaOJ1a/\n6u7uhkwmc3nOtk2cSCSakOfmcDhYuHAhKioqoNfrERcXBy6Xi/7+/gn7A0IICYyysjKUlZWN+z4c\nxovm3unTpzF//nycOHECc+bMwYYNGyCVSrFly5ZfbszhhFxL0j7mF154ARKJBA8//LDTdVqtFp2d\nnSgoKJiwWJqbm/Hss89i69atyMrKQmtrKzgcjtuNMoh7ofizSSKHtz+fXnWzZGRkICMjA3PmzAEA\nrFmzBhUVFd7cKmj19vYiISHB5bmJ7GKxUSqVsFqtaG1tBQC2q4UQQlzxKpmnpKRAqVSitrYWAHDo\n0CFMnTrVp4EFmlarHXHvT38kc1srvK6uDsDwui0WiwUmk2lCn5cQEpq86jMHgNdffx3r1q3D4OAg\n8vLy8O677/oyroCwX7NDr9cjLS3N6RrbNnHZ2dkTHk9OTg4aGxsBDCd32wSi0WalEvdobRYSjrxO\n5sXFxfjhhx98GUvA2QZszWYzBgYGXPZP27aJG6n+3JeKiorw008/sV/LZDK0t7dTMh8nKk0k4Yhm\ngLrQ3d0NkUjksivFH10sNlOnToVWq0VfXx8AQCKRwGQywWw2++X5CSGhg5K5C62trRCLxeBwOE7n\njEbjhMz8dCU+Ph48Hg/19fUAAC6XC6lUSgOhhBAnlMxdGG27OKPROOqeoL6WmpqKqqoq9mu5XE7J\nnBDihJK5C11dXSNOGDKZTD5fKXE02dnZuHTpEvu1bWlcWuOcEGKPkrkd28BYV1eXy00prFYrBgcH\nJ2zmpytTpkzB5cuX2a95PB67PC7xDg2AknBEydyObc2Ozs5OKJVKp/MmkwkikchlX/pEmTFjBtRq\ntUNLnCYQjQ+tzULCESXzKzAMA7VajdzcXKdzRqPRr10sAJCYmAiBQMBOHgKAuLg46HQ6WK1Wv8ZC\nCAlelMyvYNsIwtWEIZPJ5NfBT5vk5GScO3eO/VogECA6Ohp6vd7vsRBCghMl8ys0NTUhJibG5fK3\n/h78tMnMzGTLE23i4uKo35wQwqJkfoW6ujrIZDLweDync/4uS7TJz89HS0uLwzGpVErbyRFCWJTM\n7WzcuBGNjY2jVrIEIplPmzYNKpXKIXFHR0eDYRgMDAz4PZ5QR2uzkHBEydzOpk2b0Nra6nLtk0BU\nstgolUpwuVyH1rlt4S3qahk7Kk0k4YiSuR2GYdDZ2YnMzEync4HqLweGE3dCQgJ+/vlnh+NSqRQ6\nnS4gMRFCggslcztqtRo8Hs9lN0ug+stt0tPT2fXjbWg2KCHEhpK5nebmZshkMkgkEqdzgWyZA0Bu\nbi6ampocjvF4PIjFYipRJIRQMrfX3NwMsVjssiwx0C3zwsJCdHZ2OlWv2CYQEUIiGyVzO2+//TZE\nIpFTyzwQa7JcKScnBwCgUqkcjtsGQalE0XM0AErCESVzOwcOHIBMJoNQKHQ4PjAwAJFIBC43cC+X\nWCyGVCrF+fPnHY5HRUVRieIY0dosJBxRMv8nW8s2NTXV6Vygu1iA4YqWlJQUp2ROJYqEEICSOaun\npwfA8O4+Vwr04KdNdnY2Ll686HSckjkhhJL5P3V0dACAyx2GgqFlDgwvh9vW1ua0/K1EIkFfXx+V\nKBISwSiZ/1N7ezsA18k8WFrmKSkpSE9Px3fffedwnMfjISYmhkoUCYlglMz/qb29HUuXLnVK5lar\nlR0ADTS5XI6EhAQcPXrU6Zxt4S3iHq3NQsIRJfN/6ujowHXXXeeUzAcGBiAUCgNayWITHR2N5ORk\nqFQqh63kgF/qzalE0T0qTSThKPAZKki0tbWx3RX2ArG70Gji4+Mxc+ZMHDt2zOG4rU/fZDIFIixC\nSIBRMsdwV0pHRwcSEhKc1jEP1O5CI0lJSYFSqcSRI0ccWuEcDocW3iIkglEyB6DRaCASiVyWJQZb\nyzwzM5OdINTa2upwjkoUCYlclMwxPPgpl8tHrGQJppZ5YmIijEYjlEolmpubHc5RiSIhkWtcydxi\nsaCkpAQrV670VTwB0dHRAYlEgo8//tjhuK2SJZiSOZfLRUZGBiQSictVFCUSCbXO3aABUBKOxpXM\nX3vtNRQVFQVk9x1fam9vh1gsxttvv+1wPJgqWexlZmaCz+c7tcwBQCaTOU0qIo5obRYSjrzOUq2t\nrfjqq69w3333hXw5XHt7OwQCgdPxYJksdCWlUgmz2ezUMgd+KVG0Wq0BiIwQEih8bx/46KOP4uWX\nXx61esL+v7OlpaUoLS319ukmVEdHBzIyMpyOB8s0/iuJxWIkJyfj+++/dxqgFQgEiIqKgl6vR1xc\nXACjJIR4oqysDGVlZeO+j1fJ/B//+AeSkpJQUlIyahCh0DdpK0ssKChwOmcymSCTyQIQlXvZ2dmQ\ny+VoaWlxit3W1ULJnJDgd2VD19tuQK+6WU6cOIF9+/YhJycHa9euxeHDh3HXXXd5FUCg2coSFQqF\n07lgbZkDw3uCxsbGuuxqsSXzUO/+IoR4zqtk/vzzz6OlpQUNDQ345JNPcM011+D999/3dWx+0d7e\nDoVCAalU6rBmh23Dh2BN5vHx8RAKhS6TeVRUFPh8Pvr6+gIQWfCjtVlIOPJJmUYoV7N0dXVBKpVC\nKpU6dAuZTKagrGSxiY6ORlxcHBoaGlyel8lkVKI4glDo/iNkrMadqa6++mrs27fPF7EEhEqlglAo\ndOpfDtZKFnuZmZlobm522Z1CJYqERJbgbHb6kUqlAp/Pd0rmwdxfbpOWlgZguN//SmKxGBaLhRbe\nIiRCUDL/5273rpJ5sLfM4+PjoVAoXE4e4nA41DonJIJEfDLv7u4OiaVvXVEoFBCLxaP2m1MyJyQy\nRHQyZxgGarUaCQkJ4HK57MCY1WrF4OBgUOwuNBqFQoG4uDgcPnzY5YzP2NhYmEwmmM3mAEQXvGgA\nlISjiE7mBoMBHA4HCQkJAH4p1jeZTBCJREFbyWIjFAqRlpYGoVCIM2fOOJ3ncrmIi4uj1vkVaG0W\nEo6CO1tNMLVaDYlEEpL95Tbx8fFYuHAh9u7d6/I8JXNCIkNEJ3OVSgWxWBySZYk2CoUCSqUSarUa\ntbW1TudwPiWTAAAaqUlEQVTj4uJgMBhojXNCwlzEJ3ORSBSSZYk2CoUCvb29+PWvf+2ydc7j8RAb\nG0sTiAgJcxGfzHk8Xkh3sygUCvT09GDp0qWoqKhwOdgpl8vR09MTgOgIIf4S0cncVpYYGxsLYHjN\nDovFArPZHPSVLDZyuRxarRYikQjJyckua87lcjkMBgNVtfwTrc1CwlFEJ/POzk4oFAq2amXTpk3s\nnp+hst4Mn8+HXC5HV1cXJk2ahLq6OqdreDweZDIZ1Gp1ACIMPlSaSMJRRCdztVqNpKQkh2Oh1MVi\no1Qq0dLSMmIyB4arXtRqNS2LS0iYithkbrFY0Nvbi9TUVIfjoZ7M6+vrXV4TGxsLhmFoWVxCwlTE\nJnOtVouoqCinTSlCqSzRJiUlBRqNBklJSWhvb3e5uJZtchR1tRASniI2matUKnZNcHuhVJZow+Px\nkJaWhs7OTiiVyhHXalEoFNBoNFRzTkgYiuhkfmWN+XPPPQeLxQKhUBjAyLzjSb+5UChEbGysyyVz\nIwkNgJJwFLHJvLu7G3w+HxKJhD22detWREdHh0wliz2lUonW1lbk5eWNmMyBXwZCIxmtzULCUcQm\n8/b2dkilUvB4PIfjodbFYhMXFwcul4vExMQRB0Ft15lMJtq0gpAwE7HJvLOzE/Hx8U7HQ23w04bD\n4SAzMxODg4PQaDQwGAwur+NyudQ6JyQMRWwyv3z5MjIyMpyOh2oyB8CWJubl5aG8vHzE66jmnJDw\nE5HJvK+vDzqdDtnZ2U7nQrWbBQBSU1MxODiI1atX44MPPkB1dbXL66KjoyEUCqHT6fwcISFkokRk\nMr906RIUCgVkMhl7zGw244EHHoBAIAhgZOPD4XBQUFAAk8mEDRs24JVXXkFLS4vLa+Pj49n9TyMN\nrc1CwlFEJvO6ujpIpVKHskSTyYRHH300JCtZ7OXn56Ourg7Tp0/H2rVr8eabb7q8TqFQQK/XY3Bw\n0M8RBh6VJpJwFJHJvL6+HtHR0ZBKpeyxUJws5IpCoUB0dDTa2tqwdOlSdHZ2oq2tzek6Ho+HxMTE\nEVvuhJDQEpHJ/OLFi0hISHDoUgnFNVlGUlBQgJqaGvD5fCxevBhHjhxxeV1qaiqMRmPETyIiJBxE\nXDLXarXo6+tDWlqaw/FwSub5+flobGzE4OAgrr76apSVlcFqtTpdx+VykZWVhZaWFgwNDQUgUkKI\nr0RcMr906RJSU1Md+ssZhgnJBbZGIhaLkZaWhvr6euTk5CAmJmbEyhaJRAKZTIbW1lY/R0kI8aWI\nS+Z1dXVITEx06C83m83gcDj405/+FMDIfKuwsBDnz58HAJSWlqKsrGzEa9PT06HT6SKmVJEGQEk4\n8iqZt7S0YOnSpZg6dSqmTZuGbdu2+TquCVNfX+9UyWLrYgmnNTsyMzOh1+uh0WiwePFilJeXY2Bg\nwOW1PB4PmZmZaG5udtkdE27C6X0mxMarZC4QCPDqq6+iqqoKJ0+exP/8z/+wrcBgV1dXh6ioKIeW\neTh1sdhwuVxMnjwZ58+fh1wuR05ODs6ePTvi9TKZDGKx2GXlCyEk+HmVzFNSUjBz5kwAwzvYFBYW\nhkQSUKvVsFgsGBoacmqZh0NZ4pWmTJmC2tpaWCwWzJo1C2fOnBn1eqVSCbVaTbsRERKCxt1n3tjY\niB9//BFz5871RTwTqrq6GgUFBeBwOA7JW6/XIzY2NoCRTQyZTAaZTIbm5mbMnj0bZ86cGXU9FoFA\ngIyMDDQ1NUVEdwsh4YQ/ngcbDAasWbMGr732mstkaD/QVFpaitLS0vE83bhVV1cjPT3dYRr/wMAA\nrFZrWLbMgeGa89raWlx//fUQiURoaGhAbm7uiNfbdiNqb29Henq6HyMlJDKVlZWNWqDgKa+Tudls\nxi233IJ///d/x+rVq11eE2xVA9XV1SgtLXWoMdfr9ZBIJOBwOGG5ZkdeXh6+//57mEwmzJ49G6dP\nnx41mXM4HGRlZaG6uhpxcXFh+T+WcHyfSei6sqHr7QC9V90sDMPgN7/5DYqKirBhwwavntjftFot\nNBoNhoaGHJK5TqdjB0OD7Y+PL4hEIiiVStTX17NdLe4IBAJkZmaisbExLLtbwvF9JsSrZH78+HHs\n3r0b3377LUpKSlBSUoKDBw/6OjafsvWX9/b2IikpCcDwHyVbyzycTZ48GTU1NZgyZQra29s9mr4v\nl8sRExNDk4kICRFedbMsWrQo5Fps1dXVyMjIgEwmA58//G2bTCbweDyIRKIARzexMjIy8O2338Jg\nMKC4uBinTp3C8uXL3T5OqVSiuroaMpnMoZSTEBJ8ImYGqK0P+MoulnBvlQPDk4Ly8/NRU1ODVatW\n4bPPPvNotiefz0dWVhaamppgsVj8ECkhxFsRkcwNBgM6OjpgtVpdDn5GgqKiIpw/fx45OTlYvHgx\n3n33XY8eFxcXB6lUSkvlEhLkIiKZnz9/Hvn5+ejp6UFycjKA4f5yg8HgkMzDeWBMLpcjKSkJtbW1\nuP3221FTU+PRYCgw3E1jWxogHITz+0wiV8Qkc6VSifj4eHYN876+PgiFQoc1zcN9zY4ZM2bg7Nmz\nEIlEePDBB7F9+3b09va6fRyPx0Nubi6am5thNBr9EOnECvf3mUSmiEnmMpnMqYsl0gb10tPTweVy\n0dLSgunTp+Oaa67Bf//3f3u0lnlMTAwyMjJQX19Pa58TEoTCPpkPDAygqakJXC4XKSkp7PFIGfy0\nx+FwUFxcjMrKSgDArbfeCqFQiN27d3v0+Pj4eMTFxaGhoWHUZQEIIf4X9sm8vr4eSqXSob/carWi\nv78/LGc3ujNp0iRotVpcvnwZPB4PjzzyCMrLy3HgwAGPHp+RkQGGYXD58uUJjpQQMhZhn8xtFRzR\n0dHsMrcGgwFisRg8Hi/A0fkfj8fD/Pnzcfz4cVitVkgkEmzcuBH79u3Dvn373D6ew+EgNzcXGo0G\nPT09foiYEOKJsE/mFy5cQHx8vEddLJGyZkdubi6EQiG7Bn1ycjK2bt2Kr7/+Gl999ZXbx/P5fOTl\n5aGlpQX9/f0THa7PRcr7TCJLWCdzi8WCCxcuICoqyiGZj1RfHiklaxwOBwsXLsQPP/zA7j6UkJCA\nZ599Fn/7299gMBjc3kMsFrNrvoRahUukvM8ksoR1Mm9paYFMJoNer2f7y4eGhmAymRATExPg6AIr\nMTERWVlZ+Omnn9hjKSkp+Jd/+Rd8+eWXHt1DoVAgLS0NtbW11OVCSICFdTKvqanBpEmTYDQaoVAo\nAPyyEQWXG9bfukdmzZqFqqoqmEwm9tgtt9yCAwcOeLzbUHx8PPLz89HW1hYxe4gSEozCOqOdP38e\nqampSE5OBofDARCZ9eUjkUqlyM3NZUsVgeHW+axZszyubgGGu1wKCwthNptRU1Mz4sbRhJCJE7bJ\nnGEYVFdXIzY2lu1iASKzvnw0V111lcvW+Zdffgm1Wu3xfWyzRBUKBS5cuODRzFJCiO+EbTJvaGiA\nUCjEwMAAO/g5ODgIi8XCliheKRIHxmyt89OnT7MTgdLS0rBy5Uq89NJLY2plczgcJCcns5Uuwbra\nYiS+zyT8cZgJmsrH4XACOkvw448/htFoBMMwuOuuuyAQCNDe3o6BgQFkZ2e7fEygYw6U/v5+7N+/\nH6mpqVi0aBG4XC4YhsG2bdtgsVjw6KOPst1UnrJYLGhpaYFer0dOTk5QTdCK1PeZhAZvfz7DtmV+\n6tQppKSkICsrCwKBABaLBV1dXQ5dLmSYWCzG6tWrodPpcODAAVgsFnA4HDz44IPo7u7Gq6++OuZq\nFR6Ph+zsbLZ8sbW1lQZHCZlAYZnM29raoNfr0d/fj4KCAgBAd3c3JBLJiF0skU4kEuGGG26AxWJB\ndXU1AEAoFOK5555DcnIyHnvsMXzxxRdjHtyUyWQoKiqCyWTChQsXoNPpqFVMyAQIy2ReXl6O4uJi\nGAwGZGRkwGq1orOz02HiEHHG4/GwcOFCnDlzhh0QjY6Oxrp16/D888+jvr4eDz30EPbv3z+mvnCB\nQIC8vDykpKSgubkZNTU10Gq1lNQJ8aGw7DP/4x//iBkzZiAjIwOLFi1CV1cXdDodJk2aNOrjqC91\n2NGjR8HlcrFo0SKnc42Njdi+fTtKSkpw2223jfneDMNAo9Ggvb0dXC4XqampiIuLG3Of/HjQ+0yC\nGfWZ/5NGo0FraysGBwcxefJkWK1WdHR0IDU11e1jac2OYXPmzMHFixddliZmZ2fj8ccfx4EDB9DY\n2Djme3M4HCgUChQVFSElJQVtbW04f/48urq6/LZOOr3PJByFXcv8L3/5C0wmE9LS0vBv//ZvUKlU\n6O3tRX5+vt9jCWU1NTU4ceIEZs2ahWnTpjnNmD18+DC++uorvPjii+BwOOjv7/eqfp9hGOj1eqhU\nKuh0OkilUiQkJEAikfi1tU5IsPA2d4ZVMj916hR27dqFefPmYcmSJcjMzERVVVXQlcaFCo1GgyNH\njqCvrw8xMTEQCATIyclBQUEBeDwenn/+eXR0dKCnpwdWqxW//e1vsXTpUq+fb2hoCD09PVCr1Rga\nGoJcLodMJkNMTAwldhIxIj6Z9/T04IknnsDKlSsRFRWF5cuXQ6VSoaenh61oIWPHMAxUKhXMZjOM\nRiMuXLiA7u5uTJ06Fbm5uejq6kJmZiY0Gg22bNmCW2+9FcuWLWMfb7FY0NzcjMzMzDGtH9/f3w+N\nRgOtVguz2QypVAqZTAapVBqR69CTyBHxyfyFF15Aeno6LBYLbr31VsTExKCqqgpZWVk0fd/Henp6\ncO7cOdTV1SE3NxcLFiyASCRCe3s7Nm/ejMTERCQnJ4NhGFRUVIDL5SI3NxePPfaYV6WhAwMD0Gq1\n0Gq1MBgMEIlE7GYjtg+BQECtdxIWIjqZV1RUYOfOnViyZAmmTJmCGTNmQK1WQ6VSYfLkyX6JIRKZ\nTCaUl5ejtbUVy5cvh0KhgMFgwKVLl9Dd3Q2z2YyrrroKCoUCO3bsQH19PZ544olxTdyyWq0wGo1O\nHwzDQCwWIyYmBjExMRCLxRAKhT78bgnxj4hN5mazGY899hiKi4uRlZWFa6+9FkajEfX19cjKyhrT\nCombNm2idTu8UFNTg++//x4zZsxAUVERoqKinK5hGAZ79+7F3//+dyQlJWHOnDlYvny5z1awNJvN\n6OvrQ39/P/uZw+EgNjYWcXFxiIuLA5/PB0DvMwluEZvM9+/fj6NHj2Lu3LlYtWoVurq6oFarkZ6e\njoSEhDHdi+qPvafRaPDjjz+isbER2dnZSElJQXJyMhQKhUP3h8ViQU1NDY4ePYrvv/8e1113HRYv\nXoz09HSf9oUzDIPBwUHodDpotVro9XpER0cjLi4OaWlpDrsj2cdne/8ZhmE/hEIh+Hw+deMQv/B7\nMj948CA2bNgAi8WC++67D0899ZRPAvKU2WzGkSNH8N5772HhwoW46aab0N3djZiYGCiVSggEAo/u\nU1ZWhtLSUr/EPB72cQazgwcPIiMjA52dnejo6IDVakVubi6ioqKg1+sxMDAAsViM2NhYSKVSfPPN\nNzh37hzUajUyMjKgVCqRnp6OwcFBaLVaJCcn49prrx13NZLVaoVer4dWq0VWVhbOnTvHnrN/z20J\nm8PhsP8eHBwEAERFRSE6OhpRUVEQiUTsNfaP4fP5EAqFPtn8JFTec4rTt7zNQ3xvnsxiseDhhx/G\noUOHkJ6ejjlz5mDVqlUoLCz05nZjotfrUVZWhn/84x8Qi8WYN28eSktL0d3djczMTMTFxY3pfqHy\nBodKnCdPnsSmTZswbdo0MAyDnp4eNDQ0wGw2IzExESKRCP39/dDpdKisrER6ejoWLFgAq9UKjUaD\ngYEBdHR0QCQSITMzE3V1dXjooYcwa9Ys8Pl8mEwm6PV66HQ6REdHY9myZViwYAGEQiFbcWNbaiAu\nLo5Nqlwul+1uAYBp06Z5/D0xDMNuN2gymWA0GqHX6x1a7/bXDQ4Ogs/nQyQSQSgUQiQSgc/nsx88\nHo/9PJrDhw87zMLl8XhB+b+DUPnZDJU4veVVMi8vL8ekSZPYpWRvv/127N2716fJnGEY9gdXq9Xi\nzJkzKC8vx7lz56BUKlFQUIC0tDQUFBRAIpEgNTWVStaCDIfDQXx8POLj412enzdvHi5evMgOlhoM\nBqhUKmRmZkIoFKKvrw9ZWVlQKpW4fPkygOG1YhQKBWJjY2E0GvHll19ix44dbCIViUQQCASwWq0w\nmUxQKBRISkpCeno60tPT2ViOHz+O9vZ2MAwDuVyO+Ph4KJVKSCQSCAQCcLlcWK1W9PX1sXX2MTEx\nHlVG2bp4BgcHMTAwgMHBQZhMJgwNDWFoaAgWi4X9PJquri78/PPP7NcWiwVcLhc8Hs/hg8PhgMvl\nsp9t/7Y/PtK/bZ9t79eVn13978PVHzCj0Tjitbbr7D/b4rSPiYyPV8n88uXLUCqV7NcZGRk4deqU\nVwFUVlay08IZhoHFYoHVasWJEyfQ3NzM/sDFx8cjOTkZy5YtQ0ZGBnJycpCens7+d5eEHj6fj8LC\nQodGQH9/PxobG2G1WpGVlQU+nw+j0YhJkyaBy+VCIBCAYRiYzWaYzWbk5+dDrVbDarWCw+FgaGgI\nAwMDMJvNAH4pa2xsbMS5c+fYBcQ++OADdi9Yo9GI/v5+GAwGtkFgtVoxNDQEPp8PgUAAs9mMoaEh\n8Hg8NgExDAOr1cp2rfB4PFgsFjAMwyYq271sx2w/q1arlf1+bGWVtkRnu76iogLvvPMOGwuHw2ET\nuO25AbDPY7FY2IRvG+w1m81sjLbEb3se+1iujNF2f9vvn42tDNT2OIvFghMnTrBjEDfffDP7B8/+\neVz9MbB/3iuvufIPi81Iv+vucgDDMOjs7ERlZaVDF8aVf9RGel777+fKz64e42rsxf57u/IPmC9y\nmFd95nv27MHBgwexY8cOAMDu3btx6tQpvP766z4NjhBCIpHf+szT09PR0tLCft3S0oKMjIxxB0MI\nIcQ7Xg25z549GxcvXkRjYyMGBwfxt7/9DatWrfJ1bIQQQjzkVcucz+fjjTfewPXXXw+LxYLf/OY3\nfqlkIYQQ4prXxbA33HADampq8MYbb2DXrl3Iz8/HSy+95PLasrIylJSUYNq0aQErDTp48CCmTJky\nYpyvvPIKSkpKUFJSgunTp4PP56O3tzeoYlSpVFi+fDlmzpyJadOm4b333vNrfDbu4tRoNLjppptQ\nXFyMuXPnoqqqyu8x3nvvvUhOTsb06dNHvOY//uM/kJ+fj+LiYvz4449+jO4X7uK8cOEC5s+fj6io\nKPzXf/2Xn6P7hbs4P/zwQxQXF2PGjBlYuHAhzp496+cIh7mLc+/evSguLkZJSQlmzZqFw4cP+znC\nYZ78fALADz/8AD6fjy+++ML9TZlxGBoaYvLy8piGhgZmcHCQKS4uZqqrqx2u0Wg0TFFREdPS0sIw\nDMN0d3eP5yknLE57+/fvZ371q1/5MULPYty4cSPz9NNPMwwz/DoqFArGbDYHXZyPP/44s2XLFoZh\nGObChQt+fy0ZhmGOHj3KVFRUMNOmTXN5/ssvv2RuuOEGhmEY5uTJk8zcuXP9GR7LXZxdXV3MDz/8\nwDzzzDPMK6+84ufofuEuzhMnTjC9vb0MwzDMgQMHgvb1NBgM7L/Pnj3L5OXl+Ss0B+7iZJjh37Wl\nS5cy//qv/8p8/vnnbu85rmlq9vXmAoGArTe399FHH+GWW25hB0jHOsXeFzyJ095HH32EtWvX+jFC\nz2JMTU2FTqcDAOh0OsTHx7MlaMEU5/nz59l1zSdPnozGxkZ0d3f7Nc7FixdDLpePeH7fvn24++67\nAQBz585Fb28vOjs7/RUey12ciYmJmD17tsczmieKuzjnz5/PTsiaO3cuWltb/RWaA3dxxsTEsP82\nGAwByUeA+zgB4PXXX8eaNWuQmJjo0T3Hlcxd1ZvbJnfYXLx4ET09PVi6dClmz56NDz74YDxP6RVP\n4rTp7+/H//3f/+GWW27xV3gAPIvx/vvvR1VVFdLS0lBcXIzXXnvNrzECnsVZXFzM/rewvLwcTU1N\nAfvlHomr7yPYYgxVf/3rX7FixYpAhzGi//3f/0VhYSFuuOEGbNu2LdDhuHT58mXs3bsXDz74IADP\nSr3H1azz5AnMZjMqKirwzTffoL+/H/Pnz8e8efP8uo3bWGre9+/fj0WLFkEmk01gRM48ifH555/H\nzJkzUVZWhvr6eixbtgyVlZV+Xa/dkziffvppPPLII+z4Q0lJSVDOzmWuKJ+luRHj9+2332Lnzp04\nfvx4oEMZ0erVq7F69WocO3YMd955J2pqagIdkpMNGzawWzIydpOORjOuZO5JvblSqURCQgK7icCS\nJUtQWVnp12TuSZw2n3zyid+7WADPYjxx4gSeeeYZAEBeXh5ycnJQU1OD2bNnB1WcEokEO3fuZL/O\nyclBbm6u32L0xJXfR2trK9LT0wMYUeg7e/Ys7r//fhw8eNBtF0IwWLx4MYaGhqBWq0dcciJQzpw5\ng9tvvx3AcOHDgQMHIBAIRi0BH1c3iyf15jfeeCO+++47WCwW9Pf349SpUygqKhrP005InMDwGjBH\njx7FjTfe6Nf4PI1xypQpOHToEACgs7MTNTU1fk+SnsSp1WrZlQZ37NiBq6++Ouj2YF21ahXef/99\nAMOLg8lksnFtmjHRPGmZBVJzczNuvvlm7N69G5MmTQp0OCOqr69nX8uKigoACLpEDgCXLl1CQ0MD\nGhoasGbNGrz55ptu5/KMq2U+Ur359u3bAQAPPPAApkyZguXLl2PGjBngcrm4//77/Z7MPYkTGO5L\nu/76673a2swfMf7xj3/E+vXrUVxcDKvViv/8z/+EQqEIujirq6txzz33gMPhYNq0afjrX//q1xgB\nYO3atThy5AhUKhWUSiU2b97MrtfywAMPYMWKFfjqq68wadIkxMTE4N133/V7jJ7E2dHRgTlz5kCn\n04HL5eK1115DdXW13/84uotzy5Yt0Gg0bB+vQCBAeXm5X2P0JM49e/bg/fffh0AgQGxsLD755BO/\nx+hJnN6YsM0pCCGE+M/4V9AnhBAScJTMCSEkDFAyJ4SQMEDJnBBCwgAlc0IICQOUzAkhJAz8fxIF\n68DHQL2FAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x108a66710>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for this in [rrmse_1k, rrmse_2k, rrmse_4k]:\n",
      "    isfin = this[np.isfinite(this)]\n",
      "    print \"The proportion of voxels with rRMSE<1.0 is %s\"%(100 * len(np.where(isfin<1)[0])/float(len(isfin)))\n",
      "        "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The proportion of voxels with rRMSE<1.0 is 98.1198808637\n",
        "The proportion of voxels with rRMSE<1.0 is 99.9084760486\n",
        "The proportion of voxels with rRMSE<1.0 is 99.8665922065\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}