{
 "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: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvals\n",
        "Loading from file: /home/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%******************]  68693 of 68694 complete "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/home/arokem/usr/local/lib/python2.7/site-packages/numexpr/necompiler.py:716: DeprecationWarning: using `oa_ndim == 0` when `op_axes` is NULL is deprecated. Use `oa_ndim == -1` or the MultiNew iterator for NumPy <1.8 compatibility\n",
        "  return compiled_ex(*arguments, **kwargs)\n",
        "/home/arokem/usr/local/lib/python2.7/site-packages/numexpr/necompiler.py:716: DeprecationWarning: using `oa_ndim == 0` when `op_axes` is NULL is deprecated. Use `oa_ndim == -1` or the MultiNew iterator for NumPy <1.8 compatibility\n",
        "  return compiled_ex(*arguments, **kwargs)\n",
        "/home/arokem/usr/local/lib/python2.7/site-packages/numexpr/necompiler.py:716: DeprecationWarning: using `oa_ndim == 0` when `op_axes` is NULL is deprecated. Use `oa_ndim == -1` or the MultiNew iterator for NumPy <1.8 compatibility\n",
        "  return compiled_ex(*arguments, **kwargs)\n"
       ]
      },
      {
       "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+/AAAIABJREFUeJzt3Xt4U3WaB/DvyaVJ2zRN07Tp/X6/cpWrUHBQcQVvzOww\nrCiOjrqPrrrrii4+CjoqLrouOI4XdlAZFJ3RUcEVxkEsOhZEKRQpUEppS1uaNm2StmlubXr2j5Kz\nDU3bNElzfT/P04f2nJOTl6R9++v7uzEsy7IghBAS0Hi+DoAQQoj7KJkTQkgQoGROCCFBgJI5IYQE\nAUrmhBASBCiZE0JIEBg3md91111QKpUoLS3ljv37v/87CgsLUV5ejltvvRU9PT1THiQhhJDxjZvM\n161bh/3799sdu/baa1FbW4uamhrk5eXhhRdemNIACSGETGzcZH711VcjJibG7tiyZcvA4w0/bM6c\nOWhtbZ266AghhDhF4M6Dd+zYgdWrV486zjCMO7clhJCQ5eqkfJc7QJ977jmEhYXhV7/61ZgB+fvH\n008/7fBzf/vw59gCLc7Fixf7PIZgeS0pTs9/uMOlZP7OO+/giy++wHvvvefWk/uTTZs2+ToE4gWH\nDh3ydQiETIlJl1n279+PLVu24NChQxCLxVMREyGEkEkat2W+evVqzJ8/H3V1dUhNTcWOHTvw4IMP\nQq/XY9myZZg+fTr++Z//2VuxelxFRYWvQ3AKxRl6AuW1pDj9B8O6W6hxdFOGcbv+422BGDOZPHqf\niT9z5/vTrdEshJDgI5fLodVqfR1GUIuJiYFGo/HoPSmZX/b000/7OgTiBfQ+T0yr1dJfL1NsKoZv\nU5mFEGKHfn6n3livsTuvPS20RQghQYCSOSGEBAFK5oSQgJCRkYGvvvrK12H4LUrmhJCAwDCMWx2H\nAwMDWLVqFTIzM8Hj8RzOBl6/fj0UCgUUCgUef/xxu3NNTU1YsmQJIiMjUVhYOOoXy/vvv4/09HRI\nJBLccsstXh8RRMn8so0bN/o6BOIF9D6HtkWLFmHXrl1ISEgY9YvhzTffxGeffYaTJ0/i5MmT2Lt3\nL958803u/OrVqzFz5kxoNBo899xzWLVqFbq6ugAAtbW1uO+++/Dee++ho6MDERER3p9QyU6BKbrt\nlArEmMnk0fs8MX99jTIyMtgXXniBLSoqYmNiYth169axJpPJpXulpKSwhw4dsjs2b948dvv27dzX\nO3bsYOfOncuyLMvW1dWxIpGI1ev13PlFixaxb7zxBsuyLPvEE0+wa9as4c41NDSwYWFhdtePNNZr\n7M5rTy1zQkhAYFkW77//Pr788ks0NDTg3Llz+O1vf4uWlhbIZDLExMQ4/Pjggw+cuv/p06dRXl7O\nfV1WVoba2loAwy3vrKwsREZGcufLy8vtzo98bFZWFkQiEc6dO+eJ/7pTaNIQIWRSjh075pH7zJw5\nc1LXMwyDBx54AMnJyQCADRs24MEHH8Szzz4LnU7ndjx6vR7R0dHc11KpFHq93uE5AIiKikJ7ezsA\noL+/f9R5qVSKvr4+t+NyFiVzQsikTDYJe1Jqair3eVpaGi5duuSxe0skEvT29nJf9/T0QCKRODxn\nOx8VFcWdv3I/5JHnvYHKLISQgHHx4kW7z5OTk9HS0gKJRIKoqCiHH7t373bq3sXFxThx4gT3dU1N\nDUpKSrhzFy5c4FrqtvPFxcXc+ZqaGu5cQ0MDLBYL8vLy3Pr/TorL1fZxTNFtp9TTTz/t6xCIF9D7\nPDF//flNT09nS0tL2dbWVra7u5tdsGABu2HDhkndw2QysUajkU1JSWG//PJL1mg0cufeeOMNtrCw\nkG1ra2NbW1vZoqIi9s033+TOz507l3300UdZo9HIfvzxx6xMJmO7urpYlmXZ2tpaViqVst9++y2r\n1+vZ1atXs6tXrx4zjrFeY3dee0rmhBA7/vrzm5GRwW7evJktKipiZTIZe+edd9olY2ekp6ezDMOw\nPB6P+7e5uZk7/9hjj7FyuZyVy+Xs+vXr7R7b1NTEVlRUsOHh4WxBQQH71Vdf2Z1///332bS0NDYy\nMpK9+eabWa1WO2YcU5HMaaEtQogd+vmderTQFiGEEIcomRNCSBCgZE4IIUGAkvlltGZHaKD3mQQr\n6gC9LBBjJpNH7/PE6DWaetQBSgghxCFK5oQQEgQomRNCSBCgZE4ICQi0bdz4KJlf9vTTT/s6BOIF\n9D4HLne3jTty5AiWLVuG2NhYxMfH4xe/+AVUKpXdNbRtXBCgIWuhgd7n0KXT6XDfffehubkZzc3N\niIqKwrp167jztG2cA1N0W0KIF/jrz68nt41jWZY9duwYGxUVxX0d1NvG3XXXXVAqlSgtLeWOaTQa\nLFu2DHl5ebj22ms9ssMHIYRMhPXwtnHffPMNt145EOTbxq1btw4PPvgg1q5dyx3bvHkzli1bhsce\newwvvvgiNm/ejM2bN095oIFgcHAQLMtCKBT6OhRCpszrr7/ukfvcf//9k7rek9vGnTx5Es8++yz2\n7NnDHQvqbeOuvvpqNDU12R3bs2cPDh06BAC44447UFFRQcn8ssOHD6OxsRGLFi1CRkaGr8MhZEpM\nNgl7kie2jTt//jxuuOEGbNu2DQsWLOCOB/q2cZPeA7SjowNKpRIAoFQq0dHR4fC6kR1NFRUVqKio\ncClAb9m4caNbnWMsy6KxsRGzZs3Cd999h/Pnz2PBggUIDw/3XJDEbe6+z8S3xto2rrCwcMyRLm+9\n9RZWr14NAGhubsayZcvw1FNPYc2aNXbX2baNmzVrFoCxt42zJfiamhrcfvvt3HlXto2rrKxEZWXl\nJF6BcUxUVG9sbGRLSkq4r2Uymd35mJgYjxbxfcXdmFUqFbt7926WZVnWYrGwhw4dYj/55BNPhEY8\nKBC/N73NX18jd7eNa21tZbOystiXXnrJ4fmg3zbuymSen5/Ptre3syzLspcuXWLz8/M9GpCvuBvz\n4cOH2SNHjnBfW61WdteuXWxbW5u7oREPCsTvTW/z19fI3W3jNm7cyDIMw0okEu5j5GgWlg3ybeOa\nmpqwYsUK/PTTTwCAxx57DLGxsVi/fj02b94MnU43qmYeiKuuuRvz7t27sXTpUq4EBQz3jl+4cAE3\n3nijJ0IkHhCI35veRq/R1PP6qomrV6/G/PnzUVdXh9TUVLz99tt4/PHH8be//Q15eXk4ePDgqFlS\noUir1WJgYADx8fF2x/Pz86HRaKBWq30UGSEkVNB65pe5E3N1dTX0ej0WLVo06lxNTQ1UKhWuu+46\nd0MkHhCI35veRq/R1KP1zKeQO2t2NDY2Iisry+G5oqIitLe3Q6PRuHx/4jm0NgsJVtQyd5Ner8ef\n/vQn3HHHHeDz+Q6vOXbsGHQ6Ha655hovR0fI5IXSz6+vUMvcDzU1NSEtLW3MRA4AhYWFaGpqgtVq\n9WJkhJBQQsncTa2trUhPTx/3moiICMhkMm7qLyGEeBolczdpNBooFIoJr0tPT7ebvUYIIZ5EydwN\ng4OD0Ov1kEqlE16blpaG5uZmL0RFCAlFlMwvc2W9Dq1Wi+jo6HHr5TZxcXEwm82jFush3kXrsgQu\n2jZufJTML9u0adOkH6PRaCCXy526lmEYpKWlUanFx1x5n4l/cHfbuJGeeeYZ8Hg8HDx40O44bRsX\noq5M5hMNKaJSCyG+19DQgI8++ghJSUl2xwN92zhK5m7QarV2yXz37t14/fXXx0zqqampaG9vx8DA\ngLdCJCSoHD16FMXFxZDL5bjrrrtgNpsnfY8HHngAL7744qhNZN599108+uijSEpKQlJSEh599FG8\n8847AIBz587h+PHj2LRpE0QiEW699VaUlZXh448/BgC89957WLlyJRYuXIjIyEg8++yz+Mtf/oL+\n/n63/8/OomTuhpEt84GBAfztb39DbW0tvvzyS4fXi0QiKBQKlxbUJyTUsR7YNu7Pf/4zxGIxli9f\nPur+Qb1tHBmbxWKB0WjkdhI5evQo0tLScO+992LDhg3IzMx0uDC9bYjiRGPTCfFXt912m0fuY2vV\nOsvdbeP6+vqwYcMGHDhwwOH5oN42LpRMds0OrVYLmUwGHm/4j5uvvvoK11xzDZKSknD//ffjlVde\nwbZt20b9KZeWloZ9+/Zh4cKFHuvMIc6jtVncN9kk7EnubBu3ceNG3H777UhLS+OOjSyJBvq2cVRm\nuWyyQ9ZGlljUajUaGhowZ84cAMBVV12FlJQUh8Oo5HI5hoaGJr0BLfEMGpoY2MbaNk4ikSAqKsrh\nx+7duwEABw8exLZt25CYmIjExES0tLTgF7/4BbZs2QLg/7eNsxlr27iR54uLi7nzrmwb51Eub2sx\njim6rV/5+9//zlZXV7Msy7Iffvgh+9Zbb9mdr6+vZ++++27WbDaPemxlZSV74sQJr8RJyGT568+v\nu9vGdXd3sx0dHWxHRwerUqnY1NRU9qOPPmL7+/tZlg38beOoZe4i20gWlmXx9ddfj1oRMScnB1lZ\nWQ7rczTenJDJYxgGa9aswbXXXovs7Gzk5ubiySefdPrxcrkc8fHxiI+Ph1KpBJ/PR0xMDCIiIgAA\n9957L1asWIHS0lKUlZVhxYoV+M1vfsM9/oMPPsCPP/4IuVyODRs24OOPP0ZsbCyA4aWu33jjDaxZ\nswZKpRJGoxG///3vPfsCTICWwHXRzp07ccstt6C9vR2/+93vsHXr1lE18AsXLuD555/Ha6+9BpFI\nxB0fGBjAu+++izvuuGNUTZ0QXwuFn19foyVw/YTZbIbFYoFEIkFVVRXmz5/vsDMzKysL+fn52Lt3\nr91xoVAIpVKJ1tZWb4VMCAlylMwvm0zHmEajQUxMDFiWxeHDhzF//vwxr127di327t3LzRSzoVKL\nb1AHKAlWlMwvm8yaHbaRLPX19RCLxXZDna6kVCqxfPly7Ny50+64LZnTn7PeRWuzkGBFydwFts5P\nW4llIrfccgvq6upw6tQp7phMJgPDMLQ3KPE7MTEx3KJW9DE1HzExMR5/3yiZu0Cj0UAmk+Hw4cOY\nN2/ehNeLRCKsWbMGH330EXeMYRjasIL4JY1GA5Zl6WMKP6aiEUfJ3AUajQZarRYRERHjllhGmjNn\nDhoaGuxmkVHdnBDiKZTMJ8loNMJqteL06dOYPXu2048TiUQoKyvDjz/+yB1LSkqCWq2mVRQJIW6j\nZH6Zs2t22OrlJ06cwLRp0yb1HHPmzMH333/PfS0UChEVFUW7D3kRrc1CghUl88ucHbKm0WggFovR\n3t4+6XUXZs6cidraWhiNRu4YJXPvoqGJJFhRMp8kjUYDjUaDkpKSSc/ejIyMRH5+Po4fP84d8/Yy\nmYSQ4ETJfJK0Wi1aWlomXWKxmTt3rl2phVrmhBBPcDmZv/DCCyguLkZpaSl+9atfubR9U6BhWRbd\n3d04d+6cy8l89uzZOH78ONfpScmcEOIJLiXzpqYmbN++HdXV1fjpp59gtVrttmYKVkajET09PYiI\niIBSqXTpHjKZDAkJCTh//jwAKrMQQjzDpWQulUohFAphMBgwODgIg8HAbeUUqJzpGNNoNOjv73e5\nVW5TXFzM7R0YFRWFvr4+mtbvJdQBSoKVS9vGyeVy/Nu//RvS0tIQHh6O6667Dj/72c/srhn5Q1NR\nUYGKigp34pxymzZtmvAHXavVQq1W4/rrr3fruYqLi/H5559j1apVEIlE4PF4MJlMCA8Pd+u+ZGLO\nvM+EeEtlZSUqKys9ci+XknlDQwP++7//G01NTYiOjsbPf/5zvPfee1izZg13TTD+wHR0dKCzs5Pb\nKmoiFosFzc3NYFkWPB4PDMOAx+NBIpGgvr4eAwMD3Fjzvr4+SuaEhJgrG7ruLATnUpnlxx9/xPz5\n8xEbGwuBQIBbb70VVVVVLgcRKOrq6pCamup00u3q6oJAIEBCQgIUCgViYmIglUoRGRmJ2NhYu7o5\ndYISQtzhUjIvKCjAkSNHYDQawbIsDhw4gKKiIk/H5ldYlkVjYyOmT5/u9PUajQbx8fGQSqWQyWSQ\ny+WIjY1FcnIy0tLSuA1gbS1zQghxlUvJvLy8HGvXrsWsWbNQVlYGAHZ75QUjg8GArq4uzJw506nr\n+/v7AYDbX3AkHo+H0tJSLplTy5wQ4i6Xx5k/9thjqK2txU8//YR333034PeynGjNjqamJpjNZmRn\nZzt1v+7ubsTGxjrcTg4YHm9uuye1zL2H1mYhwYpmgF42UYdtdXU10tLSwOfzJ7zX0NAQtyDXWORy\nOeLi4lBdXU0tcy8Kxo55QgBK5k47c+YMCgsLnbq2p6cH4eHhEIlE415XXFyMEydOQCKRQK/X01hz\nQojLKJk7gWVZNDc3Y8aMGU5dr9FoEBsbO+F1ZWVlaG5uhtlsRlhYGFdnJ4SQyaJk7oTOzk4AQE5O\nzoTXDg4Ooq+vz6k9/goLC9Ha2gqVSkXT+gkhbqFk7oT29naEh4c7Nb5co9FAKpU6VVuXyWSQSqW4\ncOECJBIJ1c0JIS6jZH7ZeB1jbW1tiIqKGnNkykjOllhsCgsL0d3dDT6fTy1zL6AOUBKsKJlfNt40\n2vb2dkRHR094D5PJBLPZDKlU6vTzFhQUoL29HVarlVrmXuDOdGlC/Bklcyd0dnY61drWaDSQy+VO\nteBt8vPz0dDQAIlEAq1W606YhJAQRsncCV1dXVAoFONeY9u4YjIlFgBISUlBX18foqOjqWVOCHEZ\nJXMn6HQ6JCQkjHtNf38/eDzepFc+5PF4yM/Ph06ng9lsttvsmRBCnEXJfAJDQ0Po6+tDUlLSuNdN\nNH1/PPn5+aivr4dIJEJ7e7uroRJCQhgl88vGWrNDp9NBKBROOG68p6cHMpnMpecuKCjA2bNnER0d\nDZVK5dI9iHNobRYSrCiZXzbWkLWuri6IxWJIJJIxHzswMIChoaEJp++PJTs7G83NzVAqlejq6oLV\nanXpPmRiNDSRBCtK5hNQqVQQiUTj1sL7+/sRGRnpUokFAMLDw6FQKMDn82E0GqkjlBAyaZTMJ3Dp\n0iVIpdJxE7XBYHC4bvlkZGRkQK/Xo7e3Fzqdzq17EUJCDyXzCahUqglr4baWuTsyMzPR2tqKiIgI\ntLe30wqKhJBJoWQ+AbVaPe7YcZZl0d/f73bLPDMzE42NjUhISIDBYKAVFAkhk0LJ/LKxOsa6u7sR\nHx8/5uMsFgt4PB7CwsLcev7MzEw0NTVBqVTCZDKhp6fHrfsRx6gDlAQrSuaXOVqzg2VZ9PT0jDth\nyBOtcmB4BUWBQAChUIienh6qm08RWpuFBCtK5uMwGAxgWXbcqfwGg8HterlNVlYWtFotjEYjjEYj\nzGazR+5LCAl+lMzHoVarERERMe4Yc0+1zIHhUktzczPi4+PBsiy1zgkhTqNkPg61Wj3uhCGWZT3a\nMs/IyEBjYyPVzQkhk0bJfBydnZ0QCoVjThgymUwQCAQQCAQeeT7biBalUgmtVov+/n6aDUoIcQol\n88scrdnR1tY27hZwnmyVA4BSqeTKNmq1GpGRkdQ69zBam4UEK0rmlzkastbe3g65XD7mYzwxWWgk\nHo+HjIwMqFQqhIeHg8fjUTL3MBqaSIIVJfNxdHZ2jjvG3BPT+K+UlZWFxsZGxMfHc3Vzmg1KCJkI\nJfMx2HYOGmuM+dDQEIxGo8eTeWZmJi5cuIDY2Fj09vYiLCwMer3eo89BCAk+lMzHoNPpIBAIxiyz\nmEwmhIWFjVlPd5WtEzQmJgYajQYymYxKLYSQCbmczHU6HVatWoXCwkIUFRXhyJEjnozL5zo6OhAV\nFTXmsERP18ttUlJSuPHtWq0W0dHRlMwJIRNyOZk/9NBDuOGGG3DmzBmcPHkShYWFnozL667sGFOp\nVIiMjBwzYXtystBIAoEAaWlp6O7u5oY+Wq1WmEwmjz9XKKIOUBKsXErmPT09+Pbbb3HXXXcBGE5A\n0dHRHg3M265cs0OlUiEsLGzMZO7pYYkj2UotMpkMOp2OWuceRGuzkGDl0myXxsZGxMXFYd26daip\nqcHMmTOxdetWu5bqyBZQRUUFKioq3I3Vq1QqFYRCocOEPTQ0BJPJNO7uQ+7IyspCfX09CgoKoNVq\nkZiYCJVKBaVSOSXPRwjxjcrKSlRWVnrkXgzrwri3H3/8EfPmzUNVVRVmz56Nhx9+GFKpFM8888zw\nTRkm4IbTXRnz+vXrERcXh0cffXTUtXq9Hi0tLVNWWjp//jxef/11rFmzBiaTCXPnzkVNTQ1KS0s9\nNts0VAXi9yYJHe58f7pUZklJSUFKSgpmz54NAFi1ahWqq6tdCsBfdXR0jDnGfKrq5TZpaWncdnVa\nrRY8Hg9SqZRKLYSQMbmUzBMSEpCamopz584BAA4cOIDi4mKPBuZLRqMRJpNpzKVvp7JeDgBhYWFI\nTExEf38/NBoNACA6OppWUSSEjMnlv9lfffVVrFmzBhaLBdnZ2Xj77bc9GZfXjVyzo6OjAzKZDFKp\n1OG1/f39425Y4QmZmZno6OiAyWSCxWJBdHQ0WltbMTQ0BB6Ppge4itZmIcHK5WReXl6OH374wZOx\n+NTIDtvxxpgPDg5iYGAAYrF4SuOxTeuXy+XQarVQKpUQi8Xo6+sL+JFDvkRDE0mwoiaeA7Yx5o6S\nuW09FoZhpjSGoqIi/PTTT4iJiYFWqwUAmg1KCBkTJXMHVCoVRCIRoqKiRp2bisW1HMnIyOAmCtnq\n5rZx5zQagxByJUrmDnR0dIDP5zvs5JzK8eUjMQyDGTNm4NKlS1zLXCQSgcfjwWg0TvnzE0ICCyVz\nB1QqFQQCgcMWuNFonPJ6uc2MGTNw/vx5LpkzDMO1zgkhZCRK5pfZOsbMZjM0Gg0UCsWoUSMsy3qt\nZQ4AZWVlaGhoQF9fHwYGBgCApva7iTpASbCiZH6Zbc2OtrY2KBQKh8MSLRYL+Hy+x5e9HUtERARy\ncnK4XzAAIJFIYDabYbFYvBJDsKG1WUiwomR+hZaWFigUCocjWbzZKreZMWMGNBqNXamFWueEkCtR\nMr9CS0sLZDKZw5Es3qyX28yYMQMtLS3o7u7mjlHdnBByJUrmV2htbYVEIvGblnlKSgoAoLm5mTsm\nlUqh1+thtVq9GgshxH9RMr9CS0sLxGKxw2Tui5Y5wzDIyclBU1MTd4zP50MikaC3t9ersRBC/Bcl\n88uefvpprqORx+ONSubeHskyUmFhIbq7u2E2m7ljVDd3Da3NQoIVJfPLNm7ciLa2NiQkJMBgMIxK\n5t4eyTJSbm4u9Hq9Xd3ctiQuzQadHBqaSIIVJfMRWlpakJycDIZhIBKJ7M75qlUOANnZ2dDpdFCr\n1dwxsVgMPp9Ps0EJIQAomdtpaWlBXFyc39TLbSIiIhATE4P6+nq741RqIYTYUDIfoaWlBXK53G9G\nsoyUmZmJCxcu2B2jZE4IsaFkPsLFixfHXMfcly1zYHhJ3La2NrsauUQigdFoxODgoM/iIoT4B0rm\nlz355JPQ6XQICwvzq5EsNsXFxejp6UF/fz93jMfjISoqilrnk0AdoCRYUTK/7LnnnkNSUhL6+/tH\nzf705UgWm7S0NBiNRly6dMnueHR0NI03nwRam4UEK0rmIyQlJUGv149qmfu6VQ4AAoEAcXFxOH36\ntN1xW92chigSEtoomY8QHx/vMJn7ul5uk5GRMWpES1hYGMLCwuzKL4SQ0EPJfASFQgGDwTBqhyGj\n0ejzljkAlJSU2K3RYkOjWgghlMxHkEqlCA8PH1UbN5lMftEynzVrFtRqNbdRhY1tNighJHRRMr/s\nqquuQnh4uF+OZLFRKBQQi8Wj6uYSiQQWi4U2rHACrc1CghUlcwwn7KysLAiFQkRHR9ud84eRLCMl\nJyfjxIkTdscYhoFUKqVRLU6goYkkWFEyB6DT6SAWi2E0GkdtF+cv9XKb3Nxc1NXVjTpOdXNCQhsl\ncwBqtRrx8fHo6ekZ1TL3l3q5ja0T9MqhiFKpFH19fRgaGvJRZIQQX6JkDqCzs3PMZO5vLfOCggIw\nDIO2tja740KhECKRiIYoEhKiKJljOJnHxcUFRMtcIpFAoVCgurp61DkqtRASulxO5larFdOnT8eK\nFSs8GY9PqNVq7N+/HwzD2CVufxrJYsMwDLKyslBTUzPqHCXziVEHKAlWLifzrVu3oqioCAzDeDIe\nn1Cr1di9e7ffj2SxKS0tHbUcLjC87vng4CBMJpMPogoMtDYLCVYuJfPW1lZ88cUXuPvuu4NiTZDO\nzk4A8PuRLDbFxcUwGo2jWuEMw0Amk1HrnJAQJHDlQY888gi2bNky7rjmkX/OVlRUoKKiwpWnmnIs\ny3Lbsfl7vdxGoVBAKpXizJkzmDt3rt05mUwGlUoFpVLpo+gIIc6qrKxEZWWlR+416WT++eefIz4+\nHtOnTx83iECpTfb09HAJ29FIliuXw/UHPB4PSUlJOHny5KhkHhUVhQsXLmBwcBACgUu/qwkhXnJl\nQ9edMuCkyyxVVVXYs2cPMjMzsXr1ahw8eBBr1651OQBfs41kAQKnZQ4MTx66cgVFYDjR01othISe\nSSfz559/Hi0tLWhsbMQHH3yApUuXYufOnVMRm1fYJgytXLnSLpn740iWkcrKynDp0iWHfRbR0dHQ\n6XQ+iMr/0dosJFi5Pc480EezdHZ2Qi6XY+XKlXatcH8dyWKTn5+PoaEhqFSqUedsuw/RbNDRAqX8\nR8hkuZXMFy9ejD179ngqFp9Qq9WQSCSQSqV2v5j8dSSLTVhYGBQKBY4fPz7qnFAoREREBPr6+nwQ\nGSHEF0J+BqharYZYLPb7afyOpKenj1oO14ZKLYSElpBP5p2dnQ6XvvXnzk+bwsJCNDU1OTxnG28e\nDPMACCETC+lkPt4Y80Bomc+YMQNqtRpWq3XUObFYDD6fD4PB4IPICCHeFtLJvLe3F0KhEEajEf/z\nP//DHbeNZPH3lnliYiKEQqHD9c0BKrU4Qh2gJFiFdDK3DUvs7e3Fyy+/zB03m80QCoV+O5JlpJyc\nHOzbt8/hOZlMRsn8CrQ2CwlWIZ/MY2NjR22QHAglFptly5bh2LFj0Gq1o85FRkZicHAQZrPZB5ER\nQrwppJMSpAHDAAAZ2ElEQVR5Z2cnpFKpw87PQEnmOTk5SE9PxyeffDLqHMMwVGohJESEfDKPiIhw\n2Pnp7/VyG4VCgeTkZBw6dMhh65xKLYSEhpBO5l1dXRCJRAE5ksWGz+cjKSkJs2fPxt69e0edl0ql\nMBgMGBwc9EF0hBBvCelk3tnZCT6fj+joaG7NDpZlYTabA6ZlDgBKpRIZGRk4ceLEqHM8Hg9RUVG0\n8NZltDYLCVYhm8yvHGNuG7JmMpkQFhYGHi9wXhqlUgmGYdDR0eFwQ+eYmBhoNBofROZ/aGgiCVaB\nk7E8TK/Xg2EYGI1Gux2GAqnz00apVEKtViM7O9vhsrgxMTHo7++HxWLxQXSEEG8I2WSuVquhUChg\nsVgQGRnJHQ+kzk8biUQChmGQmZmJs2fPjjrP4/Egl8vR1dXlg+gIId4Q0sk8JiYm4FZLdIRhGCiV\nSigUCofJHBge9dLV1UVrtRASpEI2mXd2dkIikQT0GPORlEolwsPDUV9f73CtloiICAiFwnH3bSWE\nBK6QTeZqtRrh4eFcvXzjxo0YGhqC2WyGSCTycXSTp1Qq0dPTg7i4uDFXUrS1zkMZdYCSYBXSyTws\nLIxrmW/atAkmkwkikSigRrLYxMXFQafTIScnZ8xSi1wuR19f36jlC0IJrc1CglXgZS0PUavV3HR3\nm0AtsQDDuwulpaVBJpONmcz5fD5kMhm6u7u9HB0hZKqFbDLv7OzE0NCQXTIPxJEsI+Xn58NqteLs\n2bNjdnTGxcVRRyghQSgkk3l/fz+sViusVuuoYYmB2jIHgJSUFAwNDUEoFOLMmTMOr4mIiACPx6P9\nQQkJMiGZzBsaGpCamgqpVGpXHw/kMgswPJ48Ly8P06ZNc7iKIjA8jJE6QgkJPiGZzM+dO4ekpCS7\nEstTTz0Fi8USkCNZRsrLy4NIJMKFCxfQ3Nzs8Bq5XI7e3t6QXHyL1mYhwSokk3l9fT3i4uLskvn6\n9eshFovtJhAFIoVCAZFIhEWLFuGzzz5zeI1AIEB0dHRIdoTS0EQSrEIumbMsi3PnziE6OjqoOj9t\nGIZBXl4eEhMTcezYMW4xsSvRjFBCgkvIJfOOjg4IhUJYrVa7BbYCvfNzpNzcXLS2tmLJkiVjts4l\nEgkA0NK4hASJkEvm586dQ15eHnp6eoJmjPmVoqKiEBsbi2nTpuHbb791uAMRwzBITU1FS0uLw+n/\nhJDAEpLJPCcnBwaDAVFRUdzxYCmz2OTl5UGlUmHx4sXYs2ePw2ukUikkEgna29u9HB0hxNNCMpkn\nJydDIpFwwxKtVitee+21gB/JMlJWVhba2tpw3XXX4eDBg2OWU1JSUtDd3Q2DweDlCH2DOkBJsAqp\nZG42m9HS0sItfWtjMpnw5ptvBvxIlpFEIhFSUlLQ09OD+fPn4/PPP3d4nVAoRFJSEi5evBgSnaG0\nNgsJVi4l85aWFixZsgTFxcUoKSnBtm3bPB3XlGhsbERqaioMBsOokSzBKC8vD+fOncPKlStx4MAB\nDA0NObxOoVAAAE0kIiSAuZTMhUIhXnnlFdTW1uLIkSN47bXXxpw+7k/OnTuH3NxcdHd3IyYmhjse\nrMk8LS0NWq0WkZGRkEqlaGhocHgdwzBIT0/HpUuXQnpFRUICmUvJPCEhAdOmTQMwPMStsLAQly5d\n8mhgU6Gurg4FBQXo6OiAUqnkjgdrvZjP53P7gs6cORPV1dVjXhseHo7Y2Fi0trZ6MUJCiKcI3L1B\nU1MTjh8/jjlz5tgdH9nRVFFRgYqKCnefyi0sy+Ls2bNYvXo16urqEBsbCwAYGhoK2mQODI85r6ys\nxPTp0/HHP/4R//iP/zjmtYmJiTh9+jR6e3vt+hQIIVOjsrISlZWVHrmXW8lcr9dj1apV2Lp1KzcJ\nxcbfRg10dnaCx+PBarUiLi6OG8mi1+sRERERtGt2JCQkwGq1QqFQoL29HTqdDjKZzOG1fD4faWlp\naG5uRlFREfh8vpejnXrB+j6TwHRlQ9edDnqGdXEIw8DAAG688UYsX74cDz/8sP1NGcbvRkYcOnQI\nP/zwAxYvXoyhoSHuL4m2tjYwDIOkpCQfRzh1Tpw4ga6uLhw/fhyzZs3CkiVLxr2+ubkZLMsiIyPD\nOwESQgC4lztdqpmzLItf//rXKCoqGpXI/VVdXR3y8/OhUqns6uW9vb12k4eCUWFhIS5evIji4uJx\n6+Y2KSkp0Ov10Gg0XoiOEOIJLiXz7777Drt27cLXX3+N6dOnY/r06di/f7+nY/Oos2fPIj8/H52d\nnYiPjwcwPFnIZDLZbVARjEQiEbc0bk1NzYRL3/L5fGRmZqKlpQUWi8VLURJC3OFSzXzhwoVjjln2\nRwaDAR0dHZDJZBCLxYiIiAAA9PX1ITIyMiA3cJ6ssrIyfPzxx0hOTsa33347YaklMjISSqUSjY2N\nyMvLC6oJVYQEo+DPYhhevzwrKwtdXV1ISEjgjvf19QV9icVGKpUiKSkJS5cuxc6dO8dcGnckWzmq\no6NjqsMjhLgpJJK5rcRy5fjyvr4+bgiev42+mQplZWVQq9VYuXIltm3bNuFqiQzDIDMzEx0dHejv\n7/dSlFMrFN5nEppCIpk7miw0MDAAi8XClVxCYc2OhIQEiEQilJWVgWEY7N27d8LHhIWFITU1FY2N\njUGxVG4ovM8kNAV9Mrdaraivr0dGRgb0ej03Waivrw8SiSSkasEMw6CsrAy1tbW477778Omnn0Kv\n10/4OLlcDolEEjKLcRESiII+mTc3N0Mul8NgMNhNFgqlevlIWVlZ0Gq1EIvFmD17tlOtc2B4nRej\n0YjOzs4pjpAQ4oqgT+Znz57lSiyh2vk5Ep/PR1FREU6dOoVVq1Zh//796O3tnfBxPB4P2dnZUKlU\nTl1PCPGukEjmhYWFaG9v5+rlFosFVqs1aLaJm6yioiLU19dDJpNh/vz5Y+4TeiWRSISsrCw0NjbC\nbDZPcZSEkMkI6mTOsizOnDmDnJwcdHR0cFP2ba3ykfXyUFqzIzIyEmlpaTh79ixWrVqFAwcO4Ikn\nnsD27dsnLKNERUUhMTERDQ0NAdkhGkrvMwktQZ3M1Wo1hoaGwLIsYmNjERYWBsDxFP5QG7JWWlqK\nU6dOQS6X44033sDtt98Os9mMP/3pTxM+Ni4uDhEREdwaLoEk1N5nEjqCOpmfOXOGW2s9OTkZwHBr\nPVTr5SMplUqIRCI0NzcjPDwcRUVF+Kd/+iccPXp0ws06GIZBWloazGYzTSgixE8EdTK3dX62trYi\nJSUFALhabzBt3uwKhmEwY8YMVFVVcbsLyWQyFBcXo6qqasLH2zpEOzs70dbWFnAtdEKCTdAn8+zs\nbGi1Wq7z01G9PFRlZWVBqVTi8OHD3LElS5bg66+/durxYWFhKCwshF6vR319/YQLeBFCpk7QJnO9\nXo/Ozk6EhYVBqVRCIBheU2zkFH4yvGjaxYsX0dTUBACYMWMG2tvbnd4GUCgUIi8vDxEREThz5kzQ\nTPsnJNAEbTKvra1Ffn4+2tvbnaqXh2rHmEgkwtKlS3Ho0CEYDAYIBAIsWrTI6dY5MFyySUlJQWpq\nKs6fPw+1Wu23ZZdQfZ9J8AvaZF5VVYU5c+agra2Nq5cbjUbw+XxuVMtIobxmR1JSEgoKClBZWQmW\nZbFs2TIcOHAAp0+fntR9ZDIZ8vPzoVar0djYyNXi/Ukov88kuAVlMjebzaiurkZZWRn6+/uhUCgA\nDA9VjImJ8XF0/mnWrFnQ6/Woq6tDUlIS/vVf/xVbtmyZdEIXi8UoKChAWFgYTp8+je7ubr9tpRMS\nTNza0NlfnThxAtnZ2dDr9UhKSgKPx4PFYoFWq0VJSYmvw/NLfD4f11xzDfbs2YPk5GSUlpbikUce\nwZYtW5CamorIyEjMnj0bS5cunfBePB4PKSkpiImJQXNzM7q7u5Genh7yI4gImUoub+g87k19vKHz\nK6+8gpKSEgiFQsTGxqK0tBQXL17kkowjvo7ZX1RXV6O1tRUrVqwAwzDo7OyESqWCXq/Hzp078fOf\n/xzXXHON0/djWRYdHR1QqVRISEhAfHy8T3d2oveZ+DOvb+jsz2wllpkzZ6KpqQmpqamwWCzQaDR2\nG1MQx6ZNm4bBwUH89NNPAID4+HiUlZVh/vz5eOqpp7B7926nxqHbMAyDhIQEFBQUoK+vD6dOnUJ7\ne7tf1tMJCWRBl8yrq6uRk5ODS5cuQalUQiaToaOjA7GxsRAKhWM+jtbsGMbj8bB06VIcP34cx44d\ns2slJCUl4cknn8T27dvR0NAwqfuKxWLk5uYiNzcXFosFtbW1aGpqgsFg8PR/YVz0PpNgFXRllhde\neAEzZ85EZ2cnfvaznyE2Nha1tbUoLi4eN5kTe3q9Hl9++SUiIiKwdOlSuxFA3333HXbv3o0tW7a4\nvPLk4OAgurq6uLkAsbGxiImJ4eYDEBKKqMxy2enTp9HU1ISUlBRERUUhISEBHR0dkMvllMgnSSKR\n4KabbkJERAT+/Oc/o6amhmtFL1iwAMXFxXjrrbfQ19eHyspK7Nu3b1LfhAKBAAkJCSgtLUVCQgJ6\ne3tx6tQpXLhwATqdjurahExS0LTMWZbFE088geXLl0OlUmHevHlITExEbW0tioqKHI4tJ85pa2tD\nXV0dGhsbkZCQgKysLCQkJGDTpk3o6upCaWkpuru7kZubi7vvvtvlpRIGBweh1WrR3d0Ns9mMmJgY\nyOVyREZG0vILJCS4kzuD5m/aqqoqWK1WpKWloaOjA6mpqbh06RJiYmIokbspOTkZycnJGBgYQFNT\nExobG1FVVYWFCxdixowZyMrKQn9/P37729/irbfewtq1a7nyi1arRUtLC0pLSydMyAKBAHFxcYiL\ni4PJZIJWq0VzczOGhoa4xB4eHk6JnRAHgqJlPjAwgIceegj33XcfGhoaUFJSgszMTJw6dQqFhYU0\nvnkKWK1WNDU14bvvvkNOTg6uuuoqWCwWvPrqqzh58iTS09MBAK2trYiMjERZWRl+85vfgM/nT+p5\nWJaFyWSCRqOBRqMBwzCIiIiAWCyGWCyGSCSCWCye9H0J8Ufu5M6gSOaffPIJzp49i8WLF6OzsxMr\nV66ESqWCxWJBRkaGU/fYuHEjrdvhAqPRiG+++QY6nQ6LFi1CYmIizGYz6urqMDQ0hOLiYgwODuLl\nl18Gj8fDunXrEB8f71LyZVkWRqMRRqMRJpMJZrOZ+5fP5yMyMhIRERHcv446U+l9Jv4spJO5VqvF\nI488gjvvvBMqlQo33XQT+vr6cOnSpUm1ymkyietYlkV9fT1++OEHREREYMaMGUhLS7MrhwwODmLn\nzp04evQoenp6kJubi7Vr1yInJ8cjz2+xWNDf34/+/n4YDAYYDAYIhUJIJBJERUUhKioKYWFh9D4T\nvxbSyfzVV18FAERHR+PGG29EV1cXhoaGkJGRAbFYPO5jKysrUVFRAcC/k/nIOP3ZwYMHkZqaiuPH\nj8NisSA9PR3p6elITEyEQCDgkrvJZMLhw4exa9cuzJw5E9dffz0yMzM9Wgu3lWf6+vqg1+vR19cH\nHo+HsrIyqNVq8Pl8MAwDhmHA4/G4z0d+D9g+DwsL434ReEugvOcUp2f5pAN0//79ePjhh2G1WnH3\n3Xdj/fr1rt7KJX19fTh69CiOHTuGBQsWYMGCBWhpaUF8fDwSEhKc+sELlDc4UOL85ptvsHHjRuTk\n5ECj0aC5uRnHjh3j9mIVCAQQi8VQKBSIj4/H+vXr8f333+O//uu/MDAwgOTkZOj1ehgMBkilUsjl\ncpSXl2PhwoWTHs/OMAzCw8MRHh6O+Ph4LrkDw987tr1hbf/aPkZ+39g+N5vNsFqtCA8Ph1gsRnh4\nOEQiEXd+5A+fQCCASCSy++XlikB5zylO/+FSMrdarXjggQdw4MABJCcnY/bs2Vi5ciUKCws9HZ+d\noaEhHD16FJ9++ilaWlogl8tx1VVXYe7cuTCZTMjNzUVERMSUxkAmxjAMYmNjERsbixkzZgAYfu8G\nBwdhMBigVquhVqtx8eJFmM1mzJkzB3w+HwKBACkpKUhKSoLBYEB3dzeqqqqwa9culJSUwGq1wmQy\nQafTcbtHrVixAvPmzYNAIIDVasXQ0BCA4ZmsI+vytuQOAJmZmZP+Pw0ODsJkMnE1+97e3lH/55HX\nsSzLddCKRCIIhULu/2j74PF4dr9IbL9YgOFO/ZGzYxmG4R7vy7VtiP9yKZkfPXoUOTk5XOfiL3/5\nS3z22WceS+Ysy6K3txcSiQQ8Hg/t7e04fPgwvv76a7Asi8LCQmRlZSEtLQ0JCQmQy+VITEykb3I/\nxuPxuHKFTCZDbm4ugOH3uqenB2q1mlvUq7a2lnuc7Rd2X18fRCIRJBIJkpOTIRKJoNPp8OGHH+L3\nv/89WJbF4OAg9z0wNDQEHo8HoVAImUzG/XIBgH379sFkMqGnpwcsy0IkEiE6OhpZWVnc6o4Mw8Bs\nNkOn03F/KchkMkgkEqf+v4ODg3YdtCaTCYODgxgcHITVauX+HVnisX0OAD09PdzuT7bXyfaYkYnd\nVi7i8Xjc40feEwD3+cjrbL/srrzWdv2Vj7M91haL7WNgYAD9/f121418rKO/fmzfD1d+jHwsmTyX\nauYfffQR/vrXv2L79u0AgF27duH777/n6tf0ZhBCiGu8WjOfKFn7a0ciIYQEK5fqEsnJyWhpaeG+\nbmlpGXOdcEIIIVPPpWQ+a9Ys1NfXo6mpCRaLBR9++CFWrlzp6dgIIYQ4yaUyi0AgwO9+9ztcd911\nsFqt+PWvfz3lI1kIIYSMzeXhH8uXL8fWrVshEAiwY8cOvPjiiw6vq6ysxPTp01FSUuKzcZ779+9H\nQUEBcnNzHcb50ksvYfr06Zg+fTpKS0shEAig0+n8Ls6uri5cf/31mDZtGkpKSvDOO+94PUZg4ji1\nWi1uueUWlJeXY86cOXajU7zlrrvuglKpRGlp6ZjX/Mu//Atyc3NRXl6O48ePezG6/zdRnGfPnsW8\nefMgFovx8ssvezm6/zdRnO+99x7Ky8tRVlaGBQsW4OTJk16OcNhEcX722WcoLy/H9OnTMXPmTBw8\neNDLETr3vQkAP/zwAwQCAf7yl784d2PWRYODg2x2djbb2NjIWiwWtry8nD19+rTdNVqtli0qKmJb\nWlpYlmVZtVrt6tO5zJk4R9q7dy97zTXXeDHCYc7E+fTTT7OPP/44y7LDr6VcLmcHBgb8Ls5HH32U\nfeaZZ1iWZdmzZ8/65PX85ptv2OrqarakpMTh+f/93/9lly9fzrIsyx45coSdM2eON8PjTBRnZ2cn\n+8MPP7AbNmxgX3rpJS9H9/8mirOqqorV6XQsy7Lsvn37/Pb11Ov13OcnT55ks7OzvRUaZ6IYWXb4\n52zJkiXsP/zDP7AfffSRU/d1uWU+cqy5UCjkxpqP9P777+O2227jOkcVCoWrT+cyZ+Ic6f3338fq\n1au9GOEwZ+JMTEzkJqv09vYiNjbW6zvzOBPnmTNnsGTJEgBAfn4+mpqaoFarvRrn1VdfjZiYmDHP\n79mzB3fccQcAYM6cOdDpdOjo6PBWeJyJ4oyLi8OsWbN8vrnKRHHOmzcP0dHRAIZfz9bWVm+FZmei\nOCMjI7nP9Xq9T3LSRDECw8uUrFq1CnFxcU7f1+Vk3tbWhtTUVO7rlJQUtLW12V1TX18PjUaDJUuW\nYNasWfjjH//o6tO5zJk4bQwGA/7617/itttu81Z4HGfivOeee1BbW4ukpCSUl5dj69at3g7TqTjL\ny8u5Pw2PHj2K5uZmn/1wj8XR/8PfYgxUf/jDH3DDDTf4OowxffrppygsLMTy5cuxbds2X4czSltb\nGz777DPcf//9AJyft+Nys86ZJxgYGEB1dTW++uorGAwGzJs3D3PnzuVm/3nDZCYw7d27FwsXLoRM\nJpvCiBxzJs7nn38e06ZNQ2VlJRoaGrBs2TLU1NQgKirKCxEOcybOxx9/HA899BDXBzF9+nS/XG+c\nvWI+BE12c9/XX3+NHTt24LvvvvN1KGO6+eabcfPNN+Pbb7/F7bffjrq6Ol+HZOfhhx/G5s2buRm0\nV36fjsXlZO7MWPPU1FQoFApuwaNFixahpqbGq8l8MmPiP/jgA5+UWADn4qyqqsKGDRsAANnZ2cjM\nzERdXR1mzZrlV3FGRUVhx44d3NeZmZnIysryWozOuPL/0draiuTkZB9GFPhOnjyJe+65B/v375+w\njOAPrr76agwODqK7u5tb6sEfHDt2DL/85S8BDA962LdvH4RC4YTDv10uszgz1vymm27C3//+d1it\nVhgMBnz//fcoKipy9SmnLE5geC2Mb775BjfddJNX47NxJs6CggIcOHAAANDR0YG6ujqvJ0ln4uzp\n6YHFYgEAbN++HYsXL3Z6TRNvWblyJXbu3AkAOHLkCGQyGZRKpY+jGpuzrTNfuXjxIm699Vbs2rXL\nI2vUT5WGhgbutayurgYAv0rkAHDhwgU0NjaisbERq1atwuuvv+7UPB6XW+ZjjTV/8803AQD33nsv\nCgoKcP3116OsrAw8Hg/33HOP15O5M3ECw3W06667btJLrXozzv/4j//AunXrUF5ejqGhIfznf/4n\n5HK538V5+vRp3HnnnWAYBiUlJfjDH/7g1RgBYPXq1Th06BC6urqQmpqKTZs2YWBggIvxhhtuwBdf\nfIGcnBxERkbi7bff9nqMzsSpUqkwe/Zs9Pb2gsfjYevWrTh9+rTXfzlOFOczzzwDrVbL1XmFQiGO\nHj3q1RidifPjjz/Gzp07uY1LPvjgA7+L0VVTsjkFIYQQ76I1YwkhJAhQMieEkCBAyZwQQoIAJXNC\nCAkClMwJISQIUDInhJAg8H8jD4qPpQSbtAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1af25550>"
       ]
      }
     ],
     "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.0231170117\n",
        "The proportion of voxels with rRMSE<1.0 is 99.8646169971\n",
        "The proportion of voxels with rRMSE<1.0 is 99.8486039538\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}