SharedVarianceGLM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Does a union contrast of two colinear regressors in a GLM capture only their unique variances or also their shared variance?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
 
   ],
   "source": [
    "import numpy as np\n",
    "from nipy.modalities.fmri.glm import GeneralLinearModel\n",
    "from nipy.modalities.fmri.hrf import spm_hrf_compat as hrf\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate some data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1b8979cc048>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABlhElEQVR4nO19d7glRZn+W93nnJvm3omXYQIwMAxJyUNWckZlV2EXIwr80F3MEdxlXQOKawBdFUVwMWNCRUSQNASJQw5DmATD5HTv3HhCd/3+6P66q6urT6juc2O9z3Ofe06fPqc6VL/11fuFYpxzGBgYGBhMDlijfQAGBgYGBiMHQ/oGBgYGkwiG9A0MDAwmEQzpGxgYGEwiGNI3MDAwmETIjfYBVMOsWbP4ggULRvswDAwMDMYVHn/88S2c827VZ2Oa9BcsWIClS5eO9mEYGBgYjCswxl5N+szIOwYGBgaTCIb0DQwMDCYRDOkbGBgYTCIY0jcwMDCYRDCkb2BgYDCJYEjfwMDAYBLBkL6BgYHBJIIhfQMDgwmLFzfswGOrt432YYwpjOnkLAMDA4M0OP3q+wEAq688a5SPZOzAWPoGBgYGkwiG9A0MDAwmEQzpGxgYGEwiGNI3MDAwmEQwpG9gYGAwiWBI38DAwGASoW7SZ4zZjLEnGWO3+O9PYow9wRh7ijH2AGNsT397C2PsN4yx5YyxRxhjC4TfuMzf/hJj7LTMz8bAwMDAoCoasfQ/BmCZ8P4aAO/mnB8E4FcA/tPffiGA7ZzzPQFcBeDrAMAY2w/AeQDeAOB0AD9gjNmpjt7AwMDAoCHURfqMsfkAzgJwnbCZA+jyX08FsM5/fTaAn/qvfw/gJMYY87ffyDkvcs5XAVgO4PB0h29gkA2Gyw7Kjjvah2Fg0HTUm5F7NYDPAugUtl0E4FbG2BCAHQCO9LfPA7AGADjnFcZYL4CZ/vaHhe+/7m+LgDF2MYCLAWDXXXet9zwMDFJhn8tvw35zunDrx9482odiYNBU1LT0GWNvAbCJc/649NEnAJzJOZ8P4P8AfDuLA+KcX8s5X8w5X9zdrVzX10DCq1sHsGbb4GgfxrjHC+t3jPYhGBg0HfVY+scAeBtj7EwArQC6GGN/BbAP5/wRf5/fALjNf70WwC4AXmeM5eBJP1uF7YT5/jaDlDjuG0sAmPoiBgYGtVHT0uecX8Y5n885XwDPEXs3PH1+KmNsL3+3UxA6eW8GcL7/+hwAd3POub/9PD+6Z3cAiwA8mtmZGBhownX5aB+CgcGIQavKpq/V/z8Af2CMuQC2A7jA//h6AD9njC0HsA3eQAHO+fOMsd8CeAFABcAlnHMn7QkYGKTFYNl0Q4PJg4ZIn3O+BMAS//UfAfxRsc8wgHMTvn8FgCsaPUgDg2ZioFgZ7UMwMBgxmIzccQjH5XhhnXE6ZoV+Q/oGkwiG9Mchvnf3cpz53fvx3Nre0T6UCYH+YUP6BpMHhvTHIZ5csx0AsKlveJSPZGKA5J28zUb5SNLjvGsfwonfXDLah2EwhmGWSxyHoGATi41/khoL6PNJvy0//quCPLzSrAerAucczDwvAIylPy7hRcAa0s8KZOm3F4wNNFFhonJDGNIfh3BcQ/pZgki/rTD+LX2DEGL+hcsN6xMM6Y9DBKRv7l4mmEjyjkEIRyB6x5j6AQxtjENQX3ZNUchMEDhyc+ZxmEgQid4Y+iFMLx+HoKlq2bB+Jhgs+Rm5hhkmFCpG3lHCeK7GIWjaWnFMR06Lz/7+aSx91QuBdQwxTAi4LodlMTjC82HubQhj6Y9DkIPKMZZ+KnDO8dulr2Pl5gEARi6bCHhk5Vbs8flb8fir21ARbig39zaAIf1xCJq1lo2lnwry9WumBPDKxj786N4VTft9Aw8PLN/i/X9la0TTN/JOCCPvjENQB64Y0zQV5OURm0kMb7/mQfQNV3DBm3ZH3m6+reW4HLY1+UJ6KYzZ4Tyi6Rt5J4Sx9MchyIIxln46yD6RZob1Dfvlm0fKDzNZ1/ulgc51ubH0E2BIfxyC+q9MUtx07Lpxz0ubcMTX7oxsa+blIwt0pCKuKpM0Lp1IX7b0zaQ4hCH9cYgwekeWJ0bjaMYnvnjz8xguj5y8E5B+ZWTYZ6TaGW28tnUQd76wMXhP19mz9MNrYCz9EIb0xyGCOH2HR6x7o/HXD9UA2Uzdl+T1kbLAG5F3PvLrJ/GN219s4tE0Dyd/+15c9LOlwXtylzguN3H6CTCO3HEICtmsuG6EvAzn1w8VCTTz+pEFWhopS7+BweXF9TsCn8N4Q0ka3CKOXMfIOyoYS38cQgzZFHV9Y+lXx7X3rcArG/sAqPX7pso7vqk/YpZ+A4OLy/m4XxyeZrzVHLl/eXodHnhly6gc31iCIf1xCCInx+URojJFpZLhuBxfvfVFnP39fwTvZTRX0/f+j1RUTSPtuHz8yx9k8QdBDpIj9+WNffjIr5/Ee65/ZDQOb0zBkP44wI2PvoZP/+7p4D1NWyuOK1n64/vBbSboOlGdHRXJNZOPyQIdKdKXZY9qcDnHeI/+pfBlmu06bnRgv/jnj4/KcY1FGNIfB7j0pmfx+8dfx/JN/QDCB7rs8ojzcbxP0ZsJmeRVV6qZIa+0alPZ4fj5w6/ispueAQD84fHX8een1mbeniofIOn8HHf8yzvkKyHyl+UdgxCG9McBDlswHQBws08O1MErjht5WI2lnwyZ9FUk10yJw/ZJv+K4uPxPz+HXj64BAHzqd0/jYzc+lXl78ozi0VXbsPtlt+KpNT2xfTkf/9IgnS8Ndg43pJ8EQ/oZ4If3rsCy9Tua9vtT2woAgM39JQBAseJnd7o8Er1jOnky5Gujlnear+k3IrukgdzOA69sBgAseWlTbF/ZNzQeERhCbmgQmcAGNQzppwTnHFf+7UXc8sy6prVBSSbDZQeuy0P9Uore0SWt/mIF3/r7SxM6dV9+/lWXqhm8t6W/iFVbBsLonRErwxBtp8VfFWxIEZrp8glA+k5U3ilWXGMEJcDE6acEdaxmPswk2wyXnYgF58Xpp5d3rr7jZVz3wCrMn96Gfz1s13QHO0YhJ14pLf0mEN+RX70LFZdj1xntAEbOkStna7f6pF8sx9t3x6m8I/ooQnnH+1+suEbuTICx9FOCOlYzi5/RgDJUdtDvL+1HbWZh6RcrbuT/RETMkTtCcfrUP8LonbCNZiZqyYNLa9571FVJWOM1ekfsr6G8451I/3AFHzQRO0oY0k+JsOJl8x5gR7D0xeSSeMim3jEERaomsGUkO26rZeRu7S9mTshMEaffN1zOtA0RJYnFyZGcRPrjsVhfsMwlBEvfv4krNvcrvzMJq03HYEg/JciyaKbTqBJo+i7ueGEjdupswR6zOnxHbvpU8yB1fQKTvizdVEvOOvQrd+LDv3oi0/aD6B1XJP1K0u6pIcs7RIpykTnAuxYjce+LFQd/enJtZgPMgDDrLQYRbdWNsI6CUbQN6acEPSylSvMeGtHSX7N9EPvN7ULetmKOXH1L3/s/3p151SBzWpK8Q2T5d6Fyo3abQqOBpS/0k96h0NLP2tKWSY9IcbgSt/RHKmTzG7e9hI//5incn1EphKilH5VZxc9EtLfYmbQ9nlE36TPGbMbYk4yxW/z3jDF2BWPsZcbYMsbYR4Xt32WMLWeMPcMYO0T4jfMZY6/4f+dnfzojD2cELH3qyMNlB2WHI29byNks5sjVfXCDcrQTl/Prk3c4MJyhrNMjkDrNCPsE63T7YCn2eVrkfP1ClnfEPiTDcXlT1xIgrN46kHgMOhgoCf4tOWQz4Xp2tBhLvxFL/2MAlgnv3w9gFwD7cM73BXCjv/0MAIv8v4sBXAMAjLEZAL4A4AgAhwP4AmNsepqDHwtohqY/XHbwtVuXBdPX0NJ3UXZcFGwLOdvyHbnxY2kU1iTQ9GMLziTsNygSiePik799CisT9OFa2NpfDF6Tj6BXIHqR9IfLDl7e2Ictwnd0EDiMK/XLO54jdyTkHa9tiiRKiyHBmi9JyVlJaMuo7fGMukifMTYfwFkArhM2/xuAL3HurTPPOaesj7MB/Ix7eBjANMbYHACnAbiDc76Nc74dwB0ATs/oPEYNZFlkGb3zy0dew4/uWxkspE1tDJUdlB0XOZshZ7GYFqtL2raw8MREhUhqXJGtSQ4+kUieWtODm55Yi8/8/hmtNrf0h6ROpCtKOtsGwtfFiosP/N9j+N7dy7XaItD6u/LMkwYdUQcnjFSVTQoXbcmlU5XLjotVWwYizvZyEKcfPe8PHbcwVVsTEfVe/asBfBaAeEUXAvhXxthSxtjfGGOL/O3zAKwR9nvd35a0fVyjGZb+IFn4QjVNwJd3Kq4n71gMZSebOH1LWGJuokLUzFWhqQWfiEQtmK67bsTH1oHQaiejQCT9HsnS3z5YioTk6iCs5inLO945q5zHLh+Ze0/+hFzKheH/8vQ6nHbVfZFZUbGilnUKdnjzjt+7e0JLmPWi5tVnjL0FwCbOuRz02gJgmHO+GMCPAfwkiwNijF3sDyRLN2/enMVPNgUbdwxjwaV/xS3PrAeQbXIWLYBBVpuYZVhyXORt5jlyM7L0iSgm8gMhjslK0rfjpE/7UbG0RhGxRCtxS39tz1DwerjsYqjsxKJuGgURqnyO9F4VJjpSZRhIy0/rtO4dKqPkuMqZlGx8FYRZRd62xmVoataoZ8g9BsDbGGOr4en2JzLGfgHPUr/J3+ePAA7wX6+Fp/UT5vvbkrZHwDm/lnO+mHO+uLu7u4FTGVm8sM6rtfOLh18FkG1NFXrwifRFMu8broSOXMeNWGhG3kmGeG2KCkdiIedpvaKmf/5PHgWgb+mLl5OsXNG5e+9LoVHTO1QG542teFUNg9KMISBF2bfh95+RKFNDA09a3xF9XxxAy1LIJiEvzCpIEp3sqEn6nPPLOOfzOecLAJwH4G7O+XsA/AnACf5uxwF42X99M4D3+VE8RwLo5ZyvB3A7gFMZY9N9B+6p/rbxCZ8IwjIMGZK+/5s5xWpLxUDe8Ry5WVTZnAzyjmjJqpyZeV8GUIX6WZqWvtimSt7ZOlDC4t28WIbtA57VmrYfUZuyTESkL5LeJ37zFO7xC7CNBBmSpZ+2n9HXe4dCSz9w5Eqjl0j6lsUmdFhyvUgTv3QlgF8yxj4BoB/ARf72WwGcCWA5gEEAHwAAzvk2xtiXATzm7/clzvm2FO2PLvy+E2r6Gco7fge+56VNWDR7ChzXRSFnBXKB6MjNosqmuMTcRIX4sKuyNYnYhzIkfZWUsHrLQOT9IbtNx9JXt2PbIJF+untA91DW7kuClc05x3DZxR+fXIs/PulNtkdE0y9HV7fSBd3LnkHB0nfUz6Eo71iMjUho6lhHQ6TPOV8CYIn/ugdeRI+8DwdwScL3f4KMtP+xgrD2TpbyjvebD6/chodXeuPizI4CtlY8Yij48k5ZKsNwya+ewFELT8GMjkJD7YWafvyJeGVjHzpb89h5aqvOqYwZiNfpqTU9sFhUfiFep9jvfz9+IX6wxIuesjT1HdUYKm+b2pYHIFj6KQde+nZfzNKPzjo27hiOfq+JbFhxXPQXK0FJ8NTyjoL0ixW1pV8QLX02sWez9cJk5OqCiLIZpK8QWKe0huNzzrLQmrexcssAbnhwVWQ/nbr+YRmG+Gcf+fWT+NrfwvSMPz+1Fvtc/jd87+5XGm5nNCEOaE+t6cGinTojn9Nshyz9OcIgp0uISVJCRyGMFe8i0vcJLG2SHzUpO2xFx67jcmzwSb/TT1Zqprzzxb+8gIO+dEdg6Wcl74j+EXkRFYJs6Rt5x5C+NgKyr5EBqAOVVCTWDMnnGP7teC/++L6Xoynt09sbs/IBQd5RPBADpQrWbg+jTJ5a04Phsovn1zVv0ZhmQBzQVmzux75zJNL3B74NvR4Z7jy1LfgsKaW/FpK6xOIFM4LXXf5gTpZ+WpmQBqj1PcORpDLRKCm7bmDpd7ZmR/qOy3HtfSswVHJw9Z0vY5UvZf1FWmsirYwYOHJ9ScxiwPUPrELPYKlq9I7F2Ig4rMc6DOlrQnaMyRmQWfy2iClC+njesrCwewpO3nd2bFEMHfm5WpVN1wU29oVSAI0Lza5VvrZnKCDgLCAOaL1DZUyTB0f/ul33gDdz2rkrtPR1ST9phrD3zuGAE8g7g1k5cr3/G3YM48Rv3RtsF/uUI8g7NIPMwgC++em1+OqtL+LSm57B1Xe+ggtu8Nx3OStKM2m7Dt1Lcoq73PNZfOvvL8f6ZV6Sd0zIpiH9mhguO5EkGoKcDJJVqB2gduaJ8g5Fmuw2sz22X5rpq2ra7bgcm3YUg4eFBoas5YCy4+LBFeGs5Zgr78aJ31qS2e+L1mXfcCWQVQiys1b0YYhhnLptitize0rwuksi/bTXVb7/5MCN5Ay4Ljb0eolNDNlFbtHguM2ftVC0jpggBWRwjv73B0pOJJy2vcWuKe8YTd+Qfk2868cP46Av3RHbTtNw6kNZavo1LX2/I6tIX6dP03dUJOVwjmLFxY5hdZZwVvjm31/Cu378CJ58bTuWb+oDoG9hqyA/7F2t0RgG2Vc7U3CGa1v6CdvF+9bV6pH+tqzkHek9GSxin6o4PJi9keM6i/tJAwgZRDSQyhm4aXV18VBFS37nrlavNlUkISu8sV7IZqqmJwQM6dfAE6/1KLfLi2xkmZGrkk46JHkHQLAEnwidB4q+o15YxNu22ScJt0mW/opN/X47RSzxk5Z2mdFW7SsNQT7errY8DtxlWvBetvTFiB050aleJF2ilrzoyPU1/cwcuTwSsUKhoGLVTcflQQ0eGtCycHDSJSQLn97ns7b0hWMVSZ9KjYtF1eToHSPvGNLXRkmqS55lRq5q1abO1qgjFwBmd8XDKHWeJ/qO6hRoANq0o+jvw/3t2XrEwvLOHDt8rbZ7Sktmvy8/612tefz5kmNw4Zt2BxD6NVQYLDtaZKH6zvzpbUHSHeDN4GyLBRZ5WuOBcwT+CSCcQYj9tey4QfEzSuJKw8PFioO7lm0MmqWBhO5pPmNL34mQPsPPLjgcgPcMlh03Svqx6J1UTU8IGNLXhDwNzzIjV1VvXJR3yDEmbjv/qN0A6FkyvA5Lf1Nf0d/H2561pZ+zw7BRskqzXLc1bunL8k6c9J+4/BRccsJCcK7O4q0Fup5E8ucdtguWfPr4CAkWbAstOSu4ruWUg6nLecRo2O5X8iw74Qyg4vIgZr6UQWmEb//9ZVz406V4bPV2AGHYa2jpZ0v64tfztoVj9pwFwBswKy5HmxASKztyTRkGQ/p1Q9a7Zcve5dl1KNXqRhF5x+/I4rbdZ3UEx9EoeBUiJ6uKoj3cJmn6lrCcIA2gWWYIy5o+Rc0Q1asM/RkdBezU6c2mBjScuXT4VD8+76+DkBPkjpxtRerLO2kzcjlw8K7TgveBvFNxg8XRKw5XFp3Tvd7r/Cgr6iPk+A41fVne0WpG+H5U3rEtBosBG3YMoVhx0dGSQPpNLMPAOcfNT69r6mL3WcGQfp2QLTDVQ5OVM1dVCqAtbwfERBqp2LnJWaZj6VfT9B3J0m9W9I6YKxCuO5xdGzKhkQP1nUfsimntebzlgLnK75E8oHNv6XoSqdM55qUQRrG+fJooMLr3xy7qxjP/fSoAMf7fDYyEiusq+69uZEubP5hQ+DD9X7VlAD+8d0XsfNM7cqPyDuD1/18/ugYWA845ZH7wecsIlWFY8tJmfPTXT+Lbd7xce+dRhiH9OpFUn7zWNh2opISczSIWIwC05ESLRn/JQ/qOeglBifR59oQMhMlRnryTvaUv/xSFSi7snoKn/utUzJ8eOo0/cuKeseNKM4MiC5ssX1uyfEVLP41MSO0x5g1qna25UNN33ED2EOUdEbpkTBo6GSti/73yby/GfjftfXUlSx8IHbb7zunC3jt3xT4H4JfeaA7r9/jF39b3DtXYc/RhSL9OyMlXqmlcVhE8Kk0/Z8VJP/q5ty1N9I7Keg8dud7UXY7XzwpkBVecUN7JMqZa/q1Oaa1UMVrnU6fuHbwmXVqHqOg7oaXvbc9LWpJojabpQ/RNGlxmdBSC+P9SxUU7kb7DA0du9Hj12m3zs8WTJDB5VpH2vqpCNklCaslZkWihWJx+kzR9uubjITjIkH6dkK14FelnZemrpt452wosKjkEzvs8faeTD5/zcMHszRnLO1/72zKcd+1DwXsi/WLFDYgvywdUJm25iFpSJc00D7MrW/qWOm5dDOFM48ilwZtObUqLZ+k/t7Y3EtVScTx5pyAtW6hLxnR+qlW5gPizktrSV8g7gcWfsyLXNyfF6WdFyr1DZTy6KiwSzITos7EOQ/p1QnbcKuWdjEhKRXY5i6HFf7hUlj6Rppam76o1ffE4Qk0/+RgbwY/uXRlUEAVCQhwuO8G1zpL0xd86Ye/44jxJK/iRHJ1mBtXqy3B2gmNTtPTTnDO1RwTUUcjh/le24C3/+wBcHs44SN6ZJmUl6xIWDYyqVbkAxKSky//8PC7/03NabQHqOH2aPRX8pUQJtjCYZynvXLNkBf7lRw/h7hc3AggDAsYB5xvSrxeypq+09DPy3Ks6pm2xgDxk0gBEeUenvbDdTX3D+PWjr3mLh/vH0dmaQ3+xgoFiJTi2Zmn6TbP0/eN+24Fz8d13Hhz7PGlJRCuFBcc5B2MIBmsamHPSLEPU9MsO104gEjV9IOroF9shR+60don0Na83yXFJ2cSqZ+Xn/opzOhAnQwHp50RLXyB9SyT97KJ3HP8gbnjQOw+65jwxD3vswJB+nYjJOwpLP6uEJVXHzNsscMQVFGZp6MhNp+l/+ZZluOymZ/Hkmp6AdLs7vSSp7YOlpoVs0sMyXHaC69gM0v/cGfugszUf+5wGHZn7WQpHrss9oiFiCkIYq0TvAMBnf/8Mnn29t+H2AtL37c52yW9BpD9YcsA5MK0tWnRO93qXavghihUXp71hNu7+1HFavy8jkpzlXzsaSAs5O3J9xfvJ/OSsLLJy6Sdo6U265sbSn0CQrRWV9VKqZHPHVc+ebVmBdqp05Nr6jlz6huPyoBPf+cLGgASoXPOOoUrmIZv0AJb9a1esuIHFmOVsgsZou4ZFL4cXVltgphZczmGxkOSDkM0q0TsA8LvHX8dX/vpCw+3RwEnHLNbtB8LQysGid4/lBDVdTb+WL6tUcTFvWrsyg1wHEXknuKZhFE+k3o5wv+0Mna3kd5P9KIb0JxAasfR7h8rY2l/UbitJ0ydHnFre8Tu0RntEvBWHB864h1duDabRRPq9Q+XM5R16eMiBOVx2gmudpVOMCM1K0u79SyqXY7BTyDsu96xLIiH6bVlKIktf3KwzqNJXiOg6Eix9Kr0gf657uWuFmQ6WHRRyVuzaDmjWNJIzcsX/LfmoI1fW9IFs+hX5Keiaswx/u9kwpF8nGonTP/TLd+DQr9yp1U7S1NO2WBDloZJ3cmkcuQGRu3ht2yAAb1AjopzR4ckhO4bLARll1bmHyw5+89hruOWZ9f775mj6PLDIEiz9BL09kHc0lDvOORhCsk9qm0hfDCPVGVRDR673Xlx4BxDlHTXpN0vecVyOlpwVk850Y9ojGbk5Iv3QkSuGxIpt0j3OolsFq4AFP6ZvdI00DOnXAHWa+kI201vBSQ9e3hYduVXkHS1y8v4Pld3gQaSKhQAwvUO09OF/no3/YttACZ/7w7PB9SxWnMDqb0b0Ti15Jx7K6f3Xl3dYMJAkRQgRGU8VHKs65x46cknTVzty+315Z0pGpF9PqHIhZ8Wu/boevUVyVPJOLkhYjM4oxFlVltZ4aOkbeWfCgTqqLOc0K06/2nPXVohaNQDwrXMPxLxpbYH1rytDAMDr2wdDUne5QtMvZ15a+dZn10feD5fd2Kpk1yxZgX8s3xL7biOg30pa5JwIWbb008bpWywkpFqWfpfgYNYxHGg2Q63IpE4+IZJV5M91ybAeA0AmY8CbOepAGbJJlr4Upy8ijVQnoygVqmNBPxn7rG9IvwaIJB5esTXSuVVT2iwycpM6pOOG8d6iI/cdh87HPy49UYgnb7xN6qik5+dtL3OR5J2pbXkw5pF+sIhKRp37SWm9gmLFick7X7/tRbz7ukdStUPXNamEMj208udp4/Sjln70t887bBcAoQUukr6jMWWjQ6Rm2mV5x+8/Awnyju44Xs/CL568Ez1/XSPJdYWkLL/MOEXPyI5cEWH4rVazEVDWfKDp+9vHPuUDudq7TG6QdfCj+1aiJWfhk36KvlxPH8impr5MLrblE7BQMla2RgHRitTXggmteRsV1w2s+rxtobMlhx3DlcwtfTmLsyha+hlaTa5EiDKCxKkkTV8rTt+TFFSa/uorzwpeB5p+azaaPhkqUyR5h/pPaOlHP9fX9Oux9O3YtrJmtJvDOTpb89g2UAoMIDr3Qs6KhcQCwHuP3C1jeSdaH6pa0cKxBkP6NSCSwLINfcFr1RJ6WVj68oPXkrMwWHLgco6zD5qH6e0F5fQ1jRUj99O2vB3R9G3LK1AW0fQzIn15ij9ccQLL0XGjiUplxwXniJUPqAeBvFND05evbdrkLMsKLf2k3wgsfSFDVqvWj2R1xix9Ctn0+27c0m+evKO6Z7pGEuc8KDFBpE+H7pF+9B7TAPsTf9H7LAr5Ue0iJyB7OrbUP910GHmnBkQNmF65LsfmviIWdndE9m2Gpr+TnxhlWwx77jQFF/grPcWOM5XDMfq+rWB7mr4Q8dLVmvfkHX9nztM9PHS8KktfJBFxcDnuf+7BEV/Vi4pygwEsSd7x/suEkabKJiVn2b7lmWRJk6U/tS2lpo+oviyH9rbVCNnUJf165B0l6WtmsLs89EfISYktObuK3yadvFNxXDy8ciuAcM0LapeMk3HA+Yb0a0FFEtsHS6i4HHvv3BnZngnpSz3y6+84AJ8/cx8csfuMqt9Llzka/VJb3sbW/iI+87unAXjXYGpg6Yf7ppFfiARoaUTCcMWJ1DAaEiqOrusdDtaSbRR0rEnRO3aC7h4Mpppyi+fI9X4kicjVmn6a6J3oewKF/JIfJavoHZXF/v6jF0TbVpC+7vPiuBx5m+GLb3sD/umgeQBCsq02C0wbp3/1na/gvGsfxuOvbgss/UDulMh/LMOQfg2o5AAqPrbXbJn0ozc8Tcw8YVZnCy4+dmFibRhC+KDrR30Q2go2XB4uCm8zhhlTCtg6UIoQQxpdn6KN+qQEnaGSEyED1YIyvRrELyfRyKBLIJN+2jIMTHDkJl2vQ3adjuP26sZes6cE2/Sid7z/1GcXzIrORFsk6SoWvaNps6jknX127sQ17z4keK8iY21Hri+bnX/0Aizyn0Hqw9VIP8y50Ou3L2/05N3NfcUgZNPIOxMQqugvIv29JdKXa+/UM+2VIVvPclmAJKQJLZS/0iaVBbAthjldrVjXMxQ5vn0uvw1bNDOPCwrHHuCRvOgbUflOXtnUF9tWC67rWd1Jg6eTIP9YKQdTrwyD9yNJJLfrzHb89ILDMUVw5FYbUM/67v24QlGmQY4ZnzetDX/96JuCz2V/hVxwTb8MQ/x7na35SHsqR26tpK4kUFRUdJv3X5W4SEgbvSPOoEJL328/46TFZsKQfg2oIgFoQZFFEunLGqWOJSP3GVXJBRXSJhGJUJL+tDYUK26wEhPhtuc2NNweoJ7uA17KfslxA61WtPRp28sb+xtuz+E8Uc8HwgE7FqefQgcmcqql6QfHIJBgNefo8+t24Mf3r1K2B4ThiwDQ2RISu9iXLj1jH8zujNbC0df03Vif6WzNRa5lppa+G4/CCjX9ZEqj8SD1wuyI194Jtf1UPz0iMKRfAyLnk75Mlr64xB4Qn5LrdGqZGOon/XQyhIhWqVCXZTHMmeoRRI8krTzkO7YaRdI03HG5t4i3bxkOlUP5Z2pbmBncKFTWodwuUEXT13SQW4zV1PQJojaehaYPRPuPSMJn7T8n5vDUlT3KjhsJNwV80hfaVpGxviM3fi/F6J0kpF3ohAZTx+WxNR8CeWccuHIN6deA2D8o0mRzXxGdLblYdcTYkoopFtMm1CvvJH2/HsjSRbts6bOQ9GWs2NS41a1qUwbVoB8qidfQd5ppiM+evFPN0vd+W5ZAiCh0pA+Xe/X0a2n6BJEE02j6ooQVJf3w3GZ0RMsq13N8SSg7XEH6+Uh7mWv6MdKvbekHhlHKeAtxZbvA0pfIfyyjbkZhjNmMsScZY7dI27/LGOsX3rcwxn7DGFvOGHuEMbZA+Owyf/tLjLHTMjmDJkO0fijUbcdQGdM64jXZ5ZWzdCwZuUPWbemnKLgWi9MvKOSdqdFZDUHXaqpFoi2BpR/KO3R5tXwlbnK4JgAcutt0vP2QefifdxwQ2R6W49W7rp68U13TJ8yd5g2shZylWSI7qukDUaNBzFRtL8Q1dv2MXDe2RkFXW21LX5/04/eSDj2pBAOQnbwzJKwFHCtAOJFIH8DHACwTNzDGFgOYLu13IYDtnPM9AVwF4Ov+vvsBOA/AGwCcDuAHjDG1N28MQXwQ+n1Lv3eoHAmvAzwHktyJdcgpZulX6cQi0hR8qkfT7+5sCWQD8YHTtQ5rWVthIlH4gMk1eRpqz3eqJiFvW/j2vxwUi3hJYx1Sm9P82kXyYuwyDt1tBm796JvxwWP30MzI9f5HasiLlr7Ql1QO7TRx+jFLvyVfU9PXXX/CcXksCkt2YquQJtEOCGWzHT4PeIOzf0zjKCO3LkZhjM0HcBaA64RtNoBvAPistPvZAH7qv/49gJOY18POBnAj57zIOV8FYDmAw9MdfvMhWqSBpT9cDhJp7vn08bj948ciZ7OY801L09cm/ew0fZWlb1ssslBF0nfrhUjcKgucpDPRkVtJUcXUreHITUKa1H3S9P/54Hn477fuh38/Yc+a39lvbhdylqWV/CaXVgailr6qfIeINFU2RSPoO+cdhLaCHem7yjIMKTJykzR90t2/fPYb8MP3HBrZJ034rfd97z/llrQX7JisM/Ypv/4yDFfDI3cxXOXDAG7mnK+XrIZ5ANYAAOe8whjrBTDT3/6wsN/r/rYIGGMXA7gYAHbdddc6D695kOUd1+XYMVTB7r5FSP/ztlXXOrq1IMsI9RJVGnJSlWEQEa4qxVCCJxOQL1WXKMTBrbM1F3MQi0v7EYgkdMo6OzU0/SSkG0x5UHvn/ceoM6lVEB2/hQYGqpqafg2pMM3KWaKlf7afMNUsR64qEkvOw3jvUQti38tqERUqHdKet9FT8V4HGbkTwdJnjL0FwCbO+ePCtrkAzgXwv1kfEOf8Ws75Ys754u7u7qx/vmHID8Jg2cGO4XJsqbm8zepaXatme/5XrvrXA/HwZSfV/b1wjc70jlzZQR0s8ycsPh0er668E35PThLyjsF35AqaPln4upZ+Unp+NZChrB+n33ibdp2OX1V7ACC2KFr3OcvCHt0dOOONO1f9fqOoKOQduW1V/HyWIZunvWE2AGCnrpbE76UtrUzPGMk7rQUbLuf43dI1+Moty/zf1vrpEUU9lv4xAN7GGDsTQCuALgDPAygCWO5bFe2MseW+jr8WwC4AXmeM5QBMBbBV2E6Y728b03BcjvcfvQB7ze7E5//4LPqHK9ih0PTzvqYvkpkczVMPqEO25W3snBAxo0Kg6TfcYm0/QliBUiXvpLf0VU7FwJGrSM7SGWgclyeWYKiGVJa+m5wBXA10nF7uQP1uLzpEcaCJWP0Ww92fOj7x+zoczLkXvthWUJF+WOteNeDqFlxTRe989MRFOP+oBcGCPyqkWQVNBPn2WnNeYcLP/P6Z4LNxwPm1LX3O+WWc8/mc8wXwHLF3c86nc8535pwv8LcP+oQPADcDON9/fY6/P/e3n+dH9+wOYBGARzM+n8zB/UgBypbsGSphoOREKiIC3lS27PCIFaoXZUK6bGNsEToc9R2ABFkGkBf0zsLSF7/XqlgGUmXpE+TM53qgivioB/QV3ZDNkbT0azkza8o7GveS+nte0Si1p9LzVTPjeqG6rpbFqhI+kF7eoX5Hpalb8lZs1jke5J1mlFa+HsDPGWPLAWyDN1CAc/48Y+y3AF4AUAFwCec8/kSPMTiupx9S5MW6Hm85wS5pOpu3PEvfiZC+fkZuo1ZpqtLK0vtYpclgSToi//SWvjg4kVN891kdeMmvb0LJWVmVsHYVER/1IM2KSFR7p1FUS+aqdhw0FiY1WSsoQOde0r3I5yycc+h8vHnRrOCzYGao0PPb8raWUeQdp+4Ank7eIR8EBXS0KgazccD5jZE+53wJgCWK7VOE18Pw9H7V968AcEVDRzjKcHxnHFn6a/11PWVLn+Qd0QrV0vTJWmswbY6lXOFJhPxAEVeES9NlYOkLbe42sx3nLj4YB+8yDW/+n3sAiMlZlfh3deQdzeidNDow1d5puE0redZW7dzl0soyakXv6Jwj9fGcxfDNcw+MthdY+grSL9gpNH29ATxtaWU6137B0pfB4a0BsXLLABZ2T4l9PhZgMnJrwPW14A5fs1zvW/pTFfJOxeERK1QrOSuYojfWq2lvrYJrNTT9MHonS0du+DpnWXjbgXMxd1qYANY9xXPIqeQdOQmuHqTW9HX0bjR+H4GQnFWWfjWZKQxbVKPWoKeV/0CZzIrfVvUXQnshl2kZhnqQduUsOt6BKpa+6wK/evQ1nPSte7F09TatdpoNQ/o1QPHdFJ1A8o6cgZi3LZQcV9L0kzv1ht5hfOP2F2PWHL1vtFMHVTZ1lkt0ow9t3NKPyjviQ6wdpy88ePR7YruHLvDWD1DJO7rrx+pE76SL09e19JMLtFU79VoGQy2pKc0sUXVt7SqWfmveTuHITSfv6OruxYD0vT6ptvSB59b2AtArDDgSmHSkP1CsYMGlf8WfnqwdOMQ5DxJsKKxwfa8n78hhhnnf0q9X0//4b57E9+9ZgWf8DkLtbdzhFXNrtFOn0/R5JEwzUdP3/xeykHeE78kLWe+zc2fQxrDKkavpINch4LC8RePf1db0U1r6jUqDBJ3ZDB2P6jyD/qLU9OMZ7HW3qSnv0DOiu9YRWfo0WKkGM8550HeLinW0xwImHenTqkE/f/jVmvvSM2dbLFhabrNfYVMOMyRNX+zI1eqFk7Ug9t0f3rsSl/zqCW97g506beZoq2C1JFn6RIAiSWtFtciVRAWWeuw/TsafLjkmaFMdvaOn6eslZ3n/R1LTp+usmtFUG2RVpZUbgc69rBZ4QDKhKnqnvZDTXhidc02pLmXtHVmOkvNZ6LdphbKipnzVbEy6hdGf9S3r/edNrblvuJi2Z60UclZQVlkm/ZxtYaDkRC39KjddVcr37hc3Bq91o3d0NX3xwUyK06f/EXlHM1FKhDiIdAtrAgMJ0Tua2nMaSWAkQzarWfrVrjd9omMF1/rtxO9UCROlbarErNa8viNXfwBPF70jk7jK0ndcHmzX9Vk0G5PO0n9unUf6M2vE9AJxvbKzJRd47uX6NAWb4ek1PTjnhw8G26pplir9VRwwGtWfA0vf5bjtuQ047ar76pZeeA1Ln44lcOiK8o7GAxRKAsn7BJZ+Rpp++to7DX/VzxzVj9NXyVjV5R29HI96fjsJqiJvBMYY8jZTat9tBTvVwuha2dUpDCMg/jyrZjCOy4OZq0qarIVfPvIqVm5uri9g0pE+Rd/U08GJmMnCFZeza5cyEOlGb+kPV5aqx9IXHa8iseg7coHP/O5pvLSxDwOKcEcVXF5D02ck78Q/57xxxxhxdj5wWMb3oTYGivFz0NL0NfX1NKWVqfZOo6hWf7+aNR4ScONtit9v6Dtu9QE8Z1lKS789jSNX1z/jf0d7AfiYvBM/r4rLAy1f1XergXOO//jjczj7e//QOr56MelIn+53PVNZWYIh520hZ8WsRlWWaLVOTYOOSGAisTTaqUXtudEu7fJo2d2Ypm9HLX359xt9iOjcSdZRDcDU1o7hSkROalVkQcrgnOPLt7yAZ17vCbaVKy4Kda5NoDoOHemD83SWvuq6VDv3wJGraemnkXeSZlE5iwUat4hUcfq68k4Qp58N6SdZ+iQDUY2eWugvVnD61fcF0nNfsYIv/uV5XP/AKq3jrIVJR/pEUHVZ+v49pg5GpK+qFbPBXzdXRLUHlB4wcR/xmBqVIsSysQ1b3pxHXH+Jmn5CckujsgDdAyrgpjpesWTAop3CJJfWvF2T9IsVF9c/sArn/PChYFvJcasupZeEtFU2dSJpclVCNutz5OpBb50C738SCedslqjpu1y/zXQhmw1/FYBC3kmw9EnWqXdZz0dXbcWLG/rw9dteDLYteWkznl7To3egNTB5Sb8OIyPIjvX7F2XhyssJAl7cfVJb1X5b3EecLOhGmnDB0m/EchOfoaTonaSSBI1K7HRcNLioLC/x/PcUSL8tb9fU9OnnxOMsVhylZVYLaTOds9b0qx1HkJyla+mniNNPanLO1LbYWtJAaDjprS6nG7Lp/dct2Cd+jzG1g9p1OYbLvqVfJ+nTfRajrsqOW/eqeY1i0pJ+PR1cnrpOb/dIX3biAlEtX25L+dt+Xxdrw4vHpJdxyCK/UW/flslJ1vTpM+qDb5w3Fd2dLThpn50AaFj6JO8E2nV8H3HgEauNthXsmpp+uIpRuK1UcZUPaS2ksQ514/SrFVyrWoahDgd5NWiRfo1kwpv+/Wh85MT44jGkh+sVz0tbMbXx85QHp7xlKZ3JFZdj2Nf0l766vS5rPSyyGN1Wq2yGLiYf6Sss7CQEHZpI34/4kZ24AHC4n0Eaaauapa+SdyLROzUPLwbP0g8Jql6LRowyeeO8rtg6o0Gcvr9PZ2sOj/3HyThq4cyG2gnbi8o7qodQ7PBUkgHwLP1a8o5qYC9VdOUd/zdHofaOihCrEVaS1HLW/nNw8K7TlN+ZJ5S+0JHYxVwWFVrztnLdWprl6dZRShO9oyPVbRv0jDoqvJi3mXLgcVyOYjm8kD+6b0XN307y3VRb7zcNJl2cfijv1OHIlaJ3pvtrnaqMjJ984DBsHygFBcNqtRE4cl21pa9jyXiWfhgRVK9F43IOMOCFL52GnGVho+SfoOeLHrQgmkfTyRk6cq3E74sP9cwpYXhtW95Wxu6L8g39nnj6+qSvbx1SNnejqBa9U42YVQujA8D3331I4nf+dMkxeHljH9593SPpyjA0eJrVchFqt6kpf6ZYEGfFJi+Mcr+5XXhk1TbYFlMaZhXXRbHiYI/uDqzcPFDXcqeq+1xxXGPpZ4VG5J0wOSsq76jib6e05LDLjPbItrocuYJUIe6uIwuQph/+Xn2dm/vfbS/klJFJLCB5vx0ifyHK5JZn1uHBFVvqao+u6/F7eSujnaZYyUns8DM7BEtfEfXxt2fXY+//vA0v+2WZVZZTWkeunryT1tKPNlqsOLjxsdeqtOf9b6TrdHe24PDdvVlqmuidRvtrtfpCtaCd6ZyQaDdcdiKJkSos90n/jX5SZ9LA4/iO3Nmdrdhzpyl1hRer9qkI8f5ZY9KSfj2dLaxlErX0602vruZwrCiOI2Lpa05fxd+oPzmruqYvH1MYtx9a+h/+1ZN4148fqas9uiz7zOnC6ivPwmEKaUw8nlmCvNOat2Pn9bfnNgAAXli3IzgeGcWyq8ygrIUgFFY7skXD0qdQVqnN7971Cv7vH6sTv6ebnGUnkGE9kCPc6kUaSz/r9Y6/eusyXHDDUjxVRX9fvrkfU9vygVO64rrKZ5RCNlvzFnIWqysXQSnvONw4crNCaOnXvy/N0EjTFzW76t9P/oxIRCwTHNH0Ne43A4Vseu/rfYY9h2P4PklLDCx+/2NbUyMNZLMqvU98oOToHZkoSCLLVYn7L45CyKauRZok72xVBAtE2/P+N9qkOHg3CreOe6lCGKGktwqanqbv/ZflnbXbvYRNqqulwqrNA9ijuyOIOkoaeFzuZZG35GzkbSvx/B5cvgX3vLgp+C0ZxpGbIYLIjgY0fVneqVY9ryVSa752cpa4jxshfT1LhvMweapuR65k6SfNMmgfCi2jB103Tr/aOYardDEUclZQK6ktb8cepFIl6iOQz5tzjlLFRYuGYyx9aWWd6B31edRbGlk3TFS3vlA9xyaj2upgVdsL+k5DX/O/o/bPUDZ6tbIJfcUyprcXgn0rVWYbA8UKWvOWV3k34fzedd0j+MANjwFQP6dl12j6maGx5CyyYqLyznAVS79FqA1fVdP32y8naPq6td8j8k69mj6PL6CtAs02Q4ehnoVYK4sTALpa87j0jH2Chbx/dsHh+NkFh6O9JdnSpwxfeayla6xj6TPGwJie88+zSBv+WqL0IY9Z8nVPk5FrMd0EtMbanOHPlqsloIl4besgvn3Hy8H1l0ujNIKk0sr0zFYjfZIH23zS5zy5/27pLwVRS/XkIagMCs6TZ9xpMWmjd+ohKrlD02pZbz1wTuJ3WvO2VzrAtqqH17nR4wGkjFzNNHMumPr1EpUsQ4ha4u0fPzby+94x+8dYJcqkGgLZrMY5fui4hcHr6R0FHLtXN+5/ZXOsPXKEBUQinTfpqjqkD5CvpPHvcc71HPIJpZXlkskO57CEbbUSpaq2yVjmVTZVuPOTx6FnsBQsMFLL0fn+Gx7Fys0D+NfDdsG8aW3CcqL60Tvyc1lPKeSS45O+kKOTxMm0b95mVQ1EQpCcJd04Hb9ePTCWfh37EunnbAtPXH4Krvjn/RO/Q6nZFqveocOQzXCfNLV36DjFkM165VKXR3VgMWpg7507I78PxGdA2vKOxknm7HjtHYrmSRqEyNrSycgFvHuhW2s+y9LK8uW67v5VAIBP/OYpLLj0r6lKK9sWS7VcYr1GyoyOAvbonlI1LFUEZbXSLC7MOm74UBNXzqJEsVqWfkGw9L1jiB8EJQC25m3krGRNXwT1X/nX5MWFssKkJf0/P7UOl930TLB920AJX/jzcxG9XiVDzOgoVI29Je25RRFlojqOaEZu+LmWvANJ3qk3egdRfZKa/vAJ0UxKsaib975xeYdzjm/f8TIAvdlMzmKxB4kemqRwXLqn6Sx9XU2/8faSBi+ZZL5+24vYuGMYf/RXgQuidzQihmymq+mrj61me4GmX50UA+ubU3vp5R25u7bWYelTHohYjVZ1DJQEV3G5v4RqtLHewXKsXHiSBGSbkM1sID5Iv350TZCI9PDKrfjpQ6/i5Q39sX0bkda+ee6B+MWFR2CX6W0NzSbk17plGDhCa6ju5Cw3+tAyxrD6yrPw6dP2juwnF1zTsfRXbx3E3X7Ugs701ba82Yw40FSk2Zs8w6KHSqcMAxA6yBtF2uQs+TzU0SJxn5AOV3h+i8a/xxuUdwj1hmzSvZMz6dMkZ8ltVluak1CsqOSd+DF86lTvmelszflLqEYJ/cAv/R1v+d/7I9tIfqy2uFCWmHykL11YIiCyFocFS1+O3qkH7YUc3rRoVt3TZbEDps3I1U3OqtcipYFBvi61zrNvuIwt/V44nOgk1pnNBBE6wrmRo9Zx1JZ+QPralr5+OKOu1ALEr6vqcom7pHHk6so7uhp7taJyIogQK8G9hVZ7QLh+g0zEdN60hKkKxYqLlnxU3lFd58N3n4HbP34sPnjsQqUUCQArNg9E3lP/lJMOjaafATjnsY5NNdepU4mjPdcgfYJtscQOLXa6iLwjHBvTuDOeM67xkE0veqf2frYk5wSWfw3Z8k1fvweLv3Kn4njrOrzoMSjIgq5hUl2lYqDpj7QjV6/v0MBGq7Q9uGILKo6rJDpRGkhTWlk/ZNP73+h51lt7R561pgnZDMJEpeeSiLm/qK6KWXFcfxlEWyJ9dTt779yJtoKNvM3qWjOABjZZ5jEhmylQdlwMFCvKB5divElbFL3tcpRKI6j2EIlZeon19LXD7rgQ3lbf92RNv9rve78blb1qkYVYV1w8X11N3/ud+DV0JZmHkDZ6Rw6FrRe6mn5r3saB86firhc34uGVW/GuHz+C792zXDkwi0RB10QnYogxht8+tgZX+f6WeqFbe6daUTkVgntcR7hvEmigKUttktEgW/rUn6j/1CPvRNqzrPpI37+HcvauKcOQAhf+dCne8IXb1UkQDk2t4pZ+Gv0wZ1mJVoyY0RupvZOynj6TLNL65Z362gtXHooeYyOygBiGqBW9o5A+ypKlL0sxdL215R0KhW0QuslZAPDWA+fiubU78NiqbQCA1VsGlL8lkgqRh5akxLy8ku/c9UpD39OdDdcbvUOQB3SdgS3JVxJa+pXg/68eeQ17fP5WvL59MNJ/xNlirf6bs5Nn+wRKHASAciW6b7PKMEyKOP37Xt4MAIGuLIIsDRrtRdJPY1VYVTTSYiVupXLOIyO9tjNOWMywkYzcep6hpJDNeheLACRLX2cGRdaa/zB96rdP43U/jT6prpJoqenAShHZorugyT47dwEAXt02CMALVVXNjMS+RJEiaRZuaRTUZRttM6moXGI7nJ4Taq+h5qJtShY1cQCR/hu/cHvw2csb+7DfnKkAvJBf8X7WOue8XdvSr7jhc28s/QxB5RNUCxoQeVDnG1ZopLr1U4h8bnridSy49K/oHfTIUVywnCzfch2RGrUgl2Gou8pmneR0tF8//4D53kNAJETp5PVAtHyyKDv8hydeDz5LysEIo3d04/QZfvHwa3hpQ19D39OtvQMA0/w+KzrAVb8lkkrZP0+tHA9NJgiekQa/X29GLqEiOel1S4/nbRapdyUeg2ohc87DkF/ZaKh1DJ6mX/38yo4byjuSpm8cuSlACUZPKkk/SrpF39L/8i0v4H9uewmAngwhlmH4/j3LAQDrd3gWab+wYDK1K4/y6TR97329ixLVS04n7TsbT//XqThiD4/8k66L43KlHPLJ3z6FrQNh0TCdTq3S9MV2gfh5B8lZijVN6wEd5m+Xrmnoe7ohm0BI+lQELG9bSt1GJIowyUc/hr1R6Mo7dO/rXRydyD6N5ApAmTBFz6BqnQaP9KP95+cXHo57P3N8zYHOi94J21I9E+UKT4zeaVbI5qSQd8jj/urWgdhnZEGI8k7vYDmyEr0OAdtCWjs5h+n+ixYFdWJ5lNfp07KmX68k0Qg5TfXJCEgm7YWfvxXnHDof3zz3wMj2m55YizW+XAFoOnITIjAAgRhijlw/OStlLROZFHoHy2gr2Im+Al1HLhDWeSLSz9lqv0LE0ifS19T0daAfvdOYpu8EMqjfnuaFzSmsb/ptVSFFjtAnRBndb17krQPx4oYdVdvKW15bVI5Dda4lx02Ud4ylnwJkcatG8sDSFwj6b8+tj+yjpz2zYJQf8mcPg76s0yeQPu2zakt/5Pt6ERhysk4Dmn7DrVV/0H//+OvK7aKGqyMpUJaiSguWo3gIaR25Yf+JTv8P/NLfcfHPlwIAnlrTg1O+fS8GihWs7RnCX59Z7y/grffgtvshfyTvWExdm100Fui1Xl0aXU0/bfROY6Sv2x4hL1nfQMgBqoxcl/PAaJDlnXo0ffGY1X3WDdqNrcPbpIJrdf8qY8xmjD3JGLvFf/9LxthLjLHnGGM/YYzl/e2MMfZdxthyxtgzjLFDhN84nzH2iv93fvanowbdVDn9WfxMtPSfW9cb2UdcxKNeiJo+OYcH/PZFS58s1vte3qLdkQmW58kNoEooWt87FOtc9Wr6MlSDodjmpr7h2OeitaMzmBbsZFnAlYiBkDZkk9oaLDl4Yd0OrN4SzhiXvOQFCXz11mV4ZVM/nnm9F+dc8yAu+dUT2ot9AN79mNZeCCzp6x9YhR/duzK2nzgQkCNXK06/xnFu7S/igVfiK6OFmn6j0TuNafoycepWoMxZLBYlIz+nIjhPNhpq3ducEHRQrDj4+UOvxvYpV7ggMY89S/9jAJYJ738JYB8A+wNoA3CRv/0MAIv8v4sBXAMAjLEZAL4A4AgAhwP4AmNsepqDrxdErAMK0q9InWm44gRLoxFmCWu01guxlgnNMIZKYUgYAHQUwvo8D63cGtTt0YUlWfrxJCUHR33tbnz2909Htmsv4K0KIRSsqFc29sc+r2f5uGqgKXax4sbkjiRH7raBEhjzlrTUAf3uUMnBmd+9H8d/c0lsH7Ecwfpeb7AbrrhaUgthWlu+5j5KS19zlkhQFQk76mt34z3Xx1dG080CrsfSF0kwJH1vm27iUt62YnH6gT9P0adESzzmyK0Vp08GiuviO3e+gituXRbbpyQ4cmXZaVTLMDDG5gM4C8B1tI1zfiv3AeBRAPP9j84G8DP/o4cBTGOMzQFwGoA7OOfbOOfbAdwB4PQMzyWGTX3D2NJfDKSboVLcOx9E7wRx+i6WbxqI6L86VoVtsaAsAIGSP/p8R+609kJwbFv6i5gvrbHbKOTCYPLzROf69xei64HqOhxV8sxwSbQ81dNlgmrmVQtU8Gqo5MSm4//9lxdw0U+XRmYbnHM8v24H9pjVESmW1QiImAYU/YdATYozJs/S12oSQKjrV0MpK01fOFCZfDb3FcP6MFKnSrswulPFkSvKsXJdJV3SV8XOUwQd5/FzL1VcIXon2n9qHQId40sb+pTh4kA0ekfGaBdcuxrAZwHEjs6Xdd4L4DZ/0zwAYpjD6/62pO3y713MGFvKGFu6efPmOg9PjcOvuAuLv3JnYLlU0/TJgtjUV8SW/mKwWLQucv6qOeJDQprwQLGCvM3QXrCDDjdUctCuSUoixGdyxeb+aLKZE3WGhd9JVyNGxJDQnsqqF603nSzXoAxuxVHezzuXbYxmObscz6/txRvm6s+iKgp5MIn8ZCtZV94BwgieaiirondS3kt5sBZnvrJlrhtNY9excpYqUZL219W7cxaLafriMcjOXI/01dFfNTV9f2Zw7g8fSqzrUxYcuapjbQZqXjnG2FsAbOKcP56wyw8A3Mc5vz/h84bAOb+Wc76Yc764u7s7i58MyEdFEhUnOrUivfbAXdJKLZ7VvWM4TFwieam/WEFHS85fWYfjnhc3YaBYCdbfTNOmyKPfuP0lfPp3oZQjT2sJHJqavuI7orNTpbvTw/vPB8/DIbs2ru4FS9uVnJhjlSAOJpv7i1jXO4w3zutquK3w97z/oqUvkxW9lWcfuo5coD7Sj2j6KeQd8Tuy5SmSr0yYutE09VTZFNeslUuR6+rdXsKUJOEI7+VFT0qOK0TvqOWdpEPJC5Z6UtnmsuMmSp6juTD6MQDexhhbDeBGACcyxn4BAIyxLwDoBvBJYf+1AHYR3s/3tyVtbzqI7IgkxGeiLIVsEknPneateq/buXJ+nP72wZD0BwXSn9KSQ85iuHPZRnzghsewY7iCtkK6CFrLiscCP7hia/A6qXPpavqqB1209FWkTw/5JSfsqUWIFH47XHECy/tbUmjox258KnhNVurC7ilIC9HS7xmKLlJO1112BjZd3olk5FKcfuOwIs9ENdLPSt6p7cgVHZ+0Hz2vuoSYU5Q7Fgcy+f6VKi6KQUZ31CgjbkiadYjHmGTNlyo8sf7QqGXkcs4v45zP55wvAHAegLs55+9hjF0ET6d/J+dcPOqbAbzPj+I5EkAv53w9gNsBnMoYm+47cE/1tzUdcllWUa8P5R3vQyopsHNXKwDg3EPnQwe2ZcFxeKREwSDV9hj2SV/quFlY+rJkIg4CSYkwupq+ytIXiZFilCNtBdNzvYc2XMTaDQbRKa3JgyWVaNh5aqtWeyLEKfrG3qhGS9d9uCKTvj7rT63L0helCZJ30kViyf0kSbKrOG5QTE+39k41x/4Dy7dgYXcHAEE+C9ZD1pV34uWOxfey36ZYcYOETTl6h0456VjE7SVFDgDgXeukga9Z8k4a0/KHAF4F8JDfyW7inH8JwK0AzgSwHMAggA8AAOd8G2PsywAoZ/9LnPNtKdqvG/LIXrCt4AGpBBaE76jyr39HSw5L//PkuiIoVLAtz/mkkncGSh7py89JWtKXk7Nk0DmK9XmAZmr68Q5N11t3BkWa/lAp1PQ7q0TlUDIYDeJpIJ4bLb5DoC4mywMjaem/srEPU1pywVrOjaCavDOUIO9cetOzQT5Go/3HsrwF5+V1gEX0DZex28ypWLF5IBZlpx+9Ey93LA48YrY8EL0WSWUYkvqyaNiIv/O+o3bDzI4WXHXny568k0T6YyEjl3O+BMAS/7Xyu340zyUJn/0EwE8aOsIMINfaKOQswDfUwjj9+D468fkE27codgyFnYjkpf7hCqZ3FGJTyfaU8g5DdedoUufSrfuukncGS1F5Jyn7UXfq2irKO2XverYVbD9cNb7/69uHULAtzOhoPOy2GjYKOQgLLv0r9vFLfcj3NI2mP70OS/+J17YHr1dvHcTRC2fqFQgUSV+29MUoGuEiiwl42ktfJvZJjoGSg67WfKTdMHpH39KPkb7LUchZKFXcSOIk4F2LsuMiZ7EY6dO9TTp18RjFa/qmPWdhlxntIeknafqm4Jo+ZEtfnHYFZRjc+GwgDXKWV4aBpr/T2/MBIQ6XXbTm7Ni0ML28E98mdqfA0s8qekfxpeFyVN6RCSRMrtEjw5acBcbIkeu11V7IIWmoW7N9ELOntqQiXxU27ojKOzRzjFv6KeSdttoD1aOropPlA3eZptWWOFDc9MRa/MsPHwreD9eIyAKyX62rWPFmiV1tniEkO3J1+08+Z8WMwIrjosN/9mRLv1h2sX2whOkdhVgfqlV3SDxGMSEsZ7Pg2S85VTT9sWDpj1fInvN8Lq5fyhapbvYmwfKtGCL9nae2BZZ+2XWRs1nMImtrgqYvItmRq1mkS3GJ4pZ+9No7KafnjHkW13DFFUjfTlzfdc22Iewxq0OrrWrYKsVdk+9GrtSYSt7paFym2Xt2p1Zb4qyN6k6VHRd525LknexI37O643HxHDxIYCSpKnDkpvQJ5a24I9dxOTpactg+WA7aDY7HcbC1v4QZCqmNznmnzhZsGyjFPhcNR9H4yVlW8Fm5Us3SH/2M3HELOVRTvBkVvyJkMyx9wIvyyFkMs6YUAkeg43LkLBabvmXiyE2WSINzlLuYbmGwWo7cissVU+l0IXeAF8Fz7X0r8etHX/PeV7luW/qL2KlLX6ZLwtb+6ENO1UO3Sg+/bk0bAJhWh6UvQzfrWMWhZPUOCQl3SRq8zvKeOZvFfu/oK+/C/l/4ezB4EumHVTap/2jKO4rkrLLrBtdNpel7ln58AJ7eUcA3zjkAP73g8IS2REeuSPosMDxLVTT9sVCGYVxDvH6yrOKRU7aWPt2wnoEyutry6CjkAku/4nDYlhUbyVNr+iyu6Ytvg3OU+hjn+uWjZYhWYanixuqc0AOXppgU6fpPvtYDAOiocd10ibAakjIstw1Et6cqw1CHpi9Dd7aostQpc3xIkuxUyErT39JfQslxY5Z+GHCRNiM3XobBcXhgcMU0/YqLbQOlRJ/QuYt3weyEIAHxWRTVhpxtBf2/qKj3Qxj1gmvjHR3Cg0+ETo4Zz5kiWfoZkf62wRK6WnNob7GDGYfjci8jtyX6gKaP3onr9SKSppG6mr4yTr8UjfRI0vTTWDHiOe41e0pNotNdMasaZIs+2C7NANJo+q356ELcIm74wGHK7dqkr7gfm/uL4JwrM2Nj389Y06dZcRdZ+pIjV5cQPXlHtvR5wA+yPFdyXGwfLNcVSSUjagCFr20r1PQHq5C+sfRTQrT26GKSxVh2eMziSEv6ZIlsHyihqy2P9kJI+hWXw7ZYLNSQyunqwmIMj65OjoJNzMjV1PTVGblRqzApNyCNXik+mDtPbau5f0sG5S1kbOlTW/ryDCDtc7sgwR9x/N47KbfrGg6qe/mOax7ET/6xWsq9SJB3NM5TFTNPoHssR+84KeXBnB1fRMVxeaK8M1x20TNYwswa0V8zOwp431G7RbaJxoY4AORtFsjHw1XqTxlNPyVES586OFlRFUVoYVpNnzrl9sESulo9eYc6suN6IWByUlFbwcYjnz8ZD3zuBK02VdaWmBxFFo44fQa8uP00a46KkDNyVSRhWyxVNE2/kEDzn2ftW3P/rC39lpwVkwEIq7cORt6nsfQB4LcfPBLzp9ce2AhJM4NaSJLIb312fSThbKQsfeqfXQnyjrYjV7FcYtlxA36QHbmb+4pwuaffV8Pjl5+CL539xsi2o/aYiQ+fsGfkuAGy9L3jV5WGIWQdcUaYsKQvZ4KKpE8dlBJ9PvW7p2MZc9mRfhmdrTm0F3JBGJpn6VuY0hLVbNsLOczoKGD+dL1qm7X6iEjA4uLP2lU2lY7caO2dJNJPA7q1f/nwm7BXHdEqutU1k1DLMSz2nbQPbmdrvqFZp668kyQLdrbmYs55IB4Grbv0ZdLMIXTkes+tnJGrX09fbel3JGj6G/wkPJ08D8YY3n5IrKYk8rYF22Io5Kwgum8kMWFJX542ThH0c3oOiQyWvLQZKzaHC2PkLJYq6gIIH4L+4QraCjY6/PYHSxVUHI6czdApW/pNkCFEJJdhaI6lX3E4SpU4m2Q1bZ1Z5zoHWVv6szurZ/fu0d0RXM9MTlW4hFNacvjDvx0NAPjkKXvFdk0bDCBjSksOw2UnkD+CtSmkqpG6/aempR/IO9729I7cuKZfcThaCzYYA/qHoyRMRd90NH1A7Xugme60tnyib6iZmLikL91YMcIjtPSjJEukm1bPB0IJqeS4aC/YgQU2WHLgkKYvkP4Ru89IbQGrVv6puJ4jrlRx8ZokPRA4h5YoqzrcIT85qS1vY0t/ET2D8U493kk/ydKnxXb23rkzrPOU9YDT1YJDd5sOAPjoSYtw8K7TIp/rGg7Jln4eQyLpu9HChASdGY23cLi64S2+Q7wriNN3/f9pyzDEo3cqrou8ZaElZ8XkHYJuRrfq/lP1zalt+cAHJF++jpRBHdUwIZOzXt7Yhw/9PFoJWnTkkn7ZKtXHntKaw1DZyYb0hU7ZlreDQWegWEGFNH3/mFpyFn7zwaNStzmkIP3BkoN9Lr8N5xw6P7Zu7bqeIdz+/AYAepaa6kEfFhKm7ly2CXcu2xTbR3dqLkOueph2v1pozVsYLruJVh8R1fF7d+PPT60DACzeLd26DEA0wlamSNmS1DUcqiX1DZUdz/+0IzSm+obV5NgIbAu444WNWL1lIOaw3tJfRGs+TGJ6cUOfP0tO6ciVondcl8Pl3gygJWcn1r3XJX2lpe/r+dPa81jX48lHLTkrks39tXccoNVePZiQln6p4mKlsI4poNb0CxIZUDRNWj0fiKZQt+XtIKpioOh4ncyy0OlPXZNkl0ZRbSWq25/bENt2wQ2P4Yt/eQGAXvSOCoNCPZwkpLX0f3rB4bj0jH3q3l9e/EIXJJ10JMT9n33QXADAWw+YG2zba3b6ks4iIcvcnEVfBeKDCWG47GCo5Aaz0qAa7XB6LXqVL6l+/o/Pxj7bNlBCRyEXyKy3PLMeF/10Kcp+uLOur4RmFzzwEYQzh7xtJa7RoC/vxI8z75/T1LZCYOnTfexqzeHJy0/B2w6cG/teVpiQpK+azouVB6nDyA8MPcxZWPqik7O1YAe/3ec/LDk7tPTrXBu6JuS6LyJU4Zpre4aC17rE+MTlp0Te08BTLWEqLekft1c3PnTcwrr3z0reoYE7KSzy2/9yEJZ96XTkbAufOW1v/MeZ+2YegSFb5FnVZ0ky9IdKDoZKlcBAIXmlZzA96VPV2Z0643JZz2Apdt8eXLE1kEZ1QYRLZP/Acm91vpxtoWAzZTRNwba0HeQqS59mulPb8sKqXN7v2xarGSmUFhNS3lGRNi2KAoQp53KnImdrFqQvllhoy4eaPllIsqafBeSl3kSoMilFH0Cr5jnL096hkgOLVR9E7CYVkkpCVvIOkX1Shq9tseA+X+KH6mUBkZBl0s8ua1PN+juGyxgoOYG/gjR9la9GFzspMlp7hsrK/Iqy40ZWpGoURLgVhyNvAxfcsBSAJ2/mc5ZyhaukBVDqQZIjF4hmXBMXNSshS8SEtPRVpD1LcPqRFS5Pveh5ymLKnKTpU4iWqOlnhWryjipKQhwI0hZ7IwyUnCAkLQlpHtpa+JmiDkpaS5/6Cck7coTMvGlt2rkVjUKesKVJ5ov8rtQ9Tt53No7YfQbW93qac7dfZpz08J4MQg0pCklViqB3sKy8bxT5pgu6XmXXjYR1r9462JSyB7YVL6xIxyCu1UHnmjavox5MSNJXWXZ5oQORTihfYHJOZe7ILYSaPtXXtxXJWWmhcuTWi6xi2fuLFRRsq+qKSM20ZlT1atJq+o/9x8l49PMnBfewQyqfsdfsKdq5FfVAXvRGRLPqsxRy3qyFJMBuX4IJLf30pH/obtOxU2eL0pLuGSorn8OK66YKBBBX7NohOKMvOGZB066lPDCrLH3yLzYrC1fEhCT9amFSQBipEqy85JNvX9HryJk4cmVLX9b0hfobh0hhd7pI4xvIivQdlyPvL0iRhGaSvuo80so709oL2KmrNdHS1634qINmyTtyMmPOstCWt4P7GJI+afrZyDuFnBUsPC7CcXmypZ+i/4Tyjov1vd6A9r13HYw9d+pEQSLnORkssQnE7xFxUVebSPq+pW9IXw8q0hanhEQ6jt/R3zh3KoARsvSHydL32rjvMyfgZxcekbq9tMgyMSxvs6o6aLMsKkBtKcmhuboILH1JCmu2dSbysexwzcovJNsLOZtF+gSRfpaOXMCTNYoJfUU1WO8YLqeSd4gbihU3kK7m+PWbxH558r474cfvW4xp7Xn82/H1Bw2oIPd3IvYuhbwzEpr+hHTkKsOkhG0UUXHyvjth/3lTcf5RC/Ce6x/BJ05ZhAtuWJoJ6YsaeVveRkvOgsXCxTao4+46s3myQCPIslRB3o5a+jM6CpFFJprZsVUDSlaOXJJ12iVfTLNWOCJUc+R+6tS90VHI4Xv3LMfhC/RzAuRZYsG2IpFHFGFTDjT9Eqa357E9Jfm35Gylpe99Fr+XG3qHU/mEKKFvS38R63uI9D2LXuw717znUORtC0/916nabRHksg8EcVClcx2JEIcJaemrwuRylhXEdpNlZlsWPnrSIkxtz+MvH3kTTtxndqQCXhp0C+vrtuZtMMbQUcgF0Tsjod01gmwtfSuw9OdMbcWP3nto5PNmnruKgLMK2WzL55RrpY7kvZTJeUpLDp8+bW8sv+IM3Hjxkdq/G5N3JEuf1ot2BE0/zRrShELOCqLOXOnkVMbXpr5iKqOBrPr1vcPY0DsEi4UDmuj3y/Ke7khIZIuSfnNLsIiYkKSvQslx8aHjFmL1lWcFo6nqvrbm7EwsfTFVn6z+9hY74sgdS2grZJGb4P3P22EhrR+8+5Bg0XBCMy1j1WLSWT1QZ+y/Mz543B6x2URWGcZJ4JHkrISl9WwrUz04Z1lBn+go2JEy5IBH+vWWwaiGFsH/I5dkUA3W/cOVVNd77jTPql/XM4QtAyVMby8Ev0eafiFnNa3CpQhxdp1VAmE9mPCkT3ruQKSUsAfVfW3JZ0P6orOPRvSopT+2Ln0WxNghZKzSg9zRkoudazPPXSXtZfVAHbZgBj5z2j6xNkbW0s8oky/hd/MC8VG/ndKaCwY6x+WoOC629BcDqzkNWvJ2EBsvhxWr+uRAqZIqTHVqWx5teRvreoaxrb8UGbjoHHVzVhpFxNJvsuEgYkJq+iK+ee6BeGjFVhyz56xgGz03qtIDF715d+y9c+1yvY2Abm5bwcZ2X9tuhqV/1v5z8Ndn12t9N4s4/Y6WHPqKFUwRSL+9YMfOtZmzHJUVmFWpAkJW9W7qRbXaO1mjLW+j7FR8GcvrEwu7pwSzuIrLsa5nGBWXZ7LgfMEWLf2o9q0arL0SJvrXmzGGudNasb53KLYMYkD6Tap2+6uLjoist9AqzK6TwsibgQlP+nOmtuGKf94/sq2apd9Ien+9aCuElv6r/k1vhnX4/Xcfgu/Di5UX6+XXgyw0fVr+sbM1F0zVOwq52Lk2s4S06rpmHQYny1Mj8aASZN07K5Ah1Jq3scOXUFZt9WrjHLtXNxjzFv6oOC5WbukHACzcKX1doZZ8qOnLln7SYJ12pjh3WhvW9Q6jb7iMfXfuCrYT6TdLajl6z1k4WkjUFgeXXfw8j4+etKgpbYuY8KSvmgom6aLNAmmT7S12ULq1maUIdDJ9s7BuqF2x/fYWO0a6WWjBSWhmOCgh65lDLURCNpvcRl7Qt991+K54Yd0OnHfYLgC8Gc3dL27CD5asAOCtG5AWLULpg5imn0C+aX1Cc6a24sUNm1F23IilX8j5JddHyKkqGj9vP2QeLnrz7onF/LLE2BKWmwDV1PstB8wBABwwf1pT26akK3IKiUXIxlr0ThYSBZ1fZ2see/srWqkIcmYGUR9JEM/jnw5qTqVCeWBptqEvZuQ2y14hTZ8INWdbeOO8qfjTJcdgml9hMmdZeHFDX/CdedMy0PQFR249mj4dWxrMndaGzX1F9AyWR0Te+cJb98OnFAveiP0oZ7MRIXxgUlj68Q5y+hvnYPWVZzW97V9cdEQkjlnUzUfTkXvHJ47FKVfdl/nvtgkFyX598ZF4deuAMgpiVhMtfRFXn3cwrj7v4Mx/V7Y0mz18V4vTz6wN/z8ZIyqjRDzvXWe0ZzKrasmFjly5xHhSqG1ag2mu4IBWOXKzXnTnA8fsXnOfkYzmm/CkP5qhke2FXCSKR8zkbHZCTzVQmdysQWfU2eqt9StaUf92/EJc48sCMzuaZ+mPBOTkoGaH94k0/67Dd21qI2SMqJ4bItu8zbDk08dnMgBRnP7371keK1mdFEWXlvTnTAvLK4ykI7caRtIInPCkP5ZkFDGTc7QGo2nt+cwqM8ogClD5FD53+j64a9lGvLyxv6maPgA8+9+nNjXKZSTqo6hw/2dPyERSUYEkJOqXKn9xTljmz7IYWAYXmTT979z5Cvad2yV9pibftDMMscy6WCiP4vSztvTrwUgagROf9EfY6VYNEUt/FIjj/s+egJ2ntqaqxlkNpMkmVQ+lstLNlneaNZMZLZBB3ZLLNgFL1QYZBKooIRoQ6PpmMcNpyVng3EuelEssy+Rb8DO90xotorxzwLypwWsaTEZjFj6SfDDxSX8sWfqF0bX0p7R4STaq2vpZgOKskxxSRPrNkHd+dsHheHZtb+a/Ww9GLGKzie1Qj5CLEYogMsxy8R9RwhmWSF+Wd1pyHumnLUneVrBx4C7TcNI+O0UG0fwIFj2TMSblHcaYDWApgLWc87cwxnYHcCOAmQAeB/BeznmJMdYC4GcADgWwFcC/cs5X+79xGYALATgAPso5byyYXANji/RH15FLIXDi9Pg3Fx8ZC5XTBaXoJ02P37yoG3e8sDGydGVWOHavbhy7V3fmv1sPslpfOAm7z2rHlv5iUxefoTBm6pcqwyBnZU/6ooRDy32SRS8Ppi15C31FYEpL+v7z50uOiW2j52IkSjDIGEmeaqQXfQzAMuH91wFcxTnfE8B2eGQO//92f/tV/n5gjO0H4DwAbwBwOoAf+ANJUzGWyh2MtqZP4ZNi20fsMTOSrayDcw+dDyAkiiTN9X/feTAe+NwJo6aJNwvN5ohr37sY171vcVPXTn37Id493N3PslWTvl8HPkP5TDQQSHYk4yS+EE/2g44I0vTtUSD9kXwm6mJExth8AGcBuM5/zwCcCOD3/i4/BfBP/uuz/ffwPz/J3/9sADdyzouc81UAlgOIr22XMUYzSkaGqOk3y5kq485PHhu8bpZ/43/OOQArvnomDtplGoBohVERrXm7qStMjRaafSendxRw8n6zm9rGu4/YFSu+eiZm+2WGlaQ/QvJOWNzNs/wfvPRE3PeZEwL5MOtlRgl5hVE0EVEvC1wN4LMAKJB2JoAezjlVMXsdwDz/9TwAawDA/7zX3z/YrvhOAMbYxYyxpYyxpZs3b67/TBIwlm7gaGj68ipPzQBj3jqgl56xD/760TdhQQY1WQxGFnQPSWZQkf5sfwHzLB3lorxD8fpkHNExzJ3Whl1ntqNcaS7pk7U9Cob+iKLm1WOMvQXAJs7544yx45t9QJzzawFcCwCLFy/WFpuntuXRO1QeU5q+uCZms2Wnez59PBhGdqaTty28Ye7U2jtOMEwkknj/MQuwfFM/LnpzPKHoqD1m4u4XNwULAWWBfed0Yv95U1F23CDb99Iz9sVdyzbi7IOiNiH5jLJeW5pAfo3RkHdGEvUwzzEA3sYYWw3PcXsigO8AmMYYo6s/H8Ba//VaALsAgP/5VHgO3WC74juZ46Z/PxpfeOt+Yypkc3p7qMk229LffVYHFszqaKrzz8DDaDj+moWu1jy++86Dg9ILIo5aOBNAtgXJ9uiegr985E2R9Wh36mrBN849MFb5lRbm6WySpU8JwSOpDryzWcl2VVDz6nHOLwNwGQD4lv6nOefvZoz9DsA58AaC8wH82f/Kzf77h/zP7+acc8bYzQB+xRj7NoC5ABYBeDTTsxGwsHsKFnanrwKYJaKW/sh0rHxCJM2Vb9/fyDAZYeJQfnW8cd5U/N8HDsOhu03P/LdF46zWs9EsS58yjEdyEP/a2/fH196+f+0dM0Saq/c5ADcyxr4C4EkA1/vbrwfwc8bYcgDb4EXsgHP+PGPstwBeAFABcAnnvDlZQmMUYnp3M6tsikh6gM4bBQtjwmKysD6AE/beqSm/KwY21LK0m6XpE+mPIXGgKWjo6nHOlwBY4r9eCUX0Ded8GMC5Cd+/AsAVjR7kRMSIWfoTvQcbTAhEKk7WkCSbZemT49ho+gZNwUjlD4yl6KWJimYnZ00GiM9DrT7bmUFylgqUSa7yZ0wkTPgyDGMVYymqyCAdjhulTOCxgKMXzoz4qnQhyju1no3WJq1s9Y5D5mOwWME7j5jY0qch/VHCRMtKnaxY9bUzJ1T0TqP41f87MpPfydWh6R+z50z8Y/nWpl1v22J4fx2178c7DOkbGGjg1o++GT1DpUlN+FlC1PST/FDXn38YhkqTKvajKTCa/gjj/UcvGO1DMMgA+83twtEL09UsMgghEn2Spd+at5taf2iywJD+COO/3/aGEVmqUcbJ+zYn1M7AIAuIOr7xdzUXRt6ZBFh95VlBirmBwViEmJxl/F3NhbH0JwmM9mwwllEQHLnNis4x8GCuroGBwahDtPQLJqGwqTBX18DAYNQh6vhmVtpcGNI3MDAYddBiKiO1uNBkhiF9AwODUQeVYTDSTvNhrrCBgcGogzJyk0qBG2QHc4UNDAxGHSTrmKqwzYe5wgYGBqMOInsj7zQf5gobGBiMOkjTN47c5sOQvoGBwajDyDsjB3OFDQwMRh2UnGVIv/kwV9jAwGDUkTfROyMGc4UNDAxGHZafhdtiLP2mw1xhAwODUUfZcQEA+Zxx5DYbhvQNDAxGHQHpG0u/6TBX2MDAYMxgSotZ4qPZMFfYwMBg1HHsom586LiFuPjYPUb7UCY8DOkbGBiMOnK2hUvP2Ge0D2NSwMg7BgYGBpMIhvQNDAwMJhEM6RsYGBhMIhjSNzAwMJhEMKRvYGBgMIlgSN/AwMBgEsGQvoGBgcEkgiF9AwMDg0kExjkf7WNIBGNsM4BXU/zELABbMjqciQZzbarDXJ9kmGtTHWPh+uzGOe9WfTCmST8tGGNLOeeLR/s4xiLMtakOc32SYa5NdYz162PkHQMDA4NJBEP6BgYGBpMIE530rx3tAxjDMNemOsz1SYa5NtUxpq/PhNb0DQwMDAyimOiWvoGBgYGBAEP6BgYGBpMIE5L0GWOnM8ZeYowtZ4xdOtrHMxpgjP2EMbaJMfacsG0GY+wOxtgr/v/p/nbGGPuuf72eYYwdMnpH3nwwxnZhjN3DGHuBMfY8Y+xj/nZzfQAwxloZY48yxp72r88X/e27M8Ye8a/DbxhjBX97i/9+uf/5glE9gREAY8xmjD3JGLvFfz9urs2EI33GmA3g+wDOALAfgHcyxvYb3aMaFdwA4HRp26UA7uKcLwJwl/8e8K7VIv/vYgDXjNAxjhYqAD7FOd8PwJEALvH7iLk+HooATuScHwjgIACnM8aOBPB1AFdxzvcEsB3Ahf7+FwLY7m+/yt9vouNjAJYJ78fPteGcT6g/AEcBuF14fxmAy0b7uEbpWiwA8Jzw/iUAc/zXcwC85L/+EYB3qvabDH8A/gzgFHN9lNemHcATAI6Al2Wa87cHzxmA2wEc5b/O+fux0T72Jl6T+fCMghMB3AKAjadrM+EsfQDzAKwR3r/ubzMAZnPO1/uvNwCY7b+etNfMn24fDOARmOsTwJcvngKwCcAdAFYA6OGcV/xdxGsQXB//814AM0f0gEcWVwP4LADXfz8T4+jaTETSN6gD3DM9JnW8LmNsCoA/APg453yH+Nlkvz6cc4dzfhA8q/ZwAGbVcgCMsbcA2MQ5f3y0j0UXE5H01wLYRXg/399mAGxkjM0BAP//Jn/7pLtmjLE8PML/Jef8Jn+zuT4SOOc9AO6BJ1lMY4zl/I/EaxBcH//zqQC2juyRjhiOAfA2xthqADfCk3i+g3F0bSYi6T8GYJHvTS8AOA/AzaN8TGMFNwM43399Pjwtm7a/z49SORJAryBzTDgwxhiA6wEs45x/W/jIXB8AjLFuxtg0/3UbPH/HMnjkf46/m3x96LqdA+Buf6Y04cA5v4xzPp9zvgAet9zNOX83xtO1GW2nSJMcLWcCeBmeDvkfo308o3QNfg1gPYAyPI3xQnha4l0AXgFwJ4AZ/r4MXsTTCgDPAlg82sff5GvzJnjSzTMAnvL/zjTXJ7g+BwB40r8+zwH4L3/7HgAeBbAcwO8AtPjbW/33y/3P9xjtcxih63Q8gFvG27UxZRgMDAwMJhEmorxjYGBgYJAAQ/oGBgYGkwiG9A0MDAwmEQzpGxgYGEwiGNI3MDAwmEQwpG9gYGAwiWBI38DAwGAS4f8DZdjI2lrJtnAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "baseline = np.array([np.random.normal(np.random.randint(4100, 4200), 100, size=26)])\n",
    "rest = [np.random.normal(np.random.randint(4100, 4200), 100, size=20) for x in range(10)]\n",
    "tasks = [np.random.normal(np.random.randint(4300, 4600), 100, size=20) for x in range(10)]\n",
    "Y = np.append(baseline, np.hstack(np.concatenate(list(zip(tasks, rest)))))\n",
    "\n",
    "plt.plot(Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 1: Single regressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1b897aa01d0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI4AAAD8CAYAAAChF5zCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAO9ElEQVR4nO2dfYxcdbnHP8/MvtHSbqWF3lKKNFg11UTU5hLCTUTUG6jESoJANQqG+BYaNb5RNPdKvJJo4vVKosFbIwJGqRUlIqnWNxo0EQUUtYAvtVZprVZeurSW292dee4f52wzLLszZ86cmfM8M88n2ezs7Oycp5lPz5xz5vd8H1FVgqBdKmUXEPgkxAlyEeIEuQhxglyEOEEuQpwgF10TR0QuEJHfichuEdncre0E5SDduI4jIlXg98BrgH3AfcBGVX248I0FpdCtPc6/ArtVdY+qTgJbgQ1d2lZQAkNdet6VwKMNP+8Dzp7vwSPVBXrC8HiXSuktp6w5VHYJhXFw/xQTT0zLXL/rljgtEZG3A28HGBtazDlnXFFWKYXy7m/dVXYJhfGeDX+c93fdeqvaD6xq+Pm09L7jqOoWVV2nqutGqgu6VEbQLbolzn3AGhFZLSIjwOXAnV3aVlACXXmrUtVpEdkE7ACqwE2q+lA3thWUQ9eOcVR1O7C9W88flEtcOQ5yEeIEuQhxglyEOEEuQpwgF6VdOW5Eq0J94VjZZQRtYEKcqUVV9r+qPz6rGhRMiDM8PsmpF/6l7DKCNjAhzorRCT6yuj8+HDxaHy27hJ5gQpyFopw9OlV2GYVw99MhTs+oIIzKcNllBG1gQpw6ytH6ZNllFETscXrGtNZ5rG/EGQxsiINwqG6ilCAjJl6tOsLh+kjZZQRtYEKcYzrM3qmTyy6jEJZWj5RdQk8wIc6R2ij3TDy/7DIK4eKTflF2CT2hI3FEZC9wGKgB06q6TkROAr4GnAHsBS5V1SebPc/hiQXs3HFWJ6WY4eKNIU5WXqmqjzX8vBn4oap+Im393Qxc0+wJRp6qc/qOpwsoxQAbyy6gN3TjrWoDcF56+xZgJy3EkacnGd71py6UEnSLTsVR4HsiosD/quoWYLmqHkh//zdg+Vx/+IyGPBZQOzTRYSlBL+lUnH9T1f0icgrwfRH5beMvVVVTqZ5FKtkWgMVykiJzdpoGRulIHFXdn34/KCJ3kIQN/F1EVqjqARFZARxs9TxSqVBZEN2cnsgtjogsBCqqeji9/e/Ax0g6Nq8APpF+/1bLJ6tUkBDHFZ3scZYDd0jyFjMEfFVVvysi9wHbROQq4M/ApS2fqSLI2GB8ONgv5BZHVfcAL5nj/seBV7X1ZCIwbOJaZJARE69WfWyII2tPKbuMoA1MiDO5BP7yuhgN4AkT4py88DDvPGdn2WUEbWBCnJOqR3nT+C/LLqMQdk0uLbuEnmBCnBGpctrQiWWXUQi7BmQhowlxFKWm9bLLCNrAhDg1lCfrffLpOCeUXUBPMCHOtMITscNxhQlxagiHYs2xK0yIU1cZmNbZfsGEOE/UFnLb4/MGr7si1hz3kFhz7A8T4lSPwZI/lF1F0A4mxBk6PMWyH/+17DKCNjAhjk5NUnt0f+sHBmYwIQ4KOj1ddhVBG7QUR0RuAi4CDqrqi9P75my6k2Q54A3AeuAocKWqtj5aFEGG4zqOJ7LscW4GPgvc2nDffE13FwJr0q+zgRtpMuBsBqlUqJy4sL3Kg1JpKY6q3iMiZ8y6e76muw3ArZoM+rxXRJbMdDw03UisOXZH3mOc+Zru5hqpuBJoLo5UYDTeqjzR8cFxs6a7Zjyjk3N4nPp4vFV5Iq848zXdtRypOENjJ+cJ/7JK97/6OTlLCcogrzjzNd3dCWwSka0kB8UTLY9vSAKyl69/tNXDAkNkOR2/jeRAeJmI7AM+SiLMXE1320lOxXeTnI6/NUsRK0Yn+I/V3267eIsMyqf8Wc6q5kt8eVbTXXo2dXW7RURAtj9MXDmOgGx/mBAnArL9YUKcCMj2hw1xIiDbHSZerQjI9ocJcSIg2x8mxDlwdJz/+tX6sssohM+8bFvZJfQEE+JUJyqM39UfveO8rOwCeoMNcY7VWbynX1qABwMT4kRAtj9MiKO1WgRkO8OEOAARkO0LE+JEQLY/TIgTAdn+MCJOLFb3hg1xIiDbHSZerQjI9kfeTs7rgLcB/0gf9mFV3Z7+7lrgKpJxi+9W1R2tthEB2f7I28kJ8D+q+qnGO0RkLXA58CLgVOAHIvJ8Va0120AEZPsjbyfnfGwAtqrqMeBPIrKbZIbVT5v9UQRk+6OTY5xNIvIW4H7g/emk35XAvQ2PmenkfBaNDXmnrxyKgGxnVHL+3Y3AmcBZJO29/93uE6jqFlVdp6rrli2tUNN6X3wNCrn2OKr695nbIvIF4K70x8ydnI1EQLY/cokzK4HiYmBXevtO4Ksi8mmSg+M1wM9bPV8EZPsjbyfneSJyFsn46L3AOwBU9SER2QY8DEwDV7c6o4IIyPZI3k7OLzZ5/PXA9e0UEQHZ/jBx5TgCsv1hQpwIyPaHCXFGnqpz+o4+OauaL6KhzzAhTqw59ocJcbReo37kn2WXEbSBCXEiINsfNsSJgGx3mBAnArL9YUKcWHPsDxviREC2O2yIU61EQLYzTIgzubhKBGT7woQ4EZDtDxPiREC2P0yIEwHZ/jAhTgRk+yPLCsBVJD1Vy0lW/G1R1RuKHK8YAdn+yLLHmSZpf/mFiCwCHhCR7wNXUtB4xQjI9keWpaMHSCfcqephEXmEpFeqsPGKEZDtj7ZerbSj86XAzyhwvGIEZPsjszgiciLwDeC9qvqUNESv5Rmv2NjJufTUkQjIdkYmcURkmESar6jqN9O7Oxqv2DhacezMlRoB2b7IclYlJO0wj6jqpxt+Vdh4xQjI9keWPc65wJuB34jIg+l9H6bA8YpDh6dY9uO/tld5UCpZzqp+AsyXJVvIeEWdmqT2aMsW88AQNs6BY82xO2yIAxGQ7QwT4kRAtj9MiBMB2f4wIk4sVveGDXEiINsdJl6tCMj2hwlxIiDbHybEiYBsf5gQJwKy/WFCnBGpRkC2M0yIo+hAhUv3AybEiYBsf5gQJwKy/WFCnAjI9ocJcSIg2x8mxImAbH900sl5HQWNV4yAbH900skJBY1XjIBsf3TSyTkfbY9XjIBsf3TSyXkuHYxXbGzIG2MBtUMTeeoPSiLzaMXZnZx0OF6xcbTiMKPJmpx++BoQcndyFjleMdYc+yN3J2eh4xVFkJG4AOiJTjo5NxY2XjHWHLujk07O7U3+pr3xihGQ7Q4TV44jINsfJsSJgGx/mBAnArL9YUKcCMj2hwlxIiDbHybEiYBsf5gQJwKy/WFCnAjI9ocNcSIg2x0mXq0IyPaHCXGO6XAEZDvDhDgHjo4TAdm+MCFOBGT7w4Q4EZDtDxPiREC2P0yIEwHZ/siydHQMuIfkkugQcLuqflREVgNbgaXAA8CbVXVSREZJGvheDjwOXKaqe1tsBBmO03FPZNnjHAPOV9Uj6aL1n4jId4D3kTTkbRWRz5N0bt6Yfn9SVZ8nIpcDnwQua7YBqVSonBgLuTyRZemoAjMXJ4bTLwXOB96Y3n8LcB2JOBvS2wC3A58VEUmfZ25izbE7srbHVEnejp4HfA74I3BIVWcOTBqb7o6PVlTVaRGZIHk7e6zJBmDIxuFWkI1Mr1bapXCWiCwB7gBe2OmGn9HJObQoArKd0darpaqHRORu4BxgiYgMpXudxqa7mYa8fSIyBIyTHCTPfq7joxUXLTlNIyDbF1nOqk4GplJpTgBeQ3LAezdwCcmZ1ezRileQhAxcAvyo6fENEZDtkSx7nBXALelxTgXYpqp3icjDwFYR+TjwS5JuT9LvX05TKp4giTxpSgRk+yPLWdWvSRIqZt+/hyS+ZPb9/we8oZ0iIiDbHyaOSCMg2x8mxImAbH+YECcCsv1hQpwIyPaHCXEiINsfJsSJgGx/mBAnArL9YUKcCMj2hwlxIiDbHybEiYBsf5gQR2u1CMh2hglxgIEKl+4HTIgTAdn+MCEOlQoS4rjCiDixWN0bNsRBoJp5HklgABviREC2Ozrp5LwZeAUwcx59pao+mA4NuQFYDxxN7296OTUCsv3RSScnwAdV9fZZj7+QZGLMGuBskia9ph9ERUC2Pzrp5JyPDcCt6d/dKyJLZo0oehYRkO2PXJ2cqvozEXkXcL2I/CfwQ2BzOofzeCdnykyX54FZz3m8IW/VymoEZDsjVyeniLwYuBb4GzBC0lh3DfCxrBtubMhb95IxjYBsX+Tt5LygYWz0MRH5EvCB9Oe2RytGQLY/cndyzhy3pGdRr+eZoxU3ichWkoPiiWbHNxAB2R7ppJPzR6lUAjwIvDN9/HaSU/HdJKfjb221gQjI9kcnnZznz/N4Ba5up4gIyPaHif/mEZDtDxPiREC2P0yIEwHZ/jAhTgRk+8OEOBGQ7Q8T4kRAtj9siBMB2e4wIU4EZPvDhDix5tgfNsSRCozGW5UnbIgTa47dYUKc+kiFo6sWlV1G0AYmxImAbH+YECcCsv1hQpwIyPaHCXEiINsfmcVJVwDeD+xX1YuKHK0YAdn+aGeP8x7gEWBx+vMnKWi0YgRk+yNrX9VpwGuB64H3pQvUCxutGAHZ/si6x/kM8CFg5mLLUgocrRgB2f7I0h5zEXBQVR8QkfOK2nBjJ+cppw4NTOtsv5Blj3Mu8DoRWQ+MkRzj3ECBoxVPXrtUIyDbF1naY64lafcl3eN8QFXfJCJfp6DRihGQ7Y9OruNcQ0GjFSMg2x/t9o7vBHamtwsbrRgB2f4wceU4ArL9YUIcIAKynWFCnAjI9ocJcSIg2x9GxInF6t6wIQ4CQ9WyiwjawIY4AhrJ6q6wIQ7EWZUzpMWnAb0pQuQfwJ97sKllNPmU3hm9+Lc8V1XnTLwyIU6vEJH7VXVd2XUUQdn/ljiwCHIR4gS5GDRxtpRdQIGU+m8ZqGOcoDgGbY8TFMRAiCMiF4jI70Rkt4hsLrueThCRm0TkoIjsav3o7tH34qSNhJ8jGcC2FtgoImvLraojbgYuKLuIvheHZJXiblXdo6qTJGukN5RcU25U9R6SJbmlMgjizDd4LeiAQRAn6AKDIE7bg9eC1gyCOPcBa0RktYiMkLTr3FlyTe7pe3HSTtNNwA6StI1tqvpQuVXlR0RuI2l2fIGI7BORq0qpI64cB3no+z1O0B1CnCAXIU6QixAnyEWIE+QixAlyEeIEuQhxglz8P6DM6XFqg/5XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X1 = np.array([np.convolve(np.repeat(10 * [[0] * 20 + [1] * 20], 1), hrf(np.array(range(len(Y)))))[:len(Y)],              \n",
    "              [1] * len(Y)]).T\n",
    "model_1 = GeneralLinearModel(X1)\n",
    "\n",
    "plt.imshow(model_1.X, aspect=0.01, interpolation='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "task > rest: 327.57 (z = 20.94)\n"
     ]
    }
   ],
   "source": [
    "model_1.fit(Y)\n",
    "con = model_1.contrast(np.array([1, 0]))\n",
    "print(f\"task > rest: {con.effect[0,0]:.2f} (z = {con.z_score()[0]:.2f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 2: Two regressors, 100% colinearity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1b899b17ef0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAD8CAYAAADOpsDvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARy0lEQVR4nO3de4xc5X3G8e8ze7ExGLsYcI1xihVMK4IECRYIUbUUSgUuioOUgElEgFq5VFgQJWkxVG1QEktEIhAkKlJH0OAosHFJEC5y6yaARSIB4drWmCZ1jIntGMzFXuya+LL76x/nLBrbu7OzOzvzvvF5PtJqd2fn8lrMw5lzZs7vUURgVmW11AswS80hsMpzCKzyHAKrPIfAKs8hsMprWwgkXSLpF5I2SFrarscxa5Xa8T6BpC7gl8DFwBbgWeCqiFg/4Q9m1qJ2bQnOATZExMaI2Af0AQvb9FhmLelu0/3OBjbX/b4FOHekK/d2TYmjeqa1aSnWihPn7Uy9hAmxfet++t85oOH+1q4QjErSZ4HPAkzuPpbzTrkm1VKsgRseeTT1EibEjQt/NeLf2vVyaCswp+73k8vL3hcRyyNifkTM7+2a0qZlmI2uXSF4Fpgnaa6kXmARsKpNj2XWkra8HIqIA5KWAGuALuC+iHi5HY9l1qq27RNExGpgdbvu32yi+B1jqzyHwCrPIbDKcwis8hwCq7xk7xjXiy4xePTk1MuwisoiBPundrH1In92yNLIIgQ90/Zx0qW/Tr0Mq6gsQjBrUj9/N/fI+KDWkWbP4KTUS2i7LEJwtIJzJ+1PvQwbxhPvOQQdUUNMUk/qZVhFZRGCQYI9g/tSL8OG5S1BRxyIQd5yCCyRPEKA2DmYxVKsgrJ45g0idg32pl6GVVQWIdgbPWzaf0LqZdgwZnTtTr2EtssiBLsHJvFk/2mpl2HDuPy4F1Ivoe1aCoGkTcAuYAA4EBHzJR0H/AA4BdgEXBEROxrdz67+Kaxdc1YrS7E2ufwqh6AZfxYRb9X9vhR4LCJuK8cvLgVuanQHve8O8oE1703AUmzCXZV6Ae3XjpdDC4ELyp/vB9YySgj03j561r3ahqWYja7VEATwH5IC+KeIWA7MjIht5d9fB2YOd8ODhm8xhYGd/S0uxWx8Wg3BH0fEVkknAj+W9D/1f4yIKANymDIwywGO1XGBhp2QZ9Z2LYUgIraW37dLephiEO8bkmZFxDZJs4Dto92PajVqUzyFztIYdwgkHQ3UImJX+fNfAF+lmDR3DXBb+f2RUe+sVkMOgSXSypZgJvCwipcx3cADEfHvkp4FVkpaDLwGXDHqPdWEJh/5H9SyPI07BBGxEThzmMvfBi4a051J0JPF+3ZWQVk88wYnd7P79BNTL8MqKosQ7JsOv/7oxNdGmTUjixCccPQuPn/e2tTLsIrKIgTHde3hU9NeTL0MG8a6fTNSL6HtsghBr7o4ufuY1MuwYayrwAl/WYQgCAZiMPUyrKKyCMEAwY5Bf4o0T0elXkDbZRGCAwHveENgiWQRggHETp9jbIlkEYLBUCXG/VmesgjBOwNH8+DbIxbeW0I+x7hDfI5xvnyOcYd07YXp/5t6FVZVWYSge9d+jv/pb1IvwyoqixDE/n0MbN6aehlWUVmEgIA4cCD1KqyiRg2BpPuAy4DtEXFGedmwA7ZUnGZ2F7AA2ANcGxGj71lJqMfvE1gazWwJvgvcDayou2ykAVuXAvPKr3OBe8rvDalWo3bM0WNbudkEGTUEEfGkpFMOuXikAVsLgRUREcDTkqYPTZ5o+CA+x9gSGu8+wUgDtmYDm+uut6W8rHEIVINJfjlkabS8Y9xowFYjB02g65nG4DS/HLI0xhuCkQZsbQXm1F3v5PKyw9RPoDvq9+fE1j//vXEuxaw14w3BSAO2VgFLJPVR7BD3j7o/QFHmPXPB5tGuZtYWzRwifZBiJ/h4SVuAr1A8+YcbsLWa4vDoBopDpNc1s4hZk/r5+7n/OubFW/tV4dO9zRwdGmlC/WEDtsqjQtePdREu886Xy7w7xGXellIWIXCZd868JegIl3lbSnmEwGXellAWzzyXeVtKWYTAZd75cpl3h2zbM42v/eeC1MuwYXzrIytTL6HtsghBV3+NaY96FmmWPpJ6Ae2XRwj2DnLsRo9htDSyCIHLvC2lLEIQAwMu87ZksggBgMu8LZUsQuAyb0spixC4zNtSyiQEPtHe0skjBC7ztoSyeOa5zNtSGu8EuluBzwBvlle7JSJWl3+7GVgMDAA3RMSa0R7DZd6W0ngn0AHcGRG3118g6XRgEfAh4CTgJ5JOi4iBRg/gMm9LabwT6EayEOiLiL3Aq5I2AOcATzW6kcu88+Uy78aWSPo08BzwpYjYQTFt7um66wxNoDtM/fCtD8zudpl3plzmPbJ7gK8BUX7/JvBXY7mD+uFbZ585KVzmbamMKwQR8cbQz5K+Azxa/tr0BLp6LvPOmcu8h3XIpOnLgXXlz6uAByTdQbFjPA/4+Wj35zJvS2m8E+gukHQWxcuhTcDnACLiZUkrgfXAAeD60Y4Mgcu8La3xTqC7t8H1lwHLxrIIl3lbSlm8Y+wy73y5zLtDXOadL5d5d0jvu4N8YI2PDmVppHHMR5AsQuBzjC2lLEIQgwMM7v6/1MuwisoiBC7ztpTyCIHLvC2hLELgMm9LKYsQ+BxjSymPELjM2xLKIwRdNZd5WzJZhGDfsV24zNtSySIELvO2lLIIgcu881WFT/dmEQKXeefLZd4d4jJvS6mZM8vmUMwcmklxJtnyiLhL0nHAD4BTKM4uuyIidkgScBewANgDXBsRDT+P6zLvnHlLAMVpkl+KiBckTQWel/Rj4FrgsYi4TdJSYClwE3ApxbnF84BzKSZTNDxjxmXellIzp1duA7aVP++S9ArFLKGFFOceA9wPrKUIwUJgRUQE8LSk6YecmH8Yl3lbSmN65pWT6D4MPAPMrHtiv07xcgmKgNQf7xwawDViCFzmbSk1HQJJxwA/BL4QEe+qrl4pIkLSmCbq1k+gm3FSr8u8M+Uy75KkHooAfD8iflRe/MbQyxxJs4Dt5eVNDeCqn0A3+YOzw2XeeXKZN1Ae7bkXeCUi7qj70yrgGuC28vsjdZcvkdRHsUPc32h/AFzmnTWXeQNwPnA18N+SXiovu4Xiyb9S0mLgNeCK8m+rKQ6PbqA4RHrdqIvYtZ/jf/qbsa3cbII0c3ToZ8BI/aoXDXP9AK4fyyJi/z4GNo86stSsLfI4LulzjC2hPEIALvO2ZLIIgcu8LaUsQuAyb0spkxD4RHtLJ48QuMzbEsrimecyb0spixC4zNtSyiIELvO2lLIIgcu88+Uy7w7pVZfLvDPlMu8OCQKXeVsqWYTAZd45c5l3R7jM21LKIgQu87aUsgiBy7wtpSxC4DLvfLnMm4YT6G4FPgO8WV71lohYXd7mZmAxMADcEBFrGj2Gy7zz5TLvwkgT6ADujIjb668s6XRgEfAh4CTgJ5JOi4iBkR7AZd4Zc5l3wwl0I1kI9EXEXuBVSRuAc4CnRrqBy7wtpVYm0J1PMVrl08BzFFuLHRQBebruZkMT6A69r/eHb01mCgM7+8ezfrOW1Zq94qET6CgG7X4QOItiS/HNsTxwRCyPiPkRMb+HScU5Bf7K76sCxj2BLiLeqPv7d4BHy1+bmkB30P37HGNLaNwT6A6ZNH05sK78eRXwgKQ7KHaM5wE/H+VBUK/fLLM0WplAd5WksygOm24CPgcQES9LWgmspziydH2jI0OAzzG2pFqZQLe6wW2WAcuaXoXLvC2hLN4xdpm3pZRFCFzmbSllEQKXeVtKWYTAZd75qsKne7MIgcu88+Uy7w5xmbellEUIXOadM28JOsJl3pZSHiFwmbcllMUzz2XellIWIdgbPS7zzpTLvDtk255puMw7Ty7z7hCXeWfMZd6d4TJvSymLELjM21LKIgQu87aUmjm9cjLwJMVbh93AQxHxFUlzgT5gBvA8cHVE7JM0iWJY19nA28CVEbFplAdBPT5Eamk0syXYC1wYEbvLE+5/JunfgC9SDN/qk/Rtiolz95Tfd0TEqZIWAd8Armz0AKrVqB3jk2osjWZOrwxg6GBxT/kVwIXAJ8vL7wdupQjBwvJngIeAuyWpvJ/h+RxjS6jZkStdFC95TgX+EfgVsDMihl7I1w/Ymg1sBoiIA5L6KV4yvdXgAaA7j90Tq56mnnnltIizJE0HHgb+qNUHPmgCXfdUl3lbMmN65kXETklPAOcB0yV1l1uD+gFbQ8O3tkjqBqZR7CAfel/LgeUAU6efHC7ztlSaOTp0ArC/DMBRwMUUO7tPAB+nOEJ0DfBIeZNV5e9PlX9/vOH+AC7ztrSa2RLMAu4v9wtqwMqIeFTSeqBP0teBFymm1FF+/145jfodijHtDbnM21Jq5ujQf1FMoj708o0UI9cPvfy3wCfGsgiXeefLZd4d4jLvfLnMu0Nc5m0pZRECl3nnzGXeHeEyb0spixC4zNtSyiIELvO2lLIIgcu88+Uy7w5xmXe+XObdIS7zzpjLvDvDZd6WUhYhiIEBl3lbMlmEAKhMcbTlJ4sQuMzbUsoiBNRqyCGwRDIJgU+0t3TyCAGCrlrqRVhF5RECl3lbQq1MoPsu8KfA0LHNayPiJUkC7gIWAHvKyxu+7egyb0uplQl0AH8TEQ8dcv1LgXnl17kUA7kafjDIZd6WUisT6EayEFhR3u5pSdMlzYqIbSPdwGXe+arCp3vHNYEuIp6R9NfAMkn/ADwGLI2IvdRNoCsNTafbdsh9vj98a87sLpd5Z8pl3qVDJ9BJOgO4GXgd6KUYonUT8NVmH7h++Nb8MyeHy7wtlfFOoLskIm4vL94r6Z+BL5e/D02gG1I/nW5YLvPOmbcEI06gG3qdXx4N+hiwrrzJKmCJpD6KHeL+RvsD4DJvS6uVCXSPlwER8BLw+fL6qykOj26gOER63WgP4DJvS6mVCXQXjnD9AK4fyyJc5m0pZfG/X5d558tl3h3iMu98ucy7Q1zmnTGXeXeGy7wtpSxC4DJvSymLELjM21LKIwQu87aEsgiBy7wtpSxC4HOMLaU8QqAaTPLLIUsjjxD4HGNLKIsQDPbW2DNnauplWEVlEQKXeVtKWYTAZd6WUhYhcJl3vlzm3SEu886Xy7zrlGeWPQdsjYjLJM0F+oAZFJMoro6IfZImASuAs4G3gSsjYlOj+3aZt6U0li3BjcArwLHl798A7oyIPknfBhZTDNpaDOyIiFMlLSqvd2WjO3aZd85c5g2ApJOBvwSWAV8sT66/EPhkeZX7gVspQrCw/BngIeBuSSpPuxyWy7wtpWa3BN8C/hYYOpg/A9gZEUMf/RwasAV1w7ci4oCk/vL6b4105y7ztpSaGblyGbA9Ip6XdMFEPXD9BLoTT+quxLg/y1MzW4LzgY9KWgBMptgnuAuYLqm73BrUD9gaGr61RVI3MI1iB/kg9RPoTjh9RrjMO08u8wYi4maKkYuUW4IvR8SnJP0L8HGKI0TXAI+UN1lV/v5U+ffHG+0PgMu8c+Yy78ZuAvokfR14Ebi3vPxe4HuSNgDvAItGuyOXeWfMZd4Hi4i1wNry543AOcNc57fAJ8Zyvy7ztpSyeMfYZd6WUhYhAFzmbclkEQKXeVtKWYTAZd6WUiYh8In2lk4eIUDQ3ZV6EVZReYRAEG60t0TyCAH46JAlo1E+0dCZRUhvAq914KGOp8GnWX/H+N8yNn8QEcM2wWQRgk6R9FxEzE+9jongf8vE8QtxqzyHwCqvaiFYnnoBE8j/lglSqX0Cs+FUbUtgdphKhEDSJZJ+IWmDpKWp19MKSfdJ2i5pXeq1tELSHElPSFov6WVJNyZby5H+cqgcGvZL4GKKqRjPAldFxPqkCxsnSX8C7AZWRMQZqdczXpJmAbMi4gVJUykGuH0sxX+XKmwJzgE2RMTGiNhHcU70wsRrGreIeJLitNXfaRGxLSJeKH/eRTHYbXbjW7VHFULw/hykUv2MJMuApFOADwPPpHj8KoTAMibpGOCHwBci4t0Ua6hCCIbmIA2pn5FkCUnqoQjA9yPiR6nWUYUQPAvMkzRXUi/FCJhViddUeeU823uBVyLijpRrOeJDUE7IWwKsodj5WhkRL6dd1fhJepBisNkfStoiaXHqNY3T+cDVwIWSXiq/FqRYyBF/iNRsNEf8lsBsNA6BVZ5DYJXnEFjlOQRWeQ6BVZ5DYJXnEFjl/T/Ny5Ji0Ihl2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X2 = np.array([np.convolve(np.repeat(10 * [[0] * 20 + [1] * 20], 1), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "               np.convolve(np.repeat(10 * [[0] * 20 + [1] * 20], 1), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "              [1] * len(Y)]).T\n",
    "model_2 = GeneralLinearModel(X2)\n",
    "\n",
    "plt.imshow(model_2.X, aspect=0.01, interpolation='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "task > rest: 327.57 (z = 20.93)\n"
     ]
    }
   ],
   "source": [
    "model_2.fit(Y)\n",
    "con = model_2.contrast(np.array([1, 1, 0]))\n",
    "print(f\"task > rest: {con.effect[0,0]:.2f} (z = {con.z_score()[0]:.2f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 3: Two regressors, ~50% colinearity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1b899b85978>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAD8CAYAAADOpsDvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAASTklEQVR4nO3de4xc5X3G8e8zu2sbDNgFAjXgBBqctknUkIBAiKql0FTgRnFQEy6JCKRWLhUIoiQthqoNCkEiEoEgESV1BA1ESRyXBOEit24uUBIJEgihrTENdcAEOwZzsRcTgy+7v/5xzqKx8c7O7Mzs+545z0da7e7sXF6LeThzzsz5PYoIzOqskXoBZqk5BFZ7DoHVnkNgtecQWO05BFZ7fQuBpLMk/VLSeknL+vU4Zt1SP94nkDQEPA68G9gIPAhcEBHrev5gZl3q15bgZGB9RDwREbuAFcCSPj2WWVeG+3S/RwNPN/2+EThlsivPGjowDhiZ16elWDeOWLQt9RJ6Ysum3Yy+uEf7+1u/QjAlSR8DPgYwZ/gQTj32olRLsRYuu+vu1EvoicuX/GrSv/Xr5dAmYGHT78eUl70mIpZHxEkRcdKsoQP7tAyzqfUrBA8CiyQdJ2kWcD6wqk+PZdaVvrwciog9ki4F1gBDwK0R8Wg/HsusW33bJ4iI1cDqft2/Wa/4HWOrPYfAas8hsNpzCKz2HAKrvWTvGDeLITE+d07qZfRE47evpl6CdSiLEOw+eIhNZw7GZ4cWrnIIqiaLEIzM28VRZ/869TJ6Y9Xs1CuwDmURggWzR/n74wbjg1rX8Vepl2AdyiIEcxWcMnt36mVYTWURggZitkZSL8NqKosQjBPsGN+VehlWU1mEYE+M87xDYInkEQLEtvEslmI1lMUzbxyxfXxW6mVYTWURgp0xwobdb0i9DKupLELw8ths7ht9S+plWE11FQJJG4DtwBiwJyJOknQo8B3gWGADcG5EbG11P9tHD+TeNSd0s5RsvJktqZdgHerFluDPIuL5pt+XAT+MiOvK8YvLgCta3cGsl8Z545pXerAUs8714+XQEuD08ufbgHuZIgR6ZRcja5/sw1ISOOLw1CuwDnUbggD+Q1IA/xQRy4EjI2Jz+fdngCP3d8O9hm9xIGPbRrtcSh6GHILK6TYEfxwRmyQdAXxf0v82/zEiogzI65SBWQ5wiA4NtN8JeWZ911UIImJT+X2LpDspBvE+K2lBRGyWtACm3lNUo0HjQE+hszSmHQJJc4FGRGwvf/4L4HMUk+YuAq4rv9815Z01GsghsES62RIcCdyp4mXMMPCtiPh3SQ8CKyUtBZ4Czp3ynhpCcwbjZBS3QlfPtEMQEU8A79jP5S8AZ3Z0ZxKMZPG+ndVQFs+88TnDvPzWI1IvoyfmPv5C6iVYh7IIwa758Ov3DsYLiT+8PvUKrFNZhOANc7fziVPvTb2MnvhP/ij1EqxDWYTg0KEdfGjeL1IvoyccgurJIgSzNMQxwwelXobVVBYhCIKxGE+9DKupLEIwRrB13J8itTSyCMGegBe9IbBEsgjBGGKbzzG2RLIIwXiIHeOD8bEJq54sQvDi2Fy+/cKkhfdmfZVFCHyOsaWURQiGdsL8/0u9CqurLEIwvH03h//4N6mX0RPhT8NWThb/xWL3Lsae3pR6GT3R+L03pV6CdSiLEBAQe/akXoXV1JQhkHQr8B5gS0S8vbxsvwO2VJxmdhOwGNgBXBwRD0+5CgmN+H0CS6OdLcHXgZuB25sum2zA1tnAovLrFOAr5feW1GjQOGhuZys365EpQxAR90k6dp+LJxuwtQS4PSICeEDS/InJEy0fxOcYW0LT3SeYbMDW0cDTTdfbWF7WOgRqwGy/HLI0ut4xbjVgq5W9JtCNzGN83mC8HIpGI/USrEPTDcFkA7Y2AQubrndMednrNE+gO+B3F8amP/+daS4lL+FBepUz3RBMNmBrFXCppBUUO8SjU+4PUJR5H7n46amuVgkN7xVUTjuHSL9NsRN8uKSNwGcpnvz7G7C1muLw6HqKQ6QfaWcRC2aP8g/H/WvHi7f+q8One9s5OnTBJH963YCt8qjQJZ0uwmXe+brnFYdgRrjM21LKIgQu886ZtwQzwmXellIeIXCZtyWUxTPPZd6WUhYhcJl3vg4bejn1EvouixBs3jGPa/5rcepl2H586V0rUy+h77IIwdBog3l3exZplt6VegH9l0cIdo5zyBMew2hpZBGCgSrztsrJIgQxNjYwZd5WPVmEAMBl3pZKFiFwmbellEUIXOZtKWUSgsE50d6qJ48QuMzbEsrimTdIZd5WPdOdQHc18FHgufJqV0XE6vJvVwJLgTHgsohYM9VjDFKZt1XPdCfQAdwYEXv1t0t6K3A+8DbgKOAHkt4SEWOtHmCQyryteqY7gW4yS4AVEbETeFLSeuBk4P5WNxqkMu9Bs3bXYamX0Hfd7BNcKunDwEPApyNiK8W0uQearjMxge51modvvfHoYZd5Z2ptDU74m24IvgJcQzF68xrgi8Bfd3IHzcO3TnzH7HCZt6UyrRBExLMTP0v6GnB3+WvbE+iaucw7ZwekXkDfTSsE+0yaPgdYW/68CviWpBsodowXAT+b6v5c5m0pTXcC3emSTqB4ObQB+DhARDwqaSWwDtgDXDLVkSFwmbelNd0JdLe0uP61wLWdLMJl3pZSFu8Yu8w7X+ccOnXbVtVlEYJBKvMeNOdc4BDMiFkvjfPGNT46lKXJxjEPkCxC4HOMLaUsQhDjY4y//NvUy7CayiIELvO2lPIIgcu8LaEsQuAyb0spixD4HGNLKY8QuMzbEsojBEONgSnzturJIgS7DhliUMq8rXqyCMEglXlb9WQRApd556sOn+7NIgQu886Xy7xniMu8LaV2zixbSDFz6EiKM8mWR8RNkg4FvgMcS3F22bkRsVWSgJuAxcAO4OKIaPl5XJd558xbAihOk/x0RDws6WDg55K+D1wM/DAirpO0DFgGXAGcTXFu8SLgFIrJFC3PmHGZt6XUzumVm4HN5c/bJT1GMUtoCcW5xwC3AfdShGAJcHtEBPCApPn7nJj/Oi7ztpQ6euaVk+jeCfwUOLLpif0MxcslKALSfLxzYgDXpCFwmbel1HYIJB0EfBf4ZES8pKZ6pYgISR1N1G2eQHfYUbNc5p0pl3mXJI1QBOCbEfG98uJnJ17mSFoAbCkvb2sAV/MEujlvPjpc5p0nl3kD5dGeW4DHIuKGpj+tAi4Criu/39V0+aWSVlDsEI+22h8Al3lnzWXeAJwGXAj8j6RHysuuonjyr5S0FHgKOLf822qKw6PrKQ6RfmTKRWzfzeE//k1nKzfrkXaODv0EmKxf9cz9XD+ASzpZROzexdjTU44sNeuLPI5L+hxjSyiPEIDLvC2ZLELgMm9LKYsQuMzbUsokBD7R3tLJIwQu87aEsnjmuczbUsoiBC7ztpSyCIHLvC2lLELgMu98ucx7hszSkMu8M+Uy7xkSBC7ztlSyCIHLvHPmMu8Z4TJvSymLELjM21LKIgQu87aUsgiBy7zz5TJvWk6guxr4KPBcedWrImJ1eZsrgaXAGHBZRKxp9Rgu886Xy7wLk02gA7gxIq5vvrKktwLnA28DjgJ+IOktETE22QO4zDtjLvNuOYFuMkuAFRGxE3hS0nrgZOD+yW7gMm9LqZsJdKdRjFb5MPAQxdZiK0VAHmi62cQEun3v67XhW3M4kLFto9NZv1nXGu1ecd8JdBSDdt8MnECxpfhiJw8cEcsj4qSIOGmE2cU5Bf7K76sGpj2BLiKebfr714C7y1/bmkC31/37HGNLaNoT6PaZNH0OsLb8eRXwLUk3UOwYLwJ+NsWDoFl+s8zS6GYC3QWSTqA4bLoB+DhARDwqaSWwjuLI0iWtjgwBPsfYkupmAt3qFre5Fri27VW4zNsSyuIdY5d5W0pZhMBl3pZSFiFwmbellEUIXOadrzp8ujeLELjMO18u854hLvO2lLIIgcu8c+YtwYxwmbellEcIXOZtCWXxzHOZt6WURQh2xojLvDPlMu8ZsnnHPFzmnSeXec8Ql3lnzGXeM8Nl3pZSFiFwmbellEUIXOZtKbVzeuUc4D6Ktw6HgTsi4rOSjgNWAIcBPwcujIhdkmZTDOs6EXgBOC8iNkzxIGjEh0gtjXa2BDuBMyLi5fKE+59I+jfgUxTDt1ZI+irFxLmvlN+3RsTxks4HvgCc1+oB1GjQOMgn1Vga7ZxeGcDEweKR8iuAM4APlpffBlxNEYIl5c8AdwA3S1J5P/vnc4wtoXZHrgxRvOQ5Hvgy8CtgW0RMvJBvHrB1NPA0QETskTRK8ZLp+RYPAMN57J5Y/bT1zCunRZwgaT5wJ/AH3T7wXhPohg92mbcl09EzLyK2SboHOBWYL2m43Bo0D9iaGL61UdIwMI9iB3nf+1oOLAc4eP4x4TJvS6Wdo0NvAHaXATgAeDfFzu49wPspjhBdBNxV3mRV+fv95d9/1HJ/AJd5W1rtbAkWALeV+wUNYGVE3C1pHbBC0ueBX1BMqaP8/o1yGvWLFGPaW3KZt6XUztGh/6aYRL3v5U9QjFzf9/JXgQ90sgiXeefLZd4zxGXe+XKZ9wxxmbellEUIXOadM5d5zwiXeVtKWYTAZd6WUhYhcJm3pZRFCFzmnS+Xec8Ql3nny2XeM8Rl3hlzmffMcJm3pZRFCGJszGXelkwWIQBqUxxt+ckiBC7ztpSyCAGNBnIILJFMQuAT7S2dPEKAYKiRehFWU3mEQBDDQ6lX0RMxd07qJViHuplA93XgT4GJY5sXR8QjkgTcBCwGdpSXT/2244AcHXIpefV0M4EO4G8j4o59rn82sKj8OoViINcUHwzSwITApeTV080EusksAW4vb/eApPmSFkTE5q5XWwGDVkpeh0/3TmsCXUT8VNLfANdK+kfgh8CyiNhJ0wS60sR0us373GfT8K1Duv13ZGPQSsld5l3adwKdpLcDVwLPALMohmhdAXyu3QduHr41b86CgRk65FLy6pnuBLqzIuL68uKdkv4Z+Ez5+8QEugnN0+kG3uCVkntLMOkEuonX+eXRoPcBa8ubrAIulbSCYod4tC77A4BLySuomwl0PyoDIuAR4BPl9VdTHB5dT3GI9CM9X3XGXEpePd1MoDtjkusHcEn3S6sml5JXj/+31WODVkruMm/r2KCVkrvM2zo2cKXkLvO2TrmUvHocgh5zKXn1OAQ95lLy6nEIesyl5NXjEPSYS8mrxyHoMZ8rXT0OQa/N9suhqnEIemx8nl8OVY1D0GM7Fh6cegnWIYegx1xKXj0OQY+5lLx6HIIeG7RScpd5W8cGrZTcZd5NyjPLHgI2RcR7JB0HrAAOo5hEcWFE7JI0G7gdOBF4ATgvIjb0fOWZcil59XSyJbgceAyYmI/yBeDGiFgh6avAUopBW0uBrRFxvKTzy+ud18M1Z23wSsld5g2ApGOAvwSuBT5Vnlx/BvDB8iq3AVdThGBJ+TPAHcDNklSedjnwXEpePe1uCb4E/B0wcRD8MGBbREx8ZHJiwBY0Dd+KiD2SRsvrP9+LBefOpeTV087IlfcAWyLi55JO79UDD+oEujqMLRw07WwJTgPeK2kxMIdin+AmYL6k4XJr0Dxga2L41kZJw8A8ih3kvQzqBLpBKyV3mTcQEVdSjFyk3BJ8JiI+JOlfgPdTHCG6CLirvMmq8vf7y7//qC77A8DAlZK7zLu1K4AVkj4P/AK4pbz8FuAbktYDLwLnd7fEahm4UnKXee8tIu4F7i1/fgI4eT/XeRX4QA/WVkkuJa8ev2PcYy4lrx6HoNcGpHGnThyCHnMpefU4BD3mUvLqcQh6zCfaV49D0GsD0sdcJw5Bj8VQI/USrEMOQa/56FDlKIdPNEh6DnhqBh7qcAbn06z+t3TmTRGx3waVLEIwUyQ9FBEnpV5HL/jf0jt+AWu15xBY7dUtBMtTL6CH/G/pkVrtE5jtT922BGavU4sQSDpL0i8lrZe0LPV6uiHpVklbJK1NvZZuSFoo6R5J6yQ9KunyZGsZ9JdD5dCwx4F3U0zFeBC4ICLWJV3YNEn6E+Bl4PaIeHvq9UyXpAXAgoh4WNLBFAPc3pfiv0sdtgQnA+sj4omI2EVxTvSSxGuatoi4j+K01UqLiM0R8XD583aKwW5Ht75Vf9QhBK/NQSo1z0iyDEg6Fngn8NMUj1+HEFjGJB0EfBf4ZES8lGINdQjBxBykCc0zkiwhSSMUAfhmRHwv1TrqEIIHgUWSjpM0i2IEzKrEa6q9cp7tLcBjEXFDyrUMfAjKCXmXAmsodr5WRsSjaVc1fZK+TTHY7PclbZS0NPWapuk04ELgDEmPlF+LUyxk4A+Rmk1l4LcEZlNxCKz2HAKrPYfAas8hsNpzCKz2HAKrPYfAau//Ab8lbhK69RVqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X3 = np.array([np.convolve(np.repeat(8 * [[0] * 20 + [1] * 20], 1), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "               np.convolve(np.append([0] * 80, np.repeat(8 * [[0] * 20 + [1] * 20], 1)), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "              [1] * len(Y)]).T\n",
    "model_3 = GeneralLinearModel(X3)\n",
    "\n",
    "plt.imshow(model_3.X, aspect=0.01, interpolation='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "task > rest: 355.26 (z = 18.88)\n"
     ]
    }
   ],
   "source": [
    "model_3.fit(Y)\n",
    "con = model_3.contrast(np.array([1, 1, 0]))\n",
    "\n",
    "print(f\"task > rest: {con.effect[0,0]:.2f} (z = {con.z_score()[0]:.2f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 4: Two regressors, unique variance of model 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1b899be7d30>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAD8CAYAAADOpsDvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQgUlEQVR4nO3dfYwc9X3H8ffnHnzGBuzwWPOQ4ILTlqDGoYgHUbUURAMuwqAmxC4iDrLkUIFK1KYFGqmhaZBAIiGp0tI6gmJHKY5LgrCQW0oBl0Qqj4ZQYwpxwDw4BheMjanB9t19+8f8jixw59u9nfXM7e/zkk63OzO7+z3Bh5md3fmgiMAsZz1VD2BWNYfAsucQWPYcAsueQ2DZcwgsex0LgaRzJD0raYOkqzv1OmbtUic+J5DUCzwHnA28AjwKLIyI9aW/mFmbOrUnOBnYEBHPR8RuYAUwv0OvZdaWvg4975HAyw33XwFOGWvjKb3TYr/+GR0axdpx2JxtVY9Qii2b9rB966BGW9epEIxL0hJgCcDUvgM57ZhFVY1ie/End91d9QiluHL+z8dc16nDoU3A0Q33j0rL3hMRSyPipIg4aUrvtA6NYTa+ToXgUWCOpNmSpgALgFUdei2ztnTkcCgiBiVdAdwD9AK3RsTTnXgts3Z17D1BRKwGVnfq+c3K4k+MLXsOgWXPIbDsOQSWPYfAslfZJ8aNolcMT59a9Ril6Pm/d6sewVpUixDsOaCXTWd1x3eHjl7lEEw2tQhB/4zdHHHuS1WPUY5VA1VPYC2qRQhmDWznK7O744ta1/OHVY9gLapFCKYrOGVgT9VjWKZqEYIexID6qx7DMlWLEAwT7BzeXfUYlqlahGAwhnndIbCK1CMEiG3DtRjFMlSLf/OGETuGp1Q9hmWqFiHYFf1s3HNo1WNYpmoRgreHBnhw+8erHsMy1VYIJG0EdgBDwGBEnCTpIOAHwDHARuCiiHhzb8+zY/s01twzt51RauNYtlQ9grWojD3B70XE6w33rwbui4jrU/3i1cBVe3uCKW8N89F73ilhFLPWdeJwaD5wRrq9DFjDOCHQO7vpX/dCB0apwGGHVD2BtajdEATw75IC+MeIWAocHhGb0/pXgcNHe+D7yreYxtC27W2OUg+9DsGk024IfjsiNkk6DLhX0v80royISAH5kBSYpQAH6qBAozbkmXVcWyGIiE3p9xZJd1IU8b4maVZEbJY0C8Z/p6ieHnqmuYXOqjHhEEiaDvRExI50+/eBr1E0zS0Crk+/7xr3yXp6kENgFWlnT3A4cKeKw5g+4J8j4t8kPQqslLQYeBG4aNxn6hGa2h0Xo/j/Cj35TDgEEfE88MlRlr8BnNXSk0nQX4vP7SxDtfg3b3hqH28ff1jVY5Ri+nNvVD2CtagWIdg9E146vzsOJH7jxqonsFbVIgSHTt/BZaetqXqMUvwnv1n1CNaiWoTgoN6dXDzjiarHKIVDMPnUIgRT1MtRfftXPYZlqhYhCIKhGK56DMtULUIwRPDmsL9FatWoRQgGA7Z6R2AVqUUIhhDbfI2xVaQWIRgOsXO4O742YZNPLUKwdWg6t78x5v/w3qyjahECX2NsVapFCHp3wcyfVT2F5aoWIejbsYdDfvyLqscoRfjbsJNOLf6JxZ7dDL28qeoxStHzqx+regRrUS1CQEAMDlY9hWVq3BBIuhU4D9gSESekZaMWbKm4zOzbwDxgJ/CFiFg77hQS6vfnBFaNZvYEtwHfAZY3LBurYOtcYE76OQW4Of3eK/X00LP/9NYmNyvJuCGIiAclHfOBxWMVbM0HlkdEAA9JmjnSPLHXF/E1xlahib4nGKtg60jg5YbtXknL9h4C9cCAD4esGm2/Md5bwdbevK+Bru9Aoq+33VHMJmSiIRirYGsTcHTDdkelZR/S2EA3Y+qvuIHOKtMzwceNFGzB+wu2VgGfV+FUYPu47wcAUFG70g0/Nuk0c4r0doo3wYdIegX4KkW73GgFW6spTo9uoDhFemkHZjYrVTNnhxaOsepDBVvprNDl7Q5lti9N9HDIrGs4BJY9h8Cy5xBY9hwCy55DYNlzCCx7DoFlzyGw7DkElj2HwLLnEFj2HALLnkNg2XMILHsOgWXPIbDsjRsCSbdK2iJpXcOyayVtkvRk+pnXsO4aSRskPSvp050a3KwszewJbgPOGWX5TRExN/2sBpB0PLAA+ER6zN9LcpeK1dq4IYiIB4GtTT7ffGBFROyKiBcoLrg/uY35zDqunfcEV0h6Kh0ufSQtG6uB7kMkLZH0mKTHdg/tbGMMs/ZMNAQ3A8cCcykqFr/R6hNExNKIOCkiTprSO22CY5i1b0IhiIjXImIoIoaB7/LLQ56mG+jM6mJCIUjViyMuBEbOHK0CFkgakDSboqL9kfZGNOusiTbQnSFpLkUT+UbgiwAR8bSklcB6YBC4PCKGOjK5WUkm2kB3y162vw64rp2hzPYlf2Js2XMILHsOgWXPIbDsOQSWPYfAsucQWPYcAsueQ2DZcwgsew6BZc8hsOw5BJY9h8Cy5xBY9hwCy14z5VtHS3pA0npJT0u6Mi0/SNK9kn6Wfn8kLZekv00FXE9JOrHTf4RZO5rZEwwCfxYRxwOnApenkq2rgfsiYg5wX7oPcC7FtcVzgCUUzRRmtdVM+dbmiFibbu8AnqHoEpoPLEubLQMuSLfnA8uj8BAw8wMX5pvVSkvvCSQdA3wKeBg4PCI2p1WvAoen200XcJnVQdMhkLQ/8EPgSxHxVuO6iAiK5ommuYHO6qKpEEjqpwjA9yPiR2nxayOHOen3lrS8qQIuN9BZXTRzdkgUFSvPRMQ3G1atAhal24uAuxqWfz6dJToV2N5w2GRWO+P2DgGnA5cA/y3pybTsL4HrgZWSFgMvAheldauBeRSN1DuBS8sc2KxszZRv/QTQGKvPGmX7AC5vcy6zfcafGFv2HALLnkNg2XMILHsOgWXPIbDsOQSWPYfAsucQWPYcAsueQ2DZcwgsew6BZc8hsOw5BJY9h8Cy5xBY9tppoLtW0iZJT6afeQ2PuSY10D0r6dOd/APM2tXMNcYjDXRrJR0APC7p3rTupoi4sXHj1E63APgEcATwH5I+HhFDZQ5uVpZ2GujGMh9YERG7IuIFigvuTy5jWLNOaKeBDuCKVLp760ghL0020Ll8y+qinQa6m4FjgbnAZuAbrbywy7esLibcQBcRr0XEUEQMA9/ll4c8TTXQmdXFhBvoPtA0fSGwLt1eBSyQNCBpNkVF+yPljWxWrnYa6BZKmktRxLsR+CJARDwtaSWwnuLM0uU+M2R11k4D3eq9POY64Lo25jLbZ/yJsWXPIbDsOQSWPYfAsucQWPYcAsueQ2DZcwgsew6BZc8hsOw5BJY9h8Cy5xBY9hwCy55DYNlzCCx7zVxeOVXSI5J+msq3/jotny3p4VSy9QNJU9LygXR/Q1p/TIf/BrO2NLMn2AWcGRGfpGiWOEfSqcANFOVbxwFvAovT9ouBN9Pym9J2ZrXVTPlWRMTb6W5/+gngTOCOtHwZcEG6PT/dJ60/K12sb1ZLzVau9KaL7LcA9wI/B7ZFxGDapLFg673yrbR+O3BwiTOblaqpEKR+obkUHUInA7/e7gu7gc7qoqWzQxGxDXgAOA2YKWmkraKxYOu98q20fgbwxijP5QY6q4Vmzg4dKmlmur0fcDZFKe8DwGfSZouAu9LtVek+af39ERElzmxWqmbKt2YByyT1UoRmZUTcLWk9sELS14EnKFrqSL+/J2kDsJWipt2stpop33qKoon6g8ufZ5TK9Yh4F/hsKdOZ7QP+xNiy5xBY9hwCy55DYNlzCCx7DoFlzyGw7DkElj2HwLLnEFj2HALLnkNg2XMILHsOgWXPIbDsOQSWvWauLLMWDE+fWvUI1qJxQyBpKvAgMJC2vyMivirpNuB3KSpVAL4QEU+mjqFvA/OAnWn52k4MX0ebzppR9QjWomb2BCMNdG9L6gd+Iulf07o/j4g7PrD9ucCc9HMKcHP6nYUjzn2p6hGsRc1cYxzAaA10Y5kPLE+Pe0jSTEmzImJz29NOAl+ZfXfVI5Rq5/BA1SN0XFPvCVLTxOPAccDfRcTDkv4YuE7SXwH3AVdHxC4aGuiSkXa6zR94ziXAEoCpfQe2+3fUxikDe6oeoVQPvOMQAEUDHTA39Q/dKekE4BrgVWAKsBS4Cvhasy8cEUvT45gxdVbX9BINqL/qEaxFLZ0diohtkh4AzomIG9PiXZL+Cfhyuv9eA13S2E7X9XYO7656hJJ5T4CkQ4E9KQAjDXQ3jBznp7NBFwDr0kNWAVdIWkHxhnh7Lu8HAF7vuhB0v3Ya6O5PARHwJHBZ2n41xenRDRSnSC8tfeoa2zbsj14mm3Ya6M4cY/sALm9/tMlpx/CUqkewFvk/WyXbuOfQqkco1cG9b4+/0STnEJTsb346r+oRSvWtE1dWPULHOQQlm3H3/lWPUK4Tqx6g8xyCkh3y419UPYK1yCEo2dDL2Xwk0jUcgpLF4OD4G1mtOAQlU79PkU42DkHJevafXvUI1iKHoGSa2v3ftek2DkHZBnw4NNk4BCUbnuHDocnGISjZzqMPqHoEa5FDULKXzu+a64Oy4RCU7LLT1lQ9grXIISjZxTOeqHqEUq3bfXDVI3ScQ1Cyo/q66wt06zK4UK7pEKQryx4DNkXEeZJmAyuAgymaKC6JiN2SBoDlwG8BbwCfi4iNpU9eU0MxXPUI1qJW9gRXAs8AI/0oNwA3RcQKSf8ALKYo2loMvBkRx0lakLb7XIkz19qbw+9UPULJ9qt6gI5rtnfoKOAPgOuAP00X158J/FHaZBlwLUUI5qfbAHcA35GkdNll19vqHcGk0+ye4FvAXwAjJ8EPBrZFxMhXJkcKtqChfCsiBiVtT9u/XsbAdbfN1xhPOs1UrpwHbImIxyWdUdYLd2sDXQ61hd2mmT3B6cD5kuYBUyneE3wbmCmpL+0NGgu2Rsq3XpHUB8ygeIP8Pt3aQHf7G93VPXzhQd1fKN5M5co1FJWLpD3BlyPiYkn/AnyG4gzRIuCu9JBV6f5/pfX35/J+AGDNPXOrHqFUFy50CPbmKmCFpK8DTwC3pOW3AN+TtAHYCixob8TJ5aP3dNnZoYVVD9B5rXaRrgHWpNvPAyePss27wGdLmG1S6l/3QtUjWIv8iXHJhrZtH38jqxWHoGxS1RNYixyCkvVMm1b1CNYih6BkcggmHYegZL7QfvJxCMrW11v1BNYih6Bk0dtT9QjWIoegbD47NOmoDt9okPS/wIv74KUOoXu+zeq/pTUfi4hR/w8qtQjBviLpsYg4qeo5yuC/pTw+gLXsOQSWvdxCsLTqAUrkv6UkWb0nMBtNbnsCsw/JIgSSzpH0rKQNkq6uep52SLpV0hZJ66qepR2Sjpb0gKT1kp6WdGVls3T74VAqDXsOOJuiFeNRYGFErK90sAmS9DvA28DyiDih6nkmStIsYFZErJV0AEWB2wVV/HPJYU9wMrAhIp6PiN0U10TPr3imCYuIBykuW53UImJzRKxNt3dQFLsdufdHdUYOIXivBylp7EiyGpB0DPAp4OEqXj+HEFiNSdof+CHwpYh4q4oZcgjBSA/SiMaOJKuQpH6KAHw/In5U1Rw5hOBRYI6k2ZKmUFTArKp4puylPttbgGci4ptVztL1IUgNeVcA91C8+VoZEU9XO9XESbqdotjs1yS9Imlx1TNN0OnAJcCZkp5MP/OqGKTrT5Gajafr9wRm43EILHsOgWXPIbDsOQSWPYfAsucQWPYcAsve/wP4IwIBXG7zvAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X4 = np.array([np.convolve(np.repeat(2 * [[0] * 20 + [1] * 20], 1), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "               np.convolve(np.append([0] * 320, np.repeat(2 * [[0] * 20 + [1] * 20], 1)), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "              [1] * len(Y)]).T\n",
    "model_4 = GeneralLinearModel(X4)\n",
    "\n",
    "plt.imshow(model_4.X, aspect=0.01, interpolation='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "task > rest: 409.02 (z = 8.31)\n"
     ]
    }
   ],
   "source": [
    "model_4.fit(Y)\n",
    "con = model_4.contrast(np.array([1, 1, 0]))\n",
    "\n",
    "print(f\"task > rest: {con.effect[0,0]:.2f} (z = {con.z_score()[0]:.2f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 5: two regressors, shared variance of model 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1b899c506a0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAD8CAYAAADOpsDvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQwUlEQVR4nO3de4xc9XnG8e+zF9tc18GAawypUXBaEaQ4iYuLqNoUSgVuFIOUEGjErVZIJBBERC2XqgpSikSkBEKVKpUTKKYKGIsE4VpuKeEiGgkI5lJqIEk3xgS7DgZsL6YmXnv37R/nLBmWvczu7MzvF37PR1rtzJkzc16LeThnzux5X0UEZiXrSl2AWWoOgRXPIbDiOQRWPIfAiucQWPHaFgJJZ0r6maR+Sde0aztmrVI7vieQ1A38HDgD2Ao8CZwfES/M+MbMWtSuPcHJQH9EbI6IQWANsKJN2zJrSU+bXnch8ErD/a3AsvFWntV9cBzU29emUqwVRy/enbqEGbFj234Gdh7QWI+1KwSTknQpcCnAnJ7DOWXRRalKsQlccd/61CXMiCtX/GLcx9p1OLQNOK7h/rH1sndExKqIWBoRS2d1H9ymMswm164QPAkslnS8pFnAecC6Nm3LrCVtORyKiAOSLgfuB7qB2yLi+XZsy6xVbftMEBEbgA3ten2zmeJvjK14DoEVzyGw4jkEVjyHwIrnEFjxHAIrnkNgxXMIrHgOgRXPIbDiOQRWPIfAiucQWPEcAiueQ2DFcwiseC1dWSZpC7AHGAIORMRSSUcAdwOLgC3AuRGxq7UyzdpnJvYEfxoRSyJiaX3/GuDBiFgMPFjfN8tWOw6HVgCr69urgbPbsA2zGdNqCAL4D0lP1c20AOZHxPb69q+A+WM9UdKlkjZK2jg4tLfFMsymr9VuE38UEdskHQ08IOmnjQ9GREgas+NvRKwCVgH0zVng6YGWTEt7gojYVv/eAdxL1Yj3VUkLAOrfO1ot0qydph0CSYdIOmzkNvDnwCaqTnMjjUUvAu5rtUizdmrlcGg+cK+kkde5MyL+XdKTwFpJK4GXgXNbL9OsfaYdgojYDHx0jOVvAKe3UpRZJ/kbYyueQ2DFcwiseA6BFc8hsOI5BFY8h8CK5xBY8RwCK55DYMVzCKx4DoEVzyGw4jkEVjyHwIrnEFjxJg2BpNsk7ZC0qWHZEZIekPQ/9e8P1Msl6R8k9Ut6TtLH21m82UxoZk9wO3DmqGXjNdg6C1hc/1wKfGdmyjRrn0lDEBGPAjtHLR6vwdYK4I6oPA7MHek8YZar6X4mGK/B1kLglYb1ttbLzLLVavOtCRtsTaTuWHcpwJxZfQwfMqfVUsymZboheFXSgojYPqrB1jbguIb1jq2XvUdjB7qDfue42HZ63zRLMWvNdEMw0mDrRt7dYGsdcLmkNcAyYKDhsGlcvX2DHHPWL6dZillrJg2BpLuATwJHStoKfJXqzT9Wg60NwHKgH9gLXNJMEQtmD/C3x6+fcvHWfnuHZ6cuoe0mDUFEnD/OQ+9psBURAVw21SIOUbBs9v6pPs064OG3HYKO6ELMVm/qMqxQWYRgmGDv8GDqMmxM3hN0xIEY5nWHwBLJIwSI3cNZlGIFyuKdN4zYMzwrdRlWqCxCsC962bL/qNRl2Bjmdb+VuoS2yyIE2/f28bX/Wp66DBvDtz6+NnUJbZdFCLoHuuhbf2jqMmwsBVwRkkcI9g1z+Oa3U5dhhcoiBHp7kN5NL6UuwwqVRQhiaIih3QOpy7BCZRECAKopmGYdl0UI1NVF18EHpy7DCpVFCOjqQg6BJZJJCITmvP//UMvylEcIJOjNoxQrTxbvvOE5Pbx14tGpy7BCNXN55W3Ap4AdEXFSvex64AvAa/Vq10XEhvqxa4GVwBBwRUTcP9k2BufCLz895YYVZjOimT3B7cC3gTtGLb85Ir7RuEDSicB5wEeAY4AfSfpwRAxNtIGjDtnDl055pNmazWZUM9cYPyppUZOvtwJYExH7gJck9QMnA49N9KQjuvfy+b5nmtyEddKmwXmpS2i7Vj4TXC7pQmAj8JWI2EXVbe7xhnXG7UDX2Hzrgwt7OLbHf0CXo00FXPA33RB8B/gaEPXvbwJ/NZUXaGy+9YmPzo6hGJ5mKWatmVYIIuLVkduSvguMNA1qugNdoyGCXcP+K9I8HZS6gLabVghGWjDWd88BRmYXrAPulHQT1QfjxcBPJnu9AwE7vSOwRKbbge6TkpZQHQ5tAb4IEBHPS1oLvAAcAC6b7MwQwBBit68xtkSm24Hu1gnWvwG4YSpFDIeKaPdnecriG+OdQ4dw1xvLUpdhYzjniKdTl9B2WYRgz8DBPHL/ktRl2BjOOd8h6IhZbw7zwft9dihL47Vjfh/JIgS+xthSyiIEMTzE8Fv/l7oMK1QWISAgDhxIXYUVKo8QSKjX3xNYGlmEQF1ddB16SOoyrFBZhMDXGFtKeYRAXTDbh0OWRh4h6O5iuM+HQ5ZGFiEYPLybbX/2gdRlWKGyCEFv3yDzl7+SugwrVBYhWDB7gL87/l9Tl2FjKOGve7MIgYd558vDvDvEw7wtpWauLDuOqufQfKoryVZFxC2SjgDuBhZRXV12bkTskiTgFmA5sBe4OCIm/HtcD/POmfcEUF0m+ZWIeFrSYcBTkh4ALgYejIgbJV0DXANcDZxFdW3xYmAZVWeKCa+Y8TBvS6mZyyu3A9vr23skvUjVS2gF1bXHAKuBR6hCsAK4IyICeFzS3FEX5r+Hh3lbSlN659Wd6D4GPAHMb3hj/4rqcAmqgDSe7xxpwDVuCDzM21JqOgSSDgV+AHw5It5Uw3iliAhJU+qo29iBbt4xszzMO1Me5l2T1EsVgO9HxA/rxa+OHOZIWgDsqJc31YCrsQPdnA8tDA/zzpOHeQP12Z5bgRcj4qaGh9YBFwE31r/va1h+uaQ1VB+IByb6PAAe5p01D/MG4FTgAuC/JT1bL7uO6s2/VtJK4GXg3PqxDVSnR/upTpFeMmkRe/Zz5H/+79QqN5shzZwd+jEw3nzV08dYP4DLplJE7B9k6JVJW5aatUUe5yV9jbEllEcIwMO8LZksQuBh3pZSFiHwMG9LKZMQ+EJ7SyePEHiYtyWUxTvPw7wtpSxC4GHellIWIfAwb0spixB4mHe+PMy7Q2ap28O8M+Vh3h0SBB7mbalkEQIP886Zh3l3hId5W0pZhMDDvC2lLELgYd6WUhYh8DDvfHmYNxN2oLse+ALwWr3qdRGxoX7OtcBKYAi4IiLun2gbHuadLw/zrozXgQ7g5oj4RuPKkk4EzgM+AhwD/EjShyNiaLwNeJh3xjzMe8IOdONZAayJiH3AS5L6gZOBx8Z7god5W0qtdKA7laq1yoXARqq9xS6qgDze8LSRDnSjX+ud5ltzOJih3QPTqd+sZV3Nrji6Ax1Vo90PAUuo9hTfnMqGI2JVRCyNiKW9zK6uKfBPfj8FmHYHuoh4teHx7wLr67tNdaB71+v7GmNLaNod6EZ1mj4H2FTfXgfcKekmqg/Gi4GfTLIRNMtfllkarXSgO1/SEqrTpluALwJExPOS1gIvUJ1ZumyiM0OArzG2pFrpQLdhgufcANzQdBUe5m0JZfGNsYd5W0pZhMDDvC2lLELgYd6WUhYh8DDvfJXw171ZhMDDvPPlYd4d4mHellIWIfAw75x5T9ARHuZtKeURAg/ztoSyeOd5mLellEUI9kWvh3lnysO8O2T73j48zDtPHubdIR7mnTEP8+4MD/O2lLIIgYd5W0pZhMDDvC2lZi6vnAM8SvXVYQ9wT0R8VdLxwBpgHvAUcEFEDEqaTdWs6xPAG8DnImLLJBtBvT5Famk0syfYB5wWEW/VF9z/WNK/AVdRNd9aI+mfqDrOfaf+vSsiTpB0HvB14HMTbUBdXXQd6otqLI1mLq8MYORkcW/9E8BpwF/Wy1cD11OFYEV9G+Ae4NuSVL/O2HyNsSXUbMuVbqpDnhOAfwR+AeyOiJED+cYGWwuBVwAi4oCkAapDptcn2AD05PHxxMrT1Duv7haxRNJc4F7g91vd8Ls60PUc5mHelsyU3nkRsVvSw8ApwFxJPfXeoLHB1kjzra2SeoA+qg/Io19rFbAK4LC5x4aHeVsqzZwdOgrYXwfgIOAMqg+7DwOfoTpDdBFwX/2UdfX9x+rHH5rw8wAe5m1pNbMnWACsrj8XdAFrI2K9pBeANZL+HniGqksd9e9/qbtR76Rq0z4hD/O2lJo5O/QcVSfq0cs3U7VcH73818Bnp1KEh3nny8O8O8TDvPPlYd4d4mHellIWIfAw75x5mHdHeJi3pZRFCDzM21LKIgQe5m0pZRECD/POl4d5d4iHeefLw7w7xMO8M+Zh3p3hYd6WUhYhiKEhD/O2ZLIIAVDM4GjLTxYh8DBvSymLENDVhRwCSySTEPhCe0snjxAg6O5KXYQVKo8QCKKnO3UVVqhWOtDdDvwJMHJu8+KIeFaSgFuA5cDeevnkXzv67JAl0koHOoC/joh7Rq1/FrC4/llG1ZBrkj8MkkNgyUx6IB6VsTrQjWcFcEf9vMepWrMsaL1Us/Zo6tOopG5JzwI7gAci4on6oRskPSfp5roRLzR0oKs1dqdrfM1LJW2UtHFwaO/0/wVmLWoqBBExFBFLqJpsnSzpJOBaqk50fwAcAVw9lQ1HxKqIWBoRS2d1+zsCS2dK5yUjYjdV060zI2J7fcizD/hnftN+ZaQD3YjG7nRm2Zk0BJKOqnuQ0tCB7qcjx/n12aCzgU31U9YBF6ryh8BARGxvQ+1mM6KVDnQP1S0aBTwLfKlefwPV6dF+qlOkl8x41WYzqJUOdKeNs34Al7Vemlln+G8VrHgOgRXPIbDiOQRWPIfAiucQWPEcAiueQ2DFcwiseA6BFc8hsOI5BFY8h8CK5xBY8RwCK55DYMVrOgR1x4lnJK2v7x8v6QlJ/ZLuljSrXj67vt9fP76oTbWbzYip7AmuBF5suP914OaIOAHYBaysl68EdtXLb67XM8tWs32HjgX+AvhefV/AacBI97nVVBfbQ9V8a3V9+x7g9Hp9syw1uyf4FvA3wMjc+XnA7og4UN9vbLD1TvOt+vGBen2zLDXTcuVTwI6IeGomN+wOdJaLZlqunAp8WtJyYA5wOFXX6bmSeur/2zc22BppvrVVUg/QB7wx+kUjYhWwCqBvzoKJepuatVUzDXmvjYhjI2IRcB7wUER8nqoT3Wfq1S4C7qtvr6vvUz/+UN2GxSxLrXxPcDVwlaR+qmP+W+vltwLz6uVXAde0VqJZe01pUk1EPAI8Ut/ezG/6jzau82vgszNQm1lH+BtjK55DYMVzCKx4DoEVzyGw4jkEVjyHwIrnEFjxlMNfNEh6DXi5A5s6Eni9A9vpBP9bpuZ3I+KosR7IIgSdImljRCxNXcdM8L9l5vhwyIrnEFjxSgvBqtQFzCD/W2ZIUZ8JzMZS2p7A7D2KCIGkMyX9rO6F9Ft9kY+k2yTtkLQpdS2tkHScpIclvSDpeUlXJqvl/X44JKkb+DlwBlVXjCeB8yPihaSFTZOkPwbeAu6IiJNS1zNdkhYACyLiaUmHAU8BZ6f471LCnuBkoD8iNkfEILCGqjfSb6WIeBTYmbqOVkXE9oh4ur69h6qx28KJn9UeJYTgnT5ItcYeSZaBulXnx4AnUmy/hBBYxiQdCvwA+HJEvJmihhJCMNIHaURjjyRLSFIvVQC+HxE/TFVHCSF4Elhcd9GeRdU7aV3imopX96e9FXgxIm5KWcv7PgR1h7zLgfupPnytjYjn01Y1fZLuAh4Dfk/SVkkrJ3tOpk4FLgBOk/Rs/bM8RSHv+1OkZpN53+8JzCbjEFjxHAIrnkNgxXMIrHgOgRXPIbDiOQRWvP8HcoOUaEfqGjwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X5 = np.array([np.convolve(np.append([0] * 80, np.repeat(6 * [[0] * 20 + [1] * 20], 1)), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "               np.convolve(np.append([0] * 80, np.repeat(6 * [[0] * 20 + [1] * 20], 1)), hrf(np.array(range(len(Y)))))[:len(Y)],\n",
    "              [1] * len(Y)]).T\n",
    "model_5 = GeneralLinearModel(X5)\n",
    "\n",
    "plt.imshow(model_5.X, aspect=0.01, interpolation='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "task > rest: 256.72 (z = 12.68)\n"
     ]
    }
   ],
   "source": [
    "model_5.fit(Y)\n",
    "con = model_5.contrast(np.array([1, 1, 0]))\n",
    "\n",
    "print(f\"task > rest: {con.effect[0,0]:.2f} (z = {con.z_score()[0]:.2f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- testing the effects in model 1 (1 regressor) and model 2 (2 regressors, 100 colinearity) is statistically identical\n",
    "- effects in model 3 (2 regressors, 50% colinearity) are slightly worse\n",
    "- but effects in models capturing only the unique variance (model 4) or the shared variance (model 5) are significantly weaker"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "**Testing a union contrast of two colinear regressors seems to capture their shared variance in addition to their unique variances**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}