{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conditional Average Treatment Effects (CATE) with DoWhy and EconML\n",
    "\n",
    "This is an experimental feature where we use [EconML](https://github.com/microsoft/econml) methods from DoWhy. Using EconML allows CATE estimation using different methods. \n",
    "\n",
    "All four steps of causal inference in DoWhy remain the same: model, identify, estimate, and refute. The key difference is that we now call econml methods in the estimation step. There is also a simpler example using linear regression to understand the intuition behind CATE estimators. \n",
    "\n",
    "All datasets are generated using linear structural equations.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import logging\n",
    "\n",
    "import dowhy\n",
    "from dowhy import CausalModel\n",
    "import dowhy.datasets\n",
    "\n",
    "import econml\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "BETA = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         X0        X1   Z0        Z1        W0        W1 W2 W3         v0  \\\n",
      "0  1.898922 -1.493600  1.0  0.093164  0.097720  0.440277  3  2  22.528813   \n",
      "1  0.885180 -2.027130  0.0  0.790633  0.203576 -0.210274  2  2  16.619325   \n",
      "2  0.379053 -0.216526  0.0  0.049098 -0.453768 -0.439534  0  3  -1.908925   \n",
      "3  1.384133  1.860552  0.0  0.224985  1.216297  0.075126  1  0   7.027797   \n",
      "4  0.803883 -1.275331  1.0  0.511890  0.302380  1.060172  2  0  24.172841   \n",
      "\n",
      "            y  \n",
      "0  357.902396  \n",
      "1  199.083575  \n",
      "2  -21.889854  \n",
      "3  117.426939  \n",
      "4  296.879249  \n",
      "True causal estimate is 11.13817451301838\n"
     ]
    }
   ],
   "source": [
    "data = dowhy.datasets.linear_dataset(BETA, num_common_causes=4, num_samples=10000,\n",
    "                                    num_instruments=2, num_effect_modifiers=2,\n",
    "                                     num_treatments=1,\n",
    "                                    treatment_is_binary=False,\n",
    "                                    num_discrete_common_causes=2,\n",
    "                                    num_discrete_effect_modifiers=0,\n",
    "                                    one_hot_encode=False)\n",
    "df=data['df']\n",
    "print(df.head())\n",
    "print(\"True causal estimate is\", data[\"ate\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = CausalModel(data=data[\"df\"], \n",
    "                    treatment=data[\"treatment_name\"], outcome=data[\"outcome_name\"], \n",
    "                    graph=data[\"gml_graph\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAD7CAYAAABwkSklAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydeVxU1f//37Mww7DNgCACrqCI4AJSueAuiQumqZh+akxbsG/WZJhRH8uxsk9Y1ofSFqw0qk+aViKmqdiigJq5IbiguKSAINvAIDAwM6/fH/5mYmQRcGYuA+f5eMxDuXPmnNd93Tt37vucc9+HBwDEYDAYDAaDwWAwGAyGleBzLYDBYDAYDAaDwWAwGJ0LFogyGAwGg8FgMBgMBsOqsECUwWAwGAwGg8FgMBhWRci1AEbHpLa2lm7dukVERFqtltRqNRER6fV6Ki8vNylbU1ND1dXVd62zpeWIiJycnMjOzq7V5Xg8HslkMuPfMpmMeDweERG5uLiQQCBoUfsMRmejqqqKcnNzqaysjNRqNel0OqqoqCAiIldXVyK6/X1yd3cnLy8vEovFXMrtcDD/uYX5zy3Mf25h/nOLLfvPAtFOQGVlpfFVXl5OFRUVpNFoSK1Wk0ajoaqqKqqurqaamhq6desW1dbWUkVFBel0OlKpVKTX60mlUhkDyvpBpqE8ERnr6Aw4OjqSSCQiIiKJREL29vZE9E+wKpPJSCAQkFQqJaFQSM7OziQSicjR0dFY3lCHs7MzCYVCkslkJBaLycnJiWQyGTk7O5OTkxNJJBIud5XBMCE3N5f+/PNPOnPmDGVlZdG5c+fo2rVrxh+9luLu7k49e/akoKAgGjhwIA0cOJCGDx9Obm5uFlLeMWD+cwvzn1uY/9zC/OeWjug/j2XNbb9otVoqKytr9FVRUUEqlYrUarUxyKyoqKDy8nKTwLOsrKzZNuzs7MjJyYnEYjE5ODgYgyTDSGFzgRURkb29vTFQMgRaRGRShs/nk1QqNbZp6J0xcOcoZFO0tFxjo64tLVe/FwkAqVQq43sqlYoMXxdDoE50O9Cvq6szlmkqcDeM6BqCd7VaTVqtlsrLy0mv1zepUyAQkIuLC0mlUnJycjK+XF1djf93dnYmFxcXkslk5Orq2ujLMLLLYLQGlUpFu3fvpv3799PBgwfp0qVLxOfzydfXlwYNGkSBgYHUq1cv8vLyou7du1OXLl3IycnJeN4SkfE6VF5eTjdv3qSCggK6fv06/f3335SZmUlnzpyh69evE5/Pp6CgIBo3bhxNmjSJwsPDjZ08nRXmP7cw/7mF+c8tzH9u6Qz+s0DUSpSVlVFhYSEVFxdTUVERFRQUUGlpaZOBpmF4/U6EQiG5uroagw7DqJmjoyNJpVJycXExCU5kMplJ8GL4jEgkMgkOGdxiCHoNo9KGkWtDh4JarSaVSmXSyXDn3+Xl5aRSqaisrIwa+1o3F6S6urpS165dyd3dnTw8PKhr167k6elJTk5OHLjB4JrKykr6/vvvadu2bfT7778TABoxYgSNGzeOxowZQ8OHDzd2OpmL0tJSSktLowMHDtDBgwfpxIkT5OjoSFOmTKF58+ZRZGRki6bbdwSY/9zC/OcW5j+3MP+5pbP5zwLRNqLT6aiwsJByc3ONAWZBQQEVFRUZX4a/i4uLjdNXDbi5uZG7u3uzgUFjL8MoI4PRHIaAtCUvQ4dIUVFRg84PiURC7u7u1K1bN/Lw8DC+6v/t5eVFXl5e5OnpydHeMszFmTNn6JNPPqFvv/2WamtracqUKTRr1iyKjIxs0WwEc5Kfn09JSUm0fft2+u2338jT05OefPJJWrx4MXXv3t2qWqwF859bmP/cwvznFuY/t3RW/1kg2gjV1dV048YNys/PN/n38uXLxv9fu3aNtFqt8TP29vbGYNHb25u8vLwa/N/wd/fu3Y3PFzIY7QmNRkMlJSXGINVw/t/5/7KyMsrLyzOZ2iwSiahLly7Gc76xf/38/Kx+QWXcnaysLHr33Xfpu+++I19fX3ryySfpySefJHd3d66lEdHtH8VvvvmGPv74YyooKKB58+aRUqkkPz8/rqWZBeY/tzD/uYX5zy3Mf27p9P6jE1JQUIDDhw9j8+bN+M9//oPo6Gg8+OCD8Pf3h4ODA4jI+JJIJOjbty9Gjx6N+fPnIyYmBv/973+xZcsWpKam4vLly7h16xbXu8RgcEJFRQXOnTuH33//Hd988w3WrFkDhUKB2bNnY+TIkejRowdEIpHJd8rFxQVBQUGIjIzEc889h7Vr1+LHH3/E8ePHUVpayvUudSpyc3PxyCOPgMfjYejQoUhKSoJer+daVpNoNBokJCSgd+/eEIlEiImJQUVFBdey2gzzn1uY/9zC/OcW5j+3MP9v0yED0ZqaGmRkZOCHH37A2rVrsWTJEkybNg2BgYGQSCTGG2KhUIjevXtj/PjxWLRoEd544w1s3LgRv/zyCzIzM1FSUsL1rjAYNo9er8eNGzdw4sQJ7Ny5E5999hlee+01PProowgLC4O3tzd4PJ7xeymVShEcHIyZM2di6dKl+PDDD7Fz505kZ2ejrq6O693pEGi1WnzwwQdwdnaGn58ftm/f3q5/AO+ktrYWn376Kbp06QJvb298//33XEtqFcx/bmH+cwvzn1uY/9zC/DfFpgPR0tJSpKamIiEhAbGxsYiKikJgYCAEAoHxptbV1RWhoaGIiopCbGwsEhISkJKSgkuXLrGbWgajnaDRaHDp0iWkpKQgISEBSqUS0dHRCA8Ph6+vrzFQFQqF8PX1RXh4OBQKhfH7nJ+fz/Uu2AzXr1/H2LFjIRaLsXLlSlRVVXEtqc0UFRXhySefBI/Hw2OPPWYTvePMf25h/nML859bmP/cwvxviE08I1pYWEinTp2iU6dOUWZmJmVnZ1N2drYxsYpMJiN/f38KCAig/v37G1/9+vVrV4u2MhiMtlFZWUnZ2dl04cIFOn/+vPEacOHCBaqqqiIioi5dulBAQAAFBATQoEGDKDg4mIYMGcKeSa3H3r176dFHHyUPDw/asmULDRkyhGtJZuGXX36hhQsXklQqpZ9++okGDhzItaRGYf5zC/OfW5j/3ML85xbmfxOYNTy+R7RaLc6dO4ctW7YgNjYWkydPRrdu3Yyjmz169MC0adOwbNkybNiwAX/88QcKCwu5ls1gMDhCr9fj6tWr2LdvH9atW4clS5Zg4sSJcHd3N143+vTpg4cffhirVq1CUlISrly5wrVsTti0aROEQiHkcjkqKyu5lmN28vPzMWbMGMhkMvzxxx9cy2kA859bmP/cwvznFuY/tzD/m4bTQDQvLw9bt27F0qVLMXz4cGOiIKFQiEGDBkEul2Pt2rXYv38/iouLuZTKYDBsjOvXr+Pnn3/G6tWrMWfOHPTt29c4xdfV1RXjx4/HihUrsGvXrg6fJGnt2rXg8Xh47bXXbOpZlNZSU1ODqKgoiMViJCcncy3HCPOfW5j/3ML85xbmP7cw/5vHaoGoVqtFRkYGPv74Yzz66KPo3bs3iAgCgQAhISFYsmQJvvjiCxw7dgw1NTXWksVgMDoRFRUVSEtLw/r167Fo0SIEBASAx+OBz+cjKCgI0dHR+Oqrr3DhwgWupZqNL7/8EjweD/Hx8VxLsQo6nQ5PP/00JBIJDh48yLUc5j/HMP+5hfnPLcx/bmH+3x2LBqKXLl3CunXrMHXqVLi4uBiXboiIiMAbb7yB/fv3Q61WW1ICg8FgNEtRURGSk5MRGxuLUaNGwd7eHkQET09PzJ07F1999ZXNPgKwZ88eCIVCrFy5kmspVkWr1WL27NmQyWTIzs7mTAfzn/nPBcx/bmH+cwvzn1ta679ZkxXV1tZSamoq7d69m3bv3k3nz58nqVRKDz74IE2YMIHCwsIoKCiIBAKBuZpkMBgMs1JbW0vHjx+nQ4cOUUpKCh04cIBqa2spNDSUpk6dSlOnTqX77ruP+Hw+11KbpaCggIYMGUKTJk2ib775hms5Vkej0dDo0aNJp9PR4cOHSSQSWbV95j/zn0uY/9zC/OcW5j+3tMr/e418q6qq8P3332PWrFlwdnYGESEwMBDLly/Hb7/9htra2nttgsFgMDijsrISycnJeOaZZ9CzZ08QEbp27YpFixZhz5490Gq1XEtslClTpqBv3742kdLeUuTk5MDZ2RnLly+3etvMf+Y/1zD/uYX5zy3Mf25pqf9tDkRTU1OxYMECODs7QyAQYNKkSfj444/bXUbK69evG7Nn1n9t377dpNyKFSsalDl37pxFtb333nvGtnx8fCzali2xefNmoy9isbjVnz969Cgef/xx9O7dG/b29nB1dUVQUBBmzZqFTz75BDk5ORZQ3TRbtmzBkCFDjFM+iQiZmZlW1dBa7vUYdGROnz6NNWvWYNiwYcYpvC+88ALOnj3LtTQjycnJ4PF4SEtLu+e6PD0973ptnD17NogIubm5JtvvvK6+8847xvdOnjyJqVOnQiqVwsnJCRMnTjSL3jtJSEiAnZ2dVZ/7tQX/AWDXrl3o168fBALBPetsClv3/06CgoIavae488Xj8XDjxo0Gn7eG5/Vh/v+Dta459bF1/819/SktLcWnn36K8ePHw9XVFfb29ujbty/+9a9/4dSpU/es906Y/6b+6/V6pKWl4dlnn0W/fv0gEong4eGBsLAwfPPNN2ZPptQS/1sViNbU1CAhIQGDBg0CESE0NBQffvghCgoK7lmspTHcWMfGxjZbbuzYsfj888+tpOo2Q4YMYYFoI0ycOLFVQZBOp8NLL70EoVCI5cuX49y5c6ipqUFBQQH27duH8PBw4xeyrq7Ogsr/IS0tDTweD8uXL4darUZOTg66d+/e7gNRA609Bp2Nixcv4o033oCvry+ICOPGjcNPP/3EaWa8uro6DBgwAFFRUWarc8OGDSAiLF26tMF7Wq0WMpkMRIQvv/yywfslJSWQyWQms2OOHDkCiUSCRx55BPn5+SgqKsLTTz8NoVCIvXv3mk23QV9QUBBmz55t1nqbwhb8z8nJwfTp0zF48GC4uLhYNCjqCP7XJygoqEFHtoHq6moEBgaCiBATE2PynjU9rw/z/zbWvObUpyP4b87rz5NPPgmhUIj4+HjcuHEDt27dwsGDBxEYGAiBQNDksW0rzH9T/8+dOwciQnh4ODIyMlBdXY1Lly5h/vz5ICIsW7bMbLoN+u7mf4sC0draWqxbtw7du3eHWCzGokWLcPToUbMJtQYsELU9WhsE/fvf/wYRYcOGDY2+r9VqMWXKFKsGoi+88EKjPVW2AgtEW4ZOp8Mvv/yC6dOng8/nY9CgQdi6dSsnAemOHTvA5/PN2gN87do1EBECAgIavJeenm5cequxH98tW7ZgxowZxr91Oh2CgoLg5eWFqqoq43atVov+/fujR48eZs+c/tNPP4HH4+HSpUtmrbcx2rv/ADB//ny88847qKurg4+Pj8WDIlv3vz7NBUIKhQJEhIEDBzY4h63teX06u/9cXHPqY+v+m/P68+STTyI6OrpBuVOnToGI0K9fP7PpNsD8/8f/c+fOQSgUNliyTqPRoEuXLhCLxVb//b1rIPrHH38gMDAQEokECoXCZm+oWSBqe7QmCDp37hz4fD5CQ0ObLXfo0CGrBqKzZs0CEaG6utoq7ZkbFoi2nszMTMybNw98Ph/jx4+3+pTdGTNmIDw83Oz1DhgwAESEv//+22T7ypUrERMTAycnJ7i6ujZ4ZnbRokX45JNPjH///vvvICI8//zzDdpYtWoViAg//PCDWbXX1dXB29sbr7/+ulnrbYz27j8Ak5txawRFHcH/u7Fv3z7weDyIxWJkZGQ0eN/antens/vPxTWnPh3Bf3Nef5pCIpGAz+ebvQOX+d8y/4ODg0FEUKlUZtMN3N3/JtM+AqA1a9bQhAkTyNfXl86cOUMffvgh+fj4NPURBoMzNmzYQHq9nqKiopotN2LECAJAQqHQKrp0Op1V2mG0HwYOHEibN2+mv/76i6qqqigkJIQ+//xzq7StVqtp9+7d9Pjjj5u97smTJxMR0Z49e0y279mzh6ZPn04TJkygsrIyOnr0qMn7+/btM36WiOi3334jIqL77ruvQRuGbb/++qtZtQuFQpLL5bR582az1nsntuA/EZFEIjG7vuboCP43R2lpKS1cuJAA0Ntvv02DBw9uUMbantens/vPxTWnPh3Bf3Nefxrj1q1bVF1dTQMHDiQej2c+4cT8b4n/KpWKLl68SCEhISSVSs0nnO7uf6OBqE6no7lz59KqVatow4YNtHPnTurTp49ZhdkKSUlJxOPxjK+rV6/SI488QjKZjLp06UKRkZF06dKlBp8rKSmhmJgY8vPzI5FIRK6urjRlyhT6/fffm2zr/PnzNG3aNJJKpeTg4EDjx4+n9PR0kzIajYZWrlxJAQEB5ODgQG5ubjR9+nRKTk5uEPQUFRWRQqGg3r17k0gkIg8PD5o1axadOnWqyf3Lzs6muXPnUpcuXUy213+tXr2aiIi0Wq3J9jlz5rSq7fr7PXPmTJJKpeTo6EijR4+mtLS0lh2g/8/BgweJiBq9AbgbLTlWrT0PDOV37NhBRLdvQng8Hg0fPrxV7a5evdrY5qhRo4zb9+zZY9zu7u7eZp0GWnsMzHVuFRcXt+qctiWGDh1K6enptGzZMlq8eDEplUqLt3nkyBGqq6ujiRMnmr3uiIgIIiLau3evcVtpaSmdP3+eRo4c2ej7WVlZ5ODgYPL7cf78eSIi6t69e4M2DB2dFy5cMLv+iRMnUk5ODt24ccPsdRuwBf+5wtb9b47FixdTfn4+jRs3jmJiYqzadkvpzP5zdc2pj637b+nrz7Zt24iIaMWKFeaUbYT53zgVFRWUnp5ODz30EHXr1o2+/vprs2snuov/jQ2TLl++HBKJBKmpqWYdnuWSe52aO2PGDBARZsyYgUOHDqGyshIpKSmQSCS4//77TcreuHEDffr0gaenJ3bu3Iny8nJkZ2dj1qxZ4PF4DeofMmQIpFIpxo8fj7S0NKjVavz1118YPHgwRCIR/vjjD2PZp556ClKpFPv27UNVVRUKCgrw0ksvgYjw+++/G8vl5+ejV69e8PT0xK5du6BWq5GVlYWxY8fC3t4ehw4danT/xo4di99//x23bt3CkSNHIBAIUFRUhMmTJ4PP5zeacXbEiBH47rvv2tT2xYsXIZPJ4OPjg3379kGtVuP06dOYNGkSevfu3eJpoV5eXiAi/Pnnny0qb6C1x6o150H98ndOzW1tu46OjggLC2tQf2hoKLp06dJkuy3R2dpjYO5zq6XntC2zceNG8Hg8JCYmWrSdN954A3369LFI3dXV1ZBIJJBKpcap7Vu2bMH06dMB3E7GQkQYNmyY8TNr167Fc889Z1LPgw8+CCLCkSNHGrRx8eJFEBGGDh1qdv0VFRUQCAT48ccfzV63AVvw/06sNU3U1v1visTERBARZDIZrl271qLPWHtqLtC5/efqmlMfW/ffUtcfACgoKICnpyeeeuopi2gHmP+N8dZbbxkTeI4bNw6nT5+2iHagef8bBKIXLlyAnZ1do9mXbBlzBaI7d+402T5nzhwQEYqKiozbFi5cCCLC5s2bTcrW1NTA29sbEonEJNPwkCFDQEQ4fPiwSfnTp0+DiDBkyBDjtj59+mDkyJEN9Pn7+5vctD/++OMgIvzvf/8zKXfjxg2IxeIGz1Ia9m/37t0N6gaA/fv3g4jw7LPPmmxPS0tDz549TZ65bE3bUVFRjT6fkZeXB7FY3OpAtLVJtFp7rFpzHtQvf2cg2tp22xqItkRna4+Buc+tlp7Tts6LL74IT09PVFZWWqyNJ554ApMnT7ZY/RERESAiYyflokWLsH79euP7fn5+4PP5KCkpAXD7BvDnn382qaO5m8ILFy4YM7JbAh8fH3zwwQcWqRuwDf/vxJpBka37fydXrlyBi4sLiMikM/ZucBGIGtrtjP5zec2pj637b4nrT3FxMYKDg/HII49YfE1u5n9DNBoNzp07h2eeeQYCgQBvvvmmxfQ35X+Dqbm7du0id3d3WrhwYauGXds7AoGAiO7+zJ5OpzOWbYz777/f5O8ePXoQEVF+fr5x2/bt24mIaNq0aSZlxWIxTZw4kaqrq02Gz4mI7O3tadiwYSbbBg0aRN7e3pSRkWEczp48eTIdOnSIoqOj6ciRI8b9yc7OpnHjxhk/m5SURHw+nyIjI03q7NatGwUFBdHx48cpNze3wf498MADje73xIkTKSQkhL766isqKSkxbn/vvfdo6dKlJs9ctqZtw3x3w7QCA97e3uTv79+olsbw9vYmIqLi4uIWf4aobceKqGXngSXabS0t0dnaY2Duc6ul57St88orr1BhYSEdOnTIYm2UlJRQly5dLFa/4RwxnDN3Pn8yefJk0uv1lJKSQtXV1XT06NEGx1AmkxHR7WeC7sSwzVDG3Li7u5tcv8yNLfjPJbbuf330ej0tWLCAKioqaP78+TR//nyrtHsvdFb/ubzm1MfW/Tf39efWrVsUERFBgYGB9L///a/Ze29zwPxviEgkooCAAPr000/poYceopUrV9L+/fstor8p/xsEojdu3CBvb2/i85vMY2STODk5EdHt+dDNoVKpyMXFpcn373yIVyQSEdHtiyLR7Wc4y8vLyd7enpydnRt83tPTk4iICgoKTLYbnpu7k65duxIR0c2bN4mI6OOPP6avv/6aLl++TBMnTiQXFxeaPHmyMbCpr0Gv15NUKm3wjOeJEyeIiOjixYsN2nN0dGxy35ctW0ZVVVX0ySefENHtZyoOHjxITz31VJva1mg0pFaryd7e3nh8Gtv3ljB27FgiIjp9+nSLP9PWY0V09/PAUu22lpacr605BpY4t1pyTncEPDw8yMHBwaLPqFRVVVk0KYrhR2/v3r2UmZlJ9vb25OfnZ3y//nMqBw4coPvuu6/BcQ8ICCAiarSzIi8vj4ioVZ1QrcHR0ZEqKystUjeRbfjPJbbuf33effddSk1NpR49ehh/E9s7ndV/Lq859bF1/815/dFqtRQVFUU+Pj6UmJho8SCUiPl/N6ZPn05ERD///LMZVf9DU/43iDaDgoLo3LlzrR5Zau8YLjJnzpxpsoxGo6GcnBzq169fm9sRi8UklUqppqaG1Gp1g/cLCwuJ6PboUX3Ky8sbrc8QgBoCAh6PR3K5nPbv308qlYqSkpIIAM2aNYs++OADowaZTEZCoZDq6uoIt6dgN3iNHz++Vfv2yCOPUI8ePWj9+vWk0Wjo/fffp6efftokmGpN22KxmJydnammpqbRk7O0tLTF2hYvXkxCoZB++OGHZsu9/PLLxOfz6fz5820+VvdKW9rl8/lUW1vboKxKpbpnLa05BpY4t1pyTncEjh07RlVVVTRw4ECLteHq6kplZWUWq3/AgAHUo0cPOn78OH377bcNRtHHjx9PIpGI9u7dS3v27GnwvqEMEdHx48cbvGfYZqlkJ6WlpeTm5maRuolsw38usXX/DZw8eZJWrlxJPB6PEhMTrTKaZg46q/9cXnPqY+v+m/P6s3jxYtJoNLR161aTGXV9+/alI0eOWEQ/8795xGIxEbXu3rs1NOV/g0B0zpw5JJPJaNmyZRYRwhV+fn4UEBBAR44caXS0hoho69at5OHhcc83ig8//DAR3Z7mXB+NRkO//vorSSSSBidIZWUlZWRkmGzLzMyk/Px8GjJkCHl5eRHR7ekjhgxwdnZ29OCDDxqzk9Zvb9asWaTVahtk3SUiWrNmDfXs2ZO0Wm2r9ksoFNILL7xAN2/epPfff5+2bNlCCoWiQbnWtD1lyhQiapiSuri4mLKzs1uszd/fn5RKJR07dow2btzYaJns7GxKSEiguXPnGntI23KszEFr2/Xy8jL23BooKCiga9eu3bOW1h4Dc59bLT2nbZmamhpaunQpjRo1ioYOHWqxdtzd3amoqMhi9RPd7nUFQB999FGDtPBOTk4UFhZG+fn59NVXXzX63Rk7diwFBgbSDz/8QDU1NcbtOp2OtmzZQj169GgwZd1cFBUVmWSZNje24D+XdAT/a2pq6NFHH6W6ujqKiYlpstNtzpw5tHbtWotqaS2d1X8urzn16Qj+m+P6s2rVKjpz5gzt2LHDGPxYA+Y/0UsvvUSPPfZYo3X/8ssvRNTwkS5z0aT/jT1QumvXLgiFQsTExFj84WFr8ssvv8DOzg5+fn748ccfUVJSAq1Wi7y8PHz88cdwcXHBtm3bGv1sU0lnYmNjQUQ4efKkcdudGVErKipMMqJu2LDBpI4hQ4bA0dERo0aNwpEjR1BZWdlk1lypVIqxY8ciIyMDNTU1KCwsNC7IvHr1amO5wsJC+Pn5wdfXF7t374ZKpUJJSQk+++wzODg44Pvvv2/R/t1JRUUFpFIpeDweFixY0GiZ1rSdk5MDNzc3k4ytZ86cQUREBLp27driZEUGXnnlFdjZ2SE2NhbZ2dnQaDTIzc3FF198AS8vL4waNcokWUxrj1VrzoPmyre23eeeew5EhHXr1kGtViMnJwdz586Fj49Ps8mKWqKztcfA3OdWS89pW6WyshLTp0+Hq6srzpw5Y9G2DNcxS163f/jhBxAR7OzsoFarG7wfFxcHIoKXl1eTdRw+fBj29vaYN28ebty4geLiYixevBhCoRB79uyxiO6rV6+CiHDgwAGL1A/Yjv/1sVbinI7i//PPPw8iwqBBg1BTU9NkudmzZ+O9995r9D0ukhV1dv+5uObUp6P4f6/Xn02bNhkztTb1ujNxpzlg/t9m2bJl4PF4eOONN3DlyhXU1NTgypUrePnll41Ju6qqqsyuuzn/Gw1EAeC7776Dvb09Jk2ahLy8PLOL4orjx4/jscceMy5LIRKJ0L17d0RFRSE9Pb1B+cOHDzf4kqxYsQIAGmyfNm2a8XPFxcVYunQp+vTpAzs7O0ilUkRERODXX381lnnvvfeMn/Xx8cHRo0cxfvx4ODk5QSKRYOzYsUhLSzPRc+rUKSxevBgDBgyAg4MD3NzcMHz4cHz++efQ6/UmZUtKShATEwNfX1/Y2dnBw8MDkyZNQkpKSrP715w7OXcAACAASURBVET/hJHly5eDiJCRkdFkmZa0bSA7OxszZ86Ei4uLcXmRn3/+GRMnTjTqefLJJ5vVVJ+jR49CLpejR48esLOzg7OzM4YPH44PP/wQGo2mQfmWHKvWngfbt2+/6wW2Je0aUKlUeOqpp+Dl5QWJRIJRo0bhr7/+QmhoqLHu2NjYNp+vrT0G5jy3WnNO2xqnT5/GwIED4eHhYZEf1zs5efIkiAiZmZkWa0OlUkEoFGLcuHHNali4cGGz9Zw4cQJTpkyBi4sLnJycMGHChAbXO3Py3Xffwc7ODrdu3bJYG7bi/86dO5u8CWwsa7w56Aj+X79+HTwe76430oZX/UCIC8/r09n9B6x/zalPR/AfuPfrz7Rp0zgJRJn/tykvL8cXX3yBiIgI9O7dGyKRCE5OTggNDcU777xjkSAUaN7/ZiOOY8eOwc/PD87Ozli7dm2zvU8MBoPBuE1paSliYmJgZ2eHESNGtHh9wXtFq9XCw8MDb731llXasyXmzp2LMWPGWLQN5n/TMP+5hfnPLcx/bmH+c0tz/jebGjc0NJQyMzPpxRdfpNdee438/Pzoo48+oqqqquY+xmAwGJ2S4uJieu2116h3796UmJhI69evp7S0NOOyOZZGIBDQY489Rhs3bmxRBufOQklJCe3YsYMWLVpk0XaY/43D/OcW5j+3MP+5hfnPLXf1v6XRbG5uLl544QVIJBK4urrixRdfRHZ2ttmiZQaDwbBVDh8+DLlcDnt7e7i7u+Ptt99GeXk5J1qysrLA4/GQlJTESfvtkTfffBNSqdTk+XBLwfxvCPOfW5j/3ML85xbmP7fczf8WB6IGCgsL8c4776B3797g8XgYMWIE1q9fj5s3b96zWAaDwbAVLl++jNWrV2PAgAEgIoSEhGDDhg1W+bG7G3PnzkX//v1RW1vLtRTOuXnzJlxcXKya+Ir5/w/Mf25h/nML859bmP/c0hL/Wx2IGtDpdNizZw/kcjmcnJwgFAoxduxYxMXF4fTp022tlsFgMNolWq0W6enpWLFiBYYOHQoej4euXbvi+eefx9GjR7mWZ0JOTg5EIhHef/99rqVwzqJFi+Dt7W3RJBV3wvz/B+Y/tzD/uYX5zy3Mf25pif9tDkTrc+vWLXz//fdYsGABunbtCiJCjx49sHjxYuzYsaNdjBAwGAxGaykqKsK3336L+fPno0uXLiAi9OnTB0uWLMGePXtQV1fHtcQmWb16NcRiMU6cOMG1FM7Ytm0bZ9OkmP/Mf65h/nML859bmP/c0lL/eQBgzodS9Xo9HT9+nHbv3k27d++mY8eOkZ2dHT3wwAMUFhZGYWFhNGLECOrSpYs5m2UwGIx7Jjc3l9LT0+nQoUOUnp5Op06dIoFAQKNHj6YpU6bQtGnTKCAggGuZLUKv19PEiRPpxo0blJ6e3umuuWfPnqWwsDB69NFHaf369VZvn/nP/OcS5j+3MP+5hfnPLa3x3+yB6J3cvHmT9u7dSwcPHqRDhw7RuXPniIgoICCARo4cSaNGjaIRI0ZQ//79LSmDwWAwTNDpdJSZmWkMPNPS0ujatWskFAopODiYwsLCaNy4cTRx4kRydnbmWm6buHHjBo0cOZK6detG+/fvJ0dHR64lWYXr169TWFgY9erVi1JSUsje3p4THcx/5j8XMP+5hfnPLcx/bmmt/xYPRO+ktLSUDh06ZBxx+Ouvv6i6uprc3d0pJCSEgoODaciQIRQcHEz9+/cnoVBoTXkMBqMDUlNTQ2fOnKGTJ09SRkYGZWRk0KlTp0itVpNUKqWRI0fSyJEjKSwsjB544IEO9YNx4cIFGjVqFA0ePJi2b99us0F1S/n7778pIiKCRCIRHThwgFxdXTnVw/xn/lsT5j+3MP+5hfnPLW3x3+qB6J3U1dXR8ePH6c8//6RTp07RqVOn6OzZs1RbW0v29vY0cOBAk+B08ODB5OLiwqVkBoPRjikqKqKMjAyToPP8+fOk1WrJ0dGRBg0aRMHBwRQSEkIjRoygoKAg4vObXVLZ5jl16hRNnTqVunXrRrt27SIvLy+uJVmEjIwMmjp1Krm7u9OePXvazX4y/7mF+c8tzH9uOXXqFE2aNIk8PDxo//797UaXuWnP/kdERFDXrl1p37597UaXuWmz/xZ+VrVN1NXVISsrC1u3boVSqURkZKQxCRIRwdXVFWFhYYiOjkZcXBy2bt2KrKwsaLVarqUzGAwrUFdXh0uXLiElJQXx8fGIjo5GeHg4vLy8GlwnFAoFEhMTO/014sqVK+jfvz+6d++OAwcOcC3H7Hz77bdwcHCAh4cH0tLSuJbTgM7gv7OzMyZMmACVSsW1nAYw/7mF+c8Nhw4dwkMPPQQejwdXV1fmv5U5c+YM5HI5+Hw+HBwcmP+N0C4D0aa4evUqfv75Z6xduxZPP/00xowZYxKg2tvbY/DgwYiKisKKFSvw9ddfIy0tDbm5udDr9VzLZzAYraCurg6XL1/Gb7/9hi+++AIvv/wyZsyYAX9/f9jZ2YGIwOPx0LNnTzz44INYsmQJ1q9fj5SUFBQWFnItv11SUlKCmTNnQiAQYNWqVR1inbPy8nIsXLgQPB4PUVFRCAkJAY/Hw6xZs9rdUmId3f+lS5dCo9FwLalJmP/cwvy3HqmpqYiMjAQRYejQoUhMTMTNmzeZ/1bi1KlT+Ne//gU+nw+RSAQiwrhx4zBjxgzm/x3YVCDaFKWlpThy5Ag2bdqEV155BbNmzUJgYKDx4BMRxGIx/P39MWnSJERHR+M///kPNm/ejMOHD+PGjRtc7wKD0enQarW4du0aDh48iMTERKxatQoLFy7EuHHj0Lt3bwiFQuP319HRESEhIZg3bx6USiU2b96MEydOsKWh2oBer8e6detgb2+PgQMH2nTv7ObNm+Ht7Q13d3ckJycbt6ekpOC+++4Dj8dDZGRku0qh3xn8b88w/7mF+W85dDodkpOT8cADD4CIEBYWZqKrqqoKy5Ytg0QigVgsZv5bgNTUVEyePBk8Hg92dnYQCASYM2cOjhw5AoCd/43RIQLRptDpdLh+/TpSU1Px9ddf44033sCiRYsavdGVSCQYMGAAHnzwQTz++ONYsWIF1q1bhx07duDPP/9EXl5ep57Wx2C0Bo1Gg7///hvp6enYtm0b4uPjsXz5cjz22GMYM2YM/Pz8jKOa9TuKIiIisHjxYrzzzjvYvHkzjhw5goKCAq53p0OSk5ODqVOngsfjYd68eTh79izXklrMgQMHMHbsWPD5fDz99NMoLi5utFxKSgpCQ0PbZUDaGfxvzzD/uYX5bz40Gg0SExMREBAAPp+PyMhI/PnnnyZl0tPTERAQAKlUioSEBFy8eJH5b0ZSUlIwdOhQEJFxGq5CocC1a9caLc/O/3/o0IHo3airq8OVK1fw22+/YePGjVi5ciWeeOIJTJkyBQMHDjQuYG94CQQCeHt7Y9iwYZgxYwaee+45vP3229i4cSN27tyJI0eO4MqVK2yUhtFhKS8vR3Z2NtLS0pCUlIQNGzZg1apViI6ORmRkJIKDg+Hp6WnyveHxeOjWrRuGDh2K6dOn45lnnsFbb72Fb775BmlpacjLy2NT5zkkKSkJAwcOBJ/Px/z583H8+HGuJTWKXq/H3r17MX78eBARxo8f3+Bmq6nPJScnY+jQocabtJMnT1pBccvo6P63d5j/3ML8bztqtRrx8fHo3r07RCIR5HI5zp07Z1KmqqoKsbGx4PP5mDJlSoPAiPnfdnQ6HXbs2AF/f3/jvU737t2xfv36FscBzP9OHoi2hOrqauTk5CA1NRWbN2/GBx98gBdffBHz5s3DqFGj4OvrC0dHR5MbbyKCg4MDevXqhWHDhiEyMhILFy7Eyy+/jPfffx9ff/01du/ejb/++gs5OTkoKSnhejcZnZC6ujrcvHkT2dnZOHLkCHbu3IlNmzYhLi4OMTExeOyxxzBgwAD0798fPj4+EIvFDc5zFxcXBAQEYPz48ZDL5YiNjcWHH36IH3/8EYcOHcL169c7xHMQHR2dTodt27Zh8ODBICLcd999+Pzzz9tF0ocbN25gzZo16Nu3L4gIEyZMaNN0JkNAGhISAj6fj6ioqAY3bVxhC/736tXrnvxvz9iC//d6/rdnVCoVRowYYbyhZ/43T1FREZRKJdzc3ODk5ASFQoHr1683KJeamop+/fpBJpMhISGhyfrY+d86amtrsX79epMcNSEhIdi+fTt0Ol2r6+vs/nO+fEtHoaqqioqKiqigoICKi4upqKiIbt68SYWFhca/CwsL6ebNm1RcXEw1NTUN6nB1dSVXV1dyc3Mz/r+pv11cXEgqlZKTk5Pxxeh8lJeXU2VlpfFVWlpKZWVlxlf9v+98r6KiokF9jo6O5OHhQd26dSMHBwdKS0uj2tpasrOzI19fXxo0aBCNGDGCJkyYQAMGDCCxWMzBXjMsSVpaGn322Wf0ww8/EBFReHg4zZo1y7j8gjW4dOkS7dy5k3766SdKT08nFxcXWrBgAUVHR1NQUNA91Q2Afv75Z1q5ciWdPn2aZs+eTW+99Rb179/fTOrvjfbqv0gkIi8vL0pLS+vQvzft1X9znf/tjZs3b9LUqVMpLy+P9uzZQ2q1mvnfBAUFBfTZZ5/Rf//7XxKJRLRkyRJSKBTk5uZmUk6tVtPy5ctpw4YN9PDDD9PHH3/cYu/Y+d80Go2G3nzzTfroo4+osrKS7Ozs6OGHH6ZXX32VgoODzdJGZ/SfBaIcoVarqbi4uNnAobFt5eXlTdYpk8lMAtM7/5ZKpeTs7ExOTk4kkUjI2dmZhEIhSaVSEggEJJPJSCAQGG86HB0dyd7eniQSiRWd6XhUVlZSXV0dVVRUkE6no7KyMtLr9VReXk51dXVUWVlJNTU1VF1dTZWVlVRUVEQ6nY7Ky8tJrVYbg8yKigqqqKgw/q1Wqxttj8fjtbhTo/42d3d3cnBwaFDf5cuXKS0tjY4fP07p6el08uRJ0uv15OXlRaNGjaKwsDAKDQ2lBx54gEQikaXtZFiJsrIySk5Opp9++on27dtHNTU11L9/fxo9ejSNHDmSBg0aRIGBgY2eM61BpVJRVlYWnT59mtLT0+nAgQOUl5dHMpmMgoKC6OTJk5SWlkYhISFm2rPb6PV62rVrF73++uuUmZnZ7gLS9uD/9OnTadasWRQREUG5ubk0duxY8vX1pb1795Kjo6OZ9rR90t7874i/w1evXqWIiAiqq6ujvXv3Ur9+/YzvMf//4dKlS/TRRx/Rhg0bSCaT0eLFiykmJoZcXFwalD106BA9/vjjVFFRQe+99x4tWLCgTW0y//8hLy+PYmJiaPv27VRXV0fu7u700ksv0ZIlSyzWKdeZ/GeBqI2h1+uptLSUKioqGoyGlZWVGf9/69YtUqlUDQIZw2c0Gg2Vl5eTXq9vUbuOjo4kEonIxcXFGLTyeDwiuj2Sa0AqlRKfzyciMga6REQODg7G0bOmgtv65ZtCLBbf9YunVqtJq9U2W0aj0VBVVVWD7Yagkej2KLdGo2lQ3hA8EhHpdDrjyKJWqyW1Wk21tbV069Ytqq6ubnTkuzHs7OzIycmJRCIR3bx5k1xcXMjPz4969uxp0pHg4uJCjo6O5OTk1OiouJubG8lksha12VbUajVlZGRQeno6paWl0eHDh6mkpIQcHR0pODiYQkNDadSoUTR27Fjq2rWrRbUwrMOtW7coPT2dDh48SAcOHKDjx49TdXU18fl86tOnD/Xq1Yt8fHzIy8vL2LFBdPvaYPiOGDpXbt68SQUFBZSbm0tXr16l69evE9Hta8eIESNo9OjRNGbMGBo2bBgR3e4Rvn79Oh09epTc3d3Nvm96vZ5+/PFHWrlyJV24cIFmz55Nq1evJn9/f7O31Va48t/Ozs5ER2ZmJo0fP56GDh1KO3fu7DQzItqL/x2JrKwsmjx5Mrm5udGePXvI29u7ybKd1f+//vqL4uLiKCkpifz9/enll1+mRx99tNEO35qaGlq1ahWtXbuWIiIi6IsvviAvLy+z6Ois/u/evZtWrFhBGRkZREQUFBRk9NeadHT/WSDKoLKyMuPJeucInSEYMwR3KpWK9Ho9qVQqIjINxAx1Gagf6NYPDm/dukW1tbUmGu6spylaEmS2JFglMg2gDUgkErK3t29QjyFQJCLjqPGd9fD5fJJKpSQUCsnZ2dn4eUMQ7uTkRHZ2dsZg3dXVtUFdtbW1tGPHDtqwYQPt37+f/P396YknnqDo6OhG9bYXWjpq2h5+XBj3jk6no8uXL1NmZiadPXuWcnNzKS8vj/Lz80mlUlFZWRkBIJVKZfxeGM71rl27Urdu3ah79+7Uo0cPGjhwIAUFBVGvXr0abauwsJDuv/9+CggIoN27d9+1s6qtGALS119/nS5evEizZ8+mt99+22SUpr3QVv8NnYQjR45ssf93curUKZowYQKNHTuWtm3bZrHj0Z65m/8lJSWk0+mosrLyns//jsiBAwdoxowZNHToUEpKSmp0ZK85rHn94YKUlBRas2YN/frrr3T//ffTK6+8QjNnzjR+f+8kMzOTFixYQJcvX6b33nuPoqOjLaqvI/t/9epVWr9+PX355ZekUqlIKBRSREQEffLJJ9SzZ0+u5RFRB/TfrE+cMhgMs3Hu3DkoFAo4OjrC2dkZ0dHRyMjI4FpWiygvL0dKSgqUSiUiIyPh6uoKIoKTkxPCwsKgUCiwdetWFBUVcS2VYQOcOHECDg4OePHFFy3elk6nw9atW+Hv7w87OzvI5XJcvHjR4u1ag02bNkEsFqOmpuae6jl06BCcnJwwZ84ctqxZI7z77rvo1asX1zLaJUlJSZBIJJg5cyaqq6u5ltNuMKwBOmzYsEbXAG0MrVaLuLg4iEQihIWFIScnx0pqOxbl5eVITExEWFiYMfmQIQnUrVu3uJbX4WGBKIPRzikvL0dCQgKCgoJARAgNDUViYqJNZaPVarXIyspCYmIioqOjERgYCB6PByKCr68v5HI54uPjcezYsTZlnWN0fLZt2wYej4fPP//cKu11xIA0OzsbRITDhw/fc12//vor7O3tsXDhQrb80h28+uqrCAkJ4VpGu2PTpk0QCoV49tln2XX+/2NYA3TAgAHG5aWOHDly189dvnwZo0ePhr29PeLi4pifrUSn0yE1NRXR0dGQSCQQCATg8Xjw9PTE+++/f8+ddYyWwwJRBsOGSE1NRVRUFIRCIbp164bY2NgmF0xu76hUKpNRU5lMBiKCs7MzwsLCEBsbi+Tk5HaxWDWjffDKK6/Azs4OBw8etFqbhoC0X79+xoDUVkce9Ho9PDw88P7775ulvr1790IsFuP55583S30dhejoaISHh3Mto10RFxcHIkJsbCzXUtoFja0Bevbs2bt+Tq/X45NPPoGjoyNCQ0Nx5swZK6jtOJw5cwZKpdK4HJVUKgWPx0NgYCASExPZDA8OYIEog2GD5OXlQalUomvXrhAIBIiMjERKSopNj0wYRk0TEhIgl8vZqCmjATqdDtOnT0e3bt0aXTfPktTW1iIxMRF9+/Y1BqSXLl2yqgZz8NBDD2H27Nlmq++nn36CUChETEyM2eq0debMmYO5c+dyLaNdoNfrERMTA4FAgE8//ZRrOZzT0jVAm/rsQw89BKFQiNdff92mZkVxSV5eHj744AMEBweDiNCtWzf069fPZAq0Ld872TosEGUwbBiNRoOtW7ciPDwcRIT+/fsjLi4OZWVlXEszCwUFBUhOToZSqUR4eDgkEgmICC4uLggPD4dSqURycjJKS0u5lsqwEhUVFRg4cCBCQkI4eX6nfkAqEokQHR1t9aD4XoiLi4Onp6dZ6/z666/B5/Px9ttvm7VeW2XChAl45plnuJbBORqNBvPmzYNYLMb333/PtRxOuXr1KhQKBRwcHODu7g6lUomSkpIWf/63336Dj48PevbsadUZIbaKSqXCxo0bMXHiRAgEAkilUoSHh2PAgAEtfgaXYR1YIMpgdBCOHz+O6Ohok+RGp0+f5lqWWamrq2swakpEEAgECAwMhFwuR0JCArKyslgPZwfm8uXLcHd3x5w5czg7zoaA1M/PzxiQ5ubmcqKlNRw8eBBEhMuXL5u13i+//BI8Hg/vvvuuWeu1RYKDg/Hvf/+baxmcolarERERAScnJ+zbt49rOZyRmZkJuVwOOzs79OrVC/Hx8a3qQKurq4NSqYRAIMCsWbNaFbx2NmpqapCcnAy5XA5HR0eIxWJMmzYN0dHR6N+/v/EZ3D///JNrqYx6sECUwehgqFQqJCQkGIO0sLAwbN26tcNO48nPzzcZNbW3tzc++1F/1LSjjBIzbrNv3z4IhULExcVxqsMQkPr6+tpEQFpVVQWRSIRvv/3W7HV/+OGH4PF4nX4KZo8ePbB27VquZXBGcXExhg8fjm7duuHEiRNcy+GE48ePQy6Xg8/no2/fvkhISGj1b/DVq1cRFhYGiUSC+Ph4Cym1bQxJhxQKBdzd3cHn8xEWFob33nsPa9asQa9evWBnZ4fHH38c586d41ouoxFYIMpgdFB0Oh1SUlKMyY28vLwQGxtrU9MI20JdXR2OHTuG+Ph4yOVy9OnTh42adlDi4+PB5/ORlJTEtRRoNBokJCQYk49ER0cjLy+Pa1mNMmzYMCxZssQidb/55pvg8Xj44osvLFK/LeDo6IiNGzdyLYMTrly5An9/f/Tp0wcXLlzgWo7VSU1NRWRkJIgIQ4cORWJiYptyGmzduhUymQxBQUEdbmaTOcjKyoJSqTT+vgcGBkKpVCIrKwtr166Fl5cX7O3t8eyzz+Lq1atcy2U0AwtEGYxOQG5uLpRKJTw8PCASiRAVFYWUlBSuZVmNvLw8JCcnIzY2FmFhYRCLxSAieHp6IjIyEkqlEikpKaiqquJaKqOV/N///R8cHBzw119/cS0FwD8BqY+PD8RicbsMSF988UUEBwdbrP5XX30VAoEAW7ZssVgb7RWNRgMiwo4dO7iWYnUyMzPh4+ODQYMGtbtz3pIY1gC977777vn5w4qKCkRHR4PH40GhULBlROpx9uxZKJVKBAQEgIjg5+eH119/HefOnUNFRQXi4+PRrVs3ODo6QqFQtOuZKYx/YIEog9GJqKmpMUluFBAQgPj4eKjVaq6lWZXa2lqTUVNDKnehUIjAwEBER0cjMTERWVlZXEtl3IXa2lpMnDgR3t7e7erGo35A6uDgAIVCgfz8fK5lAbi9JqtAIIBKpbJYG8uWLYOdnR127txpsTbaI/n5+SAipKamci3Fqvzxxx+QSqUYN24cysvLuZZjFQxrgPr7+xufPzx69Gib6zt79iwCAwPh4eHR6b43TXHx4kWsXr0agwYNAhHB29sbL7zwgnEtZEMWYplMBmdnZygUCty4cYNj1YzWwAJRBqOTcuzYMURHR8PBwQEuLi6Ijo7u1IFXXl4etm7dCoVCgbCwMIhEIhARvLy82KhpO6e0tBT+/v4IDQ3lJJNucxgCUm9vb2NPPdc3SsXFxeDz+Ra92dXr9YiOjoZYLMaePXss1k57IzMzE0TUqdZ3TEpKgkQiwcyZM1FdXc21HItjGH3z9vY2rgGanZ19T3Vu2rQJEokEY8aMaTcdVlxx7do1xMfHIywsDDweD25ubpDL5UhOTkZdXR2A2xn1Y2NjTbIQszwQtgkLRBmMTo5KpUJ8fDx8fX1NkhsZLvidlcrKSqSmpiI+Ph5RUVHo2rWrcdQ0NDQUCoUCiYmJZs8+ymgbly5d4jyTbnPU1NQgISEBXl5e7SIgHTRokMXX/tTpdPjXv/4FBweHTrPkxIEDB0BEnHc2WItNmzZBKBTi2Wef7fDrOxcWFkKpVMLV1dU4+navszCqq6uhUCiMU3E7alLBu5Gbm2sSfLq6uhqDz/qeXLlyBQqFAvb29vD09ERcXFy763xktA4WiDIYDACmyY0EAgG8vb2hVCpRWFjItbR2w52jpnZ2diajpnFxcUhNTe0UowLtkQMHDkAkEkGpVHItpUlu3bqF+Ph4k4C0oKDA6jpeeOEFiz4nakCr1SIqKgpSqbTdPMdrSbZv3w4igkaj4VqKxYmLiwMRITY2lmspFuXy5ctQKBSQSCTo2rUrlEqlWdauvnr1Kh544AG4uLhg27ZtZlBqW9y4cQMff/wxxowZAz6fD6lUigULFmDXrl0NAvKcnBxER0dDKBSid+/eiI+PZ7+zHQQWiDIYjAbk5OQgNjYW7u7unTK5UUtRq9Umo6bu7u4gItjZ2ZmMmrKsfdbDsJ6lJZYnMSeVlZXG5BpOTk5QKBRW7fRJSkoCn89HUVGRxdvSaDSYNm0a3N3dkZmZafH2uOTLL7+Es7Mz1zIsil6vR0xMDAQCQYdeqicjIwNyuRxCoRB9+vRBfHy82R7NSE5OhqurK0JCQpCTk2OWOm2BvLw8rFu3DmPHjgWfz4ejoyPmz5+P7du3NxpYGo6BQCCAn58fEhISOv1srY4GC0QZDEaTGJIbjRw5EkSEkJAQJCQkoLKykmtp7ZZLly4hMTERCoUCoaGh4PP5jY6asmyIluPFF1+Evb09Dh06xLWUu3JnQBobG2uVRevLysogEAjwww8/WLwt4Pb6pePGjUPXrl079Hp+7777Lnr16sW1DIuh0Wgwb948iMVifP/991zLsQiGJVh4PB4GDx6MxMREswU/dXV1UCqV4PP5kMvlnWJa6fXr1xEfH4/w8HAIhUI4ODggMjISiYmJTSZKTE9Pb3AMtFqtlZUzrAELRBkMRoswJDeSSCSQSqWIjo7uVAk52oph1DQuLg6RkZHo0qULiAgODg4ICwuDQqHA1q1b2RRoM6LT6fDQQw/B09MTf//9N9dyWoQhIPX09DQGpOaY/tccoaGhFltPtDEqKysxatQodO/evcM+W/3qq68iJCSEaxkWQa1WIyIiAk5OTti3bx/XcsyKenSK5AAAIABJREFUXq9HcnIyRowYYbIEizmfN8/Ly8Po0aPh4OCAxMREs9XbHrkz4VD94LO5juz667COHDnS7MeA0f5ggSiDwWgVZWVliI+PR+/evcHn8xEeHs6SG7WS5kZNo6KiEB8fj9TU1E7xnJmlqKiowODBgxEYGGjRZUrMjVqtRlxcHNzc3ODs7IzY2FiLZYNcvnw5AgMDLVJ3U6hUKoSGhsLPz69DrjUZHR2N8PBwrmWYneLiYgwfPhyenp44ceIE13LMRm1tLRITExEUFAQej4fIyEikp6ebvZ0jR47A29sb/v7+OH36tNnrbw/8/fffJsGnTCYzJhy62wyglJQUDB8+/J7XYWXYHiwQZTAYbcKQ3MgwfcbHxwdKpRI3b97kWprNUVFRgZSUFCiVSkRGRsLV1RVEBEdHR5NRU2s8z9eRuHr1Kjw9PTFlyhSbm9ZlCEgNGTotEZDu3r0bRGT1gLCoqAhBQUHw9/fvcNll58yZg7lz53Itw6xcuXIF/v7+6NOnDy5cuMC1HLNgyGLds2dP2NnZQS6XW2z5su+++w4SiQSTJ0/ucEuMnD59Gm+++SaGDh0KIkKXLl3wxBNPYPfu3XftSNXpdEhOTsZ9990HIkJ4eLhxfVBG54EFogwG4565ePEiYmNj0aVLF4jFYpbc6B7RarXIyspCYmIioqOjERgYCB6PByKCr68v5HI54uPjcezYsQ6/ZMK9cujQIdjb22Pp0qVcS2kTlgxI1Wo17Ozs8N1335mlvtZQUFCA/v37Y/DgwVZ5JtZaTJgwAc888wzXMsxGZmYmfHx8MGjQoA4xgl1eXm58JlssFkMul+PixYsWaUur1SI2NhZEBIVCYXOdYY2h0+mQmpqKZcuWwc/PD0QEb29vPPPMM9i7d2+LZkYZRqEDAgLA5/MRGRnZKTJqMxqHBaIMBsNsVFdXIzExEcHBwSAihIaGIiEhoVMkZLA05eXlJqOmMpkMRAQnJyeEhYUhNjYWycnJKC4u5lpqu2Pz5s3g8Xj46KOPuJbSZioqKhAXFweZTAY3NzcolUqzTDkeOXIknn76aTMobD3Xrl1D7969MWzYMFRUVHCiwdwEBwfj3//+N9cyzMIff/wBqVSKcePGoby8nGs590RBQQGUSiWkUilcXFygUCiQn59vsfbKy8sRGRkJe3t7m38etKamBikpKVAoFPDy8jJ2iCoUCqSmpra4M1Sj0SAxMRH9+vUzjkJ35MRljJbBAlEGg2ER7kxupFAocOnSJa5ldRjYqGnr+M9//gM+n4+ffvqJayn3RP2AtEuXLvcckK5YsQK+vr5mVNg6Ll68CG9vb4wcObJDZOPu0aMH1q5dy7WMeyYpKQkSiQQzZ8606fUac3JyoFAoYG9vD09PT7N14DRHdnY2AgIC4OPjg6NHj1q0LUtRWVmJ5ORkyOVyuLi4gIgQGBgIpVKJY8eOtbqu+Ph4+Pj4QCQSQS6Xd5gp3ox7hwWiDAbDohQWFiIuLg69evUySW7UEaYptTcKCgqQnJwMpVKJ8PBwODg4gIjg7OxsMmrakaZCtoZnn30WEomkQzyHVF5e3iAgbcuo1YEDB0BEyM7OtoDKlnH+/Hl4enoiPDzcpoMeAHB0dMTGjRu5lnFPbNq0CUKhEM8++6zNdmKdPHnSZP3J+Ph4q5xbu3fvhkwmw8iRI23u+eebN28iMTERkZGREIlEEAgECAsLQ3x8PK5fv97q+ioqKozToB0dHaFQKJCbm2sB5QxbhgWiDAbDKtyZ3MjPzw9xcXEsAY8FqaurQ1ZWFhISEiCXyxEYGAgigkAgQGBgIORyORISEpCVldUpUuRrtVpMnz4dHh4eFnsuzNqUlJQYpxy6u7u3OiDVarVwc3PDf//7XwuqvDunTp2Cm5sbHnroIdTW1nKqpa3U1NSAiLBjxw6upbSZuLg4EBFiY2O5ltIm6i//ERwcbLX1J/V6PeLi4sDn8/HYY4/ZTIfKpUuXjJlu+Xw+JBIJIiMjkZCQ0OYlxYqKiqBUKo3PtSsUCpsLyhnWgwWiDAbD6ly4cAGxsbFwc3MzJjeyRMp8RkPy8/NNRk0lEgmICC4uLggPD4dSqURycnKHy+5ooKKiAsHBwRgwYIDF1+m0JsXFxVAqlXBxcYG7uzvi4uJa/Gz2I488ggcffNDCCu/O4cOH4ezsjNmzZ9vkclD5+fkgIqSmpnItpdXo9XrExMRAIBDg008/5VpOqzBkXx02bBgny39oNBrI5XIIhUKsW7fOau22laysLCiVSoSGhhoz3crlcmzduhVqtbrN9RYUFCA2NhaOjo7GTrGOdI1lWAYWiDIYDM4wJDcaMmSISXKjqqoqrqV1Gurq6nDs2DHEx8dDLpejT58+HX7UNC8vDz179sSYMWPuur6drWEYjXBxcYGHhwfi4uLu+n366quvIBKJ2kXCoLS0NDg6OuLxxx+3uWmhmZmZICKcOXOGaymtQqPRYN68efh/7N15XJRV+z/waxYYBhj2RUQEURFQREFNJVFxSC3U0lBLMTMdS2hyySZ9MswWp7LEXBKffJS0VDBFzBZxQcGlxAVFU3FFWVV2kGWYz++PvsxPBBR0Zu4Bz/v14h9muK9rDozOdZ9zriMSibBt2zau02m2uuY3D3ZfPX78uF5zKCoqQlBQECwsLLB37169xm4ulUqF5ORkyOVydOjQAUQEV1dXyOVyJCYmPvUKhJs3b0Iul0MsFsPR0RFKpbJN7Pdm9IMVogzDGITU1FSEhYXByMgIVlZWkMvluH79OtdpPZOysrKQkJAAhUKBgIAAmJiYgIjg6OiIkJAQREZGIjExsVV3Q05PT4eVlRUmTpzYZgrsB7WkIM3Pzwefz0d8fLyes2zc3r17IRKJEB4eznUqLVK33zY3N5frVJqttLQUw4cPh7m5ucEWUg8rLS1FVFQUOnTooGl+c+HCBb3nkZWVhV69esHJyQmnTp3Se/xHKSoqwrZt2zBp0iRNh/XevXtj8eLFOHPmjFZiXLt2DTNmzICxsTFcXV3x/ffft7kbe4zusUKUYRiDkpOTA6VSiY4dO2qaGyUkJLTJYqG1qK6urjdr6urqCiKCUCiEt7c3ZDIZYmJidHYgvK4cOHAAxsbG+Pjjj7lORWfu3LkDhUIBU1NTODg4NFmQ9u3bFzNnzuQgw8bt3LkTQqGwVZ3/unPnThARqqqquE6lWe7evYv+/fvD0dHR4AqpxtTdXLGxsYG5uTnkcvkTNdHRhnPnzsHFxQXdu3fHzZs3OcnhYZcuXcI333yDoKAgGBkZQSgUYsiQIVi+fLlWb+peu3YNMpkMRkZGcHNzQ1RUFCtAmSfGClGGYQySSqVCQkICpFIpeDweunTpAqVSyc7JNBAPz5qKRCIQEdq1a1dv1tTQl1mvX78eRNTq9sW1VH5+PhQKBcRisaYgfbChSmRkJDp06GBQN3zi4uIgEAjwySefcJ1Ks/zwww+QSCRcp9Es169fh4eHBzp16mTwR2ncuHEDcrkcpqammr2HXHb+3rdvHywtLREUFMTpXvq6JbcKhULTiM7GxgahoaGIjo7W+sz81atXIZPJIBQK0alTJ0RHR7fKvdyMYWGFKMMwBu/SpUtQKBSwtraGiYkJwsLCcPr0aa7TYh5QXl6O5ORkREVFITQ0FA4ODo3Oml67do3rVBtYuHAhjIyMWs3SxKfxYEHq4uKiOdbir7/+AhEhLS2N6xTr2bBhA3g8HpRKJdepPNZXX30FV1dXrtN4rHPnzsHZ2Rk+Pj7IysriOp0mnTt3TrNdw9XVFVFRUZxvB9i4cSOMjIwQGhrKSWfcu3fvIjY2FmFhYZolt+7u7lrb79mY9PT0ekfhsAKU0SZWiDIM02qUlJQgOjoaPj4+rLlRK5CVlYXY2FjI5XIEBATA2NgYRAQnJyeEhIRAqVQiOTmZ86MO1Gq15uB2be2fMnR5eXn1CtLly5fDwcEBS5cu5Tq1Br777jvweDysXr2a61QeacGCBejduzfXaTxSUlISLC0tMWTIEBQVFXGdTqNSU1MRGhoKHo+HHj16ICYmxiCO9FEqleDxeJDL5XptpFV3xIpUKoVQKIRQKERAQACUSiX++ecfncU9e/aspgDt3r273o7CYZ4trBBlGKZVerC5kYODAxQKBW7cuMF1WswjlJWV1Zs1tbe3BxHByMgI/v7+kMvliImJ4eT3WFlZicDAQHTs2NGgZ4m0ra4gNTExgZmZGTp37myQ+72++eYb8Hg8/Pe//+U6lSbJZDJIpVKu02hSfHw8xGIxXn75Zc5v/jTmwTNA/fz8EBMTYxCdk2tqajBt2jQIhUJER0frPN79+/eRmJgIuVyOjh07gohgZ2eH0NBQxMTE6PwGwpkzZzQ3Anx8fFgByugUK0QZhmnV6pobubi4sOZGrdDVq1cRExMDuVwOf39/8Pn8RmdN9VEc3bt3D15eXvDx8Wmz56g2JTMzEyNGjAARoUOHDgbZgGThwoUQCATYsmUL16k06tVXX8X48eO5TqNRGzZsgFAoxKxZswyiuKtTdwZonz59ODkD9HGqqqowbtw4mJqa4tdff9VZnLy8PMTExCA0NBQSiQREBG9vbygUCiQnJ+vld3b69GlNAerr64vY2Fj2/yijc6wQZRimTaiqqkJsbKymuVHXrl2hVCo5bWrBtFxpaSmSk5OhVCoREhICW1tbEBFMTU0REBAAuVyO2NhY5OXl6ST+7du3NWeMGuKskS4VFBRAKBRixIgRMDExgaurq8HtB5s/fz6MjIw4L1ZycnKQn59fb6YoKCgIb7/9NodZNa5uSalCoeA6FY3KykrExMTAw8NDcwbo33//zXVa9VRWVmLMmDEwNzfH/v37tX799PR0KJVKBAQEgMfjQSwWQyqVIioqSq/dgFNSUjQz0QMGDGA3chm9YoUowzBtzj///AO5XA5zc3NNc6NnZe9fW/SoWdPQ0FBERUUhOTlZa8dmpKenw8bGBqNHjzaoIkwfgoKCMH78eM0h9SKRCG5ubgZTkKrVarz99tswNjbGb7/9xlkeCoUCRAQigpmZGZydnWFpaYkuXbrgtddeQ3h4OBYtWoTly5frdB/fo6jVasydOxcCgcBgukKXlJQgKioK7du315wBeunSJa7TaqCsrAzBwcGwsrLCsWPHtHLNwsJCxMXFYdq0aXB0dAQRoWPHjnjnnXfw22+/6f3G14NLoQ1tJpp5drBClGGYNqu4uBjR0dHo0aOHprmRoTS+YJ5cSUlJvVlTGxsbTUHw4Kxpfn7+E8c4duwYTE1NIZPJtJi54VuzZg1MTU1RVlYGAAZZkNbW1mLy5MkwNTVFUlISJzns27dPU4g+/MXj8TTnOAqFQk7OuqyqqsLEiRMhEomwbds2vcd/WF5eHiIjI2FlZQWJRAK5XI7bt29znVajioqKMHDgQDg4ODzVDcza2lqcOHECn332GZ5//nkIhUIIBAIMHDgQX3zxBWcdqpOTkxEUFMQKUMYgsEKUYZhnQnJyMkJDQyEUCuHo6AiFQmEwB5EzT0elUiE9PR0xMTGQyWTw9vYGj8drMGuampraor1WCQkJEAqFreYcS23Izc2FQCDA9u3b633/xo0bBnWGoEqlwoQJE2BhYcHJks7KykqYmJg0WYzWNeHiYs9oaWkphg8fDnNzc86PJLp27Rrkcrnm/NrIyEgUFBRwmtOjFBQU4LnnnkO7du1w7ty5Fv/8nTt3EBsbC5lMBmdnZxARHBwcNI2GuNwqkpiYiP79+2sKUF0sN2aYlmKFKMMwz5Ts7GwolUp06NBBszcpMTGR7YlpY4qLi5GYmIjIyEiEhIRoztwzNzdHQEAAFAoFEhIScPfu3UdeZ926dSAirFy5Uk+Zc2/IkCGYMGFCo49dv369QUHKVUfN6upqze/25MmTjT6npqYGJSUlOokvlUo1y8Sb+jp69KhOYjfl7t276N+/PxwdHXHq1Cm9xn5QWloawsLCNH8nUVFRBn/MVm5uLnr27AlXV1dkZGQ062dUKhVSU1OhVCo1x6sIBAL4+/sjMjISqampnP/fkpiYiH79+oGIIJVKtbbUmGG0gRWiDMM8kx5sbkRE8PDwgFKpNOi79cyTe9Ssqbu7O8LCwpqcNY2MjGx0lrCtWr16NUxNTVFeXt7kcx4sSL28vDg74qGiogJDhw6Fvb09zp8/X++xqqoqvPzyy4iMjNRJ7GXLlkEoFDZagPL5fPj4+Gg95qNm9K9fvw4PDw906tQJly9f1nrs5qjbd8jj8dCzZ0/ExMRwvpS7OTIzM+Hh4YFu3bohMzPzkc/Ny8tDbGwswsLCNNsC2rVrh7CwMMTGxhrE+axqtbpeN2KpVIq//vqL67QYpgFWiDIM88y7cOEC5HI5zMzMIJFIIJPJcPbsWa7TYnQsNzcXCQkJiIyMhFQqhampKYgIEomk3qzpvXv38O6778LExASHDx/mOm2dq1ueu2PHjsc+99q1a5wXpOXl5Rg0aBCcnZ1x9epVAP8WqMOHDwcRwcrKSiezcWfPnm1yJpTP52PTpk1ajXfr1i2MGDGi0aY26enpcHZ2ho+Pj97Pwa0regYMGFBv3yHXM4HNde3aNbi5ucHHxwe5ubkNHq+b9YyMjIS/vz94PB5MTEwglUqhVCqRmprKQdaNq62txZYtW9C9e3fw+XyMGzeONepjDBorRBmGYf5PXXMjb29v1tzoGVRTU4P09HRER0cjLCxM83cgEAjg7e2Njh07QiwW45dffmk1H7KfVGBgIF577bVmP//q1auQyWSasdJ3QVpUVIQ+ffqgY8eOuHDhAoYMGaKZrRQIBFizZo3WY6rVatjZ2TVaiFpbW2v9HNbJkyeDiDB69Oh6Y5uUlARLS0sMGTJEr7Nx1dXViImJ0awuCAkJwZEjR/QWXxtu3rwJNzc3+Pn51Vumn5OToznX09LSUrNyQiaTITY2VmfLvZ9U3c2AXr16abacNLVcnWEMCStEGYZhHqJWq5GYmKhpbtSuXTsoFIrHLtli2p6cnBzNrGlQUJBmT6C5uTmkUikiIyORkJCAwsJCrlPVqpUrV8Lc3LzFM4kXLlxAWFgYBAIBunfvjpiYmBY1iHoad+7cQc+ePeHh4QEjI6N6XWxdXFx0UhhPnjy5Xqy6JkWLFy/Wapy0tDTNUnKBQIC33noLABAfHw+xWIyXX35Zb8d/lJWVISoqCh07doSRkRHCwsKQnp6ul9jalJubC09PT/j4+OD27dv4448/MHv2bHh6emq6cI8aNQqrV6/WzLQbGrVajV27dsHX1xd8Ph8TJkzAhQsXuE6LYZqNFaIMwzCPkJWVhcjISNjb28PIyAihoaGsudEz7O7du+jWrZumE6a7u3u9WdOwsDBER0cjPT29Vf+N5OTkQCAQYOfOnU/08+fPn9cUpD169EBsbKzOx6OgoAC9evVqUBjWFaO//PKL1mP++OOPDRoWCYVC5OTkaDVOcHBwg+J68uTJEAgEmDVrll6K/eLiYkRFRaFdu3YQiUQICwtrdlMfQ5Ofn48uXbrA3t4egYGBmg7IdbOeCQkJWp/R1rbExET07dtXsweUzYAyrRErRBmGYZrh4eZG3bp1Q1RUFEpLS7lOjdGzrKwsuLq6ol+/figrK0NWVhYSEhKgUCgQEBCg+VBraWmpmTVNTEx8ZPMfQzRo0CBMmjTpqa6hr4I0NzcXXl5ejRahdTcK/Pz8dBK3bqaybjZ08uTJWo2RlJTU5F7UsLAwrcZqTE5ODiIjI2FpaQkLCwvI5XJkZ2frPK62ZWZmYv369XjllVc0y7bt7Ozw+uuvY8OGDQZ7runDkpOTMWTIEE0BeuLECa5TYpgnxgpRhmGYFjp58iRkMhlMTU1hYWEBmUz2RGfOMa3XhQsXYGtri1GjRjXoClpdXY3U1FRERUUhLCwMbm5umpkyb29vyGQyxMTEGPxyxhUrVkAikWhlyWd6ejrCwsI03WS1WZBmZmbCzc2tySL0wS9d7GGsW8pZ96XNc03VajX8/f2b7M7L4/EQExOjtXgPunLlCuRyOUxMTODo6IjIyEiD6AjbXGVlZUhMTIRCodA0GRKLxbC0tIS1tTV+/fVXvS0b14ajR49qboQGBAQgKSmJ65QY5qmxQpRhGOYJFRUVITo6Gl5eXpoPB7Gxsa3iuALm6R09ehRmZmaYNGnSYz/QPjxrKhKJNMc+hISEaGZNDemsxezsbPD5fOzatUtr1zx37hxCQ0M1x3u0pCBtqhvs7NmzwePxmizWHpytHD16tNZeS5158+bB2NgYAoEA/fr10+q1t27dWm/GtanZ3j/++ENrMU+dOqWZxe7cuTOioqL0tv/0aTx8pqexsXG95bZbtmxBcHAw7Ozs8M8//3CdbrNlZmYiLCwMRISBAwdi//79XKfEMFrDClGGYZinVFtbq2luJBAI4OTkBIVCgVu3bnGdGqNj+/btg4mJCd58880WzfCVl5cjOTkZUVFRCA0NhYODQ6Ozplw3SQkICHjq5bmNOXv2rKYg9fX1fWxBev36dVhZWTU5o3n27FmMHTsWPB7vkTOjPB5P681c/vzzT821t27dqrXrVldXw9XVtcEe1MZek6mpaYMlmmq1ukVnR9adAUpE6NWrF2dnw7bE1atXER0djdDQUFhbW4OI4OjoiNDQUERHR2uW26rVakydOhWmpqY4evQox1k3T0lJCRQKBUQiEby9vfH7779znRLDaB0PAIhhGIbRimvXrtG6deto/fr1VFJSQmPGjCGZTEZSqZTr1Bgd2bVrF7366qsUERFBy5cvf+LrZGdn05EjRyglJYVOnjxJJ06coOrqanJyciJ/f396/vnnKSAggPr06UMmJiZafAVNW7NmDc2fP59yc3NJIpFo/frnzp2jTz/9lLZv3049e/ak//znPxQaGtrgedOnT6f169eTmZkZJSUlUZ8+fRq9XlpaGn3yySe0c+dOEgqFpFKp6j1uZGRE06ZNo7Vr17Y415qaGsrPz6e8vDwqKioitVpNxcXFVFlZSVOnTiWJRELx8fFkZ2dHzs7OZGlp2eIYD1q1ahXNnj2bamtrm3yOkZERqVQqGjx4MC1atIiCgoKIiEitVpNMJqOEhATKzMxs8u9FrVbTnj176LPPPqO///6bAgICSKFQ0KhRo54qd10pKyuj48eP0+7du2n37t10/fp1MjU1pYEDB5JUKiWpVEp+fn7E4/Hq/dzcuXNp9erV9Ouvv1JwcDBH2TcPANq4cSMtXLiQqqurafHixfTOO++QUCjkOrVG3blzh/Ly8qiwsJCqqqro/v37VFlZSUZGRmRubk5CoZCsra3JwcGBHB0dic/nc51ym9Lqx5/jQphhGKZNqqysRGxsLAICAkBE8PLyQlRUFMrKyrhOjdGBTZs2gc/n4/PPP9faNcvKyurNmtrb22uWmPr7+0MulyMmJgbXr1/XWsyH3bt3DyKRCBs3btRZDODf40nqZkj79++PhIQEzWM3btyodyaoubn5YzuEHjt2DCNGjNDMMtNDS3Qf1dW2tLQU+/fvx/LlyzF9+nT0798f7dq1e+wS2Ye/TE1N4eHhgTFjxmDhwoXYsmULLl261KzxKC0thY2NTZMzoAKBAGKxGDNmzGiwP12lUuGNN94An88Hn8/H2rVrG1y/qqoKMTEx8PT01Jw7efz48Wblpk8VFRXYt28fFi5ciH79+kEgEEAgEOC5557DRx99hEOHDj32nOfly5eDz+djy5Ytesr6yRUXF2PcuHEQCAR49913ce/ePa5TAvDvqp/Tp0/jf//7H+bNm4cXXngBbm5umi0Gzf0SCoVwdnbG0KFDER4ejrVr1+Lo0aOoqqri+iUatLY8/mxGlGEYRsdOnjxJ69ato82bN5NQKKSJEyfSe++9R97e3lynxmjRmjVrKDw8nL755huaO3euTmI8OGt65MgROn36NKnV6gazpn379iWRSKSVmOPGjaPi4mLat2+fVq73KGlpafT555/T9u3bqX///rRgwQLauXMnbd68mWpqaoiISCgUkkQioZSUlMe+h44cOUIfffQRJSUlkUAgoNraWhIKhfSf//yHFi9eTEREtbW1lJycTL///jsdPnyYUlNTSaVSkb29Pfn4+FD37t2pS5cu5OzsTE5OTuTo6EjW1tbE4/HIwsKCBAIBLVu2jF555RUSi8VUUFBAt2/fptzcXLp58yZduHCBzp8/T5cvX6aamhpycnKiwYMH09ChQ2n06NHUrl27Bnl//PHHtHTp0nozunUzvB07dqRZs2aRTCYja2vrej9XW1tLU6dOpZ9//pnUajXxeDzq0KEDXbt2jYRCIZWVldH69etp2bJllJ+fTxMmTKAFCxaQl5fXU/7mtEOlUtGJEydo//79dODAATp27BhVVlZS165dadiwYRQcHExDhw5t8LqbsmfPHhozZgwplUp6//33dZz90/nnn3/o1VdfpTt37tBPP/3E+cztlStXKCEhgZKSkig5OZmKiopILBaTt7c39ejRgzw9PcnZ2VnzvrCxsSFjY2MyMTEhsVhM1dXVVF5eTiqVigoLCykvL49u375NOTk5dOnSJUpPT6fz589TcXExmZqa0oABA2jw4MH00ksvkZ+fH6ev3RA8M+PPWQnMMAzzjCksLERUVBQ6deoEHo8HqVTKmhu1MV988QV4PB5++OEHvcQrLS1FcnIylEolQkJCYGtr22DWNDY2Fnl5eU8cY+fOneDz+cjMzNRi5o9WN6MpFAohEAgavbNvY2PT7KYz+/btQ79+/TQzipaWltizZw/eeust2NnZgYjg6emJd955B1u2bGmyMVJTmvMerqqqwpEjR7B06VKMHDkSZmZm4PP5CAgIwLfffos7d+4AAPLy8iAWi+vN4PJ4PAwdOhQJCQlN7qVVqVR4/fXXG4wXj8fDunXrEBkZCRsbG5ibm0MulxvMHva60yKBAAAgAElEQVQH93laWVk12Of5pDP+qampMDMzw/Tp07WbsA5s3LgRpqamCAwM5PRonKtXr+Ljjz+Gj48PiAi2trZ49dVX8d133yEtLU0nXYYvX76M9evXY8qUKXBxcQERwc3NDXPnzsXp06e1Hs+QPYvjzwpRhmEYPXu4uVH79u0RGRmJ/Px8rlNjtODDDz+EQCDAtm3bOIl/9epVxMTEQC6Xw9/fX9PsxsnJCaGhoYiKikJycnKzl2NVV1fDzs4OS5cu1XHmDY0dO7bJ5kNCoRCOjo64du1as68XFxcHZ2dnzTX69u0LpVLZ7CWz2lRRUYH4+Hi88cYbsLKygkgkwuuvv44xY8ZoCkhzc3PMmzfvsU2rqqur8fLLLzdatPP5fJiZmcHe3h6ffvopCgoK9PQKG1dXeIaFhWl+FxKJBFKpFEqlEqmpqU99tM/169fRrl07jBgxwqBv9KlUKigUCvB4PMjl8scuM9aF2tpaxMfHY8SIEeDz+Wjfvj0iIiKwf/9+TppVpaamYuHChejWrRuICM899xw2bNjQKjo3P4lnffxZIcowDMOhK1euQKFQwNbWFsbGxggNDUViYiLXaTFPac6cOTA2NsaePXu4TgUlJSX1Zk3r9h6amZkhICBAM2v6qBshERER8PLy0mPW/743Htcx1sjICE5OTrhx48Yjr1VcXAylUgkrKytYWFjgpZdewgsvvGAwXWHv37+P2NhYDBw4EEQEsViM8PBwlJeXP/Znq6qqMGrUqEaL0Ae/tHkMT0vk5OQgNjYWMplMc6aumZlZvcJTmzM9RUVF8Pb2hq+vL0pKSrR2XW0rKSlBSEgITExMsGnTJr3Hr62tRUJCAnr37g0+n69ZocNFMdyU1NRUyGQyiMViODg4QKlUGtQRV0+Djf+/WCHKMAxjAOqaGw0YMABEBD8/P0RHRzfrgyhjeNRqNd566y2IxWIcOnSI63TqUalUSE9PR0xMDGQyGby9vTWNeB6cNX2wQDhx4gSIqMnjU3Rh8uTJjzyK5cFi1M3NrdEljdXV1fjyyy8hkUhga2uLL774QlOcqNVqgylE62zatAnfffcdpFIpiAhDhw5Fenp6k8+vrKzESy+99NgzVIVCIQYNGqSX13Dnzh3Nmbn+/v6aM179/f2hUCiQmJios+YoKpUKI0eORPv27Q1m6XFjMjIy4OXlhfbt2+Pvv//We/z9+/fD29sbAoEAkyZNMvhzVbOysvDee+9BLBbD2dkZP/30E9cpPRU2/v8fK0QZhmEMzIN3IS0tLSGTybR+9iGjeyqVCqGhobCwsGhwxqOhKS4uRmJiIiIjIxESEqLZq2dubo6AgAAoFAq4u7vjtdde00n8hwvCy5cvP3Y29OFCy9PTU7PPEgAOHz6M7t27QywWY/HixQY9O9aYlJQU9O3bF0ZGRnj//fcbdNwuLy9HUFDQY4vQB7/qztDMyMhARkaGVvIsKCjArl27MHv2bPTq1Qt8Pl/T2XbBggVITEzU2yzWBx98ABMTE4PsAFzn2LFjsLOzQ+/evfW67xoAcnNz8frrr4OIMHr0aFy8eFGv8Z9WdnY2ZDIZ+Hw+hg4davAF3MPY+DfEClGGYRgDlZeXB6VSCTc3t3pLdwxtFodpWlVVFUaOHAk7OzucP3+e63Sa7VGzpm5ubggLC2swa/o0Dh06BKlUisuXLwMA1q9fj549e0IikdQrpIyNjZucJTUyMoKPjw/u3LmDhQsXgs/n48UXX3zs/kpDVltbi++//x7W1tbw8vJCWloagH+P9hk0aFCLilAej4fRo0fjq6++gkgkwpdffvlEORUWFiIhIQFz5szRLCvk8/nw9fXFe++9h4SEBBQVFWlzGJolLi5Or43CnkR8fDzEYjHGjBmj99Uuv//+OxwcHODm5sbZMm1t+euvv+Dn5wdTU1OsW7eO63SahY1/41ghyjAMY+DqmhuFhISAx+PB3d0dSqWy3uwPY7jKy8sxaNAgODs7t6ixjqHJzMyESCTCSy+9BKlUClNTU02jmbpZ04SEhCc6+1CpVGqKySVLlqCyslLzWFFREU6fPo0dO3Zg2bJlCA8Px/Dhw+Hu7g5jY+MGBZe5uTnEYjGio6O1+fI5devWLQQGBsLExAQrV67ULOF/uBAXiUQNilOBQAB7e3t069YNlpaW4PP54PF4GDFiRLNil5aWIjExEQqFAgEBAZobAe7u7pDJZIiNjcXdu3d1PAKPdurUKZiammLu3Lmc5vEoq1evhkAgwDvvvKPXm4m1tbWahkiTJk1CcXGx3mLrUk1NjeaG04QJEwx2Gwsb/0djhSjDMEwrcvnyZSgUCtjY2EAkEiE0NBQpKSlcp8U8RmFhIXr37o3OnTu3+GgQQzJlyhT4+voC+PeDSHp6uqYDqre3t6b4cXd3R1hYGKKjo5Genv7YWdOQkBDNUlyhUAhXV1fs3bv3sfmo1WpkZWUhOTkZX375JaysrGBpaYm33377qTuvGhqVSoX58+fXm+G0sbGBt7c3Ro4ciRkzZiAyMhLR0dHYvXs3Tp48iezsbFRWVkKpVEIoFNYrUs3MzBotiJpTeBrSTbDc3Fy4uLggODjYYDvk1h3r9Pnnn+s1bmVlJUJDQyESibBhwwa9xtaXffv2wdbWFgMGDOD8hsjD2Pg/HitEGYZhWqH79+8jJiYGvr6+ICL4+/uz5kYGLi8vD56envD09OT0rMCnkZycDCLCyZMnG308JycHCQkJiIyMhFQq1ZyHaWFhAalUisjISCQkJDQ4QqSuk++Ds3hEhHHjxjXrWKNTp04Z7IdRbbp06RIWLVoEoVCI6dOnP7bYPnr0KDw8PJrsqHv69OlWV3g+qKamBkOGDEHnzp2faCZeH5RKJXg8HlasWKHXuBUVFRg6dCisrKyQlJSk19j6dvHiRXTq1Amenp7IycnhOh0AbPybixWiDMMwrVxdcyMTExNYWVlBLpe36iWgbVlubi68vb3h4eHRamdGvby8IJPJmvXcqqoqHDt2DMuXL8f48ePRoUMHTaHZq1cvzJo1C19//fUjO+JKJBJERUU1OauakZEBR0dHBAcHt5mjHR4nISEBRkZGUCgUjT5eXl4OhUIBPp/f5D5SIyMjTedOHo+HHj16ICIiAtu3bzfYwvNhdc2JmroxwrX58+dDIBBg/fr1eo1bU1ODUaNGwdbWFmfPntVrbK5kZ2ejW7du6NWrFyd7lB/Exr/5488KUYZhmDYiNzcXSqUSrq6urLmRAcvNzYWXlxe6devWKmdGV65cCVNT0yeegcrKytIc7xEQEAATE5PHdsjl8Xjo379/g6NMCgoK4O7ujr59+6K0tFQbL6/V2LRpE3g8HlatWlXv+0lJSXBzc3tsIyOBQAA/Pz/ExcU1a9bZ0CQkJIDH4+F///sf16k0asmSJRAIBJwcNTJz5kyYmZnh2LFjeo/NpRs3bsDZ2RlBQUGc/r/Hxr/5488KUYZhmDam7qBsqVQKHo+Hzp07Q6lUtukli61NTk5Oqy1GS0pKYGFhgWXLlmnlejNnzmy06VBjM3gCgQByuVxzlMm4cePQoUOHVllIacOnn34KkUiE06dPo7CwENOnTwePx2v20TdWVlatci9tRkYGrKysmj0zr29r1qwBj8fjpGHWtm3bwOPxEB8fr/fYhuDMmTMwMTHB4sWLOYnPxr9l488KUYZhmDbs0qVLUCgUsLa21jQ3qjtLkOHWrVu30KVLF4Pa19RcERERcHNz08qsQ/fu3Zt9BEndTF7Hjh0xb9488Pl8HDx48OlfUCtVW1uLoUOHwtnZGba2ti0ax7qv1nZG8f379+Hn5wdfX1+DXIr9888/g8/nQ6lU6j32rVu3YGVlhfDwcL3HNiSrVq2CUCjU+/91bPz/1ZLxZ4UowzDMM6CkpATR0dHo2bNnveZGhvhB7lly69YtdO7cudUVo5cvXwafz8fu3buf6jqlpaVNNtJ58IvP50MkEjU4Q9THxwf379/X0qt6OqdPn8aLL74IS0tLmJubY9iwYXrpaK1UKsHn8zWNoZpzpMuD47p27Vqd56hNM2bMgLW1tUHug9+7dy9EIhHkcjkn8adMmYLOnTtr/d/1xm4WTZgwQfP4sGHD6j3m7+9f7+f1/d5Qq9UIDg5Gv3799Drjb6jjDwB79uxB165dIRAItJpbY1oy/qwQZRiGecakpqYiLCwMRkZGcHBwgEKhwPXr17lO65mVmZmJzp07o2fPnq2mSQwAvPDCCxg+fPhTXePAgQOaWU4jIyPweLx6RZK9vT38/f0xceJEfPDBB1i1ahUSEhIwffp0WFtbN+i+y5Xjx49DLBZjwoQJyM7Oxp07dzBjxgwIhUL8+eefOo+/aNEiWFtbIycnB5mZmTh69Cji4+OxatUqLFq0CNOmTcOIESPg5eUFGxubeuP8+uuv6zw/bYmNjQWPx8OOHTu4TqWBY8eOwczMDFOmTOFkufOpU6fA5/MRGxurk+vn5uZqZt03b97c4PG6Qufhzu1cvTfOnDkDPp+PrVu36izGgwx1/K9cuYJRo0ahZ8+esLCw0EshCjR//FkhyjAM84zKycmBUqmEi4uLprlRQkJCq9wz1tq1xmK0rlnMxYsXn/gaP/74I1588UW88847+OKLL7Bp0yYkJyfjxo0bTZ4JWVFRASsrKyxduvSJ42pTbW0tunfvDicnp3ozISqVCt26dYOLiwsqKyt1mkNJSQmsra3x1VdfNev5NTU1yMrKwokTJ3D48GGd5qYtV69ehaWlJSIiIrhOpYHLly/D2toaY8eO5axJTlhYGPz8/HT67/dPP/0EIoKtrS1yc3M13y8oKICLiwuOHDlS7/lcvzdee+019OnTR2fXf5Ahjj/w7xgsXboUNTU1cHZ21lshWhf7cePPClGGYZhnnEqlqtfcqEuXLlAqlQZ7Ll9blZmZCXd3d/j6+raKYrS2thZdu3bVe8OYzZs3QygUGsxS5oMHD4KI8O677zZ4bPHixSAibN++Xed5zJo1C926dWuTN5Kqq6sxYMAA+Pj4GNx2gpKSEnTv3h1+fn6c5VZcXAxTU1O9LLMePXq05ozfOpMmTcIHH3zQ4LlcvzeSkpJARDhz5ozOYgCGO/4A6v1N6rsQbc7484lhGIZ5pgkEAho1ahQlJibSP//8Q+PGjSOlUknOzs40ZcoUOnPmDNcpPhNcXFzo4MGDVFpaSlKplO7du8d1So/E5/Np7ty5FBMTQzk5OXqLu3XrVnrxxRepXbt2eov5KAcOHCAioj59+jR4rO57+/fv13ke06ZNo0uXLrXJ9+uCBQvo7NmzFBsbS2KxmOt0NNRqNU2aNInu3btHu3bt4iy3Xbt2kVqtpokTJ+o81tq1a8na2pp++eUX2r59O8XHx9OZM2doyZIlDZ7L9XsjMDCQOnfuTFu3btVZDCLDHX8i4vT90pzxZ4UowzAMo9GtWzdSKpWUmZlJK1asoDNnzlDv3r2pT58+tG7dOqqsrOQ6xTatY8eOdPDgQSopKWkVxeibb75J1tbWtHr1ar3EA0BHjx6l4OBgnVy/qKiIeDxeva/PPvuMiIhUKlW977/66qtERHTx4kUiIurQoUOD6zk7OxMR0eXLl3WS74P8/PzIxsaGUlJSdB5Ln/bu3UvffvstrVmzhjw9PblOp56PPvqI/vjjD4qNjW30968vycnJ1LdvX7K0tNR5LCcnJ1q+fDkREYWHh1NERATFxMSQSCRq8Fyu3xs8Ho+GDRum8/eEoY4/15oz/qwQZRiGYRqQSCQkk8no7NmzlJycTO7u7hQREUFubm704Ycf0s2bN7lOsc3q2LEjJSUlUXFxMUmlUiooKOA6pSaJRCKaNWsWrVmzhkpLS3Ue7+LFi1RQUEADBgzQyfWtrKwIAI0YMYL4fD5duXKFPvroIyIiEgqFBIAGDBhAP//8M23fvp2I/i1eiYjMzMwaXM/c3JyIiAoLC3WS74N4PB7179+fjh49qvNY+lJYWEjTpk2jCRMm0JQpU7hOp54dO3aQUqmkNWvW0KBBgzjN5fjx4zp7TzTmjTfeoBdffJHy8/OpS5cu5O/v3+jzDOG9MXDgQEpNTaWamhqdxTDU8TcEjxt/VogyDMMwj/T8889TbGws3bx5k+bMmUM//fQTubu706hRo2jfvn0EgOsU25y6mdGioiKDL0bDw8OppqaG1q9fr/NYdTdAunXrptM477//PqnVavr222/rff/IkSOUlZVFoaGhzbpO3XuDx+NpPcfGeHh4tKmbRLNmzSIAeptxb660tDQKCwujiIgImj59Otfp0I0bN3T+nnhY165diYjo0KFDtGvXrhb/vL7eGx4eHlRZWUl5eXk6i9Eax19fHjf+rBBlGIZhmsXJyYkUCgVdvXqVtm7dSpWVlRQcHEyenp705Zdf6mXW51ni6upK+/bto7t379ILL7xgsMt0bWxsaNq0afTNN9/odNaBiOju3btkbGysmU3RlWHDhlHv3r1p48aN9cb966+/ptmzZ5NQKNR8z8rKioiIysvLG1yn7nt1z9E1W1tbg/07aaktW7bQtm3b6IcffiAbGxuu09EoLi6mV199lfr27dvgRgUXqqurqaysjGxtbfUWMzk5mXbs2KF5/W+//Xaj//4bwnujblzu3r2rk+sb8vgbgseNPytEGYZhmBYxNjam0NBQTXOjESNG0GeffUaurq40c+ZMOnv2LNcpthmdO3em5ORkKioqokGDBlF2djbXKTVqzpw5lJeXRz/99JNO49y/f59MTU11GqPOvHnzqKKigtasWUNE/+5lO3z4cIMZsLp9i7dv325wjaysLCL6d1ZAH8zNzamsrEwvsXQpOzubIiIiaNasWTRy5Eiu09EAQG+99RaVlpbSzz//XO+GBFcqKioIgN6a0pSVldHUqVNp3bp1NGfOHBo5ciTl5ubSe++91+C5hvDeqLtp1VgxrA2GPP6G4HHjzwpRhmEY5ol5enrSihUrKCsri5YtW0ZHjhwhX19f6tOnD/344486nyF7Fri6utLBgwdJpVJRUFBQox/quObm5kZhYWH02WefkUql0lkca2trKikpodraWp3FqDNhwgRycXGhVatWUVVVFX3zzTc0Y8YMkkgk9Z43dOhQIiI6efJkg2vUfW/YsGE6z5eIqKCgQK8zM7oAgKZPn07W1takVCq5Tqee5cuXU3x8PG3evJnat2/PdTpERGRpaUlCoVBvM2Lz5s0jqVRKI0aMICKi6OhosrCwoE2bNtGvv/5a77mG8N6oWyGgq1l1Qx5/Q/DY8df1GTIMwzDMsyU5ORmhoaEQCoVo164dFAoFMjMzuU6r1cvNzYWPjw9cXV1x5coVrtNp4OrVqxAKhYiJidFZjAMHDoCIkJeXp7MYD1q2bBmICJ9//jksLCxw+/btBs+pra2Ft7c32rdvj/v372u+r1Kp4OXlBRcXl3rf16WZM2diyJAheomlK6tXr4ZQKMTx48e5TqWe48ePw9jYGEuXLuU6lQbs7e2xcuVKncf5448/0KlTJ5SUlNT7/rp160BEaN++PQoLCzXfN4T3xv79+0FEyM/P11kMQx3/h+n7HFHg8ePPClGGYRhGJ7KyshAZGQkHBwcIBAKEhIQgMTERarWa69RarYKCAvTr1w9OTk44f/481+k08Oabb6JLly6oqanRyfXv3LkDHo+HP/74QyfXf1hJSQksLS3B4/EwZcqUJp937NgxmJiYYOLEicjJycHdu3cxc+ZMCIVCveUKAH379sW7776rt3jaduXKFZibm2PRokVcp1LPvXv34Orqipdeeskg//0aMmQIpk2bptMYhYWFcHFxwcGDBxt9XCqVgogwderUet/n+r3x9ddfw8HBQacxDHn8H8RFIfq48WeFKMMwDKNTVVVViI2N1fxH2a1bNyiVykfeuWWaVlhYiAEDBsDR0RFpaWlcp1PPlStXIBQK8eOPP+osRteuXfHxxx/r7PoPmz9/PojosWN96tQpjBw5EhYWFjA3N0dQUBBSUlL0lCVQUVEBIyMj/Pzzz3qLqU01NTXo378/evfujaqqKq7T0aitrcXw4cPRsWNH3L17l+t0GrVw4UJ4enrq7PrOzs4gIs3XmDFjNI8VFhbWe6zua/ny5ZrncPneGDt2LF5++WWdxjDk8d+9e3ejjxMR/vvf/+os5zqPG39WiDIMwzB6c/LkSchkMpiZmUEikUAmk+Hs2bNcp9XqlJSUIDAwEHZ2dkhNTeU6nXreeOMNnc6Kzpw5Ez179tTJtVuzuLg4CAQCZGVlcZ3KE/n0008hEolw7tw5rlOpJzIyEkZGRjh69CjXqTSpbvnjP//8w3UqBqWsrAwSiUTny2bZ+DeuOePPmhUxDMMweuPn50fR0dGa5kYpKSnUs2dPev755ykuLo41N2omiURCv//+O/n7+1NQUBAdOnSI65Q0PvroI7px44bOOui+8cYbdPbs2UYboDzL1q9fT8OHDzeYJjotcfr0afr0009JqVRSjx49uE5H48CBA/TZZ5/RihUraMCAAVyn06QhQ4aQq6srbdiwgetUDMq2bduoqqqKJk6cqNM4bPwb15zx5wHsJHKGYRiGG2q1mg4cOEDr1q2jnTt3kr29PU2ZMoUiIiKoQ4cOXKdn8Kqrq2nKlCkUHx9PP//8M40dO5brlIiIaObMmbRnzx66fPmyTo5b6dGjB/Xq1Ys2b96s9Wu3RufPnydfX1+Ki4ujV155het0WqSmpob8/f3Jzs6O9u3bR3y+YcyR5ObmUu/evWnw4MG0detWrtN5rE8++YRWrVpFV65cIUtLS67T4ZxarSY/Pz/y8vKiLVu26DweG//6mj3+epufZRiGYZhHuH37NiIjI2Fvbw9jY2OEhoYiMTGR67QMnkqlgkwmg0AgwPr167lOBwCQl5cHiUSCr776SifXj4uLA4/HM7hlyVx56aWX4Ovri9raWq5TabFPP/0UpqamBtUJuqamBs8//zw8PDxQXFzMdTrNUlhYCFtbWyxYsIDrVAxCTEwMBAIB0tPT9RKPjX99zR1/VogyDMMwBqWysrJecyNPT09ERUWhtLSU69QMmlKpBI/Hw7Jly7hOBQCwaNEiWFlZ6aTBi1qtxoABAxAQEACVSqX167cme/bsARFh7969XKfSYpcuXYKJiQm+/vprrlOpZ+7cuTAzM9NbEaMt3377LcRiMS5evMh1KpwqKChAhw4dMH36dL3GZeP/r5aMP1uayzAMwxiskydP0rp162jz5s0kFApp4sSJJJfLqXv37lynZpC+++47mj17Nn3wwQekVCo5zaWsrIy6dOlCYWFh9PXXX2v9+ufPn6e+ffvSggULaNGiRVq/fmuQn59Pvr6+FBgYSNu2beM6nRYBQFKplO7du0cnTpwgIyMjrlMiIqKtW7fS66+/Tj/++CNNnjyZ63RaRKVS0aBBg6i6upqOHTtGxsbGXKfEiddee40OHjxIaWlp5OjoqLe4bPz/1aLx13VVzDAMwzBPq6ioCFFRUXB3dwcRISAgALGxsTrrzNqa/fDDDxAIBIiIiOB8qebKlSthYmKCGzdu6OT63333HYRCYaucDXxaVVVVGDp0KLp06dLggPvWIDo6GkKh0KCWV586dQqmpqaYO3cu16k8sYyMDEgkErz99ttcp8KJqKgo8Pl87N+/n5P4bPxbNv6sEGUYhmFajdraWiQmJiI0NBQCgQDt27dHZGQk8vLyuE7NoGzfvh0ikQjjx49HZWUlZ3lUVVWhS5cumDRpkk6ur1arMXnyZEgkEpw4cUInMQxRbW0txo8fD0tLS5w5c4brdFosJycH1tbWmD9/PtepaOTm5sLFxQXBwcGt/gbXjh07IBAIsGTJEq5T0autW7eCz+fjyy+/5DQPNv7NH39WiDIMwzCt0pUrV6BQKGBnZ8eaGzXi4MGDsLKywpAhQ1BYWMhZHjt37gSPx0NSUpJOrl9dXY2RI0fC3t7+mShGa2pqMHXqVJiYmODgwYNcp/NExo4dCzc3N5SVlXGdCoB//4YCAwPRqVMn3Llzh+t0tCI6Oho8Hk9nDcMMzfbt22FsbIw5c+ZwnQoANv7NxQpRhmEYplWra240cOBAEBF69+6N6Ohog/mQy6X09HS4uLige/fuyMzM5CyPF198ET169NDZTFNZWRlGjBgBc3Nz/P777zqJYQjKysrw0ksvwczMDHv27OE6nSfy66+/gojw559/cp2KxvTp0yGRSFpdc6LHWbFiBfh8Pt577z3Ol+nr0sqVK8Hn8xEeHm5Qr5ON/+OxQpRhGIZpM1JTUyGTySAWi2FpaQmZTIbz589znRansrKy4Ovri/bt23O2jDMjIwMikQgrV67UWYzq6mpMnToVQqEQX3zxRZv74PfPP//A19cX9vb2+Ouvv7hO54kUFxejQ4cOeOONN7hOReObb74Bn8/H7t27uU5FJ2JjYyESiTBy5Mg2t4WhoqICMpkMPB4PS5cu5TqdBlJSUtCvXz8YGRmx8W8CK0QZhmGYNqewsBBRUVFwc3MDn8+HVCp9ppsbFRQUIDAwEFZWVjh06BAnOXz44YewsLBATk6OzmKo1Wp88803MDY2hlQqxe3bt3UWS1/UajXWr18PMzMzPPfcc7h27RrXKT2xt99+G3Z2dsjPz+c6FQDA77//DoFAAKVSyXUqOpWSkgJbW1s4ODi0mcZeZ8+eRffu3WFtbY0dO3ZwnY5GbW0tfvnlFwwYMABEhP79+2PVqlXo1KkTnJyc2Pg/hBWiDMMwTJtV19woJCQEPB4Pzs7OiIyMNJgPwvp0//59jB07FiYmJti+fbve45eXl8PV1RVvvfWWzmOdOHECXbt2hYWFBb799ttWewPi/PnzGDJkCPh8PubPn4/q6mquU3pix44dA5/Px08//cR1KgCAv//+G+bm5gY1O6sLf/75J7y8vCASidC/f38QEdDQ9BsAACAASURBVF5//XVkZ2dzndoTKS0txbx582BkZISAgACddeRuqcrKSsTExMDT0xN8Ph8hISH1ehYUFRVh/PjxbPwfwgpRhmEY5pmQkZEBhUIBW1tbiESiZ7K5kUqlQnh4OPh8PpYtW6b3+HVdFZOTk3Ueq6KiAosWLYKJiQl69OiBHTt2QK1W6zyuNty+fRsREREwMjJCnz598Pfff3Od0lOprq6Gt7c3Ro4cyXUqAP5tdObo6Ihhw4ahqqqK63R0IiMjA6GhoSAihISE4MqVKwCAXbt2wc3NDZaWlvj0009RXFzMcabNU1VVhe+//x7Ozs6wsbHB2rVrDWL5fX5+PpRKJZycnCASiRAWFoYLFy40+Xw2/vWxQpRhGIZ5pty/fx8xMTHo1asXiAj+/v6Ijo5GeXk516npTd1Zb2+99ZbeZ9lGjx4NDw8PVFRU6CVe3QdyPp8PX19fbNu2zWBnFjMyMhAeHg6RSAQXFxesW7fOID5sP61ly5ZBJBLh8uXLXKeC/Px8dO3aFX379kVpaSnX6WhdWVkZIiMjIRKJ4Onpid9++63BcyoqKvDJJ5/AysoKNjY2WLJkicHuXywpKcGqVavg4uICkUiEiIgIg+hsfP78ecycORNisRh2dnb4+OOPmz2GbPz/P1aIMgzDMM+sh5sbyeVyXL16leu09OK3336DRCKBVCrV6/EuWVlZsLKygkKh0FtMADh37hzGjx8PgUCAdu3aYcGCBcjIyNBrDo25f/8+4uLiIJVKwePx4OrqijVr1rSZmbqcnBxYWloiMjKS61RQUlICPz8/dO7cGbm5uVyno1VqtRoxMTFo164drK2tERUV9dgl6YWFhVi8eDFsbGxgbGyM8ePHY9++fVCpVHrKumknTpzAzJkzIZFIIBaLER4ejlu3bnGak0qlwo4dOxAUFAQej4euXbtizZo1T3wTk40/K0QZhmEYBnl5eVAqlXB1da3X3MgQPhDoUlpamuZ4l+vXr+st7tq1ayEQCDhZcpqZmYlFixahffv2ICL4+vpiyZIlSEtL09vS3cLCQsTGxmLixImQSCQQCAQICQlBQkJCm/ubmzRpEjp27Mj5cUrV1dV44YUX4ODgYBAzs9r0999/Y8CAARAKhZDJZC3eA19RUYGNGzdqGuw4ODhAJpPhjz/+0NvKhZqaGqSkpOD999+Hu7s7iAheXl6IiopCQUGBXnJoSlPN77T1Xn2Wx58HAMQwDMMwDKnVajpw4ACtWLGC9uzZQ+7u7jRjxgx66623yM7Ojuv0dCI7O5tGjx5NN27coJ07d9KgQYN0HhMADR8+nG7fvk2nT58mkUik85gPU6lUlJSURDt27KD4+HjKyckhW1tbev755ykwMJB8fX3Jx8eHHBwcnipOdXU1Xbp0idLT0+n48eN06NAhOnfuHPF4PBo8eDCNHTuWXnnlFWrfvr2WXpnhSElJocDAQNq5cyeNGTOGszwA0BtvvEG7du2ipKQk6t27N2e5aFNWVhYtWLCANm/eTEOGDKGoqCjq2bPnU13z4sWLtGPHDtqxYwedPHmSjI2NqW/fvhQYGEh9+vQhHx8fcnd3J4FA8FRxbt68SefPn6eTJ09ScnIyHT16lMrLy8nDw4PGjh1L48aNoz59+jxVjKd16tQpio6Ops2bN5NQKKSJEyfS7NmzycvLS2cxn7XxZ4UowzAMwzQiIyOD1q9fT//973+pvLycRo8eTbNnz6aBAwdynZrWlZWV0WuvvUb79++njRs30vjx43Ue8+rVq+Tr60vz5s2jTz75ROfxHkWtVlNaWhodPnyYDh06REeOHKH8/HwiIrK3t6euXbuSk5MTOTs7k6OjI1laWpKRkRGZmZmRsbExlZSUUG1tLZWUlFBRURHdvn2bcnNzKTMzkzIyMkilUpGRkRH17NmTBg0aRIMHD6bAwECysbHh9HXrUm1tLfn7+5ODgwPt3buX01zmzJlD33//Pf32228UFBTEaS7acP/+ffruu+/o888/JysrK/rss89oypQpWo+TlZVFSUlJdPjwYTp8+DBdvnyZ1Go1icVi8vDwIBcXF3JycqL27duTpaWl5v1gZmZGlZWVdP/+faqsrKTi4mLKzc2l27dvU05ODl26dIlKSkqIiMjV1ZUGDRpEgYGBNHjwYPLw8ND662iJyspK+uWXX2jVqlV0/Phx6tmzJ4WHh9PkyZPJ1NRUr7k8C+PPClGGYRiGeYTKykqKjY2lb7/9ltLS0sjf359kMhmFhYWRWCzmOj2tqa2tpXnz5tF3331HCxcupCVLlhCfz9dpzBUrVtD8+fPp0KFDNGDAAJ3Gaqm8vDxKT0+n9PR0un79OmVnZ1N2djbl5eVRaWkpVVdXU1lZGdXU1JBEIiGhUEiWlpYkkUjIxcWFHB0dycXFhRITE6myspL++usvMjIy4vpl6c2KFSvogw8+oLNnz1K3bt04y+OTTz6hJUuW0JYtW/Ryg0XXdu/eTe+99x7l5+fT+++/Tx9++CGZmJjoJXZFRQVduHCBzp07R5cvX9a8J7Kzs6m0tJTKysqourqaysvLycTEhMRiMYnFYpJIJJr3Q7t27ahr167Uo0cP6t69O1lZWekl98dJT0+nH374gTZt2kQlJSX0yiuvUHh4OA0ePJjr1DTa5PjrbNEvwzAMw7QxqampCAsLg5GREaysrCCXy/W6t1IfNm3aBLFYjJEjR+q8iZFarcaoUaPg4uKCe/fu6TQWV/bu3QsieuSRDm1NXl4erKyssHDhQk7zWL58OXg8HtauXctpHtpw+vRpBAYGgsfjITQ0FDdv3uQ6pVbv/v37iI2NhVQqBRHBxcUFCoWCja0e6fZWJ8MwDMO0If7+/vTjjz/SrVu36MMPP6T4+Hjq3LkzBQcH0+7duwltYJHR5MmTKSUlhc6fP099+/al8+fP6ywWj8ej9evXU21tLc2YMUNncbg0bNgwcnV1pQ0bNnCdit4oFAoyNzenBQsWcJbDqlWraM6cOfT111/TzJkzOcvjaRUUFNB7771Hffr0oYqKCkpJSaHY2Fjq2LEj16m1WidPnqSZM2eSg4MDhYWFkbW1NSUmJtLNmzdJqVSysdUnrithhmEYhmmtamtrkZCQoDl6o0uXLlAqlbh79y7XqT21/Px8DBkyBObm5vjll190GispKQkCgQCrV6/WaRyuREZGwtHR0WDPL9Wmo0ePgsfjIS4ujrMcNmzYAD6fj6VLl3KWw9Oqrq5GdHQ07Ozs0L59e0RHR7eJM2W5UlRUhOjoaPTu3RtEBE9PTyiVyhZ3GGa0ixWiDMMwDKMFly5dgkKhgLW1NUxMTBAWFobTp09zndZTqaqqwowZM8Dj8aBQKHT6Qfjjjz+GSCTCqVOndBaDK5mZmRAIBIiPj+c6FZ2qra1F3759MWzYMM5yiIuLg0AgMIhzS59UYmIiunfvDmNjY8jlcpSUlHCdUqtVd1a0qakpTExMEBoaisTERK7TYv4Pa1bEMAzDMFpUWlpKW7ZsoVWrVtG5c+faRHOjdevW0bvvvkvBwcG0efNmnTS4UKvVFBwcTLdu3aKTJ0+SRCLRegwuvfDCC2RiYkIJCQlcp6Iz0dHR9O6771JaWppOj7hoys6dO2n8+PEUERFBy5cv13v8p5WRkUH/+c9/KC4ujkJCQmjFihXk7u7OdVqtTl5eHm3dupV++OEHSk9P1/wb/Nprr7W5f1daO1aIMgzDMIyOnDx5klasWEFbt24la2trevPNN+mdd94hV1dXrlNrsZSUFAoNDSVra2vauXOnTjqhZmVlUa9evejFF1+kmJgYrV+fS9u2baPJkyfTzZs32+SZoaWlpeTh4UETJ07kpAjcu3cvjR49mmbMmEErV67Ue/ynUVZWRsuWLSOlUkmdOnWi5cuX04gRI7hOq1WpOwN63bp1FB8fT6ampjRhwgR6++2328y5sW0StxOyDMMwDNP25eTkQKlUwsXFBXw+H1KpFAkJCVCr1Vyn1iK3bt1Cv379IJFIsHXrVp3E+O233yAQCPDtt9/q5Ppcqaqqgp2dXavet/goCxYsgLW1NSf7o/ft2wcTExNMnTq1Ve2jrK2tRUxMDBwdHWFjY4OoqCjU1NRwnVarcvv2bSiVSnTq1AlEBH9/f0RHR6O8vJzr1JhmYIUowzAMw+hJVVWV5rgAHo+Hrl27QqlUtqqjSyorKyGXy0FECAsLQ0VFhdZjfPXVV+Dz+di9e7fWr82ld999Fx4eHq3uBsTj3Lp1C6ampli+fLneYx89ehTm5uZ49dVXoVKp9B7/Sf3111/o378/hEIhZDIZ7ty5w3VKrUZFRQU2b96M4OBg8Pl8tGvXDgqFApcvX+Y6NaaF2NJchmEYhuHAxYsX6fvvv6f//e9/pFKpKDQ0lObNm0e+vr5cp9Ys8fHx9Oabb1KnTp0oNjaWunTpotXry2Qy2rp1Kx09epR69Oih1Wtz5dy5c9SzZ086fPgwDRo0iOt0tGbSpEl0/PhxunDhAolEIr3FPX36NA0bNowGDx5McXFxJBQK9Rb7SWVlZdGCBQto8+bNNHToUIqKiiIfHx+u02oVjhw5Qhs3bqS4uDiqqKigkSP/X3v3Hpfz3f8B/HWdKoeI0IiZsZzCKIfmNGQOd90UMSOaTfttKLPRztybbRe7EdscMofCRplDYRvJTZgzuathDRsVpqPq6nB1vX9/7KefJofq6vqG1/Px6J9c38/n9e3x2B69+ny/n88QvPzyy/jHP/4BnU6ndDyqABZRIiIiBWVnZ2PDhg348ssvSzbWCAgIwJgxY6r9L1e///47Ro0ahV9++QUrVqzA6NGjzTZ2UVERBg4ciD/++ANHjhxBw4YNzTa2klxdXeHs7Iw1a9YoHcUsjh49ih49emDTpk3w9va22LxnzpxB//794eLigsjISIsW4IowGAxYvHgx5syZgyeeeAKfffYZfHx8lI5V7SUnJ2PdunVYuXIlfv31V7Rr1w7jx4+Hn58fHBwclI5HlcQiSkREVE0cOHAAixcvxpYtW2Bvbw8/Pz+88cYb1fqA9YKCAsycOROLFy+Gv78/Fi9ebLZSkJaWhm7dusHR0RHR0dGwsrIyy7hKWrJkCd5++22kpKRUye7DltanTx8YjUYcPHgQKpXKInOeP38effv2Rfv27bF9+3bY2NhYZN6KioqKQkBAAP7880+8/fbbePfdd6t9cVZSfn4+oqKiEBYWhh9//BG2trbw8fGBr68vevXqpXQ8MiMWUSIiomomNTUVYWFh+Oqrr5CSkoKhQ4ciMDAQAwYMsNgv++W1efNmTJw4ES1btkR4eDhatmxplnETExPh5uaGESNGYNWqVWYZU0lZWVlo0qQJgoODMWnSJKXjVMqmTZswatQoHDlyBF27drXInElJSejbty+aN2+OXbt2oXbt2haZtyJOnTqFwMBAHDhwAOPGjcO8efPwxBNPKB2r2jpx4gTCwsKwbt06ZGVloV+/fvD398ewYcMeiT9C0Z3USgcgIiKi0ho3boygoCD89ttv2LBhA/Lz8zFw4EC0adMGc+fORUZGhtIR7+Dt7Y0jR47AaDTC1dUVW7ZsMcu47dq1Q1hYGEJDQzFv3jyzjKmkunXrwsvLCytXrlQ6SqUUFhbivffew9ixYy1WQi9fvoyBAwfCwcEBO3bsqLYlNC0tDYGBgejatSvy8/Nx8OBBhIWFsYSWITk5GXPnzoWTkxNcXV0RHR2NmTNnIjk5Gbt374aPjw9L6KNMuX2SiIiI6EElJiZKQECA1KpVS2xtbcXf31/OnDmjdKw75OXlyauvvioAxN/f32zHKCxatEhUKpWsXLnSLOMpKSYmRgBIXFyc0lEq7N///rfY2NjIpUuXLDLf1atXpXXr1tKxY0dFjoh5EIWFhRIcHCx169aVJk2aSGho6CO3Q7I5GAwGCQ8PFw8PD9FoNFKvXj3x9/eX48ePKx2NLIyP5hIRET1Ebm1utGjRIiQmJlbbzY02b94Mf39/NGjQAN9++y26dOlS6TE/+ugjfPbZZ9iwYQNGjhxphpTKEBE888wz+Oc//4kFCxYoHafcMjIy0KpVK7z++uuYM2dOlc/3559/4vnnn4fRaMS+ffuq5cpidHQ0pk2bhgsXLiAgIADvv/8+bG1tlY5VrZw4cQIhISH47rvvkJeXV/Lo7fDhw6vV/7vIcvhoLhER0UOkTp068Pf3R3x8PHbv3o2nn34ar7zyCp588km88847uHz5stIRAfz1qG58fDxatGiB7t27Y/bs2TCZTJUa8+OPP8aUKVMwduxY7Nq1y0xJLU+lUsHPzw9hYWEoKChQOk65zZkzB1qtFkFBQVU+V2ZmJgYPHozCwkLs3bu32pXQ8+fPw9PTEwMHDkSLFi2QkJAAvV7PEvp/zp8/j1mzZqFVq1ZwdXXFsWPHMGfOHFy9erXk0VuW0McXV0SJiIgecikpKQgJCcGSJUuQnp6OIUOGVJvNjUQEixcvxsyZM9GzZ0+EhYWhadOmFR7PZDJh7Nix2LFjB/bu3QsXFxczprWc5ORkNG/eHBs3bsSIESOUjvPAkpOT8cwzz2DevHmYMmVKlc6VnZ0Nd3d3XL16Ffv378dTTz1VpfOVR2ZmJvR6PYKDg/H0009j4cKFGDRokNKxqoXr169j48aNWL9+PY4cOYImTZrgxRdfxIQJE9CxY0el41E1wiJKRET0iCgsLMS2bdsQEhKC6OhotG7dGq+//jpeeeUVxTd2OXHiBMaOHYurV69i2bJlePHFFys8VlFREYYNG4Zjx44hNjYWbdq0MWNSyxk6dCgAYOfOnQoneXATJ07E3r17cfbs2So9giQvLw9DhgzBuXPnsG/fPrRu3brK5ioPk8mEdevWYebMmSgqKsJHH32EKVOmQKPRKB1NUfn5+di9ezfWrl2LrVu3QqvVwsPDA76+vhgyZAi0Wq3SEakaYhElIiJ6BJ08eRLLly/HunXroNVq8eKLL2Lq1KlwdnZWLJPBYMA777yDxYsXw9fXF0uWLKlwQc7NzYW7uztSU1Oxb98+NG/e3Mxpq97333+PUaNG4eLFi9X6rNhbzp07B2dnZ6xevRrjxo2rsnkKCgowfPhwnDhxAv/5z3/Qrl27KpurPPbt24dp06YhISEBL7/8Mj799FM0aNBA6ViKMZlMOHToENauXVvqvU9fX194e3sr/scvqv5YRImIiB5hWVlZ2LhxI4KDg/HLL7+gZ8+eCAwMhJeXl2KrFFu2bMGkSZNgb2+PsLAwdO/evULjpKWlYcCAAcjOzsbevXsfujJaWFiIpk2bIiAgAB988IHSce5rxIgROH/+POLi4qBWV802I4WFhRgxYgQOHDiAPXv2mGWTq8q6cuUK3nvvPaxbtw79+/dHcHCwon/QUVpCQgLWrl2LsLAwpKamol27dhg/fjz8/Pzg4OCgdDx6iHCzIiIiokdY3bp1S21u1KRJE4wZM6Zkc6MrV65YPJOXlxfi4uLw1FNPoVevXnj//fdRWFhY7nHs7e2xd+9e2Nvbo0+fPrhw4UIVpK06VlZWGDduHFatWlXpjZyq2rFjx7BlyxZ89tlnVVZCi4qKMHr0aOzfvx+7du1SvITm5eVh7ty5aNu2LQ4fPoyNGzciOjr6sSyhly9fxqJFi9C5c2c4Ozvju+++w/jx43Hu3DkkJCQgKCiIJZTKjSuiREREj5kLFy4gJCQEK1euRHZ2NoYNGwZ/f3+4u7tbPEtYWBimTJmCJ598EmFhYRUqH5mZmXjhhRdw7do1xMTEoGXLllWQtGokJCTA2dkZMTEx6Nevn9Jx7srd3R05OTn4+eefq2QDrOLiYvj6+iIyMhI7d+5Enz59zD5HeURFRWHq1Km4ceMG3n77bbz77rtV+k5sdZSZmYnIyEisXbsWe/bsgZ2dHXx8fODr64uePXsqvhEaPfy4IkpERPSYefrpp6HX63HlyhWsW7cOKSkpGDhwINq1a4dFixYhNzfXYlnGjx+PM2fOoFGjRujevTveeecdFBUVlWsMOzs77Nq1Cw4ODujXrx9+++23Kkprfu3bt0f37t2xcuVKpaPc1a5du7Bnzx7o9foqK6ETJkzAtm3bsH37dkVL6MmTJ9G7d28MHz4cffr0wW+//YbZs2c/NiW0oKAAUVFRGD9+PBwdHfHaa6/BxsYGGzduxNWrV7F8+XL06tWLJZTMgiuiREREVHLY/O2bGwUGBlpsoxgRwYoVKzB9+nS0atUKoaGh6NSpU7nGyMzMxKBBg5Camoq9e/c+NCujISEhCAwMREpKCurVq6d0nFJEBG5ubmjQoAG2b99u9vFNJhP8/PywadMmbN++Hf379zf7HA8iLS0NH3/8Mb7++mu4uLhg0aJF6NGjhyJZLK24uBgxMTH47rvvsGXLFty8eRP9+vXDuHHj4O3tzTNRqeoIERER0f/JyMiQ4OBgadGihahUKnF3d5fw8HApKiqyyPxJSUnSq1cvsbGxEb1eL0ajsVzXp6eni4uLizRr1kzOnTtXRSnN6+bNm1K7dm1ZsmSJ0lHusHHjRlGr1XLy5Emzj20ymeS1114TKysr2bFjh9nHfxCFhYUSHBwsdevWFUdHRwkNDRWTyaRIFksymUxy8OBBmTJlijg4OAgAcXV1lfnz50tycrLS8egxwRVRIiIiuoPJZEJMTAxCQkKwefNmODg4YNKkSZg8eTIaNmxYpXMbjUbo9Xp88skn6N69O1avXl2u1c2MjAwMHToUFy5cwA8//KD4pjcPws/PD/Hx8Th+/LjSUUoUFxejffv2cHFxwfr16806tohg8uTJWLlyJb7//nt4eHiYdfwHER0djcDAQFy8eLFk5+JH/ciRhIQEREREYP369UhKSkLbtm0xatQojBkzptqc1UqPEYWLMBEREVVzSUlJEhQUJPb29mJlZSU+Pj6ye/fuKp/31KlT0rFjR6lZs6YsWLCgXKujOTk5MnjwYKldu7bs2rWrClOax/79+wWAnDp1SukoJcLCwkSr1cr58+fNOq7JZJLJkyeLTqeTbdu2mXXsB3H27FkZOnSoABAPDw+5ePGixTNY0qVLl0Sv10vbtm0FgDRr1kwCAgIkNjZW6Wj0mGMRJSIiogeSn58v4eHh4ubmJgCkS5cusnz5csnNza2yOQsLC0Wv14u1tbX06NFD4uPjH/jagoICGT16tFhbW0tERESVZTSXNm3ayNSpU5WOISIiRqNRWrduLX5+fmYfe+bMmaLRaGTDhg1mH/teMjIyJCgoSKysrKRz586yb98+i85vSVeuXJHg4GDp2bOnqFQqsbe3F39/f4mNjX0sHj2mhwOLKBEREZXb8ePHxd/fX2rUqCF169YVf39/SUxMrLL54uPjpXv37qLT6SQoKEgKCgoe6DqTySTTpk0TjUYjK1asqLJ85vD5559L/fr1xWAwKB1F1qxZIxqNxuzv2b777rui0Wjk22+/Neu491JcXCyhoaHSqFEjsbe3l+Dg4HK/e/wwSEtLk9DQUPHw8BCNRiN2dnbi6+srkZGRUlhYqHQ8ojuwiBIREVGFXbt2TfR6vTz11FOiVqtLNjeqil/0i4uLZfny5VKrVi3p0KGDHD169IGv1ev1olKpZO7cuWbPZS6pqami0+ksvlL4d0ajUZycnGTixIlmHff9998XjUYja9euNeu497J3717p1KmT6HQ6CQgIkMzMTIvNbQm5ubkSHh4uHh4eotPpxMbGRjw8PCQ0NLRKn1QgMgcWUSIiIqq04uJi2b17t3h4eIhKpRJHR0eZNWuW/Pnnn2af67fffpMBAwaIVquVgIAAycnJeaDrFi5cKCqVSt566y0pLi42ey5z8PT0lIEDByqaYdWqVaLT6eS3334z25gfffSRqFQqCQkJMduY9/LHH3+Ir6+vABB3d/dyPdJd3eXm5kpERIR4e3uLjY2NWFlZiaenp3z77bcP/N8CUXXAIkpERERmdf78eQkKCpL69euLtbW1+Pj4yIEDB8w6h8lkkqVLl0qdOnXEyclJ9u/f/0DXrV+/XqytrcXb27tarhht3bpVVCqVWUtgeRQWFsrTTz8tr776qtnG/OKLL0SlUsmyZcvMNubd5ObmyqxZs8TGxkacnJwkKiqqyue0hLy8PPn+++9l9OjRUqtWLdFoNDJgwAD55ptvJD09Xel4RBXCIkpERERVwmAwSGhoqHTq1EkAiIuLi9k3N0pJSZHhw4eLSqUSX19fuXHjxn2vOXTokDRs2FC6du0qqampZstiDkVFRdK4cWOZNWtWqe8fPnzYrEWuqKhIrl+/fsf3v/nmG9HpdHLhwgWzzLNgwQJRqVTy9ddfm2W8uzGZTBIeHi7NmzcXOzs70ev1kp+fX6VzVrX8/HyJjIwUX19fqVOnjqjVaunZs6cEBwdLSkqK0vGIKo1FlIiIiKrcrc2NbGxsxM7OTgICAsxWdkREIiMjxdHRURwcHCQ0NPS+n09KSpI2bdqIo6NjtToyRURkxowZ0qxZM0lNTZUFCxaIk5OTABA3NzezzXHp0iWpUaOGTJs2TS5fviwi/78a+tprr5lljuDgYFGpVPLll1+aZby7OX78uPTs2VPUarX4+vrK1atXq3S+qnSv8pmcnKx0PCKzYhElIiIii7l69aro9Xpp3ry52Tc3ysjIkICAAFGr1dKvX7/77vialpYmffv2FVtbW9m5c2el5zcHo9EoISEhYm1tLRqNRrRarahUKgEgnTp1Mts8sbGxAkC0Wq1otVp55ZVX5NNPPxWdTmeWczVDQkJEpVKJXq+vfNi7SElJEX9/f1Gr1dK9e3c5fPhwlc1VlYxGo8TGxkpAQIA0aNBAAEi7du1Er9fLlStXlI5HVGVYRImIiMjiiouLJTIyUtzd3UWlUknLli1Fr9c/0KO193PgwAFp37691KhRQ2bNmnXPo17y8/PlpZdeEp1OJ99841EPbwAAGJxJREFU802l566oy5cvi16vF0dHx5KCCKDUV+vWrc0237fffltScAGITqcTlUolTk5Ocvz48UqN/c0334hKpZLPPvvMTGlLKywslODgYKlTp440bdpUQkNDH7qzMW8vnw0bNiwpn7NmzZKkpCSl4xFZBIsoERERKercuXMSFBQk9erVK9nc6NChQ5Ua81ZZuXXUy73GM5lM8uGHH4pKpZKpU6dKUVFRpeYur6VLl4pKpRKdTndH+bz968knnzTbnPPmzRMrK6s75rhVSIcMGSI///xzmdfGx8ff9f3LVatWiVqtlk8++cRsWW8XGRkpLVu2lJo1a0pQUJDcvHmzSuapCreXz0aNGpUqn7/++qvS8YgsjkWUiIiIqoXs7GxZvny5dOzYsdTmRnl5eRUe8/z58zJgwABRq9UyefLke54jGR4eLrVq1ZLevXvLtWvXKjxneRUWFkrfvn3vW0QdHBzMNmdgYGCZRfTW160V2VGjRt2x2ujm5iYDBw4Ug8FQ6vtr1qwRtVp9x0ZL5nD27FkZMmSIABAPDw+zPD5sCUVFRbJr1y6ZNGmS2NvbCwBxdXWVefPmPTT3QFRVWESJiIio2jl+/Lj4+vqKTqeTRo0aSVBQUIV/cTeZTBIWFiYNGzYUBwcHWbNmzV0f5YyLi5MWLVpIs2bNKv2IanlkZWWJk5NTmY/k3vqys7Mz23ze3t6lHs39+5dGoxFbW9s7fga33i3VaDTSv3//kj8ShIeHi1arlbfeestsGUVE0tPTJSAgQLRarXTp0uWBj+lR0u0rnw4ODqVWPu/33jLR44RFlIiIiKqt1NRU0ev10qxZs5LNjSIjIyv0TuCtzYw0Go107dpVjhw5Uubnbty4IQMGDBAbGxtZu3ZtZW/hgV24cEHq168vGo2mzHJYo0YNs8317LPP3rOE1q5dW44ePXrHdYMGDSopyzqdTnr27Cnr1q0TrVYrb775ptnyFRcXS2hoqDRs2FDs7e0lODjYLBtaVZX8/HzZvXt3me98nj17Vul4RNUSiygRERFVe0ajsdTmRq1atRK9Xi9paWnlHuvUqVOljvv4888/7/hMYWGhTJ48WVQqlbz//vtSXFxsjtu4r9jY2JL3NMsqiOZy6zHRsuaoVatWmSX99OnTd+TSarVSo0YNmTx5stmyxcTESMeOHUWn00lAQMA9H6dWksFgKDlqpW7duqXK5/nz55WOR1TtsYgSERHRQ+Xs2bMSFBQkdnZ2YmNjI76+vuU+C9RkMkloaKg4ODhI/fr1JTg4uMyyuWLFCrG2tpZBgwaVWVhvWbZsmdl2bt2wYcNdH5s1x6pgYWGhqNXqcpVQEZGRI0eW+R6rVquVHj16SHZ2dqVy/fHHH+Lr6ysAxN3dXRISEio1XlXIy8srKZ+2tralzvnkUStE5cMiSkRERA+lW5sbdejQodTmRn/fROdebj2uq9VqxcXFpcyzKE+cOCEtWrSQpk2blrn7blhYmACQr776qlL3c7sPP/ywzLKYk5NT6bEvXbpU7hKalJRUZp7bd9v9exl90GKek5Mjs2bNEhsbG3FycpLt27dX+h7NKTc3t6R81q5dWzQaTUn5TE5OVjoe0UOLRZSIiIgeerGxseLj4yM6nU4cHBwkKChILl269MDXnz59Wnr16lXyuO7169dL/fuNGzdkyJAhotVqRa/Xi8hfBSUpKUlq1qwpAMTa2loSExPNcj8mk0nGjBlzx+ZF5jpntawSWlYJv2XixIn33GX37yujJpNJ/Pz85MyZM/e8x/DwcHnyySfFzs5O9Hr9Pc98taSMjAwJDQ0VHx8fqVWrVqnymZqaqnQ8okcCiygRERE9MlJSUkSv10vTpk3LvbnRrWLUtGlTqVev3h0b5JhMJtHr9aLRaMTLy0sGDRokTz/9dKnNe5ydnc1WpgwGg3Tr1q3U47DmePzz9kd/H6SEXr58+Z67+d5eaOvXry+RkZHy4YcfCgDp169fmWMeO3ZMnnvuuZLib8njcu7m2rVrEhISIoMHDxYrKyuxtrYWDw8PWb16dYXeRSaie2MRJSIiokdOQUGBhIeHi7u7uwAQJycn0ev1kp6eft9rbz0qamVlJV26dLnjcdxdu3aVbE5TVhl77733zHYf165dK9kxGIAkJSVVesz58+eLTqd7oBIqIjJt2rR7nnF6q4Dq9XrJzc2V0NDQUu+4bt26tWSs5ORk8ff3F7VaLc8//7ycPn260vdTGRcuXJD58+dL7969Ra1WS82aNcXLy0vWrl1bbTdJInpUqEREQERERPSIOnv2LJYuXYpVq1YBAF566SVMnjwZHTt2vOd1586dw9SpUxEdHY1x48Zh7ty5aNy4MYqLi9G8eXMkJyeXeZ1KpcLevXvRt29fs+RPTExEt27dkJubi/j4eLRv3x5FRUW4fv06rl27hszMTJhMJmRlZcFkMqFWrVqwsrJCzZo1YWdnB0dHR9StW7dkvDfffBPBwcGoWbMm9uzZgx49etx17rS0NDRt2hT5+fl33KNKpUKjRo3wzjvvwN/fHzVq1MD+/fvh7u6OoqIiAIBarYajoyMSEhKwatUqfPTRR6hTpw4+/fRT+Pr6QqVSmeVnVB4JCQnYvn07oqKicOjQIdjZ2cHd3R0eHh7w8vKCra2txTMRPY5YRImIiOixkJ2djQ0bNmDx4sVISEiAi4sLAgICMGbMGOh0urteFxUVhcmTJyM1NRWBgYF44oknMGPGjLt+XqPRoHHjxkhISECdOnUqlTknJwdHjx5FREQEli9fDmdnZ/z555+4du0ayvMrXM2aNdG0aVO0bdsW586dw8WLFxEaGorRo0ff87oPP/wQer0eRqMRAEqKY5MmTTBjxgy89tprsLGxAfBXce/atStyc3NhMplKxtBoNKhXrx4MBgPeffddvPXWWyXXWILJZMKpU6cQFRWFjRs34uzZs2jYsCEGDx4MHx8fDBo0CFZWVhbLQ0R/YRElIiKix86BAwewePFibNmyBQ0aNMCECRMwefJkNGvWrMzPb9u2DcOHDwfw1yrf7UWrLDqdDmPGjEFoaGi5chUXFyM2NhY//PAD9u/fj+PHj8NoNKJhw4awt7eHs7MzevfuDUdHRzRu3BgODg6oV68eVCoV6tSpA41Gg5ycHBQVFcFgMCA9PR1XrlzB1atX8fvvvyMxMRE7duxAQUEBjEYjGjdujL59+6Jfv3745z//iSeeeKIky82bN+Ho6IibN29CrVZDRNCiRQvMmjULL730ErRabcln09LS4OrqiitXrpSU1ttZWVnh559/RpcuXcr186io4uJi/Pzzz4iIiMCmTZuQkpKCFi1awNPTEz4+PnjuueegVqstkoWIysYiSkRERI+tlJQUhISEYOnSpUhLS8OQIUMQGBiIAQMGlHps1NPTEz/++GOZJeteNm3ahBEjRtzzMyKCmJgYfPfdd9i2bRtu3LiBNm3aoF+/fujTpw/69OmDJk2aAADy8/MrvZr466+/onnz5jh+/Dj2799f8mUwGODm5oYRI0bA19cXq1atQlBQEACgTZs2mD17Nnx8fO4ocPn5+ejbty9OnTpV8kju31W0mJeHwWBAdHQ0IiIiEBUVhczMTLRr1w4+Pj7w9PSEi4tLlc1NROXHIkpERESPvcLCQmzbtg0hISGIjo5G69at8fLLL+O1115DdnY2WrRocd9V0L+7tUqZmJhYUiRvl5GRgdWrV2P58uU4f/48unbtihEjRsDLywtOTk7murUHYjAYsGvXLmzZsgXbtm2DwWCASqWCo6MjvvjiCwwfPrzM9zlFBGPHjkVERMR9S7pKpcLhw4fRrVu3ku+ZTKZKrUymp6dj+/bt2L59O3744Qfk5eXBzc0Nnp6e8Pb2xjPPPFPhsYmoarGIEhEREd3m5MmTWLJkCb777jtotVo4Ozvj2LFjd13tuxedTodevXphz549JUUuOzsbS5cuhV6vh8lkwosvvojXX38dzz77rLlvpULy8/OxYMECRERE4PTp0+jZsyeCgoLg6el5x2fff//9kvu4H7VajS5duuDo0aMwmUz45JNPoFKpMGvWrHLlu3z5Mn744QdERUXhp59+glqtRu/eveHh4YHRo0eXeryYiKovFlEiIiKiMmRkZCAkJAT/+te/YDAYyn29SqUq2VBo8eLF+J//+R8sXLgQc+bMgZWVFd566y1MmTKlWu/SeuDAAfzrX/9CdHQ0+vXrhy+//BLt27cHAKxevRoTJ0685/U6nQ4mkwnFxcXQaDRo3bo1FixYgE8++QQHDx5E27ZtkZiYeN8cFy5cQFRUFCIiInDo0CHUqFED/fv3h4+PD4YPH17pTaGIyPJYRImIiIjuYu3atZgwYUK5dqgF/n93WbVajeLiYuh0upIjX4KCgjB9+vRqXUD/7uDBg3jzzTdx+vRpBAYGlmxuVFxcXPIZrVaL4uJiiAhsbGzQtm1buLm54dlnn0Xnzp3h7OyM6Oho+Pr6Ijc3t2SF+eLFi3jqqafumDMhIQERERGIiIhAYmIi7O3tMXToUPj4+OCFF16AtbW1pW6fiKoAiygRERHRXbi4uCAuLq5U4boXjUYDnU6H/Px8qFQq2NraokePHti9ezcaNGiA2NhYtG7duopTVw2TyYSQkBAEBQUhJyen5HFcOzs7dO7cGV27dkXnzp3x7LPPwsnJqdS7n0ajEZ988knJ47i3rtVqtZg/fz4CAgJK7XS7efNmXLlyBc2bN8ewYcPg6emJ559/vtROvUT0cGMRJSIiIirDyZMny7XTqlarha2tLYYOHYqkpCScOXOmZNOff//735g2bVrJ6ujDqqCgAIGBgfjpp5+QkpKCTz/9FG+//fY9r7l06RJGjhyJuLi4OzY0UqvV6NSpEzp27Ijt27cjLS0NHTp0gJeXF7y8vKrNe7NEZH78sxIRERFRGZYsWQLgr7Kk0WhKVvhMJhOMRuMdj+sajUZkZmZi586dCAkJwYwZM6BWqzFt2rSSFcKH/exKa2trLFu2DMXFxZg9ezZmzpyJrKwsfPzxx2Xuqvv999/Dz8+v5NzSvzOZTIiLi4NOp0NQUBC8vb3RsmVLS9wKESmMK6JEREREZTh27BiuX7+OrKwsZGdnIysrCxkZGcjKykJWVhbS09ORkZGBzMxMZGdn4+bNm8jNzQXw1zuibdq0QWxsLOzt7RW+k6qzZs0aTJo0CRMnTsSyZctKyqjBYMDMmTPx1Vdfldq0qSxqtRphYWEYO3aspWITUTXAFVEiIiKiMnTt2rXc1yQlJeG5555DmzZtsGjRoke6hAKAn58f7O3tMWLECNSrVw96vR6//PILvL29kZSUBAD33ehJrVZjy5YtLKJEjxmuiBIRERGZQUZGBlxdXWFvb4+YmBjUrl1b6UgWs27dOowfPx6jR4/G5s2bYTQaH+hs0Vtq1qyJ9PR07oRL9BhhESUiIiIyg5EjR+LIkSM4efIkGjZsqHQcizIYDOjWrRvi4+Pv+9lb79yqVKqSHXSLiorwww8/YPDgwRZIS0TVAR/NJSIiIqqk1atXY8uWLdizZ89jV0IBoEaNGoiLi8OAAQOQnJyMn376CSKCmzdvoqioCJmZmSgqKkJOTg7y8vJQUFCArKwsGI1GZGVloaCgAAaDQenbICIL4oooERERUSVkZ2fjmWeewUsvvYSFCxcqHUdRV65cQZs2bTBr1izMmDFD6ThEVI093HuIExERESls7ty5MBqN+Oijj5SOorimTZti+vTp+Pzzz5Genq50HCKqxlhEiYiIiCrIYDBgyZIlmDFjBurVq6d0nGrh1kroypUrFU5CRNUZiygRERFRBW3evBk5OTnw8/NTOkq1YWtrizFjxmDlypX3PbqFiB5ffEeUiIiIqII8PT2hVquxbds2paNUKydOnICrqytOnjyJzp07Kx2HiKohrogSERERVYCI4NChQxg4cKBF583MzCw5+uTW15w5cwAARqOx1PdHjhxp0Wy3dOnSBfXr18eBAwcUmZ+Iqj8WUSIiIqIKOHv2LNLT0+Hm5mbRee3s7CAiGDRoENRqNZKSkvDBBx8AALRaLUQEbm5uWL9+PTZt2mTRbLeoVCr06NEDhw4dUmR+Iqr+WESJiIiIKuD3338HALRu3VqR+adPnw6TyYQFCxaU+v7Bgwfxxx9/wMfHR5Fctzg5OZX8jIiI/o5FlIiIiKgCbty4ASsrK9SuXVuR+V944QV06NABa9asQVpaWsn3v/jiC0ydOhU6nU6RXLfY29uXykVEdDsWUSIiIqIKMBgMqFmzpqIZpk2bhry8PCxZsgQAcP78ecTExMDf31/RXABQu3Zt5OTkKB2DiKopFlEiIiKiCqhXrx6ys7NRXFysWIaxY8fCwcEBX331FQoKCjB//nxMmDChWpxpmp6eDnt7e6VjEFE1xSJKREREVAH29vYwmUyKPn5qbW2NN954A9evX8f8+fOxfv16BAYGKpbndtevX2cRJaK7YhElIiIiqoAOHTpApVLh1KlTiuZ44403UKNGDXzwwQdwd3dHq1atFM1zy8mTJ9GhQwelYxBRNcUiSkRERFQBDRo0QKtWrRQ/oqRBgwYYN24cRATTp09XNMstBoMBp0+ftvjRNkT08GARJSIiIqqg/v37Y+vWrUrHgJubG1xcXNCnTx+lowAAduzYAZPJhL59+yodhYiqKRZRIiIiogqaMGECzpw5gxMnTiiaY9myZdVmNRQAVq5ciUGDBqFJkyZKRyGiaopFlIiIiKiC3Nzc0L59eyxcuNCi837zzTfw8vJCTk4Oli1bhoyMDIwaNcqiGe4mISEBu3fvxquvvqp0FCKqxlhEiYiIiCph9uzZ+Pbbby2+Krp161bUq1cPS5cuxYYNG6DVai06/90EBQXB2dkZw4YNUzoKEVVjKhERpUMQERERPaxEBD179oRarca+ffug0WiUjqSYnTt34h//+Ad27dqFgQMHKh2HiKoxrogSERERVYJKpcKKFStw8uRJfPbZZ0rHUcz169fxyiuvYNSoUSyhRHRfXBElIiIiMoMvv/wS06dPx86dOx+7IlZYWIjBgwfj8uXLOHnyJGxtbZWORETVHIsoERERkRmICMaPH49t27YhJiYGrq6uSkeyCJPJhDFjxuCnn37Cvn370KlTJ6UjEdFDgI/mEhEREZmBSqXCqlWr0KtXLwwdOhTHjx9XOlKVMxqNeOWVVxAZGYmtW7eyhBLRA2MRJSIiIjITnU6HiIgIuLi4oF+/fvjxxx+VjlRlcnNzMXz4cEREROD777/H888/r3QkInqIsIgSERERmVGtWrUQGRmJkSNHwtPTE59//jlMJpPSsczq7Nmz6NmzJ44ePYqYmBgMHTpU6UhE9JBhESUiIiIyM51Oh1WrVmHu3LmYPXs2Bg0ahOTkZKVjVZqIYNWqVXB1dYWNjQ2OHDmCbt26KR2LiB5CLKJEREREVUClUmH69Ok4ePAgfv/9d7Rr1w4LFy6E0WhUOlqFJCYmon///pg0aRLeeOMNxMbGokWLFkrHIqKHFIsoERERURVydXVFXFwcAgMD8d5776Fz587YsmULHpaDC5KTkzF16lQ8++yzyMnJweHDhzFv3jzodDqloxHRQ4xFlIiIiKiK1ahRAx9//DH++9//om3bthg5ciQ6d+6M8PBwFBUVKR2vTElJSZgyZQpatmyJbdu24euvv8aRI0fQtWtXpaMR0SOARZSIiIjIQlq1aoXw8HDExcWhdevWeOmll/Dkk0/ivffeQ1JSktLxkJ+fj02bNmHgwIFwcnLC9u3bsXDhQiQlJWHSpElQq/mrIxGZh0oeludCiIiIiB4xly9fxooVK7By5UqkpKSgU6dOGDFiBIYNG4YOHTpApVJVeYbMzEzs3r0bmzdvxo4dO5CXl4chQ4bA398fQ4cOhUajqfIMRPT4YRElIiIiUpjRaMR//vMfbN68GVu3bkVqairs7e3Rq1cv9OnTB506dUKHDh3QqFGjSs1TWFiIc+fOIT4+HocPH8a+ffvw3//+FyqVCn379oW3tze8vLzQpEkTM90ZEVHZWESJiIiIqhGTyYS4uDjs378f+/btw8GDB3H9+nUAQMOGDfHMM8+gcePGcHR0hIODA+rWrQudTodatWrBysoK2dnZKC4uRnZ2NjIzM3HlyhVcvXoVf/zxB5KSklBUVASdToeOHTuid+/e6Nu3L/r06YP69esrfOdE9DhhESUiIiKq5q5du4b4+HjEx8fj4sWLSElJQUpKCq5du4abN2+isLAQOTk5KCoqgq2tLbRaLerWrQtbW1s0a9YMDg4OaNasGdq2bYv27dujTZs23PWWiBTFIkpEREREREQWxa3PiIiIiIiIyKJYRImIiIiIiMiiWESJiIiIiIjIorQAIpQOQURERERERI+P/wWyU6Ces42lJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.view_model()\n",
    "from IPython.display import Image, display\n",
    "display(Image(filename=\"causal_model.png\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                    \n",
      "─────(Expectation(y|W3,W2,W1,X1,X0,W0))\n",
      "d[v₀]                                  \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W2,W1,X1,X0,W0,U) = P(y|v0,W3,W2,W1,X1,X0,W0)\n",
      "\n",
      "### Estimand : 2\n",
      "Estimand name: iv\n",
      "Estimand expression:\n",
      "Expectation(Derivative(y, [Z1, Z0])*Derivative([v0], [Z1, Z0])**(-1))\n",
      "Estimand assumption 1, As-if-random: If U→→y then ¬(U →→{Z1,Z0})\n",
      "Estimand assumption 2, Exclusion: If we remove {Z1,Z0}→{v0}, then ¬({Z1,Z0}→y)\n",
      "\n",
      "### Estimand : 3\n",
      "Estimand name: frontdoor\n",
      "No such variable found!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "identified_estimand= model.identify_effect(proceed_when_unidentifiable=True)\n",
    "print(identified_estimand)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Model \n",
    "First, let us build some intuition using a linear model for estimating CATE. The effect modifiers (that lead to a heterogeneous treatment effect) can be modeled as interaction terms with the treatment. Thus, their value modulates the effect of treatment. \n",
    "\n",
    "Below the estimated effect of changing treatment from 0 to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                    \n",
      "─────(Expectation(y|W3,W2,W1,X1,X0,W0))\n",
      "d[v₀]                                  \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W2,W1,X1,X0,W0,U) = P(y|v0,W3,W2,W1,X1,X0,W0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+W3+W2+W1+X1+X0+W0+v0*X1+v0*X0\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 11.138323583866043\n",
      "### Conditional Estimates\n",
      "__categorical__X1  __categorical__X0\n",
      "(-3.742, -1.16]    (-3.816, -0.439]      5.350980\n",
      "                   (-0.439, 0.172]       8.452173\n",
      "                   (0.172, 0.677]       10.231967\n",
      "                   (0.677, 1.246]       11.996779\n",
      "                   (1.246, 4.411]       14.906734\n",
      "(-1.16, -0.589]    (-3.816, -0.439]      5.981074\n",
      "                   (-0.439, 0.172]       8.985925\n",
      "                   (0.172, 0.677]       10.909247\n",
      "                   (0.677, 1.246]       12.577942\n",
      "                   (1.246, 4.411]       15.509218\n",
      "(-0.589, -0.1]     (-3.816, -0.439]      6.582177\n",
      "                   (-0.439, 0.172]       9.335388\n",
      "                   (0.172, 0.677]       11.190280\n",
      "                   (0.677, 1.246]       12.909711\n",
      "                   (1.246, 4.411]       15.754069\n",
      "(-0.1, 0.472]      (-3.816, -0.439]      6.810293\n",
      "                   (-0.439, 0.172]       9.691918\n",
      "                   (0.172, 0.677]       11.559779\n",
      "                   (0.677, 1.246]       13.262256\n",
      "                   (1.246, 4.411]       16.144426\n",
      "(0.472, 3.296]     (-3.816, -0.439]      7.313792\n",
      "                   (-0.439, 0.172]      10.285809\n",
      "                   (0.172, 0.677]       12.091017\n",
      "                   (0.677, 1.246]       13.846789\n",
      "                   (1.246, 4.411]       16.792039\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "linear_estimate = model.estimate_effect(identified_estimand, \n",
    "                                        method_name=\"backdoor.linear_regression\",\n",
    "                                       control_value=0,\n",
    "                                       treatment_value=1)\n",
    "print(linear_estimate) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EconML methods\n",
    "We now move to the more advanced methods from the EconML package for estimating CATE.\n",
    "\n",
    "First, let us look at the double machine learning estimator. Method_name corresponds to the fully qualified name of the class that we want to use. For double ML, it is \"econml.dml.DML\". \n",
    "\n",
    "Target units defines the units over which the causal estimate is to be computed. This can be a lambda function filter on the original dataframe, a new Pandas dataframe, or a string corresponding to the three main kinds of target units (\"ate\", \"att\" and \"atc\"). Below we show an example of a lambda function. \n",
    "\n",
    "Method_params are passed directly to EconML. For details on allowed parameters, refer to the EconML documentation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                    \n",
      "─────(Expectation(y|W3,W2,W1,X1,X0,W0))\n",
      "d[v₀]                                  \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W2,W1,X1,X0,W0,U) = P(y|v0,W3,W2,W1,X1,X0,W0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+W3+W2+W1+X1+X0+W0 | X1,X0\n",
      "Target units: Data subset defined by a function\n",
      "\n",
      "## Estimate\n",
      "Mean value: 15.15260356848542\n",
      "Effect estimates: [15.33955778 15.84879232 14.66910411 ... 18.4715256  14.39605376\n",
      " 17.57646996]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LassoCV\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "dml_estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     control_value = 0,\n",
    "                                     treatment_value = 1,\n",
    "                                 target_units = lambda df: df[\"X0\"]>1,  # condition used for CATE\n",
    "                                 confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\":LassoCV(fit_intercept=False), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=False)},\n",
    "                                               \"fit_params\":{}})\n",
    "print(dml_estimate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True causal estimate is 11.13817451301838\n"
     ]
    }
   ],
   "source": [
    "print(\"True causal estimate is\", data[\"ate\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                    \n",
      "─────(Expectation(y|W3,W2,W1,X1,X0,W0))\n",
      "d[v₀]                                  \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W2,W1,X1,X0,W0,U) = P(y|v0,W3,W2,W1,X1,X0,W0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+W3+W2+W1+X1+X0+W0 | X1,X0\n",
      "Target units: \n",
      "\n",
      "## Estimate\n",
      "Mean value: 11.037725784237395\n",
      "Effect estimates: [15.27390981 11.5084465  11.0123769  ...  7.07033426 11.24713532\n",
      " 10.45835842]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "dml_estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     control_value = 0,\n",
    "                                     treatment_value = 1,\n",
    "                                 target_units = 1,  # condition used for CATE\n",
    "                                 confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\":LassoCV(fit_intercept=False), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=True)},\n",
    "                                               \"fit_params\":{}})\n",
    "print(dml_estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CATE Object and Confidence Intervals\n",
    "EconML provides its own methods to compute confidence intervals. Using BootstrapInference in the example below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                    \n",
      "─────(Expectation(y|W3,W2,W1,X1,X0,W0))\n",
      "d[v₀]                                  \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W2,W1,X1,X0,W0,U) = P(y|v0,W3,W2,W1,X1,X0,W0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+W3+W2+W1+X1+X0+W0 | X1,X0\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 11.03107854161636\n",
      "Effect estimates: [15.23154757 11.4914362  11.00644985 ...  7.09605105 11.23996513\n",
      " 10.45151404]\n",
      "95.0% confidence interval: (array([15.20814535, 11.43427515, 10.97728669, ...,  6.92179496,\n",
      "       11.21199051, 10.40377708]), array([15.53654796, 11.69017454, 11.10728956, ...,  7.17041884,\n",
      "       11.35021865, 10.56899233]))\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LassoCV\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from econml.inference import BootstrapInference\n",
    "dml_estimate = model.estimate_effect(identified_estimand, \n",
    "                                     method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     target_units = \"ate\",\n",
    "                                     confidence_intervals=True,\n",
    "                                     method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\": LassoCV(fit_intercept=False), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=True)},\n",
    "                                               \"fit_params\":{\n",
    "                                                               'inference': BootstrapInference(n_bootstrap_samples=100, n_jobs=-1),\n",
    "                                                            }\n",
    "                                              })\n",
    "print(dml_estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Can provide a new inputs as target units and estimate CATE on them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[11.09338659 12.13565105 11.3207331  12.64976079 13.0817584  12.11993925\n",
      " 11.57045914 11.44465176 13.48497597 11.05031631]\n"
     ]
    }
   ],
   "source": [
    "test_cols= data['effect_modifier_names'] # only need effect modifiers' values\n",
    "test_arr = [np.random.uniform(0,1, 10) for _ in range(len(test_cols))] # all variables are sampled uniformly, sample of 10\n",
    "test_df = pd.DataFrame(np.array(test_arr).transpose(), columns=test_cols)\n",
    "dml_estimate = model.estimate_effect(identified_estimand, \n",
    "                                     method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     target_units = test_df,\n",
    "                                     confidence_intervals=False,\n",
    "                                     method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\":LassoCV(), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=True)},\n",
    "                                               \"fit_params\":{}\n",
    "                                              })\n",
    "print(dml_estimate.cate_estimates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Can also retrieve the raw EconML estimator object for any further operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<econml.dml.dml.DML object at 0x7f516ff479d0>\n"
     ]
    }
   ],
   "source": [
    "print(dml_estimate._estimator_object)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Works with any EconML method\n",
    "In addition to double machine learning, below we example analyses using orthogonal forests, DRLearner (bug to fix), and neural network-based instrumental variables. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binary treatment, Binary outcome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            X0        X1   Z0        Z1        W0        W1        W2  \\\n",
      "0     0.607055 -1.614792  1.0  0.455616  1.221222 -0.648731  1.326302   \n",
      "1     2.993224 -0.200924  0.0  0.311552 -0.804113  0.837046 -0.219313   \n",
      "2     1.991809  1.339323  0.0  0.177437  0.541844 -1.252206 -0.018509   \n",
      "3     1.375978  0.476735  1.0  0.549200 -1.052208 -1.242290  0.363920   \n",
      "4     0.802443 -0.904535  0.0  0.348527  1.370334 -0.536766 -0.011914   \n",
      "...        ...       ...  ...       ...       ...       ...       ...   \n",
      "9995 -0.376158 -0.719342  0.0  0.213180 -0.463657 -1.072515 -1.974591   \n",
      "9996  0.116404  2.395842  0.0  0.153815  1.254356 -1.060713  1.624262   \n",
      "9997  0.514493 -0.307350  1.0  0.048421 -0.141551  0.897432 -0.602098   \n",
      "9998  0.715935 -1.660456  0.0  0.895173 -0.653106  2.064832 -0.931449   \n",
      "9999  1.145919 -1.137355  1.0  0.068377  1.188003 -0.917376 -0.878206   \n",
      "\n",
      "            W3  v0  y  \n",
      "0    -0.496089   1  1  \n",
      "1    -0.455882   0  1  \n",
      "2     0.276731   1  1  \n",
      "3     0.027972   1  1  \n",
      "4    -0.133733   1  1  \n",
      "...        ...  .. ..  \n",
      "9995 -1.383749   0  0  \n",
      "9996 -0.553284   1  1  \n",
      "9997 -0.602913   1  1  \n",
      "9998  0.302548   1  1  \n",
      "9999 -0.008667   1  1  \n",
      "\n",
      "[10000 rows x 10 columns]\n"
     ]
    }
   ],
   "source": [
    "data_binary = dowhy.datasets.linear_dataset(BETA, num_common_causes=4, num_samples=10000,\n",
    "                                    num_instruments=2, num_effect_modifiers=2,\n",
    "                                    treatment_is_binary=True, outcome_is_binary=True)\n",
    "# convert boolean values to {0,1} numeric\n",
    "data_binary['df'].v0 = data_binary['df'].v0.astype(int)\n",
    "data_binary['df'].y = data_binary['df'].y.astype(int)\n",
    "print(data_binary['df'])\n",
    "\n",
    "model_binary = CausalModel(data=data_binary[\"df\"], \n",
    "                    treatment=data_binary[\"treatment_name\"], outcome=data_binary[\"outcome_name\"], \n",
    "                    graph=data_binary[\"gml_graph\"])\n",
    "identified_estimand_binary = model_binary.identify_effect(proceed_when_unidentifiable=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using DRLearner estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                    \n",
      "─────(Expectation(y|W3,W2,W1,X1,X0,W0))\n",
      "d[v₀]                                  \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W2,W1,X1,X0,W0,U) = P(y|v0,W3,W2,W1,X1,X0,W0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+W3+W2+W1+X1+X0+W0 | X1,X0\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 0.48099941411381597\n",
      "Effect estimates: [0.34133643 0.58394521 0.69402234 ... 0.46721726 0.34146445 0.41183014]\n",
      "\n",
      "True causal estimate is 0.3833\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegressionCV\n",
    "#todo needs binary y\n",
    "drlearner_estimate = model_binary.estimate_effect(identified_estimand_binary, \n",
    "                                method_name=\"backdoor.econml.drlearner.LinearDRLearner\",\n",
    "                                confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{\n",
    "                                                    'model_propensity': LogisticRegressionCV(cv=3, solver='lbfgs', multi_class='auto')\n",
    "                                                    },\n",
    "                                               \"fit_params\":{}\n",
    "                                              })\n",
    "print(drlearner_estimate)\n",
    "print(\"True causal estimate is\", data_binary[\"ate\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instrumental Variable Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/25\n",
      "10000/10000 [==============================] - 2s 158us/step - loss: 6.0581\n",
      "Epoch 2/25\n",
      "10000/10000 [==============================] - 1s 90us/step - loss: 2.5423\n",
      "Epoch 3/25\n",
      "10000/10000 [==============================] - 1s 88us/step - loss: 2.3626\n",
      "Epoch 4/25\n",
      "10000/10000 [==============================] - 1s 87us/step - loss: 2.3099\n",
      "Epoch 5/25\n",
      "10000/10000 [==============================] - 1s 87us/step - loss: 2.2638\n",
      "Epoch 6/25\n",
      "10000/10000 [==============================] - 1s 99us/step - loss: 2.2351\n",
      "Epoch 7/25\n",
      "10000/10000 [==============================] - 1s 94us/step - loss: 2.2353\n",
      "Epoch 8/25\n",
      "10000/10000 [==============================] - 1s 83us/step - loss: 2.2118\n",
      "Epoch 9/25\n",
      "10000/10000 [==============================] - 1s 84us/step - loss: 2.2095\n",
      "Epoch 10/25\n",
      "10000/10000 [==============================] - 1s 83us/step - loss: 2.1964\n",
      "Epoch 11/25\n",
      "10000/10000 [==============================] - 1s 83us/step - loss: 2.1899\n",
      "Epoch 12/25\n",
      "10000/10000 [==============================] - 1s 99us/step - loss: 2.1810\n",
      "Epoch 13/25\n",
      "10000/10000 [==============================] - 1s 110us/step - loss: 2.1768\n",
      "Epoch 14/25\n",
      "10000/10000 [==============================] - 1s 107us/step - loss: 2.1662\n",
      "Epoch 15/25\n",
      "10000/10000 [==============================] - 1s 97us/step - loss: 2.1600\n",
      "Epoch 16/25\n",
      "10000/10000 [==============================] - 1s 81us/step - loss: 2.1619\n",
      "Epoch 17/25\n",
      "10000/10000 [==============================] - 1s 89us/step - loss: 2.1534\n",
      "Epoch 18/25\n",
      "10000/10000 [==============================] - 1s 105us/step - loss: 2.1529\n",
      "Epoch 19/25\n",
      "10000/10000 [==============================] - 1s 96us/step - loss: 2.1530\n",
      "Epoch 20/25\n",
      "10000/10000 [==============================] - 1s 104us/step - loss: 2.1509\n",
      "Epoch 21/25\n",
      "10000/10000 [==============================] - 1s 91us/step - loss: 2.1388\n",
      "Epoch 22/25\n",
      "10000/10000 [==============================] - 1s 88us/step - loss: 2.1440\n",
      "Epoch 23/25\n",
      "10000/10000 [==============================] - 1s 87us/step - loss: 2.1422\n",
      "Epoch 24/25\n",
      "10000/10000 [==============================] - 1s 94us/step - loss: 2.1302\n",
      "Epoch 25/25\n",
      "10000/10000 [==============================] - 1s 92us/step - loss: 2.1334\n",
      "Epoch 1/25\n",
      "10000/10000 [==============================] - 2s 197us/step - loss: 9634.7153\n",
      "Epoch 2/25\n",
      "10000/10000 [==============================] - 1s 124us/step - loss: 6117.5493\n",
      "Epoch 3/25\n",
      "10000/10000 [==============================] - 1s 132us/step - loss: 6080.0172\n",
      "Epoch 4/25\n",
      "10000/10000 [==============================] - 1s 144us/step - loss: 6000.0626\n",
      "Epoch 5/25\n",
      "10000/10000 [==============================] - 1s 140us/step - loss: 5903.0761\n",
      "Epoch 6/25\n",
      "10000/10000 [==============================] - 1s 114us/step - loss: 5971.7064\n",
      "Epoch 7/25\n",
      "10000/10000 [==============================] - 1s 109us/step - loss: 5927.9321\n",
      "Epoch 8/25\n",
      "10000/10000 [==============================] - 1s 107us/step - loss: 6042.7108\n",
      "Epoch 9/25\n",
      "10000/10000 [==============================] - 1s 107us/step - loss: 5935.9766\n",
      "Epoch 10/25\n",
      "10000/10000 [==============================] - 1s 108us/step - loss: 5899.0229\n",
      "Epoch 11/25\n",
      "10000/10000 [==============================] - 1s 108us/step - loss: 5814.6864\n",
      "Epoch 12/25\n",
      "10000/10000 [==============================] - 1s 110us/step - loss: 5938.3598\n",
      "Epoch 13/25\n",
      "10000/10000 [==============================] - 1s 110us/step - loss: 5783.8304\n",
      "Epoch 14/25\n",
      "10000/10000 [==============================] - 1s 110us/step - loss: 5861.9487\n",
      "Epoch 15/25\n",
      "10000/10000 [==============================] - 1s 122us/step - loss: 5884.3656\n",
      "Epoch 16/25\n",
      "10000/10000 [==============================] - 1s 145us/step - loss: 5875.7650\n",
      "Epoch 17/25\n",
      "10000/10000 [==============================] - 1s 131us/step - loss: 5842.4607\n",
      "Epoch 18/25\n",
      "10000/10000 [==============================] - 1s 117us/step - loss: 5833.5643\n",
      "Epoch 19/25\n",
      "10000/10000 [==============================] - 1s 123us/step - loss: 5883.4645\n",
      "Epoch 20/25\n",
      "10000/10000 [==============================] - 1s 110us/step - loss: 5888.1516\n",
      "Epoch 21/25\n",
      "10000/10000 [==============================] - 1s 111us/step - loss: 5883.6870\n",
      "Epoch 22/25\n",
      "10000/10000 [==============================] - 1s 111us/step - loss: 5917.5608\n",
      "Epoch 23/25\n",
      "10000/10000 [==============================] - 1s 114us/step - loss: 5954.1094\n",
      "Epoch 24/25\n",
      "10000/10000 [==============================] - 1s 110us/step - loss: 5871.8566\n",
      "Epoch 25/25\n",
      "10000/10000 [==============================] - 1s 117us/step - loss: 5842.7111\n",
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: iv\n",
      "Estimand expression:\n",
      "Expectation(Derivative(y, [Z1, Z0])*Derivative([v0], [Z1, Z0])**(-1))\n",
      "Estimand assumption 1, As-if-random: If U→→y then ¬(U →→{Z1,Z0})\n",
      "Estimand assumption 2, Exclusion: If we remove {Z1,Z0}→{v0}, then ¬({Z1,Z0}→y)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+W3+W2+W1+X1+X0+W0 | X1,X0\n",
      "Target units: Data subset defined by a function\n",
      "\n",
      "## Estimate\n",
      "Mean value: 3.3913748264312744\n",
      "Effect estimates: [4.111145  2.1927567 4.0210876 ... 0.9122124 4.053177  2.2825317]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import keras\n",
    "from econml.deepiv import DeepIVEstimator\n",
    "dims_zx = len(model._instruments)+len(model._effect_modifiers)\n",
    "dims_tx = len(model._treatment)+len(model._effect_modifiers)\n",
    "treatment_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_shape=(dims_zx,)), # sum of dims of Z and X \n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(64, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(32, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17)])                \n",
    "response_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_shape=(dims_tx,)), # sum of dims of T and X\n",
    "                                    keras.layers.Dropout(0.17), \n",
    "                                    keras.layers.Dense(64, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(32, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(1)])\n",
    "\n",
    "deepiv_estimate = model.estimate_effect(identified_estimand, \n",
    "                                        method_name=\"iv.econml.deepiv.DeepIV\",\n",
    "                                        target_units = lambda df: df[\"X0\"]>-1, \n",
    "                                        confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{'n_components': 10, # Number of gaussians in the mixture density networks\n",
    "                                                              'm': lambda z, x: treatment_model(keras.layers.concatenate([z, x])), # Treatment model,\n",
    "                                                              \"h\": lambda t, x: response_model(keras.layers.concatenate([t, x])), # Response model\n",
    "                                                              'n_samples': 1, # Number of samples used to estimate the response\n",
    "                                                              'first_stage_options': {'epochs':25},\n",
    "                                                              'second_stage_options': {'epochs':25}\n",
    "                                                             },\n",
    "                                               \"fit_params\":{}})\n",
    "print(deepiv_estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Metalearners"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            X0        X1        X2        X3        X4   Z0        Z1  \\\n",
      "0    -1.511404 -0.796794  0.381156  0.546913  0.988627  1.0  0.945576   \n",
      "1     1.084781 -1.922791 -0.461289 -0.208160 -1.362516  0.0  0.631043   \n",
      "2    -0.993879  0.264821  2.495885  0.768880 -0.200069  1.0  0.564980   \n",
      "3    -2.752826 -1.278350  2.013659  0.225430 -0.296333  1.0  0.729429   \n",
      "4    -0.176540 -0.143080 -1.994730  1.445114 -1.310322  0.0  0.660578   \n",
      "...        ...       ...       ...       ...       ...  ...       ...   \n",
      "9995 -0.615514 -0.054405 -1.403155  0.047062  0.668092  0.0  0.049970   \n",
      "9996 -0.079075 -1.970365  0.590612  0.129645  1.702936  0.0  0.006828   \n",
      "9997 -0.612788 -1.583195 -0.955393  0.139663  0.633461  0.0  0.622620   \n",
      "9998 -2.649140 -1.605903 -1.236836 -1.253810 -1.396432  0.0  0.394677   \n",
      "9999 -1.006696  0.406869 -0.687812  0.247259  0.680919  0.0  0.069656   \n",
      "\n",
      "            W0        W1        W2        W3        W4  v0          y  \n",
      "0     0.083088 -0.950762  0.474662  1.346305 -1.790676   1   1.968666  \n",
      "1    -0.082310 -0.659660  1.051546  1.581745  0.195329   1  10.434861  \n",
      "2     1.221314 -0.894577  0.308504 -0.166613 -1.599631   1  13.105967  \n",
      "3     1.151428 -1.390558  2.628862  0.928396 -0.477143   1   5.665644  \n",
      "4     0.299112 -0.227596  1.117443  0.333261  0.822112   1   9.385864  \n",
      "...        ...       ...       ...       ...       ...  ..        ...  \n",
      "9995  0.688683 -1.733233  1.168689 -0.755424 -0.560150   0  -1.992805  \n",
      "9996 -0.396265  1.375946  0.622630  1.124654  0.143208   1  14.094986  \n",
      "9997  0.726669 -0.896359  2.831088  1.560724 -1.672692   1   8.784501  \n",
      "9998 -0.738313 -1.374355  1.775307  0.559349 -1.305286   1 -10.959509  \n",
      "9999  1.175807 -0.233245  2.094218 -0.437447 -0.219961   1  13.028082  \n",
      "\n",
      "[10000 rows x 14 columns]\n"
     ]
    }
   ],
   "source": [
    "data_experiment = dowhy.datasets.linear_dataset(BETA, num_common_causes=5, num_samples=10000,\n",
    "                                    num_instruments=2, num_effect_modifiers=5,\n",
    "                                    treatment_is_binary=True, outcome_is_binary=False)\n",
    "# convert boolean values to {0,1} numeric\n",
    "data_experiment['df'].v0 = data_experiment['df'].v0.astype(int)\n",
    "print(data_experiment['df'])\n",
    "model_experiment = CausalModel(data=data_experiment[\"df\"], \n",
    "                    treatment=data_experiment[\"treatment_name\"], outcome=data_experiment[\"outcome_name\"], \n",
    "                    graph=data_experiment[\"gml_graph\"])\n",
    "identified_estimand_experiment = model_experiment.identify_effect(proceed_when_unidentifiable=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                                \n",
      "─────(Expectation(y|W3,W2,X4,W1,X3,X1,X2,W4,X0,W0))\n",
      "d[v₀]                                              \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W2,X4,W1,X3,X1,X2,W4,X0,W0,U) = P(y|v0,W3,W2,X4,W1,X3,X1,X2,W4,X0,W0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+X3+X2+X1+X4+X0+W3+W2+W1+W4+W0\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 6.287323303624039\n",
      "Effect estimates: [ 4.61465613 12.53217671  8.69905637 ...  6.04112821 -3.265252\n",
      "  6.54000794]\n",
      "\n",
      "True causal estimate is 5.535602451644429\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "metalearner_estimate = model_experiment.estimate_effect(identified_estimand_experiment, \n",
    "                                method_name=\"backdoor.econml.metalearners.TLearner\",\n",
    "                                confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{\n",
    "                                                    'models': RandomForestRegressor()\n",
    "                                                    },\n",
    "                                               \"fit_params\":{}\n",
    "                                              })\n",
    "print(metalearner_estimate)\n",
    "print(\"True causal estimate is\", data_experiment[\"ate\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Refuting the estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding a random common cause variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Add a Random Common Cause\n",
      "Estimated effect:11.995163235677346\n",
      "New effect:12.028609904102334\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_random=model.refute_estimate(identified_estimand, dml_estimate, method_name=\"random_common_cause\")\n",
    "print(res_random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding an unobserved common cause variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Add an Unobserved Common Cause\n",
      "Estimated effect:11.995163235677346\n",
      "New effect:12.046656228321476\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_unobserved=model.refute_estimate(identified_estimand, dml_estimate, method_name=\"add_unobserved_common_cause\",\n",
    "                                     confounders_effect_on_treatment=\"linear\", confounders_effect_on_outcome=\"linear\",\n",
    "                                    effect_strength_on_treatment=0.01, effect_strength_on_outcome=0.02)\n",
    "print(res_unobserved)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Replacing treatment with a random (placebo) variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Use a Placebo Treatment\n",
      "Estimated effect:11.995163235677346\n",
      "New effect:0.014174758181060527\n",
      "p value:0.40645666707426864\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_placebo=model.refute_estimate(identified_estimand, dml_estimate,\n",
    "        method_name=\"placebo_treatment_refuter\", placebo_type=\"permute\",\n",
    "        num_simulations=10 # at least 100 is good, setting to 10 for speed \n",
    "        ) \n",
    "print(res_placebo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Removing a random subset of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Use a subset of data\n",
      "Estimated effect:11.995163235677346\n",
      "New effect:11.998974744372351\n",
      "p value:0.46530895233272207\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_subset=model.refute_estimate(identified_estimand, dml_estimate,\n",
    "        method_name=\"data_subset_refuter\", subset_fraction=0.8,\n",
    "        num_simulations=10)\n",
    "print(res_subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More refutation methods to come, especially specific to the CATE estimators."
   ]
  }
 ],
 "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.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}