{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Doing ordinary least squares and AR(1) modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The usual boilerplate\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=4, suppress=True)\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "# - set gray colormap and nearest neighbor interpolation by default\n",
    "plt.rcParams['image.cmap'] = 'gray'\n",
    "plt.rcParams['image.interpolation'] = 'nearest'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Imports from nipy\n",
    "import nipy\n",
    "from nipy.modalities.fmri.design import block_design, stack_designs\n",
    "from nipy.modalities.fmri import hrf\n",
    "from nipy.algorithms.statistics.api import (OLSModel, ARModel, make_recarray)\n",
    "from nipy.modalities.fmri.design_matrix import _cosine_drift as dft_basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we get the raw functional data and drop the first four scans, as dummy scans."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(64, 64, 30, 169)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_dummies = 4\n",
    "func_img = nipy.load_image('ds114_sub009_t2r1.nii')\n",
    "func_img = func_img[..., n_dummies:]\n",
    "func_img.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the times at which each remaining volume was acquired relative to the start of the first volume.  The first volume is the first of the four dummy scans we dropped:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TR = 2.5\n",
    "n_vols = func_img.shape[-1]\n",
    "volume_times = np.arange(n_vols) * TR + n_dummies * TR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We make the volume times into a structured array with field name `t`, so that nipy knows this is an array of time values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([(10.0,), (12.5,), (15.0,), (17.5,), (20.0,), (22.5,), (25.0,),\n",
       "       (27.5,), (30.0,), (32.5,), (35.0,), (37.5,), (40.0,), (42.5,),\n",
       "       (45.0,), (47.5,), (50.0,), (52.5,), (55.0,), (57.5,), (60.0,),\n",
       "       (62.5,), (65.0,), (67.5,), (70.0,), (72.5,), (75.0,), (77.5,),\n",
       "       (80.0,), (82.5,), (85.0,), (87.5,), (90.0,), (92.5,), (95.0,),\n",
       "       (97.5,), (100.0,), (102.5,), (105.0,), (107.5,), (110.0,), (112.5,),\n",
       "       (115.0,), (117.5,), (120.0,), (122.5,), (125.0,), (127.5,),\n",
       "       (130.0,), (132.5,), (135.0,), (137.5,), (140.0,), (142.5,),\n",
       "       (145.0,), (147.5,), (150.0,), (152.5,), (155.0,), (157.5,),\n",
       "       (160.0,), (162.5,), (165.0,), (167.5,), (170.0,), (172.5,),\n",
       "       (175.0,), (177.5,), (180.0,), (182.5,), (185.0,), (187.5,),\n",
       "       (190.0,), (192.5,), (195.0,), (197.5,), (200.0,), (202.5,),\n",
       "       (205.0,), (207.5,), (210.0,), (212.5,), (215.0,), (217.5,),\n",
       "       (220.0,), (222.5,), (225.0,), (227.5,), (230.0,), (232.5,),\n",
       "       (235.0,), (237.5,), (240.0,), (242.5,), (245.0,), (247.5,),\n",
       "       (250.0,), (252.5,), (255.0,), (257.5,), (260.0,), (262.5,),\n",
       "       (265.0,), (267.5,), (270.0,), (272.5,), (275.0,), (277.5,),\n",
       "       (280.0,), (282.5,), (285.0,), (287.5,), (290.0,), (292.5,),\n",
       "       (295.0,), (297.5,), (300.0,), (302.5,), (305.0,), (307.5,),\n",
       "       (310.0,), (312.5,), (315.0,), (317.5,), (320.0,), (322.5,),\n",
       "       (325.0,), (327.5,), (330.0,), (332.5,), (335.0,), (337.5,),\n",
       "       (340.0,), (342.5,), (345.0,), (347.5,), (350.0,), (352.5,),\n",
       "       (355.0,), (357.5,), (360.0,), (362.5,), (365.0,), (367.5,),\n",
       "       (370.0,), (372.5,), (375.0,), (377.5,), (380.0,), (382.5,),\n",
       "       (385.0,), (387.5,), (390.0,), (392.5,), (395.0,), (397.5,),\n",
       "       (400.0,), (402.5,), (405.0,), (407.5,), (410.0,), (412.5,),\n",
       "       (415.0,), (417.5,), (420.0,), (422.5,), (425.0,), (427.5,), (430.0,)], \n",
       "      dtype=[('t', '<f8')])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "volume_times_rec = make_recarray(volume_times, 't')\n",
    "volume_times_rec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we make the definition of a block design, in a format nipy understands.\n",
    "\n",
    "Here are the onsets of the blocks and their durations, in the standard OpenFMRI format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  10.,   30.,    1.],\n",
       "       [  70.,   30.,    1.],\n",
       "       [ 130.,   30.,    1.],\n",
       "       [ 190.,   30.,    1.],\n",
       "       [ 250.,   30.,    1.],\n",
       "       [ 310.,   30.,    1.],\n",
       "       [ 370.,   30.,    1.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "block_defs = np.loadtxt('ds114_sub009_t2r1_cond.txt')\n",
    "block_defs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make this into a structured array with fields `start` (onset time) and `stop` (offset time):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[(10.0, 40.0)],\n",
       "       [(70.0, 100.0)],\n",
       "       [(130.0, 160.0)],\n",
       "       [(190.0, 220.0)],\n",
       "       [(250.0, 280.0)],\n",
       "       [(310.0, 340.0)],\n",
       "       [(370.0, 400.0)]], \n",
       "      dtype=[('start', '<f8'), ('end', '<f8')])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "onsets, durations, amplitudes = block_defs.T\n",
    "block_spec_arr = np.column_stack((onsets, onsets + durations))\n",
    "block_spec = make_recarray(block_spec_arr, ('start', 'end'))\n",
    "block_spec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feed this block definition into the nipy utility for constructing the design matrix for a block design.\n",
    "\n",
    "We'll use the SPM HRF definition, but you can explore the `hrf` model for other HRF definitions we could have used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_exper, cons_exper = block_design(block_spec, \n",
    "                                   volume_times_rec,\n",
    "                                   hrfs=(hrf.spm,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because this is a single condition, this design is a single regressor - the block on/off step function convolved with the HRF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10de1c090>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYVOWV/z+nm61FUdAJTkBwwSXjmrjFvXFFM0aDUUEn\nmTGjIZkhy2R+UWPiI0wmk+SJyWQmjhGMMmYSxahRCSoalzbuYhRXVFxAQUOCoCwK9PL+/nir7Oqi\n6lY19d733tN9Ps/TD11Vt+/9curc97znvMsV5xyGYRhG/6YpawGGYRhG9lgwMAzDMCwYGIZhGBYM\nDMMwDCwYGIZhGFgwMAzDMAgUDETkKhFZLiLPVPn8LBF5uvDzoIjsHeK6hmEYRhhCZQazgBMSPn8N\nONI5ty/w78CVga5rGIZhBGBAiJM45x4UkbEJnz9a8vJRYFSI6xqGYRhhyGLM4FzgjgyuaxiGYVQh\nSGZQLyIyHjgHODzmdQ3DMIxkogUDEdkHmAlMcM6tSjjONksyDMPoJc45aeTvQ5aJpPCz6QciY4Cb\ngM85516tdSLnXK5+Lrnkksw1mKa+oymvukyTXk0hCJIZiMi1QCuwrYi8AVwCDAKcc24mcDEwArhc\nRARod84dFOLahmEYRuOEmk10Vo3PzwPOC3EtwzAMIzy2ArkOWltbs5awCaapPvKoCfKpyzTVRx41\nhUBC1ZtCISIub5oMwzDyjIjgcjSAbBiGYSjFgoFhGIbRt4LByy/D5z6XtYre8V//Bb/+ddYq6qer\nC049Ff70p6yV1M8TT8A//VPWKnrH9Olw221Zq6if9evhU5+CNWuyVlI/d98NF12UtYr80KeCwTXX\nwK9+BQsXZq2kPpyDn/wE/vM/s1ZSPw89BLfeCrNnZ62kfq68EmbM0BPANmyAn/4UfvazrJXUzx13\nwO23e9/QwuWX+87YunVZK8kHfSYYOAfXXw/jx8O112atpj4efRRaWmDZMnjppazV1Mfs2d7GWrKZ\n9na46SY48kjvHxq46y4YN877x/LlWaupD21+sXo13HMP7LefrgCWJn0mGPzxjyACP/qRDwYaJiTN\nng1nnQWTJukIYB0dcOONcMUV8OabviyXd+6+G3bbDS68UIeNwfvFOefAySfDb36TtZrarF0L8+bB\nrFnwyCPw5z9nrag2t94KRx0FX/6yHr9Imz4TDGbP9o3qJz4BAwfC449nrSiZzk644QY480wfEDQE\nsLY2GDPGN65nngnXXZe1otpcf73XeswxsHgxvPJK1oqSef99P1Zw2mndfpF35s6Fww6DsWP9uMEN\nN2StqDbF9uLUU+HBB2HFiqwVZU+fCAZdXd03vYiOm+iBB2DkSNh9dzjgAK/7iSeyVpVM8QYCOPvs\n/Nt4wwaYMwdOPx0GDNARwG6/HQ480PvGscfCq6/Ca69lrSqZ2bO9bUHHvbdypQ8AJ58MW24JJ57o\nM97+Tp8IBqtW+VR1r73861NO8Wlrnnn6aTi8sJG3iA7NCxZ0az7wQF93ff31bDUl8frrsO228NGP\n+tfabDxwoO9p33lntppqsWABHHGE//3YY+GZZ/w9mVdeeAH23BO22sq/1uAXMegTweC992Cbbbpf\n7703vPMOLF2anaZalGs++mi4997s9NRDqWaR/Gsut/Ghh/ognOfpj9r9YvBgOOQQ+MMfstWURLmN\nx4+H++/3pdv+TJ8JBltv3f26qcl/wffdl52mWpRrPvxwmD8fPvggO021KNec94aqXG9Li89oHnww\nO021KNdc9OOuruw0JeGcD67DhnW/p80vRo6EUaPgqaey05QH+mQwAH8TaXLIrbaCffaBhx/OTlMt\nKjVU996b34HvvuAXo0fDiBHw7LPZaUpi7VoYMsSPyRTRZmPIv+YY9IlgsHr1pl9usXeS14YqSXMe\n2bDB23LIkO73dtrJlwVefDE7XUloszHo01xJ7/77wxtv5HeKqTYbx6JPBIP33uuZpoKfpdPent+Z\nGJU0H310fktbRb1Ssi9icdwg75pLOeggWLQovwOcWv2ilAED/CK/trZMJNWkkuajjvKr6zduzEZT\nHugzwaA80hcbqnvuyUZTLSppPuQQPxNj9epsNCVRSS/os/GgQX4gOc8NVbnm1lY/INvRkYmkRPqK\nX4wY4dfPPPZYNpryQJ8NBgATJsDvfhdfTz1U0tzS4geS8zjNrZqNjzvOp9d5HPjuK37xkY/47Sny\nOEMnycZz5+Zz4FujX8SgTweDT3/a30ArV8bXVItqmidNyufCqGp6R470NeLbb4+vqRbVNJ9+Otxy\nix8HyRPO9R2/+NjHYLvt/OLKvJFk49mz8xnAYtCng8GwYb7n+tvfxtdUi2qaTz3V97Tfey++piQq\nDboVKd5EeaOa5lGj/MytvGVg69dDc7MflC/njDO8H+etpq3RL6rde3vt5Wf1PfJIfE15oE8HA4DJ\nk/PnkO3tvlc6dOimn22zja+33nxzfF1JJNl44kS/02bexjpq+UXeetpJeseM8b3tu+6Kq6kWSZon\nTfLbPLS3x9VUi2qaRfLZXsSiTwSD1as3nR1Q5KST/J4/edrLfs0a3wORKk8szWOPqtIMjCIjRvjZ\nI3PmxNVUiyTNp53m9+DP0172SXpBn1/stBPsskv+BpKT2oszz/Qb7eVxsD5t+kQwSOqdtLTAZz7j\nt13OC0l6Af72b/1q5Dw9pKeW5rPPhv/5n3zVW5M0b7ed309n1qy4mpKoZePTT/djM3naZqWW5r/7\nO7jssnyt90nSvOuusOOO/hkY/Y0gwUBErhKR5SLyTMIx/y0ii0RkgYjsF+K6RWo55PTp3iHzsv9+\nLb1Dh8K0afClL+Wnca2noXIOfvGLeJpqUUvzD37gfeOtt+JpSqKW3pEjYepU+NrX4mmqRS3N553n\nd1695ZZ4mpJwLjkzALj0UvjGN/I3bpc2oTKDWcAJ1T4UkROBXZxzuwJTgKD99FoOOWYMfOc7vnHN\nQw+lll7wz+z94IP89FxraW5uhpkz4dvfzk9JrpbmvfaCL34xP41rPX5x0UV+a4q8lORqaR482D9y\n9KtfzceY0rp1XtPAgdWPOfxwn51/61vxdOWBIMHAOfcgkLSm8xTgl4VjHwO2FpGRIa4N9d1EU6f6\nB4d8/vPZO2U9eouN6wUXwP/+b/ZBLGnWSJF99oEpU/x87axLXBs2+F0oW1qSj/vOd3zjOnVq9msl\n6rHxkCG+cT333HzMkqtH85FH+q24TzrJP2AoS+q598BnjXPn+s5N3gbA0yLWmMEo4M2S18sK7zVM\nZ6eP9sW9yasxYIAfyNpiC99oXXyxfyTiSy/5IBGTeh1yv/38NNMf/xiOP97X5B9/HJYsiV8+qlfz\nd7/rHyV4xBH+35tu8o3tO++kr7GUYiNVbZC+SEuLn0q4YoW397/9m9/6YdEiP9UzJvXaePx4nxmc\nf75fSzNzpp8k8eab8TsN9Wq+/HL/3IADD/SZ2K23wvPPw7vvpq+xlHr1Dh/u77UFC/zDp77/fb9m\n4tVX8ze9NxQDah8Sn2nTpn34e2trK62trVWPXbPGP62oqY6wNnSo71U98oiP+t/9rn8Y/RZb+G0g\nYlGvQ4IPXPPn+17g3XfD1Vf7YHD++f4nFvVqFvHZwQkneM2zZvne4JIl/lkCO++culSgdzYePtxP\nM73vPj/D6OKL/SDtzjvH3bysN5o/+UnfUN14o/eLK67w+3D95CfwhS+kq7OUejU3NcE3v+mD1803\nw89/7jeze+st38Buu236WqF3Nt5+e99OzJvnHzB0/vk+4B52mH+yYpa0tbXRFnpPFedckB9gLPBM\nlc+uAM4sef0iMLLKsa43LFni3OjRvfqTHnR2OjdihHPLlm3+OXrL977n3IUXbv7f33abc0cfHU5P\nPeyxh3PPPbf5fz95snNXXx1OTy2eeMK5/fbb/L9//33nttjCubVrw2mqxde/7tyll27+319zjXNn\nnBFOTz1st51zy5dv/t8fd5xzc+aE01OLefP8NTeXt992bvhw327kiUK72VAbHrJMJIWfSswBPg8g\nIp8E3nXOLQ9x0d5E+ko0NflIH/OBJ41qPuQQn8LGrGU2qvnww3XZuKUF9t3X2zkWoWwcq1SUtH1G\nvWjzi+2391nM88+H05QXQk0tvRZ4GNhNRN4QkXNEZIqIfBHAOXc78LqIvALMAP4pxHWh8S8X9Dnk\n8OF+QU/MJzOFuOlj7lMTyi80ad5pJ99Ax3ou9fr1vjNVafuMetFmY4jfXsQi1Gyis5xzH3XODXbO\njXHOzXLOzXDOzSw5Zqpzbpxzbl/n3JMhrgthvtwjjtAVDCCu5o4Of+NvueXmn2OvveAvf4HlQfLB\n2tQzy6UWsf2iUc0icTWHsPHBB/uxpFgzubTdezFRvwI5xJf7iU/4BWmxppxq650UF+nUmpmTRFOT\nf47AQw+F05VECBsfeig8+mi8rQm0+UUIvUOH+o7C/PlhNNVCY8YYCwsG+DR3//3j7VYY8qaPUR8O\noRf0NVTbbgs77BBvpll/DAagT/Ouu/pM+Y03wmjKC30iGCQtLa+X2A7ZqOYddvCDnIsWhdGURH+1\nMejTvM8+frr0ihVhNCXRX20s0jfHDdQHgxB1S/CZQaweYEjNTz/d+HlqEaoHGEsvhNWsKTNobvaz\noJ57LoymJDTaWGN7EQv1wSCUQw4dqmsQC+JpDqW3pcVPh+3sbPxctdBm440bvW222KLxc2nzC7v3\n8oEFgwJDhsTZfqCrC9aurb19Rj3E0hyqNyXix2diPG4ylObYNm5kkL7I4MG6/CKWXtDXXsTEgkGB\nWI3UmjW+V9Hc3Pi5YmkOZWPQp1mbXvANlSbNsfSCPr+IiQWDAhpvetNcHfOLMOdKIpTmAQP8rLgY\nU3i1+UVM+kQwCDGjIVaqGkovmOYkzC/CnCuJ0JrTblyL22do8ouY9IlgoClV7c/lANCnWZteMM3V\nWL/ej8kMGdL4uWKWtmKhPhisW9fYNglFYqV9ofSCaa6Gc/4ZFeYX6aJNsza9sVEfDDZsaGyjrCKx\n0r5QesE0V6O93Q/Q1/OMi1qYjaujTbM2vbGxYFAgVtoX0iFNc2W06YXwDZVp3hSNfhET1cHAOb9Y\nZ9Cgxs8V6wYKpRdMczW06YWwmmM1VNo0a/SLmKgOBu3tflpaiHLAgAF+QVja09u09aZAn2aN5QBt\nNgZ9mrXpjY3qYBDyyy3OMtDkkLFWQWrTrLEcYAFMn1/YmEGOCJn2gZUwqqFNcxp6094qXFvJBfRp\n1ubHsVEdDEJGerBUtRraNIfU29zss0YrH26KNs3a9MZGdTDYuDF8MEg79QupOVY5QJvmNPwiRq9V\nk41BX5lIo41jojoYbNgQtkwUa8xAU2oN+jT3d7+w8mFlQtt448Y4TxqMhfpgoK0HqDFV1aa5v/uF\nxkFvbZM3mpv9T3t7mPPlAdXBwMpEViaqhPlF+nq7unxDOHBgmPNpszH0vVKR6mDQ38sBViaqTH/3\nixiZTHu71xviYTygr0wEfW8VcpBgICITRORFEXlZRC6o8Pm2InKHiCwQkWdF5B9CXLe/lwOsTFSZ\n/u4X2kouoFNzX5tR1HAwEJEm4DLgBGBPYLKI7FF22FRggXNuP2A88GMRGdDotW2dQZwdQEOXAzTZ\nGPSVMDT2srX6hQWDnhwELHLOLXHOtQOzgVPKjvkTUHzq71bAO865hmdup9E7sVWQPdm40QeCEFt+\ngD4bg84ykabgBfpWIEPfW4Uc4hYfBbxZ8npp4b1SrgT2FJG3gKeBrwW4rsq0L3TJJe3pbf3dxqBP\ns8aSi0bNfS0zaLhUUyffAp52zo0XkV2A34vIPs65tZUOnjZt2oe/t7a20traWvGkGtO+kJqbmrqn\nt4W0Qyn93cagr6dtZaLKbNwILS3hzpdlMGhra6OtrS3oOUMEg2XAmJLXowvvlXIY8D0A59yrIvI6\nsAfwRKUTlgaDJNKI9NpS1aLmtIKB2Vhf46oteEE8v9hmm3Dny3JqaXknefr06Q2fM0SZaD4wTkTG\nisggYBIwp+yYhcCxACIyEtgNeK3RC1uqmr5mbXrBNA8aZOXDSmj0i5g0nBk45zpFZCpwFz64XOWc\nWygiU/zHbibwfWCWiDwNCHC+c25lo9e2Ekb6mrXpBZ291pCam5q6A0JIO5QS2i+07VoKNmZQEefc\nPGD3svdmlPy+Ajg5xLVKsRJG+pq16QWvediwcOfTViaCbjunFQy0+oU2zTFRvQI5dM8nVu9Ek2Zt\nesE0Q/oBTGOZSKNfxER1MNA4o0GbZm16Ib1edppoq2drLBNp9OWYqA8G2non2jRr0wt6NWtqqLTa\nWJvmmKgOBmmkfZoGCiF9zdr0gr5yQFcXdHaG2/IDbCypEhp9OSaqg4G23hTo06xNL+jTXCy5hNoB\nFGyWWSW0+UVs1AcDbWmfNs3a9IK+XmvoRgps/UklNPpyTFQHgzQGsdK86bu6/IPVQ5YDYpSJNNkY\n9A1uprEeQFsAi1Um0ubLMVEdDLRFeo3lgNA2HjjQ18e7usKdsxxtfhFaL8TxZU02Bn1+ERsLBiXY\nTb8poTWL+ICoSbO2XjZYmagS2tqL2KgOBtoGsULrBdNcCY1+obGToMnGoM8vYqM6GKTRO9E0HQ9M\ncyW09VrTyhg1TTkeMMCXDzs7w52zHG3tRWwsGJSgreQCprkSVibSZ2MRfZotM8gR2tI+K7l4tGnW\nWCaKkc1oG+fQ5hexUR0MtKV9VnLxaNOssWHVViYCfaumrUyUI7SlfVZy8WjTrK2RAn02hnQ1d3T4\nfwcEfNCvZQY5QlvaZyUXjzbNVibaFG1lIo1+HBvVwUBb2mclF482zRobVisT9USjH8dGdTBIY3pb\ncQfJNNA4n1yb5uLq5tDlAI0Nq5WJutFo49ioDgahe1RpT2/TOoVQk2bb8sNjftETbXqzQH0w0HQT\nadML+jSnqde5sOctoq3+DvrGObT5cRaoDQZplAMg3TpgWjeQthKGNhsPGOAzjeKMlNBoq7+DvlXT\n2vw4C9QGgzTKAWCpajnaNKehF/RmM2lhfmGZQW5I4wYCu+nL0aZZq19YmUhvmSit8mFsggQDEZkg\nIi+KyMsickGVY1pF5CkReU5E7mv0mmnMG4b0Sxhp3PSaHm4D+mwM6c+B11RyAX3TYdPwi+Zm/9Pe\nHva8WdFwxV1EmoDLgGOAt4D5InKrc+7FkmO2Bv4HON45t0xEtmv0ulp7gJp62aBPc5p+keYc+G23\nDXtO84uepN1epNEBiU2IzOAgYJFzbolzrh2YDZxSdsxZwE3OuWUAzrkVjV7UgoHHbvqeaPULbfVs\nbfP2NfpFbEIEg1HAmyWvlxbeK2U3YISI3Cci80Xkc41e1MpEHisT9cTKRJ4Yq7zNL/rWjKLAEzMT\nr/MJ4GhgKPCIiDzinHul0sHTpk378PfW1lZaW1s3OUZjpN+wAUaMCHtOywx6orVMpMnGzunTrLG9\nSKKtrY22trag5wwRDJYBY0pejy68V8pSYIVzbj2wXkT+AOwL1AwG1dD45Wq7gUCfZq1+oalM1NkJ\nTU1+8DQkgwfDunVhz1lEo18kUd5Jnj59esPnDFEmmg+ME5GxIjIImATMKTvmVuBwEWkWkS2Ag4GF\njVxUazkgrZs+reltaZWJNNkYdJaJNAUvML/ImoYzA+dcp4hMBe7CB5ernHMLRWSK/9jNdM69KCJ3\nAs8AncBM59wLjVzXygGe0ultaTi7tpWm5hcebXrBNGdNkDED59w8YPey92aUvb4UuDTE9SCd3hTo\n2zkR0pveVty9NfSWH5ptnAba9trXamNtmmOjegWyxm0HNGnWphf0ata0OlarjbVpjo3qYJBGpE9z\nqpg2zdr0gmku0tQEAwf6HnFozMbd9KWppWqDgca0T5tmbXpBr2ZNvVatNtamOTZqg4HGtE+bZm16\nQa9mTQ2VVhtr0xwb1cFAW9qnTbM2vWCaSzG/6Eaj5tioDQZppdaDBqVTZwV9mrXpBdNcivlFNxo1\nx0ZtMNCWWoM+zdr0gmkuxfyiG42aY2PBoIzBg9OL9No0a9MLprkUjQPIGm1smUHGpJn2aZs1kpZm\nbXrBNJeSVgkjrcFYrTa2zCBjtPWmQJ9mbXpBn+a0dgAF84tSNGqOjQWDMixV7UabXtCnuaPD7y3V\nlMKdmJZmKxN1Y2WiHKAx7dOmOS29xa2POzrCn9ts3E1amq1M1I2ViXKAtrTPygE90aZZm17Qp9nK\nRNmiNhhoS1XTeiAI6CsHgD7N2vSCPs1pllw0ao6N2mCgLVVNSy+Y5lLML7rRprm5Gbq6urdOD4k2\nv8gC1cFAU9pn5YCeaNOsTS/o0yyiT7OViXJAWmlfmkv407rpTXM35hfdaCsTgT7Nth1FDtC2C2Ga\n5QDT3I35RTfaykSgT7NlBjlAW9pn5YCeaNOsTS+Y5lK0+UUWqA0G69f77WNDk1bal5ZeMM2lmF90\nk+ZiRG2atflFFqgNBtoivfWmeqJNsza9kF7JZf16XZqLM5QGDgx7XrDMIBekdRMVnTH0g8StYe2J\nNs3a9IJpLlIcLxAJe16wYJAL0kpVm5v9T+itEtJMrdPcndI0e7TpBSsTFdHox1kQJBiIyAQReVFE\nXhaRCxKOO1BE2kVkYqPXTDNVTSPaa9ML+jR3dvqfAQPCnhd02tjKRB5tfpwVDQcDEWkCLgNOAPYE\nJovIHlWO+wFwZ6PXhPRT1TR6J5r0gj7NRb1plAOKPUArH+r1izSw7Sh6chCwyDm3xDnXDswGTqlw\n3FeAG4E/N3rBNDd9g3R6J9r0gj7NaZYDittMp1E+1NZQaQtg2vw4K0IEg1HAmyWvlxbe+xAR+Shw\nqnPu50DD/baODt/7S6McAOmVidKss6ZVDtCkOc1yAOjzizTLRJo0a/PjrEipOd2EnwKlYwmJAWHa\ntGkf/t7a2kpra2uPz9OM9GCpahFtms0vemJlIo82vfXQ1tZGW1tb0HOGCAbLgDElr0cX3ivlAGC2\niAiwHXCiiLQ75+ZUOmFpMKhE2je9lYk82jSbX/QkzRlQ5hee0gc1pVWpqER5J3n69OkNnzOE/PnA\nOBEZC7wNTAImlx7gnNu5+LuIzAJ+Vy0Q1EOaaR/oKwdYmcij1S+GDQt7ziJpZgbmF90UNccMBmnQ\nsHznXKeITAXuwo9BXOWcWygiU/zHbmb5nzR6zRg9QE2pqtYeoCYbgz7NGqfDaisTQbdfDB2a3jVi\nECSWOefmAbuXvTejyrFfaPR6MWrDaaSqLS1hz1lEa2141aqw59TqF5oCLliZqJy+MoiscgWyxi9X\nYw9Qm+Y0yxdgfgHpT+vWZmOwYJApadcA0+hR2Q6gPUnLxtrKRNrGZWJM69ZkY+g7W1KoDAYaI722\n3hTo02x+0RNtYxxgZaIssWBQAa03fcitErSWA6xM1I02vWCas0RlMNCY9qWpuanJp+3t7eHO2d7e\nvYNrGliZyKOtTBQj4GqyMViZKFM0RnptmrXpBdNcjtaAq8nGYJlBpti2A5sSWrM2vaCz1xpjq4SQ\n5UOtfqFNcxaoDAYx0j5tqyBDa9amF3T2WtMuHzY3hy0fahyX0ejLWaAyGGhM+7Rp1qYXTHMlQpeK\nNAZcjX6RBRYMKmCpqj69YGWiSmjsJGi0sZWJMkJj70SbZm16wTRXQmMw0GZjKxNliMa6pTbN2vSC\nzoYqbTuHLhPF0KvNxlYmyhBtdVbQp1mbXjDNlQjdUMXIZLTZ2NYZZIi2Gwj0adamF/T1AJ3zjcig\nQeHOWY75hU7NWaAyGGh9iIkmzdr0gr76+8aNMHCgnwKaFlYm0unLWaAyGGhM+7Rp1qYX9GlOWy9Y\nmQj0+UVWWDCoQOgbqKvLL/yxckA3VibSGQy0+oU2zVmgMhhoS/uKzigS7pzlaCsTDRrkA2RXV7hz\naisTpW1jsDIR6GsvskJlMNCW9sXoAWrTLOLr5aG3SjC/6InGzMDKRNlgwaAC2m4gMM1gZaJKaBwz\nsDJRNqgMBtocMm29YJrB/KIS2rKZgQPD77SqzS+yQmUwsFWbm2Ka9ZUDYthYW/bV1JRO+VBTe5EV\naoOBpkivsRygVbOViXqiMZvRZmfLDEoQkQki8qKIvCwiF1T4/CwRebrw86CI7N3I9bQNYsW66U2z\n+UU52rIvCKs51rRuywwAEWkCLgNOAPYEJovIHmWHvQYc6ZzbF/h34MpGrqntYRWxphCaZl27U8aw\nsbbsC8JqjjGt23Yt7eYgYJFzbolzrh2YDZxSeoBz7lHn3HuFl48Coxq5oLa0T2M5QKtmKxP1RGM2\nk0YwSBMrE3UzCniz5PVSkhv7c4E7GrmgxoFCTak16NPc2el/BgwIc75KaLWxpuwLwtpZmx9nSYq3\nzqaIyHjgHODwpOOmTZv24e+tra20trb2+FzbikKN5QBtmm2Vd2UGD4bVq8OdT1uZSJsf10tbWxtt\nbW1BzxkiGCwDxpS8Hl14rwcisg8wE5jgnFuVdMLSYFCJGHOdOzr84FOIHSU1pqraNGvTC1YmqkZ/\n94t6KO8kT58+veFzhigTzQfGichYERkETALmlB4gImOAm4DPOedebeRiHR3+3zTLASKWqmrTbOsi\nKmNlIl16s6ThJtU51ykiU4G78MHlKufcQhGZ4j92M4GLgRHA5SIiQLtz7qDNuV4MZ4TuaB+igdE4\nN1ubZm16QadmbT1tjTbOiiD9a+fcPGD3svdmlPx+HnBeiGvF6AFC+FRVW91Sm+YYjVQa5cMYNu7P\nq6a16c0SdSuQY9z0YKmqNs0xbnqt5UNtmYFGG/eFMpEFgypo67X293KA+UVl+ntpS6ONs0JdMIgx\nVQxseps2zbHHkkJgD7epTMhsJpaNQz+oKQvUBQMrE1XGykRxOgna/EJrxqjJxmk8qCkLLBhUob+n\nqto0m19UJqTejg7/nIE0p3WDPhtD3ygVqQsGMctEoXonsUouITMDbZpjlok0+UV/z740thdZoS4Y\nxCwTaeqdhNTrnD7NMRsqTX7R33vZGtuLrLBgUAVtDhm6HNDUBM3NYc5XDW02Bn2atekF05wV6oJB\nrHJA6FRzHQqqAAASRElEQVRV02Cs2bg62jRr0wumOSvUBQNbgVwZbXrBNJtfVKa/a84KlcFAW9oX\nQ/OAAX6ec2dn4+cyG1dHm2ZtesE0Z4W6YKAx7YuhOeRWCWbj6mjTrE0vmOasUBcMYqV9W2wB778f\n5lzaNGvTC6a5Gi0tvkEMsTrWbFydkJqzQmUwiBHphw2D996rfVw9aNOsTS+Y5mo0N/uAsHZt4+cy\nG1cnpOassGBQha23Dve4QG2atekF05xEKM2xSi792cZZoi4YxHTIUJFem2ZtesE0JxFKc8yGtb/a\nOEvUBYNYNcCQX642zdr0gj7NzvkBR00NlTYbg07NWaEyGGi6gUBf3TKW3i23hHXrdE2HDWXj9nY/\nHTjEE9Nqoc0vNI4ZWDDIgFhpX0iH1JaqxtLb3AxDh8KaNY2fy2xcHW2ai3qda/xcGtuLrFAXDF55\nBcaMSf86oW6gN96A4cN1NVSvvAJjxzZ+nnoIqVmTX8TSC/o0Dxrks6YPPmjsPM7p84ssURUMNmyA\nJ5+EQw5J/1qhZgc88AAccYRfFJY2oTXHIITmFStg2TLYZ58wmpLorzYGfZpfftlnnjvsEEZTEjab\nKDJPPAG77w5bbZX+tUJF+gcfhMMPb/w89RBCc2cnPPIIHHZYGE21CKH5oYd8ByHtXVah//rF2rXw\nwgtwwAFhNNUihGZtNs4aVcEg5pe75Za+3tjR0dh5tDnkc8/B9tvDX/1VGE216I83vXP6ND/2GHz8\n43Fm5kAYzQ88oMvGWRMkGIjIBBF5UUReFpELqhzz3yKySEQWiMh+m3OdmDeQiM9AGkn9Vq2CxYth\nv8363/aeEINYMW0M+jQXfaKRwc3XXvNZTKxxGW02hnCaY5a1+n0wEJEm4DLgBGBPYLKI7FF2zInA\nLs65XYEpwBW9vU5Xly8HxHTIRr/ghx+Ggw9O/5mxRbT1sqFxze+/D88+CwcdFE5TEgMH+t5xI9s7\nFG0cYxwJ+qdfvP02rFwJH/tYOE1JWDDwHAQscs4tcc61A7OBU8qOOQX4JYBz7jFgaxEZ2ZuLLFwI\nI0bAX/91AMV10mjvJOaAGzTukM7p0/z447D33n6jsFg0qlmbjTs6fJko1jgSNK75oYe83hjrOMB3\nEDo7dW9jHaLPOgp4s+T1UnyASDpmWeG95fVeJHbPBDZvhsCyZXDffb4U8JvfwJVXpqOtEpujt7MT\nFizwWczixf71zjunIq8im6N5zRr4wx98RnDPPTr8YvFi7xdLlsBtt8FXvpKKtIpsjt6ODj9h49FH\n4aWXfElr+PB09FViczSvWgVtbb7jeNttcEp5lzRFRLo11zveVpz6ev/98OabPpP52c/S1ZlEpAJG\n75g2bdqHv7e2ttLa2soBB8SrvRfpTe9k6VL4x3/0N9D48bDHHnDhhXDUUelqLKW3vam5c2HKFJ8B\nHXmkn4J3/fXxyhfgNb/1Vn3HbtgA3/42XHGFLwvtvz8cfzx8/vPpaiynN3Z+5RX4+7/3/44fD7vt\nBj/8YZxpsEV66xfXXQdf/Sp89KM+0O6wA1x1VXr6KtEbzevWwTe+AddeC4ce6tuJ006DL3whXY3l\nFDXXEwyeecb7xV/+4v1il118Saurq75spq2tjba2toY1lxIiGCwDSpd1jC68V37MDjWO+ZDSYFBk\n//03W99mU69DPvAAfPazvrd3223xxgjK2WorX8uux6EuvBBmz/Y3/pFHxtFXiXptvGIFHHss7LST\n72Vvt13q0qpSr+a5c+Gcc+Dii2Hq1Hgli3Lq1escfPGL3p/nzcvmniuy9dZ+wWYt3njD+8Whh/qs\nfNiw9LVVo147/+pX8C//Apde6jsym9P5KnaSi0yfPr33JykjRLM1HxgnImOBt4FJwOSyY+YA/wxc\nLyKfBN51ztVdIsqKer7c99/3Ef4Xv4CTT46jqxrNzb52vmaN116N3//eB4Knnoqb+lei3nGZf/1X\n30v92c/iZi6VqEfzO+/AuefCLbfErbVXoqjXuWTbzZ7tx2CefDLuGEwl6rGxc/DlL8PZZ8Mll8TR\nlUQ97cUbb8DXv+5LQ3vuGUdXvTQcDJxznSIyFbgLPyB9lXNuoYhM8R+7mc6520XkJBF5BVgHnNPo\ndWNQz5c7fTp88pPZB4IiRc3VgsEHH8CXvgSXX559IID6bHzPPb4W/Pzz2QcCqE/zN78JZ5yRfSAA\nvxVKc7P/7qs18itX+lLLzTdnHwigPhvfcIMfg7n55jiaalFLs3M+Q/za1/IXCCDQmIFzbh6we9l7\nM8peTw1xrZhsvbUflKrGwoUwa5YfyMwLtRzyhz/0q0hPOimepiRq6e3s9L2/yy7zCwHzQC3NDz/s\ns6/nn4+nqRZFzdUa+osvhokTfccmD9Sy8bp1vtRyww1+L6M8UEvz3Ll+7OiGG+Jp6g25HEDOC8OG\n+Z5HNa6+Gs47D0b2apJsuiTNwujs9LOb7r47rqYkas0aaWvzYyF5ybygtuYrr/S97Czr1+UUNVea\nmv3++/DrX8OLL8bXVY1aNv7d7/yU4kMPjaepFsOGJWueMQMuuijebrW9RdV2FLFJivRdXX7mzeTy\n0ZGMSdL84IN+pkOshTj1UKs3NXu2Lhtv2AC33gpnnhlXUy2SNN9+u88Wt98+rqYk+ppfrFzpB+Zj\nTnftLZYZJJD05T78sP98r73iaqpFkubrrsvfDVTsTVUa3Ny4EX77W78OIk8k2fiOO2Dfff20zDyR\npFlbw/ruu37NxjXXxNVUi623huVVpsXcdBOccEKcTTY3F8sMEqh1A02aFFdPPVSbhdHe7h0ybz3W\ngQN9zXfduk0/u+su+Ju/ibMFcW9ImumizS9Wr/bjGxMnxteURJKNb74ZjjkmecZcFmhsL0qxYJBA\ntS+3o8MPAuWtYYXqmu++G3bdFXbcMbqkmlTTnMdMBqrrXbvWZwannRZfUy2qab71Vr8wMg8zy0pp\nafH32caNm3523XX5bFir2fjtt/103RNPjK+pN1gwSKDal/vYY34gbty4+JpqUU3znDn5bKSgsuaO\nDl/LzluPFarb+J57/MroLBfEVUObXxS3dyjXvGaNL9F+6lPZ6Eqimo1vuw0mTPABLs9YMEig2oyG\ne++F446Lr6cekjQfe2x8PfVQSfNTT8Ho0fka1CzSV2zc1eVr75o0P/CAD7hDh2ajKQmN7UUpFgwS\n2Gor3xPp6ur5/r33wtFHZ6OpFpV6J0uX+hWxe++djaZaVNKszcagT/Ozz8K228KoUdloqkVf8Avn\n8q25FAsGCRS3dyjdu/6DD2D+/Pg7ZdZLpYG3++7zm2FltTdOLSppzvMNVLq9Q5Hly33Q/fjHs9OV\nhDYbgz7NlfQuXOjbkDyO1ZWT0+YhP5RH+4cf9lMH8zpFTFtvCjbVvHGjt3PMHV97w5Ahvqa9fn33\ne21tfsO/rDYprEVf8It33vEreA88MDtNSWi0cSkWDGpQ/gXn/cst16shTS3X/NhjfgvwbbbJTlMt\ntPtFR4evv5dsfJk7yjXff7/PyAcOzE5TEpWem553vyjFgkENtN/0r73mnXO33bLTVAttNgZ9msv1\n/vGP/oE19T6IJQu02bj8uemdnT5jHD8+U1l1Y8GgBjvs4J/0BD5Nff55OOSQbDUlMXKk3/t/zRr/\net48v0AnD7t9VqPUxtCtOc+Uan71Vd8A5HEnyiKjR8Prr/vFh6DPxs7BnXfq0jx/vp+CHvNRvY1g\nwaAGn/lM9y6DN9zgd/scMiRbTUlsuaWvtc+Z419fd53fSjnPfPrT/kZft843WK++mv/e1MSJ3X5x\n3XVw+un5HaAHP2No3Djfu3ZOh1985jN+O5LOTv8EQYj/tMPeUtpeaLBxKTl233xw8sl+MHPFCv9Y\nvbPPzlpRbSZN8svfFy/2vZQTTshaUTLbbeezrblzvY1PPz2/deEin/2s1/vBB37HT01+8eSTPkM4\n+OCsFSWz665+j6f77+++9/Kc4YLfleD66/0kiOuvh7POylpR/YgrnR+XA0TE5U3TmWf6XtWMGf55\nvXnZP70aq1f7dHXKFF8u+vnPs1ZUm1mzfDbz0kt+C+g8PBSmFscc47dQ/vWvfTaT94Zq2TK/1uSs\ns/z2E9/9btaKavOjH/mttW+/3QeFPI99FdlnH7876Z13+ifHxUBEcM415oHOuVz9eEn54re/da6p\nybkpU7JWUj8TJ3rNDzyQtZL6WLXKuZYW53bc0bmurqzV1MfMmd7G3/521krq58gjveYXXshaSX0s\nXuxcc7NzBxyQtZL6+d73vI1/+tN41yy0mw21vVYmqoMTT/SzBDSUAopMnuyzgzw9/COJbbaB44/3\nvda897CLTJzoy1maSgGTJ/u6e56eaZHE2LH+6Wua7r0zz/TrTTSNF4CViermnXf80n0tOOcf2Tli\nRNZK6mf1ar+ZV97HC0pZuVKXjbu6/HTNvO1SmsS77/rVvXkeoC8ndnsRokxkwcAwDEM5IYKBolhr\nGIZhpIUFA8MwDKOxYCAiw0XkLhF5SUTuFJFNHkQnIqNF5F4ReV5EnhWRrzZyTcMwDCM8jWYGFwJ3\nO+d2B+4FvlXhmA7gG865PYFDgH8WkT0avG5U2traspawCaapPvKoCfKpyzTVRx41haDRYHAKcE3h\n92uAU8sPcM79yTm3oPD7WmAhkNPHaVQmj1++aaqPPGqCfOoyTfWRR00haDQYfMQ5txx8ow98JOlg\nEdkR2A94rMHrGoZhGAGp+SgOEfk9MLL0LcAB36lweNU5oSKyJXAj8LVChmAYhmHkhIbWGYjIQqDV\nObdcRLYH7nPObbK2UUQGAHOBO5xz/1XjnLbIwDAMo5c0us6g0Yf0zQH+Afgh8PfArVWOuxp4oVYg\ngMb/Q4ZhGEbvaTQzGAH8BtgBWAKc4Zx7V0T+GrjSOfe3InIY8AfgWXwZyQEXOefmNazeMAzDCELu\ntqMwDMMw4pObFcgiMkFEXhSRl0Xkggx1LBaRp0XkKRF5vPBezcV1Kei4SkSWi8gzJe9V1SEi3xKR\nRSKyUESOj6jpEhFZKiJPFn4mRNZUcVFjlraqoOkrhfczs5WIDBaRxwp+/byI/Efh/SztVE1Tpj5V\nuE5T4dpzCq8zvfdKND1VoimsnRrdAzvEDz4ovQKMBQYCC4A9MtLyGjC87L0fAucXfr8A+EEEHYfj\np+E+U0sH8DfAU/gxoB0LtpRImi7BLyosP/ZjkTRtD+xX+H1L4CVgjyxtlaApa1ttUfi3GXgUOCwH\nPlVJU6Z2KlzrX4BfAXMKrzO1UxVNQe2Ul8zgIGCRc26Jc64dmI1f0JYFwqYZU83FdaFxzj0IrKpT\nx6eB2c65DufcYmAR3qYxNIG3WTmnRNJUaVHjaDK0VRVNxYWWWdrq/cKvg/E+vorsfaqSJsjQTiIy\nGjgJ+EXZtTOzUxVNENBOeQkGo4A3S14vJbtVyg74vYjMF5FzC++NdL1YXJci1Rb5ldtvGXHtN1VE\nFojIL0rS5+iapHtR46NU/86i6pJNF1pmZqtimQH4E9DmnHuBjO1URRNk61P/CXyTnuumsvanSpog\noJ3yEgzyxGHOuU/go/A/i8gRbPoF5GXUPQ86Lgd2ds7th7+hf5yFCNl0UWPm31kFTZnayjnX5Zz7\nOD5zOkJEWsnYTmWajhSRo8jQTiLyKWB5IbNLmuYezU4JmoLaKS/BYBkwpuT16MJ70XHOvV349y/A\nLfj0armIjAQQv7juz1loS9CxDD+9t0g0+znn/uIKhUrgSrrT0WiaxC9qvBH4P+dcca1LpraqpCkP\ntiroWA3cDhxATnyqoOk24ICM7XQY8GkReQ24DjhaRP4P+FOGdqqk6Zeh7ZSXYDAfGCciY0VkEDAJ\nv6AtKiKyRaE3h4gMBY7Hr48oLq6D5MV1wSXRsydQTcccYJKIDBKRnYBxwOMxNBVujCITgecy0FRp\nUWPWttpEU5a2EpHtimUEEWkBjsMPMmZmpyqaFmRpJ+fcRc65Mc65nfHt0L3Ouc8BvyMjO1XR9Png\ndkpj1HszR8on4GddLAIuzEjDTviZTE/hg8CFhfdHAHcX9N0FbBNBy7XAW8AG4A3gHGB4NR347cNf\nwQ9WHh9R0y+BZwp2uwVfW42p6TCgs+R7e7LgS1W/s7R1JWjKzFbA3gUdTwFPA/+vlm9nqClTnyq5\n1lF0z9zJzE4JmoLayRadGYZhGLkpExmGYRgZYsHAMAzDsGBgGIZhWDAwDMMwsGBgGIZhYMHAMAzD\nwIKBYRiGgQUDwzAMA/j/kqL9VL5wUakAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b23f310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(volume_times, X_exper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nipy made the contrast for this part of the design.  In this case it's just a single value == 1 (more or less!) but if we had included other columns in the block definition, these would have been factor levels, and there could be more contrasts here, for main effects and interactions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'constant_0': array(0.9999999999999998)}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cons_exper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model the low-frequency drift with the discrete cosine basis.  This is the default model used by SPM. See https://en.wikipedia.org/wiki/Discrete_cosine_transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 128 seconds wavelength cut-off high-pass filter.\n",
    "drift = dft_basis(128, volume_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10ded0090>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHQAAAD/CAYAAAAg7fI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADEJJREFUeJztnV+MH1UVxz9nt5V0uwSIllJaKQ9orQlJbUxfSmQVg6gJ\nNT4QjA/+CcZEicYHgzzRGmPURHgwITFCiCiKaEIoMdGCsDHVgBhbqVBoDRZKW7aQUgr0Zbt7fPhN\n4bfb3+zM3N/87p177/kkk52d39y9s7/vnHPuvXPvHFFVjHQYC30BRruYoIlhgiaGCZoYJmhimKCJ\nMTJBReQ6EXlORA6IyC2jqsdYiIyiHyoiY8AB4BrgKPAUcKOqPtd6ZcYCRmWhW4CDqvqiqs4C9wPb\nRlSX0ceoBF0LHO77/eXimDFiloWqWERszHFIVFUWHxuVoEeAy/p+X1ccW8DmzZsHFj569CiXXnpp\nowp9lakqNz8/P/D4sWPHWLNmzcDP9uzZM/D49u3b2b59+8DPRM7REhidoE8BV4jIeuAYcCPwhRHV\nVUnZl6yqpZ8thWs5H4xEUFWdE5GbgV304vTdqrp/FHUZCxlZDFXVPwEbXMqef/75A48vZRUrV65s\nbDWTk5ONzh+mnEuZqampxmVG0g+tVbGIlsXQMrrq5oalLIYuhYh4bRRlS+gJAyZoCa7CzM3NtXwl\nzQgqaJddaGhLc8WetiRGFi7XxdpcvUdor2OCtowJ2lFchQkdey2GJkZUFup697tYm6uFZt1tSZHQ\nLjcLQX22ckNbqMXQxIjKQl3dmfVDE8P1IbavutokKkF9flmusTC0oBZDEyMqC3XFxX3a47MOY42i\njuKzoRJro8hiaGJEZaGu5DRSFJWgNrBQTVSC+vySTVAHfD2ZyMlCrVGUGFG5XFdystCoBI2hH2oP\nuBtgrdxqLIYmRlQW6opPC83a5Tb95326wVjHcs1CS4jVQi2GJkZUFhrDROvQZCGozxkLoV1uFgt+\nfQoaGouhiRGVy42h2xKaqAR1JSeXG5WgMTSKQmMxNDGysFCLoR0lBpdb9tpTX5igLdc1NhY2ilkM\nTYyoHp/FMBxnLrcBJmg1JmjLWAw1WsViaMuYy21ADIKGdrlRCepKToJaDE2MqGYsmMutJiqX61NQ\n18aNCdoAs9BqhhJURA4BbwDzwKyqbhGRi4DfAeuBQ8ANqvrGkNdp1GRYC50HplT19b5j3wMeVdWf\nFKmaby2ODU0M/cnx8fGWr6QZwwoqnNtS3gZcXez/EpgmsKAu5VxdZ9QuF1DgERGZA36uqncBq1V1\nBkBVXxGRi0sLexopciHXRtFWVT0mIquAXSLyPD2R+ylV4bXXXntnf2JigomJiSEvJ12mp6eZnp6u\nPK+17IQichvwFnATvbg6IyKXAI+r6sYB5+uGDc2yUc7MzDhd28mTJxuXWb58uVNdq1atalzmyJFz\nkh9XUpad0Nk/iMiEiEwW+yuBa4F9wE7gy8VpXwIecq1jMarqtLkgIk7b+Ph4461NhnG5q4EHi+To\ny4D7VHWXiPwTeEBEvgq8CNxQ9ge6HEOzaxSp6v+ATQOOnwA+OcxFGe5ENVLkE9dWrj0PbUAMLjf2\ngYWh8PW0xYVYR4rseWhimMstwdXlWgxtgAlaTVSC+sRiqNEJspiX64K5XA/4XLOZ3dBfCMxCq7EY\nmhhRxVCfxGqhWbhcnzE0NFGN5bpignaUGAbnQ98Icd6GRilZWKi53MQwQTuKz5dBuQoa+tXmcd6G\nRilZDCyE7uz7JAuX6yJo6OearpiFtlgGwg9nWgxNjKhcris+uy2hLTQql+szhrq63NDdliws1MXa\nLIYanSAql+uKz35oaAuNyuXGMDh/5swZp3JtYRbaYhkI3yiyGJoYUblcV3y2ckNbqLnclslaUF/k\nJKjF0MTIwkJ9TieZm5vzVtcgLIaWEOsUFLPQlgktqMXQxDCX2zJZx1BfuAjq6jpDu9wsLNRiqBEt\n5nJbJusYai63fbKwUBdcb7bQgloMTQxzuSXEaqFRuVyfb5l2FSbrSWIpjhSFtlCLoYkRlct1xWJo\nR/EZQ2PIVToIi6ElxNooshiaGJWCisjdIjIjIk/3HbtIRHaJyPMi8mcRuaDvs1tF5KCI7BeRa0d1\n4aPGNc+ar9xsZdSx0HuATy06djaL7wbgMXpZfBGRD9PLdbYR+DRwpyzh75r+464J5sbGxhpvc3Nz\nTlvnBVXV3cDriw5vo5e9l+Ln54r964H7VfWMqh4CDgJb2rlUd1xuglhxjaEX92fxBc5m8V0LHO47\n70hxzPBEW61cJ79x+vTpdy9k2bLKJKwxDP2NiroZfl0FnRGR1X1ZfI8Xx48A7+87b11xbCArVqxw\nrH70dE3Qqakppqam3vl9x44dA8+rK6gU21nOZvH9MQuz+O4E7hORO+i52iuAf5T90aYNAp/vsO2a\noHWpFFREfgNMAe8VkZeA24AfAb9fnMVXVZ8VkQeAZ4FZ4BsauqedGa0lVW9csYheeOGFjcr0x9wm\nbNx4Tk73Sk6cOOFU1+HDh6tPWoSLBmVJ1bMY+nPB1eWG7vLY4HwJycbQXInVQm1wPjGiiqExTCUJ\nbaFRudwYBA2dfMAstIRYLdRiaGJE5XJd8TlJLGuX25QYXG7WgqYYQ0MLajE0MczlluAaQ0Pne4nK\n5cawcLdq1sWoicpCXXERJ9ZWrsXQxIjKQq2VW01UMdQaRdWYhZZgMdToBOZyS4jVQrNwuS7iWAz1\ngFlojfqD1m60ThYx1Of839AWmoXLzWnoLwtBfRJaUIuhiRFVDI2B0Baahcu1RlFHMUFr1B+0dqN1\nooqhsSZp9Ym53JYJ7XJN0JYJ3Ve2GJoYFkNbriu0hWbhcrteV5tkIajra2NixGJoYmQRQ10SzPl8\nDV2bROVyfRKry41KUNe736eFhsZiaGJEZaGu+OyHhh43zqJRlFO3JSoL9Wk1oVurrsR51UYpUVmo\nKy4u19VCQ8+TikpQn19yDI/qBhFVo8gVnzE0dCvXYmhiROVyfbpBi6EeiKEf6jLM2CZRCeqKzxga\n2kIthiZGFhbqM4aGdrmuGX5vE5GXReRfxXZd32e1M/w2TZzqkth1bGyM+fn5xptrXaETwtax0HuA\nnwH3Ljp+u6re3n9ARDbybobfdcCjIvKBtvKfxdDKDW2hlYKq6m4RWT/go0Hf7jaKDL/AIRE5m+H3\nyeEuczhy6rYM0yi6WUT2ishdfUnVLcNvYFwbRXcC31dVFZEfAD8Fbmr6RxbfzVUuNed+6Egz/Krq\nq32//gJ4uNhvlOHX10Nkny53dnbWqVwVI83wKyKXFMnUAT4P/KfYb5Th1xc2BaWPkgy/HxeRTcA8\ncAj4OliG3y4QNMNvUyu48sorneo6cOBA4zJr17q15U6dOtW4zPHjx6tPWkTWGX59Noo673K7RAwD\nC6EFtcH5xIjKQl3JyUJN0JYxQRvgczZ7rBZqMTQxorJQVyyGJobFUCNaC7UYmhhZWKitbTFM0C7j\n80sOLajF0MSIykJjWCZvFtoSb775ZifqanOidZ1JYYsxQUtYanb822+/XfpZm2QtqNHDBE2MoJPE\nglScEIMmiQUT1BgN5nITwwRNDBM0MTolqIhcJyLPicgBEbmlZplzVpjXKLNORB4TkWdEZJ+IfKtG\nmfNE5EkR2VOU+2GD+saKle47G5Q5JCL/Luqrvz7IZWRjFBu9m+u/wHpgObAX+FCNclcBm4CnG9R1\nCbCp2J8Enq9Z10Txcxx4Athas77vAL8Gdja4xheAi5p+j12y0C3AQVV9UVVngfvprQhfElXdDbze\npCJVfUVV9xb7bwH7qbEwWVVPF7vn0bsBK+sVkXXAZ4C7mlwjvdV+jfXpkqCLV3+/jIfV3yJyOT0L\nr3xtQOE69wCvANOq+myNKu4Avgs07R8q8IiIPCUiX6tbqEuCekdEJoE/AN8uLHVJVHVeVT9CbyHz\nx0Tk6oq//1lgpvAGC9bY1mCrqm6mZ93fFJGr6hTqkqBHgMv6fl9y9fewiMgyemL+SlUfalJWVU8B\nfwQ+WnHqVuB6EXkB+C29dbWL3yZTVsex4uerwIP0QlKtgp3Y6DU0zjaK3kOvUbSxZtnLgX0N67uX\n3qt56p7/PuCCYn8F8Ffgmgblr6ZmowiYACaL/ZXA34Bro2oUqeoccDOwC3iG3utx9leVK1aY/x34\noIi8JCJfqVFmK/BF4BNFt2DBy7NKWAM8XsTQJ+iJ85equhxZDezuq+thVd1Vp6CN5SZGZyzUaAcT\nNDFM0MQwQRPDBE0MEzQxTNDE+D8Ya4FWPkIZywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10de7fdd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(drift, aspect=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Stack the experimental and drift designs, keeping track of\n",
    "# the columns in the contrast matrix.\n",
    "X, cons = stack_designs((X_exper, cons_exper), (drift, {}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10e12c410>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIEAAAD/CAYAAADWiLXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD0hJREFUeJztnW2MXNV5x3//xcXCNkHILq/2uq0IMUhNaFS5Wtm1qZJG\nqI2C4AMQf2hSlKhSGiUKosJEQqilKk2keoWiFrWFukmLg3FUy0aVLBulfDDaEqJgJ61xcJUYjN8a\n08jGgPyyfvphrp3N7tyZe8/OnHPP3ucnjXbmzDlzzs595jnPeftfmRlOuxlJ3QAnPW4EjhuB40bg\n4Ebg4EbgMEQjkHSHpP2SXpf00LDqcWaPhjFPIGkEeB34GHAEeAW4z8z2D7wyZ9YMyxOsBA6Y2Rtm\ndg54FrhzSHU5s2RYRnAjcGjK67eKNKeBzEtVsSSfr06AmWl62rCM4DAwOuX10iLtl7j55pu7Fj5x\n4gRLliyZkb5o0aLSCo8cOcINN9xQq5EhZfqVu3DhQtf0o0ePcv3119eqp1+Z559/fkbahg0beOCB\nB7rmX7ZsWdf0YRnBK8BNkpYDR4H7gE9Pz3TXXXd1Lbx7925Wr149I33Xrl2DbWUPyi4mgJn1fD83\nhmIEZjYp6YvATjpxx9Nm9tr0fOvWreta/uTJk13fCzWCsgs21y5mKEOLCcxsB/ChXnlWrFjRNf3u\nu+8ufa+MK6+8slZ+6N29DLpcrDJjY2O1ywxlnqBSxZKdOXOmVpmQfxB6u/ac6RYT9GLZsmVRA8NK\nXH755SmrbwypN/YkNYKzZ8/Wyp/DLzrkgrbaCN5///2U1Q+FkAua2rh9FdFJ6wnqBoYxCXXRIb/q\nVncH58+fT1l9Y2i1EaTuC3sRemE8JnCyJKknaDKhv84cPUErjCBmnxtyQVttBO+9917K6htD6sDQ\nYwInrSfYunVryup7EnOeYHJyMqiuQZHUCDZt2lQr/8hImOOKOZ/vk0U12b+/3g70W2+9dUgtmUmb\nRgceEzh5LSWHEvOXlqMn8E0lJcSMCVptBFdccUXK6nviawdOq0jqCebPn5+y+p6E/jq9O6hb+bw4\n1fvaQW+SGoE0Y/fzUIg5WZRjTNAKIwihTTOGHhg6ee0nyGHs7t1BTRYsWJCy+p7kYHCDIqkRlB1N\nL2PHjh1DaslMYgaGHhM4yUnqCcr0CcoI9QS+n6A3SY2grgZBTNrUHWS1ipjDBE6OnsBjAsc3lZTh\n3UEkmqxPENMIUuNH00tokyfwmMBphz5B05eSU5OVPkEOizpuBA2l6aOD1PsqPCZw2uEJQvCYIBJN\n1ieIeTFTdwetOJre9NFBaiPwmMDJS58glKbP4qX2BFnpE5TdLqcfbgS9yWoVMeaFiRkThCqwDAqP\nCZy8dhaF4t1Bb1qhT9D0CZzU3UErjqbH9AQhv+rUnsBjAicvfYK5utun1d1BajfYCzeCSDRZpCKU\nkP8payOQdBA4CVwAzpnZSklXA5uB5cBB4B4zOznLdjpDZLae4AJwu5n9fEraeuAFM/u6pIeAh4u0\nWRNzy1fM0UHWngAQM0cYdwJri+ffBF6kxAhi6RP4PEFvZmsEBuySNAn8vZk9BVxrZscBzOyYpGvK\nCtfVJ9i2bVtYIyN6gpALmrsRrDKzo5J+Fdgp6cd0DGMqpd/m7t27Lz0fHR1ldHR0ls1xpjIxMcHE\nxETffAO7a7qkR4HTwOfoxAnHJV0H/IeZ3dIlv+3du7dWHffee29Q295+++3aZU6dOhVUV4g24+LF\ni4Pqeumll2rlH/hd0yUtAEbM7LSkhcAngD8HtgOfBb4GfAYo9eGx9Al8iNib2XQH1wJbJVnxOc+Y\n2U5J3week3Q/8AZwT9kHNHkVMZRWGYGZ/RS4rUv6/wEfn02jnLi0YmdR00cHqafPXZ9gwOS4lNwK\nfYKmxwSXXXbZEFpSHd9P4OSlT5DD+UCPCWqSWtN3GLRqiJgTMc8H5rh24DGBk5cnyEFCJnX/HoLr\nEwyYHLsD1ycoIWZMkNp7eEzgtEOfICY+bVyTuvoEoTuPmt4dpKYVq4gxSf2rDsGPppcQejFznDHM\nz3c5A6cV+gQxSb0sHEJW+gQ5aAt6TFC38oDt2bGIeTFTr6Z6TOD40fQyQtvmMUFN6n7RMd2mxwRz\njKbPGKaeBPOYwEnrCerqE8S6mSbE7Q5Se4KkRlBXn2Dz5s1DaslMYhpB6iFiUiNYt25drfw5GEGO\neEzgpPUEdfUJcpg2DhkdTE5OBtU1KFqxlBxCzJig1UbQ5E0lMWOC1KMDjwkc1ycow+cJIhFLnyCE\nNi0l+9H0EtpkBB4TOK5PUEbMHcCpv4eslpJzGCKGGI8bQUPxmMBpFa5PUEKbPEFW+gQ5xAQhF7TV\nRtDko+ltGh14TODkpU+wZMmSoHqavksotSfIaik5Jm3qDnxTSQmh3qPpt+TthscETl76BDkMEUNo\ndXdQV58gB0IMtdVGUFefIPTLiqkjFNJGjwmc5GR1ND0mbdpt7EZQQswhYqtjgrqk/sVUIYc2Tsdj\nAqe/J5D0NPBJ4LiZfbhIuxrYDCwHDgL3mNnJ4r2HgfuB88CXzWxn2WfX1Sc4efJkrfwX8dFBb6p0\nBxuBbwDfmpK2HnjBzL4u6SHgYWC9pFvp3Bv5FmAp8IKkD1rJf1lXn2Djxo218s+GNk0b9zUCM9st\nafm05DuBtcXzbwIv0jGMTwHPmtl54KCkA8BK4OVun11XnyDUCGLKz6c+XBpCaExwjZkdBzCzY8A1\nRfqNwKEp+Q4XaU6DGdToIMifbdmy5dLzNWvWsHbt2h655y7D6g4mJiaYmJjomy/UCI5LutbMjku6\nDvjfIv0wsGxKvqVFWlcee+yxWpWGflkxA8PU/ftUxsbGGBsbu/R6fHy8a76qRqDicZHtwGeBrwGf\nAbZNSX9G0jidbuAm4HtlH9rkTSWh5BgTVBkibgJuBxZLehN4FPhrYIuk+4E36IwIMLN9kp4D9gHn\ngC+UjQyc5lBldFAWwn+8JP/jwONVKq+rTxCzO4g5REw9fe76BAMmR8eXlT5BKDE9QciMYas9Qd0v\nLIdfWY7dgS8gOXktJYfS9E2jqT1BK4wgJm4ENal7ND3mEDEUjwmcLMlKnyCUpscEqW+T6/oEJcRc\nQErdHWR1NH3hwoVB9bgn6FN/0tqdRpCVPkEOniCkO0h9Q03XJxgwHhPUJNbR9Jgqox4TOFmSlT5B\n6Ha0pscEqT1BVvoEbgTDIatTyXN1P0FqI/CYwMnLEzS9HsjTE7RiP0HqcXg/Wm0EdY+mnzp1Kqie\nmLepy9ETeEzgpPUEdfUJnnzyyaB63BP0JqkR1NUnCDWCptNqI1ixYkWt/DH3GMZUFEttBB4TOO1Y\nSm76PEFqstpUEkrT72Te6k0lTT6aHpPU7fOYwHF9gjJy0EcaFK5PUELMIWLq7iArfYJQUn/J/Ujd\nPo8JnLyWkmP203N1JNKNrIwgJjlO+oSSlT5BTGLedCs1HhM47dAnCBmH+zxBJOrqE3hgOByy0ieo\ne1jlIiEXJlSouukTU93wmMDJayk51BPExLuDmjRZn8BjgkjU1ScI3YTS9NFBajwmcFyfoIzQiD2m\n1xkUWekTxCSHOYlB4aeSS8hBH2lQeEzgtMMThODdgRO1O2i1EdTVJzhx4kRQPTH7aY8JnCzJSp/g\niSeeGFJLZuLdwRQkPQ18EjhuZh8u0h4FPs8vbpT9VTPbUbz3MHA/cB74spntLPvsuvoEoUbQ9P0E\nqe+vXMUTbAS+AXxrWvoGM9swNUHSLXTum3wLnTumvyDpg2X3S66rTxBK0+cJGu8JzGy3pOVd3ur2\nzd4JPGtm54GDkg4AK4GXu332XFxFDDGC1J5gNoHhFyXtkfSUpKuKtBuBQ1PyHC7SnAYTGhj+HfAX\nZmaS/hL4G+BzdT/kkUceufR8zZo1rF27NrA5g2cuTBZNTEwwMTHRN1+QEZjZz6a8/Efg+eL5YWDZ\nlPeWFmldefDBB3/pdT+9gph9Z8zu4Ny5c0F19WNsbIyxsbFLr8fHx7vmq2oEYkoMIOk6MztWvLwb\n+K/i+XbgGUnjdLqBm4DvlX1orKPpTd9UknqyqMoQcRNwO7BY0pvAo8DvSboNuAAcBP4EwMz2SXoO\n2AecA75QNjJwmkOV0UG3wfzGHvkfBx6vUnksfYIQ3BNEosn6BG3aT+CriCW0yQh8AclxT1BGzHg2\ntSdohT5BzB1Mc3LtYJjEOpoeQptGBx4TOHnpE8SkTZ4gK32CJu9Ong2tjgnqHiure2xtNvg8gdMq\nstInSO02q9Ck/QRVyUqfoO79EWbDXNhUUpWsjqbHNIKYpDYCjwkc1ydoAq32BJK6Ps6ePds1vRe9\nhlkjIyNdH++++27pe6F1lfHOO+/ULhNST5WNpdNppBGcOXNmoEZQRq8LY2alj8nJydL3yjh9+nTt\n9rXCCJxm4J3ygCnrSiSVvlfmQfp5l0GhVEGJpObP/MxBzGxGv5rMCJzm4DGB40bguBE4NMwIJN0h\nab+k1yU9VLHM05KOS/phjXqWSvqupP+W9CNJX6pQZr6klyW9WpT7qxr1jUj6gaTtNcoclLS3qK/0\nPOe0MldJ2iLptaKNv1Opsl6TIjEfdAzyf4DlwK8Ae4AVFcqtBm4DflijruuA24rni4AfV6xrQfH3\nMuA/gVUV6/sK8K/A9hpt/Alwdc3v8J+BPy6ezwM+UKVckzzBSuCAmb1hZueAZ+kon/TEzHYDP69T\nkZkdM7M9xfPTwGtUENMws4t75OfTMdq+9UpaCvwB8FSdNtI5BV75+kj6APC7ZraxaOt5MztVpWyT\njGC6yslbRFA5kfRrdDxJV0mdaXlHJL0KHANeNLN9FaoYB/4MqDsWN2CXpFckfb5C/l8HTkjaWHQ9\n/yCp0oaNJhlBdCQtAr5DR2Wt7+S+mV0ws9+iI76xRlJPaRVJf0hH9W0P0zQeKrDKzD5Kx4v8qaTV\nffLPAz4K/G1R7j1gfZWKmmQEh4HRKa97qpzMFknz6BjAv5jZtjplCzf778Bv98m6CviUpJ8A36aj\n6zBdBa6sjqPF358BW+l0l714CzhkZt8vXn+HjlH0pUlG8Apwk6Tlki4H7qOjfFKFur8ygH8C9plZ\nJXFESUsuCnQVbvb36QSvpZjZV81s1Mx+g87/810z+6MKdS0ovBSSFgKf4BdqMGV1HQcOSbq5SPoY\nHbGQ/qQeFUyLbu+gE6kfANZXLLMJOAKcAd6kiI77lFkFTNK5iK8CPwDu6FPmN4t8rwJ7gQdr/m9r\nqTg6oNO/X2zbj2p8Fx+h82PaA/wbcFWVcr524DSqO3AS4UbguBE4bgQObgQObgQObgQO8P/9RpVo\nDX5BBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dc82550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(X, aspect=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we do an OLS fit of this model to the data, rearranged to a time by voxel array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_0 = OLSModel(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<nipy.algorithms.statistics.models.regression.RegressionResults at 0x10dfbc090>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = func_img.get_data()\n",
    "data_2d = data.reshape(-1, n_vols).T\n",
    "result = model_0.fit(data_2d)\n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7, 122880)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "betas_2d = result.theta\n",
    "betas_2d.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We rearrange the beta matrix from number of parameters (==P) by voxel, to (I, J, K, P)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vol_shape = data.shape[:3]\n",
    "betas_4d = betas_2d.T.reshape(vol_shape + (-1,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the center slice from the first parameter volume (p=0).  This is the parameter for the HRF regressor.  It shows some activation in the left frontal lobe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10e2aab10>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD+CAYAAAAalrhRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnVuMXWd1x//LGEIuxPg4eGziXJxMbBqnCUE0FEWQlASC\nQEr6FEGrKCHijYpUtIiENx6KSF8QUvtSFVAUQQullLiC1iY1BtEKCIoBi2BnZnJzLjO5HEgwhBDH\nXx/mzMl//z37P/vMzDlnyF4/yfI5sy/fbX9nr/Wt9a0VpRQkSdIu1o27AkmSjJ6c+EnSQnLiJ0kL\nyYmfJC0kJ36StJCc+EnSQlY08SPivRFxKCLuj4hPrFalkiQZLrFcO35ErANwP4CrADwO4B4AHyil\nHFq96iVJMgxW8sa/DMBUKeXhUsqLAP4VwHWrU60kSYbJ+hVceyaAI/T9Ucz/GFSIiHQNTJIxUUqJ\nxf6+konfmPPPPx/dbhdnnHEGNm3ahE2bNgEAjh8/3j/nd7/7XeWak046qfZ+L7zwQu2x1772tf3P\nL730UuXYq171KgDA1NQUzjvvvMoxPnf9+mq3rFv3smB07Nix2mNa3imnnNL//Pvf/75y7NixY5iZ\nmcH555/fr9di93n1q199wnV19eT+VBVOy2BKKZiamsIFF1xwQnm//e1v+59PPvnk2roofB9tO7dP\n67lu3bp+v+gzwNdp2dx27Reui17Hz5L20fr16/v9ovXk8tz4KVw3rUvEy3NUy4sITE9PY3Jy8oR7\nLlzX7XbR7Xb7f5+ZmamvR+2RpXkMwNn0fVvvbycwOTmJ6elp7NixYwXFJUni6HQ66HQ6/e/Dmvj3\nAJiMiHMAPAHgAwA+6C7gX0mF35x6rv6iurcXo7+8C7+MpZQTynP15F9i/ZXmt4nec6lf/nXr1p3w\ndtL76D3cYixfp21ouoj74osv1h7jftB76jHG9YMbB2W5C9F8nbuHa8NKzl1tFpMGBmXZE7+U8lJE\n/BWAvZhfJPx8KeUXdefzL9G4WSt1WSv1ANZWXTZu3DjuKvRZS/2ymnVZkY5fSvlvADubnPtK7cCV\nsLDWsRZYS3VZK+MDvHL7ZSSLews4MVJFXhYPF1t0qYOvU3GXxXRdxOIy3KKV1sWJjk7EdepD3XmA\nX2h8zWte0/+sfc33caqTLrLydbqo6sR0XtCrW2Rd7B7cJm2DUxUdXBcdL+4z7WsuT4/xdbp4yedq\n27l97tnRtnI/rYaakS67SdJCcuInSQvJiZ8kLWSkOr7qlqyrqN7inBmc7syoHs/3cfqi6s5cb7cW\nofVsajpy+pweYwcl1cf5u7bdOfewjqrt43NV52Y9V+Hy1BHH6bncdl1T4Lo5JyuFjzmzsfY1j7s+\nc7/5zW/6nwd5zvhcfea5DO3bpibJpuQbP0laSE78JGkhIxX1VcxikcWJ18584XzSnRjpfL1VJWnq\nk65iHYvQLKIDVbFZxTou35lAFSe2OvOhE8sZrSdfp2PL3wfdN1B3Hn937dO6uPKcCuT6k8fTqXg6\n7nwfHVvXPvcsLcfMmW/8JGkhOfGTpIXkxE+SFjJSHd/tZVfczi+n0zR153XlOZOd02W1faeddlr/\nM+9rB7w5j9G6sBlJdT1eb3j++edry1PcfnW3huF00qZ6pxsTZyZT/Zj7SceIz3X1dH2k48D31PKa\nrqfoPV2sBWdaXI4Lb77xk6SF5MRPkhYyUlFfcWKWC7PkxGsnfrJ4piYeDiulYqTbGcXeci4UltLU\nK+zUU0+tHFPxl+G+0DBZ3AZnRlKTHasMbiedXsf1dOYut3NPTYusLg2yY7NpSDLnParPmfOS5Ho7\ncV7hMdN+aaqSNCXf+EnSQnLiJ0kLyYmfJC1kpDq+0/UGcQdlHUeva+ri6txyBwkA6fRc3T3HcNu1\nPBcinHVGd3+FQ32rnsllON1S+5brOYguy6ge78pzOzab7j4cxH3Y7erjeqvOzW3X+7udl85s7NYi\n0mU3SZJG5MRPkhYyUlFfzR7OxMSilAZkcDuxGCcmO/HImey0Dc5sxaK/y7bi6qlwXzjTnt6Tr9M+\n4/YNYiZzplOXG8BlLeI+G2T8XH9yXVS14LHWZ5Dvo8+E2wXqdvVx+1QdYpVE1UYXfHY55Bs/SVpI\nTvwkaSE58ZOkhYzVZVd3rDGslzldVvUd1pNcplmn4w+SdIFxOqkLHOlwO+Jccg9tA5evei67BbsE\nEG4cHK6ebj2lqalN76P35L7WtnO/6Ji4wKYuMKbbcejWWtx6lXtelhN8c8k3fkR8PiLmIuJn9LeN\nEbE3Ig5HxJ6I2DBwyUmSjI0mov4XAVwjf7sVwN2llJ0A9gG4bbUrliTJ8FhS3iylfL+XCpu5DsAV\nvc93ANiP+R8Di4sx7kR2F0Pdea+puWS5ce75mAuQ4DwTnWebiqYsDmqwDZfjj81RKpY7M1ndPfQ6\n7Wu3w9CJ6dxevSeL4i4GvtI0aKaql+46pyI4laTpTjqXynyQXYSj9NzbXEqZA4BSyiyAzcu8T5Ik\nY2C1Fvfs6sL09HT/c6fTWVNpkJPklUK320W322107nIn/lxETJRS5iJiC4An3cmTk5PLLCZJkqbo\nS3VmZqb23KYTP3r/FtgN4CYAtwO4EcBdjW5iIpxojnHWW1SX5Z1mTXPQAz6BB+tzzhykOrBLAMGo\nzu30Y9XrGafPsVnO6bJOd3ZJT7QNPGbOJdmVp7os69Kq/7skJE2Ds7pgm7oOw+1zrrf67HI9dW3A\n7cp0Ofecm/NyaGLO+zKA/wOwIyIeiYgPAfgMgHdHxGEAV/W+J0nyB0KTVf2/qDl09SrXJUmSETHW\nQBwuXx2LlSq2sgjovJ2cOO/E6aaBGwEvnrkceC4uu1NfXCAHFnGdiKmiqTPZuUAjzlvOeUlyv6j6\nwIFG9Z58rgts6jwFFWcCbYrbDejUWy3P1WWQHIpNSF/9JGkhOfGTpIXkxE+SFjJWHd/t/OJzXUBG\nl0NddXyXIMG5ATtTEetemsSiadQUl8DD6ae6NuBy4DlzEPenJvDgMnRtwOXVY1wbnEuyjkPTHINu\njURxZjm+pzNlan/yM+L0fz3m1ovq7rHUubX3GPiKJEn+4MmJnyQtZKSivtvBNYgJhkUpF+BCxXLn\naebUAD7XxTvX9rlcfU1NRy5whPLMM88sq55cF2c2cjvNnFrl1BwVr5uaMt3OS312uE1OtXC77FyQ\nFZdT0AUh1b5uulPQ5RtoSr7xk6SF5MRPkhaSEz9JWshIdXzVoVhXUZ2GdSE1l7jkCXX3ALwZ0Jm7\nnKnIJc1wCSfcrjfW49Wk5dxyuTw1yznTossbyNepbukCm7ocfy7yEpfnTGhut5w75iI9ORdotwbl\n9H+nq+s6hbpE1zFIApY68o2fJC0kJ36StJCRivrOw+i0006rfGcRzO0mc2Kdmr5YNHaBP5x5xIl8\nvLMMqAYMUbH1+eef7392pikV617/+tf3Pz/77LOVY08++XIgpC1btlSOcRl8D6DadpenUEVKVTWY\nX//61/3PLlffAw88UDm2Y8eO/mfuI8Cn+mac2qjj7oKsNhWh3c5EF5xVj7kdm1w3l1K+KfnGT5IW\nkhM/SVpITvwkaSEj1fHVXdJFk3GBMV3CCZePzwWVrCtbr1N9jnVwrcvRo0f7n3/1q19Vjk1MTPQ/\nb9hQzUDG+rHq+FwXXTfYuHFj/7Pq/6xbqp7rXJmdfsxlbN26tXKM++L000+vHOO1EDU78g5HXW9g\nHXiQoJnOLMf9q+tMXB6PJeBdrp3rLZfvAobqM8/fX/e611WO6VpIE/KNnyQtJCd+krSQse7OY5HP\n7RhjM47ivJ9c4A+3a0rFQRY51XuNxa6pqanKMTZ3sRiu93HeeRrcg3fgqZjM4qiKf+eee25tedxn\nqiKw+Ot2Hx45cqRyjMdM1SPulze84Q2VY6xObNq0qXLsl7/85aJlA9UxU3XPBQxxuQ9doBj+7uL4\n6/PJ4+ICvrrnzO12bEq+8ZOkheTET5IWkhM/SVrIWF12Wbd0UUxcYgyXk96ZS1RfdZFtuDw1pTz8\n8MP9z2yGA4Bt27Yteo/Fzq2rp7aP26R6oFtTcLnsWK/WNQVnCnO57FkH1uzIvMaweXM1wzrr8aoD\nsxuyPhNsTtS1FtardRxcRBw+pjo+t12fHbdrsalJUtvO9XYm66Y0yZ23LSL2RcTPI+JgRHy09/eN\nEbE3Ig5HxJ6I2LDUvZIkWRs0EfWPAfhYKWUXgLcD+EhEvAnArQDuLqXsBLAPwG3Dq2aSJKtJk6SZ\nswBme5+PRsQvAGwDcB2AK3qn3QFgP+Z/DGpRcZDFF/WMarozSkUpFmld+mJXFxd0kUV7oCpWnnHG\nGZVjTz/9dP+z7ohjEVpNbyxuO5OPinws/rpgou4+rt/VrMpjpsdYRVERmtUcVXnYRKn9wmqAjju3\nnVUsoDpm+kzU3QNonpdR+8yZlFlFcIFNVU1ldWI5wTWVgRb3IuJcAG8G8AMAE6WUOaD/47C5/sok\nSdYSjRf3IuI0AF8DcEvvza/xr2vjYU9PTy/cA51O5wTHjCRJVk6320W32210bqOJHxHrMT/p7yyl\n3NX781xETJRS5iJiC4An666fnJwE4EWnJElWRqfTqVhQZmZmas9t+sb/AoD7Simfo7/tBnATgNsB\n3AjgrkWuq6CuoqzPqZmM9U79wXA/IKzD6XmsL7tdU2qe4bpolB2WXtSExtdp2/ncp556qnKMdUbd\nMebcjrluuiOO9WVtO+vqKo25nIIuqYQLlsp1m52drRxzwVOfe+65/uczzzyzcozXCnTtg9unayaP\nP/54/7M+gy7CED9nztXXRdJxOr6WzW3QZ2I5LDnxI+JyAH8J4GBEHMC8SP9JzE/4r0bEzQAeBnD9\nimuTJMlIaLKq/78A6nYFXL261UmSZBSM1HNPcTHU2aSlIhF/dwER1KSlomrduS72uprlXCAH9jRT\nM5ILfsmiowsAoWKkM4Xxd62nyxvA91Q1h8VRlxuQzZpaho4t101NhBywRHcRulTmF110Uf+zqh0H\nDx7sf77wwgsrx1gN0Hs671Gut6odfK6qY6yuqBrHfa1j60y+deRqW5K0kJz4SdJCcuInSQtZMzq+\ni4zicq+7Y85cotFrXB6/7du39z/Pzc3V3lP1QL6PmmdcEEu3S4t1RtWBnY7v9EDWnV2UHZfbzUWo\nUbdcbrvu3OM+dMd0jNjUp33NfabPGbtZq47PAVLdc6Z95nZCurUBHj8XbHPkLrtJkrwyyImfJC1k\npKK+imDsIacBIJrmRnN5y1ScZ5FP789inYqt7M3mUlMrTdN5u9xnWk+u2yAx253plMtw46CmMFZR\nVK3icXF1UfhczgUIVNunplkXnIIDgeo9zz777P5n3dXHuQ80x5/LZcfj4lK8q4rnPAV5zHSM3I7D\nOvKNnyQtJCd+krSQnPhJ0kJGquO76C5qfnJ52ZuaAVXvZD1Jd+DxFkbNLc9mJZf/T/VObpOa3tzu\nNUb1N26DutDyMW07r3e46C46DnxPl7xEdWznluv6jMdP12+4ntovPC5aHve9JvDgsd69e3fl2IED\nB/qfL7300soxt77hTNHcn/osOZdyfkbc+ldT8o2fJC0kJ36StJCx5s5jUc55I7lcdipC666tuvJU\nRWDTzTnnnFM5xqKpy/+nZjIO7OACf+oxFqFdcAg9xnVTkZ2PqUjJIqce4/K0z3jnoIr6PC4uSIeK\nu9wXWhfezaamLzfu3J+qxtXVS1EvwkceeaT/WUVvl3bcpWp35jwuQ5/B5US2yjd+krSQnPhJ0kJy\n4idJCxnr7jzWf1R/ZLdHZ7ZyASDVTMY6PidnAKqum2ruYv1KzYBu3YD1bNWBXT2bJhPRe7Kup+7K\njJqfXK43p5OyHj9IHjiOsuPWDZzu6iLUuJz0jz32WOUY94XW8+qrX44sp+spHCxVQ1qzDq51ce7Z\nzCA5IZdj3ss3fpK0kJz4SdJCRirqq2jq8oix2KXmC1YRVCTi61Qs5/uo+dB5YrH4q8Ev+ToV43gn\nmAZWXEgyApwoKnLgCjX1OTGSRUA1bzmx3KWRZhFXr+NztX18TNUxHjPeFQlUnxEOEKI4FU9x6dG5\nr5955pnKMVZJ2NwLVPtF1QdWs5xKovOBnzN9Pl3gj+WQb/wkaSE58ZOkheTET5IWMlIdX01FrEur\n6y3rRqrvOB2H7+mSPLiEglqX+++/v/+ZXTUBYNeuXf3PuvPLma2cqy+fq3o8X6fmJ9ZzNUkH68Sq\nH3P5mjuvaeQg55Ks8DqJmmp5rDmAppan/cn1VNdXdV9mNm9+Obu71pnL0GeJ1wa0Dbze4Vy89Z5N\nzdu6rrWc4JtLvvEj4qSI+GFEHIiIn0fEp3t/3xgReyPicETsiYj6lZgkSdYUS078UsoLAP6slHIp\ngIsBvKuXSPNWAHeXUnYC2AfgtqHWNEmSVaORqF9KWbA7nIT5H4tfArgOwBW9v98BYD/mfwxqcYEp\nXWAFFW3YTKdindvxxzHU3/nOd1aO3XvvvYuWDVTFQfbwA6pqgYuBr15oLMaqiZDvo+oDi/5OvFZ1\nxZkrWUVQsbWpmUzbx8dUtdBgkQy33e1W03tynjv1WuRxV683Fzuf+1DNjiyKO1VCcSY7Zxp2uySX\nQ6PFvYhY10uRPQtgfynlPgATpZQ5ACilzALY7O6RJMnaoekb/ziASyPidAB7IuJKAEVPq7t+enq6\n/7nT6ZywtzlJkpXT7XZPcAarY6BV/VLKcxHxLQBvBTAXEROllLmI2ALgybrr2EstSZLhoC9VZ7la\ncuJHxBkAXiylPBsRJwN4N4BPAdgN4CYAtwO4EcBdS93L7VRSXZL1Yz3m3E8dU1NT/c9ve9vbKsce\nffTR/uetW7dWjnH56mbJup/mnedjqpPyuforzYOn17G5TfVxl8CD+0nXTFi3VL3TuVWz26qOEZu7\nnDlWd725pCCMrotwf2r7uLzzzjuvcuzgwYP9zxxcE6hGYnrqqacqx7gMTcTBbdf1Bl4L0fnAbdf1\nMLcO4/qpjiZv/K0A7oj50VsH4M5Syv/0dP6vRsTNAB4GcP3ApSdJMhaWnPillIMA3rLI37sArj7x\niiRJ1jqhYtGqFxBRrrnmGgAnimBNxXQXsNCZm/T+btcbqw9qCnvooYf6n9UjjkW3QYJtsGiqnnts\nftK03C6ngDMVuZxt3L9qQnPmQxern01cbrejegryGKmZjMvQezq1iseMvTCBqqlWzbGshuzYsaNy\njBestTyup96z6fjpOLA6qH1dZybfs2cPSimLuvWlr36StJCc+EnSQnLiJ0kLGenuPNVFnL7D+pUL\nDqm6LOtJatJivVrrwrq6up9u3769/5l1UKC6HqAmJj5Xk3RwG974xjdWjrHeq66iHCRUzUG8jqBm\nJNZDtQ2M6tW800zXPtis5HahqU7KY+Si0GhduN56HX/XseU+1LrccMMN/c/f/OY3UYdzKXfmNB0j\n95zxGphep2ZPV0YT8o2fJC0kJ36StJCx5s5j1NzFO/BUnOdz1RTm1AC+p5rluG4qKrLIqaI3i2dq\numHRWNUVrpuabrgMNYHyPbV9LoAHf1eR3eUidAE1uL2qqrFapSoQt0lFWBZ/tS4sUmv7XCrzp59+\nuv/5rLPOqhz73ve+1/+swVlZxWOvT6Dapp07d1aOsZrjgnu4vAFaF5e/0e16rSPf+EnSQnLiJ0kL\nyYmfJC1kpDq+6i2sn7udSqrHs048SK5w1gNVt2RTkeqrrDOq/shwAg2gasJzrpua/IL1Qufu6hJc\nqBup0/HZfKhmK+4XTZrBOreWx+haBNdFcxhyGzg/HeAj8riIOC43II+D7q7kummQVXb11XUfboPT\nv934ued4OTq9km/8JGkhOfGTpIWMVNRXkY9FFufFpLD4q6I+i9QqfrIIqLutuDz1CuP7aNgwvm7L\nli219VQRmvtCgzxwXzgzp4qt/F3VKjZfOrOqeqFxf+r4uZx7PLa821Drpl6E2qa66/Q8HiMVoVnl\n0rZze7XPeMxYtNdjCvfTIDshuXyXOt2ZWJuSb/wkaSE58ZOkheTET5IWMlYdn78711vVk5zbI5+r\nehJfpwEuWQc+88wzK8dY9+JoPFqeJttg1023S2sQ8wz3k9NznVunywOnY8T6pPYnt0FNaG73mstz\nx2Ok+erZBOvapzsT2dzm9GNdE2Jz4sTEhD2X0bUCxu3A4za4NS5dw3CmvzryjZ8kLSQnfpK0kLGm\nyWZxxpk9FD6mIh9/d+mEVWRnrzQNAOmCSrCpShOH3HPPPf3P6qHmcvy52OvOrMMipnoKchvUXOl2\nxHF5KmLyd+0XF6TTBWBxKolLB825CLXPWPTn/AlahstFoOrfE0880f+sYr8T07k/nbqifcbfVT1a\nTsDcfOMnSQvJiZ8kLSQnfpK0kLFG4GF9R00grGuqDsM6nOqWLgoNoyYfNmlpPbl8DYzJ+eu1nm6H\noXMV5TbodaznDqJb8rnOhKZt51z2uqbA/avHuDy+h6JrEa5fOOqNXsfjose4X9Tlmtd6zj333Mox\n7vvZ2dnauuiuPhc9ygXbdIlinPl3qME2I2JdRNwbEbt73zdGxN6IOBwReyJiw8ClJ0kyFgYR9W8B\ncB99vxXA3aWUnQD2AbhtNSuWJMnwaCTqR8Q2AO8D8HcAPtb783UAruh9vgPAfsz/GNTigiCo+Mkm\nJxWlWNRREZNFIid6b9hQFVC4LiqasqlIYQ/AI0eOVI6xCU3NVmyecaKaesRxH6qoyCKm9gu3z12n\nAUqchxqL5RpQ0wXG4LF2gT80ICpfp2oAi/cqzrPXnT5L7J2n43748OHa8lhVdOPngp7qM+/SnHN5\n+jw6VaqOpm/8zwL4OACeSROllDkAKKXMAti82IVJkqw9lnzjR8T7AcyVUn4SEVeaU2tXljizaKfT\nOcFBJkmSldPtdk/Yg1JHE1H/cgDXRsT7AJwM4HURcSeA2YiYKKXMRcQWAE/W3WDBo815XiVJsjI6\nnU5FzZmZmak9d8mJX0r5JIBPAkBEXAHgb0opN0TE3wO4CcDtAG4EcNdS93LmPDVROBMT/4CovsrH\nVPdh/VgDY7LLrtPH9ToXIJHXGAbJk6Zuswz3hV7nyuP+VX2c1zt0baBpsEgXmUjdld1OOgf3i+q5\n3Be6lsTutRowlNcwdD2D15n0meDy3G48Ncvxuog+1/ycab+7Y8thJQ48nwHw7og4DOCq3vckSf4A\nGMiBp5TyXQDf7X3uArh6GJVKkmS4jNRzT8VPFnVUzGIRSdcGnGnKBUGsuz9QFeV0BxeLdS5wpHqM\n8Xc1k3F5bPZTVF1xAUpYbNZ7soipdXHipwsOwdexByNQNW2qKMxedqpacI6BXbt2VY6x6U1VLufV\nx9+///3vV47x2LoALNpnusORYVXD5W/UZ9d59S03cEsd6aufJC0kJ36StJCc+EnSQsaq4zfdVaT6\nlNu551x2WX9VPYnXGNQc9OCDD/Y/q87Gup+Wxzu/tm7dWjnGawWqk3LdVFdnM5aapvg6t6ag6wZu\nRyO722ouOzYDqlmO9X912OJz1dTHuruWx1xyySWV72yz1ueMTY26LsLjoME9XS5Cbp+6JLs1KH5+\ndH3DRWXKhBpJkqyYnPhJ0kJGKuqr95PDmTZY1HFebm4Hl4rlXDcVky+88ML+Z/UYc3XZvn17/7Pu\nNGPzk6oWLicdt0lz7rFJzaUPVzGSRXbdMcbt1WPO645Fag2WwuK9HuM+U7Gcxe1rr722cox34H3n\nO9+pHNu/f3//s5odWbXR9vAz4oKcqKrmduDxeKoK60x9ThXNYJtJkjQiJ36StJCc+EnSQkaq4yus\nG6lJhHXupjvE9FzVvVzOPRfkkXViNTGxqU/L41x6vEMMqJrUdG2AdW52YdU2qI541lln9T+razHr\n1WoOUjMWw/2p53E0G9WdeZ1E10VcTkG+zvXn17/+9cox7rODBw9WjnFiEzVlss7tzMvObKzPJz8/\nzhVd0fWOpvXM3HlJkjQiJ36StJCxeu6x+KQmCheIw+U7c7nXnKjv0kE7j6qLL764//nxxx+vHGOR\nlsVwoNr2xx57rHKMA05qG1gU1sCRLAKq+sDitlNlnBip/XLo0KH+Zx0H9pbjNNVA1fSmoaJY5VNT\nGPfv3r17K8f4XJfDQIOQcJsGGXd+XlR14vL0OvcMchtUjXM595ZDvvGTpIXkxE+SFpITP0layEh1\nfNWF+LszpahZh69THYrv46KYqBmQ9V41+fB1ql+xzqjXsW6rbWfdVnfgsf6vbqQXXHBB/7Oa11iv\n135xSSxc4E8+V/uMdW51c2bTm9bFrWEw2i9TU1P9zy5HvPY168daT5f/j8da10y4TYOYovlZ1meJ\nn0EdBxdhaBBX+AXyjZ8kLSQnfpK0kJGK+oPsKmJRygV5dLHlFRbltC4sPqko7EQp523FIq2Kg+xJ\np+Ig11PFQRZ/nRipZjLeyaftcZ5fXE8NqLF588tZ07Tf+TqXf28QnNcbPwfaPjaNqYrgcjS44Bfu\nOWu6q0/7hZ9JFeedyrUc8o2fJC0kJ36StJCc+EnSQtZsBB6nzzGqY7NOrGYkl1vORahx5kOX4MJF\nW+HdbIPstuLdeqoHssnH5cDT9Q0+V3VgDhLqdtLpWgSbGnWMuF+0z7iftJ5chkuIono130fdlZu6\neDvXcK0nf9e6sF6v5bn1AD5Xn4/l6PyNZldEPATgWQDHAbxYSrksIjYC+AqAcwA8BOD6UsqztTdJ\nkmTN0FTUPw7gylLKpaWUy3p/uxXA3aWUnQD2AbhtGBVMkmT1aSrqB078kbgOwBW9z3cA2I/5H4Na\nVFxSEYlhEdDtpHPimdvxp+KRM5M50ZvLcCKYxsB3sez5mIrzHHBC4fuoeM1ehSoms4rgchEuN9eb\nivPOS9OlK3c5E5x3JY+nE4udGqdw/2pdnCmanxetp9sh6nZJNs0XWblfw/MKgG9HxD0R8eHe3yZK\nKXMAUEqZBbC59uokSdYUTd/4l5dSnoiINwDYGxGHMf9jwNSuTExPTwOY/7XduHFjZa92kiSrQ7fb\nPcFxq44RHY03AAAKLElEQVRGE7+U8kTv/6ci4hsALgMwFxETpZS5iNgC4Mm66ycnJwGsTgCBJEkW\np9PpVF6qnFZMWXLiR8QpANaVUo5GxKkA3gPgUwB2A7gJwO0AbgRw11L30onvdqGxTqM6YdOdZqpf\nsa6n7pKsAztXSq1LU3OQXueitLDJaRD3WheZiNvkdGC33qBrJlxvZ8rUfuHytTy+p9bTuUfzWKtJ\n0u2W43rr88LrKzoObq3AmRbdWoS7P/eTS5bSlCZv/AkA/xERpXf+l0opeyPixwC+GhE3A3gYwPUD\nl54kyVhYcuKXUh4E8OZF/t4FcPUwKpUkyXAZa7BNFt3U04xFMBeswYl1ek+3+8mJ8yx+uh2Gg6Q9\ndh5cTc1PKkY6z0Q+5syqziyn5fGOQw1UweOnbXf9yWqOyymoYn/TPnNehO4Z1DYwrn36LLkdok4N\ncEFkMq5+kiSNyImfJC0kJ36StJCR6vjOHOPMJapDOVdRZ9Jyu5/4XLem4Ha9Ob1aj3F5gyQa4WMu\nWYO6dfKuN9XH+To1y7GZU/uP7+OOuR1xqgO7xBFuzYTLd/kUtZ7Ojdvhctlx3XQcnLnZmaldf7rn\nuo584ydJC8mJnyQtZKxpsllEUtHNid7OtNFUdNNjXJ7ztnLx6jWgpgsm0vSYC9LhPL/U+9CpL9xe\nVcdUrKyri6orbgeeU+NYpNU2NA2M6UxozivT7QJ1Oz1d4BbtT7dzT9VIhp8t93w2Jd/4SdJCcuIn\nSQvJiZ8kLWSswTadvsoutarTuICajIvw4+rmTG8K64xaT5fHr65soKqPq2sx69wu4o+2neuipiLu\nQ10LcIExnfmwaX+qy7VLHOGeF+57rQv3hQY9Zb3a5aRT/Zvr7cbPud42TSij5bu8gU3JN36StJCc\n+EnSQtZMXH0VMZ3ZwwVWcJ5RTQMyDOJJ19QM6Ly7XAAPNac1jWKk4q7bMcYmJxe40QVEUdG06W45\nNa85laipaKxmVR4zFZNd0BMub5Dnk49pee554TFzbdB6OnW3jnzjJ0kLyYmfJC0kJ36StJCxRuBh\nfdK5sKpOw+c6/dGZu9wuO9WjXQ411ntVn3M7zZzLpzMD8j3d7jzV1fk+bpekHuP7aPuauoo692h3\nrjPHunz1rjxnrnSBVJ3Z0Zky3U66QXZeugQezj26jnzjJ0kLyYmfJC1krOY8FnUGiUnvRGFn2mCR\nzO1e03qyR5fbLedioTvvLj3GYqVTEZy50onzijMxudTUbveaC5rh1A4el0FUBD6m9+Q2ud2GTmTX\ncednQr0B3W5Ovo8+89x27U/O7ej6uin5xk+SFpITP0laSE78JGkhjXT8iNgA4J8BXATgOICbAdwP\n4CsAzgHwEIDrSynPLrciTn9zbpbONVX1eD5Xd2I1dfV1awpuN5kzH7q1AReZyNXFtcHdU+vCwTbd\nGDnTojOTOVOtc2l1x9xahDPxOpdd5ejRo/3Prj/12WX9XE2uLlef2606zN15nwPwrVLKHwG4BMAh\nALcCuLuUshPAPgC3DVx6kiRjYcmJHxGnA3hHKeWLAFBKOdZ7s18H4I7eaXcA+POh1TJJklWliai/\nHcDTEfFFzL/tfwzgrwFMlFLmAKCUMhsRm5csTEQwZ+5yQR5YfHJx2V3MfbejSuvivNeaxoFX1cJ5\nW/ExNU01jcfvUoS73HkuTbYzSSrOxOQCqXCbBjGvMW4npFOrXF49Z8Z1O0R1jNj05wJ4OjXA5RRo\nShNRfz2AtwD4x1LKWwD8BvNivpY2eOlJkoyFJm/8RwEcKaX8uPf93zE/8eciYqKUMhcRWwA8WXeD\n6elpAPO/VJ1OB5s2bVphtZMkUbrdLrrdbqNzl5z4vYl9JCJ2lFLuB3AVgJ/3/t0E4HYANwK4q+4e\nk5OTAJYX/ztJkmZ0Oh10Op3+95mZmdpzm7rsfhTAlyLi1QAeAPAhAK8C8NWIuBnAwwCuX+omqluy\nbqT6HOs4Tr9yepIz66iexPdxuqUzk7mde87N0u0wdKYwxyBmQKdzD2LCY7hNznTq7unWaNRUy88L\nmyCVQdZ26spWdGy5LmrGdbvsBgnIulIaTfxSyk8B/Mkih65e1dokSTISUvZOkhYy0t15zgzhPLEU\nJ3o7nOeX885zIpgTk51Y50w+zgPPiZxNY7a7YI2DlOfE5KbHXL/oTjq+bpCU64zLw+hw/ekClA7i\nZedyGK42+cZPkhaSEz9JWsjIJn5T++IoWCt1WSv1ALIudbxS67JmJn4ppfIvIvr/BoHvoSzcr9vt\nVu6/nAgmi6H3dPcvpaDb7S65RuHu6crQ/lyKhfEZ5LphjBHXZdDrBi17WP0yLP4gJ36SJGuHnPhJ\n0kJi2GJLROTmnSQZE6WURfWwoU/8JEnWHinqJ0kLyYmfJC0kJ36StJChT/yIeG9EHIqI+yPiE8Mu\nT8r+fETMRcTP6G8bI2JvRByOiD29CMKjqMu2iNgXET+PiIMR8dFx1SciToqIH0bEgV59Pj2uuvTK\nXRcR90bE7nHWo1f2QxHx017f/Ghc9YmIDRHxbxHxi94YvW016zHUiR8R6wD8A4BrAOwC8MGIeNMw\nyxS+2CubGVd04GMAPlZK2QXg7QA+0uuLkdenlPICgD8rpVwK4GIA74qIy8dRlx63ALiPvo8zgvNx\nAFeWUi4tpVw2xvoMN7K1eiSt5j8Afwrgv+j7rQA+McwyF6nDOQB+Rt8PYT5QKABsAXBolPWhenwD\n8/EMxlofAKcA+BGAC8dRFwDbAHwbwJUAdo97jAA8CGCT/G2k9QFwOoCZRf6+avUYtqh/JoAj9P3R\n3t/GyeZC0YEBLBkdeLWJiHMBvBnADyDRikdVn554fQDALID9pZT7xlSXzwL4OKrBWsfSJz0KgG9H\nxD0R8eEx1acf2bqnAv1TRJyymvXIxb0RRweOiNMAfA3ALaWUo4uUP5L6lFKOl3lRfxuAd0TElaOu\nS0S8H8BcKeUnAJzD/yjH6PIyH036fZhXx96xSPnDrs/QI1sPe+I/BuBs+r6t97dxMhcREwAQS0QH\nXm0iYj3mJ/2dpZSF4KRjqw8AlFKeA/AtAG8dQ10uB3BtRDwA4F8wv9ZwJ4DZcfVJKeWJ3v9PYV4d\nuwyj75fFIlu/ZTXrMeyJfw+AyYg4JyJeA+ADAHYPuUwlUH2b7MZ8dGBgiejAQ+ALAO4rpXxunPWJ\niDMWVoQj4mQA7wZwYNR1KaV8spRydinlPMw/G/tKKTcA+M9R1mOBiDilJ5EhIk4F8B4ABzH6fpkD\ncCQidvT+tBDZevXqMYLFkvcCOAxgCsCto1qk6ZX9ZQCPA3gBwCOYjw68EcDdvTrtBfD6EdXlcgAv\nAfgJ5ifZvb2+6Yy6PgD+uFf+AQA/BfC3vb+PvC5Upyvw8uLeWOqBed16YXwOLjyvYxqjSzD/4vwJ\ngK8D2LCa9Uhf/SRpIbm4lyQtJCd+krSQnPhJ0kJy4idJC8mJnyQtJCd+krSQnPhJ0kL+H0HC/Xia\nM1xvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dfbc750>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(betas_4d[:, :, 15, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we make an estimate for the AR(1) parameter from the residuals of the OLS model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(169, 122880)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.resid.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Estimate AR(1) coefficient\n",
    "lag_1 = (result.resid[1:] * result.resid[:-1]).sum(axis=0)\n",
    "resid_ss = (result.resid**2).sum(axis=0)\n",
    "nz_ss = resid_ss != 0\n",
    "ar1_estimate = lag_1\n",
    "ar1_estimate[~nz_ss] = 0\n",
    "ar1_estimate[nz_ss] /= resid_ss[nz_ss]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(122880,)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar1_estimate.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We round AR(1) to nearest one-hundredth and group voxels by their rounded\n",
    "# AR(1) value, fitting an AR(1) model to each batch of voxels.\n",
    "ar1_estimate = np.round(ar1_estimate * 100) / 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "white_betas_2d = np.zeros_like(betas_2d)\n",
    "# Loop across bins estimating model for each\n",
    "for ar_val in np.unique(ar1_estimate):\n",
    "    ar_mask = ar1_estimate == ar_val\n",
    "    m = ARModel(X, ar_val)\n",
    "    d = data_2d[:,ar_mask]\n",
    "    results = m.fit(d)\n",
    "    white_betas_2d[:, ar_mask] = results.theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the estimates for the first parameter (HRF parameter) after AR(1) modeling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "white_betas_4d = white_betas_2d.T.reshape(vol_shape + (-1,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10e35ae50>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD+CAYAAAAalrhRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnW2sXeV15//rQkJ4CY6P8QvBMWAupsQJBNRhpkIR0BAa\npUqZTyidSQVF/dZRMspMFci3fpiqzJcq0syX0aQRQu1MM51JcaVkcBwLpGnTAooDJsTG5oIxL74G\njsNLXgjgZz7cc0/++++7//c5x/eec8NeP8nyPmef/bzt/dy91rPWs1aUUpAkSbeYmXYDkiSZPDnx\nk6SD5MRPkg6SEz9JOkhO/CTpIDnxk6SDnNbEj4jPRMSBiHgqIr6yUo1KkmR1iXHt+BExA+ApAJ8C\n8CKARwB8vpRyYOWalyTJanA6b/zrABwqpRwppbwN4H8CuHVlmpUkyWpy5mlcexGAo/T5eSz8MWgQ\nEekamCRTopQSS31/OhO/mu3bt6Pf7+OCCy5Ar9fDhg0bTvnNL37xi8bn97///cPjiGbb33rrrda6\nPvCBDwyP33nnnca5M844AwBw6NAhbN++vXGOf/u+972vcW5mZmbJ3+m5d999t3Hu3HPPbW3zO++8\ng7m5OWzfvn3YrkVOnjzZ2pa33357eHzmmc3bx9epCqe/1esOHz6M2dnZU+r7+c9/PjzmsV3sQxt8\n/7TvPE7cZmDhHj399NO47LLLcNZZZ7Vep3Vzf3U8uU963S9/+cvW684444zhuOh48md+BrSdCt8H\nbQs/51rfzMxMa1sWr+v3++j3+8Pv5+bm2tvRemZ5XgCwjT5vHXx3CrOzszh8+DAuv/zy06guSRJH\nr9dDr9cbfl6tif8IgNmIuBjASwA+D+D3RynALSzyOX3j61/YNvRtslhOKeWUMvit4P66619zLmfU\nv/wzMzNLvom5v+7NpvAbS+vmsdDxXPwcEQ2JwrVL26J953N6H9rafDpwHU6CWo1NaTouq8Fiu0sp\nrW/8URh74pdS3o2IfwdgNxYWCb9eSvlx2+/5L9G0WSttWSvtANZWW9avXz/tJgx5r47Laen4pZT/\nC+CKmt+upQFcK21ZK+0AsOS6y7RYS+PyXm3LRBb3FuGFFIUXg4CmqKpipFuocgtH/FnrY/HQLVo5\nMVJx5bB45lQCxS00srqiZXJ92gcuUxfi+JzeP7cYxb916pG2ha9TtWNc9Y/L0XbymGmZTj3i50fb\nyW3Rc27MGHeuVtV1pMtuknSQnPhJ0kFy4idJB5moju90S4X1eNV3nE7MOpRbC9AyWV9W3Zn1Oafj\nO33fmWC0Pm63lnn22WcPj9m5Bmjqx865R2E9VB14uN2qr9bqx+qIw/3VceHx1fUGrk/7554lZ3J1\npj5nHv3Zz37W2pZa/Vyv43FxZa6ESTLf+EnSQXLiJ0kHmaqo78RkFomcZ9JSPs1tcDnOtOg88Jw5\nSK/jOlTcZbFZfeCdOciJeW48nWrhvN4YpwYo3Ad3bx36O7cvwj0vtSqfUyG1r3zPtH/8Wx0zrkP3\npzj11j1LTo1rI9/4SdJBcuInSQfJiZ8kHWSiOr7bf6w490y388vpXrXmNtUJWS9THdjpzrwf/6c/\n/WnjnNNXGa2PzUiq63Fb1NTnYHOltsW5ATs9vtYd1bkWO5NWW6yFpa7jdRJtp3MNbysfaOrcbi1C\n12i4PhdrwcWEcGNWS77xk6SD5MRPkg4yUVHfBbhQcYXFGRWXnMnJiWtsPlFxsHa3lbaTy3SmG6XW\nK+ycc85pnON2a1u43ezhB/hddtwW7QO3Rced2627HZ34WWvGVVWG1SW3Y9Ophgpf554dVS3Yq1Cf\nJW63jqcbF77XTk1dicAf+cZPkg6SEz9JOkhO/CTpIFPV8d0OLudGyjqO6pYu6g3j9NVat9/l2qku\nmYwLVMljsVT46ZryFdYftQ+1Y6ZrA9zOUQKNMnr/uL8ulLnTgfU5c+HY3fqNO+dctbndbkeqC4iq\nLBcefVTyjZ8kHSQnfpJ0kImK+s60oWIPi4AuyKOL5+680EaJ9e7MgC5oBnsROlFY2+lMU858yIyS\nMYbLXCqbTBtO/HQirctaxPddx9ON9XLZeZYqX9upzxmX6e67CxTjVAS9ju+Dy1q0ErkI8o2fJB0k\nJ36SdJCc+EnSQSaq4zt91f1W9TLG7Zpyu+yc+7Digkq2lQ/4RBVOP2ecvqp9dxljuX4dT+cq6rLz\nOmp32blnYhTXbLfuw/Vr3/l50XvCbs/6rLLO7XZJKjyeo0QKqs3tWMuyb/yI+HpEzEfE4/Td+ojY\nHREHI+KBiFg3cs1JkkyNGlH/GwB+R767C8CeUsoVAPYCuHulG5YkyeqxrOxWSvl/g1TYzK0Abhgc\n3wvgQSz8MbA485qKZyyauqAEKoLVij0uzr0LdOBMks4z0QVdUJHPBVZ0Of6c2OpMmVyO885TdYVF\nU5cbQHH3j81Yo+Syc2ocn+PdhlqmU6tcoFGXw8D1wYnsbtei1ucCx7Yx7uLeplLKPACUUo4B2DRm\nOUmSTIGVWtyzr9nDhw8Pj3u93ppKPZwk7xX6/T76/X7Vb8ed+PMRsbmUMh8RWwAcdz+enZ0ds5ok\nSWrRl+rc3Fzrb2snfgz+LbILwB0A7gFwO4D7qwoxZg63g0oj4nAQS6ePuzxpLre8M7OoPsVl1uq1\ngI/Ao3o943TZ8847b3iswT2dCa3W9OZ0S7dDbZQEJWxucyY0p3OPole7+8f9U52b69dnwq3RMKMk\ninH5Bsehxpz31wD+EcCOiHguIv4QwJ8D+HREHATwqcHnJEl+TahZ1f83LaduXuG2JEkyISbquae4\n3U8sHqrYyiLgKAESuQ7d/cQ4MdKZfJwJRttSKwq74BCqAvFvnYg5ismO++48KF3gz1E86d58882q\n61weRmdec0E6xvVMdOqmE+ddmmztX22wlFrSVz9JOkhO/CTpIDnxk6SDTFTHd1FTnBlJ9XGXy87p\nwM7FlK8bJe8c16HJL2qjpri2jJJXj9cmtO+uLXxftA8uQo3LuddWPtC8t7oW4XY0uqhJ3PdR9GF+\ntlziFl33cZGQXL5BtybkzHS16wa15Bs/STpITvwk6SBTFfVZrHNBJRQnfrpcaM7EVJtbTtvp1Af+\n7HL1KSzKaQ48F7zk5Zdfbq3P7SJ0YnJtunIXhMSpOaOI804F4rZo/1xOwVpPT20X3z+Xp0Db4p4z\nRq/j+1mbs8CRb/wk6SA58ZOkg+TET5IOMlEdX3U91r2c+6kLmulMIi7Zhsud55InuMQfzsXU5Ulz\nu8LUtMjjomZO7pOa5bg+5yLsXEVV/+f1Bl2z4D7pugT/1vXdmdAUF4yS63PBWfVZqjWhuYQoSm3U\nIrfGpXOlNnArk2/8JOkgOfGTpINMVNR3ZggOrgE0RTAVbVh0U5GIPcGciKll1qbJdl5vLnCkeqix\nCK8in9tNtm7dryKZv/76641zx4//KhDSpk3NMIjcvw996EOt9al5zZnluC3aTm6by9X3wgsvNM7t\n2LFjeKxqDj8jOmZOFOdzzoTmdl4qzszpzLi1QTOdKupSyteSb/wk6SA58ZOkg+TET5IOMlEdX81P\nrO84s5wzI2mZmjChrUynxyu1iRxcW37yk580zrEOzroyALzxxhvDY+fSqn1Yv3798Fj1fxcw1EWh\ncfrxK6+8Mjy+8MILG+fYnPjBD36wcY4jKnGAUKDpUqt95+fAuey6HZsKj6euM/E6kEaBciY75z7M\nbVMTKP/WrUG59aJa8o2fJB0kJ36SdJCpmvNYFHbeSOqFxqgawKKUy6HmrnN5/PQ6FlUPHTrUOMei\nHIvhQLO/LsCl7iZ79dVXh8cqmnKgSlV5Lr300uGxmh1ZdFSVhPun949F/6NHjzbOudTbrBJdcMEF\nrWVu2LChce7EiRPDYxe8RMVy59lWG6xVy2BR38Xx13byfXHPtc4HVpdGUVPbyDd+knSQnPhJ0kFy\n4idJB5moju8iv6iJgnUjp/+7qD5qLnGmMNYtVZfl+s4///zGuSNHjgyP2QwHAFdeeeXwWM0/bG5z\nySgU7pPqj25NgX+r9fEag+qdbpek2/XGY6Z6PK8xqGsx6/FqHt2yZcvwWPvO5kTOzgw09WrV1bl/\nOu7uWeLPaubkcvQ6Pud2bOp94HaPY75TanLnbY2IvRHxo4jYHxFfHHy/PiJ2R8TBiHggItYtV1aS\nJGuDGlH/HQBfLqXsBPBbAP44In4DwF0A9pRSrgCwF8Ddq9fMJElWkpqkmccAHBscvxkRPwawFcCt\nAG4Y/OxeAA9i4Y9BKy5VtJpgandGqejG4pOKbm73kwuoyfU988wzjXNswtu4cWPjHAe/VNGbURWB\nxUP1bOMxZPMd4PPxsUir/WO1Q8V5VsFcvkE9x2Kr3gdut/aBRVwVadnUqO3kOrZu3do49+yzzw6P\nVSx3nonONOzyFNTuHnUqnZbJ6tjEzXkRcQmATwD4JwCbSynzwPCPw6b2K5MkWUtUL+5FxHkA/hbA\nlwZvfo0N1BoraHGxZWZmBr1e7xTHjCRJTp9+v49+v1/126qJHxFnYmHS31dKuX/w9XxEbC6lzEfE\nFgDH266fnZ0F4MX3JElOj16vh16vN/w8NzfX+tvaN/5fAniylPI1+m4XgDsA3APgdgD3L3FdA3UV\nZZdT1WXZBOP0VcVFP2H3V2cGVF2P1x9U72TpRSPbcH+176zzaxQat4PLRTHiMdMdf9xu1Tu5f04a\n0z44N2cee70P7H567Nixxjnun7aT10J0NyCf03vE6wjqAv3iiy8Oj/UZ5HULfeb4GXGu4S6waa3Z\nFvA7Gsdh2YkfEdcD+LcA9kfEPiyI9F/FwoT/ZkTcCeAIgNtOuzVJkkyEmlX9fwDQJqPfvLLNSZJk\nEkzUc09hsVzNQSySuZTPLiCCilkuLrsz53H9KkKz2OqCSjjvLvVsY0YJAOFMYfxZ28n9VVWCy1TT\nm4uPz2ZW3lEINO+fjgurf7r70AUadTviPvaxjy3ZLgDYv3//8PijH/1o4xyrASrqc39Vbaz1slMv\nUFZXnBlQd146M2sb6aufJB0kJ36SdJCc+EnSQaaq47sceE7vbCsDqM8/rrufWOdW/Wrbtm3DY05a\nAdTnV9fdh9w/bSdf5wJ4qg7Murvqlly/9o91Teeaqn1gVMfn/qk+zmOm6xt8jm3Ses6Z+pwpTOH6\neTcl0HQRdqY3l3PPuYZrH/iZ1Ofa7a4ch3zjJ0kHyYmfJB1koqK+imAs4rp49YpLIcwikYrzLl8d\ni3Vq8mGRUz2/nBuyU2UYd05FTGfO4zHUseb+Om9AVS1cCm2X39C1xXnn8dirWsX1q2mWVSAt8/nn\nnx8e845JoKnG6a4+DhLy9NNPt7bT7Ux0OfBcTkiF77s+g06lbSPf+EnSQXLiJ0kHyYmfJB1kojq+\ny8vmXBudHu122blglKoX8RbGzZs3N86xju/Mcqp38me3NqA6KX/WdnIfXOII1R/Z9Ke6JLvi6n2o\nDe6pbXERePgeuR1/Liiolsn3RdvpEniwHr9r167Guccff3x4fPXVVzfOcbtdchY9x/dIn2vnUs7X\nufWvWvKNnyQdJCd+knSQqebOcx5ILp47i3kqQr/22mutZbKIqW1hkY9NPIAXy1nEVTMge9KpaMpi\nnebOY7FVU0yzeK3nWKxUkZ3763LnaVsYvV/cdxX11auwrRw1hbFqo+1kD0NVVzTnH8MBUjg2v7bF\n5cDTICucK9DlFHTenFofj4WqsG734TjBN/ONnyQdJCd+knSQnPhJ0kGmujvP6dys02jkF9ZJVQ90\n5hKuj3O0AU29XvVO1s/VvMZlql7GbqSqc7tglLXRiDUSC+uMTsdWHZH1em0Lr2mMopPy/VTXaV77\n0PvO97rWHRrwEZvcfeC+632/6aabhscaeYnb6UJaj2J6c7/l8XVJZGrJN36SdJCc+EnSQSYq6qto\nWuttpeecCaY2fbGaQFjk0/pY/FW1g+vQ61555ZXhsQa4vPzyy4fHqnZwHSq2stjsPOk0+AWPp5ok\nWV1S8ZpNhnod/1YDR/KYqTrG5Wg7+RlRExrj8tzpM8H1qVjMATxUZOe+qxmQ76fmRWA1y6UPV7WD\nz2kfnEo5DvnGT5IOkhM/STpITvwk6SBTjcDjEmqwXjaK26PLW8461GIG36Xq07YcOXJkeMzRXABg\n586dw+ONGzc2zrkcZ6zHq3nG7WxjnVTNZNxfNT85F2GuX3Pn1QY61TJdVBjn6sv3U92vuT7XLh1P\n7rs+L6y7u6CnaubkduvaFY+9rm+4dSa+f6rHu6hT47DsGz8izoqIf46IfRHxo4j4s8H36yNid0Qc\njIgHImLdcmUlSbI2WHbil1LeAnBTKeUaAFcB+O1BIs27AOwppVwBYC+Au1e1pUmSrBhVon4pZVGW\nOQsLfyxOALgVwA2D7+8F8CAW/hi0oqINiy9qmnJx4FkkU/XBmXV4Bx6nqQaAxx57bHisO/5YHLz4\n4osb51jcdjHwVe1gM5bzTFT1weWdYzFS2+KCWHDbtC0sNruAIdo/FlVV3OXx1TL5nAs+qarFSy+9\ntGQZQHMMNZCKi53P91Z3/7GqoaI+48R5F2jE5c5zcfxrqVrci4iZQYrsYwAeLKU8CWBzKWV+0Mhj\nADa5MpIkWTvUvvFPArgmIs4H8EBE3AhAVxtaVx94Ia3X652SISVJktOn3+/bfQPMSKv6pZTXI+Lb\nAH4TwHxEbC6lzEfEFgDH266bnZ0dpZokScZAX6ocR1JZduJHxAUA3i6lvBYRZwP4NIA/BbALwB0A\n7gFwO4D7lyvL7cBz0V3UPDNOPnAAOHDgwPD42muvbZxjM50G23Quu6z7aUQcPqd6Lpv6NH883zzV\nZfmcmsJqE3io/sh1OHdXXYdx5kpew3DJUtTVl3Vi1ce5LXof3I4/HotLL720ce6JJ54YHu/bt69x\njtdz2P1a+chHPtL4XGty1TUFt57Cz6Ded7cW0kbNG/9CAPfGwt2bAXBfKeV7A53/mxFxJ4AjAG4b\nufYkSabCshO/lLIfwLVLfN8HcPNqNCpJktUlVtoj6JQKIsott9wC4FTxxYn6zlTkznE5LsiDmsJY\nfVAR9plnnhkeqxnQxYjnc+r55eLcs/fc/Px845yK24zb3cXXOXOey3OnY8Z9UHXM5SnksdA49ywm\nq0jLor/2gb389D7wLr+nnnqqcY5NvKpCsgfeZZdd1jjHufRcSnKXi1CfTx5DVXNYHdQy24Jt7t69\nG6WUJW196aufJB0kJ36SdJCc+EnSQSa6O0/1caevumQULumDyx/PZbpAjlrG9u3bh8caMYbNNXqO\ndWJ19eVzH/7whxvnWMdXV1F2D1Udkfuk6xSs96opjNdMVF9lE6W6wrI+rue4nbpu4ExT7v5xu9VN\n1rlH8xhqW77whS8Mj7/zne80zrldklyOM6e5vIGjuLCrqZjJYJtJklSREz9JOshUc+cxKp45bzL+\nrcurp6IUi7sayJHbptexWK6BMV0ADxa3Xb4zlwdOTaBsQnMecTrWXL96k7F4r/H4eTy1Phbv1VOQ\ny1HVgn+rIizfW20L90nHk82seh/YFKZedg899NDwWNUc9vJT91c2H+7YsaNxjtUXZ8LW54zHV02n\nzkxdm4eByTd+knSQnPhJ0kFy4idJB5mojq9mONZtVWdjvUXdQdmdUfUdlyuczWSqP7Je5vKyqX7M\nOpu617KOqC6YXIeaAWvdaxXnPsxlqumNd/xpfdw21ce5T6ofM24tQs2VvGaj7tEuKhN/VlMfj4U+\nS9xffSZ4PYfdtoFmVB/tA98jd7+0DzxO7jkeR6dX8o2fJB0kJ36SdJCJivrOLOdi5ysstrqceyru\nsgioUYE4WKN6k3E5Kn6yuHbRRRe1tlPVHO6vBnlwu/pYNNVzbvcamy+dWVVNi6wW6P3j8dS2cB90\nBx6bqtj7b6lyGB5DFztfueSSS4bH2nc3ZvxZg57qbxnuu6pA/My7gKgudbqru5Z84ydJB8mJnyQd\nJCd+knSQier4zkThzCyqH7NerXon64+qJ/G6gbresg6su+W4/ueee65xjnVLdQdlHdi57I6yu4rH\nSfvHbdEy+bdqDmI9e5RchE6XZf3VRfVxZs6XX365cY7NiboGxP3TnYkuaQajfeD1HI7UA/g1BX4G\ndTzdjj+XO8+dc/OqjXzjJ0kHyYmfJB1koqK+iqYssutuJLcDz+1UqhVNOfU10BQP1fzkUlqzN6AG\nZHz00UeHx5rhpNYDT8fFxV7n36qpyAWxcCYm/jyKGZDFWBVp+R65nWY61k7cZQ9DrY898jTNOY+1\nlslt03j8bP51Yr/CfXIiu/aB763LN1BLvvGTpIPkxE+SDpITP0k6yER1fGeiUF3W5VBnHVH1QNaT\n1ITGn93uPG2nc8tlM6Drn+rA/FnNlVyOXudMaNxOlx9PdUSuT/vgXFr5nun94/pcrkM958bFBfe8\n8MILq+pTl2te69m2bVvjHD9bx44da5zjNRM1HzrXW75/us7DawpuXUvNd6sabDMiZiLiBxGxa/B5\nfUTsjoiDEfFARKxbrowkSdYGo4j6XwLwJH2+C8CeUsoVAPYCuHslG5YkyepRJepHxFYAnwXwnwB8\nefD1rQBuGBzfC+BBLPwxaK/MeNKpSMtirEs17MyAzlyi6Zn5tyoqsqlIy2SvsBdffLFxzqkPztzF\nqKmIx1C9wlyZLpgoi6Z6H5ypysXHZ7VKy2AR2uU35Nx1QPO+qBrAKgkHFgGaac9VxXNBOjnPntbH\n5Wgf+L44U6bzaNTnheeA7mgcJ2187Rv/LwD8CQBWHDeXUuYBoJRyDMCmpS5MkmTtsewbPyJ+F8B8\nKeWHEXGj+WmrF8Hhw4cXy0Kv1zvlL3KSJKdPv98/xVGsjRpR/3oAvxcRnwVwNoAPRsR9AI5FxOZS\nynxEbAFwvK2AxaAXKpomSbJy6EtVcwEwy078UspXAXwVACLiBgD/oZTyBxHxnwHcAeAeALcDuH+5\nslw+cBeBR01TziTCf1x0bYB1Ut35xfqkmsl4HUEDatb2wZ3TcdF2MzwWoyQo4fo0wQWvd7g1Ewe7\nLmt9uhOS74Pq3Hz/9L47Pde5arN7rQYM5fUH1ePZtVl1bhddifugbXEuuy4ZjDs3DqdTwp8D+HRE\nHATwqcHnJEl+DRjJgaeU8hCAhwbHfQA3r0ajkiRZXaYaiINFVTX5OBGTxSwnEiksOqrHH7dFd3Bx\nmbpzj8U89SZjE5MLIqmmMEbFfm6LiuVcn5bJpiMVd/mc2y2noimbkXRcjh492tpODnSi94tz0u3c\nubNxjr0kNUApi/5uh+H3v//9xjl+llwAFjX/srqka1e1Jlc9x5/duYyrnyTJWOTET5IOkhM/STrI\nRHV81U1Y33FRRNS85nbuOdzuJ5dfjfOmqe7lAkDyzq8tW7Y0zrFOrDopj5NGy3G57FzASW6brkXw\nOd3RyGOvO9u4fjXL8bqBXse/VYeT48d/5Q6i1zFXXXVV4/OhQ4eGx3pvuRxdS+Igmq+++mrjHD8H\nLpio3iNnWmz7HdBcr9K1K+7TuCbXRl0jX5Ekya89OfGTpINMVNR3nmYur55exyqC83JTkYjFMxXL\n+bOK0FdeeeXwmM1N2hYVdzlnm4qtbH5yeeBUjGMTkxNb3Y4xHWv2WtQgHaxa6DntL8Njr/1j8V7L\n5KCWev+4v5/73Oca5zgQx/e+973GuYceemh4rPHx3f1jXHAWVdV47EdR49QDkOG+666+DLaZJEkV\nOfGTpIPkxE+SDjJRHV9hs5LqOy5iDOOi7Li8bFqmC/LIOrHuQmNTn+pznEuPd4gBzb6reY1dUznC\nD9A0Aam7K/dXXWh5h5xex+6vOi78Wc1dbDJU3ZldXHVdxOUU5PUVHc+tW7cOj7/1rW+11vfEE080\nzj388MPDYxfNya37qB7Nz4SuYfD6jXNF17UW1d0ZrkPHbJzt7vnGT5IOkhM/STrIVD33WJRScckF\nqqgNqOlSaLtdU9oWF9P84x//+PBYxXneIacptFl0fOGFFxrnWFTUMWPTm8s3oCItm+VUlXG5CHgs\ndFwOHjw4PNb7wOZLDcTBaoGeYw85VdU4tv2ePXsa5/i32k5WGXTXosvR6DzpXP4/Ftldmmwtk++Z\nego68/Y45Bs/STpITvwk6SA58ZOkg0xUx1cXTEb1FucSyeXodU7Hd2ZBLlPNa3ydSwCh7WRTnOrO\nHOzT5cDTXXaLEYuBU81r3G7tq8uB53bn8W+175xARN2ceU1D7wOPmQuWqvo478BzOeK1f1yHjif3\nV/vO/XWu4WqKZvQ+8DOi48lrL9oH/qzPma6v1JBv/CTpIDnxk6SDTFTUV88kJ6I4sZzFPBXPuEwX\n+EPbwmWqSuLa6bytWKRVcZ7VAO0r1691s/ir/WMRUM1krFq4dN7aFr5u48aNjXP8Wcvk+t2us1FQ\nMx3D99bFq9f75Uy1Lkej2xHH59yzo+PC/VN1kz0AxxHtlXzjJ0kHyYmfJB0kJ36SdJA1E4FHcQEE\n3c4o5xLp9DKXA88l6WDdzwVBVH2O9WMt361vsOut6qu8pqA6vtu1yOskOp6cZELvH++I03O840/L\n5HHRNRpew3BJT9S91u2kc4EqXR5GZ+bkMt06j7bFRVdy61Ncn54bx4W3auJHxLMAXgNwEsDbpZTr\nImI9gL8BcDGAZwHcVkp5rbWQJEnWDLWi/kkAN5ZSrimlXDf47i4Ae0opVwDYC+Du1WhgkiQrT62o\nHzj1j8StAG4YHN8L4EEs/DFoxe1iUlg0VrGOxTUnzo9iPuQdcW6XluJy4jHq+cVj4UyS6qWlOdwY\nZ7bi+rU/zkToRFOXz43VCVVznFjO/XM7L13qdH3O+B45dc/t2NTralVRlxremR3dM6jnxhH1a9/4\nBcB3I+KRiPijwXebSynzAFBKOQZgU+vVSZKsKWrf+NeXUl6KiI0AdkfEQSz8MWBa/5QePnx44Qel\nYP369TZDSpIk49Hv90/JTNRG1cQvpbw0+P/liPg7ANcBmI+IzaWU+YjYAuB42/WLG0tWIoBAkiRL\n0+v10Ou/v0WFAAAJ+UlEQVT1hp/n5uZaf7vsxI+IcwDMlFLejIhzAdwC4E8B7AJwB4B7ANwO4P7l\nylKdhvVH3TXldDa308y5S7JepkEQWQfWc6xTuWhAaiZz5iDW9fQ657LL9aneyWVqO9ll2PXBjafq\ntXyPdMy4Lc785AJHOh3f4daOnHu0e85cIE6Fx8ytRYzy7LpdkrXjwtS88TcD+FZElMHv/6qUsjsi\nHgXwzYi4E8ARALeNXHuSJFNh2YlfSnkGwCeW+L4P4ObVaFSSJKvLVOPqO5HdmfMYl0fMecSpmcx5\nRtV6Cmp9bucXi5hOhHY58Jy4qaYpHhdnYnI71FRs5QAUupuM75+L1a/j6XLusbqi973WrKrXOS9C\n1wfGBXzRZ4LrcG126thKiPrpq58kHSQnfpJ0kJz4SdJBJqrjq9nKuTYyzqzjXCddUEmXIMGV6Xbg\njaIHcluce6ZzO3aum871VvXxtt8BTTOn3gcuR/V/zjWvLrvcbrcO43RZHWtum7bTuew60xuj9682\nwYXT47U+vtc6LnwfdDxdu9vIN36SdJCc+EnSQSYq6qu4pOYTxomtTsxyJq6VSL2tfWCzoMZXZ3OU\n87JTnImJ63cBQkYJ7uF2A7Jq47zQnOo0Sv44rl/VjtrAmDpmzjvPmTL5t848qrhdkq4MN2YarJXJ\nNNlJklSREz9JOkhO/CTpIFMNtuncT9n0p7pQ7fbeWjdOLdPtJnO61yg6G+PMT27nnjNp6XWsd7q1\nFS2T26I6N/dX++fWU/icc9l1bqsuapHed65DdWUuR9c33C5C93y6NSgeJ73OPSNcv9sFWku+8ZOk\ng+TET5IOMlFRf5RceU4kcmKdq89dx6Kwi3eu5hkW5ZwZ0KVudiK7E2ndeGr/eDy1TO6TivOMC2Kp\ncP+0nXw/tT43Lk614N+qWZX758yoToV05lBn2lOxnOtw8f+1Dy4YTJrzkiSpIid+knSQnPhJ0kEm\nquOr7qw73doYJSmBq49/63bZOVdYrY/bpvocl+l20o3i7splutxyuoPLXce6ppq0WAfX/tVGfnFr\nNM5U69Yp3M5LxY21M725QKouUhCXo8+ZWy9yrsz8W/cM1pJv/CTpIDnxk6SDTNWcVxt7XUUZ50nn\nTBssOqoo7PKysfirpjC+rlbcBOrj4ztTmDNXqmjqcCZCZ0Jz6bVrg4KqasH3ZZRY9rUitFMvXept\nVS14fNX05lQg91w79YHrcCbJWvKNnyQdJCd+knSQnPhJ0kGqdPyIWAfgvwP4GICTAO4E8BSAvwFw\nMYBnAdxWSnlt3IaovuPy3DtdiPUr3Ynl8tUxTs8dJXKPc5N1ZqTaXW+jrG84/Zj1V+0fB9R07sOK\n29nm1kVYzx4lbyBf59yxnXu0C8SpvPHGG1Vl6vPJbdNx4fUOrdvtBlzNYJtfA/DtUsqVAK4GcADA\nXQD2lFKuALAXwN0j154kyVRYduJHxPkAPllK+QYAlFLeGbzZbwVw7+Bn9wL416vWyiRJVpQaUf9S\nAK9ExDew8LZ/FMC/B7C5lDIPAKWUYxGxabmC1GTHoo7zPtLrWPR3sdfdDi7FicmsWqja4erjz2q2\nYvFTr+Nzeh2Pmcsb6HYK6jkuR8VPt/uQxdZx1Q53j1ygSudd6Z4J7Z/Lp8Biul7HfXe7CLVMVj9V\nJeHr1OzIor4LbFJLjah/JoBrAfzXUsq1AH6KBTFfFYvRFY0kSaZCzRv/eQBHSymPDj7/byxM/PmI\n2FxKmY+ILQCOtxVw+PBhAAt/3Xu9HjZs2HCazU6SROn3++j3+1W/XXbiDyb20YjYUUp5CsCnAPxo\n8O8OAPcAuB3A/W1lzM7OAhgvnW+SJHX0ej30er3h57m5udbf1rrsfhHAX0XE+wDMAfhDAGcA+GZE\n3AngCIDblivEmWfU7MG7wpy7q3OzdIEc9Y9QrfuwM50412K308wlqlBdtjYXm0tiMW5gTKVWj1fT\nqRtDlxeR61NTLd9blxvQjfUoEaIY95y5aEfO/KuME1DTUTXxSymPAfgXS5y6eUVbkyTJREjZO0k6\nyER35yksdqlIySKRC+QwSk4zl8vOBUGo3W3lrnPeVs6k5URTxXm21aaAdjHiFScmuzJdu5zHH5fp\nTH2uLdqfcYOJMG7H3ygiOtexEjvwHPnGT5IOkhM/STrIxCZ+rX1xEqyVtqyVdgDZljbeq235tZz4\npZTWf46IQETgxIkTw+Ol/rVdN2riguXqOHHixFjXjcNyZbS1ZSUY9X61PSu193klyYmfJMl7hpz4\nSdJBYrXFpojIzTtJMiVKKUvqdqs+8ZMkWXukqJ8kHSQnfpJ0kJz4SdJBVn3iR8RnIuJARDwVEV9Z\n7fqk7q9HxHxEPE7frY+I3RFxMCIeGEQQnkRbtkbE3oj4UUTsj4gvTqs9EXFWRPxzROwbtOfPptWW\nQb0zEfGDiNg1zXYM6n42Ih4bjM3D02pPRKyLiP8VET8e3KN/uZLtWNWJHxEzAP4LgN8BsBPA70fE\nb6xmncI3BnUz04oO/A6AL5dSdgL4LQB/PBiLibenlPIWgJtKKdcAuArAb0fE9dNoy4AvAXiSPk8z\ngvNJADeWUq4ppVw3xfasbmRr51V1uv8A/CsA36HPdwH4ymrWuUQbLgbwOH0+gIVAoQCwBcCBSbaH\n2vF3WIhnMNX2ADgHwMMAPjqNtgDYCuC7AG4EsGva9wjAMwA2yHcTbQ+A8wE8vcT3K9aO1Rb1LwJw\nlD4/P/hummwqFB0YwLLRgVeaiLgEwCcA/BMkWvGk2jMQr/cBOAbgwVLKk1Nqy18A+BM0g7VOZUwG\nFADfjYhHIuKPptSeYWTrgQr03yLinJVsRy7uTTg6cEScB+BvAXyplPLmEvVPpD2llJNlQdTfCuCT\nEXHjpNsSEb8LYL6U8kMAbiPCJO/R9WUhmvRnsaCOfXKJ+le7Pase2Xq1J/4LALbR562D76bJfERs\nBoDlogOvNBFxJhYm/X2llMXgpFNrDwCUUl4H8G0AvzmFtlwP4PciYg7A/8DCWsN9AI5Na0xKKS8N\n/n8ZC+rYdZj8uCwV2fralWzHak/8RwDMRsTFEfF+AJ8HsGuV61QCzbfJLixEBwaWiQ68CvwlgCdL\nKV+bZnsi4oLFFeGIOBvApwHsm3RbSilfLaVsK6Vsx8KzsbeU8gcA/n6S7VgkIs4ZSGSIiHMB3AJg\nPyY/LvMAjkbEjsFXi5GtV64dE1gs+QyAgwAOAbhrUos0g7r/GsCLAN4C8BwWogOvB7Bn0KbdAD40\nobZcD+BdAD/EwiT7wWBsepNuD4CPD+rfB+AxAP9x8P3E20JtugG/WtybSjuwoFsv3p/9i8/rlO7R\n1Vh4cf4QwP8BsG4l25G++knSQXJxL0k6SE78JOkgOfGTpIPkxE+SDpITP0k6SE78JOkgOfGTpIP8\nf/gTbRJN7cdYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10df4a050>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(white_betas_4d[:, :, 15, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x112760090>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD+CAYAAAAalrhRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnWusXcV1x//LOOEZO07ANsHBbUTspuQeHmppiBNhGgIR\nlko/oZSbBor6IVKrpElbBaJIUT+0Kv1SJWq/VE0TQu6lTZ+4cQrGtYzUlrZYOJwbgw2E8rYvEPOI\nSWIenn64xydrL+7+n5nt87hk/3+S5X3u7D2z9uwzZ681a80aSylBCNEulk1aACHE+NHAF6KFaOAL\n0UI08IVoIRr4QrQQDXwhWshxDXwz+5iZ7TOzB83s88MSSggxWqypH9/MlgF4EMBHADwN4B4AH08p\n7RueeEKIUXA8b/yLADyUUnospfQqgL8FcNVwxBJCjJLlx3HtWQCecJ+fxMKPQQUzU2igEBMipWSL\n/f14Bn4211xzDbrdLjqdDqamptDpdAAAZovK9AaiOXL06NH+8QknnFApW7bsp0rM66+/Xil77bXX\nAAC33norrrnmmkqZryfK5dv3bQPA8uU/7cIjR45Uyt761rfWynLCCSfgm9/8Jj7xiU+8oczfA4PJ\nyc5drL2ZmRlMT0/jLW95S6Xs1Vdfra3T90WU2X+OfeY/L9be7OzsG54PUO3reK/sGcXvSO51vl8i\nvj+jLKys7rxBZcuXL8c3vvENfPKTn8Qrr7yy6DXdbhdzc3P9z7Ozs7X1H8/AfwrA2e7zut7f3sCx\njlusA4UQw6HT6fRfqsDoBv49AM4xs/UADgD4OIDfYBfEX7HcX0ZGfHvl/sLmvgWAN74J6ohvLw97\ni+e+4aMs7B4i/p7idb6MveHj8/Nyx3tgz4HJvWzZMpgZli1b1tfQcvCyRVmG8T0bN4tpc8f+DYPG\nAz+l9LqZ/S6A7ViYJPxqSumBuvOnpqaaNjV0loos/td50iyVPgEkSx3nnXfe0Opq7M7LbsAsbdu2\n7dhxLOsfMzlK3sDs7em1g5K3JWvP253xvGHYgUwW9uZmlMwN5F7H3vjs3uN1/hnFN76fMyn5vvg2\nSuYG2Hdp1DY+mzNhWplny5Ytk53cq8N3dEmns4kj1ulscqhOrnhdLPNfTiYzuz/2w8ImL9m5TE52\nHRvc7Ase6881X+L9+YEfTadcWSLsGeXKGU1KZjZ62Ur60z+zeN1PfvKT2vaaoJBdIVqIBr4QLUQD\nX4gWMlYbv2SShwWHeNikErOrc4OHYp2D2h/GdV7ukok/Zj8yu9CXRVs2dyKV2c5sQi1OVPn5lBJX\nrT+3JJjIw+xxL1csKwnA8ucy92icMxmGXV9pa6i1CSHeFGjgC9FCxqrqM9WbqTJMLS8xH3KvY7KU\nmA8lLsM6WJ0lfmDmc2dmVdO4AeZC83JHl11dHPogcs2VqLJ7ck0J1nash7kBhxXN2QS98YVoIRr4\nQrQQDXwhWshYbfxcN0c8N9pezHXD3HlNYW4dZvsxezJ3lV3JPAGbb2Dh0bluQGb/x37wtnt02bGy\n3FwL7B6YO4/NmcT2cuc+4r2z585Cyll7LDy6iatPb3whWogGvhAtZKKr81g6KIZXbaK65FXHE088\nsVI2DBcaU4WZG4lFabFlpCyNFFvZxlbLlSwV9f0ZVVivsjNzLJK7GpHVWbKaM9esKvkOsii73OjA\nErOq6ZLyOvTGF6KFaOAL0UI08IVoIUsmZJedW5LIkdmdubKUuFly7eWScFAWPszseG+DRzmZ7cyy\n87AEok2z+vh7iPMGuRmNSkJaWcp1FgKd20a8h9y5pDgH5dOzs+/EMNLl6Y0vRAvRwBeihUw0cm8Y\n0WslEVUscQRT55mazK5jKi1bMeavY66+iC/zyRkHycnq9Ko3S0bB3IAl7TFX5jB2GMrNnBuvi+R+\nX9gKyrjz0jjRG1+IFqKBL0QL0cAXooVMNGSX7Q7CVjgx15u3qZhNyuyyODcwjESOJSG0bE4hd3VX\ndMM1Xano62HuwwgLaWUblDSFzRflziWV2Pi5cxGxTrbpCSM3g1IuA9/4ZvZVM5s3s6772yoz225m\n+83sDjNbWdyyEGJi5Kj6XwNwRfjbDQB2pJQ2AtgJ4MZhCyaEGB0DVf2U0n/0tsL2XAXgkt7xzQB2\nYeHHoAivsjC1cRGZFj1m50WarsSKeDXLb+oY24+mTG7e+RLYSjpmPjD1M3djzAgzO1gduckvIrmr\nD0vKWL8w84/hXXglG2Oy59eEppN7q1NK8z0hDgJYfdySCCHGxrAm9+hP0MzMTP94ampqSe0LL8TP\nCt1uF3Nzc1nnNh3482a2JqU0b2ZrATzDTp6enm7YjBAil06nU3mpzs7O1p6bO/Ct9+8YWwFcB+Am\nANcCuC2rkoL96nJXxJXsZe9tzRguOYyVgrmr1WKdTTfpiLLkbhZR4v7J3d8w4mVhq/NKNi/JDZ1m\nSTrZ5hfMjo/3wPo693vA5naahrDnkuPOmwXwXwA2mNnjZvZbAP4UwEfNbD+Aj/Q+CyHeJOTM6l9T\nU3TZkGURQoyJsUbulWwnnKvWRVcRc9N5F8kpp5xSKfOr2UoSMniVL8rp94GLrj4Pu4eSZCJNo/p8\nv8Qy5jplKntuQg3mtorX+TpLkmayRBy527GXmKm56jwz45j5x5KX5KJYfSFaiAa+EC1EA1+IFrJk\nNtSI5Lp8mtpshw8frpWlxMZndplPpsjszqbhmUzO2F6ui6kkNJVlmmGrHVkmJIbvT2b/R/w9sZWC\nLLEpe0axTvbcc8nNitQUvfGFaCEa+EK0kIkm28zdC42pS0zdLXGJsFV23i1XktyDJQVh6trb3va2\n2vNeeOGF2vZYNJlXP2Of+TbYCjUGcx9GmDnGXHb+OZQkS2WJTX1fRHW+7rwIq5NFLZaYOSyCsok5\noTe+EC1EA1+IFqKBL0QLGauNX2In5e6h3tRdEslNDhnLcl1azH6M9/ryyy/XlrHQX3/uqaeeWltn\n3LPNy83coxHmdswNf4394mVjdbLvUlN3bO4KzUH4Poz96Z9f/C7luqKH8Z3XG1+IFqKBL0QLGauq\nH1Ubr5KV7JPm1bOmWxQzYh3M1cfccr6M3V9JEstc4t55nqgqsqQZvp6TTz65UsZcbyyi0ZshUU7W\nZx6WwDOSu8V0rlkDVOUsMS2GkYClJP9/HXrjC9FCNPCFaCEa+EK0kLHa+Gwf+JIQTBbqm7s3GluF\nxurwYaOD8PMB0S5j2Xl+/OMf94+j6833xV133VUp27x5c/+YbeBR4kZiqwj9dVFO9ox+9KMf1crC\nkqWyMO7cDTVYnWz+Jtbp5xhiv/jr4lxEbuLPCHOBlmxG05ej+AohxJseDXwhWshYVf2oyjAVxbt5\nWLQaSwAZYSv3WJLH3BV4UXXz9xDl8qpxNB/uvvvu/vHFF19cK+emTZsqZf6eYqKRVatWLXpeJN6f\nV1VZ30a3nL+nKIt350UT4cknn+wfr127tra9+N1hK/c8TRNcsFWELGKz5PuZu5qzZFVfHXrjC9FC\nNPCFaCEa+EK0kIlm4PG2ZrSTvF0f7SRfxtw6LESYua1KQjd9+y+++GKl7Kyzzuof//CHP6yUPf30\n0/3jaMv6spdeeqlStnLlyv5xnPvw97BixYpKmXcBMRsx7inIXJtspZkvO+200ypl3n6N7q7Vq+t3\nXPfPpWSFWu6qt0juKtCSOpkL238/vUsXqPbnWGx8M1tnZjvNbK+ZzZnZp3t/X2Vm281sv5ndYWYr\nB9UlhFga5Px0vAbgcymlcwFcDOB3zOwXANwAYEdKaSOAnQBuHJ2YQohhkrNp5kEAB3vHh83sAQDr\nAFwF4JLeaTcD2IWFH4Na2F5v0aWVm0ufqd6xLHdvObZ6LZZ973vf6x9Hlf2BBx7oH5955pmVMp8Y\n4ytf+Uql7Pbbb+8fX3755ZWyk046qX/MIumiCu2vi33t1fS4pyBbLZcbxRhVfX9dvAefaJSt3IvP\nwdf57ne/u1I2Pz/fPy5Z2caiD9nqQ5aAhUXu5UaylqxkraPoCjP7OQDnA/hvAGtSSvNA/8eh3jgT\nQiwpsif3zOw0AP8A4DO9N3/8qaz96ZyZmekfT01NodPplMophBhAt9vF3Nxc1rlZA9/MlmNh0N+S\nUrqt9+d5M1uTUpo3s7UAnqm7fnp6OksYIURzOp1O5aU6Oztbe27uG/9vANyfUvqy+9tWANcBuAnA\ntQBuW+Q6irdpWEYVtgFEtHdYYsxcl10Ml/QuLmYDR1tr3bp1/ePonnnPe97TP963b1+l7Pzzz6+V\nxcsdbURvS0fXYu79Rpu06XwK2z/Of45ZfbztzFahxf70z+W5555DHSy5Z1OXHds0g4WlR1n89ywm\nSx1GZinPwIFvZpsATAOYM7M9WFDpv4CFAf8tM7sewGMArh6qZEKIkZEzq/+fAOqm2C8brjhCiHFg\nw0jsSBswS9u2basr6x+X7L2Wq3ZF1Y25CNmKKv85RtJ5FTO6n5ia58v8fnhA1aUVTSDmsvMurdw9\n72KdTJ2P/cJWk7Hn4GErIb3LE6j2S1SF/XNhW4TH1YBsX71c9Tq25+8htuf7N7pDfXtsD8Ncl+SW\nLVuQUlp0AwLF6gvRQjTwhWghGvhCtJCxrs6LeNvIh5QCPByUJWTMTRwZr2NuuXvvvbd/7DPEANUE\nlzE09amnnuofn3766bWyrFmzplLm+4W5JKNtyVyivh4fwrpY+3V1sv3/mJyM+By8ey/en7dt2SpC\n1i8sW8473vGOStkzz9SGptTWEduIsrBVkt5WZ/euDDxCiEZo4AvRQsaq6jP3CFtpFhmG24q5MaML\nxqvC73//+ytlzHxYv359bZn/XJIwpO48oHpPsV/857hS0CftiO5KFp3H9nNjLjzmmmJbhPvPLHKP\nEe/BP9toAjEzh+1zx8yxXNmi2ca+r9o7TwiRhQa+EC1EA1+IFjLRZJsetq8ec5eUrLLztjQLMb3v\nvvsqZRs3buwfxzBSXw/be63EBcPsf7bBBbOd2bzBD37wg0XriOeyPeJL8Dap30cvysn6M8ribX72\nbKObePfu3f3jd73rXZWyc845p3/8xBNPVMpieLaHheWy726uO7ZkXqsOvfGFaCEa+EK0kLGq+kxt\njZFYbNvj3G2QWa73qKp5FT6687x6xlZbsSSIbGVbVNl9BODzzz9fKfN9xvaPYy5C5upjxDrZdtDM\nteiJ6jyT09cZ22OJOBlnnHHGovUDwGOPPdY/ju48trLUn8siBeN3nu1TMAz13qM3vhAtRANfiBai\ngS9ECxmrjc9cdizTTElIIlul5e16FuLJEmoy2CYIcW7A76UXVyZ6F1fsM18PWy1Xsmc7CzH19bCw\n1XjvuRt/RDmZu9ITnwn7Lvl5EraPX5z38XVGe9xfF58Rc936vmfjgYVjx7Im6I0vRAvRwBeihUw0\nEYeH7U026FyPV2mZ6h3r96pczPXuiWU+ecPDDz9cW2dUFf0+e1FOrxrnJuwEqskoY3s+sq3pNuAs\nsUnsT6+KR1Xfu61YNCBTaaMs/v5imU/MybamZio6S9LJYC7JSO53Xok4hBCN0MAXooVo4AvRQpbM\n6rwIC4lkeLs32rmevXv3Vj57G+rgwYOVMm+Pb9iwoVLGQn2964jtHxfDVv39MhdThO2hnusCinMY\nbDWZdzuylYLx/ny/eLcmUH1mcQ6D7Tvv22DZeeJ3ycsdnx9zl7IEnp6SjWLYeWyup8m+egNHopmd\naGb/Y2Z7zGyvmf1J7++rzGy7me03szvMbGVx60KIiTBw4KeUjgC4NKV0AYAOgF/tbaR5A4AdKaWN\nAHYCuHGkkgohhkaWqp9SOqbTnYiFH4vnAVwF4JLe328GsAsLPwa1RJWErWxjkUrMZGCrtK688sr+\n8a233lop89tWHz58uFLmVcVY5hM0xoix3G3Ao7rr1Va24i/2C1NNWUSjJ6rJntw926Is8Xn5LbxZ\nxB9Tk1lfR5WdyeldfdGsYlGE3syJpkxuHv8SM6BpUpc6smows2W9LbIPAtiVUrofwJqU0jwApJQO\nAlh93NIIIcZC7hv/KIALzGwFgDvMbDOA+DNa+7M6MzPTP56amkKn0ymXVAhB6Xa7mJubyzq3aFY/\npfSSmX0HwC8BmDezNSmleTNbC6B2v6Hp6emSZoQQDeh0OpWX6uzsbO25Awe+mZ0O4NWU0otmdjKA\njwL4IwBbAVwH4CYA1wK4rVTQ3Ewp8Txm57LML3feeWf/+Otf/3ql7FOf+lT/OO6h5u36WKe3C6Pd\n6cuinH5uYOXKqkPEu27iSkF/f9FG9LYtW1HIQpmbbozBEkeyTEEsKSib24n37t2QLCFqyV6Lfp4i\n2vE+PDq252Hfz5LQaX8de3655LzxzwRwsy3cwTIAt6SU/r1n83/LzK4H8BiAq4tbF0JMhIEDP6U0\nB+DCRf5+CMBloxBKCDFarMm+W0UNmKVt27bVlS16DPBEjl4FK1n95CPGoroU3XSeAwcO9I/f9773\nVcr8XnNxK2yvzjftZ5aXfdWqVZUytgIvdy+CqNqzlXS+zuiS9G45dl3sF5aTPrcP3/72t1c++3t6\n4YUXaq+L24V78yHm1WerFv33lSUMYSshmUsyXlcXobplyxaklBYNK1SsvhAtRANfiBaigS9ECxnr\n6jzm2mAuO+ZSYiujohuJuVKY68bb0t/+9rcrZR/84Af7x8x+jOG03k7zYaNA1T0UE3F6ueO8BLOd\nmX3s24uy+PZiSKu3gWOfeRs12v/MvZYLc2lFm9e3H+eL/HOIbtWnnnqqf+z3FwSq9xDnFLxdX5LY\nlLlH2bxWE/TGF6KFaOAL0UImqup72IqjqPZ4tTK6fHw90ZXiVcCoLvk2SpJDxm2e64jXeVXx0KFD\ntWUlap1XHWNCDXZ/Xr2Pfe37N6rCzLRgK+n8/bGc+xH2HWF70nm5Y/3eDHj00UdrZfErOwFgz549\n/eNnn322UrZixYpamXOTZjITYRj76umNL0QL0cAXooVo4AvRQsZq47P9wFmmGbbHWMRfx+zHxx9/\nvFLmM/AwWc4999xKmbfx77nnnkrZpZdeWiunt4+jO8i7BaMd6PuixJXp+4KtUIvuQ/859jtbAejP\njS5Cb1f7VW5A1bXIXKBsJV3Eh1XH+QYWXuvb8zZ9xNv0QPW5xDpz7Xq2anEY6I0vRAvRwBeihYxV\n1WerkZhLguVCZ26kqNL6aLJzzjmnts6IryfKecYZZ/SPt2zZUinz6mdURX17UU7mXvMw1TS6GVki\nB++yi4k/oqvRwyLwfBsx4s+3wdyhMRqQrdxj/eTlZBGNsV/8M4pyRvXe459LdKtGuT1enWcrUllE\nai564wvRQjTwhWghGvhCtJCJ7p3HMvCwpJLMtvT1sOw8bD93v+EDUHW3RZvN1/P8889XynyoaNyo\ngm2skJuElF3HQpKZG867vgDe174s2p3ero92ba69ylySLKS15Ln7e4rt5YZ/R7x9zlzK8Tl71ynb\n2GRsG2oIIX620MAXooWMVdWPKp9310QVjEUqMfca27PNw1xaUZ33RJXPy83UebY6L6qN/h5iJB1L\nXsJUeE90k3niM2KRe/4zi670SU6BqhkQ752tsmP9yZ47U+d9GzFZKnM15m7jzvosuk6ZiefbiM+Z\nbaFdh974QrQQDXwhWogGvhAtZKw2PrMtI2zDAk90l/hzSxIyvvOd7+wfR1vd21cx/HTXrl3947gK\n7UMf+lD/mG0YwuSMdp+352KZD0ku2Xvdl7EQ07iSzvdh7Bc/NxCfkf8cQ1OZnMyW9vWwOtl+gzGT\njr8uzlM03QPPf2YrCll4+zBW6mW/8c1smZnda2Zbe59Xmdl2M9tvZneY2cpBdQghlgYlqv5nANzv\nPt8AYEdKaSOAnQBuHKZgQojRkaXqm9k6AFcC+GMAn+v9+SoAl/SObwawCws/Btl4FZet7mKum6jW\nMXWQRfV5NS/mgfcmSkzIePbZZ9e259tgK+lKIuKY2sru3buOolru+zCaD769GNUX93fzeFdYfEbM\n5GOJMTwsapG5Flku+6hC+xV40fxruh8Ai65kz48lpsl141bqyzzvzwH8IQAv2ZqU0jwApJQOAlhd\n3LoQYiIMfOOb2RYA8yml75rZZnJq7c/VzMxM/3hqagqdTqdERiFEBt1uF3Nzc1nn5qj6mwD8mpld\nCeBkAG8zs1sAHDSzNSmleTNbC+CZugqmp6ezhBFCNKfT6VReqrOzs7XnDhz4KaUvAPgCAJjZJQB+\nP6X0m2b2ZwCuA3ATgGsB3DaorhKXFtvPPdeVEmF7jPtwzfXr11fK7rrrrkVlBqor96K7yxPb826y\naAP7MnZ/0Vb258a+9nZ9vI6Fu/o+Y+2xBJ7RrvayMZdWnEPw9cS97A8cONA/jqsr/XNhq0BjVp24\nN2EdzN2cu9Iy1sNWnZbUWdvWcVz7pwA+amb7AXyk91kI8SagKIAnpXQXgLt6x4cAXDYKoYQQo2Ws\nkXtRRWF5xFmZVwFLXBksQYJX66JaHnPpe3w9bLXVl770pUrZFVdc0T/+wAc+UCmLK/KawFYfxjKf\nxz+aK2zfwKZ59b0aG9V5L2d8Rr5fmCsxumOZqu9hqn1UvX1fsBV40Yxj5O4byKJVc1GsvhAtRANf\niBaigS9EC5no6jxvz7FVaGzPduYGjHjbKLbnbbGYeYWFWbIQU9/eF7/4xUoZW2WXG8oc22PJGv38\nQ8yy48uY65QR7WPvGov9yfZM9GWxX3x/xtWVDz30UP/4ve99b6XM9xN77hEvS+wzZuP75xCv83MT\n8bvq74mtzhsGeuML0UI08IVoIRN153mi2urVT5bDnLkymAuGRZpFWP7/3MSYTGVnK+Ii7B48UYX1\nLq7oCvPnxjpZmZfb7yEQy+Kz9e61aAawCENP/E589rOf7R/ffvvtlTJmHrEEF6zM18OSXTIXaKyT\nrRRsuqqvDr3xhWghGvhCtBANfCFayFht/GiLeFsl2ure1mO2UEnCQrYKLXeDhGgjenuZuXXivnos\no4q3bV9++eVKmc+Cc9FFF1XK7r//p5nRou3sw2aZ/Rg3lXjuuecWlRmo9u/q1dU8LP6efEhwvAdm\nA7MNWGJY7o4dO/rHbF/EiO+LeH8s4Subv2H2vz831uk/s70Bmfs3F73xhWghGvhCtJCxqvrMpRXV\nFZakw39mkXqxLFclYtFrsU4fIRfVOq9ux3z1hw4d6h/HBBBsRdyZZ57ZP3744Ydr5WSr12KUnZeb\nlTFXpjcJoiyx31mEIYuEZAk8ciPbSr4TTbejznXvMTdudGU2SajJ0BtfiBaigS9EC9HAF6KFjNXG\nZ5QkzfR2WrSn2P7qzGZjNqkn1un3VGP7uUd33iOPPNI/vvDCCytl3pUZbXw/bxA3xvB9Efd687Z0\ndJN5mMs1utB8e9Hm9nMasU7ffizz9cR5EeYKyyV3w5XYBptDYKHhbNMTdg/RpvefWRLSXPTGF6KF\naOAL0ULGquoPy3XCVBumSnmVrOm+ZWyvvqiCefU6uux8gk2WhIRFH7KosLPOOqtStnv37v7xqlWr\naq/bsGFDpezuu+/uH7OtotnzW7duXeWzN3NYIoz4TFh0JTMDcreYjv3pTZJYxurxssX787Kw71xJ\n5F4T9MYXooVo4AvRQjTwhWghSyZkNzfrDMBX2THYPvCsPX9uDKVkNirLQsPk9vcX3WQshNaXff/7\n36+UedcYW523d+/eSpnPrJM7DxLljqHFLJFqnVyLtVFHrJP1Gdv/j2XuYeHm/t5ZnWzVae6KQqBZ\nBp6sgW9mjwJ4EcBRAK+mlC4ys1UA/g7AegCPArg6pfRibSVCiCVD7uvyKIDNKaULUkrHFoHfAGBH\nSmkjgJ0AbhyFgEKI4ZOr6hve+CNxFYBLesc3A9iFhR+DWqJK5CPdoiuMrWLKjcBjah1z9bEtilnU\nW1TPWGQbcy369mKdrH3WZyxffV3b8dySxKbM3cXITXrKvhMlKzZzZWHfl4h/7swtx/o6wq4bZbLN\nBOBOM7vHzH6797c1KaX5XsMHAayuvVoIsaTIfeNvSikdMLMzAGw3s/1Y+DHw1P7szMzM9I+npqbQ\n6XSKBRVCcLrdLubm5rLOzRr4KaUDvf+fNbN/AXARgHkzW5NSmjeztQCeqbt+eno6SxghRHM6nU7l\npTo7O1t77sCBb2anAFiWUjpsZqcCuBzAHwHYCuA6ADcBuBbAbaWCers+2sDejmH2FHPdRJjtzDYz\nYPZ4rr0cV5qxOpiryMNkKQlJ9rJFWfx1sSw3bDW6QP2zjnIy+9jjk4cC1aSksU4Pm4tgKxMZJRu3\nsPbY/JQvK3EN15Hzxl8D4J/NLPXOn0kpbTez3QC+ZWbXA3gMwNXFrQshJsLAgZ9S+j8A5y/y90MA\nLhuFUEKI0bJkIvdiBBzLd57bBlP7Y51+r76Yy56Ruz8eW2lWYj6wOuvkijRN3BjlZAknfF+w1Yds\nJV3EtxGTdNadNwhmUjLVO9cMYCsoGez70jQJiUex+kK0EA18IVqIBr4QLWTJJNscBcyeinZSrl3P\n7MeSfdJKVrrV1dkUlk2GuRZLVtIxm5TNizC7mrkPc+1eVgez20vmU9i8AevP3O9rSfhwbX3FVwgh\n3vRo4AvRQpZMss1IbtJMVifLk16iguUmeSxZFeZdatGV6euM7fkouKbJGtiqsGgGsD0MWYRabqIR\n5nItSbbJyI28LInca+pKZatAGU3vvQ698YVoIRr4QrQQDXwhWsiSceeV2IjsOm8/ss0vStxruRtH\n5Np2g2Durtxw2xIblIXe+vstsf+ZTcrs+ro6Yj25KwMXk61JGZv3ia7M3LmBkrkrlvC1SQi23vhC\ntBANfCFayERV/WFsexyvi0k7685lqlvT/OqR3Egsdh5z9UVyk5c0pWTPNmausGQbLLmnv67ElZkr\nN0tewr4TTc0/Zq6UfAeboDe+EC1EA1+IFqKBL0QLmaiNn+ueKdk8oene4ew6ltUnd8+/WH9uyG4J\nufYkS/LI7PESuXLnPlhYNbuOJbgsceflhsKW2P/+2UbXGwtzzg3HHgZ64wvRQjTwhWghE1X1mSrM\nEjJ6WFQi5nQHAAAGUElEQVRYiUo7jASGTLaSe2D94s0Cdg9xnwIW+VVXR2yjZM82dl3TPRMYXhYW\nYRjxZfE6toqQJRrx/Vuy0tM/o1F/P/XGF6KFaOAL0UI08IVoIVk2vpmtBPDXAN4P4CiA6wE8CODv\nAKwH8CiAq1NKL7J6mm6Mwephq5hYtpwSWXJdTNEtl2uzRVeRty1j26yMycJgcx+ekmxHvk4211Ky\nQo1l0vFEm57VydyAvg+ZOzbui+j74pVXXqmU+fbZPoVxHiZ3o5Fccr/9XwbwnZTS+wCcB2AfgBsA\n7EgpbQSwE8CNxa0LISbCwIFvZisAfDil9DUASCm91nuzXwXg5t5pNwP49ZFJKYQYKjmq/s8DeM7M\nvoaFt/1uAL8HYE1KaR4AUkoHzWz1oIqGtWKsacJCz7BcU15VZO6gqJr6eqLrja3cY9t5e/WzaUJG\nloySJZxgZtWwYPsG5u4pWAIzq0466aT+cUnkXu7qymgiMBdvE1dfzhXLAVwI4C9TShcCeBkLan7s\n3eH0thBi5OS88Z8E8ERKaXfv8z9iYeDPm9malNK8ma0F8ExdBTMzM/3jqakpdDqd4xBZCLEY3W4X\nc3NzWecOHPi9gf2EmW1IKT0I4CMA9vb+XQfgJgDXAritro7p6eksYYQQzel0OpWX6uzsbO25uSG7\nnwYwY2ZvAfAIgN8CcAKAb5nZ9QAeA3D1oEpKwhBzM9SwOtnKL2avMhuKbfJw5MiRShlLKsnCSHPt\n42gH+vbYvbN7KNkHPncfwRJ7n7necjPWlKz4Y67h3I1UIv7em2YtKilrQtbATyndB+CXFym6bKjS\nCCHGgiL3hGghS3bvvFz1sGmiw5K8+rkrqkqua7raiq1CqzuvpM6S5JBNVWFG0+fZ5DxgOK7h42m/\njqZ7AWTXf9w1CCHedGjgC9FCxjbwu93uuJoayFKRZanIAUiWOn5WZRnbwM8NLBgHpbKklPr/ml5n\nZpV/R48eRbfbHUlo6yCiLECz58Pu73hkm5ubO646hsmb+XvLkKovRAvRwBeihdiwI4Le0ICZFu8I\nMSFSSovaTCMf+EKIpYdUfSFaiAa+EC1EA1+IFjLygW9mHzOzfWb2oJl9ftTthba/ambzZtZ1f1tl\nZtvNbL+Z3dHLIDwOWdaZ2U4z22tmc2b26UnJY2Ynmtn/mNmenjx/MilZeu0uM7N7zWzrJOXotf2o\nmd3X65v/nZQ8ZrbSzP7ezB7oPaNfGaYcIx34ZrYMwF8AuALAuQB+w8x+YZRtBr7Wa9szqezArwH4\nXErpXAAXA/idXl+MXZ6U0hEAl6aULgDQAfCrZrZpErL0+AyA+93nSWZwPgpgc0rpgpTSRROUZ7SZ\nrX301bD/AfgAgH9zn28A8PlRtrmIDOsBdN3nfVhIFAoAawHsG6c8To5/wUI+g4nKA+AUAP8L4Bcn\nIQuAdQDuBLAZwNZJPyMA/wfgneFvY5UHwAoA31/k70OTY9Sq/lkAnnCfn+z9bZKsTi47MICB2YGH\njZn9HIDzAfw3QrbiccnTU6/3ADgIYFdK6f4JyfLnAP4Q1WStE+mTHgnAnWZ2j5n99oTk6We27plA\nf2VmpwxTDk3ujTk7sJmdBuAfAHwmpXR4kfbHIk9K6WhaUPXXAfiwmW0etyxmtgXAfErpuwBYcP44\nn9GmtJBN+kosmGMfXqT9Ucsz8szWox74TwE4231e1/vbJJk3szUAMCg78LAxs+VYGPS3pJSOJSed\nmDwAkFJ6CcB3APzSBGTZBODXzOwRALdiYa7hFgAHJ9UnKaUDvf+fxYI5dhHG3y+LZba+cJhyjHrg\n3wPgHDNbb2ZvBfBxAFtH3GbEUH2bbMVCdmBgQHbgEfA3AO5PKX15kvKY2enHZoTN7GQAHwWwZ9yy\npJS+kFI6O6X0Hix8N3amlH4TwL+OU45jmNkpPY0MZnYqgMsBzGH8/TIP4Akz29D707HM1sOTYwyT\nJR8DsB/AQwBuGNckTa/tWQBPAzgC4HEsZAdeBWBHT6btAN4+Jlk2AXgdwHexMMju7fXNO8YtD4Cp\nXvt7ANwH4A96fx+7LE6mS/DTyb2JyIEF2/rY85k79n2d0DM6Dwsvzu8C+CcAK4cph2L1hWghmtwT\nooVo4AvRQjTwhWghGvhCtBANfCFaiAa+EC1EA1+IFvL//Yt5hyximD4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e489210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(white_betas_4d[:, :, 15, 0] - betas_4d[:, :, 15, 0])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}