{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mediation analysis with DoWhy: Direct and Indirect Effects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from dowhy import CausalModel\n",
    "import dowhy.datasets\n",
    "\n",
    "# Warnings and logging\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import logging\n",
    "logging.getLogger(\"dowhy\").setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         FD0        W0        v0          y\n",
      "0  11.009406  0.521425  2.654356  33.373490\n",
      "1  -6.848972 -0.637748 -1.947278 -22.097380\n",
      "2  -0.710087 -0.192628 -0.120618  -2.832689\n",
      "3  -0.092251 -0.314083 -0.229222  -1.640303\n",
      "4  13.049299  1.138318  3.028420  41.798294\n"
     ]
    }
   ],
   "source": [
    "# Creating a dataset with a single confounder and a single mediator (num_frontdoor_variables)\n",
    "data = dowhy.datasets.linear_dataset(10, num_common_causes=1, num_samples=10000,\n",
    "                                     num_instruments=0, num_effect_modifiers=0,\n",
    "                                     num_treatments=1,\n",
    "                                     num_frontdoor_variables=1,\n",
    "                                     treatment_is_binary=False,\n",
    "                                    outcome_is_binary=False)\n",
    "df = data['df']\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Modeling the causal mechanism\n",
    "We create a dataset following a causal graph based on the frontdoor criterion. That is, there is no direct effect of the treatment on outcome; all effect is mediated through the frontdoor variable FD0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:dowhy.causal_model:Model to find the causal effect of treatment ['v0'] on outcome ['y']\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAFbCAYAAADfpZU+AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd1iTVx838G8ggGzZiKCAyq4IiIgsq0WrdVG1ilVxr6etrU/RukVc7dvW+mitUmvrrlbROqtWkCUuEJQhokzZIDthJDnvH77JKwKyktwJnM915RIzzvklufPLybnPYBFCCCiKoqju5k8FpiOgKIqiJIMmeIqiqG6KJniKoqhuis10AJTkVFRUgBACLpeLuro6EEJQUVEhur2+vh4cDqfVx3M4HNTX17d6u4KCArS1tVu9XUlJCRoaGqL/q6urQ1lZGYqKitDS0gIAaGpqgs2mhyElm/h8PoqKilBUVISKigrw+XxUV1eDx+NBTU0NKioqUFVVRe/evWFiYgIdHR2mQ26CfrIYVllZifLyctGlpqYGXC4XlZWVor+rq6tRXV0NDoeD2tpaVFZWgsPhgMvlory8XFSOQCBoM2nLqje/DLS0tKCoqAgNDQ2oqqpCU1MTWlpaUFVVhbq6OrS1taGmpgZVVVXo6OhAVVUVampq0NbWhpaWFnR0dEQXRUVFhp8ZJQ+4XC7u37+PJ0+eICkpCUlJScjIyEBRUREEAkG7y1FVVUXfvn1hbW0NBwcHODg4wNnZGba2tmCxWBJ8Bi1j0VE04lFTU4OioiIUFxejtLQUJSUlKC4uRkVFRZME/vb/W3r5hS3cNxOcpqYmVFVVoaGh0SzZKSgoQENDA0pKSi0mSmVlZairqwMAevfuLTrQWCwWevfu3epzYrPZ0NTUbPX2jv4CELZ8GhoaUFtbC6DlLybhl1ZVVRW4XK7oS43L5YLD4aCiogJcLlf0BSf8hdISTU1NUbLv3bt3k+Tfu3dv6OrqwsjICIaGhjAwMICBgQH09fVbfU5U9yAQCHD37l1cvXoVERERuH//PhoaGqCrq4v33nsPdnZ2sLa2Rp8+fWBiYgIjIyPo6upCQUFB9KtTeHzX1dXh1atXyM/PR35+PnJzc5GSkoKUlBSkpqaioaEBhoaG8PLywvvvv49JkybBzMxMGk/zT5rg36GhoQF5eXnIy8tDbm4uCgoKUFhYKErixcXFKCoqQklJCbhcbpPHamhowMDAALq6us2Syrv+r6mpCWVlZYaesfwSdj9VVVW98wv17f+/evUKpaWlTb5o2Wy2KNG/nfxNTU1hYmKCvn37ol+/fk26oCjZFx0djRMnTuDvv/9GQUEBBgwYgPfffx8+Pj7w9vZGv379xFofj8fDo0ePEBkZiYiICERERKC6uhpDhw7F1KlTMXfuXPTp00esdb6h5yZ4Ho+H3NxcZGVlITc3Fy9fvhR9+wr/LiwsFN2fzWbD2NgYRkZGMDIygr6+PgwMDGBsbCz68BsaGsLIyAgGBgbo1asXg8+O6gg+n4+SkhLRpbCwECUlJSgtLW32hZ6Xl9fky1xLSwumpqYwNTVF3759YWZmhr59+6Jv374wNzeHpaUlVFVVGXx2VHV1NY4ePYoDBw4gKSkJjo6OmDp1Kvz8/ODg4CDVWBoaGnDr1i2cP38e58+fR2VlJSZNmoRly5bhgw8+EHd13TvBl5eXIz8/HwUFBcjIyGhySUlJEX1QlZWVoaenBxMTE1haWop+lr35d79+/ejJQArA6/5a4TH15vEl/Fv4r5COjg4sLS1bvPTv35+eJ5CQmpoa/Pbbb9i5cycqKysxceJELFmyRBKJtFMaGhrw999/IyQkBLdu3cLgwYOxfv16TJs2TVz99fKf4AUCATIzM0X9XampqUhOTsbz589FfbmKioowNTVt8QNmbm4OQ0NDhp8F1d3U1tYiMzOzWcMiIyMDmZmZonMGvXr1gqWlJWxtbWFraws7OzvY2trCxsaG/grsJD6fjwMHDmDTpk3g8/n44osv8NVXX8ncCJc3PXz4EFu3bsXly5fh5uaGffv2wcXFpavFyleCz8jIwKNHj0RJ/OnTp3j69Know2JmZib6kFhZWTVpJdF+bUpWEEKQn58vSvjp6el4+vQpUlJS8Pz5czQ2NkJBQQEWFhaihG9ra4vBgwfDwcGBHsvv8PDhQyxbtgyPHz/GV199hW+++UamE/vb4uPj8dVXXyEmJgZLly7Fjh073jkUuQ2ym+Dz8/MRFxcnuty7dw8lJSUAgD59+sDe3h52dnaifwcPHiwaW01R8orH4yEnJwfJyclISUlBRkYGkpOTkZCQgNraWrDZbFhZWcHFxUV0cXZ2hpqaGtOhM0ogEOC7777Dpk2bMGLECOzfvx92dnZMh9UphBAcP34cgYGBUFVVxalTpzB8+PDOFCUbCb6oqAhRUVGIjo7Gw4cPkZiYiJqaGigpKcHBwQFOTk5wdnaGs7MzHB0de/zBTPU8fD4faWlpiI+PF10SEhJQWVkJNpsNW1tbODs7Y8SIEfDy8oKNjQ0j466Z8OrVK8yYMQORkZHYsWMHVq1a1S2ee0lJCQICAvDvv/9i586d+O9//9vRIphJ8BkZGYiKihIl9bS0NCgqKsLR0RHDhg0TJfP33nuP/hylqFYQQvDixQtRwn/w4AHu3buH2tpaGBgYwNPTE15eXvDy8sKQIUO65SCB7OxsjBs3DhwOB2fPnsXQoUOZDkmsCCH4/vvvsXbtWixfvhx79uyBgkK7V5iRToIvLy/H9evXcfXqVYSFhSEvLw+9evWCq6srvL294enpiREjRtAuForqIh6Ph/j4eERHRyMyMhIxMTEoLS2FhoYGPD09MW7cOHz00UcYMGAA06F22bNnz/D+++/DwMAAV69ehYmJCdMhSUxoaCg+/fRTTJo0CSdPnmzvyCvJJfikpCRcuXIFV69exZ07dwAAnp6e8PX1hbe3N1xdXaGioiKJqimK+n8IIUhNTUVUVBTCwsJw48YNVFRUwNraGhMmTMD48ePh5eUFJSUlpkPtkPz8fHh4eMDY2Bj//PNPV05Eyo2IiAiMGzcOs2fPxsGDB9vTDSXeBJ+UlIQjR47gr7/+QnZ2NgwNDTFu3DiMHz8eY8aMeee0eIqiJI/H4yE6OhpXr17FlStXkJKSAi0tLXz00UeYO3cufH19ZX5cfm1tLdzd3cHn8xEZGQk9PT2mQ5Kay5cvw8/PD5s2bcLGjRvbuvufIF1UUlJC9uzZQ1xcXAgAYmFhQdavX0/u3r1L+Hx+V4unKEqCMjIyyN69e4mnpydhsVjExMSEBAYGkqSkJKZDa9XixYuJnp4eyc7OZjoURvz8889EUVGRREREtHXXU51O8OHh4cTPz48oKysTTU1NMn/+fBIREUEEAkFni6QoikHp6elk06ZNxMLCggAgLi4u5NChQ6Suro7p0ETOnz9PWCwWOX/+PNOhMMrPz4+YmZmRioqKd92tYwleIBCQc+fOEWdnZwKAeHl5kWPHjpHa2tquRUtRlMwQCATk9u3bZO7cuURFRYUYGxuTnTt3kpqaGkbjqqurIxYWFmTu3LmMxiELSktLiZ6eHvnmm2/edbf2J/jbt28TV1dXoqCgQKZNm0bu37/f9Si7gTa+QeVOeXk50yHIne78muXn55O1a9cSTU1NYmRkRPbt20d4PB4jsXz//fdETU2NvHz5kpH6Zc3u3buJqqoqycnJae0ubSf4V69ekXnz5hEWi0U+/PBD8ujRI/FG2QXXr18n/v7+BAABQObMmUOSk5NFt0dERJBJkyaJfm2I62ddY2Mj2blzJ/Hw8CCKiopiKZNJXC6XbNu2jQwfPpwoKCi0eX+BQEBOnz5NPvroIzJkyBDi6+tLJk6cSFasWEF27txJVq1aJZE4BQIB2bNnD1m9ejUZOXIk8fT0JE+fPpVIXW3p6Gsm70pKSsiqVauIsrIyGTp0qNTzAI/HI6ampiQwMFCs5a5bt46oqKgQAOSDDz4gkZGRJDc3lyxfvlyUV/z8/EhYWJjoMWFhYWTYsGGExWKRzz77jDQ0NBBCXh+fhw4dItOmTSPr1q0jCxcuJCdOnBBrvG+qr68npqam72rFvzvBx8fHEwsLC2JiYkLOnTsn/gjFgMvlEgBEW1u7xZO6+fn5BADJy8sTa70cDofo6OgQoMvnqWVCe59PcXExGTlyJBkwYAC5e/eu6JwLn88nx44dI7q6umTBggUSifGnn34i6urqpLGxkZSXlxM/Pz9y7949idTVHt3tGGiP5ORk4unpSXr16kVCQkKkVu/Vq1cJi8Ui6enpYi977dq1BAD573//2+T6KVOmEADk6NGjzR5z8OBBEhAQ0OS6oKAg0r9/f/Lq1StCyOvGcf/+/clPP/0k9piFNmzYQIyNjUljY2NLN7ee4CMjI4mGhgbx8fEhhYWFEgtQHAAQa2vrFm/j8XgEgER+VlpbW3erD3dbz4fP5xN3d3eio6NDSktLW7xPeHg4mTFjhsTis7KykkjZndXdjoH2EAgEZNeuXURBQYGsXr1aKnXOnTuXeHp6SqTssrIyoqqqSvr27dskT8TFxREAZPTo0c0eExAQQGJjY0X/z8rKImw2m+zYsaPJ/bZt20bU1NRa/bx01YsXLwgA8u+//7Z086kW57wmJydj3LhxmDBhAv79918YGRl1YeQms4RjemV9bK88CA0NRWxsLL755ptWxx6PHDkS06dPl0j9ubm53WKNEXnHYrGwZs0aHDp0CN9//z327t0r8TqjoqLg6+srkbJ1dXXh5+eHvLw8XL9+XXT9kCFDoKOjg1u3buHFixei62tra/H48WO4ubmJrjtx4gR4PB5Gjx7dpOxRo0aBw+Hg0KFDEolduGJudHR0i7c3S/CEEMyePRtDhgzBsWPHutX6FYQQ/P3331iyZAlMTU1RXl6OgIAA6OnpwcHBAQ8fPhTdt7KyEoGBgfjmm2+watUqjBkzBqtWrRKtMf+m9PR0TJw4ETo6OnB1dUV4eLjotgcPHsDNzQ2LFy9GYGAgFBUVUV1dDeD1xhHffvstFi5ciKFDh+KDDz7AkydPwOfzcfv2bXz55ZcwNzdHXl4efHx80K9fP4SEhEBXVxcsFgsbNmwQ1bN//34oKCggJCTknWULcTgcrFq1CkuWLMGGDRuwdu1a0T6prQkNDQWAZgfx26ZOndqu17G978fly5exbNkycDgcFBYWYtmyZVi2bBlqamrafJ+EM/6EXwxVVVX44YcfRNd15Jhoz2vWmfe0vLz8nceJrJo/fz6CgoKwevVqPHv2TGL1FBUVITMzE+7u7hKrIyAgAACaJOLw8HDRXsZvXn/27FlMmjSpSWNDmGBNTU2blCvcezUxMVEygQMYMWIEYmNjW77x7Tb9tWvXCIvFanKyUtbhHV00wtsJef3TMjc3l6irqxMAZNu2bSQrK4scO3aMACDDhg0jhBBSVVVFBg0aRDZv3iwqo6ioiAwaNIhYWFiIRk0If56vXLmS3Lhxgxw4cICoqakRBQUFkpiYSAghZNCgQURHR0d0fmDatGmkqKiIEELIokWLSGpqqqgOX19fYmhoSIqLi0lMTAxRVVUlAMiOHTvIzZs3ycKFC0l1dTX53//+RwCQq1evih6bnZ1N/P39Rf9vrezKykrS2NhIhg0bRhYtWiTqQ3/+/DlRVFR8Z3fD0KFDCYB2jxxq63V89epVu94Pobff5/a+T5aWls2el/C69h4T7X3NOvuevus4kWU8Ho9YW1uTFStWSKyOe/fuEQASndjE4/GIiYkJYbPZpKCggBBCiL+/PwkPDyfq6urEyMhIdDLV29u72bkAR0dHAoBwOJwm19fW1hIAZPjw4RKLfevWra3lv+Z98IGBgcTZ2VliwUhCexO8kJWVVZPrBAIBMTQ0JMrKyoSQ12fWAZD8/Pwmjzty5AgBIDqTL0zwlZWVovv89NNPBIBorK6+vj4BQH755RciEAhIYmIiqaioIHfv3hWdpX/7cunSpSZxlpWVNYmjvr6emJmZkYkTJ4qu27BhA4mPjyeEkDbL3rt3LwHQ7Et80KBB70zwbm5uLb4urWnv69jW+yH09vvc0ffpTW9f11YM7XnNuvKetnacyIPNmzcTGxsbiZV/5coVAkDi4/DXrFlDAJBdu3aRsrIy4uzsTAQCAVmwYAEBQM6dO0fS09OJu7t7s8d6eXkRAITL5Ta5nsPhEAASzak///wz0dfXb+mm5n3wZWVlcreFHZvNhkAgaPE2Ho/XbCGlt/txWSwWdHR00NDQAACIiYkBAGhqaja5n7e3NwCIFk8TenMVzClTpgAAUlJSAAC//PILNDQ0sHz5cnh4eKC+vh7a2tp48OAB7OzsQAhpdpkwYUKTOHV1dZvUp6ysjJUrV+Ly5ct48eIFGhoakJaWBicnJwBos+wbN24AACwsLJqU29YypMINFFJTU995P6H2vo5tvR9dLb892oqhPa9ZV97T1o4TeWBkZCTajEcShHsnS3oLwze7aY4fP46ZM2eCxWJh0aJFAICQkBD88ccf+PTTT5s91sbGBgBQUVHR5HphV6EkV7rU0NBotXu12Sfa0tISycnJrSZMWWRubo7KysoWb3v16lWHFyMSfmizsrKaXC882fyuD57wPv369QMATJs2DQkJCRgzZgxiY2MxYsQIHD16FGVlZcjMzGzxjeHz+W3GuGjRIqirq2Pfvn24cOECpk2bJrqtrbLz8vJE9+sIHx8fAMDdu3fbdf+uvI6yUP6b2vOadeU9be04kQeJiYkYOHCgxMoXbrn3dvIUN1tbW7i6uuL58+cIDg4WJfLhw4fDzs4ON27cwOHDh/HJJ580e6y9vT2A16tcvkm4+bqnp6fE4i4rK2vWYBBqluBnzJiBvLw8nDlzRmIBiZuLiwuKi4uRmZnZ7LaIiAh4eXl1qDxhC/DKlStNrs/NzQWAd+7KLryPsMW2adMmDBgwANevX8fJkyfB4/Gwbt062NjYiE7IvSklJQX79u1rM0ZtbW0sWrQIhw8fxunTp+Hn5ye6ra2yha2Nt59fW2bPng1nZ2fs2bOn2YEsVFdXhz/++ANA117H9mhv+cJWc319PYDX27sJGwSknYuptuc168p72tpxIusKCgpw8uRJ+Pv7S6wOYQNNkr8ShISteFdXV1Grm8ViYeHChSCEYOjQoTAwMGj2uDlz5kBbW7vJAAsACAsLg5KSEmbNmiWxmEtKSlpvxLbUcbNkyRKip6dHnj9/Lta+IklJS0sjKioqxMXFRTRtt76+nly6dIkYGRmJ+qaF+vfvLzrBJmRiYkIAkIaGBlJbW0vs7e1J3759m/TvfvHFF2TEiBGiky02NjZN+lMFAgFZvnw5mTRpkqhsVVVV0cSHhoYGoqWlRVxdXQmXyxUt6rRgwQJy/Phxsn79euLr6yvq0xfGWV1d3eLzzsjIIAoKCiQ4OLjJ9W2V/ejRI6KoqEh0dXXJtWvXCIfDIbdu3SKampoEAMnIyGj1tU5JSSH9+vUjFhYW5Ny5c6IJFrW1teTWrVtk1KhRovHB7X0d23o/CHk9Vhn/b7VSofaWL5ywsmHDBvLs2TPy448/iiYoXbt2jfB4vDZjaM9r1pX3tLXjRJbV19eTUaNGkUGDBkm0f5zD4RBlZWVy8uRJidUhVFpaSpSUlMiff/7Z5Pri4mKipKRETp8+3epjd+3aRQYNGkSqqqoIIYRUVlaSgQMHkqCgIInGPHbsWPLpp5+2dFPLE51qamrI0KFDiZmZGWNTwTvq6dOnZOrUqcTCwoKYm5uT/v37k+nTp5PHjx83ud++fftEJ76Cg4NJRUUF2b17t+i6NWvWEA6HQ6qqqkhgYCDx9fUlq1atIoGBgWTr1q1NVta7ceMGmTBhAvHx8SGLFi0in3/+ebO1OgAQJycnsnPnTjJr1izy0UcfiRJoZmYmmThxItHR0SFGRkZk8eLFpLi4mNTU1JCgoCBRTIsXL272JSW0cuXKFidRtFa2UEREBBkxYgTR0NAglpaWZOfOncTLy4ssXbqU/Pvvv++cGFZVVUV27dpFxo8fT8zNzYm9vT1xdHQk69ataxZLW69je96Pe/fukaVLlxIAhMVikS1btpCEhIR2lU/I6wbAsGHDiJqaGvH19SVpaWnE09OTzJ49m5w6dYr88MMP7Tom2vOadfY9fddxIotqa2vJxIkTiba2dqvHpjgNGzaMfPbZZxKvhxBCFi5c2Gw0zLuuFxIuVTB79myybt06Mm3aNBISEiLRFXb5fD7p3bs32bdvX0s3n2p1w4/y8nKMHz8eKSkp+O2335r08VIU1XOlpaVh+vTpyM/Px+XLlzF8+HCJ17lu3TqcPHkSGRkZHdmTtNsLCwvD6NGjkZKSAltb27dv/rPVV0pHRwe3b9/Gp59+iunTp2PatGmiE00URfU89fX1CA4OxpAhQ6Cqqoq4uDipJHcAmDdvHnJycnDr1i2p1CcvfvvtNwwfPryl5A6ghZOsb1JRUcH+/ftx+/ZtJCcnw8rKCitXrkRxcbFEgqUoSvYIBAL89ddfsLe3x86dO7FmzRpERUWhf//+UovBysoKHh4e2L17t9TqlHU5OTkIDQ0VDeNsUXv7erhcLvnxxx+JgYEB0dTUJKtWreqxW2ZRVE9QW1tLfv75ZzJgwACipKRElixZQnJzcxmLJzw8nAAgN2/eZCwGWTJnzhxiYWHxrh23Or5lX1VVFfn++++JmZkZYbPZZMKECeSvv/6SqW29KIrqvNjYWLJ8+XKio6NDVFVVyfLly2VmRN348ePJe++91+PzTUxMDFFQUGg22uctrZ9kbUtjYyNCQ0Pxxx9/4ObNm9DS0sKMGTMwd+5ciS4KRFGU+OXk5ODYsWM4duwY0tLSYG9vj4CAAMybN6/Fcd9MycjIgJOTExYsWNBju2sqKyvh5OQEGxsbXLly5V0rrP7Z6QT/pvz8fJw8eRJHjhxBUlISBgwYgIkTJ2L8+PHw8fGBsrJyV6ugKErMEhMTcfXqVVy5cgWxsbHQ09ODv78/AgIC4OzszHR4rTpx4gTmzJmDM2fO9LjRfXw+H1OnTsX9+/eRmJjY1pfvn2LfqSAuLo6sWbOGvPfeewQA0dDQIH5+fuTQoUPtXqSKoijxq62tJRcvXiRLly4lZmZmBAAxMjIiCxYsIBcvXhRNDJMHn3/+OVFRUWmylV5PsGTJEqKqqkqio6Pbc/fOd9G0R3Z2tqiFEBYWhrq6Ojg4OMDb2xuenp7w8vJC3759JVU9RfVoNTU1iI2NRVRUFKKionD37l00NDTA2dkZH330ESZMmABnZ2e5HFcuEAjw6aef4urVq7h48aJonaTuSiAQ4L///S/27duHc+fOYdKkSe15mHi6aNqDy+UiPDwct27dQnR0NOLj48Hj8WBpaQkvLy94eXnB09MT1tbW0giHorqd4uJixMTEIDIyEtHR0UhISACPx8PAgQPh6ekJHx8fjBs3Tq53aHtTQ0MD5syZg7///htHjx5tcRGw7qC+vh7z5s1DaGgojhw5gpkzZ7b3odJL8G/jcDiIj49HTEwMoqOjERUVhcrKSmhpaeG9996Di4uL6GJrayuXrQyKkpTy8nIkJycjLi5OdElNTQWLxYKNjQ08PT3h4eEBHx8fqY5XlzaBQIBVq1Zh7969WLNmDbZu3dqtdqHLysqCv78/UlNTERoailGjRnXk4cwl+Lc1NjaKDtT4+HjEx8cjOTkZjY2N0NTUxJAhQ+Ds7AwnJyfY29vDxsYGGhoaTIdNURLF4/Hw4sULJCcnIzExUfTZEK7maWFhAWdnZ9HF3d1dbtaRF6dDhw5h5cqVcHR0xPHjx2Fpacl0SF125swZLF26FGZmZjh9+nSrs1XfQXYSfEvq6+vx5MkTxMfHixL/kydPUF9fDxaLhX79+sHGxkaU8O3s7GBra9vq2sgUJavq6urw9OlTPH36FCkpKUhNTUVqairS09PR0NAAFouFgQMHNknmLi4uorXSqdfLMs+cORPPnz/HunXrEBgYCBUVFabD6rAXL17giy++wLVr17B06VL8+OOPUFVV7UxRsp3gW8Ln85GRkdHkQ5CSkoKnT5+ipqYGAGBoaAhra2tYWlpiwIABop3HLS0tu03/IyV/qqurkZGR0eySnp6OrKws8Pl8KCkpwdLSslmjxcbGBmpqakw/BZnX2NiIH3/8EcHBwTA2NsbWrVsxY8YMKCoqMh1am4qLi/H9999j7969GDhwIPbt29fVk8fyl+DfJTs7W9QCSk9PF32AsrOzRVuvqaurN0n4FhYWMDMzg4mJCUxNTWFsbEz7+6lOefXqFfLz85GTk4P8/HxkZWU1SeTCDStYLBZMTExEx+DAgQNFyXzgwIF03kgXPX/+HFOmTAGHw0Fubi4GDhyItWvX4pNPPpH4tn+dkZOTg71794q2bVy7di1WrFjRbKvRTuheCb41fD4fL1++FH3QMjMzRX9nZWWhqKhIdF82mw1jY2P069cPJiYm6Nu3r+gLwMzMDAYGBjAyMkLv3r0ZfEaUNHE4HJSUlKCwsBCFhYWiBJ6Xl4fc3Fzk5+cjNzdXtHco8HqfWHNzc1Ej4u1GhSwmmu7g0qVLCAgIgLm5Oc6ePQs+n49t27bh1KlT0NLSQkBAABYuXCjaX5gpjY2N+OeffxASEoJr167B0NAQX3/9NZYtWybOX2o9I8G3pb6+vsUPbF5eHvLz85GdnY2ioiLweDzRY5SVlWFgYABDQ0MYGxvDwMAA+vr66NOnDwwMDEQXHR0d0YX+MpANVVVVKC8vR0VFBcrKylBYWIiSkhJREhf+XVRUhKKiomZ7rBoYGIi+8E1NTWFiYtKsQfD2RuCUZPH5fAQHB4v2Uj148GCTfuvCwkIcPnwYv/76K7KysmBra4upU6diypQpGDJkiFS6cKqrqxEWFobQ0FBcunQJFRUVGD16NJYuXYrJkyeLo8X+Nprg24vP56OoqKjVJFBaWori4mLRbW+25oS0tLSgo6OD3r17i5K+8G/hv+rq6tDQ0IC2tjZUVVWhpqaG3r17N/m7p6qtrQWHw0F1dTWqq6vB5XJRU1ODqqoqcDgc1DsA05AAACAASURBVNbWoqKiAuXl5aKL8P9vXv/2hvJsNlv0hWxkZARDQ8MmX976+vowMDCAsbExjI2NaetbxpSUlGDWrFmIjo7G3r1737l8rkAgQHR0NEJDQ3H+/Hnk5ORAW1tbNPHS2dkZ9vb2ov1YO4vH4+HZs2dISkrC3bt3ERUVhYSEBAgEAowYMQIff/wxPv74Y0kPYaUJXlJqampQUlLSapJ5+zrhv7W1taKTxa1RU1ODqqoqtLW1oa6uDmVlZSgrK0NdXR3A6w25FRQU0KtXL1ErRjjaQlVVtUmCYrFY7/zS0NDQaLVlUVFR0eqm1bW1taLzHsDrL8iqqioAr1vQfD4fPB4P1dXVTa5raGhAbW0tCCGoqKgAh8MBl8sVbZLdGhUVFairq7f5BdrSbbK0mBbVMZGRkZg5cyY0NDRw9uxZDB48uN2PJYQgKSkJEREROHToEDIzM0XHqK6uLqysrNCnTx+YmprCyMgIWlpaUFFRgZqaGlRUVFBdXS06hisrK/Hy5UsUFRUhOztbNPqJzWbDzs4OPj4+8Pb2hre3NwwNDSX1cryNJnhZJWyVcjicJomuoqICtbW14HK5qKqqEh1kdXV1ol8N5eXlAF73HdfX10MgEIgS5NuJV5hQWyMsqyXCL5eWvPmFA/z/LxJCCJ4/fw5jY2OYm5uLxmwLy2Kz2aLujbd/uaipqUFNTQ1aWlrQ0NCAmpqa6NcO7f7qWQgh+N///ofAwECMGzcOR44c6fSv28rKSlhYWOCrr77CihUr8PjxY6SkpOD58+coLCzEy5cvUVxcLOqmbWxsRENDg6jxo6mpCS0tLfTt2xfGxsYwMzMTDd+2tbVlcqim+Bcbo7q/b7/9lvTr16/Tj9+1axdhs9nk9u3bYoyK6ikqKyvJtGnTCJvNJrt27eryptYbN24kvXv3JuXl5e+83+jRo8miRYu6VJeUneo+c3opqbG0tMTLly9RX1/fqdbJ6tWrERcXhxkzZiAuLo4uOEe1W0JCAqZNm4bq6mpcv369o1P3m6moqMDevXuxevXqNn8BZGdnY+TIkV2qT9ro71qqwywsLCAQCJCbm9upx7NYLBw+fBh6enqYNm1aky4jimrN0aNH4eHhAVNTUyQkJHQ5uQPAd999B0VFRXz22WfvvJ/weJe3dX1ogqc6TLjOR0ZGRqfL0NDQQGhoKFJTU7Fq1SpxhUZ1Q3V1dVi6dCnmzZuHRYsW4ebNm+jTp0+Xyy0rK8O+ffuwevXqNoe1FhQUoL6+Hubm5l2uV5poFw3VYcIRKF1J8ABgbW2NI0eOwM/PDy4uLpg/f76YIqS6i+zsbEyfPh1paWk4e/YsPv74Y7GV/d1330FFRQUrVqxo875ZWVkAIHcJnrbgqU6xtLREZmZml8uZPHkyAgMDsWLFCsTFxYkhMqq7uHTpEpycnMDj8fDo0SOxJvfS0lLs378fa9eubdeqtNnZ2VBSUury+Hhpowme6hQLCwuxJHgA2LlzJ3x8fDB16lSUlpaKpUxKfvF4PGzZsgVTpkzBhAkTEBMTI/blf3fu3AkNDQ0sW7asXffPysqCmZmZXCxa9iaa4KlOsbS07HIXjZCCggJOnjwJBQUF+Pv7g8/ni6VcSv6UlJTgww8/xLfffouDBw/i6NGjnV0qt1WFhYU4cOAAvvnmm3av+5KdnS133TMATfBUJ1lYWIgtwQOvZw6eO3cOMTExCAoKElu5lPyIjIyEo6MjcnJycO/evXcuOdAVu3btgra2NpYsWdLux2RlZcndCBqAJniqkywsLERLLIiLk5MTDhw4gG3btiE0NFRs5VKyjRCCPXv24IMPPoCrqyvu37/foSUHOqKgoAAhISFYv359h34ZZGdn0wRP9RziGCrZkrlz52LhwoWYN28eUlNTxVo2JXuqqqrwySef4Ouvv0ZwcDAuXLgg0QX1tm/fDn19/Q79OiCEICcnh3bRUD1H//79oaCgILYTrW/at28fbGxs8PHHH4sWI6O6n4SEBDg7OyM2Nha3b9/GmjVrwGKxJFZfbm4uDh06hPXr13doBnZRURG4XC5N8FTPoaKigr59+4q9BS8s+9y5cygrK+tQPyklP96clfrw4UN4eHhIvM7t27fDyMiow/MthGPgaRcN1aOIc6jk28zMzPDnn3/i7Nmz+OmnnyRSByV9dXV1WLJkiWhW6r///gtjY2OJ15uTk4Pff/8dGzdu7PCWiNnZ2WCz2TA1NZVQdJJDZ7JSnSauyU6tGTVqFLZt24bAwEA4OTl1dQNiimGSnJXaluDgYJiYmGDu3LkdfmxWVhb69u0LNlv+0iVtwVOdJu6hki1ZvXo1/Pz8MGPGDLx8+VKidVGSI8lZqW3JysrC0aNHsXnz5k5taC6vY+ABmuCpLrC0tERWVpZEJya9ufLk9OnTUV9fL7G6KPGTxqzUtgQFBaFfv36YPXt2px4vr2PgAZrgqS4YNGgQGhoakJ2dLdF6NDQ0cP78ebrypJwpLi6W+KzUtjx//hzHjx/H5s2bO93Fkp6ejoEDB4o5MumgCZ7qNGtrawBAWlqaxOuysrLCkSNH8Msvv+D333+XeH1U10RGRmLIkCESn5Xalq1bt8LCwgIzZ87s1ON5PB6ys7MxaNAgMUcmHTTBU53Wu3dvGBkZSSXBA3TlSXkgzVmpbUlPT8epU6cQFBTU6dZ7VlYWGhsbaYKneiZra2upJXiArjwpy6Q9K7UtmzdvhpWVFWbMmNHpMtLT0wGAdtFQPZO0EzxdeVI2SXtWaltSUlJw+vRpbNmyBQoKnU9z6enpMDAwgLa2thijkx6a4KkukXaCB+jKk7KGiVmpbdmyZQtsbW0xderULpXz/Plzue2eAWiCp7rI2toa+fn5qKqqkmq9Tk5OOHjwIF15kkFMzUptS3JyMs6dO4etW7d2qfUOvG7B0wRP9VjCkTTCvkppmjNnDl15kiHZ2dnw9vbG6dOncfbsWezZs0dmZnpu2rQJ9vb2mDJlSpfLkuchkgBN8FQXWVhYQFlZWerdNEJ05UnpY3JWaluePHmCCxcuYNu2bV1uvfN4POTk5NAWPNVzsdlsWFpaMpbg6cqT0iMLs1LbsmHDBjg5OWHixIldLiszMxONjY1y3YKXjd9UlFxj4kTrm4QrT44dOxZubm748ssvGYuluyouLsasWbMQExODgwcPMjZx6V3i4+Nx6dIlXL58WSwjeOR9iCRAW/CUGDCd4IGmK09GREQwGkt3IyuzUtuyadMmuLi4YNy4cWIpLz09HYaGhnI7RBKgCZ4SA2trazx79gwCgYDROOjKk+L15qzUYcOGMTortS0PHz7E1atXsW3bNrGNv5f3IZIATfCUGFhbW4PD4SAvL4/ROOjKk+JTVVWF6dOni2alnj9/ntFZqW3ZsGED3N3dMXbsWLGVKe8jaACa4CkxkOaiY22hK092nXBW6t27d2ViVmpb7ty5g+vXr2Pr1q1iLZcmeIoCoK+vDz09PZlI8ABdebIrZHFWals2bdoEDw8PjB49WmxlNjY2yv0QSYAmeEpMrKysZCbBA3TlyY6S1VmpbYmJicGtW7ewbds2sZabmZkJHo8n9wmeDpOkxEIWRtK8befOnUhMTMTUqVPx8OFD6OvrMx2STGJyr9Su2rBhA0aNGoWRI0eKtVzhEMkBAwaItVxpoy14SixkMcG3tfJkTk4O/vjjD2aCkxGyPCu1LWFhYbh9+zY2b94s9rK7wxBJgCZ4Skysra2Rm5sLDofDdChN6OrqIjQ0FDExMdiyZYvo+lu3bsHR0RHz58+XuS8maZCHWaltCQoKwpgxY+Dt7S32srvDEEmAdtFQYmJjYwOBQIBnz55hyJAhTIfTxJAhQ3Dw4EEEBARgyJAheP78OdatWwcWiwUlJSWcPn0amzZtYjpMqZGHWaltuXnzJiIjIxEVFSWR8tPS0mBlZSWRsqWJRQghTAdByT8ejwcNDQ389ttv+PTTT5kOp0Xz5s3D7du3kZub22RSlqWlJV68eMFgZNITGRmJmTNnQkNDA2fPnpXZiUtt8fb2hoaGBq5evSqR8k1NTbFy5UoEBgZKpHwp+ZN20VBiwWazYWVlheTkZKZDadGzZ89w584d5OXlNZtxm5GRgSdPnjAUmXTI06zUtly7dg1RUVES6XsHXk/yys/Ph62trUTKlyaa4Cmxsbe3R1JSEtNhNHPx4kU4OzuLhr69TVlZGadPn2YgMvGJioqCmZkZUlJSmt0mb7NS2xIUFISJEyfCzc1NIuU/ffoUhBCa4CnqTfb29jLXgt+7dy+mTJkCDofTYnIHgIaGBhw7dkzKkYkPl8vF3Llz8fLlS0yePBk1NTWi2+RtVmpbLl26hPv372Pjxo0SqyM1NRW9evWCubm5xOqQFprgKbFxcHBAVlYWamtrmQ5FpF+/flBXV29zt6GcnBw8fPhQSlGJ18aNG0WLq2VnZ2P+/PkA5HNWaluCg4MxefJkuLq6SqyO1NRUWFtbQ1FRUWJ1SAtN8JTY2NvbQyAQyNT2eZMnT0Z6ejqmT58OAK3u8iOv3TT379/H7t27Rb9OGhsbce7cOYwZMwbz5s3Dl19+iVu3bsnFrNS2XLhwAQ8fPpRY37tQampqt+ieAWiCp8TI0tISqqqqMtcPb2xsjBMnTuDSpUswNDRssTXf0NCA48ePQ54GlTU0NCAgIKDZlxYhBLdu3cK3336L7du3d4uWKCEEQUFBmDp1qsSH4dIET1EtUFRUhI2Njcz1wwtNmDABaWlpWL58OVgsVrNEX1hYiNjYWIai67jt27cjPT29xXMLLBYLP/zwA0pKShiITPzOnTuHx48fS7TvHXj9pZmZmdltEjwdB0+J1Zw5c1BWViax8cnicufOHcyfPx8vXrwQLWGgrKyMZcuWYc+ePe0qg8vlIi8vD69evUJVVRUEAgEqKysBANra2lBQUIC2tjb09PTQp08fqKqqii3+x48fw8XFpdUTxwCgpKQEd3d3hIWFyXUrXiAQwMnJCXZ2djh16pRE63ry5AkGDx6MJ0+ewMHBQaJ1ScGfNMFTYrVr1y7s378fOTk5TIfSpvr6euzYsQM7duwAi8VCY2Mj9PX1UVhY2CQhFhYW4t69e0hOTsaTJ0+QkpKCnJwcVFRUdKg+XV1dmJmZwc7ODu+99x4cHBzg5uYGQ0PDDpXD4/Hg6uqK5ORkNDY2tnn/bdu2Yf369R2qQ5acPn0as2bNQmJiosST7pkzZzBr1izU1tZCRUVFonVJwZ90qQJKrBwcHPDy5UtUVlbK/EJNKioqCAoKwieffIIFCxbg/v37KC0txfXr11FbW4t///0XkZGRePr0KVgsFszNzeHg4IDx48fD3NwcJiYmMDExgYGBATQ1NcFisUTjyysqKkAIQVVVFUpKSlBQUIC8vDxkZWUhKSkJISEhyM7OBiEEdnZ28Pb2hq+vLz788EOoqam9M+4ff/wRT548abZ42pvYbDb4fD6UlZXb9SUgqwQCAbZv3w5/f3+ptKhTU1MxYMCA7pDcXyMUJUY5OTkEAImMjGQ6lA6pra0lc+fOJWw2m7DZbKKoqEhGjBhB1q5dS65du0aqqqrEXmdlZSW5cuUKWbNmDRk+fDhRVFQkampqxM/Pj5w6dYrU1dU1e0xaWhpRVlYmAJpdlJSUCIvFIsrKymT8+PHkyJEjEolbmo4fP04UFRVJamqqVOqbNm0a+fjjj6VSlxScogmeEjs9PT2yd+9epsNol2fPnpFVq1YRXV1doqSkRHx9fckvv/xCSkpKpB5LcXEx+fXXX8m4ceOIkpIS0dfXJ4GBgeTFixeEEEL4fD5xd3cnSkpK70zq1dXVUo9dEng8HrGxsSEBAQFSq9PGxoZs2rRJavVJGE3wlPj5+PiQRYsWMR3GO7148YIsWbKEsNls0rdvX7JmzRqSk5PDdFgiBQUFZNeuXcTc3JwoKCiQ6dOnk02bNhEAREFBgbBYLKKmpkb8/f3J+fPnCZfLZTpksfvjjz+IoqIiSUtLk0p9dXV1hM1mkzNnzkilPimgCZ4Svy+++IIMGzaM6TBaVFxcTBYsWEAUFBSInZ0dOXXqFOHz+UyH1Soej0eOHj1KrKysRC326dOnk4sXL7bYhdNd8Hg8YmVlRRYsWCC1OuPj4wkAkpKSIrU6JewUHQdPiZ2jo2ObJwGljRCCkJAQ2NjY4MaNGzh+/DiePHmCmTNntjq7VRYoKipizpw5SElJwZ49e6Cnp4fw8HCUlpZCWVmZ6fAk5siRI8jMzJTq6J+kpCSoqKh0i40+hGT3yKbklqOjI7hcLp4/f850KACAkpISTJgwAf/5z38wf/58pKamwt/fX6YT+9sUFRXxxRdf4OnTp5g1axYWL14MPz8/vHr1iunQxK6xsRHbt2/H/PnzpbrLVHJyMmxsbNpct0ieyM8RTskNBwcHsNlsJCYmMh0KoqOj4ejoiJSUFERGRuL777+HhoYG02F1mra2Nvbs2YOwsDDEx8djyJAhuH//PtNhidXvv/+O3NxcfPPNN1Ktt5tMbmqCJnhK7FRUVGBlZcV4gg8NDYWvry+GDRuGR48ewd3dndF4xMnb2xuPHj2Cg4MD3n//fVy+fJnpkMSisbERu3btwuLFi2FhYSHVupOSkmiCp6j2cHR0ZDTB//bbb/jkk0+waNEihIaGyvUGF63R09PDxYsXMWvWLPj5+eHEiRNMh9Rlv/76K/Lz87F27Vqp1ltVVYXc3Nxul+C7T2cTJVMcHR2xb98+RuoODQ3F0qVLsX79egQFBTESg7Sw2Wz8+uuv0NHRwfz586Gjo4Px48czHVan1NfXY+fOnViyZAlMTU2lWndSUhIIId0uwdMWPCURgwcPxsuXL1FaWirVemNjY/Hpp59i+fLl3T65v+nbb7/F3LlzMX36dMTHxzMdTqeEhISgtLQUa9askXrdSUlJ0NDQQP/+/aVetyTRBE9JhKOjI4DXqx5KS2VlJWbNmoUPPvig3StCdhcsFgsHDhyAu7s7/P39m2zbJw/q6urw7bffYvny5ejbt6/U6xeeYJXn7QxbQhM8JREmJiYwNjaWamvyP//5D+rr63H48GG5GgIpLmw2G8eOHUN5eTm+/PJLpsPpkF9++QWvXr3C119/zUj9CQkJcHJyYqRuSep5nwJKapydnaWW4GNiYnDixAmEhITAwMBAKnXKoj59+mD//v04fPgw4uLimA6nXbhcLr7//nv85z//gYmJidTrJ4Tg8ePHol+d3QlN8JTEuLi4SCXJEELw9ddfY/To0ZgwYYLE65N106ZNg7u7OwIDA5kOpV1+/vlnVFZWMhbvixcvUFVVJfGtAJlAR9FQEuPs7Ixt27ZJfG34mJgY3L17V6Ym/BBCcPjwYfzzzz+wsrJCUVERRo0ahVmzZkml/h07dmDkyJGIj4+Hs7OzVOrsjNraWvyf//N/8Pnnn3d44xNxSUhIgIKCAuzt7RmpX5JogqckxsXFBYQQJCYmwtvbW2L1HD58GM7OznB1dZVYHR0VHByMw4cP49GjR9DR0UF5eTmcnJxQUlKClStXSrx+Hx8f2Nrail4bWbV3715wOBysWrWKsRgSExNhZWUl1zOcW0O7aCiJMTMzg6GhoUS7aXg8Hs6dO4eAgACJ1dFR2dnZCA4OxtKlS6GjowMA0NHRweLFi7Fu3TqUlZVJJY558+bhzz//BJHRXTlramrw448/4ssvv2T0vElCQkK37J4BaIKnJEzSJ1oTExNRVVUFX19fidXRUSdOnACPx8Po0aObXD9q1ChwOBwcOnRIKnGMHj0aZWVlSE1NlUp9HfXTTz+hvr4eX331FaNxJCQkdMsTrABN8JSEOTs7S7QFHxsbCx0dHVhbW0usjjf99ddf0NXVBYvFwoYNG0TX79+/HwoKCggJCUF0dDQANJuNaWZmBgBSW8LB0dER6urquHPnjlTq64iqqirs3r0bX331FXR1dRmLo6ysDC9fvqQteIrqDBcXF6SlpUls4k12djYGDRoktXHv06dPF82Q9fDwEF0/YcIEzJw5E0uWLEF+fj4AiLpnhISJLDMzUyqxstlsWFpaIjs7Wyr1dcTu3bvB5/Olcj7iXRISEgCAJniK6gxnZ2cIBAKJtVrLysqgp6cnkbJbs3TpUpiZmeGXX34RXffrr7+KhvlpaWkBQLNZkcL/NzQ0SClSQF9fX2p9/u1VWVmJPXv24L///W+zL0FpS0hIgKGhIYyNjRmNQ1JogqckytzcHPr6+hLrpuFwOFBVVZVI2a1RVlbGypUrcfnyZbx48QINDQ1IS0sTzYS0sbEBAFRUVDR5XHl5OQBIdTKPurq6zC1b8MMPP4AQgs8//5zpUJCYmNgtZ7AK0QRPSZyTk5PETrQKhyBK26JFi6Curo59+/bhwoULmDZtmug24XhqYVeNUEFBAQDA09NTanG+evVK6r9w3qWiogJ79+7F6tWrZWIJ58TExG57ghWgCZ6SgqFDh+LBgwcSKVtfXx8lJSUSKftdtLW1sWjRIhw+fBinT5+Gn5+f6LY5c+ZAW1sb4eHhTR4TFhYGJSUlqU12Al5vVyhLCf67776DoqIiPvvsM6ZDAYfDQUpKikzPE+gqmuApiRs2bBhSU1ObdVmIg729PdLS0sDhcMRedlu++OIL1NTUwMnJCUpKSqLrdXV1sXbtWhw8eBDV1dUAXo8aCQkJwYYNG0SjaSStqqoKGRkZMrPGeVlZGfbt24fVq1dDU1OT6XCQmJgIHo+HoUOHMh2KxNCZrJTEubm5gRCCuLi4ZmPDu2rEiBFobGzEw4cPJTpbtiUWFhb4/PPPsXz58ma3rV69Gvr6+lixYgX69euHZ8+eYfXq1Vi0aJHU4ouNjQWfz5eZrQq/++47qKioYMWKFUyHAgB4+PAhdHR0pLqxt7TRBE9JXJ8+fWBmZoZ79+6JPcH369cPAwcOxIULF6Se4IHXk3VawmKxsHDhQixcuFDKEf1/Fy5cgIODA4yMjBiLQai0tBT79+9HUFCQzCwJEBcXBxcXl263BvybaBcNJRVubm64d++eRMoOCAjA0aNHUV9fL5Hy5RGXy8Wff/6J+fPnMx0KAGDnzp3Q0NDAsmXLmA5FRJjguzOa4CmpGDZsGO7evSuRsufNm4fKykqcPHlSIuXLo8OHD4PL5WL27NlMh4LCwkIcOHAAa9euhZqaGtPhAHh9gvXp06c0wVOUOLi5uaG4uBg5OTliL9vU1BSLFi3Cxo0bGTnZKmtqamoQHByMzz77jLEleN+0a9cuaGtrY/HixUyHIvLo0aNuf4IVoAmekpKhQ4eCzWZLrJtmy5YtqKqqwo4dOyRSvjzZsmUL6uvrsXbtWqZDQUFBAUJCQrB+/XqpT0h7l4cPH0JXVxfm5uZMhyJRNMFTUqGmpgZ7e3uJJXgjIyN8++232LlzZ7Px5z3JjRs3sHv3buzevVsmxr9v374d+vr6Uh091B5xcXEYOnRotz7BCtAET0nRsGHDJJbgAWD58uWYNGkS5s6di7y8PInVI6uysrIQEBCAGTNmYN68eUyHg9zcXBw6dAjr16+HiooK0+E00RNOsAI0wVNS5Obmhri4ODQ2Nkqsjt9++w1aWlr48MMPGVnCgCklJSUYO3YsjI2NceDAAabDAfC69W5kZCQzI3mEampqkJaWRhM8RYmTm5sbuFwukpOTJVaHrq4url+/jqqqKowbN07mVlKUhKKiIowdOxY8Hg/Xrl0TrWbJpJycHPz+++/YuHEjlJWVmQ6niUePHoHP53f7E6wATfCUFNna2kJTU1Oi3TTA61E1N2/eRHFxMTw8PKS2/joT0tPTMWLECNTU1ODmzZsys+xtcHAwTExMMHfuXKZDaSYuLg76+vro378/06FIHE3wlNQoKirC1dVVKjsMWVlZ4c6dO1BTU4ObmxuuXbsm8Tql7eLFixg+fDj09fURExMjM1Pus7KycPToUWzevFnmWu8AcPfuXbi5uTEdhlTQBE9JlYeHB2JiYqRSl7GxMSIjIzF27Fh89NFH+Prrr1FXVyeVuiWJy+Vi5cqVmDJlCvz8/BAeHs7optVvCwoKQr9+/WRiklVL7t69i+HDhzMdhnQQipKif/75hwAg+fn5Uq33jz/+IJqammTgwIHk2rVrUq1bnC5evEgsLCyItrY2OXnyJNPhNJOenk7YbDY5duwY06G0KD8/nwAgN2/eZDoUaThFW/CUVLm7u0NRURGxsbFSrTcgIAApKSlwcnLCuHHjMHnyZDx69EiqMXTFgwcPRHG7u7sjNTUV/v7+TIfVzNatW2FhYYGZM2cyHUqL7t69CwUFBbi6ujIdilTQBE9JlZaWFuzt7aXWTfMmU1NTnDlzBjdu3EB+fj5cXFwwZcoUqZwT6KzIyEh89NFHGDZsGCoqKhAeHo4TJ06gT58+TIfWTHp6Ok6dOoWgoCCw2bK5UO29e/dgZ2cHbW1tpkORCprgKamTZj98S3x9ffHgwQNcvHgRBQUF8PDwwODBg/Hzzz/j1atXjMUlVFZWhv/973+wt7eHj48PKioqcO3aNcTGxsLHx4fp8Fq1efNmWFlZYcaMGUyH0qoe1f8O0D54SvqOHz9OlJSUSG1tLdOhEEIIuX//Plm4cCFRV1cnbDab+Pr6kgMHDpDc3FypxZCVlUV+/vlnMnr0aMJms4mmpiZZsmQJiYuLk1oMXZGcnEwUFBTImTNnmA6lVTwej2hoaJBff/2V6VCk5RSLEEKY/pKhepasrCxYWFggIiKCkU06WlNVVYUrV64gNDQU165dQ21tLSwtLeHl5QUPDw84OjrCzs6uyxtWVFVVISUlBYmJiYiJiUFkZCSys7OhqamJ8ePHY+rUqRg/fjzU1dXF9Mwk75NPPkFKSgoeP34MBQXZ7Bh49OgRnJ2dkZSUJNoYvZv7UzY7yqhuzdzcHKampoiJiZGpBK+lpQV/qoR73QAAIABJREFUf3/4+/uDy+UiNjYWkZGRiIyMxNmzZ1FbWwsWiwVzc3P069cPZmZmMDY2hr6+Pnr37g0Wi4XevXuDEILKykoIBAJUVlaiuLgYRUVFyM3NRU5ODrKysgAAGhoacHNzw4IFC+Dj4wM3Nzf06tWL2RehE5KTk3Hu3Dn89ddfMpvcgdfdM1paWrC1tWU6FKmhLXiKETNmzEBtbS0uX77MdCjtIhAIkJmZiaSkJCQnJ+Ply5fIy8tDQUEBysrKRAm9vLxclOgVFBSgra0NfX19sFgs5OXlYcWKFXBwcICDgwPMzc27xWqGU6dORXp6OhISEmQ6wQcEBKCgoAA3btxgOhRpoS14ihkeHh7YsmULBAKBTCcFIQUFBQwYMAADBgzA5MmTO/z46OhoeHl5YerUqbCyspJAhMx48uQJLly4gPPnz8v8+3j37l2ZPgEsCbL9jlDdloeHB8rLy/H06VOmQ5GK4cOHQ1tbGzdv3mQ6FLHasGEDnJycMHHiRKZDeadXr14hPT29xyxRIEQTPMUIR0dHqKurMzpcUprYbDZ8fHy6VYKPj4/HpUuXsHXrVpnvahIucEcTPEVJAZvNhru7OyIjI5kORWp8fX0RFhYm0fXwpWnTpk1wcXHBuHHjmA6lTdHR0bCysoK+vj7ToUgVTfAUY3x8fBAWFsZ0GFIzZswYVFdX4/79+0yH0mUPHz7E1atXsW3bNplvvQNAVFQUvLy8mA5D6miCpxgzcuRI5Ofn4/nz50yHIhVWVlawsLDoFt00GzZsgLu7O8aOHct0KG2qr6/HgwcP4OnpyXQoUkcTPMWYYcOGQV1dHbdv32Y6FKkZPXq03A/Tu3PnDq5fv46tW7cyHUq7PHjwAHV1dbQFT1HSpKysDHd3d0RERDAditT4+vri/v37cr1f7KZNm+Dh4YHRo0czHUq7REVFoU+fPjKzIYo00QRPMWrkyJE9qh9+9OjRIITI7a+WmJgY3Lp1C9u2bWM6lHaLjo6WqRnT0kQTPMWontYPr6enB2dnZ7nth9+wYQNGjRqFkSNHMh1KuwgEAty5c6dH9r8DAJ3JSjFK2A8fHh6OgQMHMh2OVIwZMwanT59mOowOCwsLw+3bt+WqS+3JkyeoqKjokf3vAG3BUwxTUlLCiBEj5CppdJWvry9evHiBjIwMpkPpkKCgIIwZM0auujuioqKgra0NBwcHpkNhBE3wFON62nj4ESNGQENDQ666aW7evInIyEhs3LiR6VA6JCoqCh4eHlBUVGQ6FEbQBE8x7v3330dBQQHS09OZDkUqlJWV4e3tLVcJPjg4GOPGjZO7vuye3P8O0ARPyQBXV9ceNx7e19cXt27dAp/PZzqUNl27dg1RUVHYvHkz06F0SEZGBl6+fNlj+98BmuApGaCkpAQPD48e1Q8/ZswYVFRU4MGDB0yH0qagoCBMnDhR7hbqioqKgoqKCoYOHcp0KIyhCZ6SCe+//z7CwsLQU/afsbOzg5mZmcx301y6dAn379+Xu753AIiIiJDbXbLEhSZ4SiZ88MEHKCgoQFJSEtOhSM3o0aNlPsEHBwdj8uTJcHV1ZTqUDgsPD8f777/PdBiMogmekgnOzs4wMDCQ+YQnTr6+voiNjUVVVRXTobTowoULePjwodz1vQPAixcvkJWVhVGjRjEdCqNogqdkgoKCAkaNGtWjEvyYMWMgEAhk8uQyIQRBQUGYOnUqhgwZwnQ4HRYWFgY1NTW5O28gbjTBUzLD19cXERER4HK5TIciFfr6+nB0dJTJL7Vz587h8ePHctn3DrzunvH09ISKigrToTCKJnhKZowdOxZcLhd37txhOhSp8fX1lbkELxAIEBwcjE8++QSDBw9mOpwOI4QgPDy8x3fPADTBUzLE1NQU1tbWMpfwJMnX1xdpaWnIyspiOhSRM2fOICkpCevXr2c6lE5JTk5GYWEhTfCgCZ6SMWPGjJH7DTE6wtPTE2pqavj333+ZDgXA69b79u3b4e/vL7frt4SFhUFbWxvOzs5Mh8I4muApmeLr64uEhAQUFxczHYpU9OrVC15eXjLzq+XUqVNITU3Fhg0bmA6l08LDwzFy5Mgeu/7Mm2iCp2TKyJEjwWazZaZFKw3CZQsEAgGjcfD5fGzbtg2zZ8+GjY0No7F0Fp/PR0RERI8f/y5EEzwlUzQ1NeHu7i4zLVpp8PX1RVlZGeLj4xmN4/jx40hPT8e6desYjaMrHj16hPLyctr//v/QBE/JHOHIkp6ybMF7772HPn36MHrugc/nY8eOHQgICICVlRVjcXRVWFgYDA0N5fb8gbjRBE/JHF9fX+Tl5SE1NZXpUKSCxWLhgw8+YPRXy5EjR5CZmSm3I2eEhMsTsFgspkORCTTBUzJn6NCh0NXVxfXr15kORWp8fX1x584d/N/27j6s5vv/A/jznE4nJbFCN27aN1ZKEV/GpFrSsAhhskqsmJuxbzYyuvgxIzNf5TYjrUxcEatQGuXQupGbqDXL0A1RVNL9qXM+vz/89FuTu/qc8zk3r8d17fpe6ni/n+X6Pvv0Pp/P+11TUyP3uZuamvDdd99h7ty5MDMzk/v8bBGLxUhNTaX197+hgicKR0NDAx999BFOnz7NdRS5cXFxQVNTEy5cuCD3ucPDw1FcXIyVK1fKfW42paamoqamBi4uLlxHURhU8EQhubq6QiQSoaqqiusocmFkZARra2u5L9M0NTUhKCgI8+bNw7/+9S+5zs22M2fOwMLCQql/C2EbFTxRSB9//DEYhlGru2naesirubkZlZWVMptz3759KCkpwTfffCOzOeQlMTER48eP5zqGQqGCJwpJX18fI0eOxKlTp7iOIjcuLi7Iy8uDSCTCrl27MGnSJHTt2hVDhgzp8NgFBQWYPXs20tPTWz7W2NiITZs2Yf78+ejdu3eH5+DSgwcPkJOTg3HjxnEdRaEIuA5AyMu4urriv//9L6RSKfh81b0WKS8vR3JyMhISEqCnp9fysBfw7AqejdtFz58/j4MHD+LgwYNwcXHBhg0bkJmZicePHyMgIKDD43MtISEBWlpacHR05DqKQqGCJwpr4sSJ+Oabb5CVlaWS+3rfuHEDc+bMQXZ2NoBnZ9OKxWIAz4r9OQMDgw7Pdfv2bWhpaaGxsRHnz5/HiBEjYGJiAnd3d/Tq1avD43PtzJkzcHR0hI6ODtdRFIrqXhYRpWdtbY13331XZZdppFIprl+/DoZhwDBMS7n/U8+ePTs81507d9DU1AQALf9bVlaGqKgojBkzBleuXOnwHFyRSCQ4e/YsLc+0gQqeKDRXV1ecPHmS6xgyYWtri9WrV79yUyw+n89Kwf/5558v7HXz/LeEixcvYvjw4Zg6dSoePXrU4bnkLTMzExUVFfQGaxuo4IlCc3V1RXZ2Nu7du8d1FJlYs2YNBg0aBE1NzTY/LxAIoK+v3+F57ty589LPPS/6kydPorCwsMNzyduZM2fQp08fWFpach1F4VDBE4U2ZswYdO7cWWWXaQQCAQ4dOvTSR+t5PF6HC76mpuaVt1ry+XwIBAIcPXoUw4YN69BcXDhz5gwmTJjAdQyFRAVPFJqWlhacnZ1VtuABwNLSEhs3bmyz5KVSaYffZH3V1Tufz4empiZOnTqFKVOmdGgeLpSXl+Py5cu0/v4SVPBE4bm6uuLs2bOoq6vjOorM+Pv7w87O7oWlGolE0uEr+JcVvIaGBrS0tJCQkKC0j/cnJSWBx+PB2dmZ6ygKiQqeKLyJEyeioaEBKSkpXEeRGT6fj4MHD75Q8FKplJWCFwqFrT6moaEBbW1tJCcnK/XmXGfOnMEHH3yArl27ch1FIVHBE4VnbGyMoUOHIi4ujusoMvXuu+8iJCTkhaUaNpZo/v6wlEAggJ6eHlJTUzFy5MgOjc0liUSChIQEfPzxx1xHUVhU8EQpuLu74/jx460eAFJFfn5+cHNza3Ul39Er+Fu3brXc+66pqYlu3brhwoULGDx4cIfG5Vp6ejrKysqU8r0DeaGCJ0ph+vTpePz4MX777Teuo8hcaGgoOnfu3HIl/84773RovPz8fADPrty7d++O9PR0lTjxKDY2Fv3791fa82PlgbYqIErB3NwcAwcORExMjMrvN2JkZITw8HBMnToVGhoaqKurQ0FBAZ48eQKJRILq6mo0NzdDR0cHWlpa0NbWRrdu3WBiYvLCDwOpVIr79+8DAHr37g2RSIS+ffty8WWxLi4ujq7eX4MKniiNadOmISwsrM11amVXX1+PS5cuIScnB7m5ucjNzYW2tjbq6+vfaqdHbW1t9OrVCxYWFrC2toaJiQmamppgbm6O8+fPw9jYWIZfhfzcvHkT+fn5mDx5MtdRFBqPUZeTjYnSu379OmxtbZGenq7Ubw4Cz66sMzIycPr0aYhEIly6dAlisRj6+vqwsbGBlZUVzMzMoKGhgeHDh8PQ0BD6+vrg8/no0qULBAIB6urq0NjYiIaGBlRUVKCkpAQlJSUoLi5GXl4e8vLy8Mcff0AsFqN79+5wdHSEk5MT3Nzc0KdPH66/BR2yefNmbNmyBaWlpa/c6kHNHaGCJ0rFwsICkydPxvfff891lHZJTU3FoUOHEBsbiwcPHqBfv35wcnKCo6MjHBwcWF8+aW5uxrVr13DhwgWIRCKIRCJUV1dj2LBhmDZtGmbPnq2UV/V2dnawsLDAgQMHuI6iyKjgiXIJCAjA0aNHX/l0pqKprq5GZGQkQkNDkZubi8GDB2PatGmYOnWq3N/sFIvFOHfuHE6cOIETJ06gqqoKbm5uWLBgAcaOHSvXLO1VVlYGExMTxMTE0BLNqx0BQ4gSyczMZAAw165d4zrKa1VXVzPBwcGMoaEh06lTJ2bGjBnMr7/+ynWsFo2NjUx0dDQzduxYhsfjMYMHD2aio6MZqVTKdbRX2r9/P6Otrc3U1NRwHUXRHabbJIlSGT58OPr27YuYmBiuo7yURCLBrl27YGpqirVr12L+/PkoKSlBdHS0Ql0lC4VCzJgxA7/++isuXbqEvn37YubMmRg1apRC7w8fGxuLsWPHonPnzlxHUXhU8ESp8Hg8TJkyBUePHuU6SpsuX76MESNGwN/fH35+frh79y7Wr1/f4XvZZW3YsGGIi4vD5cuXIRQKMWLECCxevBhVVVVcR2ulvr4e586do6WZN0QFT5TOtGnT8Oeff+KPP/7gOkoLqVSKoKAgjBo1Crq6usjOzsbmzZsVvtj/aejQoTh//jzCw8MRExMDW1tbZGRkcB2rRVJSEhoaGuDq6sp1FKVABU+UzujRo2FsbKwwyzQVFRUYN24c1q5di02bNiElJQVWVlZcx2o3Ho8Hb29v5OTkwNLSEg4ODti6dSvXsQA8W54ZOXIkjIyMuI6iFKjgidLh8/lwc3PDsWPHuI6CwsJCjB49Grdu3cJvv/2Gr776SmUewurRowdOnTqF7777DgEBAViyZMkLx/7JU3NzM06ePAk3NzfOMigbKniilGbOnInr16/j999/5yxDfn4+Ro0aBaFQiLS0NKU8Del1eDweli9fjujoaOzfvx+zZs2CRCLhJEtKSgoePXqEadOmcTK/MqKCJ0rJ0dERvXv3xuHDhzmZv6SkBOPGjUPfvn0hEolgYmLCSQ55cXd3R2JiIuLj47Fw4cJW2w/LS3R0NIYNG4b+/fvLfW5lRQVPlBKfz4eHhwd+/vlnuZdNbW0txo8fDx0dHZw8eVJtDptwdHREdHQ0wsPDsWHDBrnO3dzcjNjYWHzyySdynVfZUcETpeXp6YnCwkKkp6fLdV5/f3+UlJQgISGhw4dxKJuJEyciJCQE69atw4ULF+Q279mzZ2l5ph1oqwKi1GxsbODg4IBdu3bJZb5ffvml5fARdd6q1t3dHZcvX0ZOTo5cfoP57LPPkJeXp1C3bCqBI3QFT5TarFmzEB0d3XJikSw1NjZi2bJl8Pb2VutyB4B9+/ahrq4OQUFBMp+rqakJsbGxmDFjhsznUjVU8ESpeXp6ory8HL/++qvM59q5cydKS0uxceNGmc+l6AwMDBAYGIiQkBAUFxfLdK6kpCRUVlZi+vTpMp1HFdESDVF6o0ePhqmpKQ4dOiSzOSQSCd59913MmjWLla2Kt2zZgu+++w5VVVXg8/lwdnaGUCgEwzBoaGjArVu3UFxcjP379+PYsWNITEwEAHz44YcAnu1QaWJiAjc3N3h7e0NLS6vV+AzD4MCBA0hMTIS5uTlKS0sxZswYfPrppx3O/pxYLEa/fv3g5eWFTZs2sTbuP/n4+ODWrVtIS0uT2RwqinaTJMpv9+7djI6ODlNdXS2zOU6fPs3weDzm1q1brI15//59BgDTv3//Fz4nkUgYV1dX5q+//mLu3bvHAGDefffdVp+PjY1lzMzMmP79+zO5ubmt/v66desYU1NTpqKigmEYhqmoqGBMTU2Z4OBg1vIzDMMEBgYyRkZGTFNTE6vjPtfY2Mi88847zLZt22Qyvoo7TAVPlF55eTkjFAqZn3/+WWZzzJ49mxk9ejSrY0qlUgYAY2Fh0ebnU1NTmYcPHzIMw7z0dSUlJYyRkRHTr18/pq6ujmEYhikoKGAEAgGzcePGVq/dsGEDo6Ojwzx+/Ji1r+H27dsMAObs2bOsjfl3sbGxDI/HYwoLC2Uyvoqj7YKJ8tPX14eLiwuioqJkNsfFixfh4uLC6piv2tIgOzsbdnZ2MDQ0fOUYxsbG+Pbbb3H79u2W/WIOHTqE5uZmODs7t3rtmDFjUFdXh/3793c8/P8xMzODmZkZUlNTWRvz76Kjo2FnZ6cyB4XLGxU8UQmenp5ISkpCaWkp62OXlpbi7t27+OCDD1gf+5/EYjFycnKwZMmSN/4706dPB5/PR1JSEgC0lO0/D+t+fg7r9evXWUr7zKhRo2TyLEJDQwPi4+Pp4aYOoIInKmHy5Mno1KkToqOjWR+7sLAQwLPzYGXhzz//BI/HA4/Hg5aWFgYNGoScnJw3/vvdunVDz549W/blKSkpAYAXtirW19cHANy9e5el5M+Ym5ujoKCA1TEBICEhATU1NfRwUwdQwROVoKOjgxkzZsjkEObHjx8DgMyeWrWwsADDMGAYBs3NzcjPz4epqelbjSEQCFrupNHT0wPw4hLQ8z+LxWIWUv8/AwMDlJeXszomAERGRsLJyUnl9/mRJSp4ojJ8fX2RnZ3N+nFz9fX1AIBOnTqxOm5bNDQ08N5772Hx4sVv/HfEYjFKS0tha2sLABgwYAAA4MmTJ61eV1lZCQCsF6auri5qa2tZHbO8vBwJCQnw8fFhdVx1QwVPVIadnR2srKwQFhbG6rjPlzr+WZiyNH/+/Dd+bXJyMpqamlreVB04cCCA/1+qee7BgwcAnj03wKby8vKW5R+2REVFQSAQYOrUqayOq26o4IlKmTt3LqKiolBXV8famM+XZh49esTamABY2QWzsbERq1atgq2tLZYuXQoA8Pb2RteuXZGSktLqtcnJydDU1GT1YSfg2feF7eWriIgIzJgxA7q6uqyOq26o4IlKmTNnDhoaGlg9lNvc3BxCoRDXrl1jbUwALcsar/th9PzzDQ0NrT5+9epVuLi4oLKyEocOHYKmpiaAZ2+mfvPNN9i7dy+qq6sBAE+fPsWPP/6IwMDAlrtp2HL16lXY2NiwNl5eXh6uXLlCyzMsEHAdgBA2de/eHZMmTUJYWBhrBaGtrQ1bW1ukpaVh1qxZrIyZlpbWspRUXFyMlStX4pNPPsHQoUNbvS41NRXh4eEAnt3N8+GHH0JLSwtaWlrQ1NTEzJkz4ePj88KV7ooVK9C9e3csWrQIffv2RX5+PlasWAE/Pz9W8j8nlUqRmZnJ6v7wP/30E0xNTeHg4MDamOqK9qIhKicxMRETJkxAXl4eLC0tWRlz1apViIqKwp07d8Dn0y++zyUnJ8PZ2Zm177VUKoWpqSnmzp2L9evXs5BQrdF2wUT1fPTRRzA1NW258mXDnDlzUFRUhHPnzrE2pioICwvDyJEjWftBmpSUhHv37sHT05OV8dQdFTxROXw+H3PmzEFERARr+8Sbm5vDzs4O27ZtY2U8VVBUVITjx4+zuuwTEREBOzs7mT1Upm6o4IlK8vPzQ3l5OeLi4lgb89tvv0VCQgLOnj3L2pjKLDAwEMbGxvDy8mJlvKdPnyIuLo7eXGURrcETlTVhwgTweDycPn2atTFdXV1RXFyMrKysF/ZgVydpaWmwt7dHVFQUZs6cycqY+/btw9KlS/HgwQN069aNlTHVHK3BE9Xl6+uLM2fOtOwlw4YdO3agsLAQK1euZG1MZVNVVQUvLy+MGzeO1Y3AIiIiMGXKFCp3FlHBE5Xl5uaG7t27s/pmq5mZGXbv3o2QkBAcO3aMtXGVhUQigY+PDxoaGhAREfHKLY/fRn5+PtLS0mh5hmVU8ERlCYVC+Pr6Ys+ePWhsbGRtXE9PT3zxxRfw8vJ64WlRVbdo0SIkJSXh6NGj6NGjB2vj7t27F3369GF9z311RwVPVNoXX3yByspKVp9sBYDg4GBMnToVU6ZMgUgkYnVsRSSVSuHv748DBw7gyJEjsLOzY23s+vp6RERE4PPPP4eGhgZr4xIqeKLiTExM4O7uzvrtjXw+HxERERg/fjzGjRsnk33oFUVjYyM8PT2xe/duHDx4EG5ubqyOHx0djadPn2LOnDmsjkuo4IkaWLp0Ka5evYq0tDRWxxUKhTh8+DAWLFiAWbNmYdWqVWhubmZ1Dq4VFBTgww8/REJCAhISEuDh4cH6HHv27IG7uzvt+y4DVPBE5Y0aNQrvv/8+tm/fzvrYfD4fwcHB2Lt3L0JCQuDg4IA7d+6wPg8XoqOjMWTIENTW1iI9PR1jxoxhfY7r168jMzMTCxcuZH1sQgVP1MSSJUsQExOD4uJimYzv5+eHrKws1NTUwNraGhs2bGD1jV15un37NlxdXeHh4QEPDw9kZmaythXBP+3atQsDBgygjcVkhAqeqIWZM2fC0NAQoaGhMpvDysoKV65cwdq1axEUFAQbGxtERUVBIpHIbE42lZWVYcWKFbC2tkZRURFSUlKwZ88eaGtry2S+6upqHDlyBAsXLmTtdkvyDwwhauJ//ud/GH19faa2tlbmcxUVFTGzZ89mBAIBM2DAACYiIoKpr6+X+bztUVhYyHz99ddM586dGUNDQyY4OJgRi8Uyn3fHjh2MtrY2U1FRIfO51NRhKniiNkpLS5lOnTox+/fvl9uc+fn5zOzZsxlNTU3GwMCAWbZsGfP777/Lbf6XEYvFTFxcHDNx4kRGQ0ODMTY2ZrZu3SqXH37PWVtbM76+vnKbTw0dpr1oiFrx8fHB1atXcePGDbkuCzx8+BAHDhzAvn37UFBQAEtLS0ybNg1TpkyBra2tXO7/rq6uRnJyMo4fP474+Hg8efIEzs7O+PzzzzF58uSWE6Hk4cKFC3B0dERWVhaGDRsmt3nVzBEqeKJWrl27hqFDhyI5ORlOTk5yn18qlSI1NRXHjx/HiRMnUFRUhK5du2L06NGwt7fH0KFDMXDgwA7fMtjc3Iz8/Hzk5uYiIyMDFy9eRHZ2NqRSKUaNGgV3d3e4u7vD1NSUpa/s7cyaNQu3b9/GpUuXOJlfTVDBE/Vjb2+P7t2748SJE5zmYBgGubm5EIlEuHjxIlJTU1FSUgLg2bmq5ubmMDY2Ru/evWFoaAg9PT1oaWlBR0cHWlpaqK6uRnNzM6qrq1FVVYV79+6htLQUhYWFuHXrFsRiMQQCAaysrODo6AgHBwc4ODigZ8+enH7djx49Qp8+fbBnzx7MnTuX0ywqjgqeqJ9jx47Bw8MDeXl5MDc35zpOK+Xl5bhx4wby8vLw119/4eHDh7h37x7Kysrw9OlTNDY2oqamBk1NTdDV1YWmpia6dOkCPT099OrVC0ZGRujTpw8GDBiAgQMHwtLSUuG2Nf72228RHByM4uJi6OjocB1HlVHBE/UjkUhgZWUFe3t77N+/n+s4b+3w4cPw9vZWyqdmGxoa8K9//Qu+vr6sHtRN2kT7wRP1o6GhgRUrViAyMpLVveLJ64WHh6OyshKLFy/mOopaoIInamn27NkwMTGhM1blSCqVYtu2bfDx8YGxsTHXcdQCFTxRS5qamvD398e+fftQVlbGdRy1cOLECfz111/4z3/+w3UUtUEFT9TWvHnzoKuri61bt3IdRS1s3rwZkydPltm+NuRFVPBEbeno6CAgIAA7duzA/fv3uY6j0k6dOoWsrCysWrWK6yhqhQqeqLVFixbBwMAAmzdv5jqKSlu/fj3c3NwwfPhwrqOoFSp4otY6deqEVatWYe/evbh79y7XcVTSqVOncOnSJQQGBnIdRe1QwRO15+fnh969e2Pjxo1cR1FJdPXOHSp4ovY0NTWxZs0a/PTTT7h58ybXcVRKfHw8srKysGbNGq6jqCUqeEIAeHt7w8bGBitXruQ6isqQSCT45ptvMH36dPz73//mOo5aooInBM/OVv3hhx8QGxuLc+fOcR1HJYSHhyM/Px/fffcd11HUFhU8If9nzJgx+Oijj7B8+XJIpVKu4yi1+vp6rF+/HvPnz8d7773HdRy1RQVPyN9s2bIFN27cwJEjR7iOotSCg4NRWVlJd85wjAqekL8ZNGgQfHx8sGrVKtTX13MdRyk9evQImzdvxldffQUjIyOu46g1KnhC/mHjxo148uQJPfzUToGBgdDV1cXXX3/NdRS1RwVPyD8YGhpi9erV+P777+nhp7d0/fp1hIWFISgoCLq6ulzHUXtU8IS04csvv4SpqSmWL1/OdRSl8p///AfDhw+Hp6cn11EIqOAJaZNQKMSOHTsQExODM2fOcB1HKRw7dgwikQg//PADeDwe13EIqOAJeamxY8di8uTJ8Pf3R2NjI9dxFFpdXR2WL18OLy8v2NnZcR2H/B8qeEJeYfv27SguLkZQUBBJ+aHzAAAQ0ElEQVTXURTaunXrUFFRgU2bNnEdhfwNFTwhr9C3b1+sW7cOGzduRF5eHtdxFFJOTg62bduGzZs3o1evXlzHIX/DYxiG4ToEIYpMIpFgxIgR0NHRgUgk4nx9+fDhw/D29kZzczOnOYBn56yOHj0aEokE6enp4PPpmlGBHKF/DUJeQ0NDA2FhYcjIyEBYWBjXcRTKnj17kJWVhb1791K5KyD6FyHkDQwePBhLly7F119/jZKSEq7jKIQHDx4gMDAQy5cvh62tLddxSBuo4Al5Q+vXr4eBgQG++uorrqMohCVLlqBbt25YvXo111HIS1DBE/KGdHR0sGvXLhw5cgTx8fFcx+FUQkICYmJisHPnTnTu3JnrOOQlqOAJeQvjx4/HrFmz8MUXX6CmpobrOJyoq6vD4sWL8emnn8LV1ZXrOOQVqOAJeUvbtm1DTU2N2h5DFxgYiIqKCmzZsoXrKOQ1qOAJeUuGhobYsmULtm/fjrS0NK7jyFVGRga2b9+OLVu2wMTEhOs45DWo4Alph88++wzjx4+Ht7e32izV1NXVwcfHB05OTvDz8+M6DnkDVPCEtNO+fftQVVWFgIAArqPIRUBAAEpLS3HgwAHOH/Yib4YKnpB2MjY2xs6dO7Fnzx4kJCRwHUemkpOTsWvXLuzatQt9+vThOg55Q7RVASEd5OHhgdTUVOTk5OCdd95hffyPP/4YxcXFLX+urq7GgwcPYG5u3up1CxcuxKJFi1ifv6qqCoMGDcKQIUPwyy+/sD4+kZkjAq4TEKLsQkNDYWNjgy+//BKRkZGsj3/37l3cvHnzhY/n5ua2+nN1dTXrcwPPDj+pr6/H3r17ZTI+kR1aoiGkg7p164awsDD8/PPPOHbsWKvP5ebmwtfXt0PlO3v2bGhqar7yNTweDzNnzmz3HDU1NVi3bh3Ky8tbfTw+Ph4RERHYs2cPDA0N2z0+4QhDCGGFr68v0717d+bBgwcMwzBMaGgoIxQKGQDMoUOH2j1uQUEBw+PxGABt/sfj8ZgRI0Z0KPv+/fsZAIyhoSFz4cIFhmEY5tGjR4yRkREzZ86cDo1NOHOY1uAJYUltbS0GDx4Mc3NzaGtr48SJE2AYBgKBAFOnTkV0dHS7xx4+fDiuXLmCtv7vqqGhge3bt3do/d3V1RVnzpxpGX/9+vXIzs5Genq6zN5bIDJ3hAqeEBaFhoZi5cqVqK2tbbVfu46ODiorKyEUCts17o4dO7Bs2bI294DX0NDA/fv3272EUldXB319/VbHEvL5fOjr62Pnzp0dWvohnKL94AlhA8MwCAkJwZIlS1BTU/NCEdfV1UEkErV7fA8Pjzav3vl8PsaMGdOh9fHExESIxeJWH5NKpXjy5AkWLFiAxMTEdo9NuEUFT0gHlZaWwtnZueUKWyKRvPAaoVCIuLi4ds/Ro0cPODo6QkND44XPeXt7t3tcAPjll18gELx4Q11zczOePn2Kjz/+GF9++SWampo6NA+RPyp4Qjpo1apVSElJgVQqfelrxGIxYmJiOjSPt7f3C1fxAoEAbm5u7R5TIpEgLi7upeUtlUrBMAy2b9+ODRs2tHsewg0qeEI66Pvvv8eMGTMA4JWP8D948ADZ2dntnsfd3b3VlbZAIMCkSZPQtWvXdo+ZmpqKqqqqV76Gx+PB0dER8+bNa/c8hBtU8IR0kIGBAaKjoxEdHQ09Pb2X3rPe0WUaPT09TJgwoaXkJRIJPD092z0e8Gx55mVv/GpqakIoFGLTpk1ITk5G7969OzQXkT8qeEJYMmPGDPz555/46KOP2rySZ2OZxsvLq2WNX1tbGxMmTOjQeMeOHXvhDVbg2Zu3gwcPRk5ODgICAuhAbSVF/2qEsMjQ0BAnT57ETz/9BG1t7Reu5nNycnDv3r12jz9x4kRoa2sDAD755BN06tSp3WPduHHjhSwCgQCamprYuHEjMjIyXtjvhigXKnhCZGD27NnIy8vD+++/3+rqV0ND463Oc5VIJCgpKcG1a9eQkpKC1NRUvP/++wCA/v374+zZs/jtt9/w+++/o7Ky8q0yxsbGtvoBpKGhgffeew+XLl1CQEBAm3fsEOVCDzoRIkNSqRRbt27F6tWrwTAMJBIJxo4di6SkpFavq6+vx6VLl5CTk4Pc3Fzk5ubizp07KC0tfeXdOf+kra2NXr16wcLCAtbW1rC2tsbQoUNhaWn5wrKRra0trl+/DoFAAIZhEBgYiNWrV7923xuiNOhJVkLkITc3F59++ilycnKgqamJx48fIzc3F6dPn4ZIJMKlS5cgFouhr68PGxsbWFlZwcLCAsbGxjAxMYGhoSH09fXB5/PRpUsXCAQC1NXVobGxEQ0NDaioqEBJSQlKSkpQXFyMvLw85OXl4Y8//oBYLEbPnj1hb28PJycnuLm5gc/no0+fPmAYBu+99x6ioqIwbNgwrr9NhF1U8ITIS1NTE+bNm4fIyEgYGBjg8ePH6NevH5ycnODo6AgHBwf07duX1Tmbm5tx7do1XLhwASKRCCKRCNXV1TAzM8OdO3ewYMEC/Pe//+3QWj5RWFTwhMhadXU1IiMjERoaitzcXFhbW+OTTz7B1KlTYW1tLdcsYrEY586dQ0xMDI4fP46amhq4ublhwYIFGDt2rFyzEJmjgidEVmpqahAWFoZNmzahqqoKkyZNwvz58xWmSMViMWJjY/Hjjz/i3LlzGDRoEFavXo3p06fTmauqgQqeELZJJBKEhoZizZo1kEgkWLp0Kfz9/RV6y93Lly9j/fr1OHnyJEaMGIGdO3fi3//+N9exSMfQbpKEsOny5csYMWIE/P394efnh7t372L9+vUKXe4AMGzYMMTFxeHy5csQCoUYMWIEFi9e/NptDIhio4InhAVSqRRBQUEYNWoUdHV1kZ2djc2bNyt8sf/T0KFDcf78eYSHhyMmJga2trbIyMjgOhZpJyp4QjqooqIC48aNw9q1a7Fp0yakpKTAysqK61jtxuPx4O3tjZycHFhaWsLBwQFbt27lOhZphxc3gSaEvLHCwkJMmDABdXV1+O2331TqXvIePXrg1KlT+OGHHxAQEICCggKEhITQvjRKhAqekHbKz8+Hk5MTevTogbNnz8LExITrSKzj8XhYvnw5+vXrB09PT5SVlSEqKoq2MVAS9KOYkHYoKSnBuHHj0LdvX4hEIpUs979zd3dHYmIi4uPjsXDhwjaPDySKhwqekLdUW1uL8ePHQ0dHBydPnuzQgRvKxNHREdHR0QgPD6fTnZQEFTwhb8nf3x8lJSVISEiAgYEB13HkauLEiQgJCcG6detw4cIFruOQ16AHnQh5C7/88gvc3d1x/PhxTJkyhes4nHF3d8fly5eRk5OjNr/BKCF6kpWQN9XY2AhLS0vY29sjIiKC6zicKi8vh4WFBebNm4dNmzZxHYe0jZ5kJeRN7dy5E6Wlpdi4cSPXUThnYGCAwMBAhISEoLi4mOs45CWo4Al5AxKJBMHBwVi8eDF69erFdRyFsGjRIhgYGGD37t1cRyEvQQVPyBtISkrC/fv3MX/+fK6jKAyhUIg5c+bgp59+QnNzM9dxSBuo4Al5A0eOHIGdnR369+/PdRSFMnfuXDx8+BAikYjrKKQNVPCEvIGLFy/CxcWF6xgKx8zMDGZmZkhNTeU6CmkDFTwhr1FaWoq7d+/igw8+4CzDzz//DB0dHfB4PAQFBbUsiRw6dAhCoZDTu3pGjRqF9PR0zuYnL0cFT8hrFBYWAgAsLCw4y+Dl5YVly5YBACZNmgSB4Nk2Uvb29nB1dYWPjw9n2czNzVFQUMDZ/OTlqOAJeY3Hjx8DAOdPrfr7+0NXVxfbtm1r+dihQ4fg6+vLYapn35fy8nJOM5C2UcET8hr19fUAgE6dOnGaw8DAAEuWLEFkZCTu378PhmFw7tw5jB8/ntNcurq6qK2t5TQDaRsVPCGv8fxUpidPnnCcBFi2bBmEQiG2bduGK1eu4P33329ZruFKeXk59PX1Oc1A2kb7wRPyGs+XZh49esT5Mk337t2xcOFChIaGoqysDGvWrOE0D6AY3xfSNrqCJ+Q1zM3NIRQKce3aNa6jAAC++uoriMViFBUVKcR9+VevXoWNjQ3XMUgbqOAJeQ1tbW3Y2toiLS2N6ygAACMjI7i4uHD+5irw7LDxzMxMTm8hJS9HBU/IG3B2dkZ8fDykUinXUVBbW4ubN29i2rRpXEfB+fPn8eTJE4wZM4brKKQNVPCEvIE5c+agqKgI586d4zoKdu3ahSVLlkBHR4frKAgLC8PIkSNhaWnJdRTSBtoPnpA3ZG9vjy5duuD06dNynzsjIwPz589HXV0dJBIJbt68CS0tLbnn+LuioiJYWFhg586dCrFcRF5A+8ET8qa+/fZbJCQk4OzZs3Kfu3Pnznj69Cn4fD6ioqI4L3cACAwMhLGxMby8vLiOQl6CruAJeQuurq4oLi5GVlaWQpQsV9LS0mBvb4+oqCjMnDmT6zikbXRkHyFv486dOxgyZAg+++yzVlsGqJOqqioMGTIEAwYMwKlTp8Dj8biORNpGSzSEvA0zMzPs3r0bISEhOHbsGNdx5E4ikcDHxwcNDQ2IiIigcldw9CQrIW/J09MTmZmZ8PLygoGBAZycnLiOJDeLFi1CUlISfv31V/To0YPrOOQ16AqekHYIDg7G1KlTMWXKFLU4zUgqlcLf3x8HDhxoOd2KKD4qeELagc/nIyIiAuPHj8e4ceMQHR3NdSSZaWxshKenJ3bv3o2DBw/Czc2N60jkDVHBE9JOQqEQhw8fxoIFCzBr1iysWrVK5Q6fLigowIcffoiEhAQkJCTAw8OD60jkLVDBE9IBfD4fwcHB2Lt3L0JCQuDg4IA7d+5wHYsV0dHRGDJkCGpra5Genk7bESghKnhCWODn54esrCzU1NTA2toaGzZsQGNjI9ex2uX27dtwdXWFh4cHPDw8kJmZSVsRKCkqeEJYYmVlhStXrmDt2rUICgqCjY0NoqKiIJFIuI72RsrKyrBixQpYW1ujqKgIKSkp2LNnD7S1tbmORtqJCp4QFmlqaiIgIAB//PEHPvjgA/j4+MDa2hqRkZFoaGjgOl6bioqKsHz5cpiZmSEyMhJBQUG4evUqHB0duY5GOoieZCVEhm7duoUNGzbg8OHD0NPTg4+PD3x9fWFlZcVprqamJiQmJuLHH39EQkICevbsia+//hoLFixQiF0qCStoqwJC5OHhw4c4cOAA9u3bh4KCAlhaWmLatGmYMmUKbG1toaGhIfMM1dXVSE5OxvHjxxEfH48nT57A2dkZn3/+OSZPngxNTU2ZZyByRQVPiDxJpVKkpqbi+PHjOHHiBIqKitC1a1eMHj0a9vb2GDp0KAYOHAgTE5MOzdPc3Iz8/Hzk5uYiIyMDFy9eRHZ2NqRSKUaNGgV3d3e4u7vD1NSUpa+MKCAqeEK4wjAMcnNzIRKJcPHiRaSmpqKkpAQAoK+vD3NzcxgbG6N3794wNDSEnp4etLS0oKOjAy0tLVRXV6O5uRnV1dWoqqrCvXv3UFpaisLCQty6dQtisRgCgQBWVlZwdHSEg4MDHBwc0LNnT46/ciInVPCEKJLy8nLcuHEDeXl5+Ouvv/Dw4UPcu3cPZWVlePr0KRobG1FbWwuxWAxdXV1oamqiS5cu0NPTQ69evWBkZIQ+ffpgwIABGDhwICwtLdV6W2M1RwVPCCEqirYLJoQQVUUFTwghKooKnhBCVJQAwFGuQxBCCGFd5v8CeBk7oV6QR6sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = CausalModel(df,\n",
    "                    data[\"treatment_name\"],data[\"outcome_name\"],\n",
    "                    data[\"gml_graph\"],\n",
    "                   missing_nodes_as_confounders=True,\n",
    "                   logging_level=logging.INFO)\n",
    "\n",
    "model.view_model()\n",
    "from IPython.display import Image, display\n",
    "display(Image(filename=\"causal_model.png\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Identifying the natural direct and indirect effects\n",
    "We use the `estimand_type` argument to specify that the target estimand should be for a **natural direct effect** or the **natural indirect effect**. For definitions, see [Interpretation and Identification of Causal Mediation](https://ftp.cs.ucla.edu/pub/stat_ser/r389-imai-etal-commentary-r421-reprint.pdf) by Judea Pearl.\n",
    "\n",
    "Natural direct effect: Effect due to the path v0->y\n",
    "Natural indirect effect: Effece due to the path v0->FD0->y (mediated by FD0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:dowhy.causal_identifier:If this is observed data (not from a randomized experiment), there might always be missing confounders. Causal effect cannot be identified perfectly.\n",
      "INFO:dowhy.causal_identifier:Continuing by ignoring these unobserved confounders because proceed_when_unidentifiable flag is True.\n",
      "INFO:dowhy.causal_identifier:Mediators for treatment and outcome:['FD0']\n",
      "INFO:dowhy.causal_identifier:All common causes are observed. Causal effect can be identified.\n",
      "INFO:dowhy.causal_identifier:All common causes are observed. Causal effect can be identified.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimand type: nonparametric-nde\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: mediation\n",
      "Estimand expression:\n",
      "Expectation(Derivative(y, [FD0])*Derivative([FD0], [v0]))\n",
      "Estimand assumption 1, Mediation: FD0 intercepts (blocks) all directed paths from v0 to y except the path {v0}→{y}.\n",
      "Estimand assumption 2, First-stage-unconfoundedness: If U→{v0} and U→{FD0} then P(FD0|v0,U) = P(FD0|v0)\n",
      "Estimand assumption 3, Second-stage-unconfoundedness: If U→{FD0} and U→{y} then P(y|FD0, v0, U) = P(y|FD0, v0)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Natural direct effect (nde)\n",
    "identified_estimand_nde = model.identify_effect(estimand_type=\"nonparametric-nde\", \n",
    "                                            proceed_when_unidentifiable=True)\n",
    "print(identified_estimand_nde)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:dowhy.causal_identifier:If this is observed data (not from a randomized experiment), there might always be missing confounders. Causal effect cannot be identified perfectly.\n",
      "INFO:dowhy.causal_identifier:Continuing by ignoring these unobserved confounders because proceed_when_unidentifiable flag is True.\n",
      "INFO:dowhy.causal_identifier:Mediators for treatment and outcome:['FD0']\n",
      "INFO:dowhy.causal_identifier:All common causes are observed. Causal effect can be identified.\n",
      "INFO:dowhy.causal_identifier:All common causes are observed. Causal effect can be identified.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimand type: nonparametric-nie\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: mediation\n",
      "Estimand expression:\n",
      "\n",
      "Estimand assumption 1, Mediation: FD0 intercepts (blocks) all directed paths from v0 to y except the path {v0}→{y}.\n",
      "Estimand assumption 2, First-stage-unconfoundedness: If U→{v0} and U→{FD0} then P(FD0|v0,U) = P(FD0|v0)\n",
      "Estimand assumption 3, Second-stage-unconfoundedness: If U→{FD0} and U→{y} then P(y|FD0, v0, U) = P(y|FD0, v0)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Natural indirect effect (nie)\n",
    "identified_estimand_nie = model.identify_effect(estimand_type=\"nonparametric-nie\", \n",
    "                                            proceed_when_unidentifiable=True)\n",
    "print(identified_estimand_nie)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Estimation of the effect\n",
    "Currently only two stage linear regression is supported for estimation. We plan to add a non-parametric Monte Carlo method soon as described in [Imai, Keele and Yamamoto (2010)](https://projecteuclid.org/euclid.ss/1280841733).\n",
    "\n",
    "The estimator converts the mediation effect estimation to a series of backdoor effect estimations. \n",
    "1. The first-stage model estimates the effect from treatment (v0) to the mediator (FD0).\n",
    "2. The second-stage model estimates the effect from mediator (FD0) to the outcome (Y).\n",
    "\n",
    "For estimating the natural indirect effect, there is also an additional second-stage model that estimates the effect of treatment on the outcome, conditioned on the mediator. It assumes the same model as given for for the `second_stage_model` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:dowhy.causal_estimator:INFO: Using Two Stage Regression Estimator\n",
      "INFO:dowhy.causal_estimator:b: FD0~v0+W0\n",
      "INFO:dowhy.causal_estimator:INFO: Using Linear Regression Estimator\n",
      "INFO:dowhy.causal_estimator:b: y~FD0+v0+W0\n",
      "INFO:dowhy.causal_estimator:INFO: Using Linear Regression Estimator\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-nde\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: mediation\n",
      "Estimand expression:\n",
      "Expectation(Derivative(y, [FD0])*Derivative([FD0], [v0]))\n",
      "Estimand assumption 1, Mediation: FD0 intercepts (blocks) all directed paths from v0 to y except the path {v0}→{y}.\n",
      "Estimand assumption 2, First-stage-unconfoundedness: If U→{v0} and U→{FD0} then P(FD0|v0,U) = P(FD0|v0)\n",
      "Estimand assumption 3, Second-stage-unconfoundedness: If U→{FD0} and U→{y} then P(y|FD0, v0, U) = P(y|FD0, v0)\n",
      "\n",
      "## Realized estimand\n",
      "(b: FD0~v0+W0) * (b: y~FD0+v0+W0)\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 11.704233518070275\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import dowhy.causal_estimators.linear_regression_estimator\n",
    "causal_estimate_nde = model.estimate_effect(identified_estimand_nde,\n",
    "                                        method_name=\"mediation.two_stage_regression\",\n",
    "                                       confidence_intervals=False,\n",
    "                                       test_significance=False,\n",
    "                                        method_params = {\n",
    "                                            'first_stage_model': dowhy.causal_estimators.linear_regression_estimator.LinearRegressionEstimator,\n",
    "                                            'second_stage_model': dowhy.causal_estimators.linear_regression_estimator.LinearRegressionEstimator\n",
    "                                        }\n",
    "                                       )\n",
    "print(causal_estimate_nde)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the value equals the true value of the natural direct effect (up to random noise). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.704233518070275 11.712215459759832\n"
     ]
    }
   ],
   "source": [
    "print(causal_estimate_nde.value, data[\"ate\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameter is called ate because in the simulated dataset, the indirect effect is set to be zero. \n",
    "Now let us check whether the indirect effect estimator returns the (correct) estimate of zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:dowhy.causal_estimator:INFO: Using Two Stage Regression Estimator\n",
      "INFO:dowhy.causal_estimator:b: FD0~v0+W0\n",
      "INFO:dowhy.causal_estimator:INFO: Using Linear Regression Estimator\n",
      "INFO:dowhy.causal_estimator:b: y~FD0+v0+W0\n",
      "INFO:dowhy.causal_estimator:INFO: Using Linear Regression Estimator\n",
      "INFO:dowhy.causal_estimator:b: y~v0+W0\n",
      "INFO:dowhy.causal_estimator:INFO: Using Linear Regression Estimator\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-nie\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: mediation\n",
      "Estimand expression:\n",
      "\n",
      "Estimand assumption 1, Mediation: FD0 intercepts (blocks) all directed paths from v0 to y except the path {v0}→{y}.\n",
      "Estimand assumption 2, First-stage-unconfoundedness: If U→{v0} and U→{FD0} then P(FD0|v0,U) = P(FD0|v0)\n",
      "Estimand assumption 3, Second-stage-unconfoundedness: If U→{FD0} and U→{y} then P(y|FD0, v0, U) = P(y|FD0, v0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+W0-(b: FD0~v0+W0) * (b: y~FD0+v0+W0)\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 0.000848600067884675\n",
      "\n"
     ]
    }
   ],
   "source": [
    "causal_estimate_nie = model.estimate_effect(identified_estimand_nie,\n",
    "                                        method_name=\"mediation.two_stage_regression\",\n",
    "                                       confidence_intervals=False,\n",
    "                                       test_significance=False,\n",
    "                                        method_params = {\n",
    "                                            'first_stage_model': dowhy.causal_estimators.linear_regression_estimator.LinearRegressionEstimator,\n",
    "                                            'second_stage_model': dowhy.causal_estimators.linear_regression_estimator.LinearRegressionEstimator\n",
    "                                        }\n",
    "                                       )\n",
    "print(causal_estimate_nie)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Refutations\n",
    "TODO"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}