{
 "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')"
   ]
  },
  {
   "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 -24.766460 -1.390853 -5.320208 -106.930248\n",
      "1  30.601065  2.605898  6.336803  133.087050\n",
      "2  -2.585961 -0.190542 -0.531693  -11.210317\n",
      "3   0.085982  0.094062 -0.227000    0.469403\n",
      "4  10.951420  0.637002  2.034312   47.314086\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": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAFbCAYAAABGcY7OAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd1wU1/4//tcuSy9LFVZAKYkNEYjGUMQGSjWWq2JNNEZNcvMREzX2bhKTG5NrEk1MotcUDUaDBUUBGwiIhdBEuiC9usDS2eX8/vC38wUpguzuAJ7n47EP2dmz57zXnZ33zJmZcziEEAKKoiiKeklx2Q6AoiiKothEEyFFURT1UqOJkKIoinqp8dgOgBqYmpqaUFtbCwAQi8UQiUQAgJaWFlRVVbUp29DQgPr6+ufW2d1yAKClpQVlZeUel+NwONDV1WWe6+rqgsPhAAB0dHSgpKTUrfYp6mVTV1eH/Px8CIVCiEQiSCQSVFdXAwD09PQAPP09GRoaQiAQQFVVlc1w26CJ8CVQU1PDPKqqqlBdXY3GxkaIRCI0Njairq4O9fX1aGhoQG1tLZqamlBdXQ2JRILKykq0tLSgsrKSSWitk5y0PACmjpeBpqYmVFRUAADq6upQU1MD8P+Spa6uLpSUlMDn88Hj8aCtrQ0VFRVoamoy5aV1aGtrg8fjQVdXF6qqqtDS0oKuri60tbWhpaUFdXV1Nj8qRbWRn5+PO3fuIDk5GQ8ePEBKSgpyc3OZpNddhoaGGDJkCGxsbDB69GiMHj0ajo6O0NfXl1PknePQq0b7LrFYDKFQ2OGjuroalZWVEIlETJKrrq5GVVVVm8QnFAq7bENZWRlaWlpQVVWFhoYGs5GWHil1tWEHADU1NWZDLd3QA2hThsvlgs/nM21K9w6lnj0K60x3y3V01Nndcq33YgkhqKysZF6rrKyE9Oci3VEAnu5oNDc3M2U623GQHtFKdx5EIhHEYjGqqqrQ0tLSaZxKSkrQ0dEBn8+HlpYW89DT02P+1tbWho6ODnR1daGnp9fhQ3pkS1E9UVlZieDgYFy9ehURERHIysoCl8uFlZUVbG1tMWrUKAwdOhQCgQBmZmYwMDCAlpYWs94CYLZDVVVVKC0tRXFxMfLy8vD48WMkJSUhOTkZeXl54HK5sLGxweTJkzF9+nS4u7szO5nyRBOhggiFQpSUlKC8vBxlZWUoLi7GkydPOk100u6FZ/F4POjp6TEbPelRg6amJvh8PnR0dNpsHHV1ddtsPKXvUVFRaZOcKHZJk670qFx65C7doRGJRKisrGyzk/Ps86qqKlRWVkIoFKKjn3VXSVJPTw+DBg2CoaEhjIyMMGjQIBgbG0NLS4uF/w2KbTU1NTh16hROnz6NGzdugBACJycnTJ48GRMnToSjoyOz0ysrT548QWRkJMLDwxEREYF//vkHmpqa8PLywoIFC+Dr69ut0x0vgibCFySRSFBSUoL8/HwmwRUXF6OsrIx5SJ+Xl5cz3YdS+vr6MDQ07HLD1NFDepRFUV2RJsTuPKQ7ZGVlZe12vtTV1WFoaAgTExMYGRkxj9bPBQIBBAIBjI2NWfq0lKwkJyfj8OHD+OOPP9DU1AQvLy/MmTMHvr6+3eqNkaXCwkKcO3cOZ8+exfXr12FsbIwVK1Zg9erVMDMzk2lbNBF2oL6+HkVFRSgsLGzz76NHj5i/c3NzIRaLmfeoqakxyWrw4MEQCATt/pY+NzMzY84vUVRf0tjYiIqKCiZJStf/Z/8WCoUoKCho07WsoqICAwMDZp3v6F9ra2uFb1Cp53vw4AG+/PJLnDx5ElZWVlixYgVWrFgBQ0NDtkMD8DQp/v777zh06BCKi4uxYMEC7Ny5E9bW1jKp/6VMhCUlJcjOzkZOTk6bfx8/foz8/HzU1dUxZdXV1WFqasr0fwsEApibm0MgEMDU1BSmpqYwNjaGhoYGi5+IotghEolQUFCA4uJi5Ofno7CwEAUFBSgoKEBRURHy8vJQUlLSpkdER0cH5ubmsLS0hIWFBSwsLJi/LS0t251DpuSnoKAA69atw19//QUHBwfs2LEDb775Zp89n9zU1ITjx4/j888/R2FhIT788EPs2rWr1z1lAzIRNjY2Ii0tDRkZGe2SXXZ2NnMJPo/Hg5mZWZsfpDTJmZmZYfDgwaxcwURRAwkhBCUlJSgqKmKSZH5+fpvfZVFREXNek8/nt0uSVlZWGDZsGKysrMDj0Yvde0sikeDbb7/Fzp07MWjQIHz11VeYOXNmn02Az2pubsbRo0exbds2qKqq4ptvvsH8+fNfuL5+nQiFQiGSk5Px8OFDPHr0CI8ePUJycjLS0tKYK/r09PRgZWXV4WPIkCH0R0VRfUBTUxPy8/OZ33Hr0xGPHj1CdnY2CCHg8XgYMmQIrKysMGrUKNjY2MDKygo2NjYQCARsf4x+IT8/H0uWLEFMTAw2btyITZs29dtbdMrLy7Fp0yYcO3YMixcvxuHDh1/o6LBfJMKSkhLEx8cjPj4eSUlJSEtLQ1paGnNiX1dXF8OGDcOIESMwfPhw5vHqq6/2qZs2KYp6MTU1NUhLS0N6ejpSU1OZbUB6ejpzKsPAwAAjRozAiBEjYGtrC3t7e9jZ2dFzkq2EhIRg8eLFMDIyQkBAAOzs7NgOSSYuX76MZcuWgc/nIzAwEKNHj+7R+/tUIpRIJMjIyEBCQgLi4uKQkJCA+Ph4FBcXAwDMzc0xZswYJuENGzYMI0eOxKBBg1iOnKIoNhBCkJubi/T0dKSlpSE1NRWpqalISEhAeXk5AMDS0pJJitJ/LSws2A2cBcePH8fKlSuxcOFC/PDDDzK//YFtRUVFWLBgARITE3Hu3DlMmjSp2+9lNREWFhYiKioK0dHRiImJQWJiIurq6sDj8TBy5Mg2K6+9vT0MDAzYCpWiqH4mPz+f2ZmWPrKyskAIgZ6eHuzt7eHs7AxnZ2c4OTkN6It0Dhw4gA0bNmDr1q3Ys2dPvzkX2FONjY1YunQpLly4gNOnT2PGjBndep/CEqFEIkFycjIiIyMRHR2NqKgo5OTkQElJCWPGjIGzszMcHBxgb2+P0aNH0y5NiqJkTiQSITExEfHx8YiNjcXt27eRlpYGDoeDkSNHwsXFhUmOr776KtvhysSxY8fw7rvv4ptvvoG/vz/b4chdS0sL3nvvPfzxxx8ICQmBq6vrc98j10T46NEjBAcH4/Lly4iMjER1dTV0dHTg5OQEZ2dnuLi44I033qCjV1AUxZry8nLcvn0bUVFRiIqKwv3799HQ0ABjY2NMmjQJ3t7e8PLy6penYEJCQuDr64stW7Zg9+7dbIejMBKJBH5+frh27Rru3LmDYcOGdVlepomwqakJt27dQnBwMIKDg5Gamgo+n49p06Zh6tSpcHFxgY2NDR3Bn6KoPqupqQmxsbGIjo5GWFgYwsPD0dTUhLFjx8Lb2xve3t4YN24cuNy+PYtdcXEx7OzsMH36dPz+++9sh6NwjY2NcHV1hUQiwe3bt7scxKTXibC+vh5BQUE4deoUwsLCIBKJMGrUKPj4+MDLywsTJkyQ2/hwFEVR8lZbW4vr168zO/i5ubkYNGgQfHx84OfnB3d39z65c+/t7Y2MjAz8888/L+3QjFlZWXBwcMB7772HL7/8svOC5AXdunWLvPXWW0RbW5soKSmR6dOnk0OHDpHs7OwXrVIu8vLyCIB2j7Nnz7Ypt3Xr1nZlUlJS5Brbf/7zH6YtU1NTubbVn/z555/M/4uqqmqP33/37l3y9ttvEwsLC6Kmpkb09PSIjY0NmTNnDjl8+DDJzMyUQ9SdCwgIIHZ2dkRNTY35XElJSQqNoad6+x0MZImJieSLL74gb7zxBgFAjI2Nib+/P3n48CHboTEuXLhAOBwOiYyM7HVdxsbGz902/utf/yIASH5+fpvlz25XP//8c+a1uLg44u3tTfh8PtHS0iJubm4yifdZR44cIcrKyiQ9Pb3TMj1KhA0NDeTIkSPE1taWACBjx44lBw8eJMXFxb0OVt6kP+yNGzd2WW7SpEnk559/VlBUT9nZ2dFE2AE3N7cebYQlEglZv3494fF4ZMOGDSQlJYU0NDSQ4uJiEhoaStzd3ZkfZHNzsxwj/38iIyMJh8MhGzZsICKRiGRmZhIzM7M+nwilevodvGwyMjLI7t27iZWVFQFAJk+eTAIDA0lLSwtrMTU3N5ORI0eSefPmyazOn376iQAga9eubfeaWCwmurq6BAA5evRou9crKiqIrq4uaWpqYpbFxMQQdXV14ufnRwoLC0lZWRlZuXIl4fF4JCQkRGZxS+OzsbEh//rXvzot061E2NTURL777jtiZmZGVFVVyfLly8ndu3dlFqgi0ETY//R0I7xlyxYCgPz0008dvi4Wi4mXl5dCE6G/v3+He8r9BU2E3SORSMjly5fJjBkzCJfLJba2tuSvv/5iJSGeP3+ecLncLo+Aeio3N5cAICNGjGj3WlRUFNHQ0CAAOky+AQEBZObMmcxziURCbGxsiEAgIHV1dcxysVhMhg8fTszNzUlDQ4PMYieEkMDAQMLhcEhWVlaHrz/3bG94eDjs7e3xySefYM6cOcjKysKxY8fw+uuvd6uPlqIUITU1Ffv378fYsWOxcuXKDssoKSlh+/btCo0rLy8PAOg9sAMcl8uFp6cnLly4gISEBNjY2GDBggVwc3NDSkqKQmM5duwYpk6dKtPbP8zNzTFy5EikpqYiNze3zWshISF47733oKWlhatXrzLDW7Z+3cPDg3keERGB5ORkzJ07t83QbkpKSli4cCHy8vJw8eJFmcUOADNmzIBAIMDx48c7fL3TREgIwRdffIGpU6fCysoKycnJOHjwIExNTWUaIEXJwk8//YSWlhbMmzevy3JOTk7MmJWK8OxGgRr4Ro8ejT///BP37t1DXV0dHBwc8PPPPyukbZFIhODgYLz99tsyr9vT0xMAcOXKlTbLr1y5ghkzZmDq1KkQCoW4e/dum9dDQ0OZ9wLA9evXAQDjxo1r14Z02bVr12QaO4/Hw9KlS/Hnn392+HqHiVAikWD+/PnYtWsXfvrpJwQFBcHS0lKmgfUX586dA4fDYR45OTnw8/ODrq4uDAwM4Ovri6ysrHbvq6iowMcffwxra2uoqKhAT08PXl5euHHjRqdtpaamwsfHB3w+HxoaGpgyZQqioqLalGlsbMSOHTswYsQIaGhoQF9fHzNmzMCFCxfabXTLysqwZs0aWFhYQEVFBUZGRpgzZw7i4+M7/XxpaWmYP38+DAwM2ixv/di3bx8AQCwWt1k+d+7cHrXd+nPPmjULfD4fmpqacHV1RWRkZPe+oP9fREQEAGDMmDE9eh/Qve+qp+uBtPz58+cBPJ3Oi8PhwNHRsUft7tu3j2lzwoQJzPIrV64wy1vPGfei62tPvwNZrVvl5eU9Wqf7k9deew1RUVFYt24dVq9ejZ07d8q9zZiYGDQ3N8PNzU3mdUuP6kJCQphlT548QWpqKpydnTt8/cGDB9DQ0GiTP1JTUwGgw8l1pQda6enpMo/fzc0NmZmZKCoqav9iR/2lGzZsIOrq6uTWrVsy7adlU2/PEc6cOZMAIDNnziTR0dGkpqaGhIWFEXV1dfL666+3KVtUVEQsLS2JsbExCQoKIlVVVSQtLY3MmTOHcDicdvXb2dkRPp9PpkyZQiIjI4lIJCL37t0jY8aMISoqKuTmzZtM2XfffZfw+XwSGhpK6urqSHFxMVm/fj0BQG7cuMGUKywsJEOHDiXGxsbk0qVLRCQSkQcPHpBJkyYRNTU1Eh0d3eHnmzRpErlx4wapra0lMTExRElJiZSVlRFPT0/C5XI7vOLSycmJnDx58oXazsjIILq6usTU1JSEhoYSkUhEEhMTyfTp04mFhUW3z08JBAICgNy5c6db5aV6+l31ZD1oXb6+vr5X7WpqahIXF5d29Y8dO5YYGBh02m534uzpdyDrdau763R/duzYMcLhcMivv/4q13Z2795NLC0t5VJ3fX09UVdXJ3w+nznHHhAQQGbMmEEIISQzM5MAIG+88Qbznq+++op8+OGHbeqZNm0aAUBiYmLatZGRkUEAkNdee03m8VdXVxMlJSXy999/t3utXSJMT08nysrKHV7905/JKhEGBQW1WT537lwCgJSVlTHLli1bRgCQP//8s03ZhoYGMnjwYKKurt7mSls7OzsCgNy+fbtN+cTERAKA2NnZMcssLS2Js7Nzu/iGDRvWZqPx9ttvEwDkxIkTbcoVFRURVVVVMnbs2A4/X3BwcLu6CSHk6tWrBAD54IMP2iyPjIwkQ4YMaXPxSU/anjdvHgFAzpw506ZsQUEBUVVV7XEi7OlFXD39rnqyHrQu/2wi7Gm7L5oIuxNnT78DWa9b3V2n+7uPPvqIGBsbk5qaGrm18c477xBPT0+51e/h4UEAMAdJy5cvJ99//z3zurW1NeFyuaSiooIQ8jTpXbx4sU0dXSXC9PR05o4EeTA1NSVff/11u+XtukYvXboEQ0NDLFu2rFeHoX2N9IbX53W1SCSSLm+OffYiIXNzcwBPBxCXOnv2LADAx8enTVlVVVW4ubmhvr6+TfcBAKipqeGNN95os8zW1haDBw9GQkICczjv6emJ6OhorFq1CjExMcznSUtLw+TJk5n3njt3DlwuF76+vm3qNDExgY2NDWJjY5Gfn9/u840fP77Dz+3m5gYHBwccP34cFRUVzPL//Oc/WLt2bZtzbj1pW3q+ofXJdAAYPHjwc4dFerY8AGbGge56ke8K6N56II92e6o7cfb0O5D1utXddbq/27RpE0pKShAdHS23NioqKuR6YZZ0HZGuM8+e//P09ERLSwvCwsJQX1+Pu3fvtvsOpdNi1dbWtqtfukxeU2cZGhq22X5JtUuERUVFGDx4cJ8fPqinpOOZVldXd1musrISOjo6nb7O5/PbPJcO29PS0gLg6Tm8qqoqqKmpdTiag7GxMQAwU0tJSc+bPEs6vmFpaSkA4NChQ/jtt9/w6NEjuLm5QUdHB56ensyGtXUMLS0t4PP57c7x/fPPPwCAjIyMdu11NTXLunXrUFdXh8OHDwN42o8fERGBd99994XabmxshEgkgpqaWofjzfZkbEfplCuJiYndfs+LflfA89cDebXbU91ZX3vyHchj3erOOj0QGBkZQUNDo+NzVDJSV1cn10l2pUkvJCQESUlJUFNTg7W1NfN66/OE4eHhGDduXLvvfcSIEQDQ4c5SQUEBAPRoJ7gnNDU1UVNT0255u2xnY2ODlJSUHu9Z93XS/9jk5OROyzQ2NiIzM7NXlx2rqqqCz+ejoaGBmTi4tZKSEgBP955bq6qq6rA+aQKUbpA4HA6WLl2Kq1evorKyEufOnQMhBHPmzMHXX3/NxKCrqwsej4fm5maQp13g7R5Tpkzp0Wfz8/ODubk5vv/+ezQ2NuLAgQNYuXJlm415T9pWVVWFtrY2GhoaOlw5nzx50u3YVq9eDR6PhzNnznRZ7pNPPgGXy0VqauoLf1e99SLtcrlcNDU1tStbWVnZ61h68h3IY93qzjo9ENy/fx91dXU9njS2J/T09CAUCuVW/8iRI2Fubo7Y2Fj88ccf7XoRpkyZAhUVFYSEhODKlSvtXpeWAYDY2Nh2r0mXyeNiH+Dp+qyvr99uebtEOHfuXOjq6mLdunVyCYQt1tbWGDFiBGJiYjrcWwWAv/76C0ZGRr1eUWfPng3gaTdza42Njbh27RrU1dXbrSA1NTVISEhosywpKQmFhYWws7ODQCAA8LTLQHrVlbKyMqZNm8Zcnde6vTlz5kAsFre76hQAvvjiCwwZMgRisbhHn4vH48Hf3x+lpaU4cOAAAgICsGbNmnbletK2l5cXgPaXZJeXlyMtLa3bsQ0bNgw7d+7E/fv3cezYsQ7LpKWl4ciRI5g/fz6zV/oi35Us9LRdgUDA7C1LFRcXt7un60X09DuQ9brV3XW6P2toaMDatWsxYcIEvPbaa3Jrx9DQEGVlZXKrH3h61EcIwbffftumWxR42vPm4uKCwsJCHD9+vMPfzqRJkzBq1CicOXMGDQ0NzHKJRIKAgACYm5u3O2UgK2VlZW2usmZ0dELx0qVLhMfjkY8//piIxWJZnqtk1eXLl4mysjKxtrYmf//9N6moqCBisZgUFBSQQ4cOER0dHXL69OkO39vZRQ8bN24kAEhcXByz7NkrAqurq9tcEfjsyCd2dnZEU1OTTJgwgcTExJCamppOrxrl8/lk0qRJJCEhgTQ0NJCSkhKya9cuAoDs27ePKVdSUkKsra2JlZUVCQ4OJpWVlaSiooL8+OOPRENDg5w6dapbn+9Z1dXVhM/nEw6HQ956660Oy/Sk7czMTKKvr9/misXk5GTi4eFBBg0a1ONRTTZt2kSUlZXJxo0bSVpaGmlsbCT5+fnkl19+IQKBgEyYMKHNxQo9/a56sh50Vb6n7X744YcEAPnuu++Yodrmz59PTE1Nu7xYpjtx9vQ7kPW61d11ur+qqakhM2bMIHp6eiQ5OVmubUm3Y/Lcbp85c4YAIMrKykQkErV7ff/+/QQAEQgEndZx+/ZtoqamRhYsWECKiopIeXk5Wb16NeHxeOTKlStyiTsnJ4cAIOHh4e1e63SItZMnTxI1NTUyffp0UlBQIJfA2BAbG0uWLFnCXBauoqJCzMzMyLx580hUVFS78rdv32434OzWrVsJIaTdch8fH+Z95eXlZO3atcTS0pIoKysTPp9PPDw8yLVr15gyzw66fffuXTJlyhSipaVF1NXVyaRJk9oNQhsfH09Wr15NRo4cSTQ0NIi+vj5xdHQkP//8c7vhnCoqKsjHH39MrKysiLKyMjEyMiLTp08nYWFhXX6+TvaPGBs2bCAASEJCQqdlutO2VFpaGpk1axbR0dFhLu+/ePEicXNzY+JZsWJFlzG1dvfuXbJ06VJibm5OlJWViba2NnF0dCQHDx4kjY2N7cp357vq6Xpw9uzZDv9fW18Z3J12pSorK8m7775LBAIBUVdXJxMmTCD37t0jY8eOZereuHHjC6+vPf0OZLlu9WSd7m8SExPJ6NGjiZGRUburwuUhLi6OAPId1L2yspLweDwyefLkLmNYtmxZl/X8888/xMvLi+jo6BAtLS0ydepUuQy6LXXy5EmirKxMamtr273W5TRMsbGx8PPzQ2lpKXbu3IkPP/yQzhxPURT1HEKhEPv27cN3332HcePG4dSpU8wVu/IkkUggEAiwZs0abNu2Te7t9Sd+fn4oLi5GeHh4u9e6vDR07NixSEpKwkcffYRt27bB2toa3377Lerq6uQWLEVRVH9VXl6Obdu2wcLCAr/++iu+//57REZGKiQJAk9vE1uyZAmOHTvWrSuYXxYVFRU4f/48li9f3nGB7h5W5ufnE39/f6Kurk709PTIRx99RNLS0mR44EpRFNU/3b59myxdupSoqakRQ0ND8umnn5KqqipWYnnw4AHhcDjk3LlzrLTfF+3Zs4fw+fxOBzPo8Qz1paWlOHbsGI4cOYLHjx/D0dERixcvxvz582FkZCSL5E1RFNXnZWdn4+TJkzhx4gRSUlLg4OCA999/H4sWLeryflxF8PPzQ0JCApKSkqCsrMxqLGwrKyvDK6+8gk8++QRbt27tsEyPE6GUdPSAEydO4OzZs2hoaICLiwu8vLzg7e0NW1vbXgVPURTVl0gkEty5cwfBwcG4fPky4uLiYGRkBD8/PyxdurRPTU2XlZWFUaNG4fPPP8fHH3/MdjiseueddxASEoKMjAxoaGh0WOaFE2FrdXV1uHjxIi5duoQrV66gtLQU5ubm8Pb2hre3N9zc3FjfQ6Ioiuqp8vJyhISE4NKlSwgNDUVFRQUsLS3h7e2NGTNmwM3NTWFTevXUp59+ir179+L27dtwcHBgOxxWnDlzBvPnz8fZs2cxc+bMTsvJJBG21tLSgtjYWAQHByM4OBj379+HsrIyxo8fDxcXF7i4uMDJyYlOVEpRVJ+Tn5+PqKgoREdHIyoqCvHx8VBSUoKrqyu8vLzg4+PDDMbQ17W0tMDNzQ1FRUWIiop66ba5Dx8+hIuLCxYvXozvv/++y7IyT4TPKi0tRUhICCIiIhAdHc3M1jxixAg4OztjwoQJcHJywvDhw+UZBkVRVBsSiQRJSUlM4ouMjERubi54PB7s7e3h4uKCyZMnw83NrcMxYfuDoqIiODs7w8TEBFevXn1peuby8vLg4uKCoUOHIiwsDGpqal2Wl3sifNaTJ08QHR3N7HHdu3cP9fX1MDQ0hIODA+zt7WFnZwd7e3sMHz68z3Y7UBTVfzQ0NCA5ORlxcXFISEhAQkIC4uPjIRKJwOfz4ezsDGdnZ7i4uGD8+PEDKmGkp6djwoQJGDNmDM6ePdtvk3p3PX78GB4eHlBRUUF4eDj09PSe+x6FJ8JnNTc3IzY2Fnfu3EF8fDzi4+Px8OFDNDU1QU1NDaNHj26THMeMGdPl7BAURb3cysrKkJCQ0CbppaamQiwWQ1NTE7a2trC3t4eDgwOcnJxgY2Mz4GbbeVZ8fDy8vb1hYmKCS5cuMWMXDzQJCQnw9vaGoaEhrly50u3PyXoi7IhYLEZaWhoePnyI5ORkxMbG4u7du8xMDHp6ehg1ahRsbGxgZWUFKysrjBo1CiNGjOhyLkGKogYGsViM3NxcPHr0CMnJyXj48CHzt3SaI+l2YuzYsczjZd5G5OTkwNPTE7W1tThx4gQmTpzIdkgydeLECbz//vt4/fXXERgY2G4Ksq70yUTYmcePH+PBgwdITU1FWloa0tLSkJqayiRINTU1DBs2DMOHD2f+tbKygoWFBQYPHtzhfH8URfVNYrEYeXl5yMnJwaNHj5Ceno60tDSkpKQgOzsbzc3N4HA4MDc3Z37zI0eOxPDhwzFmzJgezWf5snjy5AlWrFiBoKAgbN++HVu2bOn39xlWV1fD398fv/76K/z9/fHFF18w8252V79KhJ0RCoVIT09HSkoK0tLSkJ6ejtTUVGRmZjJzuKmqqmLo0KGwsLBgHpaWlszfsp5zjqKorkkkEhQWFiInJwfZ2dnIzs5GTk4O88jPz2emc9LU1NAKZ2AAACAASURBVGR2bocPH44RI0YwyW8gnc9TBEIIDh06hA0bNuCVV17BoUOH+u3RYUBAANatW4empiYcO3YMM2bMeKF6BkQi7ExLS0uHPzTpv61/aOrq6rCwsICZmRkGDx4MMzMzmJiYYMiQITAxMYGZmRmMjY1f2m4ViuqJpqYmFBcXIz8/H4WFhSgoKEBBQQGKioqQm5uLgoIC5Obmorm5GcD/21FtvXPaemfV2NiY5U808GRlZeH//u//cOXKFfj5+WHHjh0YOXIk22F1S0REBHbs2IFbt25hxYoV+Pzzz3t1e8iAToTPIxaLkZ+f32ZPND8/H0VFRcjLy0NRUREqKiqY8kpKSjA2Noa5uTlMTExgbm4OgUAAgUAAIyMjGBkZwdjYGEZGRnQvlRqQqqurUVxcjLKyMpSXl6O0tBSFhYXMQ/r7KSkpYd7D4XBgbGyMwYMHw9TUlHlIE52lpSUEAgE9daFgeXl5WLhwIUQiEVpaWvDw4UP4+flh/fr1cp08+EURQhAWFob9+/fjxo0bmDJlCvbv34/x48f3uu6XOhF2R0NDA7Mn+2ySlO7tlpSUoLa2ts37NDQ0YGRkBBMTExgZGcHQ0BCDBg1iEqWhoSGMjIygp6cHPT096Ovrs/QJqZeVWCyGUChkHtLkVlJSgtLSUpSVlaGsrAwlJSXM342NjW3q0NHRweDBgyEQCJjeFGmPikAggLm5OYyNjfv9eaiB5tq1a1i8eDF0dXVx+vRp2NjYIDAwEHv37kViYiLGjRuH1atXY968eT266EQeiouL8dtvv+Hnn39GZmYmpk6dip07d8q0O5cmQhmpq6tDWVkZiouLUV5ejrKyMpSWlqKkpIR5Lt3AlJeXo6GhoV0drZOi9O/Onuvo6IDP50NLS4t5UC+fqqoq1NTUMI8nT560SW6tnz/7WnV1dbv6NDU1mR046c6adAdO+tzY2BiDBg2CoaEhnZ+0nxGLxdi3bx/27t2LRYsW4Ycffmi37YiMjMSPP/6IM2fOAADc3d0xZ84c5vYLRcjKykJQUBACAwMRFRUFHR0dvPXWW1i1ahVsbGxk3h5NhCwRiUQoLy/vcsPV0bKqqqpO69TV1W2TGJ99zufzoa2tDS0tLairq0NbWxs8Hg98Ph9KSkrQ1dWFkpISdHR0oKKiAk1NTaipqUFdXV2B/zMDT01NDZqbm1FdXQ2JRAKhUIiWlhZUVVWhubkZNTU1aGhoQH19fZukVlVVBZFIxDyvrq5GdXU181wkEnXYHofD6fZOVetlhoaGnQ5KTPV/BQUFWLhwIe7du4f9+/fD39+/y/JCoRAXLlxAYGAgQkND0dDQgOHDh8PV1RXOzs6wtbXFqFGjer3OVFZW4sGDB0hMTERUVBTCw8NRUFAAXV1dzJgxA3PmzIGHh4dct0M0EfYzLS0tePLkCaqrq9sdDQiFQubv2tpaVFZWttuQSt/T2NiIqqqqbk/eqampCRUVFejo6DBJU3pOp/XIDXw+n7k5WZpogaddxdKjh86Sa+vynVFVVX3uD08kEjEXQXWmsbGxwwmmpUkLeHqUL+0KlP5/KSsrM8kLeHrlo/TISiwWQyQSoampCbW1taivr+/wyL8jysrK0NLSgoaGBrPjoqen125HRkdHB5qamtDS0uqwV0BfXx+6urrdapN6eVy/fh2LFy8Gn8/HX3/9hTFjxvTo/bW1tYiKisLhw4dx7949CIVC1NfXg8vlwtLSEkOHDoWpqSkEAgGzYwU83TZIfyMSiQRVVVUoLS1lLqTKyclBXl4egKfbDicnJ7i6umLixIl44403FNalThMhBaFQyKyszx6hSJOBNLlUVlaipaUFlZWVANomAmldUq0TbevkVFtby9zWIvVsPZ3pTpJTVVWFsrIyCCFd3k/U0dBL6urqzLiErZNuc3Mzbt++DTs7O1hbW7cZ3UhaD5fLBZ/PB4/Hg7a2NvN+6U6AlpYWlJWVmZ0FPT095gicouRBIpFg79692Lt3L2bPno2jR4++8Dm/pqYmjBw5ElOnTsWPP/6IR48eISkpCQ8fPkR+fj4KCgpQWFiIyspKlJWVoaWlBbW1tczvQrquDxo0iLkS39zcHKNHj4aNjQ2GDh0q40/fA7KbA5ii+o4NGzaQ4cOHy6y+lpYWsnPnTgKArFmzhkgkEpnVTVHyUFJSQqZNm0bU1NTIf//7317Xd/DgQaKmpkZyc3OfW/aDDz4gEydO7HWbijKwB9ijXlqenp5IS0vDo0ePZFIfh8PBrl278L///Q8//PADFixYgPr6epnUTVGyduPGDdjb2+Px48eIiYl57vnA56mtrcVnn32GDz74AObm5s8tX1BQAFNT0161qUg0EVID0oQJE6CtrY2QkBCZ1rts2TJcvnwZYWFhcHNzQ1lZmUzrp6jeaGlpwRdffIFp06bByckJd+7cgZ2dXa/r/eabb1BbW4uNGzd2qzxNhBTVB6ioqGDKlCkyT4QA4ObmhsjISBQWFsLJyQlpaWkyb4Oieqq0tBReXl7YuXMnDhw4gL///lsmF04JhUJ8/fXX2LBhQ7fHby0sLMTgwYN73bai0ERIDVgeHh64evVqu5vAZcHGxgYxMTHQ09ODs7Mzbt26JfM2KKq7bt68CXt7e6SnpyMiIqLXXaGtff7551BSUsLatWu7VV4ikaCkpIQmQorqC7y8vFBbW4vo6Gi51G9iYoKbN2/C1dUV06ZNw8mTJ+XSDkV1hhCCL774Au7u7nB0dERcXJxMhhyTKiwsxKFDh7B169ZuX91cXFwMiURCu0Ypqi+wtLTEsGHD5NI9KqWpqYm///4bq1evxpIlS7Br1y65tUVRrZWVlbXpCg0MDJT5PaS7d++Gnp4eVq9e3e33FBYWAgA9IqSovsLT0xNXrlyRaxtKSko4ePAgfvzxR3z66adYsWIFc1M+RclDeHg47O3tkZqaivDwcJl2hUplZGTgf//7H/bs2dOjUV0KCgrA4XBoIqSovsLDwwOJiYnMXqo8rVq1CkFBQTh9+jS8vb27HA6Pol4EIQQHDx7EtGnT8PrrryMuLg5vvPGGXNrasWMHLC0t8dZbb/XofQUFBdDX12cGpugPaCKkBrQpU6ZATU1Nrt2jrXl6euLWrVtITU3FhAkTkJubq5B2qYGvvLwc3t7eWL9+Pfbu3YuzZ892ODqSLCQlJeGvv/7CZ5999txhD59VWFjYr84PAjQRUgOcuro6XF1dFZYIAcDOzg4xMTHg8XhwdHREbGyswtqmBqaIiAjY2dnh4cOHuHXrFjZu3CjX+Rs3bdoEBwcHzJkzp8fv7W+3TgA0EVIvAQ8PD4SGhkIikSisTVNTU0RERMDe3h6TJk1CUFCQwtqmBg5pV6i7uzvGjRuH+Ph4ODo6yrXNyMhIBAcH4/PPP3+hZNvfbqYHaCKkXgKenp4QCoW4e/euQtvV1tbGhQsXsHTpUsyePRvff/+9Qtun+rfy8nL4+PgwXaHnzp2TW1doa9u2bcPEiRMxbdq0F3p/QUFBvzsi7FnnL0X1Q6NGjcLQoUMREhICJycnhbbN4/Hwww8/YMSIEfD390dGRga++eYbZqoqiurI3bt34efnB4lEgoiICIWtt5cuXUJ4eHiv7r2lXaMU1Ud5eHjI/TaKrvj7++PUqVP4+eefMXfu3A7nQqQoaVfohAkTYGtri/j4eIUlwZaWFmzfvh2zZs164Tbr6upQWVlJu0Ypqi/y8PDAvXv3WB0ke+7cubh27RoiIyMxZcoUlJSUsBYL1fdUVFRgxowZWL9+PbZs2YJz585BX19fYe0HBAQgISGhV4NCFBQUAABNhBTVF7m7u0NJSQlhYWGsxuHk5ITbt2+jqqoKTk5OSElJYTUeqm+4d+8eXn/9dSQmJuLmzZvYtWuXQrvPm5ubsXPnTixZsqRXs1VI79cVCASyCk0haCKkXgo6OjpwdXXFpUuX2A4F1tbWiI6Ohrm5OVxcXHDjxg22Q6JY0ror1MbGBnFxcXBxcVF4HEePHkVubi527tzZq3oeP34MVVVVGBsbyygyxaCJkHpp+Pj44PLlyxCLxWyHAn19fYSGhsLLywvTp0/HTz/9xHZIlIJVVVVh3rx5WL9+PTZv3ozz58/DwMBA4XE0NDTg008/xerVq2FlZdWrurKzszF06NB+dzFY/4qWonph5syZEAqFcpuNoqdUVVXxxx9/YOvWrVi9ejX8/f1BCGE7LEoB7t+/DwcHB9y5cwc3btxQeFdoawcPHoRQKMSWLVt6Xdfjx49hYWHR+6AUjCZC6qVhbW2N4cOH4+LFi2yHwuBwONi1axeOHj2KH374AQsWLEBDQwPbYVFyIu0KdXFxgZWVFe7fv48JEyawFk9VVRW+/PJLrF27FiYmJr2uLycnhyZCiurrfH19++QoL++88w6Cg4MREhICNzc3lJeXsx0SJWPV1dXw8/NjukJDQ0NZP5f25ZdfoqWlBevWrZNJfTQRUlQ/4Ovri9TUVGRkZLAdSjvu7u6IjIxEfn4+nJyckJ6eznZIlIzExsbCwcEB0dHRuH79OqtdoVJlZWX47rvvsHnzZpmMWCORSJCfn08TIUX1da6urjAwMOgTV492ZPTo0YiJiQGfz4ezszMiIyPZDonqpZ9++gnOzs6wsLDA/fv34erqynZIAIA9e/ZAS0sLH374oUzqy8/PR3NzM02EFNXXKSkpYfr06X3qPOGzBAIBwsPD4ezsDHd3dwQEBLAdEvUCpF2hH3zwATZv3oywsDCZnIeThZycHPz888/YuXMnNDQ0ZFJndnY2ANBESFH9gY+PDyIiIlBZWcl2KJ3S1NTE2bNnsWrVKixatKhXo31QivfPP//gtddew82bN3H58uU+0RXa2o4dOzB48GAsX75cZnXm5ORATU2tzyT7nug73wxFKYiXlxcIIayPMvM8SkpK+Pbbb/HNN99g7969WLlyJZqbm9kOi3oOaVfokCFDEB8f/8KzOMhLcnIyTp48iU8//RQqKioyq1d6oYw850mUF5oIqZeOvr4+nJ2d+3T3aGv+/v74+++/cfLkSfj4+KC6uprtkKgOiEQiLFiwAO+//z7Wrl2LsLCwPjnU2NatWzFq1Cj4+fnJtN7+eg8hQBMh9ZLy9fXFpUuXFDpZb2/MmjULN27cQGJiIiZMmIDc3Fy2Q6JaiYuLw2uvvYYbN27g8uXL2L9/P5SUlNgOq5179+7hwoUL2L9/v8y7arOzs2kipKj+xNfXFxUVFYiJiWE7lG4bP348YmJiIJFI4OTkhH/++YftkCgAv/32G1xcXGBmZob4+HhMnz6d7ZA6tWnTJjg7O8Pb21vmdffXewgBmgipl9TIkSPx6quv9pvuUSkLCwvExMRgzJgxmDRpUr+LfyARiURYtGgRli1bhjVr1uDq1at9sitUKjQ0FNevX8e+fftkXrdYLEZBQQFNhBTV3/j4+PTLRKKtrY2goCAsWrQIs2bNwqFDh9gO6aXz8OFDODo64urVq326K1SKEIIdO3bAx8cHkydPlnn9eXl5EIvFNBFSVH/j6+uLBw8eMPc/9Sc8Hg9HjhzBgQMH8H//93/w9/dHS0sL22G9FH777Te8/vrrMDIyQnx8PDw8PNgO6bn+/vtv3L17F3v37pVL/Tk5OQD65z2EAE2E1Ets4sSJ0NXV7bOjzHSHv78/Tp06hZ9++gnz5s1DXV0d2yENWDU1NVi8eDGWLVuGd999F2FhYRg8eDDbYT2XRCLBjh07sGDBAjg4OMiljZycHKirq2PQoEFyqV/eaCKkXlrKysqYNm1av+webW3evHm4du0aIiIiMHXqVJSUlLAd0oCTkpICR0dHhIaGIjg4GAcPHoSysjLbYXXL8ePHkZGR0etJd7uSk5MDS0vLfnkPIUATIfWS8/X1xc2bNyESidgOpVecnZ1x+/ZtCIVCODk5ISUlhe2QBozffvsN48aNg4GBAeLj4+Hp6cl2SN3W0NCAPXv2YMWKFRg+fLjc2unPV4wCNBFSLzkfHx+IxeI+P8pMd7zyyiuIjo6GmZkZXFxccPPmTbZD6tfq6+uxcuVKpiv06tWrMDU1ZTusHvnhhx9QWlqKbdu2ybUdmggpqh8zMDCAo6MjLly4wHYoMmFgYICwsDB4enrCw8MDv//+O9sh9UupqakYP348zp07h0uXLvWrrlCpmpoa7N+/H2vWrIGZmZlc2+rPN9MDNBFSFGbNmoXz588PmHE8VVVVceLECWzevBlvv/02du3aBUJIh2XPnDmD+vp6BUfYt0m7QjU0NHDv3j14eXmxHdIL+eqrr9DY2IhPPvlEru3U19ejoKAAr776qlzbkStCUS+57OxsAoBcvXqV7VBk7pdffiHKysrk7bffJo2NjW1eCwkJITwej+zZs4el6BRLKBSSurq6Tl+vq6sjK1euJBwOh6xZs4Y0NTUpMDrZKisrIzo6OmTfvn1ybysxMZEAIA8ePJB7W/JCEyFFEULs7OzIv//9b7bDkIvQ0FDC5/PJ1KlTiVAoJIQ83XhpamoSDodD1NTUSF5eHstRyt+SJUvI8uXLO3wtJSWF2NraEgMDA3Lx4kUFRyZ7H330ETEyMiLV1dVyb+vMmTOEy+V2uZPR19GuUYoCMHv2bJw9e7bTLsT+bNq0abh16xYyMzMxfvx4REdHw8PDA42NjSCEQCKRyL37jG1BQUH4448/8L///a/dedPff/8d48aNg5qaGu7duwcfHx+WopSNgoIC/Pjjj9i+fTu0tbXl3l5GRgbMzc2hrq4u97bkhu1MTFF9QUJCAgFAbt++zXYocpOfn0/GjBlDhg4dSpSVlQkA5sHhcMitW7fYDlEuhEIhMTY2JlwulwAgampqJDk5mdTX15M1a9YwXaHPdh33VytWrCBDhw4lDQ0NCmnvnXfeIe7u7gppS154rGZhiuojxowZg2HDhuHs2bNwdHRkOxy5EAgEMDMzw8OHDyEWi9u8pqSkhH//+9+Ii4vrUzOpy8KHH36IJ0+eMEPQicVizJo1C6qqqigsLMT58+cxY8YMlqOUjfT0dPz66684evQoVFVVFdJmRkYGRo8erZC25GVgrfEU1QszZ85EYGAg22HIjb+/P0JCQtolQeBpcnjw4AFOnDjBQmTyc/HiRZw4caLNFcFisRjZ2dnQ0dHBP//8M2CSIPB00l1ra2ssWrRIYW2mp6f37ytGQRMhRTFmz56NzMxMPHjwgO1QZO6///0vDh061OVExIQQrF+/HrW1tQqMTH4qKyuxYsWKDo9wxWIxoqOjB9SgA7Gxsfj777+xf/9+8HiK6ewTiUQoKSnBsGHDFNKevHAIGYBXB1DUC2hpaYG5uTnee+89bN++ne1wZCY6Ohqurq4A8NwZKng8HjZt2iSTWQrq6uqQn58PoVAIkUgEiUSC6upqAICenh4AQFdXF4aGhhAIBDLvylu8eDFOnz7d5f2hKioquHv3Luzs7GTaNhs8PDwgFApx584dhY35GRsbi3HjxiEtLa1fJ0OaCCmqlQ8++AAxMTEDbvb3lJQU/Prrr/jxxx9RXV0NLpfb6dGhsrIy0tPTuz1SSH5+Pu7cuYPk5GQ8ePAAKSkpyM3NZZJedxkaGmLIkCGwsbHB6NGjMXr0aDg6OkJfX79H9QDA+fPnMWvWrOeWU1JSgpWVFeLi4qCpqdnjdvqKW7duYeLEibh27RqmTp2qsHYDAgKwZMkS1NXVQUVFRWHtyhpNhBTVSlhYGKZPn46srCxYWVmxHY7MNTY24sKFCzh8+DDCw8PB4/HaHTEpKyvjzTffxJkzZzqso7KyEsHBwbh69SoiIiKQlZUFLpcLKysr2NraYtSoURg6dChzcY6BgQG0tLSgpKQEHR0dAIBQKAQAVFVVobS0FMXFxcjLy8Pjx4+RlJSE5ORk5OXlgcvlwsbGBpMnT8b06dPh7u4ONTW1Lj/jkydPMGzYMAiFwi6PgJWVldHc3AwtLS38+OOPWLx4cU/+K/sUFxcXaGpqIjQ0VKHt7t27F7/++isyMzMV2q7MsXjFKkX1OU1NTURfX58cOHCA7VDkLi0tjWzcuJHo6+sTDodDlJSU2txSER4ezpQViUTkl19+IR4eHkRFRYUoKyuTiRMnkh07dpCrV6+SmpoamcdXUVFBzp8/Tz7++GMybtw4wuVyiba2Npk/fz4JDAzsdOQXPz+/dreHSB/S5Xw+nyxZsoRcuHBBYbcZyMu5c+cIh8MhMTExCm970aJFxMfHR+HtyhpNhBT1jKVLl5IJEyawHYbCNDQ0kICAADJ58mTC4XCYZGFra0sSExPJBx98QHR0dIiamhqZPXs2+f3335kRahSpoKCAHDp0iLi7uxMul0sEAgHZtm1bm1Fxzp8/3y75qaioEADE2NiYrFmzhoSFhZHm5maFxy8PEomE2NnZkblz57LSvoODA1m/fj0rbcsSTYQU9YzAwEDC5XJJYWEh26EoXFZWFtm8eTPR09MjAAiXyyWvvvoq2b9/PykrK2M7PEZBQQHZv38/MTc3J8rKymTp0qXk/v37xMDAgAAgPB6PACCvvPIK2bp1K4mNjWU7ZLn47bffiJKSEnn48KHC225paSFaWlrkl19+UXjbskYTIUU9o66ujmhqapIjR46wHYrC5efnEz8/P8LhcIiVlRVxd3cnIpGI7bA61djYSI4cOUIsLCyYkWNsbW3J559/TlJTU9kOT66ampqItbU1eeedd1hpPycnhwAgUVFRrLQvS/Q+Qop6hrq6Ojw8PHD27Fm2Q1EYiUSCb775BiNHjsT9+/cRGBiIzMxMhIWFQUtLi+3wOqWiooJVq1YhJCQEc+fOhZ6eHioqKmBlZSXXGdn7giNHjiAvL0/uk+52JiUlBQAGxP8zTYQU1YHZs2fj2rVrzNWNA1l+fj7c3NywefNmfPTRR0hKSsKsWbMUdi+aLAwbNgynTp1Ceno6vLy8sGDBAixduhQikYjt0OSitrYWn376KT744ANYWlqyEkNKSgqMjY1hYGDASvuyRBMhRXXA19cXHA4Hly9fZjsUuQoJCYG9vT1KSkpw584d7N69u1/PImBoaIhffvkFly5dQmhoKMaOHTsgRwo6ePAgampqsGnTJtZiSE1NxciRI1lrX5ZoIqSoDujq6mLy5MkDunv0+PHj8PX1hbe3N+7fvz8gRleR8vLyQnx8PAQCAVxdXREeHs52SDJTWVmJAwcOYN26dTA2NmYtjpSUFJoIKWqgmz17NoKDg1FXV8d2KDJ34MABvPPOO9i0aRN+/fXXfj2qSmcEAgFCQ0Mxbdo0eHh4ICgoiO2QZGL//v3gcrn4+OOPWY0jJSUFI0aMYDUGWaGJkKI6MXv2bDQ0NCAsLIztUGTq2LFj2LBhA7755hvs3bu3X50L7ClVVVUEBATgrbfegp+fH27dusV2SL1SVFSE7777Dps3b2ZG6WFDRUUFysvL6REhRQ10xsbGcHR0HFDdoyEhIVi9ejW2b98Of39/tsNRCC6Xix9++AHe3t548803kZ6eznZIL2zPnj3Q1dXF+++/z2oc0itGaSKkqJfAnDlzEBQU1OUMBv1FcXEx3nrrLSxYsAC7d+9mOxyFUlJSwokTJ/Dqq69i4cKFaGpqYjukHsvOzsaxY8f6xAVNKSkp0NbWhqmpKatxyApNhBTVhfnz50MoFOLatWtsh9Jr77zzDnR0dHD48GG2Q2GFqqoq/vzzT2RkZLB2711vbNu2DRYWFli2bBnboSApKQmjRo0aMN3qNBFSVBfMzc3x+uuv4/Tp02yH0itBQUG4cuUKjh8/Dm1tbbbDAQDEx8fDx8cHurq60NbWhru7O6KiouTaprW1Nb766iv897//RUZGhlzbkqWkpCQEBARg3759Cpt093nx2Nrash2GzNBESFHPMW/ePAQGBvbL7jTg6WzsGzduxNy5c+Hi4sJ2OACAO3fuwNnZGdra2khJSUF2djasrKwwefJkuU8ltGLFCgwbNgybN2+WazuytHnzZtja2uJf//oX26EAAJKTkwdUIqRjjVLUc+Tm5hIOh0MuXbrEdigv5Pz584TL5ZL09HS2QyGEPJ0xwcbGhggEAlJXV8csF4vFZPjw4cTc3FzuUyMFBgYSDodDsrKy5NqOLERGRhIA5MqVK2yHQgghpLCwkAAg169fZzsUmaFHhBT1HObm5njjjTf6bffosWPHMHXqVLz66qtshwIAiIiIQHJyMubOndvmog8lJSUsXLgQeXl5uHjxolxjmDFjBgQCAY4fPy7XdmRh27ZtcHV1hYeHB9uhAHjaLQoAo0ePZjkS2aGJkKK6Yd68eTh37hwaGxvZDqVHRCIRgoOD8fbbb7MdCuP69esAgHHjxrV7TbpM3hcn8Xg8LF26FH/++adc2+mty5cv4+bNm9i/fz/boTCSkpJgYmICIyMjtkORGZoIKaob5s+fj6qqqn53c31MTAyam5vh5uYml/orKyvB4XDaPPbt2wfg6bnJ1svnzp0L4OkYlQBgZmbWrj7p5fiKuNfPzc0NmZmZKCoqkntbL4IQgm3btuHNN9+Es7Mz2+EwBtqFMgBNhBTVLWZmZnB0dOx33aO3b9+GpaUlBAKBXOrX1dUFIQSenp7gcrnIzMxkbk3g8XgghMDJyQknT57EmTNnADxNngA6HNZNOuWTImb9cHR0hJKSEm7fvi33tl7EqVOnEB8fjz179rAdShs0EVLUS0zaPdrQ0MB2KN32+PFjhcwXt379erS0tODrr79uszwqKgoFBQWYN29et+ohhACAQu5P09bWhomJCR4/fiz3tnpKIpFg9+7dWLRoUZ8aDF0ikSAlJYUmQop6Wc2bNw81NTX9qnu0oqJCIfPFubm5wcHBAcePH0dFRQWz/D//+Q/Wrl3b5t43XV1dAE/n1HuWdJm0jLwZGhq2ibevOHr0KLKysrBr1y62Q2kjMzMT9fX1NBFS1MuqPOVBFgAAIABJREFUP3aP1tXVKWw4rnXr1qGuro4ZuSY9PR0RERF4991325STzliQn5/fro6CggIATyfaVQRNTU3U1NQopK3uamhowL59+7Bq1SpYW1uzHU4bSUlJ4HK5A2aMUSmaCCmqB+bPn4/z58/3m+5RPT09hZxvAwA/Pz+Ym5vj+++/R2NjIw4cOICVK1e2G8lmypQpAIDY2Nh2dUiXyevinmc9efIE+vr6Cmmru7777juUl5f3yRv+k5KSYG1tDQ0NDbZDkSmaCCmqB+bOnYuampp+M3O9oaEhysrKFNIWj8eDv78/SktLceDAAQQEBGDNmjXtyk2aNAmjRo3CmTNn2uxQSCQSBAQEwNzcHD4+PgqJuaysDIaGhgppqzuqqqrwxRdfYO3atX1yQOukpCSMGTOG7TBkjiZCiuoBU1NTTJw4sc/ffyZlY2OD+Ph4SCQShbS3atUq8Pl8bNu2DbNmzepwY87lcnH06FE8efIEy5cvR3FxMSoqKvDvf/8bGRkZ+Pnnn6Gmpib3WB8/foyKioo+dWP4V199BYlEgvXr17MdSodiY2Px2muvsR2GzNFESFE9tHDhQgQFBaG6uprtUJ7L2dkZ1dXVzPxx8qatrY1Vq1aBEIJ169Z1Ws7R0RHR0dGoqqrC8OHDYWFhgYyMDNy8eVNhI6hER0dDWVm5wxv72VBWVoaDBw9i06ZNfa67Fnh6S0teXh4cHBzYDkXmOER6vTJFUd0iFAphYmKCo0ePYsmSJWyH0yWJRAKBQIA1a9b0y6mH5MnPzw/FxcUIDw9nOxQAgL+/P06fPo3MzMw+eQ7u2rVrcHd3R1FREUxMTNgOR6boESFF9ZCenh6mTZvWL7pHlZSUsGTJEhw7dgwtLS1sh9NnVFRU4Pz581i+fDnboQB42k175MgRbN++vU8mQQCIi4uDsbHxgEuCAE2EFPVCFi5ciNDQUJSWlrIdynOtWLECOTk5CAoKYjuUPuPw4cNQU1Pr9o3+8rZz504IBAKsWLGC7VA6FRcXh7Fjx7IdhlzQREhRL2DWrFlQVVVFYGAg26E8l42NDebNm4eNGzeiubmZ7XBYV1ZWhq+++gobNmzocJg3RUtLS8OJEyewb98+qKiosB1Op+Li4gbk+UGAJkKKeiGamprw9fXtF92jAPDZZ58hOzsb3333HduhsG7jxo3Q0tLCRx99xHYoAJ5Oujty5EgsXLiQ7VA6VVdXh/T0dJoIKYpqa+HChbh161afHKvyWdbW1tixYwe2bNmCuLg4tsNhzZkzZ3D8+HEcPny4T5yLu3//Ps6dO4fPPvsMXG7f3RwnJCRAIpEM2ERIrxqlqBfU1NQEgUCALVu2dHmrQF/R0tICNzc3FBUVISoqSiFjkPYlDx8+hIuLCxYvXozvv/+e7XAAAO7u7hCJRIiJiVHIQOMv6vDhw9iyZQuEQmGfjvNF9d1dEIrq41RUVDBr1qx+0z3K5XJx8uRJNDY2wtfXt8NBrweqvLw8eHp6YvTo0fjqq6/YDgcAcPXqVVy7dg379+/v88lFen6wr8f5omgipKheWLhwIWJjYxUykawsCAQChISEICsrCzNnzoRIJGI7JLl7/Pgxpk2bBl1dXVy4cEEho9Z0x/bt2+Hl5cWMvdqXDeQLZQCaCCmqV6ZMmQITE5N+c1QIPJ3ZITQ0FA8fPsSkSZP67AztspCQkABnZ2eoqqoiJCQEenp6bIcEAAgMDMSdO3ewb98+tkN5rubmZiQnJ8Pe3p7tUOSGJkKK6gUlJSX4+fnhxIkTbIfSI/b29oiOjkZdXR3Gjx+PiIgItkOSuRMnTsDV1RUjRoxAREQEBAIB2yEBeDraz/bt2zF//vx+MW7ngwcP0NDQMGDvIQRoIqSoXluyZAkyMjJw584dtkPpEQsLC0RHR2PcuHGYOnUqdu/ePSDuM6yursby5cuxdOlSrFixApcvXwafz2c7LMZvv/2GtLS0Pjfpbmfu3LkDbW3tATcHYWs0EVJUL40bNw42Njb4448/2A6lx/T19REY+P+1d+dxNaf9/8Bf55wWiZZRKSFqpijZCi1owSSyDhKF7DfdYxnbYIYZBne2LMU3IZJsY4mhmiylpGhsRWpQKqWkSIvqnM/vj7n1G7eYdJbrnHo/H4/+mNPpul6neTx6+WzXdRJ+fn7YsGEDevbsqdBHh0eOHEHnzp1x7tw5nDlzBlu3bpWrh9SrqqqwZs0aTJ06tXaDYnmXlJSEXr16yfXjHeJqvJ+MEBmaMGECwsLCUFVVxTrKZ+PxePDx8UFKSgrat28PR0dHeHh4yGzHCkmIjY2Fo6MjJk6ciKFDhyItLQ3Dhg1jHesDu3btQl5eHn744QfWUeotKSkJvXv3Zh1DqqgICZEAT09PFBcXIzIyknWUBjMxMcFvv/2GU6dOISUlBV26dMGECRPwxx9/sI5WJ47jEBUVBWdnZzg4OIDP5yMhIQGBgYFy+YzkmzdvsH79evj4+KBdu3as49RLaWkp0tLSqAgJIf+sffv26Nu3r0KeHv1fI0aMwJ07d3D06FGkpqbCysoKvXr1QlBQEF69esU6HvLz8+Hr6wtTU1O4uLiAx+MhJiYGly5dkus/2Fu2bEF5eTkWL17MOkq93bhxA0KhUK5/r5JARUiIhHh5eSE8PBwlJSWso4iNz+djzJgxuHPnDq5evQozMzP4+PigdevWcHNzw759+5Cfny+zPI8ePYKfnx/69+8PQ0NDrF+/HkOGDEFKSgouXryI/v37yyxLQxQVFWHz5s1YvHgx9PT0WMept6SkJLRp0waGhoaso0gVLbFGiIS8fv0a+vr62LFjh1xvp9NQxcXFCA8Px8mTJxEVFYXKykqYmZmhX79+sLOzg6WlJczNzcVew7OkpAQpKSm4e/cu4uPjERMTg9zcXGhpaWHYsGEYPXo0XFxcoKamJqFPJn2LFi3CgQMH8PjxY7Rs2ZJ1nHobPXo0eDwefv31V9ZRpIqKkBAJGjduHAoLC3H58mXWUaSqrKwM8fHxiI2NRUxMDJKTk1FRUQE+n4+OHTvCyMgIhoaGMDAwwBdffFH7ILu2tjaEQiFev34NoVCIV69eoaCgAPn5+cjJyUFmZiays7MBAJqamrC1tUW/fv3Qv39/9OnTB8rKyiw/doM8e/YMX331FdatW4d58+axjvNZ2rZtCx8fHyxbtox1FKmiIiREgsLDwzFy5Eg8fvwYHTp0YB1HZoRCIR4/fox79+7h/v37yMnJQW5uLp49e4aSkhIUFxeD4ziUlJSAz+dDU1MTAoEAGhoa0NPTw/Pnz9GuXTsMHToUXbp0gYWFBYyMjFh/LImYOXMmoqKi8PDhQ6iqqrKOU2/Pnj2DoaEhLl68CGdnZ9ZxpIqKkBAJqq6uRps2bbBw4UJ8//33rOMoDE9PTxQUFCAqKop1FInKyMiAhYUFAgMDMWXKFNZxPsvJkycxduxYFBcXQ0NDg3UcqaKbZQiRIGVlZbi7uyMkJIR1FIXi5OSE+Ph4vH37lnUUiVq5ciWMjY3h6enJOspnS0pKQufOnRt9CQJUhIRInKenJx48eIDk5GTWURTGwIEDUV5ejuvXr7OOIjF3797FiRMnsG7dOigpKbGO89mawoP071AREiJhNjY2MDU1xYEDB1hHURhGRkbo0KFDo7rJaOnSpejZsydGjRrFOspnEwqFuHnzJhUhIaThpkyZgtDQ0EZ3qk+anJ2dG00RxsXFISIiQiE23a3L7du3UVpair59+7KOIhNUhIRIwaRJk/Dq1SucO3eOdRSF4eTkhOvXr6OsrIx1FLEtW7YMDg4OGDBgAOsoDRIXFwdtbW2Ym5uzjiITVISESIGhoSEGDRqE/fv3s46iMAYMGICqqirEx8ezjiKWs2fPIj4+Hhs2bGAdpcHi4+NhZ2fXqHec+Lum8SkJYWDKlCmIjIzEs2fPWEdRCAYGBujUqZNCnx4ViUT48ccfMXr0aNjY2LCO02AJCQmwt7dnHUNmqAgJkZJRo0ZBU1OzUSzELSvOzs64dOkS6xgNFhYWhnv37mHNmjWsozTY48ePkZOT02SuDwJUhIRIjYqKCtzd3bF3717QuhX14+TkhOTkZIVcuLy6uhqrVq2Cl5eXQl9bi4uLg4qKCqytrVlHkRkqQkKkyNvbG+np6UhMTGQdRSE4OTmB4zhcvXqVdZTPtmfPHmRnZyvUprt1iY+Ph7W1tUItai4uKkJCpMja2hpdu3ZFcHAw6ygKoVWrVrC0tFS464QVFRVYt24dZs+eDWNjY9ZxxBIfH9+kTosCVISESN2kSZNw5MgRlJeXs46iEBTxOuG2bdvw6tUrLF++nHUUsRQXF+PBgwdN6kYZgIqQEKnz8vJCeXk5Tp8+zTqKQnBycsLdu3dRWFjIOkq9lJSUYOPGjViwYAFat27NOo5Y4uPjwXEcbG1tWUeRKSpCQqRMT08Prq6u9ExhPTk6OkIgEODKlSuso9SLr68vOI7DggULWEcRW3x8PDp16gRdXV3WUWSKipAQGfD29salS5eQlZXFOorca9myJaysrBTiOmFeXh62b9+O5cuX124+rMji4uKa3GlRgIqQEJkYOnQodHR06KaZenJyclKIIly7di00NTUxZ84c1lHEVllZiZs3b1IREkKkQ1lZGZMnT0ZQUBCEQiHrOHLPyckJaWlpyM3NZR3lozIzMxEUFIRVq1ahefPmrOOILSEhAZWVlXB0dGQdReaoCAmRkVmzZiE3N7fR7cIuDX379oWqqqpcHxX+8MMPMDIygre3N+soEnH58mUYGxujQ4cOrKPIHBUhITJiYmICBwcH7Nmzh3UUude8eXP06dNHboswJSUFhw8fxpo1a6CsrMw6jkRcunQJzs7OrGMwQUVIiAzNmDEDZ8+epYW468HJyQnR0dGsY9Rp+fLl6NKlC8aOHcs6ikSUlZXhxo0bcHJyYh2FCSpCQmTom2++gZaWFu1eXw/Ozs54+vQpnjx5wjrKe5KSknDu3Dls2LCh0WxTFBsbi6qqqiZ5fRCgIiREplRVVeHp6Yk9e/ZAJBKxjiPXbGxsoK6uLnerzCxbtgz29vZwdXVlHUViLl++jM6dO6NNmzasozBBRUiIjM2aNQuZmZly9wde3qioqMDOzk6urhNGRkbi8uXLWLt2LesoEnX58uUme30QoCIkROY6deoEW1tbummmHpycnHDp0iW52MaK4zisWLECbm5ucHBwYB1HYkpKSnDr1q0me30QoCIkhIkZM2bg9OnTKCgoYB1Frjk7OyMvLw9paWmso+D48eO4deuWQm+6W5crV66A47hGVe6fi4qQEAbc3d2hrq6OkJAQ1lHkmrW1NbS0tJifHhUKhVi9ejU8PDzQvXt3plkk7fLly+jWrRt0dHRYR2GGipAQBtTU1ODh4YE9e/bIxWk/eSUQCNCvXz/mRbh//378+eef+Omnn5jmkIamfn0QoCIkhJlZs2bh4cOHCrkbuyy9W3eU1V22lZWV+PnnnzF9+nSYmJgwySAtBQUFSElJadLXBwEqQkKY6dq1K3r16oXAwEDWUeSak5MTioqKcPfuXSbz+/v748WLFwq/6W5drly5UnvU3ZRRERLC0KxZs3D8+HE8f/6cdRS59e76FYvHTd68eQNfX198++23aNu2rcznl7aoqCj07t0bGhoarKMwRUVICEMTJkxAixYtsG/fPtZR5BaPx4OjoyOT64QbN25EVVUVlixZIvO5ZSEqKgqDBw9mHYM5KkJCGFJTU8OUKVOwe/du2p7pE5ycnBATE4Oamho8e/YMhw4dgre3N+zs7CQyfl03LL148QJ+fn5YsmQJvvjiC4nMI0/u3r2L7OxsKkJQERLC3Jw5c5CTk4Pz58+zjiKXCgsLUVVVherqahgbG8PQ0BCTJ0/GgQMHUFJSIpE5bG1tcfDgwfduyPnll1/QrFkz+Pj4SGQOeXPhwgXo6OjAysqKdRTmeBzdu00Icy4uLuDz+bhw4QLrKHIhKioKFy5cQGRkZO3D9MrKyqiqqnrvfc7Ozrh48aJYcz1//hz6+voA/rqBydfXF507d4apqSk2b96MuXPnijW+vHJ2doahoSE9ywo6IiRELsyZMweRkZFIT09nHUUu3LlzB35+fnjw4AE4jgPHcR+UIJ/Pl8gi0X//naempmLw4MEYNGgQDA0NMWPGDLHHl0dv3rzBtWvX6LTof1EREiIH3NzcYGRkRI9S/Nd3330HGxubT256q6SkhNatW4s918OHDyEQCACg9jrtkydP8PjxYwwePBi3b98Wew55Ex0djerqagwaNIh1FLlARUiIHBAIBJgxYwb27t2L8vJy1nGY4/P5CA4OBo/H++h7eDwedHV1xZ4rPT0dSkpK771WXV0NALh69Sp69uyJiRMnyt2+iOKIiIiAlZUV9PT0WEeRC1SEhMiJGTNmoKKiAkePHmUdRS6YmZnhl19++ejmtzU1NRIpwrS0tA9Ou/59Do7jcPjwYSxdulTsueQFPTbxPipCQuSErq4uvvnmG+zYsYN1FLmxcOFC9O7du85TpEKhUCJHNCkpKZ9c71UgEMDDwwOhoaFizyUPHjx4gCdPnlAR/g0VISFyZM6cObh16xaSkpJYR5EL/3SKVNwiFAqFyM7O/uT8s2bNwqFDhz55vVKRXLhwAdra2ujduzfrKHKDipAQOWJvb4+ePXsiICCAdRS5YWZmhnXr1tV5ilTcInzy5Alqamo++v3FixfD39//o6dnFVFERAS+/vrrD66LNmWN5/8uIY3ErFmzcOTIERQWFrKOIjcWLFhQ5ylSca8RPnz4sM7XeTwetmzZgg0bNog1vrypqKhAXFwcXFxcWEeRK1SEhMiZiRMnQk1NDXv37mUdRW7w+XwcOHDgvVOkqqqqUFdXF2vc9PR0qKio1P43j8cDn8/Hvn37sGDBArHGlkfR0dGorKykIvwfVISEyBl1dXVMnz4dO3bsqL2NnwCmpqbvnSKVxPqf6enptTfK8Pl8KCkp4cSJE5gyZYrYY8uj06dPo3fv3hJZiKAxoSIkRA75+PigoKAAJ06cYB1Frrw7RQpAIg/Tp6amorq6GgKBAKqqqrhw4QJGjRol9rjySCgU4uzZsxg5ciTrKHKH1holRE6NHTsWmZmZuHHjBusociU9PR2Wlpawt7fHli1b8Pr1a1RVVaGqqgplZWUAAA0NDQgEAqipqUFDQwO6urrQ09P74KYXPT09FBYWQltbGxcvXkSPHj1YfCSZiImJgaOjIx48eIBOnTqxjiNX6LYhQuTUvHnz0K9fPyQkJMDW1pZ1HJl7+/Ytbt++jfv37yM9Pb32Ky8vD1VVVbh8+fJnFRefz4euri4MDQ1hamqKDh06oLCwELq6uoiJiUHnzp2l+GnYO3PmDDp37kwlWAc6IiREjvXp0wcdO3bEkSNHWEeRupKSEkRHR+Pq1atITEzErVu3UFVVBTU1NZiamtZ+GRoaQl9fH+np6Rg9ejQ0NTWhrKwMZWVltGjRAgDw6tUriEQiVFZW4tWrVygsLER+fj7y8/ORnZ2N9PR03LlzB1lZWeA4Di1atIC1tTVsbGzg6OiI/v37Q01NjfFvRLJMTEzg7u6OdevWsY4id6gICZFj7zagffToEdq3b886jsRlZGTg2LFjiIiIwPXr18FxHLp37w4bGxv06dMHffr0wVdfffXJNUcbKiUlBTo6Onj+/DkSExNx/fp1JCYm4v79+1BTU4ODgwOGDBmCsWPH1m7TpKhu376NHj16IDExkR6krwMVISFy7N1mtBMnTmw0z7S9fPkShw8fRmhoKK5fv47WrVtjyJAhcHFxwaBBg5jvBp+Tk4PIyEhEREQgKioKZWVlGDhwIDw9PfHNN98o5JHi6tWrsWfPHuTk5EjlHxWKjoqQEDm3du1abNmyBdnZ2WI/N8dSRkYGdu7cib1790IkEsHNzQ1eXl5wdXWV21VOKisr8fvvvyMkJARnzpyBmpoaJk+ejEWLFqFdu3as49Vb9+7dYW9vD39/f9ZR5BNHCJFrhYWFnJqaGhcQEMA6SoPcvXuXc3Nz43g8Hmdqasr5+/tzb968YR3rs+Xn53OrV6/mWrduzamoqHAzZszgcnJyWMf6R0+ePOEAcFFRUayjyC16jpAQOaejo4OJEyfCz88PIpGIdZx6y87OxpQpU9C9e3c8e/YM4eHhePDgAebMmaOQR7atW7fGqlWrkJWVhd27dyMqKgpfffUVli5dilevXrGO91GnTp2ClpYWHBwcWEeRW1SEhCiA+fPnIyMjA5GRkayj/CORSISdO3fCwsICcXFxOHToEG7evAk3N7dGsXi1qqoqvL298fDhQ6xbtw779u2Dubk5zpw5wzpanU6fPg03N7f3lpIj76NrhIQoiEGDBkEgECAiIoJ1lI/KzMyEp6cnkpKSsGTJEqxcuRLNmjVjHUuqioqKsHDhQhw8eBDu7u74v//7P2hqarKOBeCvbPr6+ggLC8OYMWNYx5Fbiv/PM0KaiHnz5iEqKgqpqamso9QpKioK1tbWKC0txY0bN7B27dpGX4IA0KpVKxw4cAARERG4evUqevXqhZSUFNaxAADh4eFQUlKiTXj/ARUhIQpi6NChMDc3x8aNG1lH+YCfnx+GDBmCwYMHIyEhAd26dWMdSeZcXFxw8+ZNtG7dGra2tjh37hzrSAgLC4Orq2vtQgOkblSEhCgIHo+H7777DocPH8bTp09Zx6n1yy+/YOHChdiwYQMOHTqE5s2bs47EjIGBAS5dugR3d3d88803+PXXX5llKSwsxOXLl+Hh4cEsg6KgIiREgXh6ekJfXx/btm1jHQXAXyX4448/YteuXVi0aBHrOHJBWVkZe/bswaxZszB+/HicOnWKSY6jR49CTU0NQ4cOZTK/IqGbZQhRMJs2bcLq1auRlZWFVq1aMctx8uRJjBkzBv7+/vjXv/7FLIc8mzNnDg4ePIj4+HiZny62t7eHsbExQkJCZDqvIqIiJETBlJaWon379li8eDGWL1/OJMP9+/fRp08fTJo0iVYr+YSamhq4uLjg0aNHuHXrFrS1tWUy79OnT9GhQwecO3cOQ4YMkcmcioxOjRKiYFq2bIl//etf2LZtGyoqKmQ+v0gkwsyZM2Fubg4/Pz+xx2vRogV4PF69voKCgj76M3w+H9ra2ujWrRvmzJmD5OTkT857+/ZtDB06FFpaWmjZsiUGDhyI+Ph4sT/P3ykpKeHYsWOoqqrCsmXLJDr2p4SFhUFbWxsDBw6U2ZyKjIqQEAU0b948lJaW4uDBgzKfOygoCImJiQgMDISysrLY47158wa3bt0CAIwYMQIcx9X59feVUer6merqaqSlpeHnn39GWloarK2t4e3tjfLy8g/mTExMhJ2dHVq2bIkHDx7gyZMnMDY2hqOjI6KiosT+TH/XqlUrbNmyBUFBQbh27ZpEx/6YsLAwjB07lh6iry8mC7sRQsQ2c+ZMztjYmKupqZHZnFVVVVzbtm25b7/9VqLj3rp1iwPAjRgx4qPvcXBw4Pbs2VPvn1myZAkHgBs+fDgnEolqXxcKhZyFhQVnYGDAlZeX175eU1PDmZmZce3ateMqKysl8Kne5+TkxA0aNEji4/6vBw8ecAC4K1euSH2uxoKOCAlRUEuWLEFWVhZOnjwpszl//fVX5OXlYf78+TKb850rV65g+vTp9X7/hg0b0KdPH4SHh7+3sXFsbCxSU1MxZsyY97ZUEggE8PDwQHZ2tlSeAVy8eDF+//133L17V+Jj/93hw4fRpk0b9O3bV6rzNCZUhIQoKBMTE4wcORK+vr4ym3Pfvn0YPnw4OnbsKLM5fXx8GlS8PB4PPj4+AICAgIDa1y9dugQAsLa2/uBn3r128eLFhkT9pMGDB8PMzAzBwcESH/vvjh49ivHjx0MgEEh1nsaEipAQBbZ06VLcvHmz9o+7NFVWViIuLg4jRoyQ+lyS8u6o6Pr166iurgYApKWlAQDatm37wfsNDQ0BAOnp6RLPwuPx4ObmhujoaImP/c7NmzeRnp5OD9F/JipCQhRYr1694OjoKJOjwuvXr6OiogLOzs5Sm+PMmTMf3A0qzuMZ+vr6AP56jOHFixcAgJKSEgCocyuod0uRFRcXN3jOTxkwYABSUlJQWFgolfHDwsJgYmICKysrqYzfWFEREqLgli5diqioKNy+fVuq8/z555/Q1NSU6s7sdd01Onfu3AaPx33mY9Lv3s/j8Ro856dYWFiA4zg8fvxY4mOLRCIcO3YMHh4eUsvfWFEREqLgXFxc0L17d6xdu1aq8xQUFEBXV1eqc0haXl4egL+WPdPR0QEAaGlpAQDKyso+eP+71969R9L09PQAAM+fP5f42NHR0cjNzYWnp6fEx27slFgHIISIh8fjYeXKlRgzZgzu3bsHS0tLqcxTXl7OZGf5nTt3Nvhn4+LiAAC2tra1zzx26tQJAJCTk/PB+3NzcwEApqamDZ7zU5o1awZlZWW8efNG4mMHBQXB3t4eZmZmEh+7saMjQkIagVGjRsHS0hLr1q2T2hytWrVCUVGR1MaXNJFIVHt98e+nV52cnACgzpVn3r02YMAAqWQqKSlBdXW1xI+si4qKEB4ejmnTpkl03KaCipCQRoDH42HFihU4duyY1DaF1dPTQ2FhIYRCoVTGl7Tvv/8eSUlJGDVqFMaOHVv7uoODA8zNzXHixAlUVlbWvi4UCnHkyBG0a9dOajs25OfnA/j/p0glJSQkBMrKyrQLfQNRERLSSIwZMwbm5ub4z3/+I5Xxu3fvjrdv39YubSZvRCIRCgoKcObMGQwYMAC+vr6YOnUqQkND37t5hM/nY+/evXj58iW8vb2Rn5+PoqIizJ07FxkZGdizZw+aNWsmlYwJCQlo1qyZxE9fBgcHY8KECbQBb0MxXNWGECJhoaGhnEAg4NLS0iQ+tkgk4gwMDLj/g3KXAAAbgUlEQVQNGzZIdFx1dXUOwHtfrVu3/uyf4fF4nKamJmdpacn961//4pKTkz85xh9//MG5urpyGhoaXIsWLThnZ2cuLi5Okh/tA56entzAgQMlOub169c5AFxiYqJEx21KaBsmQhoRoVAICwsL2NraYv/+/RIff/r06UhKSpL6MmGN0evXr9GuXTv89NNPEl2ibubMmUhISMC9e/ckNmZTQ6dGCWlEBAIBli9fjpCQEGRkZEh8/Hnz5uHevXsyWcmmsdm7dy84jsOUKVMkNmZZWRmOHj1KN8mIiY4ICWlkhEIhzM3N0b9/f+zZs0fi4zs5OUEkEuHKlSv04HY9lZaWwtzcHKNHj8a2bdskNu7+/fsxe/Zs5OTkKNwznvKEjggJaWQEAgGWLVuGAwcO4MmTJxIff9u2bbh27ZrUF49uTFatWoXy8nKsXLlSouPu3bsXI0eOpBIUExUhIY2Qp6cn2rVrJ5U1SLt27QofHx8sWrQImZmZEh+/sbl69Sp27NgBX19fiRZWeno6rl27RqdFJYBOjRLSSO3Zs6f2kQAjIyOJjl1WVoa+fftCJBIhPj6ebtv/iKysLPTq1Qv9+/fH8ePHJXoqecmSJQgLC0NmZiZtuSQmOiIkpJGaPHky2rRpI5XnCtXV1XH69Gnk5eVh/PjxePv2rcTnUHSFhYVwc3NDmzZtcODAAYmWYFVVFQ4ePIhp06ZRCUoAFSEhjZSKigpWrFiBvXv3SmW3AyMjI4SHh9fuUVhRUSHxORRVXl4eHB0dUVFRgbNnz0p8jdbjx4+jqKiITotKCJ0aJaQREwqF6NKlC3r16oWDBw9KZY7k5GS4uLigc+fOOH78eO0egE1VSkoKRo4cCWVlZURHR9du9itJffr0QceOHXHkyBGJj90U0REhIY2YQCDAqlWrEBoaijt37khlDisrK8TExCA/Px9WVlZISEiQyjyK4MiRI7CxsUGbNm1w5coVqZTgtWvXkJSUhG+//VbiYzdVVISENHLu7u7o1q0bVq9eLbU5LCwscOPGDfTs2RMODg5YtWpVk7puWFxcDG9vb0yYMAHTp0/HxYsX0bp1a6nMtX37dvTs2RN2dnZSGb9JYra4GyFEZn777TcOAHft2jWpziMSibht27ZxLVq04MzNzbmrV69KdT55cOzYMU5fX58zMDDgTp8+LdW5cnNzOWVlZS4kJESq8zQ1dERISBMwZMgQODg4YNmyZVKdh8fj4dtvv8W9e/fQvn179OvXDyNHjsT9+/elOi8LMTExsLW1hbu7O4YOHYr79+9jxIgRUp3T398f2tra720rRcRHRUhIE7F27VrExsYiOjpa6nN16NABFy5cwPnz55GVlYWuXbvC09Ozzs1wFQnHcYiIiMCgQYPg6OgIDQ0N3LhxA0FBQdDS0pLq3G/fvkVQUBDmzJkDVVVVqc7V1FAREtJE9O3bF66urvj+++/ByehmcVdXVyQnJyMkJASpqamwtrZG//79cezYMYV63OLly5fYtWsXLCws4OrqCgC4ePEiIiMjYWVlJZMMhw4dQklJCWbNmiWT+ZoSenyCkCbk7t276NGjB44fP47Ro0fLfP4rV67Az88P586dg7q6OkaPHo0JEybA0dERysrKMs/zKW/evMGFCxcQGhqKiIgI8Pl8eHh4YP78+bC0tJR5nu7du6Nbt244cOCAzOdu7KgICWli3N3dcefOHaSkpEBJSYlJhvz8fBw9ehShoaG4ceMGWrZsiQEDBsDV1RXOzs748ssvZZ5JKBTi/v37iIqKwoULFxAXF4eamho4Oztj4sSJGDVqFDQ0NGSeCwAuX74MZ2dnJCUloVevXkwyNGZUhIQ0MRkZGTA3N0dQUBAmT57MOg4eP36MCxcuICIiApcvX0ZZWRl0dHTQp08f9O7dG5aWlujUqRNMTEygoqIikTnLy8vx8OFDpKen4/bt20hMTMTNmzdRWlqKVq1aYdCgQXB1dcXgwYOhp6cnkTnFMWrUKLx48QJXr15lHaVRoiIkpAmaMWMGfv/9dzx8+FCubrx4+/YtkpOTkZiYiKSkJCQmJiIzMxMcx0FJSQkdOnRAmzZtYGBgAD09Pejp6aFFixZo3rw5+Hw+NDU1wXEcSkpKAPy1D+CbN29QWFiIZ8+eoaCgANnZ2cjOzq4d09TUtLZ0bWxsYGlpKVfrd2ZlZcHExARhYWF0t6iUUBES0gTl5OTA1NQUa9euxcKFC1nH+aSKigqkp6fXfuXl5SE/Px/Pnz9HYWEhysrKUFZWBqFQiNevXwMAtLW1AQAaGhpQV1eHrq5ubXkaGhrCzMwMZmZmMDY2lthRprR8++23OHPmDB49esTsVHZjR0VISBO1YsUKBAQEICMjAzo6OqzjSISTkxPMzc3h7+/POopE5Ofnw9jYGBs3bsTcuXNZx2m06PEJQpqo77//Hmpqavjpp59YRyEfsXHjRmhqamLq1KmsozRqVISENFEtWrTAmjVrsHv3bqSmprKOQ/7HixcvEBgYiKVLl0JNTY11nEaNipCQJszb2xvdunXD0qVLWUch/2PTpk1o1qwZpk+fzjpKo0dFSEgTxufz4evri99++w0RERGs45D/evHiBQICArBo0SK0aNGCdZxGj4qQkCbO2dkZY8aMwb///W9UVlayjkPw17qw6urqdIOMjFAREkKwdetW5OfnY+PGjayjNHmZmZnYvXs3Vq1aRUeDMkJFSAhB27Zt8cMPP2D9+vV4/Pgx6zhN2vLly2FkZIRp06axjtJkUBESQgAACxcuxJdffkmn4xi6c+cOjh49ivXr18vdIuSNGRUhIQQAoKSkhO3btyMyMhKnTp1iHadJWrx4MXr37o1Ro0axjtKkUBESQmo5Ojpi0qRJmDt3LoqLi1nHaVLCw8MRHR2NjRs3gsfjsY7TpFAREkLes23bNvD5fHz33XesozQZb9++xaJFi+Dh4YG+ffuyjtPkUBESQt6jqamJXbt2ITg4mJ4tlJH//Oc/yMvLg6+vL+soTRIVISHkA8OGDcPYsWMxffr02i2NiHRkZ2fD19cXK1euhKGhIes4TRIVISGkTv7+/qiursby5ctZR2nUFixYAAMDA8yfP591lCaLipAQUicdHR1s3boVu3fvRnR0NOs4jVJERAR+/fVX7NixQ642SG5qqAgJIR81YcIEjBgxAjNnzsSbN29Yx2lUysvLMXfuXIwbNw6DBw9mHadJoyIkhHxSQEAASkpKsGrVKtZRGpUffvgBRUVF2Lp1K+soTR4VISHkkwwMDLBx40b4+fkhLi6OdZxG4c6dO9i+fTs2bdqENm3asI7T5PE4juNYhyCEyDeO4+Di4oLc3FwkJyejWbNmrCPVycnJCebm5vD392cd5aNqampgY2MDdXV1XLlyhR6elwN0REgI+Uc8Hg+BgYHIycmhTXzFtGXLFqSmpiIwMJBKUE5QERJC6qVDhw7YuXMnduzYgXPnzrGOo5BSU1Px448/YtWqVTAzM2Mdh/wXFSEhpN68vLzg4eGBadOmIT8/n3UchVJTUwNvb2/07NkTixcvZh2H/A0VISHkswQEBKB58+bw9vYG3WJQf2vXrkVqaiqCg4MhEAhYxyF/Q0VICPksmpqaOHToEH7//XcEBASwjqMQbt26hXXr1mH9+vUwNTVlHYf8DypCQshns7e3x4oVK7Bo0SLcvXuXdRy5VllZCS8vL9jb28PHx4d1HFIHKkJCSIP88MMP6NGjByZMmICKigrWceTW4sWLkZOTg/3794PPpz+58oj+rxBCGkRJSQmhoaHIzs7GihUrWMeRSxEREfD390dAQAA6dOjAOg75CCpCQkiDdezYETt27ICfnx9+++031nHkSkFBAby9veHl5YUJEyawjkM+gYqQECKWSZMmYfz48Zg+fToKCgpYx5ELHMdh6tSpaN68OXbs2ME6DvkHVISEELHt2rULampqmDBhAoRCIes4zG3duhWRkZEIDQ2FhoYG6zjkH1AREkLEpqmpiVOnTuHatWv4+eefWcdhKj4+HsuWLcOaNWtgY2PDOg6pBypCQohEdOvWDVu2bMHatWsRERHBOg4TBQUFcHd3h4uLC63JqkBo9wlCiERNnjwZ586dQ3JyslTvlNy2bRv27dsHkUhU+1p2djZUVVWhp6dX+5pAIMCaNWswbNgwqWUBAJFIBFdXVzx8+BDJyclo1aqVVOcjkkNFSAiRqLKyMvTu3RstW7ZEbGwsVFRUpDJPfHw8+vbt+4/vEwgEyMvLg66urlRyvLNy5Ups2rQJ8fHxsLKykupcRLLo1CghRKLU1dVx6tQpPHjwAMuWLavzPc+fPxd7Hjs7O7Rt2/aT7xEIBBg4cKDUS/D8+fNYv349duzYQSWogKgICSESZ2pqisDAQPj5+eHEiRPvfS8+Ph5du3bF7du3xZqDx+PB09MTysrKH30Px3Hw9PQUa553srOzsW3btg9ef/r0KSZPnozx48djxowZEpmLyBYVISFEKtzd3TF79mxMnToVaWlpAIDAwEA4OjqisLAQhw4dEnsODw8PVFdXf/T7SkpKGD58uNjzAEBISAjmz5+PmTNn1s5ZXV2N8ePHQ19fH3v27JHIPET26BohIURq3r59C3t7e1RXV2PgwIHYsmVL7fd0dXWRl5cn9pZEnTp1wsOHDz94XUlJCSNHjsTx48fFGh/468jS2NgYmZmZEAgEsLOzw+nTp7Fy5UqEhIQgKSkJnTt3FnsewgYdERJCpEZVVRWBgYF4/fo1tm/f/t73CgsLceXKFbHn8PLyqvP0qFAoxMSJE8UeH/jrdG5mZmbtuNevX4eZmRl2796Nffv2UQkqOCpCQojU3L59G8OHD0dubi5qamre+56ysrJETo96enp+MDbw1007rq6uYo8PAPv373+vbKurq1FSUgJlZWU0b95cInMQdqgICSFSceTIEdjY2OD58+d1Xserrq7GsWPHxN7CycjICD179gSPx6t9TVlZGePGjYOqqqpYYwNARUUFjh49+sFnqKmpQXV1NYYPH/7B0S5RLFSEhBCJ27p1Kzw8PPD27ds6j9beqaiowNmzZ8Web9KkSe9da6yurpbYjg+//vrrR8ua4ziIRCLMmzcPc+fO/eRnJfKLipAQInEzZ87E4sWLIRAIPvl4g0AgwIEDB8Seb/z48e+tMNOqVSs4OjqKPS4ABAUFvXe0WReBQICQkBCJXPMkskdFSAiROHV1dfj6+uLWrVvo0aMH+Hx+nWVSU1ODyMhIvHjxQqz59PT04ODgAIFAABUVFXh5eYl9NyoAZGVlITY29qM7arybY/z48cjIyMDAgQPFnpPIHhUhIURqLC0tcf36dezfvx8aGhpQUlKq832SeMTBy8sLIpEIVVVVGD9+vNjjAcCBAwc+mpnP56NTp06Ii4vDoUOH0Lp1a4nMSWSPniMkhMjEy5cvsXTpUuzduxcCgaD2ehqPx0OvXr2QmJj4WeOVl5cjJycHxcXFKC0tRWlpKcaNGwdNTU2EhYWBx+NBS0sLOjo6MDAw+OwbZziOQ4cOHfD06dP3XldRUYGKigrWrl0LHx8fiRx5EraoCAkhMnX16lVMnz4djx49qj3lyOPx8Oeff8LY2PiD9+fk5CAxMRGpqalISUnBgwcP8PTpU7x+/fqz5tXR0UH79u1hYWGBLl26oEuXLrCxscEXX3xR5/tjYmLeu86opKRU+2zi5s2b39vhgig2KkJCiMxVVVXB19cXa9asAcdxqK6uxpo1a7By5UqUlJTg/PnziI6ORmxsLB49egQ+nw9jY2NYWlrC3NwcRkZGMDAwQNu2bdGqVSu0aNECAoEAUVFR6Ny5M9q0aQMAePXqFQoKCpCfn4/s7GxkZWXh3r17SE1NRXZ2Nvh8PiwsLODo6Iivv/4aAwcORLNmzQAA3t7eOHToEGpqasDj8dCtWzcEBgaiV69eLH91RAqoCAkhzDx69AizZ89GdHQ0dHV10bNnT1y+fBkcx8HW1haOjo7o378/bGxsoK6u/o/j1dTUfPSa3v96+fIl4uLiEBMTg9jYWPzxxx+1D+GPHDkSU6dORWVlJTQ1NbFx40ZMmzYNfD7dVtEYURESQphJTU1FQEAA9u/fj8rKSjg4OGDatGlwc3ODlpaWTLM8e/YMp0+fxqlTp3Dx4kVwHAcrKysEBwejS5cuMs1CZIuKkBAicykpKfD19cXhw4dhbGyMadOmYezYsXj58iWsra1Zx8Ps2bOhpKSE8PBw5OfnY/z48Vi1ahVMTExYRyNSQEVICJGZ3NxcfPfddzh27Bh69OiBH3/8EcOHD//HB9ZZqaqqQnBwMNavX49nz57Bx8cHq1evRsuWLVlHIxJERUgIkTqhUIjt27dj1apV0NPTw6ZNmzBixAi5LcD/VV1djb1792LlypVQVVXF1q1bMW7cONaxiITQlV9CiFTl5ORgwIAB+P7777FgwQLcu3cPI0eOVJgSBP5axHv27NlIS0uDq6srxo8fDy8vL5SWlrKORiSAjggJIVITGRmJiRMnQldXF0eOHEG3bt1YR5KICxcuYMqUKdDU1MTJkyfpZhoFR0eEhBCpCA4OhpubG4YMGYKbN282mhIEAFdXV9y+fRsGBgbo168fYmJiWEciYqAiJIRI3ObNmzF16lQsW7YMBw4cqNczgIrGwMAAUVFRGDRoEFxcXCSynRRhg06NEkIkat++fZg+fTq2bt2KefPmsY4jdSKRCLNnz8ahQ4cQGRmJfv36sY5EPhMVISFEYiIjI+Hm5obly5fjp59+Yh1HZoRCIdzd3XHx4kUkJibC1NSUdSTyGagICSESkZ+fj27duuHrr79GSEgI6zgy9/btW/Tr1w9CoRAJCQlQUVFhHYnUE10jJIRIxNSpU6GhoYGAgADWUZhQVVVFWFgYMjIysHLlStZxyGegIiSEiO3s2bOIiIhAcHBwk151xcTEBJs2bYKfnx8yMjJYxyH1RKdGCSFiqampQdeuXdGlSxccO3aMdRzmhEIhunXrhk6dOuHEiROs45B6qN9+JYQQ8hHnz5/Hw4cPcebMGdZR5IJAIMCaNWvwzTff4PHjx3VuNkzkC50aJYSIZd++fXB2dsZXX33FOorcGDZsGAwMDBAcHMw6CqkHKkJCSIOVlpbi/PnzmDx5MusockVJSQleXl4ICwtjHYXUAxUhIaTBrl+/jurqagwYMECm85aUlIDH4733tXbtWgB/XbP8++tjxoyRabZ3BgwYgD///BN5eXlM5if1R0VICGmwhIQEdOzYEQYGBjKdV0tLCxzHwcXFBXw+H3/++WftIwtKSkrgOA62trYIDQ1ldsOKjY0NBAIBEhISmMxP6o+KkBDSYFlZWTAzM2M2/8KFCyESibBly5b3Xo+Pj8fTp08xduxYRsmAli1bQl9fH1lZWcwykPqhIiSENFhRURFatWrFbP6vv/4alpaWCA4ORlFRUe3rGzduxL///W8oKyszywYAOjo67+Ui8omKkBDSYOXl5VBTU2OaYf78+SgvL69d0SY9PR2XLl3CzJkzmeYCAHV1dbx584Z1DPIPqAgJIQ2mra2N4uJiphkmTpyI1q1bY+fOnXj79i02b96MyZMnQ1tbm2kuAHj58iW++OIL1jHIP6AiJIQ0mI6ODgoLC5lmUFVVxZw5c1BQUIDNmzcjNDRUbrZ/KiwshI6ODusY5B9QERJCGszCwgK3b9+GUChkmmPOnDlQU1PDypUrMXDgQHz55ZdM8wB/3UhUVFSELl26sI5C/gEVISGkwezs7PD69Ws8ePCAaQ4dHR14enqC4zgsXLiQaZZ3rl27BmVlZVhbW7OOQv4BFSEhpMEsLS2hq6uL06dPs44CW1tbWFlZoX///qyjAABOnz4NW1tbNG/enHUU8g+oCAkhDSYQCODp6Yl9+/ZBJBIxzbJ79265ORosKirCmTNn4O3tzToKqQcqQkKIWKZNm4bMzEycPXtWpvMGBQVh1KhRePPmDXbv3o3i4mKMGzdOphk+JiAgAM2aNWP6QD+pPypCQohYLCwsMHbsWCxduhTV1dUynfv06dPQ1tbGrl27cOTIESgpsd9ZrrCwEJs2bcLixYuhrq7OOg6pB9qYlxAitkePHsHc3Bzr16+Xm9OTrEydOhWRkZHIyMig64MKgo4ICSFiMzExwY8//ojly5fj1q1brOMwc+LECQQHByMgIIBKUIHQESEhRCJEIhEGDBiAvLw8xMfHM12DlIX79+/D3t4eEydOxM6dO1nHIZ+BipAQIjF5eXmws7ODvr4+oqOjm8w1suzsbNjb28PIyAi///47mjVrxjoS+Qx0apQQIjEGBgaIjIzEo0ePMGLECJSWlrKOJHVZWVkYNGgQtLS0EB4eTiWogKgICSESZWpqiqioKNy/fx8ODg6Neof2O3fuwM7ODqqqqoiMjJSLhb7J56MiJIRIXPfu3XHt2jWUl5ejd+/eiI2NZR1J4kJDQ9GvXz906tQJsbGxMDAwYB2JNBAVISFEKjp06IBr167B2toazs7O+Omnn2T+nKE0vH79Gt7e3vDy8sK0adNw4cIFaGpqso5FxEA3yxBCpIrjOPj7+2Px4sX48ssv4e/vLzfrgX6uI0eO4LvvvkNVVRX27duHYcOGsY5EJICOCAkhUsXj8eDj44OUlBS0b98ejo6O8PDwYL5jxeeIjY2Fo6MjJk6ciKFDhyItLY1KsBGhIiSEyISJiQl+++03nDp1CikpKejSpQsmTJiAP/74g3W0OnEch6ioKDg7O8PBwQF8Ph8JCQkIDAxscs9INnZUhIQQmRoxYgTu3LmDo0ePIjU1FVZWVujVqxeCgoLw6tUr1vGQn58PX19fmJqawsXFBTweDzExMbh06RJ69+7NOh6RArpGSAhhKi4uDrt378aJEycAAAMHDsTo0aMxZMgQ6OvryyTDo0ePcPbsWZw8eRLx8fHQ0NDApEmTMHPmTFhYWMgkA2GHipAQIheKi4sRHh6OkydPIioqCpWVlTAzM0O/fv1gZ2cHS0tLmJubi72GZ0lJCVJSUnD37l3Ex8cjJiYGubm50NLSwrBhwzB69Gi4uLhATU1NQp+MyDsqQkKI3CkrK0N8fDxiY2MRExOD5ORkVFRUgM/no2PHjjAyMoKhoSEMDAzwxRdf1D7Irq2tDaFQiNevX0MoFOLVq1coKChAfn4+cnJykJmZiezsbACApqYmbG1t0a9fP/Tv3x99+vSBsrIyy49NGKEiJITIPaFQiMePH+PevXu4f/8+cnJykJubi2fPnqGkpATFxcXgOA4lJSXg8/nQ1NSEQCCAhoYG9PT0oK+vj7Zt26Jdu3bo0qULLCwsYGRkxPpjETlBRUgIIaRJo7tGCSGENGlUhIQQQpo0KkJCCCFNmhKA46xDEEIIIaz8P58z85JJGs97AAAAAElFTkSuQmCC\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",
    "\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",
    "\n",
    "**Natural indirect effect**: Effect due to the path v0->FD0->y (mediated by FD0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "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, [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": "stdout",
     "output_type": "stream",
     "text": [
      "Estimand type: nonparametric-nie\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 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",
    "#### Natural Indirect Effect\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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "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",
      "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: 19.43485415766647\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import dowhy.causal_estimators.linear_regression_estimator\n",
    "causal_estimate_nde = 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_nde)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the value equals the true value of the natural indirect effect (up to random noise). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19.43485415766647 19.42468979102858\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 direct effect is set to be zero. \n",
    "\n",
    "#### Natural Direct Effect\n",
    "Now let us check whether the direct effect estimator returns the (correct) estimate of zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "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, [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: y~v0+W0) - ((b: FD0~v0+W0)*(b: y~FD0+v0+W0))\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 0.0006101936244320427\n",
      "\n"
     ]
    }
   ],
   "source": [
    "causal_estimate_nie = 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_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.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
}