{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Figure 5 Comparing the DTM and the SFM\n",
      "\n",
      "In this notebook, we compare the accuracy of the DTM and the SFM to each other directly in each b value. \n",
      "\n",
      "First, we import some modules: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.stats as stats\n",
      "\n",
      "import nibabel as ni\n",
      "\n",
      "import osmosis.viz.mpl as viz\n",
      "import osmosis.utils as ozu\n",
      "import osmosis.model.sparse_deconvolution as ssd\n",
      "import osmosis.model.dti as dti\n",
      "import osmosis.model.analysis as oza"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We set the data-path to point to the installation-specific location where data might be found:"
     ]
    },
    {
     "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": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Point to the data-files of 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": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Get the white-matter mask: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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",
      "mask[np.where(wm_nifti==1)] = 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We set the solver parameters for Elastic Net according to the calculations in the AppendixA notebook:"
     ]
    },
    {
     "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": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Initialize the class instances of the SFM objects. Use AD/RD derived from the corpus callosum ROI distribution for each b value. Parameters are not saved, or read from a file (hence `params_file=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=mask, solver_params=solver_params, params_file = 'temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'])\n",
      "SD_1k_2 = ssd.SparseDeconvolutionModel(*data_1k_2, mask=mask, solver_params=solver_params, params_file = 'temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'])\n",
      "ad_rd = oio.get_ad_rd(subject, 2000)\n",
      "SD_2k_1 = ssd.SparseDeconvolutionModel(*data_2k_1, mask=mask, solver_params=solver_params, params_file = 'temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'])\n",
      "SD_2k_2 = ssd.SparseDeconvolutionModel(*data_2k_2, mask=mask, solver_params=solver_params, params_file = 'temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'])\n",
      "ad_rd = oio.get_ad_rd(subject, 4000)\n",
      "SD_4k_1 = ssd.SparseDeconvolutionModel(*data_4k_1, mask=mask, solver_params=solver_params, params_file = 'temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'])\n",
      "SD_4k_2 = ssd.SparseDeconvolutionModel(*data_4k_2, mask=mask, solver_params=solver_params, params_file = 'temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'])"
     ],
     "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": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Likewise, initialize the DTM objects with the same data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "TM_1k_1 = dti.TensorModel(*data_1k_1, mask=mask, params_file='temp')\n",
      "TM_1k_2 = dti.TensorModel(*data_1k_2, mask=mask, params_file='temp')\n",
      "TM_2k_1 = dti.TensorModel(*data_2k_1, mask=mask, params_file='temp')\n",
      "TM_2k_2 = dti.TensorModel(*data_2k_2, mask=mask, params_file='temp')\n",
      "TM_4k_1 = dti.TensorModel(*data_4k_1, mask=mask, params_file='temp')\n",
      "TM_4k_2 = dti.TensorModel(*data_4k_2, mask=mask, params_file='temp')"
     ],
     "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": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculate the rRMSE for each b value and each model, based on cross-prediction: a model is fit to one data-set and used to predict the other data-set from the same b value. The rRMSE is the average of the two folds of this process (1=>2 and 2=>1)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dtm_rrmse_1k = oza.cross_predict(TM_1k_1, TM_1k_2) \n",
      "dtm_rrmse_2k = oza.cross_predict(TM_2k_1, TM_2k_2)\n",
      "dtm_rrmse_4k = oza.cross_predict(TM_4k_1, TM_4k_2)\n",
      "\n",
      "sfm_rrmse_1k = oza.cross_predict(SD_1k_1, SD_1k_2) \n",
      "sfm_rrmse_2k = oza.cross_predict(SD_2k_1, SD_2k_2)\n",
      "sfm_rrmse_4k = oza.cross_predict(SD_4k_1, SD_4k_2)\n"
     ],
     "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": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In each b value, we will plot a 2D, or image histogram of all voxels and their rRMSE values in the two models:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "idx = np.logical_and(np.logical_and(np.isfinite(dtm_rrmse_1k), dtm_rrmse_1k<1.5), np.isfinite(sfm_rrmse_1k))\n",
      "fig, cbar = viz.scatter_density(np.hstack([0.5, dtm_rrmse_1k[idx], 1.5]), np.hstack([0.5, sfm_rrmse_1k[idx], 1.5]), \n",
      "                                          return_cbar=True, vmin=0, vmax=3, cmap=matplotlib.cm.gray_r)\n",
      "fig.axes[0].plot([0,100],[100,0], 'k--')\n",
      "fig.axes[0].plot([100,0],[50,50], 'k--')\n",
      "fig.axes[0].plot([50,50],[100,0], 'k--')\n",
      "\n",
      "cbar.set_ticks([0,1,2,3])\n",
      "cbar.set_ticklabels([1,10,100,1000])\n",
      "\n",
      "fig.set_size_inches([10,10])\n",
      "fig.savefig('figures/scatter_rrmse_1k.svg')\n",
      "fig_hist = viz.probability_hist(dtm_rrmse_1k[idx] - sfm_rrmse_1k[idx])\n",
      "fig_hist.set_size_inches([8,6])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "-c:1: RuntimeWarning: invalid value encountered in less\n",
        "/home/arokem/usr/lib/python2.7/site-packages/osmosis/viz/mpl.py:353: RuntimeWarning: divide by zero encountered in log10\n",
        "  imax = ax.matshow(np.log10(np.flipud(data_arr.T)), cmap=cmap, **kwargs)\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAJECAYAAAD0VzVyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X90VOWdx/HPJAGBBPmVkJAfGCCBgMYYEKstaFoWkEIV\nrUehe4SC7dpy7MLp6QrddluCXdGuttpytFZsD7Yr0F9KsW60LQ0UKaARlRJSIpoaIIIsgjL5HWb/\nYJPmzh2SCcnce+e579c5OeWZeW7uk1tn5jvf+32eJxAKhUICAAAwTILbAwAAAIgFghwAAGAkghwA\nAGAkghwAAGAkghwAAGAkghwAAGAkghwAABAzS5cuVXp6ugoLCzseO3XqlGbOnKnx48dr1qxZOn36\ndMdza9euVX5+vgoKCvTSSy91PF5RUaHCwkLl5+dr+fLlUZ2bIAcAAMTMkiVLVFZWZnnsgQce0MyZ\nM3Xo0CHNmDFDDzzwgCSpsrJSmzdvVmVlpcrKyrRs2TK1L+f35S9/WU899ZSqq6tVXV1t+52REOQA\nAICYmT59uoYNG2Z57Le//a0WL14sSVq8eLGee+45SdKWLVu0cOFC9evXT7m5ucrLy9OePXtUV1en\njz76SNdcc40kadGiRR3HdIUgBwAAQw0fPlyBQMCxn8GDB0c1ruPHjys9PV2SlJ6eruPHj0uSjh07\npuzs7I5+2dnZOnr0qO3xrKwsHT16tNvzEOQAAGCoDz74wNHznT17tsfHtAdIsUCQAwAAHJWenq73\n3ntPklRXV6eRI0dKOp+hqa2t7eh35MgRZWdnKysrS0eOHLE8npWV1e15CHIAADCYk7eronXTTTdp\nw4YNkqQNGzZo/vz5HY9v2rRJzc3Neuedd1RdXa1rrrlGGRkZuvTSS7Vnzx6FQiH97Gc/6zimK0kX\nd8kAAAC6t3DhQm3fvl0nT55UTk6O1qxZo1WrVun222/XU089pdzcXP3iF7+QJE2aNEm33367Jk2a\npKSkJD322GMdwdNjjz2mz3/+82poaNCnP/1p3Xjjjd2eOxBqn5sFAACMEst6l0hCoZC8FFaQyQEA\nwGBOBzleQk0OAAAwEpkcAAAM5mQmx2vI5AAAACORyQEAwGBkcgAAAAxDJgcAAIORyQEAADAMmRwA\nAAxGJgcAAMAwBDkAAMBI3K4CAMBg3K4CAAAwDJkcAAAMRiYHAADAMGRyAAAwGJkcAAAAw5DJAQDA\nYGRyAAAADEMmBwAAg5HJAQAAMAyZHAAADEYmBwAAwDAEOQAAwEjcrgIAwGDcrgIAADAMmRwAAAxG\nJgcAAMAwZHIAADAYmRwAAADDkMkBAMBgZHIAAAAMQyYHAACDkckBAAAwDEEOAAAwErerAAAwGLer\nAAAADEMmBwAAg5HJiaGysjIVFBQoPz9fDz74oO357373u0pMTNQll1yirKwsfec734n6WJzX3XV6\n6KGHVFxcrOLiYhUWFiopKUmnT5+O6lic1911+uCDD3TLLbeoqKhIH/vYx3TgwIGoj4W0dOlSpaen\nq7CwMOLzVVVVuu666zRgwAA9/PDDlue4vtHp7hpv2bJFRUVFKi4u1pQpU7Rt27aO57jGiFuhGGpt\nbQ2NGzcu9M4774Sam5tDRUVFocrKSsvzmZmZoRkzZtie7+5YnNfT67R169bQjBkzLupYv4rmOn3t\na18LrVmzJhQKhUJVVVVc4x7asWNH6LXXXgtdccUVEZ8/ceJE6JVXXgl94xvfCD300EMdj3N9o9fd\nNT579mzHv998883QuHHjQqEQ1zjeSQqlpqY69hPjsKLHYprJ2bt3r/Ly8pSbm6t+/fppwYIF2rJl\ni+X5rKwsDRw40PZ8d8fivJ5ep2eeeUYLFy68qGP9KprrdPDgQX3yk5+UJE2YMEE1NTU6ceIE1zhK\n06dP17Bhwy74fFpamq6++mr169fP8jjXN3rdXePk5OSOf589e1apqamSuMaIbzENco4ePaqcnJyO\ndnZ2to4ePWp5fuTIkdq1a5eKioq0adMm7d+/P6pjcV5PrlN9fb1efPFFffazn+3xsX4WzXUqKirS\nb37zG0nnPxT+/ve/68iRI1zjGOP69q3nnntOEydO1Jw5c/SDH/xAEtfYBIFAwLEfr4lpkNPdHxwI\nBDRixAjV1tbqjTfe0MyZM/XCCy/EckjG6cl/VFu3btW0adM0dOjQHh/rZ9Fcp1WrVun06dMqLi7W\nunXrVFxcrMTERK5xjHF9+9b8+fN18OBBbd26VXfeeadCoZDbQwJ6Jaazq7KyslRbW9vRrq2tVXZ2\ntuX548ePa9CgQZKk4cOHKxAI6NSpU8rOzu7yWJzX3TXubNOmTR23qnp6rJ9Fc50GDx6sn/zkJx3t\nMWPGaNy4cWpoaOAaxxD/DcfG9OnT1draynuxIfz8ZSCmmZyrr75a1dXVqqmpUXNzszZv3qybbrrJ\n8vzBgwf1zjvvqLm5WT/96U81cOBADR8+vNtjcV601+nMmTPasWOHbr755h4f63fRXKczZ86oublZ\nkvTkk0/qhhtuUEpKCte4j4VnFri+fefw4cMd1/e1116TJI0YMYJrjLgW00xOUlKS1q1bp9mzZ6ut\nrU133XWXJk6cqCeeeEKSdPfdd2vu3LkqKCiQJGVkZOjXv/51l8fCKpprLJ2/1z579mwNHDiw22Nh\nFc01rqys1Oc//3kFAgFdccUVeuqpp7o8FlYLFy7U9u3bdfLkSeXk5Ki0tFQtLS2Szl/f9957T1On\nTtWHH36ohIQEPfroo6qsrFRKSgrXN0rdXeNf//rXevrpp9WvXz+lpKRo06ZNkvhv2AR+zuQEQtx0\nBQDASIFAQOnp6Y6d7/jx456q5WJbBwAAYCS2dQAAwGB+vl3VZSanu2XAy8vLNWTIkI4tAyJtyZCV\nlaW0tDSWAwcAAI7qMpOzZMkSfeUrX9GiRYsu2OeGG27Qb3/7W8tjbW1tuueee/Tiiy9q1qxZGjJk\niDZt2qRFixbppptu6iha83N0CQDwLyfrVvz8WdtlJqe7ZcClyP9HtS8DfuLECeXn52vRokV64YUX\nIi4HHgqFdMMNNygUCvETw59vf/vbro/B9B+uMdfYhB+ucex/4JxeFR4HAoGOLRk+/elPq7KyUtI/\nlgFv/9/2ZcAjLQdeUlKi119/XSUlJXrkkUd6MxwAADznkUceUUlJiUpKSnT99dc7fn4/b+vQq8Lj\nyZMnq7a2VoMGDdL//M//aP78+Tp06FDH89H8weXl5SopKVF5eXlvhgIAgCetWLFCK1asUDAY1Ny5\nc90ejq/0KpMzePDgji0Z5syZo5aWFssy4O1LrrcvA36h5cDnz5/fm2EgCiUlJW4PwXhc49jjGsce\n1zg22gOcMWPGOH5uP2dyehXkdF70Z+/evQqFQpYtGVJTU3Xo0CH97Gc/05w5cy64HPiKFSt6MwxE\ngTeu2OMaxx7XOPa4xrHx4x//WGPHju1YDR3O6PJ2VXfLgP/qV7/S448/rqSkJA0aNMi2DPjcuXPV\n2Niojz76SLfeeivLgQMAfGn58uWSpIQE59fg9WKGxSmubusQCASoNAcA+IqTn32BQMDRXeOPHDni\nqc91tnUAfGz16tVuDwEAYoZMDuBjvAaBvhcMBtXQ0KDU1NSIzzudycnJyXHkXJJUW1vrqfcUMjkA\nAPSRYDCoefPmad26dW4PBWKDTgAA+kR7gJObm6v/+I//cHs4HfxceEwmBwCAXuoc4Kxfv16JiYlu\nDwkikwMAQK80NjZ6OsDxcyaHIAfwsW9/+9tuDwGIe5dccom+8IUvaMGCBZ4LcPyO2VUAADjI6dlV\nubm5jpxLkmpqajz1uU5NDgAAMBK3qwAAMJifa3LI5AAAEKVgMKgvfvGL+t///V+3h4IoEOQAABCF\n9mnira2tGjp0qNvDiVogEHDsx2sIcgAfY+8qIDqsgxOfmF0F+BivQaB7fR3gOD27auzYsY6cS5Le\nfvttT72nUHgMAEAXfvazn8V1BseLt5GcQiYH8DFeg0D3QqGQQqGQEhL6psLD6UzOuHHjHDmXJB0+\nfNhT7ylkcgAA6IJXi2qjFc9j7y0KjwEAgJHI5AA+xt5VgFUwGNSZM2eUmZnp9lD6jJ8zOdTkAAB8\np6qqytIuKCjomEX18Y9/XP/5n/8Zs3M7XZOTn5/vyLkkqbq62lOf62RyAAC+13ma+Jo1a9weTp/y\ncyaHmhwAgK/V19ez0J+hyOQAAHyrublZX/7ylzVp0iRjAxw/Z3IIcgAAvlNQUCDp/Bo49957r+bP\nn29kgON33K4CfIy9q+B3gUBAn/3sZ40OcPy8QSezqwAf4zUIOM/p2VXtWSsnVFVVeeo9hUwOAAAw\nEkEOAMAXgsGgFi1apOPHj7s9FEf5+XYVhccAAKPU1dVZ2qNGjVIwGNTcuXM1duxYpaWluTQyOI0g\nBwBgtM4Bzvr16/tsN/F44cUMi1P89f80AAv2roLp6uvrfR3g+B2ZHMDHmEIO0z3//PO+D3D8nMlh\nCjkAwGihUMhTH/ROTyG//PLLHTmXJB04cMBTn+tkcgAARvNSgOMGP//9/szdAQAA45HJAQAYIRgM\n6v3331dubq7bQ/EUMjkAfInCY8Sbqqoqy0+79mnijz/+uIujg9cQ5AA+Vlpa6vYQgF7rvA7O2rVr\n3R4OPITbVQCAuOX3hf6iwe0qAADiTEtLCwEOukQmBwAQl/r166d/+7d/05w5cwhwuuDnTA5BDgAg\nbhQUFHTZBjojyAF8jL2rAPP5OZPDtg4AADjI6W0drrrqKkfOJUmvv/66pz7XuYkJAPC8YDCohQsX\n6siRI24PJe4EAgHHfryGIAcA4GnBYFDz5s3TgAEDNGrUKLeHgzhCTQ4AwLPaA5zc3FytX79eiYmJ\nbg8p7ngxw+IUMjkAAE8iwEFvEeQAPsbeVfCyF154gQCnD/i5JofZVYCP8RqE14VCIU9+ePaG07Or\npkyZ4si5JKmiosJT7ylkcgAAnmVagANnUXgMAIDB/BwokskBALguGAzq0KFDbg8DhiHIAQC4qn0W\n1Y9+9CO3h2IkPxcec7sK8DH2roJT6urqLO3Dhw9LkhoaGrRy5UplZGTov/7rv9wYGgxGkAP4GFPI\n4abOAc7KlSuZJh4jXsywOIXbVQAAx7W1tWnVqlUEOIgpMjkAAMclJiZq0aJFuuqqqwhwYszPmRyC\nHACAK5xcpA7+RJADAIi58N3D2U3cOX7O5FCTA/gYhccATEaQA/hYaWmp20OADwSDQX32s5/V22+/\n7fZQ4DPcrgIAxEwwGNTcuXM1duxY5ebmuj0cX/Lz7SqCHABAr1RUVFjamZmZkqT6+nrdeeeduuyy\ny7R+/XolJHDzAM4iyAEA9LnOAc5DDz1EgOMiP2dy+K8OANDndu7cqdzcXAIcuIpMDuBj7F2FWJk1\na5ZmzZrl9jAgMjkAfIop5ABMRiYHANAj4YXGrFzsbWRyAAC4SMFgUPv373d7GIANQQ4A4KI1NDRo\n7ty5Wr9+vdtDwQUEAgHHfryGIAcAcFEaGhq0YsUKjR07Vt///vfdHg5gQ5AD+BiFx7gYU6ZMUUFB\ngb75zW/qyiuvZKE/j/NzJicQCoVCrp08EJCLpwd8j9cgLsa5c+c0c+ZMVjK+SE6+7gKBgK6//npH\nziVJO3bs8NR7CrOrAAA9kpCQoDVr1ui6664jwIGnEeQAAHrsE5/4hNtDQJS8eBvJKYTgAADASGRy\nAAAwGJkcAL7E3lXoTjAY1M0336y//e1vbg8F6DEyOYCPMYUcXQkGg5o3b55yc3OVl5fn9nBwkcjk\nXMDSpUuVnp6uwsLCiM9v2bJFRUVFKi4u1pQpU7Rt27aO58rKylRQUKD8/Hw9+OCDfTtqAEBMdQ5w\n1q9fr8TERLeHBPRYl+vk/PnPf1ZKSooWLVoUcV+SYDCo5ORkSdL+/ft1yy236K233lJbW5smTJig\nP/zhD8rKytLUqVO1ceNGTZw40Xpy1ugAAMfU1dVZ2ocPH7b1GTdunOrr67Vo0SLl5OTov//7vwlw\n+pjT6+R86lOfcuRckrRt2zZPfa53mcmZPn26hg0bdsHn2wMcSTp79qxSU1MlSXv37lVeXp5yc3PV\nr18/LViwQFu2bOmjIQMAYunVV19Vbm6uHnroIQIcxLVe1+Q899xz+vrXv666ujq99NJLkqSjR48q\nJyeno092drb27NkT8fjONQElJSUqKSnp7ZAAAL1w/fXXO7pKrunKy8tVXl7u2vn9XJPT6yBn/vz5\nmj9/vv785z/rzjvvVFVVVY+Op/ARcM/q1at5DQIxFv4FvrS01L3B+Eyfza6aPn26WltbderUKWVn\nZ6u2trbjudraWmVnZ/fVqQD0kdLSUoIcQ4XX30jSsWPHLO1x48bZ+owaNSpmYwKc1qt1cg4fPtxR\nYPTaa69JkkaMGKGrr75a1dXVqqmpUXNzszZv3qybbrqp96MFAPSp+vp6vfHGG24PAzHk513Iu8zk\nLFy4UNu3b9fJkyeVk5Oj0tJStbS0SJLuvvtu/frXv9bTTz+tfv36KSUlRZs2bTr/S5OStG7dOs2e\nPVttbW266667bDOrAADuap9FlZ+fr6KiIreHA/S5LqeQx/zkTCEHXMVr0Fzd3a5qaGjQypUrlZOT\nY5lFxe2q2HN6CvnMmTMdOZck/f73v/fUewrbOgCAzzQ0NGjFihW2AAcwDds6AD7G3lXmOnPmjO2x\n5ORkhUIhLV++XKNHj2ahP5/wYq2MUwhyAB9jZpX/BAIBrVixQhMnTiTAgfEIcgDAZ6644gq3hwAH\n+TmTQ00OAAAwEpkcAIgzkVaWLygosLSHDBkiSQqFQh3f5Jk55U9kcgAAxqmvr9edd96pyspKt4cC\nuIIgB/AxCo/N1b7QX1pamiZMmOD2cOAiP694TJAD+BgbBZopGAxq0aJFrIMD3yPIAQCDBINBzZs3\njwAHEIXHABB32ouKO2svRn7ttdeUnp7OQn/o4MXbSE4hyAEAg0yePFmTJ08mwAFEkAMAgNH8nMmh\nJgfwMfauAmAyMjmAjzGF3H11dXWWdqQF+yoqKiztKVOmSDpfZLxv3z5NmzaNhf5wQWRyAABxJRgM\nau7cufr5z3/u9lAAzyKTAwBxpj3AGTt2rB577DG3hwOPI5MDAIgLDQ0NHQHO+vXrlZDA2zhwIWRy\nACBOhEIhrVy5UhMmTCDAQdT8nMkhyAF8bPXq1RQfu+zw4cPd9klOTu7498qVKzVv3jwCHCAKvEoA\nH2Pvqvgzfvx4Ahz0CBt0AgAAGIYgBwA8KhQKuT0EIK5RkwMAHnLmzBlJUn19vf71X/9VK1as0G23\n3ebyqBDPvHgbySlkcgDAY+rr6/XlL39ZI0eO1MSJE90eDhC3yOQAPsbeVd7THuBkZWXpvvvuYzdx\n9BqZHAC+xPRxb2loaCDAAfoQmRwA8Ii///3vGjdunL7xjW8Q4KDP+DmTQ5ADADESvnv422+/bevT\neffwgoICzZ07N+bjAvyCIAcAAIP5OZNDTQ4AAIiZpUuXKj09XYWFhR2PnTp1SjNnztT48eM1a9Ys\nnT59uuO5tWvXKj8/XwUFBXrppZd6dW6CHMDHKDx2T2Njo9588023hwEfcHtbhyVLlqisrMzy2AMP\nPKCZM2fq0KFDmjFjhh544AFJUmVlpTZv3qzKykqVlZVp2bJlOnfu3MX/7SEXl9QMBAKs6Am4iNfg\nxaurq7M9dvDgQUu7f//+lvbAgQMlnZ9FtWLFCuXk5OjZZ5+N3SDhSU6+7gKBgKOLSf7qV7+K+LfV\n1NToM5/5jPbv3y/pfP3Z9u3blZ6ervfee08lJSWqqqrS2rVrlZCQoJUrV0qSbrzxRq1evVrXXnvt\nRY2HmhwAcFB7gJOZmamvf/3rbg8H6JX3339f77//fo+PO378uNLT0yVJ6enpOn78uCTp2LFjloAm\nOztbR48evejxEeQAgEM6Bzjf/OY3mSYOR8Sy8HjkyJEaOXJkR7uysrLHv6O7Hcx7M35qcgDAId/4\nxjcIcACp4zaVdP7Wb3uglJWVpdra2o5+R44cUVZW1kWfh0wOAFyE9o00OwuvwZk2bZql/Z3vfEdj\nxowhwIGjvDiF/KabbtKGDRu0cuVKbdiwQfPnz+94/HOf+5y++tWv6ujRo6qurtY111xz0echyAF8\njL2rnJWXl+f2EADHLVy4UNu3b9fJkyeVk5OjNWvWaNWqVbr99tv11FNPKTc3V7/4xS8kSZMmTdLt\nt9+uSZMmKSkpSY899livgjRmVwHARaiqqrI9dvLkSUs7PJMTaUZW5xWP4Q9Oz6664447HDmXJG3e\nvNlTn+vU5ABADPRmbQ8AfYMgBwD6WENDg7761a/q5ZdfdnsogOuLAbqJmhwAiEL47anwhf8kKS0t\nTQ0NDVq5cqUyMjI0dOhQy3EFBQUxHyeAfyDIAYA+0jnAWblyJbOo4AlezLA4hdtVgI+xd1XfaWxs\nJMABPIYgB/Cx0tJSt4dgjBMnTmjs2LEEOPAcanIAwMfCp3YfPnzY1qd9c812EydOtLWXL1/e94MD\ncNHI5AAAACORyQEAwGBevI3kFDI5ANBDDQ0N+tOf/uT2MAB0gyAH8DH2ruq5hoYGrVixQn/84x/d\nHgoQFQqPAfgSU8gjO3v2rO2x1NRU1dfX66tf/aouu+wyPfroo0pI4Hsi4GUEOQAQhfr6en3pS19S\nTk6O7rvvPgIcxA0vZlicwqsUAKKwatUqAhwgzpDJAYAo3HvvvcrMzCTAQdzxcyaHIAeA0cI31gwG\ng7Y+x44ds7QzMzNtfZKTky3tUaNG9cHoAMQSX0kAH6PwGDCfn2dXEeQAPsbeVZG1tbW5PQQAfYAg\nBwA6aWxs1LJly/TKK6+4PRQAvURNDgD8v8bGRq1Zs0ZjxozR5MmT3R4O0Ce8eBvJKQQ5AOJW+O7h\nkYqBT548aWlHWuhv2LBhamho0Le+9S3l5OTo2WefVWJiYpfnAuB9BDkAfK+hoUH33nuvRo0apZUr\nV9oCHCCekckB4EvsXXXeyZMnNX78eC1btowABzAIQQ7gY0whPy8nJ0df+cpX3B4GEBNkcgAgDp05\nc8bSPnz4sK3PBx98YGm3trba+uTm5nZ7Lhb/A+IPQQ4AAAbzcyaHdXIA+EpjY6N27drl9jAAOIAg\nB4BvNDY26v7771dFRYVCoZDbwwEcwbYOAHzJT4XHjY2Nuu+++5SWlqZly5Z58g0ZQN8KhFz8OhMI\nBPg2Bbgonl6DkRbj2717t6V9+vRpW5/8/Hw1NDRo5cqVysjI0NatW5kmDlc5+boLBAL64he/6Mi5\nJOnJJ5/01HsKmRwAxnvwwQeVkZHBQn+AzzC7CoDxli1bphEjRhDgAD5DkAPAeCNHjnR7CIBr/Fx/\nRpADwBO622zzj3/8o+2YpCTrW1ikBftSU1P7YHQA4hFBDuBjJu5d1dbWxm0poBM/Z3IoPAZ8zLQp\n5I2Njfr3f/937d271+2hAPAAMjkAjNDY2Ki1a9cqPT1dU6ZMcXs4gGf4OZNDkAPAk9avX29pnzt3\nztYnLS1N0vkA5+GHH1Z+fr7uu+8+y+2qgoKC2A4UgGcR5ACIa+1bNaSlpdkCHABkcgAgbn344YfK\ny8vTP//zPxPgALDosvB46dKlSk9PV2FhYcTnq6qqdN1112nAgAF6+OGHLc+VlZWpoKBA+fn5evDB\nB/tuxAD6jAmFxyNHjtSiRYsIcIALYIPOC1iyZInKysou+PyIESP0wx/+UF/72tcsj7e1temee+5R\nWVmZKisrtXHjRh08eLBvRgygz5SWlro9BACImS5vV02fPl01NTUXfD4tLU1paWn63e9+Z3l87969\nysvLU25uriRpwYIF2rJliyZOnNjrAQOIfzt37rQ9Fr4YYHhmJj093XZMfn6+pT1kyJBuf2+kBQMB\nmCkmNTlHjx5VTk5ORzs7O1t79uyJ2LdzurykpEQlJSWxGBIAAzQ2Nurll1/WjBkz3B4KELXy8nKV\nl5e7dn4v3kZySkyCnJ5cUBNqAgDEXmNjo9asWaOMjAx96lOf8vUbN+JL+Bd4bhM7JyZBTlZWlmpr\nazvatbW1ys7OjsWpAPhAU1OTvve97ykjI0P33HMPAQ7QA35+vfRJkBMKhSztq6++WtXV1aqpqVFm\nZqY2b96sjRs39sWpAPShWOxdFV4DI9nrYCJtu9CvXz9Lu72mr7GxUd/73vc0fvx43XfffUpI+Md8\nifAaHOptAHTWZZCzcOFCbd++XSdPnlROTo5KS0vV0tIiSbr77rv13nvvaerUqfrwww+VkJCgRx99\nVJWVlUpJSdG6des0e/ZstbW16a677qLoGPCgeLhd/PjjjysjI8MW4ACIDpmcC+gu+5KRkWG5LdXZ\nnDlzNGfOnIsfGQDo/FIWl156KQEOgB5jxWMAnjZ06FC3hwDENT9ncvhqBAAAjEQmB0CfOnbsmO2x\n8AVD+/fvb+szYMAAtbS0KCkpSYFAwLbQH7uJAxeHTA4AX/JS4XFTU5MeeeQRvfbaa24PBYAhCHIA\nH/PKomRNTU36wQ9+oLS0NBUXF7s9HMAobNAJAC5pamrSj370I6WlpWnRokXMogLQZ6jJAdArFRUV\nXbYlqb6+3tKePn26JKmhoUErVqzQpEmTtHbtWkuAQw0OgN4iyAHgmrNnz6qoqEhf+tKXyOAAMeLF\n20hOIcgB4Jq0tDQtW7bM7WEAMBRBDuBjsdi7CoC3+DmTQ34Y8DEvTSEHgL5GJgfABVVVVXXb55VX\nXrG0T506ZeuTlZWlpqYm7dq1SyUlJcrMzLT1YQdxIDbI5ABADDU1Nen73/++Dh8+rFAo5PZwAPgE\nmRwAMdXc3Kzvf//7Sk1N1dKlS5lFBTiMTA4AxEBzc7OefvppAhwAriCTA/jY6tWruyw+PnnypKX9\n+9//3tanra3N0s7IyOj49zPPPKORI0fqRz/6kRITE3s3WAAXhUwOAF+K9d5VN998sz73uc8R4ABw\nBZkcADGTnJzs9hAA+BhBDgAABuN2FQD0UmtrK9PDAXgKmRwAkqRt27bZHnvzzTct7UhBTG5urpqa\nmvTDH/5tVRg1AAAeoklEQVRQ06ZN05w5c2x9hgwZYmmz8B/gHDI5AHypL/auag9wRowYoWuuuaYP\nRgUAfYNMDuBjvd27qrm5uSPAWbx4MevgAB5EJgcAeqi5uVmbNm0iwAHgWWRyAJ+qqKiwtHft2mXr\nE77Z5vjx4zv+ffbsWV1xxRV9ttBfXV2dpU3dDtA3yOQAQA+lpKRo3rx5LPQHwLPI5AAAYDAyOQB8\n6cc//rHbQwCAmCHIAXws2iCnpaVFFRUVOnfuXIxHBKCvBQIBx368httVgA+UlZXZHtu7d68k6Xe/\n+50ke5GxJKWmpqq5uVk///nPNXToUGVmZtpqcPqqQJhCYwB9jUwOgAvqHODMnz+fImMAcYVMDoCI\nWlpaLAEO6+AA8cmLt5GcwrsWgIh27txJgAMgrpHJAQxQVVVlae/cudPSPnLkiO2Yjz76SNddd53O\nnDkjScrKyrI8f9ttt+myyy6zBDj5+fm238MifoC3kckB4EvXXXfdBZ8bMGAAGRwAcY1MDgAABiOT\nA8DXWlpaWAMHgHHI5AAeF15vM2TIEFufgwcPWtrhNTgnTpywHTNx4kRJUlNTk5544gnNnj1bM2bM\nsPS59tprLe1I9TbhNTkAvIVMDgBfag9wRowYoU9+8pNuDwcA+pQngpzVq1dr9erVF3wu0tLR9Ke/\nX/pPnDjR8pOZmamHH344Yv/Nmzfr1ltvVWlpqeXnlVdesfVtamrSt771Lb311lvas2ePbrvtNt16\n66269dZbtXnz5qjH39V4vHg96U9/t/s7LdIYYvXjNYFQKBRy7eSBgFw8PRAXorldtXv3bkv7jTfe\nsLTDb1e1tLSovLxcb731lh555BElJCQoOzvb9nsv5nYVU8iBrjn52ddVEBYLq1ev9tTnuicyOQCc\nde7cORUUFEgS08QBw/k5k0PhMeBx7cFIu61bt9r6VFRUWNrtC/y1i5T9WbJkiZ5//nkVFxdLksaN\nG2frE01WhswNAK/iKxwAADASmRwAAAzmxdtITiGTAxiupaVFu3btirjY35IlS1wYEQA4g0wO4CHh\nM6kkac+ePZb2vn37bH2am5st7REjRnQ8vnnzZo0ZM0aTJ0+2FBmfPXtWt99+u86ePSuJ2hrAVGRy\nABinPcAZMmSIvvSlLzGLCoDvkMkBDNQ5wJk3bx4BDuBjZHIAGGX79u0EOAB8j0wOYKAbbrhBSUlJ\nBDgAfJ3JIcgBHBLNbt1/+ctfbI+F7zt1+vRpW5/Ro0db2hMmTLC0x44daztmypQpWt1p37hI4+ur\nYmS2fgDgBr7mAT5WWlrq9hAAxJift3UgyAHiXGtrq9ra2tweBgB4DkEOEMdaW1u1fft2267jAABq\ncgDHHDx40PZYMBi0tCMt9Hf8+HFLOzMzU9L5lYyff/55DR8+XPPnz7cUGd94442WY9zeaJMaHMA9\nXryN5BQyOUAcag9wBg8erE996lPMogKACMjkAHGmLwOcb3/72304MgBe5OdMDkEOEIdyc3NVVFTU\n6wxO+/RxADARQQ4QZ/r166fi4mK3hwEgTpDJAdDntm3bZmm/9dZbtj7hO4x/+OGHtj7hRbuzZs2y\n9UlOTu7yGADwI4IcAAAM5udMDlMyAA9rbW3Vnj17WOwPAC4CQQ7gUa2trXr55Zd19uzZmH0To/AY\nMJ+ft3XgdhXQB7Zu3Wp7LHwV4rffftvW58SJE5Z2amqqpPPTxHfu3KlLL71UX/ziFy2zqPLz822/\nZ8iQIRc17tLS0h4FOrHcxBMA+hpBDuAxLS0t2rZtm5KTk/Xxj3+chf4A9IoXMyxO4d0T8Jj9+/cT\n4ABAHyCTA3hMUVGRAoEAAQ4A9BJBDuAxiYmJbg8BgEH8fLuKIAe4CGVlZZZ2pKLiv/71r5Z2eJGx\nJGVkZFjaV155pa1P+GMFBQVRj7M7Pd276mKLjMMLlilWBuAEghzARa2trZKkpCR3XopMIQfM5+dM\nDjf9AZe0tbVp+/btOnz4sNtDAQAjkckBXNDW1qY333xTw4cP1/jx490eDgCD+TmTQ5ADdKOiosL2\nWGVlpaUdvtGmJJ08edLSzszMlHT+FtWf/vQnDR8+XIsXL7bMohozZozt90ycOLHHY/ZaDYzb5wfg\nTwQ5gIPaA5yUlBRde+21TBMHEHN+zuTwDgs4KBAIaPTo0br22ms98cZD4TEAkxHkAA5KTEzUhAkT\nPBHgSOf3rgJgNjboBHwivFbl2LFjtj7ha94cOHDA1if8sSNHjtj6hK+BM3r0aFuf8I01I/UJH2M0\n9S3hfaqqqmx9+nK9HQDwIjI5AADASAQ5QIy0trbqwIEDamtrc3soAHzMz7erug1yysrKVFBQoPz8\nfD344IO25z/44APdcsstKioq0sc+9jFLGr+7YwFTtbW1affu3WpqamIGFQC4pMt337a2Nt1zzz0q\nKytTZWWlNm7cqIMHD1r63H///Zo8ebLeeOMNPf3001q+fHnUxwImamtr09/+9jcNGjRIxcXFnvx2\n066ne1cBiD9+zuR0WXi8d+9e5eXlKTc3V5K0YMECbdmyxbI42cGDB7Vq1SpJ0oQJE1RTU6MTJ07o\n8OHD3R4LxFp3i+Lt3r3bdszf/vY3S/vVV1/t9vemp6dLOn+LqqKiQoMHD9bNN99syeJMnTrV9nvS\n0tIs7eTkZFuf8ALh8HNH40JFxt1NIffaooIA0BNdBjlHjx5VTk5ORzs7O9u2smtRUZF+85vfaNq0\nadq7d6/+/ve/68iRI1EdK1nfZEtKSlRSUnKRfwrgvpqaGg0cOFCFhYXcpgIgSSovL1d5eblr5/di\nhsUpXQY50VyYVatWafny5SouLlZhYaGKi4uVmJgY9UVlMTKYZOzYsZ5N2wJwR/gXeNanck6XQU5W\nVpZqa2s72rW1tcrOzrb0GTx4sH7yk590tMeMGaNx48apoaGh22MB05C9AeA1fv7S1WWQc/XVV6u6\nulo1NTXKzMzU5s2btXHjRkufM2fOaODAgerfv7+efPJJ3XDDDUpJSYnqWKAvRVOr8uyzz1ra+/bt\ns/UJv6164sQJW5/Bgwdb2uG1NdL5LwmdhS/8J0nvv/++pT1t2jRbn2j+rljVylCDAyCedfm1Mykp\nSevWrdPs2bM1adIk3XHHHZo4caKeeOIJPfHEE5LO78ZcWFiogoICvfjii3r00Ue7PBYwRVtbm1pb\nW90eRq9wuxgwn59nVwVCoVDItZMHAnLx9DBMNBmP8NlUF5vJGTRokA4cOKBhw4YpJydHY8aMsfW5\n4oorLO0pU6bY+pw5c8bSvuWWW2x9YpnJ4TUIOM/J110gENBjjz3myLkkadmyZZ56T2HvKqCH2tra\ndODAAQ0YMIA6MwCe58UMi1OokgR6oK2tTe+8844GDBig/Px8X795AIDXkcmBMSLdsvnlL39paf/l\nL3+xtF977TXbMcePH7e0U1NTJZ0PcCorKzVo0CDNmjXLMpPqyiuvtP2e8FtYwWDQ1if89lRFRYWt\nT/htrkg7iof/7SziBwAEOUDUEhISlJaWpvT0dKaKA4gbfs44E+QAUQoEAsrIyHB7GH2KvasAmIwg\nB/AxppAD5iOTA3hcpDqU8BqXSLvc79q1y9J+/fXXLe1jx47Zjhk2bJilnZmZaesTntGJ1Cd8gcBx\n48bZ+oT/XZF+T3h9zYU22+yMGhwAIMgBImpra1NdXR3BAoC45+dMDtWTQJj2aeLnzp2jwBgA4hiZ\nHKCTc+fO6Z133lH//v2Vk5Pj629AAMzg5/cxvqYC/+/cuXM6ceKErwIcCo8BmIy9qxAXIhUe79ix\nw9IOLzKWpAMHDlja4UW8Q4cO7fj3yZMn1dTUpOnTp1sCnNzcXNvvDa/VmTBhgq1PUlL3idLPfOYz\nlvbOnTttfSLtTN5XeA0CznN676r169c7ci5J+sIXvuCp9xRuVwH/b8SIEZL8ndoFAJMQ5AD/j+AG\nAMxCkAMAgMH8/AWOIAeeFF6bEr6InyTt3bvX0n7jjTdsfcIX+0tJSZF0vsg4FApp9OjRtmPGjh1r\naV922WW2Pu2bdrZramqy9Qn/3dXV1bY+27Zts7SjWQww0to94TVL0SwYCACmI8iB77TPohowYIDb\nQ3Ede1cB5vNzJocp5PCV9gAnKSlJQ4YMcXs4rmMKOQCTkcmBb4RCoY4AZ8SIEb7+dgPAP/z8Xkcm\nB74QCoX04YcfEuAAgI+QyYHjwgtpDx8+bOsTvtDf/v37bX2i2VF88ODBks4HOSkpKbryyistAU52\ndrbtmIEDB1rakRb1q6+vt7QHDRpk6/Puu+9a2gsWLLD1CS8YjnQLLZpNQik0BnAhfv5SRyYHvhAI\nBDRo0CBfv9gBwG8IcgAfo/AYMF8gEHDsx2sIcgAfKy0tdXsIABAz1OTAcbt377a0I9Xk7Nu3z9J+\n5ZVXbH3Ca3uSk5Mlna+/qa+v18CBA22ba4bXt7TvV9VVn86beIafq11ra6utT1pamqUdaZPR8Bqc\naOpvAKAnvJhhcQqZHBglFArp7Nmzkvz9wgYAEOTAIO0BTkJCAkXGAABuV8EMBDgAEJmf3w8JctCn\nKioqLO2GhgZbn3feecfSfumll2x93nzzTUv7+PHjtj6d61mampo0cOBA5eTkWF7Qw4cPtxwTXl8T\nviaOJLW0tFja7be/OguvyZk2bZqtz5kzZyxtL65lw95VAExGkAMj9O/f3xbgoHtMIQfM5+f3RWpy\nYASvrtEAAHAPmRwAAAzm5y+AZHIAAICRyOTgooUvxhdJeXm57bHwxQAjLfR3+vRpS3vAgAGSzs+i\namlpUSAQUGZmpqVPeFuyL/YXXjAcXmQs2TfbTElJsfUJBoOWdqQFDSMVIwOA08jkAHGgc4ATaWdw\n9ByFxwBMRpCDuBAe4Pj5m0lfYu8qwHxs0Al4HAEOAKCnyPnjor3xxhu2x15++WVLe+fOnbY+4TU4\n4fUtkpSQ8I/4OxQKKTU1VcnJyZYAJ7zeJisrq9sx9+vXz9LOycmx9bn00kst7UgbdIZv/Bm+0SYA\neIWfvxgS5MDzAoFAxOJfAAC6wu0qAAAQM0uXLlV6eroKCws7Hjt16pRmzpyp8ePHa9asWbYZtX2F\nIAfwMfauAsznduHxkiVLVFZWZnnsgQce0MyZM3Xo0CHNmDFDDzzwQEz+doIceFIoFHJ7CL7AFHIA\nsTZ9+nQNGzbM8thvf/tbLV68WJK0ePFiPffcczE5NzU5iNqmTZss7b/+9a+2PuGL/0Va6K+5ubnb\nc2VnZ1u+FURa6G/UqFHd/p7wWp7wF1ok+fn5lnb4AoKSN3cUB4BIYll4fODAAVVWVvb4uOPHjys9\nPV2SlJ6eruPHj/f10CQR5MCj/DwbAADixeWXX67LL7+8o/3LX/6yx78jlmvsEOQAAGAwL35pTE9P\n13vvvaeMjAzV1dVp5MiRMTkPNTkAAMBRN910kzZs2CBJ2rBhg+bPnx+T85DJgaTIm22GL+T3+uuv\nW9rhG21K9hqcaOpv2iP4zt82whf2i1RLM3DgQEs7fBE/yV63k5GRYWlHqvV5//33Le0pU6bY+lRV\nVVna8Vqjs3r1aoqPAcO5nclZuHChtm/frpMnTyonJ0dr1qzRqlWrdPvtt+upp55Sbm6ufvGLX8Tk\n3AQ5cJ3bL0A/Ky0tJcgBEFMbN26M+Pgf/vCHmJ+bIAcAAIP5+YskNTkAAMBIBDkAAMBI3K7yqUiF\nxuH+8pe/WNrhhcaRCo+jWak4PT3dkj4N39FbklJTU7tsS/bi5EsuucTWJ7zQOLzPxIkTbcdEs8hg\nvBYaA/AfblcBDvLzC85r2LsKgMnI5AA+xswqwHx+/mJJJgcAABiJTI5PhdedfOtb37L1Ca+5Ca/R\niUZ4TYwkjR492tKOtJz38OHDuzxGOr+JZ1dtyV7vEwwGLe1jx47ZjommJgcA4gWZHAAAAMOQyQEA\nwGBkcgD4EoXHAExGkAP4WGlpqdtDABBjgUDAsR+v4XaVD1RUVNge27Ztm6UdvuO4JL366qs9Plc0\nC/2lpaVZ2pEKfcMLjyPtMJ6ent7l75WkmpoaS/vGG2+09QEAmIkgB33Ki5E8APiZn9+XuV0FAACM\nRJADAACMxO0qHwivv5GkP/3pT5b2rl27bH1aWlq6/L2XXXaZ7bHwGpxIiwGG19tE2iQzJSWl23P1\n79/f0k5OTrb1mTJliu0x/AN7VwHm43YVAF9iCjkAk5HJAQDAYGRyAAAADEMmx0DhNTiR6m3CN9ts\namrq9vcOGTJEo0aN6vhWkJWVZesTvlZNpE0zw+t0LrnkEluf8LVzItXtsJEmAHSPTA4Qhc4BDgAA\nXkeQg6gR4JiHwmPAfH7e1oEgB/Ax9q4CYDJqcgAAMJgXMyxOIciJcz/96U9tj4UXGu/evdvW5/Tp\n05b2gAEDbH3Ci33DC41zcnJsx4Rvmhm+8J9kXzAwPz/f1mfIkCGWNkXGAICe4nYVAAAwEpkcSJJC\noZBaWlrUr18/X6c2AcA0fn5PJ5MDSef3qfLzC8Gv2LsKgMnI5MSZqqoqS3v//v22Pjt37rS033vv\nPVuf9kX7QqGQzpw5oxEjRqigoMAS6IQv2ldQUGBpR6rJGThwoKU9dOhQW5/wmpzw3wvnMIUcMJ+f\nv8CSyfGx9gAnISHBFuAAABDvyOT4WHNzsxITE5WSkkKAAwCG8vP7e7eZnLKyMhUUFCg/P18PPvig\n7fmHHnpIxcXFKi4uVmFhoZKSkjqmJ3d3LNx1ySWXEOAAAIzVZZDT1tame+65R2VlZaqsrNTGjRt1\n8OBBS5+vfe1r2rdvn/bt26e1a9eqpKREQ4cOjepYuI8ABwDM5udtHbq8XbV3717l5eV1FIouWLBA\nW7ZsibgjtCQ988wzWrhw4UUdC7u6ujrbY6+//rqlHb6buCQdP37c0s7Ly7P1KS4utrTHjRtn61NU\nVNTl+EaOHGl7LPz/Xxbx87bVq1dTfAzAWF0GOUePHrXMoMnOztaePXsi9q2vr9eLL76oxx57rEfH\ndn6DLSkpUUlJSU/GjyidO3dOra2t6t+/v9tDgYeUlpYS5AAxVl5ervLyctfO78UMi1O6DHJ6cmG2\nbt2qadOmdUwZjvZY3mBj79y5c6qrq1P//v07po4DAJwR/gWejXGd02WQk5WVpdra2o52bW2tsrOz\nI/bdtGlTx62qnh6L2AmFQqqrq1NSUpJSU1PdHg4AwGFkci7g6quvVnV1tWpqapSZmanNmzdr48aN\ntn5nzpzRjh079Mwzz/T4WPxDeA1OpELtzZs3W9qnTp2y9bn88sslnS8cr6ysVF5enj796U8rIeEf\ndebhC/JFqtsJ3yQzMzPT0qbeBgDgZV0GOUlJSVq3bp1mz56ttrY23XXXXZo4caKeeOIJSdLdd98t\nSXruuec0e/Zsy2q3FzoWzgiFQqqsrNQll1xiC3AAAPCDbhcDnDNnjubMmWN5rD24abd48WItXrw4\nqmPhjEAgoNGjR+vSSy8lwMEFsXcVYD5uV8FI4bebgHAU/gMwGUEOAAAGI5MDTwgv5L3//vttfYLB\noKU9efJkW59PfOITlvZll11m65OcnGxpR6qXorAYABDPKNYwQGtrq/785z/ro48+cnsoAACP8fO2\nDgQ5ca61tVU7duxQ//79lZKS4vZwAADwDIKcONbW1qYdO3YoOTlZ11xzjSejaHgbhceA+fycyQmE\nQqGQaycPBOTi6T3n2WeftbRfeOEFW5/2Bfmam5u1ceNGTZgwQcuXL7dME7/xxhstx0Ta6JN6G0i8\nBgE3OPm6CwQC2rFjhyPnkqTrr7/eU+8pFB7HqXfffVfDhw+3BTgAAHTmxQyLUwhy4lReXp7y8vII\ncAAAuAA+IQEAgJHI5AAAYDBuV8Fx27Ztsz325ptvWtpTp06VJDU1NenEiRPKycnRlClTLH3CF/WL\nhCJjXAh7VwEwGberPK6pqUnr1q3Tzp073R4KDMQUcsB8fp5CTpDjYe0BzogRI3THHXe4PRwAAOIK\nt6s8qrm5uSPAWbRoEbOoAAAXxYsZFqcQ5LhkyJAhtsf+5V/+RZIUCoW0YMECFRcX66GHHlJiYmJH\nH+prAACIDkGOBwUCAd1777266qqrLAEOAAA95edMDvdAPGrKlCkEOIg5Co8BmIwgB/Cx0tJSt4cA\nIMaYXQUAAGAYanJc0r6oXzAY1B133KHvfve7FBUDAPqcFzMsTiGT46JgMKh58+YpLS1NEyZMcHs4\nAAAYhSDHJe0BTm5urtavX0+RMQAAfYzbVS4gwIFXsHcVYD5uV8FRr776qvLz8wlw4DqmkAMwWSAU\nCoVcO3kgIBdPDwCA45z87AsEAnrllVccOZckTZ061VOf62RyAACAkajJAQDAYNTkIGaCwaB2797t\n9jAAAPAdgpwYap9FtWHDBreHAkRE4TFgPj9v60DhcYwwTRzxwOTXIOBVThceV1RUOHIu6fxq/l56\nT6EmJwYIcAAAXuHFDItTuF3Vx0KhkG677TYCHAAAXEYmp48FAgGtXbtWhYWFBDgAANf5OZNDkBMD\nV111ldtDAADA97hdBfgYe1cBMBmzq3opFAr5OhUIAOgZp2dXvf76646cSzp/J8NLn+tkcnohGAzq\nxhtv1L59+9weCgAACENNzkUKBoOaO3euxo4dq6KiIreHAwBARH6+20Am5yJ0DnDWr1+vhAQuIwAA\nXkMmp4cIcAAA8YRMDqK2f/9+TZo0iQAHRmDvKgAmY3YV4GO8BgHnOT27av/+/Y6cS5IKCws99Z5C\nKgIAABiJmhwAAAxGTQ4iCgaDKi8vd3sYAADgIhDkXED7LKqNGze6PRQAAHARuF0VQedp4o8//rjb\nwwFihr2rAPP5+XYVs6vCsA4OACCWnJ5ddeDAAUfOJUmXX365pz7XyeSEueOOOwhwAADG8HMmhyAn\nzHe/+10VFBQQ4AAAEOcIcsJMmjTJ7SEAANBn/JzJIV0BAACM5Osgx0vFUYAb2LsKMF8gEHDsx2t8\nG+QEg0HNmjVLu3fvdnsogGtKS0vdHgIAxIwva3KCwaDmzZun3NxcTZ061e3hAAAQM17MsDjFd5mc\nzgHO+vXrlZiY6PaQAABADPgqk0OAAwDwGzI5PlFdXa3LL7+cAAcAAB/wVSbnqquu0rp169weBuAZ\n7F0FwGTsXQUAgIOc3ruqurrakXNJUn5+vqc+1311uwoAAPiHsUFOMBhUWVmZ28MAAMBVLAZomPZZ\nVL/61a/cHgoAAHCJcYXHnaeJP/HEE24PBwAAV3kxw+IUozI5rIMD9Ax7VwEwmVGzq2655RYNHTqU\nAAeIEjMcAec5Pbvq7bffduRckjR27FhPvacYFeS89dZbGjNmDAEOECWCHMB5BDnOMaomJy8vz+0h\nAADgKdTkAAAAGCZug5xz5865PQQAAOBhcRnkBINB/dM//ZO2b9/u9lCAuMbeVYD5/LwYYNwVHgeD\nQc2dO1djx47V+vXrlZAQl3EaAMCnnC48rqmpceRckpSbm0vh8cUiwAEAoGe8mGFxStxECQQ4AACg\nJ+Imk/Puu++quLhYDz/8MAEOAABR8nMmJ+5qcgAAiGdO1+S8++67jpxLkkaPHu2pz3VSIoCPsXcV\nYD5mV7l1cjI5gKt4DQLOczqTU1tb68i5JCknJ8dT7ymeyOQ88sgjlnYwGNSzzz7r0mjMVF5e7vYQ\njMc1jj2ucexxjc3j50yOJ4Kc5557ruPf7bOonn/+eU9Fg/GON67Y4xrHHtc49rjGMImnZld1nib+\n5JNPejIqBAAgnvj5s9QTmRzpHwHOmDFjWAcHAAD0muuFxwAA+I2ThcdHjx515FySlJWV5alSE1dv\nV3npQgAAYCI/JxS4JwQAAIzkqcJjAADQt8jkAAAAGIZMDgAABiOTAwAAYBgyOQAAGIxMDgAAgGHI\n5AAAYDAyOQAAAIYhyAEAAEbidhUAAAbjdhUAAIBhyOQAAGAwMjkAAACGIZMDAIDByOQAAAAYhkwO\nAAAGI5MDAABgGDI5AAAYjEwOAACAYcjkAABgMDI5AAAAhiHIAQAARuJ2FQAABuN2FQAAgGHI5AAA\nYDAyOQAAAIYhkwMAgMHI5AAAABiGTA4AAAYjkwMAAGAYMjkAABiMTA4AAIBhCHIAAICRuF0FAIDB\nuF0FAABgGDI5AAAYjEwOAACAYcjkAABgMDI5AAAAhiGTAwCAwcjkAAAAGIZMDgAABiOTAwAAYBgy\nOQAAGIxMDgAAgGEIcgAAgJG4XQUAgMG4XQUAAGAYMjkAABiMTA4AAIBhCHIAADBYIBBw7CeSpUuX\nKj09XYWFhQ7/5QQ5AAAghpYsWaKysjJXzk1NDgAABnO7Jmf69Omqqalx5dwEOQAAGMzJICclJcWx\nc0WDIAcAAEOFQiG3h+AqanIAAICRCHIAAICRCHIAAEDMLFy4UB//+Md16NAh5eTk6Kc//alj5w6E\n/H7DDgAAGIlMDgAAMBJBDgAAMBJBDgAAMBJBDgAAMBJBDgAAMBJBDgAAMNL/Ae7sTVoHj5SAAAAA\nAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x49a5bd0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFtCAYAAAA0zQKQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wVfWdx/HPgcQliCAouUGghVJSCQnJtbj4RHvT9KKo\nZHBgHJ2Wpops12nXQVxH3I4ap1uM47AuYGd0W6rxoViXFmQr0kLluvWpVAnqWAvWJTbQmwhiRCQh\nD/e3f1wTHhJu7vM9v5P3ayaTcB/O+c7h4cPve77nXMcYYwQAAHJmSK4LAABgsCOMAQDIMcIYAIAc\nI4wBAMgxwhgAgBwjjAEAyLG4wri7u1t+v1/z5s2TJB06dEjBYFDFxcWaM2eOWltbM1okAABeFlcY\nr1q1SiUlJXIcR5JUV1enYDCoPXv2qKqqSnV1dRktEgAALxswjPft26fNmzfrpptuUs/9QTZt2qSa\nmhpJUk1NjTZu3JjZKgEA8LC8gV5w66236oEHHtDhw4d7H2tpaZHP55Mk+Xw+tbS09HlfzyoaAIDB\nItmbWsZcGf/mN79RYWGh/H7/aXfgOM5pg9cYw1cGv+65556c1+D1L44xx9krXxzjzH+lIubK+JVX\nXtGmTZu0efNmtbe36/Dhw1q0aJF8Pp+am5tVVFSkcDiswsLClIoAAGAwi7kyXrFihZqamrR37149\n/fTT+sY3vqEnnnhC1dXVqq+vlyTV19dr/vz5WSkWAAAvSug645529PLly7V161YVFxfrhRde0PLl\nyzNSHGILBAK5LsHzOMbZwXHOPI6xuzkm1Ub36TbsOCn30AEAsEUquccduAAAyDHCGACAHCOMAQDI\nMcLYIzo7JW4RDgB2Iow9YsMG6V/+JddVAACSQRh7xCefSEeO5LoKAEAyCGOPaGuTjh3LdRUAgGQQ\nxh7R1iZ1dOS6CgBAMghjj2BlDAD2Iow9gpUxANiLMPYIVsYAYC/C2CMIYwCwF2HsEbSpAcBehLFH\nsDIGAHsRxh5BGAOAvQhjj6BNDQD2Iow9gpUxANiLMPaInjA2JteVAAASRRh7xNGj0e9dXbmtAwCQ\nOMLYI9raot9pVQOAfQhjjyCMAcBehLFHtLVJZ53FRDUA2Igw9gBjpPZ2adQoVsYAYCPC2AOOHZPy\n8qSCAlbGAGAjwtgD2tqiQfwP/8DKGABsRBh7QE8Yn3EGYQwANiKMPeDElTFtagCwD2HsAbSpAcBu\nhLEH0KYGALsRxh5AmxoA7EYYewBtagCwG2HsAbSpAcBuhLEH0KYGALvFDOP29nbNmjVLFRUVKikp\n0Z133ilJqq2t1YQJE+T3++X3+7Vly5asFIv+0aYGALvlxXpy2LBh2r59u4YPH66uri5ddtlleuml\nl+Q4jpYtW6Zly5Zlq07EcGKbmpUxANgnZhhL0vDhwyVJHR0d6u7u1ujRoyVJxpjMVoa49YRxXh4r\nYwCw0YBhHIlEdMEFF+j999/XzTffrOnTp2v9+vVas2aNHn/8cc2cOVMrV67U2Wef3ee9tbW1vT8H\nAgEFAoF01o7PtbVJn/+fiTAGgCwJhUIKhUJp2ZZj4lzifvLJJ7r88stVV1enkpISjR07VpJ01113\nKRwOa+3atSdv2HFYPWfJv/2bdOaZUiQSbVP/6Ee5rggABp9Uci/uaepRo0bpqquu0uuvv67CwkI5\njiPHcXTTTTdpx44dSe0c6cEAFwDYLWYYHzx4UK2trZKktrY2bd26VX6/X83Nzb2v2bBhg8rKyjJb\nJWLiOmMAsFvMc8bhcFg1NTWKRCKKRCJatGiRqqqq9J3vfEe7du2S4ziaPHmyHnnkkWzVi370hLHE\nNDUA2ChmGJeVlWnnzp19Hn/88cczVhAS1xPGkQgrYwCw0YDT1HC/njDu7CSMAcBG3A7TA7gdJgDY\njTD2gKNHmaYGAJsRxh7A7TABwG6EsQdwnTEA2I0w9gCuMwYAuxHGHsAAFwDYjTD2ANrUAGA3wtgD\naFMDgN0IY8t1dkrGSPn5tKkBwFaEseV6VsWOQ5saAGxFGFvuxA+JoE0NAHYijC13YhjTpgYAOxHG\nljs1jFkZA4B9CGPLnRjGeXlSd3f0oxQBAPYgjC13Yhg7DvenBgAbEcaWOzGMJYa4AMBGhLHl2tqk\n4cOP/5ohLgCwD2FsuVNXxgxxAYB9CGPL0aYGAPsRxpbrb2VMmxoA7EIYW442NQDYjzC2HG1qALAf\nYWw52tQAYD/C2HK0qQHAfoSx5WhTA4D9CGPL0aYGAPsRxpZjZQwA9iOMLcfKGADsRxhb7uhRBrgA\nwHaEseVoUwOA/Qhjy9GmBgD7xQzj9vZ2zZo1SxUVFSopKdGdd94pSTp06JCCwaCKi4s1Z84ctba2\nZqVY9MV1xgBgv5hhPGzYMG3fvl27du3SW2+9pe3bt+ull15SXV2dgsGg9uzZo6qqKtXV1WWrXpyC\nNjUA2G/ANvXwzz+5vqOjQ93d3Ro9erQ2bdqkmpoaSVJNTY02btyY2SpxWrSpAcB+A4ZxJBJRRUWF\nfD6fKisrNX36dLW0tMjn80mSfD6fWlpaMl4o+kebGgDslzfQC4YMGaJdu3bpk08+0eWXX67t27ef\n9LzjOHIcp9/31tbW9v4cCAQUCARSKhZ90aYGgNwIhUIKhUJp2ZZjjDHxvvhHP/qRCgoK9LOf/Uyh\nUEhFRUUKh8OqrKzUX/7yl5M37DhKYNNIQne3lJcnRSJSz/+H/uu/pNdfj34HAGRPKrkXs0198ODB\n3knptrY2bd26VX6/X9XV1aqvr5ck1dfXa/78+UntHKlpb5eGDTsexBIrYwCwUcw2dTgcVk1NjSKR\niCKRiBYtWqSqqir5/X5de+21Wrt2rSZNmqRnnnkmW/XiBKe2qCXOGQOAjWKGcVlZmXbu3Nnn8TFj\nxmjbtm0ZKwrxaWuTPh9278U0NQDYhztwWay/lTFtagCwD2FssdO1qVkZA4BdCGOLcc4YALyBMLYY\nbWoA8AbC2GK0qQHAGwhji7EyBgBvIIwtxjljAPAGwthitKkBwBsIY4vRpgYAbyCMLUabGgC8gTC2\nGG1qAPAGwthisdrUfHolANiDMLZYf2E8dKg0ZEj0s44BAHYgjC3WXxhLnDcGANsQxhY7erT/MGai\nGgDsQhhbLNbKmCEuALAHYWyx04UxK2MAsAthbDHOGQOANxDGFqNNDQDeQBhbjDY1AHgDYWwx2tQA\n4A2EscVoUwOANxDGFmtrk4YP7/s4bWoAsAthbDHa1ADgDYSxpYyR2ttpUwOAFxDGljp2TMrPj34o\nxKloUwOAXQhjS52uRS2xMgYA2xDGlooVxqyMAcAuhLGlBloZE8YAYA/C2FK0qQHAOwhjS9GmBgDv\nIIwtRZsaALyDMLYUbWoA8I6YYdzU1KTKykpNnz5dpaWlWr16tSSptrZWEyZMkN/vl9/v15YtW7JS\nLI6jTQ0A3pEX68n8/Hw9+OCDqqio0JEjR/TVr35VwWBQjuNo2bJlWrZsWbbqxCloUwOAd8QM46Ki\nIhUVFUmSRowYoWnTpmn//v2SJGNM5qvDadGmBgDviPuccWNjoxoaGnTRRRdJktasWaPy8nItXrxY\nra2tGSsQ/aNNDQDeEXNl3OPIkSNauHChVq1apREjRujmm2/W3XffLUm66667dNttt2nt2rV93ldb\nW9v7cyAQUCAQSEvRYGUMALkWCoUUCoXSsi3HDNBv7uzs1NVXX625c+dq6dKlfZ5vbGzUvHnz9Pbb\nb5+8YcehlZ1B//7v0UD+8Y/7Prdxo/Too9Kzz2a/LgAYrFLJvZhtamOMFi9erJKSkpOCOBwO9/68\nYcMGlZWVJbVzJK+jI9qO7g8DXABgl5ht6pdffllPPvmkZsyYIb/fL0lasWKF1q1bp127dslxHE2e\nPFmPPPJIVorFcceOSaNH9/8cbWoAsEvMML7ssssUiUT6PD537tyMFYT4HDsWDd3+MMAFAHbhDlyW\nihXGtKkBwC6EsaUGCmPa1ABgD8LYUrSpAcA7CGNLMU0NAN5BGFuKNjUAeAdhbCna1ADgHYSxpZim\nBgDvIIwtNdDKmDY1ANiDMLbUQGHc2Sn1c78WAIALEcaWijVN7TjHAxkA4H6EsaVirYwlhrgAwCaE\nsaUGCmOGuADAHoSxpeIJY4a4AMAOhLGlaFMDgHcQxpaiTQ0A3kEYWyrWNLVEmxoAbEIYW6irK/o9\nL+/0r6FNDQD2IIwtNFCLWqJNDQA2IYwtFE8Yc0tMALAHYWwhVsYA4C2EsYXiDWNWxgBgB8LYQgNN\nUksMcAGATQhjC9GmBgBvIYwtRJsaALyFMLZQvNPUrIwBwA6EsYVoUwOAtxDGFqJNDQDeQhhbiGlq\nAPAWwthCtKkBwFsIYwtxO0wA8BbC2EKsjAHAWwhjCzHABQDeQhhbiOuMAcBbYoZxU1OTKisrNX36\ndJWWlmr16tWSpEOHDikYDKq4uFhz5sxRa2trVopFVDzT1LSpAcAeMcM4Pz9fDz74oN555x299tpr\n+slPfqJ3331XdXV1CgaD2rNnj6qqqlRXV5eteiHa1ADgNTHDuKioSBUVFZKkESNGaNq0adq/f782\nbdqkmpoaSVJNTY02btyY+UrRizY1AHhLXrwvbGxsVENDg2bNmqWWlhb5fD5Jks/nU0tLS7/vqa2t\n7f05EAgoEAikVCyijh2TxoyJ/Rra1ACQWaFQSKFQKC3biiuMjxw5ogULFmjVqlU666yzTnrOcRw5\njtPv+04MY6QPK2MAyL1TF5n33ntv0tsacJq6s7NTCxYs0KJFizR//nxJ0dVwc3OzJCkcDquwsDDp\nApC4eMK4oEBqb89OPQCA1MQMY2OMFi9erJKSEi1durT38erqatXX10uS6uvre0Ma2dHRMXAYn3mm\ndPRoduoBAKQmZpv65Zdf1pNPPqkZM2bI7/dLku677z4tX75c1157rdauXatJkybpmWeeyUqxiDp2\nbOBLm4YPlz77LDv1AABSEzOML7vsMkUikX6f27ZtW0YKwsDiaVMPH87KGABswR24LBRPGNOmBgB7\nEMYWineA6+hRyZjs1AQASB5hbKF4wjgvL/rF5U0A4H6EsYXimaaWaFUDgC0IYwvFM00tMVENALYg\njC0UT5taYqIaAGxBGFso3jCmTQ0AdiCMLZTIypg2NQC4H2FsoXgHuGhTA4AdCGML0aYGAG8hjC3E\nNDUAeAthbJlIROrsjD+MWRkDgPsRxpbp6IgGseMM/Fra1ABgB8LYMvGeL5ZYGQOALQhjy8Q7SS1x\nzhgAbEEYWyaRlTFtagCwA2FsmXgnqSXa1ABgC8LYMomeM6ZNDQDuRxhbhjY1AHgPYWwZpqkBwHsI\nY8swTQ0A3kMYW4aVMQB4D2FsmUSmqTlnDAB2IIwtwzQ1AHgPYWwZ2tQA4D2EsWW4tAkAvIcwtkwi\n09QFBVJbW/RjFwEA7kUYWyaRlfGQIdHXtrdntiYAQGoIY8skMk0t0aoGABsQxpZJZGUsMVENADYg\njC2TTBizMgYAdyOMLZNoGNOmBgD3I4wtk8g0tcTKGABsMGAY33jjjfL5fCorK+t9rLa2VhMmTJDf\n75ff79eWLVsyWiSO45wxAHjPgGF8ww039Albx3G0bNkyNTQ0qKGhQVdccUXGCsTJaFMDgPcMGMaz\nZ8/W6NGj+zxujMlIQYgt0UubaFMDgPvlJfvGNWvW6PHHH9fMmTO1cuVKnX322X1eU1tb2/tzIBBQ\nIBBIdnf4HG1qAHCHUCikUCiUlm05Jo4lbmNjo+bNm6e3335bkvThhx9q7NixkqS77rpL4XBYa9eu\nPXnDjsPqOQPmzZOWLJGqq+N7/a23Sl/4QvQ7ACBzUsm9pKapCwsL5TiOHMfRTTfdpB07diS1cySO\naWoA8J6kwjgcDvf+vGHDhpMmrZFZtKkBwHsGPGd8/fXX68UXX9TBgwc1ceJE3XvvvQqFQtq1a5cc\nx9HkyZP1yCOPZKNWKLlp6gMHMlcPACB1A4bxunXr+jx24403ZqQYDIxpagDwHu7AZRna1ADgPYSx\nZbjpBwB4D2FsGaapAcB7CGPL0KYGAO8hjC3D5xkDgPcQxpZJdJqac8YA4H6EsUWMoU0NAF5EGFuk\nq0saMkQaOjT+99CmBgD3I4wtkugktUSbGgBsQBhbJNEWtSQNGxZ9XySSmZoAAKkjjC2STBg7jlRQ\nwOoYANyMMLZIopPUPWhVA4C7EcYWSWZlLDHEBQBuRxhbJJUw5vImAHAvwtgiyUxTS7SpAcDtCGOL\n0KYGAG8ijC1CmxoAvIkwtgjT1ADgTYSxRWhTA4A3EcYWoU0NAN5EGFuEaWoA8CbC2CK0qQHAmwhj\ni9CmBgBvIowtwjQ1AHgTYWwR2tQA4E2EsUUIYwDwJsLYIqlMU3POGADcizC2CCtjAPAmwtgihDEA\neBNhbJFkp6m5tAkA3I0wtkiyK2MubQIAdyOMLUKbGgC8iTC2SLLT1LSpAcDdYobxjTfeKJ/Pp7Ky\nst7HDh06pGAwqOLiYs2ZM0etra0ZLxJRtKkBwJtihvENN9ygLVu2nPRYXV2dgsGg9uzZo6qqKtXV\n1WW0QByXbBifcYbU1RX9AgC4T8wwnj17tkaPHn3SY5s2bVJNTY0kqaamRhs3bsxcdThJsmHsOJw3\nBgA3y0v0DS0tLfL5fJIkn8+nlpaW0762tra29+dAIKBAIJBwgTgu2UubpOOt6pEj01sTAAxWoVBI\noVAoLdtKOIxP5DiOHMc57fMnhjFSl+zKWGJlDADpduoi89577016WwlPU/t8PjU3N0uSwuGwCgsL\nk945EpPsNLXERDUAuFnCYVxdXa36+npJUn19vebPn5/2otC/VFbGTFQDgHvFDOPrr79el1xyiXbv\n3q2JEyfq0Ucf1fLly7V161YVFxfrhRde0PLly7NV66BHmxoAvCnmOeN169b1+/i2bdsyUgxiSzWM\naVMDgDtxBy6LpGOaGgDgPoSxRWhTA4A3EcaWiESk7m4pPz+59xPGAOBehLElOjqiLeoYl3XHdOaZ\nnDMGALcijC2RSotaYmUMAG5GGFuCMAYA7yKMLZHKJLVEmxoA3IwwtgQrYwDwLsLYEqncl1oijAHA\nzQhjS6S6MqZNDQDuRRhbgjY1AHgXYWwJwhgAvIswtgTT1ADgXYSxJdKxMiaMAcCdCGNLpDpNPXas\n9OGHkjHpqwkAkB6EsSVSXRmPGiUNHSq1tqavJgBAehDGlkg1jCXpC1+Q/va39NQDAEgfwtgShDEA\neBdhbIlUp6klwhgA3IowtkS6VsYffJCeegAA6UMYWyLVaWqJlTEAuBVhbAnOGQOAdxHGlkhHGH/x\ni4QxALgRYWyJ9vbUw3jcOOnAAamzMz01AQDSgzC2xMGD0rnnpraNvDypqEjavz89NQEA0oMwtkRz\nc3RlmyrOGwOA+xDGlgiHo6vaVBHGAOA+hLElmpvTF8ZcawwA7kIYW6C9XTp6VBo9OvVtsTIGAPch\njC3Q0hJdFTtO6tvi8iYAcB/C2ALpOl8ssTIGADcijC2QrvPF0vEwNiY92wMApI4wtkA6w3jkSGno\nUKm1NT3bAwCkLi+VN0+aNEkjR47U0KFDlZ+frx07dqSrLpwgnWEsHV8dp2MgDACQupTC2HEchUIh\njRkzJl31oB/hsOT3p297PZc3lZenb5sAgOSl3KY2nHzMuHSvjJmoBgB3SXll/M1vflNDhw7V9773\nPS1ZsuSk52tra3t/DgQCCgQCqexu0MpUmxoAkLxQKKRQKJSWbTkmhaVtOBzWuHHjdODAAQWDQa1Z\ns0azZ8+ObthxWDWnyRe/KL34ojRpUnq2t26d9Oyz0tNPp2d7AIDUci+lNvW4zz+5YOzYsbrmmmsY\n4MoAY47f9CNdWBkDgLskHcZHjx7Vp59+Kkn67LPP9Lvf/U5lZWVpKwxRH38sFRRIw4alb5uEMQC4\nS9LnjFtaWnTNNddIkrq6uvStb31Lc+bMSVthiEr3+WIp+lGMBw5InZ1Sfn56tw0ASFzSYTx58mTt\n2rUrnbWgH5kI47y86Db37ZMmT07vtgEAieMOXC6XiTCWuLwJANyEMHa5cDjaVk43zhsDgHsQxi6X\nqZUxYQwA7kEYuxxhDADeRxi7HGEMAN5HGLtcOEwYA4DXEcYu19ycuQGuDz6QIpH0bxsAkBjC2MU6\nOqTDh6Vzzkn/tkeOjF7e9Prr6d82ACAxhLGLffihNHasNCRDv0tXXSU991xmtg0AiB9h7GKZOl/c\ngzAGAHcgjF0sU5PUPS65RPq//4uGPgAgdwhjF8vU8FaP/Hxpzhxp8+bM7QMAMDDC2MUyvTKWaFUD\ngBsQxi6WjTC+4grp97+Xjh3L7H4AAKdHGLtYpge4pOi09vTp0v/+b2b3AwA4PcLYxTJ9zrgHrWoA\nyC3C2MWy0aaWomH8m99IxmR+XwCAvghjlzImGsY+X+b3VV4utbdLe/Zkfl8AgL4IY5c6fFgaOlQa\nMSLz+3IcWtUAkEuEsUtlq0Xdo6dVDQDIPsLYpbI1vNWjqir6oRHNzdnbJwAgijB2qWyvjM88U7rl\nFumf/olBLgDINsLYpbJxjfGp7r5bamqSHn00u/sFgMGOMHapbK+MJemMM6QnnpDuuEPauze7+waA\nwYwwdqlchLEklZZKy5dLNTVSd3f29w8AgxFh7EKdndH7RV9wQW72v3Rp9HKn//zP3OwfAAabvFwX\ngL7++7+lL39ZqqjIzf6HDpUee0z6x3+MfubxxRfnpg4AGCxYGbuMMdLKldJtt+W2jsmTo+eP58+X\nXnstt7UAgNcRxi7z4ovSZ59JV16Z60qiH6/42GNSdTWBDACZRBi7zMqV0rJl0hCX/M7MnSvV1xPI\nAJBJjjGZucWD4zjK0KY96y9/kb7+damxUSooyHU1J9u8Wfrud6P/Wbj+eimPaQMAOEkqueeS9Rck\n6cEHpX/+Z/cFsRRtm2/YIP3sZ9JXvhL93tGR66oAwBsIY5c4cEB65hnp+9+P/z2hUChj9fTn0kuj\n57Qfe0xav16aMiX6H4jDh7NaRlZl+xgPVhznzOMYu1vSYbxlyxadf/75mjp1qu6///501jQorVkj\nLVwoFRbG/55c/eWaPVvaskX69a+lP/4xOnn9r/8q/e1vOSkno/gHLDs4zpnHMXa3pM78dXd36wc/\n+IG2bdum8ePH68ILL1R1dbWmTZuW7vo8raMjGmgPPSR98IH0wgu5rigxF14oPf109Bz36tXR66LL\ny6Vp06Kt7PPPj14vPX68NGxYrqsFAPdKKox37NihL3/5y5o0aZIk6brrrtOzzz5LGA+grU167z1p\n925p587olPJXvhKdnq6utncoatIk6T/+Q7rnHmnHjugg2u7d0v/8j/TXv0Y/9GLUKGnChGgwjx8v\nnXde9GvcOGnsWOmcc6Rzz5VGjoze/QsABpOkpqnXr1+v3/72t/rpT38qSXryySf1xz/+UWvWrDm+\nYf5FBQAMMslOUye1FosnaLmsCQCA+CQ1wDV+/Hg1NTX1/rqpqUkTJkxIW1EAAAwmSYXxzJkz9d57\n76mxsVEdHR365S9/qerq6nTXBgDAoJBUmzovL08PPfSQLr/8cnV3d2vx4sUMbwEAkKSkrzOeO3eu\ndu/erb/+9a+68847dejQIQWDQRUXF2vOnDlqbW3t856mpiZVVlZq+vTpKi0t1erVq1MqfrCI55ru\nW265RVOnTlV5ebkaGhqyXKE3DHScn3rqKZWXl2vGjBm69NJL9dZbb+WgSrvFe3+CP/3pT8rLy9Ov\nf/3rLFbnDfEc41AoJL/fr9LSUgUCgewW6AEDHeODBw/qiiuuUEVFhUpLS/XYY48NvFGTJrfffru5\n//77jTHG1NXVmTvuuKPPa8LhsGloaDDGGPPpp5+a4uJi8+c//zldJXhSV1eXmTJlitm7d6/p6Ogw\n5eXlfY7Zc889Z+bOnWuMMea1114zs2bNykWpVovnOL/yyiumtbXVGGPM888/z3FOUDzHuOd1lZWV\n5qqrrjLr16/PQaX2iucYf/zxx6akpMQ0NTUZY4w5cOBALkq1VjzH+J577jHLly83xkSP75gxY0xn\nZ2fM7abtdpibNm1STU2NJKmmpkYbN27s85qioiJVVFRIkkaMGKFp06bp73//e7pK8KQTr+nOz8/v\nvab7RCce+1mzZqm1tVUtLS25KNda8Rzniy++WKNGjZIUPc779u3LRanWiucYS9KaNWu0cOFCjR07\nNgdV2i2eY/yLX/xCCxYs6B26Pffcc3NRqrXiOcbjxo3T4c/vE3z48GGdc845yhvgRhJpC+OWlhb5\nfD5Jks/nGzAMGhsb1dDQoFmzZqWrBE/av3+/Jk6c2PvrCRMmaP/+/QO+hqBITDzH+URr167VlW74\n0GmLxPtn+dlnn9XNN98sifsVJCqeY/zee+/p0KFDqqys1MyZM/XEE09ku0yrxXOMlyxZonfeeUfn\nnXeeysvLtWrVqgG3m9AAVzAYVHNzc5/Hf/zjH5/0a8dxYv4lOnLkiBYuXKhVq1ZpxIgRiZQw6MT7\nj5E55bpu/hFLTCLHa/v27fr5z3+ul19+OYMVeU88x3jp0qWqq6vr/Si6U/9cI7Z4jnFnZ6d27typ\n3//+9zp69KguvvhiXXTRRZo6dWoWKrRfPMd4xYoVqqioUCgU0vvvv69gMKg333xTZ5111mnfk1AY\nb9269bTP+Xw+NTc3q6ioSOFwWIWn+cSDzs5OLViwQN/+9rc1f/78RHY/KMVzTfepr9m3b5/Gjx+f\ntRq9IN5r59966y0tWbJEW7Zs0ejRo7NZovXiOcZvvPGGrrvuOknRIZjnn39e+fn5XDoZp3iO8cSJ\nE3XuueeqoKBABQUF+trXvqY333yTMI5TPMf4lVde0Q9/+ENJ0pQpUzR58mTt3r1bM2fOPP2G03VS\n+/bbbzd1dXXGGGPuu+++fge4IpGIWbRokVm6dGm6dut5nZ2d5ktf+pLZu3evOXbs2IADXK+++iqD\nRUmI5zh/8MEHZsqUKebVV1/NUZV2i+cYn+i73/2u+dWvfpXFCu0XzzF+9913TVVVlenq6jKfffaZ\nKS0tNe+j04VxAAABDklEQVS8806OKrZPPMf41ltvNbW1tcYYY5qbm8348ePNRx99FHO7aQvjjz76\nyFRVVZmpU6eaYDBoPv74Y2OMMfv37zdXXnmlMcaYP/zhD8ZxHFNeXm4qKipMRUWFef7559NVgmdt\n3rzZFBcXmylTppgVK1YYY4x5+OGHzcMPP9z7mu9///tmypQpZsaMGeaNN97IValWG+g4L1682IwZ\nM6b3z+6FF16Yy3KtFM+f5R6EcXLiOcYPPPCAKSkpMaWlpWbVqlW5KtVaAx3jAwcOmKuvvtrMmDHD\nlJaWmqeeemrAbSb1QREAACB90jZNDQAAkkMYAwCQY4QxAAA5RhgDAJBjhDEAADlGGAMAkGP/D1P/\nanAh1r1LAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1e537ad0>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "idx = np.logical_and(np.logical_and(np.isfinite(dtm_rrmse_2k), dtm_rrmse_2k<1.5), np.isfinite(sfm_rrmse_2k))\n",
      "fig, cbar = viz.scatter_density(np.hstack([0.5, dtm_rrmse_2k[idx], 1.5]), np.hstack([0.5, sfm_rrmse_2k[idx], 1.5]), \n",
      "                                return_cbar=True, vmin=0, vmax=3, cmap=matplotlib.cm.gray_r)\n",
      "fig.axes[0].plot([0,100],[100,0], 'k--')\n",
      "fig.axes[0].plot([100,0],[50,50], 'k--')\n",
      "fig.axes[0].plot([50,50],[100,0], 'k--')\n",
      "cbar.set_ticks([0,1,2,3])\n",
      "cbar.set_ticklabels([1,10,100,1000])\n",
      "\n",
      "fig.set_size_inches([10,10])\n",
      "fig.savefig('figures/scatter_rrmse_2k.svg')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAJECAYAAAD0VzVyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VOWdx/HvJAFBQe5EctEhJJKAEVPQagXJqoAUFKjr\nJbaFgt3ai33pvlpfZXftStztCl5qL3hrY7vUbgGrLkjXRtflFay2ilIvLAGJl2iICCKKcpKQC7N/\nsBMzZw7JJJk558xzPu/XK696Zp4z58kpmfnN7/k9zxOKRCIRAQAAGCbD6w4AAACkAkEOAAAwEkEO\nAAAwEkEOAAAwEkEOAAAwEkEOAAAwEkEOAABImWXLlik7O1ulpaWdjx08eFCzZs3S6aefrtmzZ+vj\njz/ufO62225TUVGRiouL9dRTT3U+vm3bNpWWlqqoqEg33HBDQtcmyAEAACmzdOlSVVdXxzy2cuVK\nzZo1S7t379ZFF12klStXSpJqa2u1fv161dbWqrq6Wt/+9rcVXc7vW9/6lh588EHV1dWprq4u7jWd\nEOQAAICUmTFjhkaMGBHz2OOPP64lS5ZIkpYsWaINGzZIkjZu3KiKigoNGDBA4XBYhYWFeuGFF7R3\n7159+umnOueccyRJixcv7jynOwQ5AAAYauTIkQqFQq79DB06NKF+7du3T9nZ2ZKk7Oxs7du3T5L0\n3nvvKS8vr7NdXl6eGhsb4x7Pzc1VY2Njj9chyAEAwFAfffSRq9c7fPhwr8+JBkipQJADAABclZ2d\nrffff1+StHfvXo0dO1bSsQxNQ0NDZ7s9e/YoLy9Pubm52rNnT8zjubm5PV6HIAcAAIO5OVyVqMsu\nu0xr1qyRJK1Zs0YLFy7sfHzdunVqbW3V22+/rbq6Op1zzjk65ZRTdPLJJ+uFF15QJBLRQw891HlO\nd7L6dssAAAB6VlFRoS1btujAgQPKz8/XrbfequXLl+vKK6/Ugw8+qHA4rIcffliSNGnSJF155ZWa\nNGmSsrKydO+993YGT/fee6++9rWvqbm5WV/84hd1ySWX9HjtUCQ6NwsAABgllfUuTiKRiPwUVpDJ\nAQDAYG4HOX5CTQ4AADASmRwAAAzmZibHb8jkAAAAI5HJAQDAYGRyAAAADEMmBwAAg5HJAQAAMAyZ\nHAAADEYmBwAAwDAEOQAAwEgMVwEAYDCGqwAAAAxDJgcAAIORyQEAADAMmRwAAAxGJgcAAMAwZHIA\nADAYmRwAAADDkMkBAMBgZHIAAAAMQyYHAACDkckBAAAwDEEOAAAwEsNVAAAYjOEqAAAAw5DJAQDA\nYGRyAAAADEMmBwAAg5HJAQAAMAyZHAAADEYmBwAAwDBkcgAAMBiZHAAAAMMQ5AAAACMxXAUAgMEY\nrgIAADAMmRwAAAxGJieFqqurVVxcrKKiIq1atSru+dtvv12ZmZk64YQTlJubq3/9139N+Fwc09N9\nuvPOO1VWVqaysjKVlpYqKytLH3/8cULn4pie7tNHH32kRYsWacqUKfr85z+vHTt2JHwupGXLlik7\nO1ulpaWOz+/atUvnnXeeBg0apLvuuivmOe5vYnq6xxs3btSUKVNUVlamqVOnavPmzZ3PcY+RtiIp\n1N7eHpkwYULk7bffjrS2tkamTJkSqa2tjXk+JycnctFFF8U939O5OKa392nTpk2Riy66qE/nBlUi\n9+n73/9+5NZbb41EIpHIrl27uMe99Mwzz0T++te/Rs444wzH5/fv3x958cUXI//0T/8UufPOOzsf\n5/4mrqd7fPjw4c7/fu211yITJkyIRCLc43QnKTJ69GjXflIcVvRaSjM5W7duVWFhocLhsAYMGKCr\nr75aGzdujHk+NzdXgwcPjnu+p3NxTG/v0+9+9ztVVFT06dygSuQ+7dy5U3/zN38jSZo4caLq6+u1\nf/9+7nGCZsyYoREjRhz3+TFjxmjatGkaMGBAzOPc38T1dI9POumkzv8+fPiwRo8eLYl7jPSW0iCn\nsbFR+fn5ncd5eXlqbGyMeX7s2LH685//rClTpmjdunXavn17QufimN7cp6amJj355JO6/PLLe31u\nkCVyn6ZMmaLHHntM0rEPhXfeeUd79uzhHqcY9ze5NmzYoJKSEs2dO1c/+9nPJHGPTRAKhVz78ZuU\nBjk9/cKhUEijRo1SQ0ODXn31Vc2aNUtPPPFEKrtknN78o9q0aZOmT5+u4cOH9/rcIEvkPi1fvlwf\nf/yxysrKtHr1apWVlSkzM5N7nGLc3+RauHChdu7cqU2bNumrX/2qIpGI110C+iWls6tyc3PV0NDQ\nedzQ0KC8vLyY5/ft26cTTzxRkjRy5EiFQiEdPHhQeXl53Z6LY3q6x12tW7euc6iqt+cGWSL3aejQ\nofrVr37VeTx+/HhNmDBBzc3N3OMU4t9wasyYMUPt7e28FxsiyF8GUprJmTZtmurq6lRfX6/W1lat\nX79el112WczzO3fu1Ntvv63W1lb9+te/1uDBgzVy5Mgez8Uxid6nQ4cO6ZlnntGCBQt6fW7QJXKf\nDh06pNbWVknSL3/5S82cOVNDhgzhHieZPbPA/U2eN998s/P+/vWvf5UkjRo1inuMtJbSTE5WVpZW\nr16tOXPmqKOjQ9dee61KSkr0wAMPSJKuu+46zZs3T8XFxZKkU045RY8++mi35yJWIvdYOjbWPmfO\nHA0ePLjHcxErkXtcW1urr33tawqFQjrjjDP04IMPdnsuYlVUVGjLli06cOCA8vPzVVlZqba2NknH\n7u/777+vs88+W5988okyMjL005/+VLW1tRoyZAj3N0E93eNHH31Uv/nNbzRgwAANGTJE69atk8S/\nYRMEOZMTijDoCgCAkUKhkLKzs1273r59+3xVy8W2DgAAwEhs6wAAgMGCPFzVbSanp2XAa2pqNGzY\nsM4tA5y2ZMjNzdWYMWNYDhwAALiq20zO0qVL9d3vfleLFy8+bpuZM2fq8ccfj3mso6ND119/vZ58\n8knNnj1bw4YN07p167R48WJddtllnUVrQY4uAQDB5WbdSpA/a7vN5PS0DLjk/H9UdBnw/fv3q6io\nSIsXL9YTTzzhuBx4JBLRzJkzFYlE+Enhzy233OJ5H0z/4R5zj0344R6n/gfu6VfhcSgU6tyS4Ytf\n/KJqa2slfbYMePR/o8uAOy0HXl5erldeeUXl5eX6yU9+0p/uAADgOz/5yU9UXl6u8vJyXXDBBa5f\nP8jbOvSr8Phzn/ucGhoadOKJJ+qPf/yjFi5cqN27d3c+n8gvXFNTo/LyctXU1PSnKwAA+NKNN96o\nG2+8UZZlad68eV53J1D6lckZOnRo55YMc+fOVVtbW8wy4NEl16PLgB9vOfCFCxf2pxtIQHl5uddd\nMB73OPW4x6nHPU6NaIAzfvx4168d5ExOv4Kcrov+bN26VZFIJGZLhtGjR2v37t166KGHNHfu3OMu\nB37jjTf2pxtIAG9cqcc9Tj3ucepxj1PjF7/4hQoKCjpXQ4c7uh2u6mkZ8EceeUT33XefsrKydOKJ\nJ8YtAz5v3jy1tLTo008/1Ze+9CWWAwcABNINN9wgScrIcH8NXj9mWNzi6bYOoVCISnMAQKC4+dkX\nCoVc3TV+z549vvpcZ1sHIMBWrFjhdRcAIGXI5AABxt8gkHyWZam5uVmjR492fN7tTE5+fr4r15Kk\nhoYGX72nkMkBACBJLMvS/PnztXr1aq+7ArFBJwAASRENcMLhsH74wx963Z1OQS48JpMDAEA/dQ1w\nqqqqlJmZ6XWXIDI5AAD0S0tLi68DnCBncghygAC75ZZbvO4CkPZOOOEEff3rX9fVV1/tuwAn6Jhd\nBQCAi9yeXRUOh125liTV19f76nOdmhwAAGAkhqsAADBYkGtyyOQAAJAgy7L0d3/3d/rwww+97goS\nQJADAEACotPE29vbNXz4cK+7k7BQKOTaj98Q5AABxt5VQGJYByc9MbsKCDD+BoGeJTvAcXt2VUFB\ngSvXkqS33nrLV+8pFB4DANCNhx56KK0zOH4cRnILmRwgwPgbBHoWiUQUiUSUkZGcCg+3MzkTJkxw\n5VqS9Oabb/rqPYVMDgAA3fBrUW2i0rnv/UXhMQAAMBKZHCDA2LsKiGVZlg4dOqScnByvu5I0ZHIA\nBBJTyIHPRGdR3XPPPV53BUlCJgcAEHhdp4nfeuutXncnqcjkAAAQUCz0Zy4yOQCAwGptbTU+wAly\nJocgBwAQWAMGDND111+vhQsXGhngBB3DVUCAUXiMoAuFQrr88suNDnCCvEEnKx4DAcbfIOA+t1c8\nLi4uduVakrRr1y5fvaeQyQEAAEYiyAEABIJlWVq8eLH27dvndVdcFeThKoIcAIDxLMvSvHnzlJWV\npTFjxnjdHbiE2VUAAKNFA5yCggJVVVUlbTfxdOHHDItbgvX/NIAY7F0F0wU9wAk6MjlAgDGFHKZ7\n5JFHAh/gBDmTwxRyAIDRIpGIrz7o3Z5CPnnyZFeuJUk7duzw1ec6mRwAgNH8FOB4Ici/fzBzdwAA\nwHhkcgAARrAsSx988IHC4bDXXfEVMjkAAonCY5giOovqvvvu87or8BEKj4EA428QJki3aeJuFx6X\nlpa6ci1J2r59u6/eUxiuAgCkrXQLcLzAcBUAAGmmra2NAAfdIpMDAEhLAwYM0E033aS5c+cS4HQj\nyJkcghwAQNqaN2+e112AjxHkAAHG3lWA+YKcyWF2FQAALnJ7dtVZZ53lyrUk6ZVXXvHV5zqDmAAA\n37MsSxUVFdqzZ4/XXUk7oVDItR+/IcgBAPiaZVmaP3++Bg0apHHjxnndHaQRanIAAL4VDXDC4bCq\nqqqUmZnpdZfSjh8zLG4hkwMA8CUCHPQXQQ4QYOxdBT974oknCHCSIMg1OcyuAgKMv0H4XSQS8eWH\nZ3+4Pbtq6tSprlxLkrZt2+ar9xQyOQAA3zItwIG7KDwGAMBgQQ4UyeQAADxnWZZ2797tdTdgGIIc\nAICnorOo7r//fq+7YqQgFx4T5AABxt5V8FrXaeJ33HGH192BYZhdBQDwRFDXwXF7dtU555zjyrUk\naevWrb76XCeTAwBwXUdHhy699NLABThwF7OrAACuy8zM1M0336yZM2cS4KSYH2tl3EKQAwDwxIUX\nXuh1F2A4ghwAAAwW5EwONTlAgLF3FQCTEeQAAVZZWel1FxAAlmXp8ssv11tvveV1VxAwDFcBAFLG\nsizNmzdPBQUFCofDXncnkBiuAgAgyboGOFVVVcrI4CMH7iKTAwBIOgIc/yCTAwBAEm3evFmFhYUE\nOPAUmRwgwNi7Cqly6aWX6tJLL/W6G1CwMznsXQUAgIvc3rvq/PPPd+VakvTcc8/56nOdTA4AAAYL\nciaHgVIAQL9YlqXt27d73Q0gDkEOAKDPorOoqqqqvO4KjiMUCrn24zcEOQCAPuk6Tfzuu+/2ujtA\nHIIcIMDYuwp9xTo46SPImRxmVwEBxt8g+uLo0aOaNWuWTjvtNAKcPnB7dtUFF1zgyrUk6ZlnnvHV\newqzqwAAvZKRkaFbb71V5513HgEOfI0gBwDQa26uvYL+8eMwklsIwQEAgJHI5AAAYDAyOQACib2r\n0BPLsrRgwQK9/vrrXncF6DUyOUCAMYUc3bEsS/Pnz1c4HFZhYaHX3UEfkck5jmXLlik7O1ulpaWO\nz2/cuFFTpkxRWVmZpk6dqs2bN3c+V11dreLiYhUVFWnVqlXJ7TUAIKW6BjhVVVXKzMz0uktAr3W7\nTs6f/vQnDRkyRIsXL3bcl8SyLJ100kmSpO3bt2vRokV644031NHRoYkTJ+rpp59Wbm6uzj77bK1d\nu1YlJSWxF2eNDgDwHQKc1HJ7nZwLL7zQlWtJ0ubNm331ud5tJmfGjBkaMWLEcZ+PBjiSdPjwYY0e\nPVqStHXrVhUWFiocDmvAgAG6+uqrtXHjxiR1GQCQSn/5y19UWFhIgIO01++anA0bNugf/uEftHfv\nXj311FOSpMbGRuXn53e2ycvL0wsvvOB4fteagPLycpWXl/e3SwCAfrj44ot18cUXe90NY9TU1Kim\npsaz6we5JqffQc7ChQu1cOFC/elPf9JXv/pV7dq1q1fnU/gIeGfFihX8DQIpZv8CX1lZ6V1nAiZp\nU8hnzJih9vZ2HTx4UHl5eWpoaOh8rqGhQXl5ecm6FIAk4c0WgMn6FeS8+eabnQVGf/3rXyVJo0aN\n0rRp01RXV6f6+nq1trZq/fr1uuyyy/rfWwBAUlmWpZdeesnrbiCFgrwLebfDVRUVFdqyZYsOHDig\n/Px8VVZWqq2tTZJ03XXX6dFHH9VvfvMbDRgwQEOGDNG6deuOvWhWllavXq05c+aoo6ND1157bdzM\nKgCAt6KzqEpKSjRt2jSvuwMkXbdTyFN+caaQA57ibzC4mCbuHbenkM+aNcuVa0nSf//3f/vqPYVt\nHQAgYAhwEBRs6wAEGHtXBU8kEtGXvvQlApwA8WOtjFsYrgKAgHnppZdUVlZGgOMRt4erZs+e7cq1\nJOmpp57y1ec6mRwACBiKjIMlyJkcanIAAICRyOQAgMEikUigv8mDTA4AwECWZWnevHl67bXXvO4K\n4AmCHCDA2LfKXNFp4tnZ2Zo8ebLX3YGHgrziMUEOEGDsXWUm1sEBjiHIAQCDEOAAn6HwGAAM8sor\nr2jixIm65557CHAgKdiFxwQ5AGCQ888/X+eff77X3QB8gSAHAACDBTmTQ00OEGDsXQXAZAQ5QIAx\nhTy9WZalZ5991utuwOeYQg4ASCvRhf5++9vfet0VwLeoyQGANBMNcAoKCnTvvfd63R34nB8zLG4h\nkwMAaaRrgFNVVaWMDN7GgeMhkwMAaSISieiKK64gwEGvkMkBEEgUHqeXUCikVatWEeAACQpFIpGI\nZxcPheTh5YHA428QcJ+bf3ehUEgLFixw5VqStHHjRl+9p/BVAAAAGIkgBwB8yk/fiIF0RJADAD5k\nWZYuueQSvfTSS153BWmOxQABAL5hWZbmz5+vnJwclZWVed0dIG0xhRwIMPau8p9ogBMOh1VVVaXM\nzEyvu4Q058cMi1vI5AABxhRyfyHAAZKLTA4A+MTOnTtVUlKin//85wQ4SJogZ3IIcgDAJ6ZNm6Zp\n06Z53Q3AGAQ5AAAYLMiZHGpyAABAyixbtkzZ2dkqLS3tfOzgwYOaNWuWTj/9dM2ePVsff/xx53O3\n3XabioqKVFxcrKeeeqpf1ybIAQKMwmPvWJalp59+2utuIAC8Xidn6dKlqq6ujnls5cqVmjVrlnbv\n3q2LLrpIK1eulCTV1tZq/fr1qq2tVXV1tb797W/r6NGjff7dCXKAAKusrPS6C4EUnUW1fv16r7sC\npNyMGTM0YsSImMcef/xxLVmyRJK0ZMkSbdiwQdKxva8qKio0YMAAhcNhFRYWauvWrX2+NjU5AOCi\nrtPE77//fq+7A/TLBx98oA8++KDX5+3bt0/Z2dmSpOzsbO3bt0+S9N577+ncc8/tbJeXl6fGxsY+\n948gBwBcwjo48EIqC4/Hjh2rsWPHdh7X1tb2+jV62hKiP/1nuAoAXHLNNdcQ4AA6lr15//33JUl7\n9+7tDJRyc3PV0NDQ2W7Pnj3Kzc3t83UIcgDAJbfffjsBDlzndeGxk8suu0xr1qyRJK1Zs0YLFy7s\nfHzdunVqbW3V22+/rbq6Op1zzjl9/t0ZrgICjL2r3DVx4kSvuwC4rqKiQlu2bNGBAweUn5+vW2+9\nVcuXL9eVV16pBx98UOFwWA8//LAkadKkSbryyis1adIkZWVl6d577+3XcFUoEolEkvWL9PrioZA8\nvDwAAK5z87MvFArpqquucuVakrR+/Xpffa4zXAUAKdCftT0AJAdBDgAkmWVZmj17tp577jmvuwL4\nsibHLdTkAAi8vXv3xhyPGzeuz6/VdZp41/U+ALiPIAcAkoR1cOBHfsywuIXhKiDA2LsqeQhwAP8h\nyAECjL2rkuftt99WaWkpAQ58h5ocAAiw/tTgRJ1xxhn62c9+loTeAEgWMjkAAMBIZHIAADCYH4eR\n3EImBwB6ybIsbdq0yetuAOgBQQ4QYOxd1XuWZWnevHnauHGj110BEhLkwmOCHCDAmELeO9EAp6Cg\nQL/4xS+87g6AHlCTAwAJ6BrgVFVVKSOD74hID37MsLiFv1IASMDixYsJcIA0QyYHABJw55136rTT\nTiPAQdoJciaHIAcAEjB+/HivuwCgl/hKAgQYhceA+ZhdBSCQ2LvKWUdHh9ddAJAEBDkA0IVlWbr4\n4ou1efNmr7sCoJ+oyQGA/xedJj5+/HjNnDnT6+4ASeHHYSS3kMkBAMUGOFVVVcrMzPS6SwD6iUwO\ngMAjwIHJyOQACCT2rjqmsbFRn/vc5whwAMOEIpFIxLOLh0Ly8PIAALjOzc++UCikZcuWuXItSfrV\nr37lq891MjkAAMBI1OQAAGAwanIAICAsy9Lvf/97r7sBwAUEOQACw7IszZ8/X0888YSv6gaAVGJb\nBwCBFKS9q6IBTjgcVlVVlS/fkAEkF0EOEGBB2bvKHuAwTRxBQiYHAAz29a9/nQAHCCBmVwEw3h13\n3KFx48YR4AABQ5ADwHh5eXledwHwjB+HkdzCcBUAADASmRwgwEzcu6q9vV1ZWby1AVFkcgAEkmlT\nyC3L0qxZs/Rf//VfXncFgA/wdQeAEbpOE7/kkku87g7gG0HO5BDkAEgLe/fujXts3LhxklgHB4Az\nghwAaY0AB+hekDM51OQASGsffPCBzj77bAIcAHG6DXKWLVum7OxslZaWOj6/a9cunXfeeRo0aJDu\nuuuumOeqq6tVXFysoqIirVq1Knk9BpA0JhQeh8Nh3X777QQ4wHGwrcNxLF26VNXV1cd9ftSoUfr5\nz3+u73//+zGPd3R06Prrr1d1dbVqa2u1du1a7dy5Mzk9BpA0Qdm7CkAwdVuTM2PGDNXX1x/3+TFj\nxmjMmDFx0zW3bt2qwsJChcNhSdLVV1+tjRs3qqSkpN8dBhBM0SJjAEhUSgqPGxsblZ+f33mcl5en\nF154wbFt13R5eXm5ysvLU9ElAAawLEuPPPKIlixZ4nVXgITV1NSopqbGs+v7cRjJLSkJcnpzQ02o\nCQCQepZlad68eSooKNDixYsD/caN9GL/As8wsXtSEuTk5uaqoaGh87ihoYEN8gD0WdcAp6qqigAH\n6IUg/70kJciJRCIxx9OmTVNdXZ3q6+uVk5Oj9evXa+3atcm4FIAk8tPeVfbF/rou9Nc1wMnIYOUL\nAInpNsipqKjQli1bdODAAeXn56uyslJtbW2SpOuuu07vv/++zj77bH3yySfKyMjQT3/6U9XW1mrI\nkCFavXq15syZo46ODl177bUUHQM+lA7Dxd/61rcIcIB+CHImJxSxp2HcvHgoFJcFAhBMx8vk7Nu3\nT2PGjCHAgTHc/OwLhUL67ne/68q1JOnnP/+5rz7X2dYBgK9lZ2d73QUgrQU5k8NXIwAAYCQyOQB8\nYdy4cTpy5IgGDhwY6G+eQLIF+e+JTA4QYH4qPLYsS3PmzNFjjz3mdVcAGIIgBwgwvyxK1nWa+KJF\ni7zuDmAUNugEAI+wDg6AVKEmB0C/HG/qdyJtCHAApBJBDgDPHDp0SOeff77+5V/+hQAHSBE/DiO5\nhSAHgGdycnL0ox/9yOtuADAUQQ4QYH7auwpAagQ5k0N+GAgwP00hB4BkI5MDoF+cCo2dNDU16ZFH\nHtFXv/pVx+cTKWAG0HtkcgAghZqamrR48WK9/PLLvtq8D4DZyOQASCnLsrR48WLl5+frzjvvZBYV\n4DIyOQCQApZlaf78+Z0BTmZmptddAhAgZHKAAFuxYkWvio/tdTNS97UzN954o8LhsKqqqnoMcKjB\nAVIjyJmcUMTDAfJQKMT4POCh3v4N9jbI+fDDDzV8+HAyOEAXbn72hUIh3XTTTa5cS5LuuOMOX32u\nk8kBkDKjRo3yugsAAowgBwAAgwV5uIrCYwBJ0dLS4qs0NQCQyQGQsOPtMB5dB+eaa67Rd77zHQ96\nBuB4yOQACKRk7F0VDXDy8/O1YMGCJPQKAJKDTA4QYP3du8q+0B+zqAD/IZMDAL3EQn8A/I5MDoCE\nbdq0qfO/P/nkE40bN04PPfQQAQ7gY2RyAKCXTj75ZFVUVBDgAPAtMjkAABiMTA6AQOpv4TEA+BlB\nDhBglZWVCbWzLEs//vGPdfTo0RT3CECyhUIh1378huEqAJKcN9+UYtfBmTNnDjU4ANIGmRwAx9U1\nwGGaOIB0QyYHgCMCHMAMfhxGcguZHACOfvSjHxHgAEhrZHIAA+3atavHNpZl6Rvf+Ia2bdsmScrJ\nyYl5fvny5ZowYUJMgOO0QScAfyOTAyCQvvGNbxz3uaFDh5LBAZDWyOQAAGAwMjkAAq2lpUUdHR1e\ndwMAkopMDpBmnNaz6alWZtiwYcd9rqmpSX//93+vBQsW6PLLL+/V6wLwPzI5AAKpqalJ3/zmN5Wf\nn69FixZ53R0ASCpfBDkrVqw47h46K1ascFw6mva0D2r7nJwc3XXXXd22Lykpiflxat/U1KTy8nK9\n+OKLeuyxxzR58uTO9qtXr/bN70t72pvW3m1OfUjVj9+EIpFIxLOLh0Ly8PJAWkpkuMo+hdw+XNXU\n1KSKigq9+OKL2rFjhzIynL/vFBcX97O3AOzc/OzrLghLhRUrVvjqc90XmRwA7mpvb9f06dMl6bgB\nDgAzBDmTQ+ExkGacioHt2Z2dO3fGHJeUlMSdM3/+fN1999167733JEkXXnhhEnsJAN7jKxwAADAS\nmRwAAAzmx2Ekt5DJAQzX0tKiBx54QO3t7XHPLV682IMeAYA7yOQAPpLIzCmnNtG6mqhTTz1VktTc\n3Kwbb7xRpaWlKikpidmLatiwYVq5cmXnsdOmnsyuAtIfmRwAxokGODk5OaqqqmKzTQCBQyYHMFDX\nAOfmm28mwAECjEwOAKPcf//9BDgAAo9MDmCgb37zmxo4cCABDoBAZ3IIcgCf27ZtW8zxoUOH4toM\nHDgw5ji6mnGUU7GyJN1111363ve+J4ki43TSl53ogSBiuAoIsONt9AnAHEHe1oEgB0hzLS0tamtr\n87obAOCiDDKgAAAfRUlEQVQ7BDlAGmtubtY//uM/6g9/+IPXXQEA36EmB7Cx1zv0tdahL6/jVG9T\nV1cXc5yXlyfpWIDzwx/+UDk5OVq2bFlMkbG9jmfq1KnHvaYbtRzJuqc4hvuH3vDjMJJbyOQAaai5\nuVk/+MEPdMopp+gHP/gBs6gAwAGZHCDNJDPAueWWW5LYMwB+FORMDkEOkIbOO+88XXnllf3O4KxY\nsSI5HQIAHyLIAdLM4MGDVVFR4XU3AKQJMjkAks6+M7hTsah95+/t27fHtRk7dmzMcWtra1wb+0J+\nx1v8zysUygLwAkEOAAAGC3Imh9lVgI+1tLTo3//93x2zNwCA7hHkAD7V0tKiVatWaf/+/SmbIk7h\nMWC+IG/rEIpEIhHPLh4KycPLA31mr3lxWsTvwIEDMcdO2Rj7Y9FzWlpa9OMf/1ijR4/WnXfe2WOQ\nk5OT021/jrf5Jn+DgPvc/LsLhUK64447XLmWJN10002+ek+hJgfwma4Bzte//nUW+gPQL37MsLiF\n4SrAZzZt2tQZ4GRk8CcKAH1FJgfwmQULFigrK4sABwD6iSAH8JmBAwd63QUABgnycBVBDpAElmXF\nPfbBBx/EHDc2Nsa1KSwsjDkePXp0XBv7DuLV1dU9tkl08b2ue1fZFyaUpGHDhvXpdQHADwhyAA+1\ntLQoFArphBNO8OT6TCEHzBfkTA6D/oBHjhw5on/+53/WE0884XVXAMBIZHIADxw5ckT333+/CgoK\ntGDBAq+7A8BgQc7kEOQg0LZt2xb3WE8L60lSXV1dj23efffdmOPc3FxJxwKce+65R+PGjdOXv/xl\nHTx4sLONvQbGqY9jxozpsY39d7BvFirF1/E4XdvOaeNP6nQA+BVBDuCiI0eO6O6779aYMWO0dOlS\npokDSLkgZ3J4hwVclJGRoalTp/omwKHwGIDJvH+XBQJkwIABmjVrli8CHEmqrKz0ugsAUizIG3Qy\nXIVAO+mkk+Ies9evNDc3x7Wxb7750UcfxbWxr3lz9OjRuDaJ7EtlX28nHA73eE5f6mScznGqwQGA\ndOGPr5MAAABJRpADpMiRI0f02GOPqbW11euuAAiwIA9X9RjkVFdXq7i4WEVFRVq1alXc8x999JEW\nLVqkKVOm6POf/7x27NiR8LmAqY4cOaLVq1frk08+UVYWo8IA4IVug5yOjg5df/31qq6uVm1trdau\nXaudO3fGtPm3f/s3fe5zn9Orr76q3/zmN7rhhhsSPhcw0ZEjR/Tggw9q1KhRWrx4sW+KjJ103bsK\ngJmCnMnp9ivm1q1bVVhY2FnoePXVV2vjxo0qKSnpbLNz504tX75ckjRx4kTV19dr//79evPNN3s8\nF/CavYBYii/0/fjjj+PaHD58OOa4qalJktTa2qrf/va3Gj58uL7yla9I+qzgeMCAAXGv09LSEnPs\ntIeVfaHBIUOGxLWxFzknWjDcdQp5Igv9OW3imUiRs/284uLihPoHAP3RbZDT2Nio/Pz8zuO8vDy9\n8MILMW2mTJmixx57TNOnT9fWrVv1zjvvaM+ePQmdK8W+yZaXl6u8vLyPvwrgveeee07Dhw/XwoUL\nfZ3BAeCempoa1dTUeHZ9P2ZY3NJtkJPIjVm+fLluuOEGlZWVqbS0VGVlZcrMzEz4prIYGUxywQUX\nKBQKEeAA6GT/As/6VO7pNsjJzc1VQ0ND53FDQ4Py8vJi2gwdOlS/+tWvOo/Hjx+vCRMmqLm5ucdz\nAdMksu4NALiJTM5xTJs2TXV1daqvr1dOTo7Wr1+vtWvXxrQ5dOiQBg8erIEDB+qXv/ylZs6cqSFD\nhiR0LtAf9hoSp9qQZ599NubYPp17+/btcefYa2A+/fTTuDb22plBgwb12GbPnj1xbaZPnx5z7FTz\nMmrUqLjH7OwLGNo323S6N/ZrOW3Qab9/EyZM6LEvThKpwempjoiNQHuHzVSBHmZXZWVlafXq1Zoz\nZ44mTZqkq666SiUlJXrggQf0wAMPSJJqa2tVWlqq4uJiPfnkk/rpT3/a7bmAKVpbW3XkyBGvu9Ev\nq1ev9roLAFKM2VXdmDt3rubOnRvz2HXXXdf53+edd55ef/31hM8FTNDa2qqHH35YBQUF+sIXvuB1\nd/rsnnvu0fXXX+91NwAgJVilDOiltrY2bdiwQcOGDdO5557rdXcAoFt+zLC4hSkgQC+0tbXpscce\n07BhwzRv3jxmUQGAj5HJQVpwKsi1e/XVV3tsYy/+dSoqti/+F13Er2uAM3369JhznXYYty8YOHLk\nyLg2b775Zsyx0wxEex+ddjwvKiqKObYXT9uvI8UvIOhUlGp/rK+7km/bti3meOrUqT1eqy8otv1M\nUH9voCuCHCBBmZmZKikp0Zlnnhm3UjEA+FWQh6sIcoAEZWRk6KyzzvK6G0n1ne98x+suAEDKEOQA\nAcbMKsB8ZHIAn3PaSNNeU+LUxv6Yvd7Gqb7FvmCg00J/9s02Tz755Lg2iWy+2XV/N0l655134tqc\neeaZcY/Z2X9P+6J9Tgv99aXexqnOI5FFGZ1qcJIhkWsHBfcCiMfUEMBBW1ubnnvuubiABwDSTZAX\nAyTIAWza2tr0xz/+Ua2trcrKItkJAOmKd3Cgi2iAM3ToUM2cOZN1cACkPT9mWNzCOzjw/9ra2vT4\n448HKsBZsWKF110AgJQhkwNf2rx5c8yxZVlxberr62OO9+/fH9fmjTfeiDlubm6OOe66iN+OHTs0\naNAgnXvuuWpra+t8/IMPPoh7XXsxstO1CwoKYo77upmnfYdxe9GzJA0ZMiTm2L74n9Pu4du2bVNl\nZaUuvfRSSVJOTk5cG/uigsniVORs/z0TKVamuPYz3AscT5AzOQQ5wP8rKSmRpEBkcAAgCAhygP9H\ncAMAZiHIAQDAYAxXAS6y12M4bR5pX59mx44dcW0++eSTmGN7TYcUv7llJBKRJLW3t+vo0aOOU8QH\nDx7c7bETp4X+7HU79roZSfrwww9jjp3qbey1Mm+99VZcmzFjxnTbv5deeinusREjRkj6rE7JqabD\nfk+Li4vj2tg333TSl9qeZ599NuZ4+vTpcW3sG7c69Y9F8oDgIshB4LS3t+svf/mLxowZo9NPP93r\n7nhq6dKlXncBQIoFOZNDEQICJRrgnHTSSZo4caLX3fHcsmXLvO4CAKQMmRwERtcAp6ysTKFQqHP4\nCgBMRSYHMFxHR4f+/Oc/xwQ4AACzkclBSjkt+mYvgnVa6G/Pnj0xx++++25cG/uO4vZj6bPC40gk\nohEjRigcDqupqanz+YEDB8adY8/uOLWxFxVnZmbGtbEvEOj0e9p3L3cqcnYqqLazF/baj+1F2pJU\nVFQUc2xfgFFyXiDQri+L9tkLhiXppJNO6vZ1nf4t2XdXT2Qn9UQ49c+pqBlIB0H+UkcmB4EQCoWU\nnZ0d6D92AAgaghwgwFavXu11FwCkWCgUcu3HbwhygAC75557vO4CAKQMNTlIKaeF/lpaWmKOnRa3\ns2+sefDgwbg2H330UcxxtNamo6NDDQ0NysvLc1zsryunWhp7vY0T+xYQTnU7XTf/TPR1nPpjr+Vx\nqq+x1+DYFx7Mz8+POyda6xP9X/sCjE7XTmRjTacNTe2LFTrV+tjrdrxcxI/6G5jEjxkWt5DJgVE6\nOjr0xhtvKBKJOAYMAIDgIMiBMaIBzsCBAxUOhwP97QUAwHAVDEGAAwDOgvx+SJCDPnOqz7DX4NTX\n18e1sa9B0tjYGNdm3759McdOdR6HDx/u/O+DBw8qKytLOTk5amtr63zcXl9j3wDTqWanvb095the\nNyPF1+A41bPY61DstUhS/No+fanjkeJreYYPH95j/yTpy1/+cudzTnVFTvVSdvb/b8LhcFwbe22P\n04ad9pqbvtTgOP2bZENOILgIcmCEESNGaODAgYH+xtIXX/nKV7zuAoAUC/L7IjU5MIJf12gAAHiH\nTA4AAAYL8hdAghyknY6OjoRqVwAAwUaQgz5zKh61F346FR7bi1mdNt+0F7M2NzdLOlaYe+DAAZ1w\nwgkaOXJkTBundXHswZD9G41T4fEJJ5wQc+xU6Bvd+DPKXujr1MbJqFGjYo67FlNHjRgxIubYXjzt\nxL7BaXZ2dlwb+wKL9vspSUOHDo05dtos1N6/7du3x7Wx/55O/3bsC/Bt27Yt5ti+gacUv0FnX3m5\n8CCQakHO5FCTg7QRDXAyMzPjdu9G32zYsMHrLgBAyhDkIC10DXBGjhwZ6G8myUSQA5iPDToBHyPA\nAQD0BTU56DOn+gx7HUVdXV1cm4aGhm6PpWPFxVGRSEQDBw7UwIEDYxbUs29C6RT82Otr7AveOS2A\nZ6/BsS8O6PS6Tgv9nXjiiTHHBQUFcW3sNUv2GhgpfsHAROqI7HU7TvUs0T5H++nUxn5/nBbbmzBh\nQsyxU12RfTFC+0KJUvy/nalTp/Z4bRYMBHoW5C+GBDnwvVAoFPchDgBATxiuAgAAKbNs2TJlZ2er\ntLS087GDBw9q1qxZOv300zV79uy4jHWyEOQAAXbFFVd43QUAKeZ14fHSpUtVXV0d89jKlSs1a9Ys\n7d69WxdddJFWrlyZkt+dIAe+cvToUR06dIjF/lxCkAMg1WbMmBG3ntbjjz+uJUuWSJKWLFmSspme\n1OTguOwFmi+99FLM8bPPPht3zu7du2OOnQqP7a/b1NQk6ViBcVtbm0KhkDIyMmK+FTgVxdoLcJ3a\n2NfTGTx4cMyx00J/9kJj+zlSfGGv0y7f9sLj/fv39/g6Tv2xv7ZTkbNTH7tyWnzPvkDfO++8E9dm\n9OjRMcf2Xd2dXtteEO7UP3tRsZsoMkbQpLLweMeOHaqtre31efv27etcpDQ7O1v79u1LdtckEeTA\nJ7oGOFlZWYGeDQAA6WLy5MmaPHly5/Hvf//7Xr9GKtfYIciB5whwACB1/Piemp2drffff1+nnHKK\n9u7dq7Fjx6bkOtTkwHMdHR0EOAAQIJdddpnWrFkjSVqzZo0WLlyYkuuQyYEk5wXSXn311Zhje03O\nrl274s6xP+a0YKB9sbhoYNPW1tb5mL32w2lTSnsNjlM9i/28SCTS4zn2+hHLsuLa2F/HqR6o6+8j\nOW/Yad9U1H6OUx+d7oW9UNt+j5026Bw2bJh+8Ytf6Bvf+IYk6cwzz4xrk5OTE3Ps9P9nIptkJmMj\nTadaGjbWBHrm9ZfHiooKbdmyRQcOHFB+fr5uvfVWLV++XFdeeaUefPBBhcNhPfzwwym5NkEOPOf1\nH2CQdQ1yACAV1q5d6/j4008/nfJrE+QAAGCwIH+RpCYHrrIP8wAAkCoEOXBNJBJxXOMFAIBUYLgq\noOwFm/YiYyl+sb9XXnkl5thpASj7rtqJZG6cNt+0F/IOHz48ro19oT+n4l97Ya/9WvbXkOIXGXQq\nBrb/Xk7Bm/1aTilje8Fw193Xj8epP/bVRO3FyU4LbV1yySWSPluYz6mQ3KnQ2M5e7Otm8S+FxkDP\nGK4CEEi33HKL110AgJQhkwME2IoVK7zuAoAUI5MDAABgGDI5AWWvt3n99dfj2tjrdHbs2BFzbK+/\nSZS9VsW+UaQUX2PitDGkfSNIp9oee23KwIEDY46bm5vjzrFfy6lux16T49TGvrHmmDFj4trY2fsn\nSUeOHIk5ti/QJzn/7l051StVV1f3eO2SkpKYY2pggPRDJgcAAMAwZHIAADAYmRwAgfTb3/7W6y4A\nQMoQ5AAB9h//8R9edwFAioVCIdd+/IbhqgDYvHlz3GN1dXUxx88//3xcm//93/+NOX7nnXd6fW2n\nYtaRI0fGHNsLiKX4XaudXse+W7jT7tz2glz74oBOiwzaF+hzWtDQ/jskskCfvRBZit9hPJFFD52K\njNvb22OOL7roophjp4Lhbdu2SfqsINqpoJlCYwDpjCAHAACD+THD4haGqwAAgJEIcgAAgJEYrjKQ\nfaPF1157La6NfbPN7du3x7V59913e31te33N6NGj49rY607s9TdSfL2Nvb5Fiq+dsW+sKcXXr5x4\n4olxbezsfXaqtxk6dGiPr2Ov5cnNzY1rc/jw4ZjjcDjc4+s6sb+OfWNNp402c3Jy9L3vfa+zFof6\nG8BMDFcBCKTvfe97XncBAFKGTA4AAAYjkwMAAGAYMjkGeumll2KO7RttOrXpS/2N06aU9roO+3oy\nUny9jVNNjn1DSad1cux1MU71NvbNNjs6Orrti9O18vPz49o0NTXFHNvrgyTnNW/szjzzzB7b2O+F\n0/3qqZ5m7969vT4HgBnI5AAAABiGIAcIsBUrVnjdBQApFuRtHQhygACrrKz0ugsAkDLU5AAAYDA/\nZljcQpCT5qKbLHb15z//Oeb45ZdfjmuzZ8+eHl/bvrjeKaec0u2xFF8k61R8ay8QHjVqVI+v41Qg\nbC/AdfpDtm+KaT8nujllV5988knMsb3IWIovRnYqjLabMGFC3GOHDh2KOS4uLo5rY1/csS8FwxQZ\nAwgihqsAAICRyOQAAGCwIA9XkclB3B5LCI5bbrnF6y4AQMqQyUkz9vqMp59+Oq6NvQbHaaG/aGAT\niUR09OhRDRw4UCeddFJMxB/duDHKvnHl2LFj417XvpGm06aZ9pqcU089Na7Nhx9+GHPstKigfTHC\nTz/9NK7NhRdeGHP81ltv9Xjt5ubmmGOnTUYty4o5tt8rKbE6GHtNjhOnOp1kYQo5YD4yOQikaIAj\nKS7AAQAg3ZHJCbBoNicjI4MABwAMFeT39x4zOdXV1SouLlZRUZFWrVoV9/ydd96psrIylZWVqbS0\nVFlZWfr4448TOhfeysjIIMABABir2yCno6ND119/vaqrq1VbW6u1a9dq586dMW2+//3v6+WXX9bL\nL7+s2267TeXl5Ro+fHhC58J7BDgAYLYgb+vQ7XDV1q1bVVhYqHA4LEm6+uqrtXHjRpWUlDi2/93v\nfqeKioo+nYvEPPnkkzHHzz//fFyb999/P+bY6R9eXl5ezLHTTtv2glt7UbHTbtj2xf9OOOGEuDb2\n85za5Obmxhw7LdpnZ1/4T4pfVNBp4UE7+++dysLfVL52IlasWEHxMQBjdRvkNDY2xnz45eXl6YUX\nXnBs29TUpCeffFL33ntvr87t+gZbXl6u8vLy3vQfCYoWGWdmZnrdFfhIZWUlQQ6QYjU1NaqpqfHs\n+n7MsLil2yCnNzdm06ZNmj59euc3+UTP5Q029SKRiD799FNlZmbGZTYAAKll/wLPxrju6TbIyc3N\nVUNDQ+dxQ0ND3DBH1Lp16zqHqnp7LlLn6NGj+vTTT5WRkRG3Pg0AwHxkco5j2rRpqqurU319vXJy\ncrR+/XqtXbs2rt2hQ4f0zDPP6He/+12vz8XxPfvss3GPvfHGGzHHe/fujWszaNAgSccCnD179mjE\niBEqKCiI+Ydu31zTacG7yZMnO75ulFMtjb0mx2ljzWidVlR9fX1cm7a2th5fZ+rUqXGP9SSRGhin\newoASD/dBjlZWVlavXq15syZo46ODl177bUqKSnRAw88IEm67rrrJEkbNmzQnDlzYj6Ijncu3BGJ\nRLRnzx4NGDAgLsABACAIelwMcO7cuZo7d27MY9HgJmrJkiVasmRJQufCHaFQSKNHj9bgwYMJcHBc\n7F0FmC/InwGseGwwanDQEwr/AZiMIAcAAIORyYEvbN68OeZ4x44dcW0aGxtjjp2yNaeddlrMsb3I\nWJKys7Njjp1247ZPN7fXVB04cCDunOnTp8cc23dNl+KLf71eEM8ukd3DAQD+xy7kBujo6ND27dvV\n1NTkdVcAAD4T5G0dCHLSXEdHh1577TVlZWU5TrMGACCoCHLSWHt7u1577TUNGjRIxcXFvoyi4W8U\nHgPmC3ImJxSJRCKeXTwUkoeX953f//73McdbtmyJa3P48GFJxxbL+5//+R+NGjVKF198ccw/rksu\nuSTmnObm5rjXsdfgOC2sZ18Uj1oV8/A3CLjPzb+7UCikZ555xpVrSdIFF1zgq/cUCo/T1P79+3Xy\nySfHBTgAAHQV5M8Igpw0lZubq9zc3ED/4wUAoDvU5AAAACORyQEAwGBBzvgT5HjEaafrrKzY/zsK\nCwslSUeOHNGBAweUm5sbt/jfpEmT4l7HPpW8L7t1SxQaBwF7VwEwGcNVPnfkyBFVVVXp+eef97or\nMBBTyAHzBXkKOUGOj0UDnJEjR2rRokVedwcAgLTCcJVPtbS0dAY4V111lTIyiEcBAL3nxwyLWwhy\nPOJU79Le3i5JikQiuu2221RSUqKbb75ZmZmZnW1ycnJ6fB0AAECQ40uhUEhXXXWVZs2aFRPgAADQ\nW0HO5DAG4lOnn346AQ5SjsJjACYjyAECrLKy0usuAEgxZlcBAAAYhpqcFHBa6M8uWjBsWZauuuoq\n3X777briiitS3TUAQMD4McPiFjI5HrIsS/Pnz9eYMWM0ceJEr7sDAIBRCHI8Eg1wwuGwqqqqKDIG\nACDJGK7yQFNTk6655hoCHHiOvasA8wV5uIogJwV6WqBvy5YtKioq0n333UeAA08xhRyAyUKRSCTi\n2cVDIXl4eQAAXOfmZ18oFNKLL77oyrUk6eyzz/bV5zo1OQAAwEgMVwEAYLAg1+SQyUkxy7L0/PPP\ne90NAAAChyAnhaLTxNesWeN1VwBHFB4D5gvytg4UHqcI6+AgHZj8Nwj4lduFx9u2bXPlWpI0depU\nX72nUJOTAgQ4AAC/8GOGxS0MVyVZJBLR3/7t3xLgAADgMTI5SRYKhXTbbbeptLSUAAcA4LkgZ3II\nclLgrLPO8roLAAAEHsNVQICxdxUAkzG7qp8ikUigU4EAgN5xe3bVK6+84sq1pGMjGX76XCeT0w+W\nZemSSy7Ryy+/7HVXAACADTU5fWRZlubNm6eCggJNmTLF6+4AAOAoyKMNZHL6oGuAU1VVpYwMbiMA\nAH5DJqeXCHAAAOmETA4Stn37dk2aNIkAB0Zg7yoAJmN2FRBg/A0C7nN7dtX27dtduZYklZaW+uo9\nhVQEAAAwEjU5AAAYjJocOLIsSzU1NV53AwAA9AFBznFEZ1GtXbvW664AAIA+YLjKQddp4vfdd5/X\n3QFShr2rAPMFebiK2VU2rIMDAEglt2dX7dixw5VrSdLkyZN99blOJsfmqquuIsABABgjyJkcghyb\n22+/XcXFxQQ4AACkOYIcm0mTJnndBQAAkibImRzSFQAAwEiBDnL8VBwFeIG9qwDzhUIh1378JrBB\njmVZmj17tp5//nmvuwJ4prKy0usuAEDKBLImx7IszZ8/X+FwWGeffbbX3QEAIGX8mGFxS+AyOV0D\nnKqqKmVmZnrdJQAAkAKByuQQ4AAAgoZMTkDU1dVp8uTJBDgAAARAoDI5Z511llavXu11NwDfYO8q\nACZj7yoAAFzk9t5VdXV1rlxLkoqKinz1uR6o4SoAABAcxgY5lmWpurra624AAOApFgM0THQW1SOP\nPOJ1VwAAgEeMKzzuOk38gQce8Lo7AAB4yo8ZFrcYlclhHRygd9i7CoDJjJpdtWjRIg0fPpwAB0gQ\nMxwB97k9u+qtt95y5VqSVFBQ4Kv3FKOCnDfeeEPjx48nwAESRJADuI8gxz1G1eQUFhZ63QUAAHyF\nmhwAAADDpG2Qc/ToUa+7AAAAfCwtgxzLsnTxxRdry5YtXncFSGvsXQWYL8iLAaZd4bFlWZo3b54K\nCgpUVVWljIy0jNMAAAHlduFxfX29K9eSpHA4TOFxXxHgAADQO37MsLglbaIEAhwAANAbaZPJeffd\nd1VWVqa77rqLAAcAgAQFOZOTdjU5AACkM7drct59911XriVJp556qq8+10mJAAHG3lWA+Zhd5dXF\nyeQAnuJvEHCf25mchoYGV64lSfn5+b56T/FFJucnP/lJzLFlWfrP//xPj3pjppqaGq+7YDzucepx\nj1OPe2yeIGdyfBHkbNiwofO/o7Oo/vCHP/gqGkx3vHGlHvc49bjHqcc9hkl8Nbuq6zTxX/7yl76M\nCgEASCdB/iz1RSZH+izAGT9+POvgAACAfvO88BgAgKBxs/C4sbHRlWtJUm5urq9KTTwdrvLTjQAA\nwERBTigwJgQAAIzkq8JjAACQXGRyAAAADEMmBwAAg5HJAQAAMAyZHAAADEYmBwAAwDBkcgAAMBiZ\nHAAAAMMQ5AAAACMxXAUAgMEYrgIAADAMmRwAAAxGJgcAAMAwZHIAADAYmRwAAADDkMkBAMBgZHIA\nAAAMQyYHAACDkckBAAAwDJkcAAAMRiYHAADAMAQ5AADASAxXAQBgMIarAAAADEMmBwAAg5HJAQAA\nMAyZHAAADEYmBwAAwDBkcgAAMBiZHAAAAMOQyQEAwGBkcgAAAAxDkAMAAIzEcBUAAAZjuAoAAMAw\nZHIAADAYmRwAAADDkMkBAMBgZHIAAAAMQyYHAACDkckBAAAwDJkcAAAMRiYHAADAMGRyAAAwGJkc\nAAAAwxDkAAAAIzFcBQCAwRiuAgAAMAyZHAAADEYmBwAAwDAEOQAAGCwUCrn242TZsmXKzs5WaWmp\ny785QQ4AAEihpUuXqrq62pNrU5MDAIDBvK7JmTFjhurr6z25NkEOAAAGczPIGTJkiGvXSgRBDgAA\nhopEIl53wVPU5AAAACMR5AAAACMR5AAAgJSpqKjQF77wBe3evVv5+fn69a9/7dq1Q5GgD9gBAAAj\nkckBAABGIsgBAABGIsgBAABGIsgBAABGIsgBAABGIsgBAABG+j+EZ7icleRMjgAAAABJRU5ErkJg\ngg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x24fce4d0>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "idx = np.logical_and(np.logical_and(np.isfinite(dtm_rrmse_4k), dtm_rrmse_4k<1.5), np.isfinite(sfm_rrmse_4k))\n",
      "fig, cbar = viz.scatter_density(np.hstack([0.5, dtm_rrmse_4k[idx], 1.5]), np.hstack([0.5, sfm_rrmse_4k[idx], 1.5]), \n",
      "                                return_cbar=True, vmin=0, vmax=3, cmap=matplotlib.cm.gray_r)\n",
      "fig.axes[0].plot([0,100],[100,0], 'k--')\n",
      "fig.axes[0].plot([100,0],[50,50], 'k--')\n",
      "fig.axes[0].plot([50,50],[100,0], 'k--')\n",
      "cbar.set_ticks([0,1,2,3])\n",
      "cbar.set_ticklabels([1,10,100,1000])\n",
      "\n",
      "fig.set_size_inches([10,10])\n",
      "fig.savefig('figures/scatter_rrmse_4k.svg')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAJECAYAAAD0VzVyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X90FPW9//HX5gciAvLTYH7U5Uc0AXMxArZeoaRFQAoi\n/qiStkLB3uNtr1XP0V7tue2R8L2tv29/UX8Ge9CeAmpbkHtttL2eYG0tFH9BCUhU0BD5KUJ1SEiy\n2e8f3E2zs0OySXZnZj/zfJyTAzM7s/PJwO6+9z3veX9C0Wg0KgAAAMNkeT0AAACAdCDIAQAARiLI\nAQAARiLIAQAARiLIAQAARiLIAQAARiLIAQAAabN06VLl5eWprKysY92RI0c0c+ZMnXvuuZo1a5aO\nHj3a8djdd9+t4uJilZSU6MUXX+xY/9prr6msrEzFxcW65ZZbkjo2QQ4AAEibJUuWqKamJm7dPffc\no5kzZ2rXrl2aMWOG7rnnHklSXV2d1q5dq7q6OtXU1Ohb3/qWYu38vvnNb2rlypWqr69XfX19wnM6\nIcgBAABpM23aNA0dOjRu3XPPPafFixdLkhYvXqx169ZJktavX6/Kykrl5uYqHA5r3Lhx2rRpk/bt\n26dPPvlEF110kSRp0aJFHft0hSAHAABDDRs2TKFQyLWfQYMGJTWuAwcOKC8vT5KUl5enAwcOSJI+\n/PBDFRYWdmxXWFioxsbGhPUFBQVqbGzs9jgEOQAAGOrjjz929Xiffvppj/eJBUjpQJADAABclZeX\np/3790uS9u3bp7POOkvSyQxNQ0NDx3Z79+5VYWGhCgoKtHfv3rj1BQUF3R6HIAcAAIO5ebkqWfPn\nz9eqVaskSatWrdKCBQs61q9Zs0YtLS3avXu36uvrddFFF2nUqFEaPHiwNm3apGg0qqeeeqpjn67k\n9O6UAQAAdK+yslIbN27U4cOHVVRUpOXLl+vOO+/Utddeq5UrVyocDuvpp5+WJI0fP17XXnutxo8f\nr5ycHD300EMdwdNDDz2kr3/962pqatKXvvQlXXbZZd0eOxSN3ZsFAACMks56FyfRaFR+CivI5AAA\nYDC3gxw/oSYHAAAYiUwOAAAGczOT4zdkcgAAgJHI5AAAYDAyOQAAAIYhkwMAgMHI5AAAABiGTA4A\nAAYjkwMAAGAYghwAAGAkLlcBAGAwLlcBAAAYhkwOAAAGI5MDAABgGDI5AAAYjEwOAACAYcjkAABg\nMDI5AAAAhiGTAwCAwcjkAAAAGIZMDgAABiOTAwAAYBiCHAAAYCQuVwEAYDAuVwEAABiGTA4AAAYj\nkwMAAGAYMjkAABiMTA4AAIBhyOQAAGAwMjkAAACGIZMDAIDByOQAAAAYhiAHAAAYictVAAAYjMtV\nAAAAhiGTAwCAwcjkpFFNTY1KSkpUXFyse++9N+Hx++67T9nZ2TrttNNUUFCg//zP/0x6X5zU3Xl6\n4IEHVF5ervLycpWVlSknJ0dHjx5Nal+c1N15+vjjj3XllVdq4sSJ+uxnP6vt27cnvS+kpUuXKi8v\nT2VlZY6P79y5UxdffLH69++vBx98MO4xzm9yujvH69ev18SJE1VeXq5JkybppZde6niMc4yMFU2j\ntra26NixY6O7d++OtrS0RCdOnBitq6uLezw/Pz86Y8aMhMe72xcn9fQ8bdiwITpjxoxe7RtUyZyn\n22+/Pbp8+fJoNBqN7ty5k3PcQy+//HL09ddfj55//vmOjx88eDD617/+Nfof//Ef0QceeKBjPec3\ned2d408//bTj71u3bo2OHTs2Go1yjjOdpOiIESNc+0lzWNFjac3kbN68WePGjVM4HFZubq4WLlyo\n9evXxz1eUFCg008/PeHx7vbFST09T7/61a9UWVnZq32DKpnztGPHDn3hC1+QJJ133nnas2ePDh48\nyDlO0rRp0zR06NBTPj5y5EhNnjxZubm5ces5v8nr7hyfccYZHX//9NNPNWLECEmcY2S2tAY5jY2N\nKioq6lguLCxUY2Nj3ONnnXWW/vznP2vixIlas2aNtm3bltS+OKkn5+n48eN64YUXdPXVV/d43yBL\n5jxNnDhRv/nNbySd/FB4//33tXfvXs5xmnF+U2vdunUqLS3VnDlz9NOf/lQS59gEoVDItR+/SWuQ\n090vHAqFNHz4cDU0NOitt97SzJkz9fzzz6dzSMbpyX+qDRs2aOrUqRoyZEiP9w2yZM7TnXfeqaNH\nj6q8vFwrVqxQeXm5srOzOcdpxvlNrQULFmjHjh3asGGDrr/+ekWjUa+HBPRJWu+uKigoUENDQ8dy\nQ0ODCgsL4x4/cOCABgwYIEkaNmyYQqGQjhw5osLCwi73xUndnePO1qxZ03Gpqqf7Blky52nQoEF6\n4oknOpZHjx6tsWPHqqmpiXOcRvwfTo9p06apra2N92JDBPnLQFozOZMnT1Z9fb327NmjlpYWrV27\nVvPnz497fMeOHdq9e7daWlr0i1/8QqeffrqGDRvW7b44KdnzdOzYMb388su64oorerxv0CVzno4d\nO6aWlhZJ0uOPP67p06dr4MCBnOMUs2cWOL+p8+6773ac39dff12SNHz4cM4xMlpaMzk5OTlasWKF\nZs+erUgkohtuuEGlpaV69NFHJUk33nij5s6dq5KSEknSqFGj9Otf/7rLfREvmXMsnbzWPnv2bJ1+\n+und7ot4yZzjuro6ff3rX1coFNL555+vlStXdrkv4lVWVmrjxo06fPiwioqKVFVVpdbWVkknz+/+\n/fs1ZcoU/f3vf1dWVpZ+8pOfqK6uTgMHDuT8Jqm7c/zrX/9aTz75pHJzczVw4ECtWbNGEv+HTRDk\nTE4oykVXAACMFAqFlJeX59rxDhw44KtaLqZ1AAAARmJaBwAADBbky1VdZnK6awNeW1urM888s2PK\nAKcpGQoKCjRy5EjagQMAAFd1mclZsmSJvv3tb2vRokWn3Gb69Ol67rnn4tZFIhHddNNNeuGFFzRr\n1iydeeaZWrNmjRYtWqT58+d3FK0FOboEAASXm3UrQf6s7TKT010bcMn5HyrWBvzgwYMqLi7WokWL\n9Pzzzzu2A49Go5o+fbqi0Sg/afy56667PB+D6T+cY86xCT+c4/T/wD19KjwOhUIdUzJ86UtfUl1d\nnaR/tAGP/RlrA+7UDryiokJvvvmmKioq9OMf/7gvwwEAwHd+/OMfq6KiQhUVFfr85z/v+vGDPK1D\nnwqPL7zwQjU0NGjAgAH63e9+pwULFmjXrl0djyfzC9fW1qqiokK1tbV9GQoAAL5066236tZbb5Vl\nWZo7d67XwwmUPmVyBg0a1DElw5w5c9Ta2hrXBjzWcj3WBvxU7cAXLFjQl2EgCRUVFV4PwXic4/Tj\nHKcf5zg9YgHO6NGjXT92kDM5fQpyOjf92bx5s6LRaNyUDCNGjNCuXbv01FNPac6cOadsB37rrbf2\nZRhIAm9c6cc5Tj/OcfpxjtPjscce05gxYzq6ocMdXV6u6q4N+LPPPquHH35YOTk5GjBgQEIb8Llz\n56q5uVmffPKJrrrqKtqBAwAC6ZZbbpEkZWW534PXjxkWt3g6rUMoFKLSHAAQKG5+9oVCIVdnjd+7\nd6+vPteZ1gEIsGXLlnk9BABIGzI5QIDxGgRSz7IsNTU1acSIEY6Pu53JKSoqcuVYktTQ0OCr9xQy\nOQAApIhlWZo3b55WrFjh9VAgJugEACAlYgFOOBzW97//fa+H0yHIhcdkcgAA6KPOAU51dbWys7O9\nHhJEJgcAgD5pbm72dYAT5EwOQQ4QYHfddZfXQwAy3mmnnaZvfOMbWrhwoe8CnKDj7ioAAFzk9t1V\n4XDYlWNJ0p49e3z1uU5NDgAAMBKXqwAAMFiQa3LI5AAAkCTLsvQv//Iv+uijj7weCpJAkAMAQBJi\nt4m3tbVpyJAhXg8naaFQyLUfvyHIAQKMuauA5NAHJzNxdxUQYLwGge6lOsBx++6qMWPGuHIsSXrv\nvfd89Z5C4TEAAF146qmnMjqD48fLSG4hkwMEGK9BoHvRaFTRaFRZWamp8HA7kzN27FhXjiVJ7777\nrq/eU8jkAADQBb8W1SYrk8feVxQeAwAAI5HJAQKMuauAeJZl6dixY8rPz/d6KClDJgdAIHELOfAP\nsbuofv7zn3s9FKQImRwAQOB1vk18+fLlXg8npcjkAAAQUDT6MxeZHABAYLW0tBgf4AQ5k0OQAwAI\nrNzcXN10001asGCBkQFO0HG5CggwCo8RdKFQSFdffbXRAU6QJ+ik4zEQYLwGAfe53fG4pKTElWNJ\n0s6dO331nkImBwAAGIkgBwAQCJZladGiRTpw4IDXQ3FVkC9XEeQAAIxnWZbmzp2rnJwcjRw50uvh\nwCXcXQUAMFoswBkzZoyqq6tTNpt4pvBjhsUtwfqXBhCHuatguqAHOEFHJgcIMG4hh+meffbZwAc4\nQc7kcAs5AMBo0WjUVx/0bt9CPmHCBFeOJUnbt2/31ec6mRwAgNH8FOB4Ici/fzBzdwAAwHhkcgAA\nRrAsS4cOHVI4HPZ6KL5CJgdAIFF4DFPE7qJ6+OGHvR4KfITCYyDAeA3CBJl2m7jbhcdlZWWuHEuS\ntm3b5qv3FC5XAQAyVqYFOF7gchUAABmmtbWVAAddIpMDAMhIubm5+s53vqM5c+YQ4HQhyJkcghwA\nQMaaO3eu10OAjxHkAAHG3FWA+YKcyeHuKgAAXOT23VUXXHCBK8eSpDfffNNXn+tcxAQA+J5lWaqs\nrNTevXu9HkrGCYVCrv34DUEOAMDXLMvSvHnz1L9/f5199tleDwcZhJocAIBvxQKccDis6upqZWdn\nez2kjOPHDItbyOQAAHyJAAd9RZADBBhzV8HPnn/+eQKcFAhyTQ53VwEBxmsQfheNRn354dkXbt9d\nNWnSJFeOJUmvvfaar95TyOQAAHzLtAAH7qLwGAAAgwU5UCSTAwDwnGVZ2rVrl9fDgGEIcgAAnord\nRfXII494PRQjBbnwmCAHCDDmroLXOt8mfv/993s9HBiGu6sAAJ4Iah8ct++uuuiii1w5liRt3rzZ\nV5/rZHIAAK6LRCK6/PLLAxfgwF3cXQUAcF12dra+973vafr06QQ4aebHWhm3EOQAADzxxS9+0esh\nwHAEOQAAGCzImRxqcoAAY+4qACYjyAECrKqqyushIAAsy9LVV1+t9957z+uhIGC4XAUASBvLsjR3\n7lyNGTNG4XDY6+EEEperAABIsc4BTnV1tbKy+MiBu8jkAABSjgDHP8jkAACQQi+99JLGjRtHgANP\nkckBAoy5q5Aul19+uS6//HKvhwEFO5PD3FUAALjI7bmrLrnkEleOJUl/+tOffPW5TiYHAACDBTmT\nw4VSAECfWJalbdu2eT0MIAFBDgCg12J3UVVXV3s9FJxCKBRy7cdvCHIAAL3S+TbxH/3oR14PB0hA\nkAMEGHNXobfog5M5gpzJ4e4qIMB4DaI32tvbNXPmTJ1zzjkEOL3g9t1Vn//85105liS9/PLLvnpP\n4e4qAECPZGVlafny5br44osJcOBrBDkAgB5zs/cK+saPl5HcQggOAACMRCYHAACDkckBEEjMXYXu\nWJalK664Qm+//bbXQwF6jEwOEGDcQo6uWJalefPmKRwOa9y4cV4PB71EJucUli5dqry8PJWVlTk+\nvn79ek2cOFHl5eWaNGmSXnrppY7HampqVFJSouLiYt17772pHTUAIK06BzjV1dXKzs72ekhAj3XZ\nJ+ePf/yjBg4cqEWLFjnOS2JZls444wxJ0rZt23TllVfqnXfeUSQS0Xnnnac//OEPKigo0JQpU7R6\n9WqVlpbGH5weHQDgOwQ46eV2n5wvfvGLrhxLkl566SVffa53mcmZNm2ahg4desrHYwGOJH366aca\nMWKEJGnz5s0aN26cwuGwcnNztXDhQq1fvz5FQwYApNOrr76qcePGEeAg4/W5JmfdunX67ne/q337\n9unFF1+UJDU2NqqoqKhjm8LCQm3atMlx/841ARUVFaqoqOjrkAAAfXDppZfq0ksv9XoYxqitrVVt\nba1nxw9yTU6fg5wFCxZowYIF+uMf/6jrr79eO3fu7NH+FD4C3lm2bBmvQSDN7F/gq6qqvBtMwKTs\nFvJp06apra1NR44cUWFhoRoaGjoea2hoUGFhYaoOBSBFeLMFYLI+BTnvvvtuR4HR66+/LkkaPny4\nJk+erPr6eu3Zs0ctLS1au3at5s+f3/fRAgBSyrIsbdmyxethII2CPAt5l5erKisrtXHjRh0+fFhF\nRUWqqqpSa2urJOnGG2/Ur3/9az355JPKzc3VwIEDtWbNmpNPmpOjFStWaPbs2YpEIrrhhhsS7qwC\nAHgrdhdVaWmpJk+e7PVwgJTr8hbytB+cW8gBT/EaDC5uE/eO27eQz5w505VjSdLvf/97X72nMK0D\nAAQMAQ6CgmkdgABj7qrgiUajuuqqqwhwAsSPtTJu4XIVAATMli1bVF5eToDjEbcvV82aNcuVY0nS\niy++6KvPdTI5ABAwFBkHS5AzOdTkAAAAI5HJAQCDRaPRQH+TB5kcAICBLMvS3LlztXXrVq+HAniC\nIAcIMOatMlfsNvG8vDxNmDDB6+HAQ0HueEyQAwQYc1eZiT44wEkEOQBgEAIc4B8oPAYAg7z55ps6\n77zz9POf/5wAB5KCXXhMkAMABrnkkkt0ySWXeD0MwBcIcgAAMFiQMznU5AABxtxVAExGkAMEGLeQ\nZzbLsvTKK694PQz4HLeQAwAySqzR3y9/+UuvhwL4FjU5AJBhYgHOmDFj9NBDD3k9HPicHzMsbiGT\nAwAZpHOAU11draws3saBUyGTAwAZIhqN6stf/jIBDnqETA6AQKLwOLOEQiHde++9BDhAkkLRaDTq\n2cFDIXl4eCDweA0C7nPzdRcKhXTFFVe4cixJWr9+va/eU/gqAAAAjESQAwA+5advxEAmIsgBAB+y\nLEuXXXaZtmzZ4vVQkOFoBggA8A3LsjRv3jzl5+ervLzc6+EAGYtbyIEAY+4q/4kFOOFwWNXV1crO\nzvZ6SMhwfsywuIVMDhBg3ELuLwQ4QGqRyQEAn9ixY4dKS0v1s5/9jAAHKRPkTA5BDgD4xOTJkzV5\n8mSvhwEYgyAHAACDBTmTQ00OAABIm6VLlyovL09lZWUd644cOaKZM2fq3HPP1axZs3T06NGOx+6+\n+24VFxerpKREL774Yp+OTZADBBiFx96xLEt/+MMfvB4GAsDrPjlLlixRTU1N3Lp77rlHM2fO1K5d\nuzRjxgzdc889kqS6ujqtXbtWdXV1qqmp0be+9S21t7f3+ncnyAECrKqqyushBFLsLqq1a9d6PRQg\n7aZNm6ahQ4fGrXvuuee0ePFiSdLixYu1bt06SSfnvqqsrFRubq7C4bDGjRunzZs39/rY1OQAgIs6\n3yb+yCOPeD0coE8OHTqkQ4cO9Xi/AwcOKC8vT5KUl5enAwcOSJI+/PBDfe5zn+vYrrCwUI2Njb0e\nH0EOALiEPjjwQjoLj8866yydddZZHct1dXU9fo7upoToy/i5XAUALvnKV75CgAPoZPZm//79kqR9\n+/Z1BEoFBQVqaGjo2G7v3r0qKCjo9XEIcgDAJffddx8BDlzndeGxk/nz52vVqlWSpFWrVmnBggUd\n69esWaOWlhbt3r1b9fX1uuiii3r9u3O5Cggw5q5y13nnnef1EADXVVZWauPGjTp8+LCKioq0fPly\n3Xnnnbr22mu1cuVKhcNhPf3005Kk8ePH69prr9X48eOVk5Ojhx56qE+Xq0LRaDSaql+kxwcPheTh\n4QEAcJ2bn32hUEjXXXedK8eSpLVr1/rqc53LVQCQBn3p7QEgNQhyACDFLMvSrFmz9Kc//cnroQC+\nrMlxC0EOAKRQ7DbxoqKiuH4fANxH4TEApAh9cOBHfsywuIVMDhBgzF2VOgQ4gP8Q5AABxtxVqbN7\n926VlZUR4MB3glyTw+UqAEiB888/Xz/96U+9HgaATsjkAAAAI5HJAQDAYH68jOQWMjkA0EOWZWnD\nhg1eDwNANwhygABj7qqesyxLc+fO1fr1670eCpCUIBceE+QAAcYt5D0TC3DGjBmjxx57zOvhAOgG\nNTkAkITOAU51dbWysviOiMzgxwyLW3iVAkASFi1aRIADZBgyOQCQhAceeEDnnHMOAQ4yTpAzOQQ5\nAJCE0aNHez0EAD3EVxIgwCg8BszH3VUAAom5q5xFIhGvhwAgBQhyAKATy7J06aWX6qWXXvJ6KAD6\niJocAPg/sdvER48erenTp3s9HCAl/HgZyS1kcgBA8QFOdXW1srOzvR4SgD4ikwMg8AhwYDIyOQAC\nibmrTmpsbNSFF15IgAMYJhSNRqOeHTwUkoeHBwDAdW5+9oVCIS1dutSVY0nSE0884avPdTI5AADA\nSNTkAABgMGpyACAgLMvSM8884/UwALiAIAdAYFiWpXnz5un555/3Vd0AkE5M6wAgkII0d1UswAmH\nw6qurvblGzKA1CLIAQIsKHNX2QMcbhNHkJDJAQCDfeMb3yDAAQKIu6sAGO/+++/X2WefTYADBAxB\nDgDjFRYWej0EwDN+vIzkFi5XAQAAI5HJAQLMxLmr2tralJPDWxsQQyYHQCCZdgu5ZVmaOXOm/ud/\n/sfroQDwAb7uADBC59vEL7vsMq+HA/gGmRwAyGD0wQHghEwOgIxGgAN0jUwOAGSoQ4cOacqUKQQ4\nABJ0GeQsXbpUeXl5Kisrc3x8586duvjii9W/f389+OCDcY/V1NSopKRExcXFuvfee1M3YgApY0Lh\ncTgc1n333UeAA5wC0zqcwpIlS1RTU3PKx4cPH66f/exnuv322+PWRyIR3XTTTaqpqVFdXZ1Wr16t\nHTt2pGbEAFImKHNXAQimLoOcadOmaejQoad8fOTIkZo8ebJyc3Pj1m/evFnjxo1TOBxWbm6uFi5c\nqPXr16dmxAAAAElIS+FxY2OjioqKOpYLCwu1adMmx207p8srKipUUVGRjiEBMIBlWXr22We1ePFi\nr4cCJK22tla1tbWeHd+Pl5HckpYgpycn1ISaAADpZ1mW5s6dqzFjxmjRokWBfuNGZrF/gecysXvS\nEuQUFBSooaGhY7mhoYEJ8gD0WucAp7q6mgAH6IEgv15Scgt5NBqNW548ebLq6+u1Z88etbS0aO3a\ntZo/f34qDgUghTJh7ip7gJOVRecLAMnpMpNTWVmpjRs36vDhwyoqKlJVVZVaW1slSTfeeKP279+v\nKVOm6O9//7uysrL0k5/8RHV1dRo4cKBWrFih2bNnKxKJ6IYbblBpaakrvxCA5GXC5eJvfvObBDhA\nHwQ5kxOK2tMwbh48FErIAgFAZwcOHNDIkSMJcGAMNz/7QqGQvv3tb7tyLEn62c9+5qvPdaZ1AOBr\neXl5Xg8ByGhBzuTw1QgAABiJTA4A3zhx4oT69esX6G+eQKoF+fVEJgcIMD8VHluWpdmzZ+s3v/mN\n10MBYAiCHCDA/NKUrPNt4ldeeaXXwwGMwgSdAOAR+uAASBfeTQB4hgAHQDpReAzAM8eOHdMll1yi\n//f//h8BDpAmfryM5BaCHACeyc/P1w9+8AOvhwHAUAQ5QIBlwtxVAPomyJkc8sNAgPnpFnIASDWC\nHACusCxLjzzyiK/mtQGCgFvIASCNLMvSvHnztGnTJoIcAK6hJgdAWsUCnHA4zG3igAf8mGFxC+82\nANLGHuBkZ2d7PSQAAUKQAwRYTwuP9+3bl/DTlVtvvZUAB/BYkGtyQlEPL5CHQiGuzwMe6ulr0Cmo\nOfvss0+5/UcffaQhQ4YQ4ACduPnZFwqF9J3vfMeVY0nS/fff76vPdWpyAKTN8OHDvR4CgAAjyAEA\nwGB+vIzkFmpyAKREc3Ozr9LUAEAmB8gw3RX7Sl3XyfTFqZ43dhfVN77xDX31q19Ny7EB9A6ZHACB\nlIq5qzrfJr5w4cIUjAoAUoNMDhBgfZ27ij44gP+RyQGAHiLAAeB3ZHKADJOuepueampq0vTp0/X9\n73+fAAfwsSBncghyAPTKiBEj+ny5CwDSiSAHAACDBTmTQ00OEGBkYgCYjCAHCLCqqqqktrMsS//1\nX/+l9vb2NI8IQKoFeYJOLlcBAdDTiTU763wXFR2Nu2Y/z34pEgeCikwOgFPiNnEAmYxMDgBHBDiA\nGfx4GcktZHIAOLrjjjsIcABkNDI5QAA41Ybs27dPt9122ykn/Lz55ps1YMAAHTx4sMvn8ZO+1B6l\ngt/PD4IpyJkcghwgwG677bZTPjZo0CAXRwIAqUeQAwCAwYKcyaEmB4CampoUiUS8HgYApBSZHMBH\nnGpKjh07FrdcUlLS7X49qQ05fvy4rr/+el1zzTWqrKzs8nm8rnkB0HNkcgAEUizAOeecc3Tdddd5\nPRwASClfBDnLli075Rw6y5Ytc2wdzfZsH4TtH3zwQZWWlsb9dLd9fn5+Us9//PhxXXjhhXr11Ve1\nZs0aFRYWKj8/X/n5+XrwwQe7fP7OP346n92N32//vmwfzO3d5jSGdP34TSjqYZ/2UChEm3igE7cu\nV1mWpZkzZ+rVV1/V3r17lZXl/H0n0y5X+X18gOTuZ19XQVg6LFu2zFef677I5ABwV1tbmyoqKiTp\nlAEOADOQyfHq4GRyYIhkMginarrXmT1rIyVmbnqb7XHS09dgMr+Dm5kTJsREJnI7k1NVVeXKsSTp\nrrvu8tXnOl/hAACAkbiFHAAAg/nxMpJbyOQAhrMsSz/84Q/V1taW8Nhdd93lwYgAwB1kcoBuJFNv\nk0wtSDL1NqeaSLM7hw8fjlt+5ZVXJJ3sZHzHHXdoxIgRmjhxYtxs4uFwWAsXLtTOnTtP+bxnnnlm\nt8f+8MMP45aTORepugvKb/U/9nOZbG1UEFA/5Z0gZ3IIcgBDxQKcUaNG6d/+7d/iAhwACAKCHMBA\nnQOcO+64Q01NTV4PCYBHgpzJoSYHMNDKlSs7AhwyOACCikwOYKAbbrhB/fr1I8ABEOhMDkEOkAKv\nvfZawrozzjgjbjmZIt5kiozffffdhHUjRozo8bHOPvtsLes0b5zT72CXqkLaTCw6TWbMFBqfWib+\nmyPzcbkga+NiAAAgAElEQVQKCDA3O6EC8EaQp3UgyAEyXHNzs1paWrweBgD4DkEOkMGampr07//+\n73r22We9HgoA+A41OUA3kmnQl5+fn7CNU+1Md9t8+umnCduMHDkybvn000+XdDLA+e53v6vCwkLN\nmDEjriGgvTmgvWans9jvYq8hkpKrMelNrYVT/c+kSZPilnvTWC9VTQaT4eWxqG9BT/jxMpJbyOQA\nGaipqUm33nqr8vPz9b3vfY+7qADAAZkcIMOkMsC57bbbUjgyAH4U5EwOQQ6QgaZOnaqvfOUrfc7g\nEOQAMBlBDpBhTj/9dF1//fVeDwNAhiCTA+CUkikwdZrJO1W3ddsLgt98882EbYYMGRK3fNlll8Ut\nO41vx44dcctO47U3FUym4NWpqNhemO1UqG3Xm0JjU2dAp9AY6B2CHAAADBbkTA53VwE+1tzcrJ/+\n9Kc0+wOAXiDIAXyqublZd911l/bv35+2W8RXrVqVlucF4B9BntYhFI1Go54dPBSSh4cHkrJhw4Zu\nt7EsK2FdcXFx3PKHH36YsE1bW1vccqzxX1NTk+644w6NGjVKS5cujQtyYs0AOxs4cGCX43PaJz8/\nX/n5+R3jcqr7eOWVV+KWnZoK9qZuJ5m6GL81wPPbeJC53PzsC4VCuv/++105liR95zvf8dXnOjU5\ngM90DnDuuOMOHTlyxOshAchgfsywuIXLVYDPPPXUUx0BDp2MAaD3yOQAPrN48WLl5OQQ4ABAHxHk\nAD5z2mmneT0EAAYJ8uUqghzAxt7Mzqlx3bFjx3q8TTKcZiF3Khq227t3b9yyvYGg03O89dZb+upX\nv6q33npLUmJzQCnx93IqsLYXHiczw3hv9KZY2Wmb3qLQGMg8BDmAh5qbmxUKhTzL3nzta1/z5LgA\n3BPkTA6Fx4BHmpub9YMf/EAvvvii10MBACORyQE80NzcrOXLlysvL09z5871ejgADBbkTA5BDjJW\nMs3Z7Nu8++67cctTp05N2KepqSluOZmamDfeeCNhXf/+/eOWI5GIJOnEiRP60Y9+pJEjR+q6667T\nxx9/3LGNU02OvdFfMnddXXDBBd1us3Xr1rjl3tacODU5tEvm38o+iWiqJui01wiloj7Ia+msPQJM\nQpADuKhzgLNkyZKOwAcA0iXImRxqcgAXZWVladKkSVqyZImysrx/+a1fv97rIQBA2nj/LgsESG5u\nrmbOnOmLAEciyAGCIMgTdHK5Ckaz1ynY60d++9vfJuzz97//PW7Z6fZu+zZOtTT79++PW7b3k5Gk\nwYMHxy13N9GmJLW2tiasswdN9p43TseOHcv+Z1fs9UpSYl2TU58c+3l36iFk78FjnxzUqTbK3g/I\n6XntNTj22h8pufofu1Q9T29QfwMkxx9fJwEAAFKMIAdIk5aWFr3wwguOmRcAcEuQL1d1G+TU1NSo\npKRExcXFuvfeexMe//jjj3XllVdq4sSJ+uxnP6vt27cnvS9gqhMnTuiXv/ylLMtiok0A8EiXQU4k\nEtFNN92kmpoa1dXVafXq1QnX+n/4wx/qwgsv1FtvvaUnn3xSt9xyS9L7AiY6ceKEVq5cqSFDhmjB\nggW+KTJ28uUvf9nrIQBIsyBncrosPN68ebPGjRuncDgsSVq4cKHWr1+v0tLSjm127NihO++8U5J0\n3nnnac+ePTp48KDefffdbvcF+sJefOlUCGoveG1paYlbbmtrS9hnwIABccvDhw9P2MY+IWY0Gu14\n/l/84hcaOnSorrrqqrgAp7GxMeF5jh8/HrdsbyAoJRbXOikvL49bthfgOhUeS9K//uu/dvzdfm4k\n6fDhw3HLTtvYC42dxptMQa79389p0lM7+/8Be7Gy0/OmqjjYrSLjVEqmeSJgki6DnMbGRhUVFXUs\nFxYWatOmTXHbTJw4Ub/5zW80depUbd68We+//7727t2b1L6StGzZso6/V1RUqKKiope/CuC9l19+\nWUOHDtU111zj9VAA+ERtba1qa2s9O74fMyxu6TLISebE3HnnnbrllltUXl6usrIylZeXKzs7O+mT\n2jnIATLdF77wBYVCIWVlZam9vd3r4QDwAfsX+KqqKu8GEzBdBjkFBQVqaGjoWG5oaFBhYWHcNoMG\nDdITTzzRsTx69GiNHTtWTU1N3e4LmIYiYwB+QybnFCZPnqz6+nrt2bNH+fn5Wrt2rVavXh23zbFj\nx3T66aerX79+evzxxzV9+nQNHDgwqX2BVLI3k3Nib3hnrzmRpI8++ihu2T6pp5Pm5uaEdbE6nZh+\n/folbGOvyXG63dxem+LUeLC7STKdmuSdqk6ns7Fjx/b4eZzOl71ux6nexl7jkkwzQDv7eBGPGhwE\nTZdBTk5OjlasWKHZs2crEonohhtuUGlpqR599FFJ0o033qi6ujp9/etfVygU0vnnn6+VK1d2uS9g\nipaWFkWjUceOyJli1apVWrx4sdfDAJBGZHK6MGfOHM2ZMydu3Y033tjx94svvlhvv/120vsCJmhp\nadGqVat07rnnavr06V4Pp9eefPJJghwAxmLuKqCHYgHO0KFDNW3aNK+HAwBdCnImx79dygAf6hzg\n2PvgAAD8hUwOfCmZpmX2Jm+HDh1K2MZeyNv5jr9T7WOfYTxWQNza2qo1a9bozDPP1Oc+97m4Wcad\nbhdPpolfXl5e3PInn3ySsE0kEun2ed5///245YKCgrjlnJzEl3qssV/sz5EjRyZs87//+79xy06/\nk70g2Ol57DOBO81U7lTU3B17kXMyhcf2/1vJSmZ89iJsCn0BbxHkAEnKzs7WhAkTdMEFFzgGIwDg\nR0G+XEWQAyQpKytLF154odfDSKmvfvWrXg8BANKGIAcIsK997WteDwFAmpHJAXzGXsvgVEdhb/7n\n1Gyvvr4+bvnIkSNxyydOnEjYx96gz+nSlL02xakmxz75p1NdjL2xn1PNS1NTU8I6u3POOSdu2X5u\nnBr/2ZseOkmmYaC9weLf/va3hG3szQqdfs/u6nacGggmU/Nifx77cZw4TfZqb1botA0Af+HWEMBB\na2urNm7c6DjrNgBkklAo5NqP3xDkADatra167rnn1NLS4ph9AQBkBt7BgU5aWlr03HPPafDgwZox\nYwZ9cABkPD9mWNzCOzjwf1paWvTMM88EKsB55plnvB4CAKQNmRxkBKdZtt977724ZafZue3NAO2z\njndu8LZ9+3bl5uZq/PjxOnDgQMd6pwk47Q0DnWRnZ3e5LCUWFTsV+tqLo50uodkLqu1jdqotGjRo\nkJ555hktXbpUknNxt73Rn1MRtL3AevDgwQnb5Obmxi07FYnbZx1PZgb0ZAqa7Q0f7ceREn9Pp+Jk\ne6GxvRDZSTJNLbs7TrLHSkZvxoPMF+RMDkEO8H9KS0slJQZGAIDMRJAD/J8gXJ4CgCAhyAEAwGBc\nrgI85DRZo53TRJqNjY1xy/Z6G6d1H3/8saSTdSSRSMSxVsVeLxKboLMzew2HU72N/Y3Fvo8kDRgw\nIG65ubm522PZG/1JiY3y+vfvH7e8d+/ehH1iTfyOHj0at9yV4uLihHX2c3jw4MGEbex1RU41Oano\nSeT0fyAcDsctO9W3JNPYrzd1Mb2peUlV/Y0TanAQNAQ5CJy2tjb98Y9/VF5enoqKirwejqeuvvpq\nr4cAIM2CnMmhCAGBEgtwzjjjjI5C4yC75pprvB4CAKQNmRwERmtra0eAM2XKlEB/uwEQHEF+ryOT\ng0Boa2vT7373OwIcAAgQMjlIq2Qam9lnCpcSm8459a6JFRHHOBXXxopgo9GoRo4cqWHDhsUVp0Yi\nkYR97M3tnGYqtxf2OgVN9oJhp8Jae5O8QYMGJWxjb/TnVORsPxf2xohDhw5N2Mf+PE7N9uzNCZ3O\ncXfPKyU2MHQqPLZfPtyyZUvcslOjv+6eI1nJzLbeGzTfgx8E+UsdmRwEQigUUkFBQaBf7AAQNAQ5\nQIAxdxVgvlAo5NqP3xDkAAFGkAPAZNTkwHX2Oh17DYyUWB+yZ8+ehG3s62KTZkYiEb3zzjsaM2ZM\nwnPbl50a9NnrR5y2sT+P0ySe9m81TvU/yUyk2d7eHrfs1LTPPmGovWbI6djHjx+P+zOZZoBOnP79\nuuN0Tnfs2BG3bG9waF924jSRq32/3k6A2Zv6Gmpw4Ad+zLC4hUwOjBKJRLR161a1t7c7ztYNAAgO\nghwYIxbg9O/fXyUlJYH+9gIA4HIVDNHW1kaAAwAOgvx+SJCDtHKaTPK9996LW3aalPKdd96JW7ZP\nxinF10gcOHBA7e3tGjx4cNx6++Sa9rqTrKzEZGYyNTmDBw/uch8psb7Gqc+LvYbE6VjDhg2LW7b3\nwJESz6G9/47TPtnZ2ZozZ05HPU8yNTkjRoxIWGevyXHqOZPM5Jv2/yv2uiynmhx7zYtTTY593aRJ\nk7odi1PdTrp66QBIH4IcGOGss87SaaedFuhvLL0xZ84cr4cAIM2C/L5ITQ6M4NceDQAA75DJAQDA\nYEH+AkiQg4zT1tbWq94sAIBgIchBSr322mtxy4cOHUrYZv/+/XHLThN0vv3223HLscLj9vZ2NTQ0\naMCAAQkN75wmfbRPFmlfthcmS4mFtE7fgmIN9GKcJsC0F6rai4GdtrE/r5RYtGtvDiglFizbf0+n\nZoD2YmWnwmN7MbDTNvZ/B6ci42SKmu2F2U5F6z19jmTZG/0l0xwQyBRBzuRQk4OMEQtwcnNzHe/y\nQc/9/ve/93oIAJA2BDnICJ0DnLPPPjvQ30xSiSAHMB8TdAI+RoADAOgNanKQUvb6msOHDydss3Xr\n1rhlp8ZrnRu4RaNRhUIhZWVl6eDBg3HrO3MKfuw1GvYJMZ2a79nXOT2vvb6mtbU1YRt7bY99Ek2n\n8ThN9Gl/ngEDBnT7PPbLeU7HHj58uKR/NNRzugTY1NQUtzxy5MiEbZK5dGhvyJfMPslMyGmvpXGS\nTH0NE2nCZEH+YkiQA98LhUKOH/4AAHSFy1UAACBtli5dqry8PJWVlXWsO3LkiGbOnKlzzz1Xs2bN\n0tGjR9NybIIcIMCuvPJKr4cAIM28LjxesmSJampq4tbdc889mjlzpnbt2qUZM2bonnvuScvvTpAD\nX4lEItq1a5djXxek3lVXXeX1EAAYbtq0aQm9xJ577jktXrxYkrR48WKtW7cuLcemJgdJsxd5btmy\nJWEb++zh9mVJ2rZtW9xyrCi1vb1dR48eVXZ2tj799NNui+UGDRoUt+zUCM7ebM/euM6piNf+PE4z\nlTsVLNvZm+I51RXZn9vexE86OfloZ07FvydOnIhbtjdhdGrGZ5853al42s5pNnP7Oqfx2Rs1OhUD\n2/9/JVMMnEzhcTJ6c2wgU6Sz8Hj79u2qq6vr8X4HDhxQXl6eJCkvL08HDhxI9dAkEeTAJzoHOIMH\nD3b8MAUA+MuECRM0YcKEjuVnnnmmx8+Rzh47BDnwnD3ACfLtjgCQan58T83Ly9P+/fs1atQo7du3\nLyFjnSrU5MBzzc3NBDgAECDz58/XqlWrJEmrVq3SggUL0nIcMjlI2o4dO+KW9+zZk7DN9u3b45ad\nGv3t3r07bjnWqK5z4zl7vUgyE0M6bWOvr7HXhtifQ0qsi3GaWNPeiNBpG/uxnGZOt9cVOU30aZ+Q\n0+l57E0X7eNxqjH55JNP9Mwzz+jLX/6ypORqmpy2sU+k6bSNvbGffSJXp22Sqbex/169raWhBgcm\n8/rLY2VlpTZu3KjDhw+rqKhIy5cv15133qlrr71WK1euVDgc1tNPP52WYxPkAAHWOcgBgHRYvXq1\n4/o//OEPaT82QQ4AAAbzOpPjJWpy4KpIJJJw+QUAgHQgyIFrIpGIdu/endA/BgCAdOByFSQ5Fwjb\nC43ffvvtuOXNmzcn7GNv9NfQ0CDpZKFua2urQqGQ2tvbE5rX2dmLV50KcmOzaMc4NaHr7nmdGvTZ\nt3EKyuzFyU4FzMkUT9sLmD/++OOEbey/u9O5izXV6smxY+tiRbdOvYns5yKZGb2TQaEv4B4uVwFp\n1DnAsX/4wltLlizxeggAkDZ84iCt7AFOkL9R+NHSpUu9HgKANAvy+y5BDtIqFAopOztbWVlZgX6h\nAQDcR5ATUPYaHHtDNymxsd/rr78et2yvv5GcJ+TsjlO9yLBhw+KWhwwZkrCNvZGevXGd1P1Emva6\nHinxXNjHIiU223NqBlhUVBS3HGt62Jm9JsjpXNi3SWbM9lqaESNGJOyTTH2NU61WKjg1+qNOB0iP\nIH/BpCYHAAAYiUwOAAAGI5MDpABN/jLPsmXLvB4CAKQNQQ5Sor29PakJFeEvVVVVXg8BQJqFQiHX\nfvyGy1UBZZ+1+qWXXkrY5o033ohb/tvf/ha33JsiYymxcd5ZZ52VsI29uNapGaB9nb0QWUos2rUX\n5DoVJtuLdJ0aBtqbATptYzdp0qSEdfaZ3JN5ntLS0oR1x44di1u2//s6FWXbORUZ24uTnWYPd/q9\nukORMQA3EOQAAGAwP2ZY3MLlKgAAYCSCHAAAYCQuVwWAUx1FbW1t3LJTY7+33norbnn37t09Pvbg\nwYMT1uXn58ctjxo1KmEbe12MU7M9e4NAp23sjfzsNS9tbW0J+2Rlxcf+TpNv2p/Hvo+U+Ls71f9M\nnjw5YZ2dvdGf/d9FSpycdOrUqd0+ryTdddddHX9Ppjlgb+pvAHiLy1UAAolbyAGYjEwOAAAGI5MD\ndINGfwCATEMmJwDs/W0kqa6uLm55y5YtCdvYe7gkw97zZvTo0Qnb2GtVnCaPtG/jNEGn/VhOgZhT\nnU5nY8aMSVjX3Nzc5T5O+9nrjHqL/jEAUo1MDgAAgGEIcoAAo/AYMF+Qp3UgyAECjLmrAJiMmhwA\nAAzmxwyLWwhyDFRdXR23/Oc//zlhG3sx8gcffNDj4xQVFSWssxfkOhXS2ouInSbfzMmJ/6/pNImn\nvUmfvWmeJOXl5cUt2yfodCoYtk926bRNMo3zAADe4nIVAAAwEpkcSJKi0aii0ahvi8cAAL0T5Pd0\nMjlQe3u7otGo18OABzrPXQUApiGTk+E2bNiQsO7111+PW37jjTcStnnvvfcknczgdG5+F8voSIn1\nK/YanNLS0oTntdfbhMPhhG3sNS7Z2dkJ29gns3RqBmjX2tqasK64uDhu2d7g0Knexv57m1x/wy3k\ngPnI5CCQYgGO0wzaAABkOjI5ARaJRJSVlaV+/fqpra3N6+EAANKATE4XampqVFJSouLiYt17770J\njz/wwAMqLy9XeXm5ysrKlJOTo6NHjya1L7yVk5Ojfv36BfoFAAAwV5dBTiQS0U033aSamhrV1dVp\n9erV2rFjR9w2t99+u9544w298cYbuvvuu1VRUaEhQ4YktS+8R4ADAGYL8rQOXV6u2rx5s8aNG9dR\nPLpw4UKtX7/eseBUkn71q1+psrKyV/siOa+88krc8rZt2xK2sRfX7t+/P2Eb+4zdo0aNSthm3Lhx\nccv2Il57oz1JOvfccxPW2dmfx4nTzOR29oLgnTt3drtNMkXEQZoJfNmyZRQfAzBWl0FOY2Nj3B01\nhYWF2rRpk+O2x48f1wsvvKCHHnqoR/t2foOtqKhQRUVFT8aPJLW3t6u9vT2hkzCCraqqiiAHSLPa\n2lrV1tZ6dnw/Zljc0uUnXk9OzIYNGzR16tSOW32T3Zc32PRrb2/X0aNHlZOTo8GDB3s9HAAIFPsX\neCbGdU+XQU5BQYEaGho6lhsaGlRYWOi47Zo1azouVfV0X6RPW1ubjh49quzsbA0aNMjr4QAAXEYm\n5xQmT56s+vp67dmzR/n5+Vq7dq1Wr16dsN2xY8f08ssv61e/+lWP90XP2OtOnIq5Y5Nttre36/33\n31d2drYGDx4c9x/dPpHmeeedl/A8559/ftxyQUFB3PJnPvOZhH0GDhwYtzx16tSEbex1RU7b9IbJ\nTfsAAD3XZZCTk5OjFStWaPbs2YpEIrrhhhtUWlqqRx99VJJ04403SpLWrVun2bNnx3WpPdW+cEc0\nGtUHH3ygfv366bTTTgt0JA8ACKZuq1DnzJmjOXPmxK2LBTcxixcv1uLFi5PaF+4IhUIaOXKkBgwY\noCNHjng9HPgUc1cB5gvyl1xutTGYfQ4mwI7CfwAmI8gBAMBgZHLgSzU1NQnr6uvr45YbGxsTtjnt\ntNPilv/pn/4pYRt74fGkSZO63aalpSVu+bLLLkvYJxnJFBrv27cvbjlIDfoAAKnB9NMGiEQi2rZt\nm44fP+71UAAAPhPkaR0IcjJcJBLR1q1blZOTE3d3GwAAQUeQk8FaWlq0detW9e/fXyUlJb6MouFv\nFB4D5gtyJicUjUajnh08FJKHh/e9NWvWJKyLzX/S2tqq3/3ud8rKylJZWVncf66xY8fG7ePUado+\nkaZTD6Mzzzwzbpm6GPPwGgTc5+brLhQK6eWXX3blWJL0+c9/3lfvKRQeZ6j9+/frzDPPVDgc9mX0\nDADwhyB/RhDkZKiioiIVFRV1TOEAAADiUZMDAACMRCYHAACDcbkKrrM3u5OkDz/8MG750KFDkqQT\nJ07o0KFDKiws1LBhw+K2cSoYLisri1v+4he/2NfhwlDMXQXAZFyu8rkTJ07oscce06uvvur1UGAg\nbiEHzBfkW8gJcnwsFuAMGzZMV199tdfDAQAgo3C5yqeampo6ApzKykplZRGPAgB6zo8ZFrcQ5Hjk\n2LFjCevOOOMMSVI0GtW3v/1tjRgxQosWLYoLcBYtWtTtc5eUlKRuoAAAZCiCHB8KhUK6+eabdfTo\nUTI4AIA+CXImh09Qn7rgggsIcJB2FB4DMBmfokCAVVVVeT0EAGnG3VUAAACGoSbHI5ZlSTp5F9V3\nv/td3XzzzZo6dWrcNpdffnnCfswEDgDoCT9mWNxCJsdDTU1NuvXWWzV06FCdc845Xg8HAACjEOR4\nJBbg5Ofn63vf+56ys7O9HhIAAEbhcpUHLMsiwIEvMHcVYL4gX64iyPHAli1bdOGFF+rhhx8mwIGn\nuIUcgMlC0Wg06tnBQyF5eHgAAFzn5mdfKBTSX//6V1eOJUlTpkzx1ec6NTkAAMBIXK4CAMBgQa7J\nIZOTZpZl6S9/+YvXwwAAIHAIctLIsizNmzdPq1at8noogCMKjwHzBXlaBwqP0yQW4ITDYVVXV3MX\nFXzJ5Ncg4FduFx6/9tprrhxLkiZNmuSr9xRqctKAAAcA4Bd+zLC4hctVKRaNRnXNNdcQ4AAA4DEy\nOSkWCoV09913q6ysjAAHAOC5IGdyCHLS4IILLvB6CAAABB6Xq4AAY+4qACbj7qo+ikajgU4FAgB6\nxu27q958801XjiWdvJLhp891Mjl9YFmWLrvsMr3xxhteDwUAANhQk9NLlmVp7ty5GjNmjCZOnOj1\ncAAAcBTkqw1kcnqhc4BTXV2trCxOIwAAfkMmp4cIcAAAmYRMDpK2bds2jR8/ngAHRmDuKgAm4+4q\nIMB4DQLuc/vuqm3btrlyLEkqKyvz1XsKqQgAAGAkanIAADAYNTlwZFmWamtrvR4GAADoBYKcU4jd\nRbV69WqvhwIAAHqBy1UOOt8m/vDDD3s9HCBtmLsKMF+QL1dxd5UNfXAAAOnk9t1V27dvd+VYkjRh\nwgRffa6TybG57rrrCHAAAMYIciaHIMfmvvvuU0lJCQEOAAAZjiDHZvz48V4PAQCAlAlyJod0BQAA\nMFKggxw/FUcBXmDuKsB8oVDItR+/CWyQY1mWZs2apb/85S9eDwXwTFVVlddDAIC0CWRNjmVZmjdv\nnsLhsKZMmeL1cAAASBs/ZljcErhMTucAp7q6WtnZ2V4PCQAApEGgMjkEOACAoCGTExD19fWaMGEC\nAQ4AAAEQqEzOBRdcoBUrVng9DMA3mLsKgMmYuwoAABe5PXdVfX29K8eSpOLiYl99rgfqchUAAAgO\nY4Mcy7JUU1Pj9TAAAPAUzQANE7uL6tlnn/V6KAAAwCPGFR53vk380Ucf9Xo4AAB4yo8ZFrcYlcmh\nDw7QM8xdBcBkRt1ddeWVV2rIkCEEOECSuMMRcJ/bd1e99957rhxLksaMGeOr9xSjgpx33nlHo0eP\nJsABkkSQA7iPIMc9RtXkjBs3zushAADgK9TkAAAAGCZjg5z29navhwAAAHwsI4Mcy7J06aWXauPG\njV4PBchozF0FmC/IzQAzrvDYsizNnTtXY8aMUXV1tbKyMjJOAwAElNuFx3v27HHlWJIUDocpPO4t\nAhwAAHrGjxkWt2RMlECAAwAAeiJjMjkffPCBysvL9eCDDxLgAACQpCBncjKuJgcAgEzmdk3OBx98\n4MqxJOkzn/mMrz7XSYkAAcbcVYD5uLvKq4OTyQE8xWsQcJ/bmZyGhgZXjiVJRUVFvnpP8UUm58c/\n/nHcsmVZ+u1vf+vRaMxUW1vr9RCMxzlOP85x+nGOzRPkTI4vgpx169Z1/D12F9V///d/+yoazHS8\ncaUf5zj9OMfpxzmGSXx1d1Xn28Qff/xxX0aFAABkkiB/lvoikyP9I8AZPXo0fXAAAECfeV54DABA\n0LhZeNzY2OjKsSSpoKDAV6Umnl6u8tOJAADAREFOKHBNCAAAGMlXhccAACC1yOQAAAAYhkwOAAAG\nI5MDAABgGDI5AAAYjEwOAACAYcjkAABgMDI5AAAAhiHIAQAARuJyFQAABuNyFQAAgGHI5AAAYDAy\nOQAAAIYhkwMAgMHI5AAAABiGTA4AAAYjkwMAAGAYMjkAABiMTA4AAIBhyOQAAGAwMjkAAACGIcgB\nAABG4nIVAAAG43IVAACAYcjkAABgMDI5AAAAhiGTAwCAwcjkAAAAGIZMDgAABiOTAwAAYBgyOQAA\nGIxMDgAAgGEIcgAAgJG4XAUAgMG4XAUAAGAYMjkAABiMTA4AAIBhyOQAAGAwMjkAAACGIZMDAIDB\nyOQAAAAYhkwOAAAGI5MDAABgGDI5AAAYjEwOAACAYQhyAACAkbhcBQCAwbhcBQAAYBgyOQAAGIxM\nDmAM9zIAAAC3SURBVAAAgGEIcgAAMFgoFHLtx8nSpUuVl5ensrIyl39zghwAAJBGS5YsUU1NjSfH\npiYHAACDeV2TM23aNO3Zs8eTYxPkAABgMDeDnIEDB7p2rGQQ5AAAYKhoNOr1EDxFTQ4AADASQQ4A\nADASQQ4AAEibyspK/fM//7N27dqloqIi/eIXv3Dt2KFo0C/YAQAAI5HJAQAARiLIAQAARiLIAQAA\nRiLIAQAARiLIAQAARiLIAQAARvr/ikO3R4fyGtkAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1ecc4310>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We define a function used for permutation testing of the statistical significance of the differences between the two models. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def permutation_test_rel(arr1, arr2, n_perms=1000):\n",
      "    \"\"\"\n",
      "    Test the difference between two nd arrays\n",
      "\n",
      "    \"\"\"\n",
      "    # Linearize and get the finite elements: \n",
      "    arr1 = arr1[np.isfinite(arr1)]\n",
      "    arr2 = arr2[np.isfinite(arr2)]  # We assume that the indices are good for both\n",
      "    orig_arr = np.vstack([arr1, arr2])\n",
      "    orig_diff = np.median(arr1 - arr2)\n",
      "    perm_diff = np.empty(1000)\n",
      "    for ii in xrange(n_perms):\n",
      "        idx1 = np.random.randint(0, 2, orig_arr.shape[-1])\n",
      "        idx2 = np.mod(idx1+1, 2)\n",
      "        perm_arr1 = orig_arr[idx1, np.arange(idx1.shape[0])] \n",
      "        perm_arr2 = orig_arr[idx2, np.arange(idx2.shape[0])]\n",
      "        perm_diff[ii] = np.median(perm_arr1 - perm_arr2)\n",
      "        \n",
      "    return orig_diff, perm_diff"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The test is applied to each b value between DTM and SFM:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "orig_diff1k, medians1k = permutation_test_rel(dtm_rrmse_1k, sfm_rrmse_1k)\n",
      "orig_diff2k, medians2k = permutation_test_rel(dtm_rrmse_2k, sfm_rrmse_2k)\n",
      "orig_diff4k, medians4k = permutation_test_rel(dtm_rrmse_4k, sfm_rrmse_4k)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We also derive from the permutations a boot-strap confidence interval on the difference of the medians:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ci_1k = stats.scoreatpercentile(medians1k, 97.5) - stats.scoreatpercentile(medians1k, 2.5) \n",
      "ci_2k = stats.scoreatpercentile(medians2k, 97.5) - stats.scoreatpercentile(medians2k, 2.5) \n",
      "ci_4k = stats.scoreatpercentile(medians4k, 97.5) - stats.scoreatpercentile(medians4k, 2.5) "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following function finds the boot-strap distribution of the individual medians"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def boots_strap_medians(arr, n_iters=1000): \n",
      "    \"\"\" \n",
      "    Produce a boot-strap distribution of the median of an array\n",
      "    \"\"\" \n",
      "    # Linearize and exclude nan/infs:\n",
      "    arr = arr[np.isfinite(arr)]\n",
      "    medians = np.empty(n_iters)\n",
      "    for ii in xrange(n_iters): \n",
      "        this_arr = arr[np.random.random_integers(0, arr.shape[0]-1, arr.shape[0])]\n",
      "        medians[ii] = np.median(this_arr)\n",
      "    \n",
      "    return medians"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compute the boot-strap confidence intervals of the median rRMSE in each model/b-value"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "medians_dtm_1k = boots_strap_medians(dtm_rrmse_1k)\n",
      "medians_dtm_2k = boots_strap_medians(dtm_rrmse_2k)\n",
      "medians_dtm_4k = boots_strap_medians(dtm_rrmse_4k)\n",
      "\n",
      "medians_sfm_1k = boots_strap_medians(sfm_rrmse_1k)\n",
      "medians_sfm_2k = boots_strap_medians(sfm_rrmse_2k)\n",
      "medians_sfm_4k = boots_strap_medians(sfm_rrmse_4k)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dtm_ci_1k = (stats.scoreatpercentile(medians_dtm_1k, 97.5) - stats.scoreatpercentile(medians_dtm_1k, 2.5))\n",
      "sfm_ci_1k = (stats.scoreatpercentile(medians_sfm_1k, 97.5) - stats.scoreatpercentile(medians_sfm_1k, 2.5))\n",
      "\n",
      "dtm_ci_2k = (stats.scoreatpercentile(medians_dtm_2k, 97.5) - stats.scoreatpercentile(medians_dtm_2k, 2.5))\n",
      "sfm_ci_2k = (stats.scoreatpercentile(medians_sfm_2k, 97.5) - stats.scoreatpercentile(medians_sfm_2k, 2.5))\n",
      "\n",
      "dtm_ci_4k = (stats.scoreatpercentile(medians_dtm_4k, 97.5) - stats.scoreatpercentile(medians_dtm_4k, 2.5))\n",
      "sfm_ci_4k = (stats.scoreatpercentile(medians_sfm_4k, 97.5) - stats.scoreatpercentile(medians_sfm_4k, 2.5))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We display the median  +/- 95% confidence interval in each condition"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots(1)\n",
      "ax.bar([1,2], \n",
      "       [np.median(dtm_rrmse_1k[np.isfinite(dtm_rrmse_1k)]), \n",
      "        np.median(sfm_rrmse_1k[np.isfinite(sfm_rrmse_1k)])],\n",
      "        yerr=[dtm_ci_1k, sfm_ci_1k], color=[0.8, 0.8, 0.8], ecolor='k')\n",
      "\n",
      "ax.bar([4,5], \n",
      "        [np.median(dtm_rrmse_2k[np.isfinite(dtm_rrmse_2k)]), \n",
      "        np.median(sfm_rrmse_2k[np.isfinite(sfm_rrmse_2k)])], \n",
      "        yerr=[dtm_ci_2k, sfm_ci_2k], color=[0.59, 0.59, 0.59], ecolor='k')\n",
      "\n",
      "ax.bar([7,8], \n",
      "        [np.median(dtm_rrmse_4k[np.isfinite(dtm_rrmse_4k)]), \n",
      "        np.median(sfm_rrmse_4k[np.isfinite(sfm_rrmse_4k)])], \n",
      "        yerr=[dtm_ci_4k, sfm_ci_4k], ecolor='k', color=[0.32, 0.32, 0.32])\n",
      "ax.set_ylim([1/np.sqrt(2), 0.8])\n",
      "ax.set_xlim([0.8, 9])\n",
      "ax.grid()\n",
      "ax.set_xticks([0])\n",
      "fig.set_size_inches([10,6])\n",
      "fig.savefig('figures/rrmse_bars.svg')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAFhCAYAAABZBV+UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFxlJREFUeJzt3W1sV+XBx/HfX0sMm9P60MytZam3gBQ1wlqyzCyxc1s6\nTcSFaYS4CA4nCTGby0wWl70weyFhY5kOZGEucVtYoUvoZPOhJGSrmhlLYMQlgwgzVitLGERaVIy4\n2vvFfbezAcqDVx+sn8+rXv9zTq/rnED49n8O/1YGBgYGAgDAB3LWeC8AAGAyEFUAAAWIKgCAAkQV\nAEABogoAoABRBQBQwEmjqqOjI7NmzcqMGTOycuXKY7YfPHgwX/3qVzNnzpxceeWV+fWvf33KxwIA\nTBaVkT6nqr+/P5dffnm2bt2a2trazJs3Lxs2bEhDQ8PQPvfff3/eeeedrFixIgcPHszll1+e/fv3\np1KpnPRYAIDJYsR3qrZt25bp06envr4+U6ZMycKFC7N58+Zh+3zqU5/K4cOHkySHDx/ORRddlKqq\nqlM6FgBgsqgaaeO+ffsybdq0oXFdXV26urqG7fOtb30r1113XT796U/njTfeyO9///tTPrZSqXzg\nEwAAGCsj/SKaEaPqVKLngQceyJw5c9LZ2ZmXXnopX/nKV/LCCy8UWRwAwERxsi4a8fZfbW1tenp6\nhsY9PT2pq6sbts9zzz2XW265JUly2WWX5dJLL82LL76Yurq6kx4LADBZjBhVTU1N2bt3b7q7u3P0\n6NG0tbVl/vz5w/aZNWtWtm7dmiTZv39/XnzxxfzP//zPKR0LADBZjHj7r6qqKmvWrElLS0v6+/uz\ndOnSNDQ0ZN26dUmSZcuW5Qc/+EHuuOOOXH311Xnvvffy4x//OBdeeGGSHPdYAIDJaMSPVBj1ySsV\nz1QBAB8KJ+sWn6gOAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIA\nKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgq\nAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFBA\n1XgvAAA+Kjo7O9PZ2Tk0HhgYSKVSSXNzc5qbm8dtXZRRGRgYGBi3ySuVjOP0ADBuDh8+nLq6uhw+\nfHi8l8IpOlm3uP0HAFCAqAIAKMAzVQDw/84///wxvx1XqVTGZJ7zzjsvfX19YzLXR5VnqgDg/1Uq\nlSxYsGBM5nr33Xfz1FNPZf78+WMyX3t7u39zPyDPVAEAjAFRBQBQgGeqAGCMHDhwIAcPHkyS9Pf3\n57333svu3btz8cUXp6amZpxXxwclqgBgjNTU1AyLpyuvvHIcV0Npbv8BABQgqgAAChBVAAAFiCoA\ngAJEFQBAAaIKAKAAUQUAUICoAgAoQFQBABQgqgAAChBVAAAFiCoAgAJEFQBAAaIKAKAAUQUAUICo\nAgAoQFQBABQgqgAACqga7wUAH26dnZ3p7Ow85vXm5uY0NzeP+XoAxouoAj6Q98fTa6+9lvb29nz7\n298e30UBjIOTRlVHR0fuueee9Pf3584778z3v//9YdtXrVqV3/3ud0mS//znP9m9e3cOHjyY6urq\nrFixIuvXr89ZZ52Vq666Ko8++mjOOeec0TmTgvzkDWdm//79+c1vfiOqgI+kysDAwMCJNvb39+fy\nyy/P1q1bU1tbm3nz5mXDhg1paGg47v6PP/54HnzwwWzdujXd3d257rrrsnv37pxzzjm59dZbc8MN\nN2Tx4sX/nbxSyQjTTwhPPvlknn766axcuXK8lwIT3o4dO3LXXXdlx44d470UOCOVSiULFiwY72WM\nivb29gn/b+5Ed7JuGfGdqm3btmX69Ompr69PkixcuDCbN28+YVS1trZm0aJFSZLzzjsvU6ZMyZEj\nR3L22WfnyJEjqa2tPcPTGD99fX159dVXx3sZcMYuuOCC9Pb2jumclUplTOaprq7OoUOHxmQugJMZ\nMar27duXadOmDY3r6urS1dV13H2PHDmSLVu2ZO3atUmSCy+8MN/73vfymc98JlOnTk1LS0u+/OUv\nH3PckiVLhqKturo6c+bMGbrFNngLbjzHu3btGlrrRFiPsfHpjnt7e7N27drs2bMnSTJz5swkGZXx\nv//97zz77LO57777xmS+Bx98MIMmyvU2/nCPBx04cCBJUlNTM6nGgybK9Z7o48Gvu7u7cypGvP23\nadOmdHR05JFHHkmSrF+/Pl1dXVm9evUx+7a1taW1tTWbN29Okrz00ku58cYb8+yzz+b888/PLbfc\nkptvvjm33Xbbfyc/w9t/F1544aT86fSCCy7I66+/Pt7LYJKpVCpDP+yMtldffTW/+93vct99943J\nfMuXL3c7g6Lc/mMkH+j2X21tbXp6eobGPT09qaurO+6+GzduHLr1lyTbt2/PNddck4suuihJsmDB\ngjz33HPDoupMHTp0KNu3b//A3+dUbNmyJU8//XQeeOCBUZ+rqalp1OcAAEbHiB/+2dTUlL1796a7\nuztHjx5NW1tb5s+ff8x+fX19eeaZZ3LTTTcNvTZr1qw8//zzefvttzMwMJCtW7dm9uzZ5c8AGFd7\n9uzJE088kSeeeCLPP/98zj777DzxxBNDt+sAPipGfKeqqqoqa9asSUtLS/r7+7N06dI0NDRk3bp1\nSZJly5YlSR577LG0tLRk6tSpQ8deffXVuf3229PU1JSzzjorn/3sZ3PXXXeN4qkA42HmzJlDzzkB\nfJSN+EzVqE9+hs9UVSqVUb39t2PHjqH/Ev7Pf/4zr7zySr70pS+lsbExjY2NozZvU1OT+90UN5bP\nVI01z1RRmmeqGMkHeqbqo+r98fTmm2/mrbfeyic/+clxXhUAMJGJqpM499xzc+655473MgCACW7E\nB9UBADg1ogoAoABRBQBQgKgCAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoAC/poaT6uzs\nTGdn5zGvNzc3p7m5eczXAwATkajipN4fT08++WR+9atfpb29fXwXBQATjNt/nJb33nsvR48eHe9l\nAMCEI6oAAApw+28SuPDCC3Po0KExnbNSqYz6HBdccEFef/31UZ8HAEoQVZPAoUOHsn379jGZ69ln\nn82mTZvy4IMPjvpcTU1Noz4HAJTi9h8AQAHeqeKkduzYkR07diRJXnnllbz66qv55S9/mcbGxjQ2\nNo7z6gBgYhBVnNT74+mdd97J22+/nerq6nFeFQBMLKKK03LOOefknHPOGe9lAMCE45kqAIACRBUA\nQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBU\nAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIAC\nRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIA\nKEBUAQAUIKoAAAoQVQAABYgqAIACThpVHR0dmTVrVmbMmJGVK1ces33VqlWZO3du5s6dm6uuuipV\nVVXp7e1NkvT29ubmm29OQ0NDZs+eneeff778GQAATAAjRlV/f3/uvvvudHR0ZNeuXdmwYUN27949\nbJ977703O3fuzM6dO7NixYo0Nzenuro6SfKd73wnN9xwQ3bv3p2///3vaWhoGL0zAQAYRyNG1bZt\n2zJ9+vTU19dnypQpWbhwYTZv3nzC/VtbW7No0aIkSV9fX5599tl885vfTJJUVVXl/PPPL7h0AICJ\no2qkjfv27cu0adOGxnV1denq6jruvkeOHMmWLVuydu3aJMnLL7+cmpqa3HHHHXnhhRfS2NiYhx56\nKB/72MeGHbdkyZLU19cnSaqrqzNnzpw0NzcnSTo7O5PkmPGgHTt2JEkaGxsnxXjwHE92/q6H8emM\nB+3ZsydJMnPmzEk1HjRRrrfxh3s86MCBA0mSmpqaSTUeNFGu90QfD37d3d2dU1EZGBgYONHGTZs2\npaOjI4888kiSZP369enq6srq1auP2betrS2tra1D72Rt3749n//85/Pcc89l3rx5ueeee3Leeefl\nRz/60X8nr1QywvQnXnSlku3bt5/2cRNdU1OT6/E+Z3o9GK5SqQz9sDPZLF++3J8RiqpUKlmwYMF4\nL2NUtLe3+/vyAZ2sW0a8/VdbW5uenp6hcU9PT+rq6o6778aNG4du/SX/965WXV1d5s2blyS5+eab\n87e//e20Fg8A8GExYlQ1NTVl79696e7uztGjR9PW1pb58+cfs19fX1+eeeaZ3HTTTUOvXXLJJZk2\nbdrQW/Rbt27NFVdcUXj5AAATw4jPVFVVVWXNmjVpaWlJf39/li5dmoaGhqxbty5JsmzZsiTJY489\nlpaWlkydOnXY8atXr85tt92Wo0eP5rLLLsujjz46SqcBAHzYdHZ2Dnt+6cCBA6mpqUlzc/PQ800f\nJiM+UzXqk3umahjPVA3nmaoyPFMFp84zVeNnYGAgZ5111oRe4wd6pgoAgFMz4u0/AOCj6/zzz8/h\nw4fHdM5KpTIm85x33nnp6+sr+j1FFUBB739GpL+/P2+88Uaqq6s/tM+I8NF2+PDhMbsdOjAwkD/8\n4Q9jNl97e3vx7ymqAAp6fzy9+OKLmT9/fl588cXxXRQwJkQVADAuDhw4kIMHDyb5v3eqPvGJT2T3\n7t25+OKLhz4N/sNEVAEfKRdccEF6e3vHdM6xekakuro6hw4dGpO5oISampph8TR79uxxXM0HJ6qA\nj5Te3t4x+4iJ/fv35xe/+EXuv//+MZlv+fLlYzIPcHyiCqCgPXv2ZO/evUmSN998M2+99VaeeOKJ\nzJgxY+iXQQOTk6gCKGjmzJlD8fTWW29l2rRpueaaa8Z5VcBY8OGfAKPk4x//uKCCjxBRBQBQgKgC\nAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoABRBQBQgKgCAChAVAEAFCCqAAAKEFUAAAWI\nKgCAAkQVAEABogoAoABRBQBQgKgCAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoABRBQBQ\ngKgCAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoABRBQBQgKgCAChAVAEAFCCqAAAKEFUA\nAAWIKgCAAkQVAEABogoAoABRBQBQgKgCAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoABR\nBQBQgKgCAChAVAEAFCCqAAAKEFUAAAWcNKo6Ojoya9aszJgxIytXrjxm+6pVqzJ37tzMnTs3V111\nVaqqqtLb2zu0vb+/P3Pnzs2NN95YduUAABPIiFHV39+fu+++Ox0dHdm1a1c2bNiQ3bt3D9vn3nvv\nzc6dO7Nz586sWLEizc3Nqa6uHtr+0EMPZfbs2alUKqNzBgAAE8CIUbVt27ZMnz499fX1mTJlShYu\nXJjNmzefcP/W1tYsWrRoaPzaa6/lySefzJ133pmBgYFyqwYAmGCqRtq4b9++TJs2bWhcV1eXrq6u\n4+575MiRbNmyJWvXrh167bvf/W5+8pOf5PDhwyecY8mSJamvr0+SVFdXZ86cOWlubk6SdHZ2Jskx\n40E7duxIkjQ2Nk6K8eA5nuz8XQ/j0xkP2rNnT5Jk5syZk2o8yPU4s+thPHw86MCBA0mSmpqaSTUe\n5Hqc2vUY/Lq7uzunojIwwltImzZtSkdHRx555JEkyfr169PV1ZXVq1cfs29bW1taW1uH3sl6/PHH\n89RTT+Xhhx9OZ2dnfvrTn+ZPf/rT8MkrlTN6B6tSqWT79u2nfdxE19TU5Hq8z5leD4arVCrDftiZ\nTJYvX37af0ZcD0ZSqVSyYMGC8V7GqGhvbz+jvy+ux3+drFtGvP1XW1ubnp6eoXFPT0/q6uqOu+/G\njRuH3fp77rnn8sc//jGXXnppFi1alD//+c+5/fbbT2vxAAAfFiNGVVNTU/bu3Zvu7u4cPXo0bW1t\nmT9//jH79fX15ZlnnslNN9009NoDDzyQnp6evPzyy9m4cWOuu+66/Pa3vy1/BgAAE8CIz1RVVVVl\nzZo1aWlpSX9/f5YuXZqGhoasW7cuSbJs2bIkyWOPPZaWlpZMnTr1hN/L//4DACazEaMqSa6//vpc\nf/31w14bjKlBixcvzuLFi0/4Pa699tpce+21Z7hEAICJzyeqAwAUIKoAAAoQVQAABYgqAIACRBUA\nQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBU\nAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIAC\nRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIA\nKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgq\nAIACRBUAQAGiCgCgAFEFAFCAqAIAKEBUAQAUIKoAAAoQVQAABYgqAIACRBUAQAGiCgCgAFEFAFDA\nSaOqo6Mjs2bNyowZM7Jy5cpjtq9atSpz587N3Llzc9VVV6Wqqiq9vb3p6enJF7/4xVxxxRW58sor\n8/Of/3xUTgAAYCIYMar6+/tz9913p6OjI7t27cqGDRuye/fuYfvce++92blzZ3bu3JkVK1akubk5\n1dXVmTJlSn72s5/lH//4R55//vk8/PDDxxwLADBZjBhV27Zty/Tp01NfX58pU6Zk4cKF2bx58wn3\nb21tzaJFi5Ikl1xySebMmZMkOffcc9PQ0JB//etfBZcOADBxVI20cd++fZk2bdrQuK6uLl1dXcfd\n98iRI9myZUvWrl17zLbu7u7s3Lkzn/vc547ZtmTJktTX1ydJqqurM2fOnDQ3NydJOjs7k+SY8aAd\nO3YkSRobGyfFePAcT3b+rofx6YwH7dmzJ0kyc+bMSTUe5Hqc2fUwHj4edODAgSRJTU3NpBoPcj1O\n7XoMft3d3Z1TURkYGBg40cZNmzalo6MjjzzySJJk/fr16erqyurVq4/Zt62tLa2trce8k/Xmm2+m\nubk5P/zhD/O1r31t+OSVSkaY/sSLrlSyffv20z5uomtqanI93udMrwfDVSqV4/6wMxksX778tP+M\nuB6MpFKpZMGCBeO9jFHR3t5+Rn9fXI//Olm3jHj7r7a2Nj09PUPjnp6e1NXVHXffjRs3Dt36G/Tu\nu+/m61//er7xjW8cE1QAAJPJiFHV1NSUvXv3pru7O0ePHk1bW1vmz59/zH59fX155plnctNNNw29\nNjAwkKVLl2b27Nm55557yq8cAGACGTGqqqqqsmbNmrS0tGT27Nm59dZb09DQkHXr1mXdunVD+z32\n2GNpaWnJ1KlTh17761//mvXr1+cvf/nL0EcudHR0jN6ZAACMoxEfVE+S66+/Ptdff/2w15YtWzZs\nvHjx4ixevHjYa1/4whfy3nvvFVgiAMDE5xPVAQAKEFUAAAWIKgCAAkQVAEABogoAoABRBQBQgKgC\nAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoABRBQBQgKgCAChAVAEAFCCqAAAKEFUAAAWI\nKgCAAkQVAEABogoAoABRBQBQgKgCAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoABRBQBQ\ngKgCAChAVAEAFCCqAAAKEFUAAAWIKgCAAkQVAEABogoAoABRBQBQgKgCAChAVAEAFCCqAAAKqAwM\nDAyM2+SVynhNDQBw2kbKpqoxXMcxxrHnAACKcvsPAKAAUQUAUICoAgAoQFQBABQgqgAAChBVAAAF\n/C8WUjRgkKEHPQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1f502c90>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We also print the values to the notebook, supporting the text in the paper manuscript:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print [np.median(dtm_rrmse_1k[np.isfinite(dtm_rrmse_1k)]), \n",
      "        np.median(sfm_rrmse_1k[np.isfinite(sfm_rrmse_1k)])]\n",
      "print [dtm_ci_1k, sfm_ci_1k]\n",
      "\n",
      "print [np.median(dtm_rrmse_2k[np.isfinite(dtm_rrmse_2k)]), \n",
      "        np.median(sfm_rrmse_2k[np.isfinite(sfm_rrmse_2k)])]\n",
      "print [dtm_ci_2k, sfm_ci_2k]\n",
      "\n",
      "print [np.median(dtm_rrmse_4k[np.isfinite(dtm_rrmse_4k)]), \n",
      "        np.median(sfm_rrmse_4k[np.isfinite(sfm_rrmse_4k)])]\n",
      "print [dtm_ci_1k, sfm_ci_1k]\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[0.77969796018681747, 0.77403867631975576]\n",
        "[0.0011072940946956766, 0.00097989156449773684]\n",
        "[0.78027668620035762, 0.75709104161134455]\n",
        "[0.0010044846549480679, 0.00077619529484007632]\n",
        "[0.78807302716479444, 0.7581707734866665]\n",
        "[0.0011072940946956766, 0.00097989156449773684]\n"
       ]
      }
     ],
     "prompt_number": 20
    }
   ],
   "metadata": {}
  }
 ]
}