{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conditional Average Treatment Effects (CATE) with DoWhy and EconML\n",
    "\n",
    "This is an experimental feature where we use [EconML](https://github.com/microsoft/econml) methods from DoWhy. Using EconML allows CATE estimation using different methods. \n",
    "\n",
    "All four steps of causal inference in DoWhy remain the same: model, identify, estimate, and refute. The key difference is that we now call econml methods in the estimation step. There is also a simpler example using linear regression to understand the intuition behind CATE estimators. \n",
    "\n",
    "All datasets are generated using linear structural equations.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import logging\n",
    "\n",
    "import dowhy\n",
    "from dowhy import CausalModel\n",
    "import dowhy.datasets\n",
    "\n",
    "import econml\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "BETA = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         X0        X1   Z0        Z1        W0        W1 W2 W3         v0  \\\n",
      "0  0.680855 -0.081822  0.0  0.363662 -1.457851  0.486804  2  3  21.904925   \n",
      "1 -1.366921  1.606660  1.0  0.505995  0.022552 -0.744353  3  3  36.342609   \n",
      "2  0.448970 -0.942918  0.0  0.790112  0.488717 -0.630550  2  0  13.715890   \n",
      "3 -0.413223 -1.072336  0.0  0.999187 -2.757718 -0.282076  2  0  13.279910   \n",
      "4 -0.133890  0.752749  1.0  0.307046 -0.829915 -1.434104  3  1  21.662710   \n",
      "\n",
      "            y  \n",
      "0  282.284051  \n",
      "1  244.554705  \n",
      "2  143.071317  \n",
      "3   78.660784  \n",
      "4  219.933669  \n",
      "True causal estimate is 10.230462989973262\n"
     ]
    }
   ],
   "source": [
    "data = dowhy.datasets.linear_dataset(BETA, num_common_causes=4, num_samples=10000,\n",
    "                                    num_instruments=2, num_effect_modifiers=2,\n",
    "                                     num_treatments=1,\n",
    "                                    treatment_is_binary=False,\n",
    "                                    num_discrete_common_causes=2,\n",
    "                                    num_discrete_effect_modifiers=0,\n",
    "                                    one_hot_encode=False)\n",
    "df=data['df']\n",
    "print(df.head())\n",
    "print(\"True causal estimate is\", data[\"ate\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = CausalModel(data=data[\"df\"], \n",
    "                    treatment=data[\"treatment_name\"], outcome=data[\"outcome_name\"], \n",
    "                    graph=data[\"gml_graph\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwAAAADPCAYAAABC17ZyAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydd1hUx/f/D9IW2EaTJh2UKih2FBWJFRsGkxhLNCqWRGMsmKgxJuYjdlLUoLEQE01Q0WBQoxgLYsdGUexK7yx9F3bfvz/87X5ZQaXc3QWzr+fZ54F7755z7tw7s+fMnJnRAABSo0aNGjVq1KhRo0bNf4J2qjZAjRo1atSoUaNGjRo1ykMdAKhRo0aNGjVq1KhR8x9CS9UGqGkd1NTUUHl5+SvPSyQSEggEzZLN4XBIS+vVr5q+vj7p6uo2S7YaNaqgoqKCRCIRlZSUkDSLksfjkba2NnG5XBVb17oRCoVUWVkpK0MiIgMDA9LR0SEej0ft2qn7pV6FtB2WliERkY6ODhkYGMjKUM2rEQgEVFNTQ6WlpURE1K5dO+LxeKSrq0v6+voqtq51I62vAoGAJBIJEb1o87S0tIjH46nYutaNSCSiiooKuTZP6veoss1TBwCtiPLyciovL6eKigoqKSmhsrIyEgqFVFpaKnPQq6urqaqqiiorK0koFFJ5ebmsQROLxTKHpKysjGpra4mI5F66qqoqqq6uJiKS+xFpLnw+nzQ0NGT/13WImgOLxSI9PT0i+r8fNiIiLS0t4nA4RPR/jbZUf7t27YjP55OmpiZxuVzS1tYmNpstkyWtaGw2W1bhuFwuGRgYEJvNVjtsahqktraWEhMT6fbt25ScnEx3796lzMxMSk9PlwuWuVyuzKEgItLT0yNLS0uysrIiFxcXcnd3J09PT+rZs+d/ysm4d+8eJSYmUnJyMqWmptKTJ08oJyeH8vPzZdfo6+tTTU0N1dTUENGLem5mZkYdOnQgR0dH8vDwIA8PD+rZsye1b99eVbeidHJycujy5cuUkpJCycnJ9OjRI8rMzKS8vDxZu66trU1aWlpUVVUl+56ZmRmZmZmRg4MDubm5kaenJ/n4+JCzs7OqbkXplJeX05UrVygpKYlSUlLo3r17lJmZSdnZ2bLfPg0NDWKz2VRWVib7HpfLJSsrK7K2tiY3Nzdyd3cnb29v8vb2fm0H1ttEbW0t3bx5k27duiVr8zIyMigjI6NeWZWVlcl+61ks1ivbPOlv+H+BtLQ0unHjBiUlJdHdu3fp0aNHlJubS3l5ebJr9PT0qLa2tl6bZ2VlVa/NMzMzU6i9GupJwMxQWVlJxcXFcp+ioiK5/0tLS6m8vJxKS0tJIBDIHP7y8nIqLi6WyZI6pmw2u0HHVtpboaenRywWi9hstqznUVNTU+aU13WUiYgMDQ1lf9eNOhvqoZfKUgQAqKSkRO7YyyMMdcujbo9D3cCmpKSEJBIJlZSUkFgsptLSUlmk/apASSAQUEVFBZWXl8s1aDweT1bmbDabDA0NZX9zuVwyNDSUfYyMjOT+l16rpu3z6NEjio6OplOnTtHFixeJzWZTly5dyMPDg9zc3MjGxoYsLCzIzMyMDA0N6/XclJSUUH5+PuXk5NDz58/p7t27lJqaSjdv3qTMzEzq3r07+fv709ixY6lr164qukvFUFhYSEePHqVjx47R+fPnqby8nLp27Uru7u7k4eFBdnZ21KFDBzI3Nyc+n19v1E/aDmZkZFB2djbdv3+fkpOTZY5cp06dqH///jR69GgaOHDgW9XbLRQKKS4ujmJiYujcuXP06NEj8vT0lDkDHTt2JHNzc7K2tiY+n1/PqRIKhVRcXEzZ2dmUmZlJT548oeTkZEpJSaEbN24Qn8+n/v3704gRIygwMJD4fL6K7pR5ANC1a9fo8OHDdObMGUpMTCRbW1vy9vYmNzc3cnV1pQ4dOpCFhQWZmprW660Wi8VUXFxMubm5lJWVRc+fP6fU1FRKTk6mGzdukFAoJF9fXxo8eDCNHTuW7OzsVHOjCuLx48d0+PBhOnnyJF28eJH09fXrtXmWlpZkZmYm80fqUlJSQgUFBZSdnU3p6elybV56ejp169ZN1ub5+PjIdRi2dYqLi2Vt3rlz56i0tFSuzbO3tycrKysyNzcnQ0PDem1eRUUFFRcXU3p6OuXk5Mi1ecnJyeTs7EwDBgygUaNG0aBBgxhv89QBQAOIxWLKz8+XfXJycqigoIDy8/Nl0VxBQYGcoy8UCklDQ6OeY1j3w+fzicPhyDmWPB5P5vBzOJx6PepqFEtJSYlcIPby/wKBgEpLS+sFd3U/YrGYdHR05J61iYkJmZqakrm5OZmamso+ZmZmsr+1tbVVfftq6EUj/uuvv9Lu3bspLS2N3nnnHRoxYgT5+fmRq6srY3rS09Pp3Llz9M8//9Dff/9NfD6fJk2aRDNmzCBra2vG9CiTmpoa+uuvv2jHjh105swZ6tq1K40ePZoGDBhA3bp1Y+wdFwgEFB8fT6dPn6bDhw+TQCCgcePGUUhICHXv3p0RHarg4sWLFBERQYcPH6b27dvT2LFjyd/fn/r27Ssb8WwpIpGIrl69SmfPnqW//vqL7ty5QwEBATRz5kwKDAxUWEePonny5Alt376dfvvtN6qsrKSRI0fSkCFDqH///mRpacmIDgCUmppKZ8+epdjYWDp9+jR17tyZpk6dShMnTmyzo8cCgUDW5qWmplJAQACNGDGC+vfvT66uroz5IJmZmbI27+jRo8ThcGRtnq2tLSM6lE1tbS3FxMTQjh076N9//yUvLy8aM2YMDRgwgLp3785Ym1daWkrx8fH077//0uHDh6moqIiCgoIoJCSEevbsyYiO/1QAIBKJKDs7mzIyMigrK4uysrIoIyODcnJyKD09Xc7pJyIyNTWVOXLt27eXOW8mJibUvn37ej3Bb1OviprG01CAkJeXJwsUpUFkfn6+7H+JREKGhoZkZmZGJiYmsh4qa2trMjc3pw4dOsiGVFkslqpv8a3k6dOnFBYWRnv37qXu3bvTjBkzaNSoUYw5Xq9DJBJRXFwc7dy5k2JjY2n48OG0fPnyNjMqUFZWRj/99BP98MMPxOVyacaMGfTBBx+QlZWVUvRfv36dIiMjae/eveTs7ExffvkljRkzpk10nkgkEvrzzz9p7dq1lJ6eTlOmTKEpU6aQl5eXUvQ/e/aM9u3bR7/88guJRCJasGABzZo1q82kp12+fJm+++47iouLo9GjR9O0adPI399fKWk6AoGAjhw5Qtu3b6c7d+7Q1KlTKTQ0VGnvfUt5/vw5rV27liIjI6lr1640Y8YMGj16tFICmZqaGjp9+jTt3LmT/v77bxo8eDAtX768zQTw5eXltG3bNgoPDyd9fX1Zm6eszpubN2/Snj17aO/evWRnZ0dffvklBQUFtWz+AN4icnNzcfnyZfzxxx9Ys2YNQkJCMGLECHh5ecHMzAxEBB0dHdja2sLX1xfBwcGYP38+1q1bh7179yIuLg537txBdnY2ampqVH07at5SxGIxcnNzkZycjDNnzmDfvn3YtGkTFixYgPfffx99+/aFg4MDWCwWiAjGxsbw8PDAsGHD8PHHH2P16tX47bffcOHCBWRmZkIikaj6ltoU+fn5CAkJgZ6eHqZMmYLU1FSV2pOVlYXQ0FBwOByMHj0a9+/fV6k9r0MkEmH9+vUwNjaGn58fTp48qdL3r6KiAlu3boWNjQ28vLxw8uRJldnSGGJiYuDq6goHBwf88ssvqKqqUpktYrEYsbGx6N27N8zMzPDTTz+htrZWZfa8iZSUFAwdOhRcLhcrVqxAbm6uSu25ffs2PvjgA+jr62P+/PkoLi5WqT2vo7CwEHPmzIGenh4mTpyI5ORkldqTnZ2NL7/8EhwOB4GBgbh3755K7XkdNTU12Lx5M0xNTeHr64sTJ06otM2rrKxEREQE7Ozs4OnpiePHjzdbVpsKAIRCIVJSUhAdHY2NGzfik08+QWBgINzd3aGvrw8NDQ1YWlrC19cXEydOxPLlyxEREYGjR4/i5s2byMnJUfUtqFHTJPLz83Hnzh3ExsZi586d+PrrrzFlyhT4+fnBxsYGmpqaYLFY6NSpE4YOHYpZs2Zh7dq1iIqKwq1bt1BZWanqW2g1SCQS/PLLLzA2NkZQUBAePHigapPkyM/Px/z582FgYICvvvoKQqFQ1SbJceHCBXh4eKBz5844ffq0qs2RQygU4qeffoKhoSHGjx+PrKwsVZskx/PnzzF69GiYmppi+/btEIlEqjZJjuPHj8PV1RVdu3bF1atXVW2OHJWVlQgNDYW+vj4WL16MwsJCVZskR2pqKkaMGAEzMzPs3btX1ebUY8+ePTA1NcXo0aNbnaNdUFCAzz//HAYGBli2bJlKA+KGuHTpEry8vODh4dHqOhdEIhG2bdsGY2NjjBs3DhkZGU2W0SoDgJycHJw9exYRERFYuHAhRowYAScnJ2hpaYHD4cDHxwfjx4/HkiVLsG3bNhw/fhz37t1DdXW1qk1Xo0apiEQiPHz4EKdOncKOHTvw5ZdfYsKECejZsycMDQ3Rrl072NnZYciQIZg3bx62bt2KuLg4pKenq9p0pVJQUIDRo0fD1tYWx44dU7U5r+XWrVvo0aMHfHx8WsVoQG1tLVauXAkOh4NNmza16tHRvLw8TJo0Ce3bt281z/nQoUMwMjLCjBkzWp3zWheRSIQ1a9aAzWYjLCwMYrFY1SYhOTkZnp6e6Nu3L1JSUlRtzms5cuQIrKys8P7776OkpETV5qCoqAjjxo2DtbU1jh49qmpzXsudO3fQu3dveHt74+7du6o2B7W1tfj222/BZrOxbt26Vhew16WgoABTp06FiYkJYmJimvRdlQYAQqEQN27cwO7duzF//nwMGDBA5rTY29tj6NChmD9/PrZt24bTp083K8JRo+a/TG5uLs6dO4ft27dj0aJFGDlyJJydnWXBdJ8+fTBnzhxERETg6tWrb+WIQXJyMmxsbBAcHNyqh+nrIhKJsGTJEvB4PJw4cUJldggEAgwaNAienp6t3gGry759+8Dj8fDtt9+qzAaJRIIlS5bA2NgYhw8fVpkdTeXmzZvo1KkTRo4ciYqKCpXZER0dDQ6Hg5UrV7bq1KS65OfnIzAwEM7Oznj48KHK7Lh79y7s7e0xduzYVh101qWmpgbLli0Dl8vF33//rTI7ysrKMGTIELi5ueHOnTsqs6OpHDhwAIaGhlixYkWjU5SUFgAIhUJcvHgRmzdvxuTJk+Hl5QVtbW0YGRnB398fn3/+OSIjI3Hr1i11T74aNQpGJBIhJSUF+/fvR2hoKIYMGQIzMzNoamrCzc0NH3zwAdauXYtz58616aAgPj4eRkZG+O6771RtSrM4cOAAOByOSlILsrOz4eXlhVGjRqnUEWwuqampsLW1xaxZs5Tem11TU4OJEyfC2dkZjx49UqpuJhAIBAgICEDPnj1V4kBu27YNXC4XsbGxStfdUiQSCb744gu0b98e169fV7r+S5cuwdjYGF9//bXSdTPB4cOHweFwsGvXLqXrzs3NhY+PD4YNG4by8nKl628paWlpsLe3x7Rp0xoVNCssACgoKEBMTAyWLl2Kfv36gcViwdzcHGPGjME333yDmJgYPHv2TFHq1ahR0wyys7Nx/PhxrFmzBuPHj4e1tTW0tbXRs2dPLFiwAAcPHkR2draqzWwUN27cAJ/Px86dO1VtSos4e/YsuFwuDh06pDSdJSUl8PLywkcffdRmel8bIjMzEy4uLpg3b57SdEokEkydOhVdunRBXl6e0vQyjVAoRHBwMHr16qVUZygyMhKGhoatbi5CUwkPD4eJiYlSU1ru3LkDQ0ND/Pzzz0rTqQguXLgAHo+HP/74Q2k6S0tL4ePjgw8//LBVpzm+iezsbHh4eGD27NlvvJaxAKCqqgonTpzAvHnz4OrqCk1NTXh4eCAkJASRkZEqHQ5To0ZN80lPT8f+/fvx6aefomvXrtDS0oKjoyNmzpyJw4cPo6ysTNUm1iMjIwNmZmbYvHmzqk1hhGPHjoHD4eDSpUsK11VbW4uBAwdi7Nixbdr5l/Ls2TN06NABGzduVIq+lStXwsnJ6a1YdEIoFOKdd97BqFGjlLLySVxcHDgcDs6fP69wXcpg1apVsLW1RX5+vsJ1ZWdnw9LSEmvXrlW4LmVw6tQpsNlsxMfHK1yXWCzG4MGDERgY2KadfykZGRmwtbXFmjVrXntdiwKAp0+fYtu2bRg5ciT09fXh6OiITz/9FLGxsW0m11aNGjVNo6ysDHFxcVi8eDHc3d2hq6uLgIAAbNq0CWlpaao2D2KxGP7+/pgxY4aqTWGUH374Afb29gqfYPjNN9/Aw8OjTad+vUxiYiLYbLbCe5XPnj0LDofTKiYyMoVAIICjoyM2bdqkUD35+fmwsLDAnj17FKpH2bz33nsYOXKkQgMoiUSCIUOGYMqUKQrToQp+/vlnWFtbKzwNLSwsDC4uLm0y7edV3L59G2w2GwkJCa+8pskBQGVlJaKiohAQEAAtLS34+voiLCxMJbluatSoUT1Pnz5FREQEgoODweFw4ObmhrCwMJWlP2zZsgUuLi5tMm/9dUgkEowYMQIzZ85UmI6kpCQYGBggKSlJYTpUxcaNG+Hi4qKwHr6qqirY2dlh+/btCpGvSq5cuQIDAwOFjuRPmDAB77//vsLkq4qSkhLY2toiMjJSYTp27NgBJyenVjka21LGjh2Ljz76SGHy7969CwMDA9y8eVNhOlTFDz/8ACcnp1cuKd3oAOD8+fOYPHkyDAwM0KNHD2zbtk2lvfxHjx4FEck+mpqacueHDBkiO6erq8u4/vnz58vkK2N4T1X0798fRARbW9tGf2f//v3o378/eDweDAwM4O7ujunTp+PUqVMKSyn48ccfYWtrC01NTRAR1q9frxA9TaU55deWqaiowK+//ooBAwaAxWIhODgYx48fV9rGKaWlpTAzM2v2yjnGxsZy7UrdNm7cuHHo1q2b3PU9e/aUXbt7927Z8aSkJAwYMAB6enqwtLTEihUrGHn3nzx5An19fYVt5DN06FCEhoY2+/tMlV9cXBzMzMywbNmyZtvyMrW1tfDw8MDWrVsZk1mXsLAw9OjRg5F3vbi4WK4cG/q8PMlTEWVWl9mzZ2P8+PEKkX3lyhVwOBxG9m9oTtkpqr5KiYqKgpWVlUI6JcrKymBhYdHkJSClMFFnb968iQ8//BAmJiYwMDBA165dGVu4ICMjA2w2G4mJiYzIe5lRo0ZhwYIFzf5+S8uvuLgYYWFhcHV1hY6ODszNzTFhwgRGluoWi8Xo2rUrwsPDGzz/xgDgn3/+QZ8+fWBiYoLPP/9c5TvIvYyVlRXmzp3b4LkDBw688hwTnDp16q0PAABgzZo1jXZgp02bBgsLC0RFRUEgEKCiogK3bt3Chx9+CCJq0a51ryI3NxdEhL1797bKnOWmlN/bxKNHj7BixQpYWVmhc+fOOHDggMIDgfXr12PgwIEtkjFixAh06tRJ7lhtbS34fD40NDTqjWyMHz9eLvUpJycHJiYmmDt3LkpLS5GYmAhzc3MsXry4RXZJmTdvHiZPnsyIrLokJiaCx+O1uGOnJeUnFouxcuVK9OrVC1paWow7s4cPH4aNjQ3jqwIJhUKYmpoytkFacXEx3N3d6x2vrKyEq6ur3AiXostMSk5ODvT09BSS5jdmzBjG7G5K2QGKr6/Ai9G7bt26ISIigjGZUr7//nv4+vq2SEZL2zwzMzOMHz8ez58/h0AgwPr16+sF9S1h0aJFChkdSkpKApvNRkFBQYvktKT8Fi5cCFtbW1y+fBlVVVW4desWnJyc0KdPnxbZJCU2NhaWlpYNjny+MgDIyMjAmDFjYGJigg0bNrTa3Ch1AKB4GuvAbt++HRoaGq9MB+vdu7dCAoCkpCQQUatNW/ivBgBShEKhLJdzwIABCs2PdnNzw59//tkiGVu2bAER4fHjx7JjCQkJcHNzAxHh999/lx2vra2Fj4+P3PdDQ0NhYmIit3nMjz/+CG1tbUbSou7duwc9PT3GR2Dnzp2LkJCQFstpSfk9ePAA69evR21tLXR1dRl3ZsViMaytrfHPP/8wKvfAgQPo2LEjYwHuq5zYmTNnQkdHRy5dQdFlVpcPPvgAS5cuZVRmbm4udHV1GVsVsCllByi+vkrZsWMHevTowZg8Kd7e3i3ubW9pm2dlZVWvPerRowd69uzZIrukPHr0CCwWi3Ffa8GCBZg6dWqL5bSk/BYuXFhv1abPPvsMJiYmLbYLeBF8Ojo6NrgZXDtqgDNnzlCXLl3IyMiI0tLSaOHChWRgYNDQpW2KDRs2kIaGBmloaNDx48dp4MCBpK+vT97e3nTt2jW5awHQxo0bycnJiXR0dMjBwYHCwsJIIpHUk5ucnEw9e/YkFotFLi4udPToUbnz586do169epGBgQFZWVnRF198QcXFxXK61q1bR/b29qSjo0OOjo60YcMGAlDP7tjYWBo0aBCxWCzS0NCgy5cvy85paGjQ0qVLiYho9erVsmPh4eGN0kNEJBQK6ZNPPiFDQ0Pi8/n0ySefUE1NTaPKd9OmTdS7d2/y8fFp8PzFixdp6NChjSrfxj6rDRs2kKenJxEReXp6koaGBv3222+N0jFr1izS0NAgFoslkzdmzBjS0NAgDw8POR2NfW8aW34teebJyclvfKdaEzo6OhQSEkJ3796l7t27U8+ePengwYOM60lLS6Nnz57R6NGjWyRn2LBhRER07Ngx2bF//vmHlixZQubm5nTixAnZ8cuXL1OvXr3kvn/8+HHq27cvaWtry44FBARQTU0NxcXFtcg2IqJOnTqRq6srnTx5ssWy6nLkyBGaMGFCi+W0pPycnJxo0aJFpKmp2WI7GqJdu3b03nvv0eHDhxmV+9dff9GECRNIQ0ODEXl8Pp+Sk5Pljh05coS2b99OYWFh5O3tLTuu6DKry4cffsh42cXGxlLPnj3JxsaGEXlNKTsixddXKePHj6ebN29SZmYmYzKfPHlCd+/epaCgoBbJaWmbl5GRQXw+X+6Yvb095eXltcguKQ4ODuTt7S1nBxO0hjZvw4YNFBISQkRE5eXlFBcXRwcPHqQvvviixXYREWloaLy6zXs5Irh48SK4XK5CJ6wwSVNHAKS9xYGBgXj48CHy8/Ph5+dXL6Jdvnw52Gw2/vrrL5SWluLvv/8Gh8ORGxaUjgCMGDEC9+7dQ0FBAebMmQMtLS1ZL6dIJAKXy8WaNWtQUVGBtLQ0eHh4yPXSLF++HO3bt8eFCxdQUVGB06dPg8fjyeWqSu0eOHAgkpKSUFJSgm7duiEpKQnV1dUNbvyxZMkSucizMXoWLlwIDoeDY8eOobS0FPv374e1tfUbe7ALCwtBRI0ecWlM+Tb2Wb1qBKAxOnbv3l1vjsiyZcvq9SA11pbGll9LnvmNGzfe+E61Zo4dOwYul8v47qg7d+6Ev78/I7I6deqEESNGyP7v3bs3cnJy8NFHH8HMzEzW07ts2bJ6u1bq6enVW3e+oqICRIRvvvmGEfs+++wzzJ8/nxFZwIu5Bbq6uqiqqmJEXkvKT4qierNjYmLg6enJqEw7OzucOXOGUZl1ycjIgJGREYYOHfraUQZFjwAIBAJoamoy2jM+bdo0rFixgjF5L/OmslNGfZXSvXt3REVFMSZv79696Nu3LyOymKizdenatSs+/vhjRmwDXvgzjVnbvrFkZmZCS0uLsYnTLS2/zZs3g4jAZrNfmbPfXE6cOFEvRQl4KQVIJBKhU6dO+OmnnxhVrkiaGwDU3VTnl19+gZaWliwvtLy8HCwWq94P7KJFi6CjoyMb6pIGAKdOnZJdU1VVBT6fL1uCMD8/H0SEixcvyq7Zu3evrJGW6np57d4FCxbAxsamnt2vajwWLVoEKysrWQ68UCiEm5ubbFflxuipqKiArq4uPvvsM7lrgoOD3xgAPHz4EETUqOHhxpZvY55V3evqBgCN1dHUAOB1tjS2/Fr6zN/0TrUFjhw5AhMTE0aXtFy4cGGLJnPV5bPPPoO+vj6qqqpQVFQkG8r+888/QUS4ceMGAKBPnz5y+cQSiQREhCVLlsjJE4vFIKIWTbCty549exAQEMCILOBFnqi3tzdj8ppbfnVRlDP7/PlzaGlpMZauU1lZCSJCUVERI/JeRiwWY+DAgWjfvv0b9xZQdAAAAE5OTjh37hxj8nr37o2DBw8yJq8ubyo7ZdVXKTNnzsTKlSsZk/fFF18wlubMRJ2VcuXKFbDZbLmUmJayb98++Pn5MSbv1KlTcHNzY0weE+Un7Qy0s7NjNHjKycmBhoaGXJob8FIK0NWrV6miooJmz57NyNCDMmjXrl2DaTlERBKJhNq1azDLiRwcHGR/GxkZUW1tLVVWVhIRUUpKClVXV5OXl5fcd7p06UIikYiSkpLkjru6usr+ZrFY1LFjR9kQpImJCX377bf07rvv0hdffEHHjh2joKAgWr16tZyu0NBQuVSezZs30/Pnz6msrExOV8eOHRu8n5CQEMrKypKlHx06dIhGjhxJurq6jdbz8OFDEgqF1LlzZznZ7u7uDeqsi5GRERFRPXsboqnl+7pnxZSOxvI6Wxpbfi195m96p9oCo0ePJkdHRzp+/DhjMgsKCsjU1JQRWcOGDaPKyko6d+4cnTp1igYPHkxERIMHDyZNTU06ceIE5ebmEpfLJX19fdn3NDQ0SE9Pj6qrq+XkSf9nKpXS1NSUCgoKGJFF9KLsTExMGJPX3PJTBqamplRbW8tYylxBQQFpaWnVS4FgirVr19LZs2dpz549ZGZmphAdTcHU1JTy8/MZk5efn89YvX2ZN5WdsuqrFEXUW1W3eS8jEAho2rRp9Msvv5C9vT0jthH9N9o8fX198vf3pw0bNtDOnTvpypUrjNhmYmJCGhoaVFhYKHdczjsuLi6m9u3bv9Jpbo1wOJxXNuRFRUXE5XIbPKelpSX7+1V5m43N53w59xKA3HeXL19Ot27doh49etAff/xB1tbW9fKxdu3aRXgxIiP34XA4r9UlxcnJiQICAmjbtm1ERM4DggkAACAASURBVLR9+3ZZXllj9eD/55+/fN+NeR8MDQ3JxcWFbt68+cZrpTS2fBvzrFqiA3XmQBDRa+c8vM6WppZfS555Y96p1o6FhQUVFRUxJk8kEpGOjg4jsvr3708GBgZ07Ngx+ueff2j48OFE9CK/uFevXnTixAk6ceKEbE5LXZydnSk9PV3u2PPnz4noRT1lAhaLVc9paQkikUjWWcAELSk/RaOjo0Pt2rUjoVDIiDyhUEg6OjqM5f/X5erVq/TVV1/RvHnzZHnGqkZXV5exsiNitt7WpbFlp4z6KkVXV5fxetsa2jwpAoGARo8eTSEhIfTee+8xYpeU/1KbJ53XePv2bUZs09TUJC0trXr1Vs4z6datG6WlpdHdu3cZUaoMPD096eLFiyQWi+udi4+Pl5vM2Vjc3NyIxWLRrVu35I7funWLdHR0ZA9Hyv3792V/C4VCevDggazX9969ezR8+HAyNTWlsWPH0q+//kpTp06lJUuWyOm6ceOGnMy4uLh6PclvYtasWXTq1CmKiYkhPT09uei7MXqcnZ1JV1eX7ty5I3dNWlpao/QvWbKELl26RImJifXO5eXlEZfLpUOHDjW5fJtDY3VwuVwSiURyIwoPHz5sls7Gll9Ln/mb3qm2QF5eHp05c4Z8fX0Zk2lkZMRYr66uri4NHDiQjh8/TomJidSjRw/ZuWHDhtGlS5fozz//bNCxGDp0KF24cEEukDx9+jRpaWlRQEAAI/Yx2fNH9KLsmAzGWlJ+iqaoqIgkEgkZGxszIs/Y2JgqKysZdU6IXoymTpgwgdzd3Wnt2rVy50pKSsjFxYVRfY2lsLCQ0Z5Tpt89oqaVnTLqq5TCwkLG621raPOIXkwE9vf3p5kzZ9Knn37KiE11eVvbPGdnZ8rOzpY7JvUpmWqjSktLSSQS1a+3L+cKffnll3B3d28zS1ueOXMGGhoaGD9+PFJSUlBZWYmnT5/i22+/hbW1db3lSxvKFz98+DCISG4yyPLly8HhcHD06FGUlZXJJi4uWrRIdo10DsCQIUPw4MEDFBYW4pNPPpGbBHz37l1oaWnhwIEDqKysxLNnz+Dr6yuXv7t8+XLo6+sjKioKZWVluHXrFlxdXfHrr7++1u6XqampgaWlJTgcToObgjRGz8KFC8HlcnHixAmUlZXh0KFDMDU1bfQylnPnzoWFhQUOHDiA0tJSlJWV4fTp0/D09JTbDr0x5dvYZ/W6ScBv0vHw4UO0a9cOmzZtQnl5OWJiYmBvb//KOQBvsqWx5deSZ96Yd6o1U1FRAT8/P8bXsl+zZg0mTJjAmDzp0m4vy0xMTAQRwcHBocHvZWdnw9jYGJ9++inKyspw8+ZNWFhYYOHChYzZtnbtWkbXxb5y5QrMzc0Zkwc0v/ykKCqf/erVq2jfvj1j8sRiMdhsNlJTUxmTCQCTJk2Cnp5eg3KLi4sbnNSn6DkAYrEYPB4PKSkpjMkMDAzEli1bGJMHNK3slFFfpYwaNQo//vgjY/I2btyId999lzF5za2zt27dgouLi9x8SODFpldM+ZKbN29GUFAQI7KAF5uXGRsbM7o3TXPKz9HREcOGDcO9e/dQVVWFhIQEODk5wdbWlrHl92/dugVDQ8N6x+sFAGKxGBMnToS9vT2uXbvGiHJF8/fff8PPzw9cLheampqyTSkePHggd92OHTvkdmzbv38/vv32W7lj0pUcJBIJ1q9fD3t7e2hpacHOzg7/+9//ZBM+6+4EHBsbiy5dukBHRwedOnWSc74lEgn279+P3r17g81mw9jYGO+++y6eP38ud82GDRvg4OAAHR0ddOzYUW5d2JftHjRo0CvLYuXKlbCxsWlwQ6w36QGA6upqzJ07FzweDxwOB1OnTsVXX30l092YNdz37duHPn36wMDAAGw2G15eXti4caPcBJQ3lW9jn5V0wxHpx8zMrNE6pGzevBkWFhbgcrlYuHAhvvjiC5m89PT0Jr03jS2/ljzzxrxTrZW0tDR4eXlh8ODBjK04I+XMmTOM7rfw+PFj2QZzdZFIJDA3N3/t5Lvbt2+jf//+YLFYMDc3x7JlyxjdpG7kyJH4/vvvGZMnEomgp6eHR48eMSazueVna2tbb+fWbdu2MWbX5s2bMWbMGMbkAcCgQYPwyy+/MCbv7Nmzb9zNtq4Tq+gyk3Lr1i3w+XxGN1JbvXo1o4F7U8sOUHx9BV6896ampq/cJ6c5JCQkMBq4N7fO8vn8V5Y1UwHAuHHjsH79ekZkAS/W42ez2YzuS9Oc8ktKSsLUqVPh6OgIFosFe3t7zJgxAxkZGYzZtWXLFgwfPrze8QY3ApM6JwYGBvj000/bzGiAGjVqWielpaVYsWIFDAwMEBoa2uCuhC2loqICbDYbt2/fZlx2a6KsrAwcDofxje8CAgIYX36uNTJo0CBGgycA+O677zB69GhGZbZGvv32W8aDp0uXLsHExES2Yt3bSkJCAoyNjRlt+6qrq8Hj8dpMZ21zqaysBJ/PR2JiIqNyhw8fjnXr1jEqszUybNiwBoOnBmcnamho0MKFC+n27dv0/Plzsre3p8WLF9OjR48aulyNGjVqGiQrK4tWrVpFdnZ2FB8fTxcuXKCwsDC5ydRMoa+vT+PHj6ddu3YxLrs1ERUVRZ06dWrW/KbXMWXKlLe+7B4/fkwJCQmMbP5Tl0mTJtGJEycoJyeHUbmtCYlEQrt376apU6cyKrdXr15kbGxMMTExjMptbezcuZMmTpzIaNunq6tLH3zwwVtfbw8dOkQ2NjbUtWtXRuVK2zy8tBDI20R6ejr9+++/NHHixPonGxM93LhxA++99x5YLBYGDBiAX3/99Y3rwKpRo+a/iVAoxKFDhzBixAjo6upi5MiRjK4b/jquX78ODoeDrKwspehTNiKRCB07dsSePXsYl11ZWQkLCwscOXKEcdmthY8//pjxuSdSxowZw+jmbK2NvXv3wsHBQSGjdz/88AO8vb0ZTS1qTTx+/Bj6+vqMpptIuX37NgwMDJCens647NZAbW0t3NzcsH37dsZlV1dXo0OHDoxuztbamD179ivnizUqAJBSUFCA77//Ht7e3uByufj4449x+PBhxnZSU6NGTdukqqoKJ06cwCeffAJTU1N07NgRa9asQWZmptJtef/99xndRKU18f3338PT05Px/GQp27dvh4uLy1uZjnHr1i3o6+vj6dOnCpGfmpoKfX19pKWlKUS+KikvL4etrS327dunEPlCoRCOjo7YtWuXQuSrmvfee0+2OagimDx5MiZNmqQw+apk27ZtcHFxUUjgCbzYVNHJyQmVlZUKka9KUlJSoK+vj4cPHzZ4vkkBQF1u3LiBRYsWwc3NDTo6OggICMDGjRsVEuGqUaOm9fH06VNs3boVgYGB0NfXh4ODA+bOnYv4+HhGV1ZoKs+ePYOhoWGDK2G1ZVJTU8HhcBQ6mlJbW4tevXoxtqNya6GyshIeHh5YtWqVQvUsXrwY3bt3r7fjZlvn448/hr+/v0Lr9bFjx8Dn8xndPbY18PvvvzdqF+eWkJmZCWNjY4XtqKwq7t+/Dy6XW291ISYRi8Xo27cvYzsqtxaqq6vh7e2N5cuXv/KaZgcAdXny5ImcI2Bra4uJEydi69atuHPnzls7rKdGzX+Je/fuYdeuXZg2bRo6deoEXV1dvPPOO9i0aRPu3bunavPkiIqKgomJCe7fv69qUxihuLgY7u7uWLFihcJ1PXr0CDweDwcOHFC4LmUgkUgwefJk+Pn5KWzkRIpQKISPjw8++eQThepRJrt27YKpqalSRvPmz5+P7t27M7b8oaq5c+cOeDwejh8/rnBdhw8fhpGR0VvTCSsQCODl5YXQ0FCF65J2GilqhEsVTJ8+HX369HntyAkjAUBdqqqqcPr0aaxatQpDhgwBl8sFj8fDsGHD8M033yAuLg5FRUVMq1WjRg2DCAQCnD9/HmvXrsWoUaNgamoKAwMDDBw4ECtWrMCxY8da/Y/0smXLYG9v3+bnA1RWVsLPzw9BQUEKd2ClHD9+HBwOB//++69S9CmS0NBQODk5KbQHti7Pnj1Dhw4d8N133ylFnyI5evQoOBwO4uPjlaJPKBRi8ODBGDZsWJsfRXny5AksLS2VusrMqlWrYGNj0+bnA1RXV8Pf3x8jR45UWOrPy8TFxYHNZuPkyZNK0adIVqxY0ajfPsYDgJepra3F7du3sWXLFnz44YdwdHSEhoYGbGxsMGrUKKxYsQIHDx7Ew4cPVZo2oEbNf5Vnz54hJiYG33zzDcaNGydXR4ODgxEeHo6rV68qrSFmColEgpCQEDg6OtbbE6StUFRUhH79+sHf35/xfRPexO+///7KTQXbAmKxGPPnz4eVldUrc2AVRXJyMkxNTbF8+fI2+7v2xx9/gMPh4K+//lKq3rKyMvTq1QuDBw9GaWmpUnUzRVJSEqytrbF48WKl6/70009hZ2fX6kZlG0tJSQkGDhyIfv36KX2xmaioKHA4HERHRytVL1NIJBIsXrwY5ubmjZqLpPAAoCGkvYs//vgjpk+fjm7duoHFYoHL5aJv376YMWMGNmzYgL///hsPHjxQWq+XGjVvKxKJBE+ePMGJEyfw/fffY/bs2Rg4cCCMjIygra0NLy8vTJ48GRs3bsTp06dRWFioapMZQSKR4Msvv4SpqWmb682+e/cu7Ozs4Ofnp3TnX8pff/0FDoeDzZs3tylHtqSkBEFBQXB1dVXYpN83cf/+fTg4OGDixImtfrSsLmKxGKtXrwaPx1No7vXrKC8vx/Dhw9GlSxdGN6dTBrGxsTA0NMT//vc/peuurKzEli1bwOfzwePx2lxvdlpaGjw9PTFmzBiVTMo9c+YMOnfuDF1dXaxbt65NtXmlpaUYP348Onbs2Oh5NCoJABqipqYGycnJ2LdvH5YvX47g4GB4eXmBxWJBR0cH7u7uCAoKwhdffIE9e/YgPj4eGRkZ6vkFatTUITs7G5cuXcLvv/+OFStWYPz48fD29oaenh50dHTg6uqKsWPHYunSpdizZw9u3LgBoVCoarMVzu7du8HhcLB8+fI2kVqwa9cucLlcjB49GiYmJvDz85PtNq1srl69Cnt7e4wcORJ5eXkqsaEpXLp0CQ4ODggMDFR5IJuTk4OAgAC4uLjgxo0bKrWlMWRlZeGdd96Bi4sL7ty5o1Jbampq8Pnnn4PP57eJ3OyqqiosXLgQPB4Pf/75p1J1l5aWIjw8HJaWlvDy8kJkZCQiIyPB4XCwdOnSNtHG//rrr+ByuVi6dKnSO31PnTqFLl26gMViQVtbGyNGjICTkxOGDRumtNTBlnDt2jU4Oztj2LBhTdq4t9UEAK9CLBbj8ePHOH78ODZv3oxZs2Zh4MCBsLW1haamJnR1deHs7Ix33nkHM2fOxP/+9z/s27cPFy9ebPO5v2rUvExeXh6uXr2KqKgorFu3DnPmzMGwYcPg6uoKFosFDQ0NWFlZoV+/fpg+fTrWrVuHmJgY3L9/v82l8DBNamoqunbtCjc3t1Y7GpCSkoKBAwfC2toap0+fBvCiNzQ8PBxmZmbw9fVFXFyc0u0qKSnBhAkTYGRkhG3btrXKjpfCwkKEhITAwMAAmzZtajW9d2KxGGvWrIGBgQE+++wzCAQCVZtUj5qaGmzevBk8Hg/Tpk1rVUt7x8bGwsLCAkOHDm21qXzHjh2Dk5MTevbsqdR0s7y8PKxcuRKGhobw9fVFTEyM7L0/cOAATExM4OLigk6dOqlsNOdN3L17FwEBAbC0tMQ///yjNL1isRgxMTFwcXEBi8WClpYW3n//fVngW1paismTJ4PP5+PHH39slZkoRUVFmDt3LgwMDJo1YtHqA4DXIRKJkJmZievXryMqKgphYWGYOXMmAgIC4ODgAE1NTejo6MDCwgI+Pj4IDAzEzJkzsXLlSkRERODUqVNITk5GSUmJqm9FzX+c6upq2bscExODiIgIhIaGYtKkSQgICICbmxs4HA6ICIaGhvDx8UFwcDBCQ0MRHh6OqKgoXL9+vU2lGqgCsViMyMhIGBsbIyAgAFeuXFG1SQBeTBicOXMm9PT0MG/evAZzn8vKyhAWFgYjIyP4+vrKAgRl8u+//8LNzQ3u7u6IjIxsFT+KdcslMDCw1S4j+fDhQwwfPhwmJiZYuXJlqwgExGIxoqKi0KlTJ3Ts2LHVpoyUl5cjNDQUenp6mDRpktLndLyKhIQEBAYGwtDQEOHh4UqrD0+ePMG8efNgYGCAwMBAXLp0SXYuNzcXwcHBaN++PQ4cOCBr80xMTODr66uykcSXefr0qewe5s2bp7T6IBKJsHv3bnTo0AG6urpgsViYM2cOMjIyGrz+7NmzcHd3h6ura6tq86SdQoMGDWr2yk9tOgB4E0KhEE+fPsWFCxcQFRWF8PBwLFq0CB9++CH8/PzQsWNH6Ovrg4jA5/Ph4uKCfv36YezYsZgzZw6++uor/Pjjj/jzzz9x9uxZJCcnIzc3t9X0LKlp3RQUFCA1NRXnzp3D3r17sWXLFnz99deYO3cugoOD4efnBzc3NxgbG4OIwGKx4ODggL59++L999/HggULsGnTJuzbtw/nz5/Hw4cPVZYL/rZRUFCAZcuWgcvlwt/fH1FRUUofJpdIJDh9+jSCg4Ohp6eHadOmNSrfWerwSnv9lD2aIRKJsGPHDtjb26NTp07YvHmzSlZ2S01Nxbx588Dn8+Hv76+03aZbSlxcHPr27QsTExMsXrxYJb3a+fn5WLduHRwdHeHs7IzIyMg2MUKYlpaG4cOHQ09PDxMmTFDJM6+ursZvv/2Gvn37gs/n4+uvv0ZxcbFSdCclJWHSpEmyQCglJUXuvHT54+Dg4HqpIEVFRfjqq6/A4/EwYMAA7N+/XyVt3pkzZ/Dee+9BT08PU6ZMUdr7X11djbVr14LH40FbWxsWFhbYunVro+Ya1NTUYNeuXbL6snHjRhQUFCjBannu3buHzz//HIaGhujfv3+LgzkNAKD/OCUlJZSZmUn5+fmUk5ND+fn5sk9OTg4VFBRQfn4+5ebmUklJCbVr145MTEzIyMiIDA0N5T6vO8bhcIjD4aj6dtU0gYqKCiorK6OioiIqLi6W+7zuWEFBAdXW1hKHwyFjY2N69uwZcTgcsrS0pI4dO5KHhwe5u7tT+/btyczMjCwsLMjExETVt/ufo7i4mHbv3k0RERFUXFxMY8eOpXHjxtGAAQNIR0eHcX0SiYSuXbtG0dHRdODAAaqqqqJp06bRzJkzydbWtkmyysvLacuWLbR27Vpyc3Oj1atX04ABAxi3+VXU1NRQdHQ0/fzzz3TlyhUaOnQoBQUFUWBgIPH5fIXofPz4MUVHR9OhQ4fo9u3bpKurS1u2bKEJEyYoRJ8iuXDhAm3bto2io6OpW7duFBQUREFBQU1+DxpLUVERHT16lKKjo+nUqVPUt29fCgkJodGjR5OWlpZCdDJNbGwsffDBB7Rp0ya6f/8+7d69m4yMjGjcuHE0btw46tq1K2loaDCuVygU0unTpyk6OpqOHDlC5ubmFBISQpMnTyYej8e4vpe5ceMGhYWF0dGjRyk4OJi++uorcnJykp3PycmhOXPm0KVLl2jr1q00duzYV8oqKSmhPXv2UEREBBUUFNCYMWNo3Lhx5O/vr5A2DwBdv36doqOj6eDBg1RWVkYfffQRhYSEkL29PeP6XkYgEFBoaChFRkaSUCgkb29vWr16NQ0dOpTatWvXJFm1tbV0+PBh+vnnn+nSpUs0ePBgCgoKopEjR5KhoaFC7H/y5AkdPnyYDh06RDdv3qR3332XZs2aRX369GmxbHUA0EREIpEsIJA6e3UdwZedwrr/S+Hz+cRms4nNZpOBgQEZGhrK/mez2cTj8YjD4RCbzSY9PT3icDikpaVFfD6f2rVrR3w+n7S0tIjD4ZCOjg4ZGBiQnp4esVgsFZaM6hCJRFRRUUFVVVVUXV1NFRUVJBKJqLS0lMRiMZWUlJBYLCaBQEBCoZCys7NJIpFQWVkZlZeXU3l5OZWWlpJAIKCKigoqLy+nsrIyKikpIWn14PF4DQZ50r9fPmZqakomJiayZyIQCOjatWt04cIFSkhIoIsXL5KWlhb16NGDfH19ycfHh/r166cwx0nN6wFACQkJdOjQIYqOjqbCwkLq3bs39evXj3x8fMjd3Z3s7OyaLDc3N5eSk5Pp9u3bdP78ebpw4QIREbm5uVFJSQldvnyZ9PX1W2R7WVkZbd26lcLCwsjd3Z2+++476t+/f4tkNpX79+/TwYMHKTo6mm7dukVdunShfv36Ua9evcjDw4OcnZ1JW1u7STLLysooNTWV7ty5QwkJCXT+/HnKysqid955h8aNG0djxoyhPXv20Jo1a+js2bPk6uqqoLtTLIWFhRQdHU3R0dH077//ko2NDfn5+VGfPn2oc+fO5OrqSmw2u0kyRSIRpaWlUUpKCl2+fJnOnz9PSUlJ5OPjQ+PGjaN3331XKc4Xk/z22280e/Zs2r17N7377rtE9MIxP378OEVHR9PRo0dJW1ub+vXrR/369SMvLy/y9PRscscKAHry5AmlpKRQYmIinTt3jq5cuUIWFhayQKNnz56KuMV6XLhwgdauXUvnz5+njz76iJYuXUoWFhZy1+zbt48++eQTGjNmDG3atKnRvyEA6NKlS7I2Ly8vj3r16kV+fn5ybV5TA6q8vLx6bZ5EIqFRo0ZRUFAQDRkyRCGBxsukpqbS/Pnz6d9//yUtLS0aOXIkbdiwoVnteEM8fPhQ1ubduHGDvLy8yM/PT9bmdezYscltXnl5eb02LzMzkwYNGkTjxo2jsWPHMhpoqAMAJVJcXEzl5eUyJ7OkpETmhFZUVJBAIKDS0lKZU1pWVtagE1tbW0tlZWUyx7cuBgYGpKOjQ9ra2rIfDU1NTeJyuUREpKGhIddA8Pl8WQXncrmkqan5Svt1dXWb7KxUV1dTVVXVK88DoJKSEtn/dZ1u6b0TvagYNTU1RERUWVlJQqFQVg51kQZC0nKQ3pM0eNLX16fY2FiytLSkHj16kI+PD3G5XDIwMJAFX1wuVxaccblc4vF4Te4peBO1tbWUlpZGCQkJdOHCBUpMTKS0tDTq1KkT+fj4UN++fcnX15fc3NwU0qOl5vXcu3ePzp8/T/Hx8XT79m1KS0sjHR0dsrGxIQsLC7K0tJQLzqXvsUgkotzcXMrIyKDMzEwSCARkb29PXl5e5OvrS35+fuTt7U01NTXUv39/6tatG23ZsoURm+sGAt26daPvvvuOevTowYjsppCTkyMru+vXr1NKSgoJhUKytrYmCwsLsrKyIh6PRywWi/T09IjoRV2vra2l7Oxsys3NpczMTMrKyiJLS0vy8PCQBWO9e/eu1wZ99dVXtHv3bjp//nybc2pfpry8nC5evEjnz5+nK1euUFJSEuXl5ZGlpSV16NCBzM3NyczMTNYBRESyjo/i4mLKysqirKwsysjIIH19fXJ3d6du3bqRn58f9evXj9q3b6/iO2we33//Pa1atYqOHDlCfn5+DV5TW1tLiYmJFB8fTwkJCZSUlERPnjwhIyMjWfmZmZmRtra2zImqqamh8vJyqqyspMzMTMrJyaFnz56RWCwmV1dXmVPXr18/cnZ2Vsq9isViOnDgAK1du5by8vJowYIFFBISUi97oKSkhD799FM6efIkbd++nUaPHt0ivWlpaRQfHy9r8+7evUva2tqyNs/Kyor09PSIzWaTtra2rM2rqamh3NxcSk9Pp6ysLCouLiYHBwfq3Lkz+fr6Ur9+/ahLly5KGWUSi8X0+++/0+rVq+nBgwfUvn17WrhwIX3++ecK1Z+bm0vx8fF0/vx5WZtXVVVF1tbWZGlpSZaWlsTn8+XavLKyMlnZ5eTkyNo8c3Nz8vDwkAVjvXv3JgMDA4XYrQ4A3gJe1fNN9KLnWSKREFHDTjRRw056Q871y7zsvBO9OYgg+r8gpS51g4vGBC/Sv18eCWkMeXl5tHv3bvr555+ppqaGpk+fTnPnziVTU9NGfV9R5OTk0LVr1ygxMZESEhIoISGBdHR0qHv37uTr60t9+/alPn36tLjHWE3TqampoYcPH1JGRobMSZUG5aWlpcThcEhDQ4MMDQ3JxMSEzM3NqUOHDtSxY0dZg/8y2dnZ1L17d1q+fDnNmjWLMVuLiorohx9+oPDwcOrevbvKAoG6PH36lNLT0yk9PZ1yc3NJIBBQdXU1VVZWkra2Nmlra1NlZSXt3r2b9u7dSxYWFuTs7ExGRkaNkh8aGkp//PEHnT9/XmEpNKqisLCQ7t+/Tzk5OZSRkUECgYD+/vtv6tSpE/H5fNLT0yM9PT3i8Xhkbm5OVlZWZGtrSzY2Nqo2vcUAoKVLl9Kvv/5Kx48fJ29v7yZ9v6Kigh48eCCrt4WFhbJOprKyMuJyuaSlpUU8Ho/MzMzI0tKSrK2tydHRUelpUdXV1RQZGUnr168nDQ0NWrx4MU2ZMoV0dXXrXXvy5EmaNm0a9erVi37++WeFpI/W1tbSw4cPKT09XdbmSTsgX27zjI2NycLCgjp06EDOzs5K/41KSkqiH374gfbv30+VlZXk5OREGzdupJEjRyrVjro8e/aMnj9/ThkZGZSTkyPX5mlpaZGOjg5xOBzi8Xiy4N7Z2ZmMjY2VZ2SLZhCoUdOGEYvFOHXqFAIDA8FisRAcHNyqlkqT7o0RERGBSZMmyVa2cnNzw8yZMxEZGYnk5GRVm6mmBVy8eBFsNlshK3MUFBRg5cqV4PF4CAgIwNWrVxnXwSQ1NTVgs9nNWtFCIpFg9uzZcHZ2RnZ2tgKsa11YWVnh8uXLqjZDodTU1GDatGlwcHBotct/MoF0DX8rKyvZGv6vmpBdWVmJefPmgcfjISIiQsmWti4KCwsREREBb29vaGpqQlNTEwMHDqw3MVrNq1EHAGrUAHjwzZpl+QAAIABJREFU4AFCQ0NhZGQEHx8fREREKH0b8saQmZmJmJgYhIaGwtfXFywWC+bm5ggMDMTKlStx6tQpleygqKb5SJclVdSOp9JAgMvlIiAgANeuXVOIHiYYMGAAdu7c2azvSiQSTJ8+HZ6enipZoUOZ6Onp4f79+6o2Q2FUVFRg+PDh8PHxQW5urqrNUQivW8O/IRISEuDo6IiAgAA8f/5ciZa2HqqrqxETE4Pg4GBoa2vDyMgIurq6mD17NtLT01VtXptDHQCoUVOHqqoqREZGwtvbGzweD/PmzWvVW9GLRCJcv34d4eHhmDRpEuzs7KClpaUeJWhjzJ8/H15eXgrdxyE/P18uELh+/brCdDWXZcuWYfr06c3+fm1tLd5//3106dJFaUszKpuqqioQkcp3OVYURUVF8PX1hb+/f6vYK4FpXl7D/+LFi6+9vra2FqtWrQKbzcYPP/zwn1uGXCwW48yZM5g+fTp4PB7s7Ozg6uoKNpuNefPm/SdG/BSFeg6AGjWvIDExkbZv306//fYb9enTh2bOnElBQUFvnOOgarKysmTzCC5cuEDXr18nIyMjucnF3bp1+8+uGtUaEYvFNHLkSGKxWHTo0CGFTvwuKCign376icLDw6lfv360atUq6tq1q8L0NYXY2FhasmQJpaSkNFtGTU0NjRs3jgoLC+mff/5p8go6rZ3MzEyysbGhmpoaxhcnUDVZWVk0dOhQcnJyon379r1VbVRSUhKtX7+eDh48SO+++y4tXbqU3NzcXvud9PR0mjhxIuXn59P+/fvJy8tLSdaqntu3b9Pvv/9O+/fvp9raWurevTs9fvyY8vLyaM6cOfTZZ5+pV81rKaqOQNSoae3k5OQgLCwMtra2cHR0RFhYWL1NVlozdUcJgoODYWZmBi0tLfj4+GDevHmIjIzEkydPVG3mf56ioiI4Ozvjm2++UYq+/Px8hIaGynoiExMTlaL3dRQVFUFTU7PFvdtCoRDDhg1DQEDAW7d53p07d2BkZKRqMxgnNTUV1tbWmDt3LsRisarNYYz4+HgEBgbKeqyfPXvWqO9FR0fDyMgIkyZNapXpqIrg+fPnCA8PR5cuXaCnp4dx48bhs88+g4uLC+zs7PDTTz+9dfVZlagDADVqGklDk4YTEhJUbVazyMzMRFRUFObNmwdfX1/o6OjAwsICgYGBCAsLQ3x8PKqrq1Vt5n+OlJQU8Hg8HDp0SGk68/Ly5AKBGzduKE13Q7i6uiI2NrbFcioqKuDn54ehQ4e+Ve/ymTNn4OTkpGozGOXKlSswMTFBaGioqk1hBIlEgpiYGPTp0wdcLhehoaGNnpcinejL5/Px559/KthS1VNUVITIyEgEBARAS0sLAQEB2LlzJ8LDw2FnZwcXFxfs2bMHIpFI1aa+dagDADVqmsH9+/frTRpuy5Nvy8vLER8fLxslMDU1hba2tmyUICoqCjk5Oao28z/BiRMnwOFwcOXKFaXqzc3NlQsEbt68qVT9Uj7++GMsX76cEVkCgQDdu3dHUFDQK1dWaWscOnQIPXv2VLUZjHHy5EnweDz8/PPPqjalxYjFYkRFRcHd3R3m5uZYuXIlSkpKGv395ORkeHp6YsCAAW/1pNbi4mLs2bMHw4cPh66uLvr374/t27cjIyMD4eHh6NChAzw9PREZGYna2lpVm/vWog4A1KhpAWVlZYiIiEDnzp3B5/Mxb968tyad5uVRAm1tbVhYWCA4OBjh4eGIj4+HUChUtZlvJZs3b4aFhYVKnABpIKCvr4/g4OBmLcvZEnbu3IkBAwYwJq+4uBhdunTB5MmT34rUkh07dmDYsGGqNoMR9u7dCzabjQMHDqjalBZRXV2NyMhIODs7w9HREeHh4U1OVdm+fTvYbDa+++67t+I9fZmSkhJERkYiMDAQurq66NmzJzZs2IDnz5/LlkI1NzdHly5dEBUV9Z+b7KwK1AGAGjUMcf36dUyaNAm6uroICAh447JubY2ysjLEx8cjLCwMgYGBMDExgYGBAXx9fWWjBHl5eao2861h9uzZ6NKli0JXBnodLwcC9+7dU4reR48eQVdXl9ERtby8PLi6uuLjjz9u83UyLCwMH374oarNaDHh4eEwNDTEuXPnVG1KsxEIBAgPD4eFhQW6dOnSrB7rqqoqTJ8+He3bt0dcXJyCLFUNFRUViImJwaRJk6Cvrw83NzesXLlStoStdGWyukuhqlEe6gBAjRqGyc7ORlhYGGxsbODk5ISw/8fenYdFWbV/AL9RtoGZYd8VBEUFwQ0jEPdQS+FVszFNx1xytNSpXBrLCq03Hbcal0rccUnFJcVSX8EVXCpcwV1wYcclFmWf+f7+6GJ+kqgsz8wz4PlcF38EzDk30wM+93POuW+lstHWJU9JSUFUVBTkcjkCAgLQpEmTZ1YJ2N7NuikvL0doaCjefvttXp8I3rt3D3K5XK+JQIsWLXD48GFOx0xPT4eXlxc++eQTTsfVN4VCgalTp/IdRp1pNBp89tlncHZ25m2bWX1lZ2cjIiIC1tbW9bpxvXHjBtq3b4+ePXsiMzOT4yj58fRNv6Wl5TM3/cA/DxeeLkn8slKojG6wBIBhdKSiogIxMTEIDQ2Fubk5pFJpg/0Hr6YKCgqqrBLY2tpCKBQiJCQECoUCMTExDaqCEt/y8/PRrl07fPHFF3yH8kwicP36dZ3NNXbsWM7OATzt7t278PDwwJw5czgfW18mTJiAiIgIvsOok4be3TclJUX7OxAWFlavczq//vorbGxsIJfLG/xDkqKiompv+v/9N+LOnTtV3j9Dbkr4KmAJAMPowbVr16BQKGBtbd0oDg3XRuUqgUwmg6+vL5o0aQIvLy9IpVKoVCokJiY2yj2vXElJSYGDgwNWr17NdygA/rmJlsvlEAgEkEgkOulIu3HjRnTt2pXzcYF/nro6OztDqVTqZHxde/vtt7F06VK+w6i1htzd9/z585BKpRAIBJBKpfU6F1NeXg6FQqH3al9ce95Nf3UrhCkpKZDJZDAzM4NEIsGVK1d4iJj5N5YAMIweFRQUIDIyEv7+/nB0dIRCocCdO3f4Dkuv8vPzERsbi4iICISFhcHGxgYikajKKkFj7XJaVydOnIBQKMTRo0f5DkWrMhGwtLSEVCrl9IluRkYGTExMUFBQwNmYT7t06RLs7Ozw448/6mR8XerVqxc2b97Mdxi10lC7+1bW8BeJRJDL5fU+lJ+Wlobg4GB06dKlQRaLePqmXygUvvCmH/jn9+zpxEkXDwuYumMJAMPwJD4+HhKJpNEeGq6piooKJCcnV1klMDIy0q4SREZGIjk5+ZVfJVi3bh3s7OwMbuvEnTt3IJPJtP/IcxVf69atsX//fk7Gqs758+dhY2ODVatW6WwOXWjfvr1O3xeuZWRkwN/fH0OGDGkQTZzUajViYmIQFBQEBwcHREREcPJA4tSpU3B2dsbEiRMbVF+KwsJCbN++He+88w4sLCzQvn17fPvtty/cAnju3DlIJBJt87PGXNK0IWMJAMPwLDMzE0qlEs2aNYO3tzeUSuUr/wQ8Ly+vyiqBlZWV9sBYREQEYmJi8OjRI77D1Lvp06fD19e3VrXF9eX27dtVEoFbt27Va7yJEydi5syZHEVXvVOnTsHKygq//PKLTufhkpubG86cOcN3GDXSkLr7lpaWIioqCj4+PvDw8IBKpeKsA+/mzZshEokazLaznJwcrF69GgMHDoS5uTk6dOiAb7755qUFAP69YtJYDjY3ViwBYBgDUVpaiujoaISGhmpvoi5cuMB3WAahcpUgMjISUqkUvr6+aNq0KXx9fausEjT2FRS1Wo1BgwbhjTfeMNiDg/9OBFJSUuo0zrZt29ClSxeOo3tWXFwcRCIR9u7dq/O5uCAQCBrEVoqG0t23sLAQKpUKzZs3h5+fH6KiojhrGldRUQGFQgF7e3scOXKEkzF15fbt21CpVAgNDYWpqSkCAgJeuL3nabGxsQgODoa9vT0iIiJeyYczDRFLABjGAF29ehVyuRxCoRABAQGIiooy2Bs+vmRlZSEmJgYRERHapMnKyqrKKoEhPimvr6KiIgQFBWHEiBEGnfCkpqbWKxHIycmBsbGxXlbDKrsvG/rWmuLiYhCRwa8QNoTuvpU16G1tbbWlPLn8fXr48CFCQ0PRvn17g93vn5ycDKVSiZCQEJibmyM0NBQqlapGT+4rt0q99tprcHR0RERERIM638GwBIBhDFp+fj4iIyPRrl07ODk5QaFQ4O7du3yHZZDKy8ufu0ogk8kQFRXVaFYJcnNz0apVK8ydO5fvUF7q6URAJpMhPT29xq/t3Lkztm3bpsPo/t+uXbsgFotx7NgxvcxXF+np6WjSpIlBb6cx9O6+T5ei1FXzqWvXrqFNmzaQSCS8NfKrTkVFBeLj46FQKODt7Q1LS0uEhYUhKiqqxjfvarUa0dHRVbZKvSoV7RoblgAwTAPx9KHhsLAwxMbGNoqbWV3KzMxETEwMFAqF9imXs7MzwsLCEBERgdjY2Ab7j9eVK1dga2uLDRs28B1KjVy5ckVbEaSmicCXX36J0aNH6yG6f2zcuBHW1tb1qu+uS5cuXYKtrS3fYTyXIXf3vXnzprYUpS5r0O/ZswdisRj//e9/DeLvc3FxMWJiYiCTyeDs7AwHBwdIpVLExMSgtLS0xuNUnpHw9vaGl5cXIiMj2ap0A8cSAIZpYDIyMqBUKuHm5obWrVtDqVSyPZc1VF5ejsTERKhUKkilUnh6esLY2PiZVYKG4vjx4xAKhYiNjeU7lBq7fPlylUQgIyPjud976tQpODg46PWJ99q1a2FjY4OzZ8/qbc6aOnr0KFq1asV3GM8w5O6+Z8+erVKKUpcN7JYtWwaxWMz7eZJHjx4hOjoaUqkUIpEIXl5ekMvliI+Pr/XvUlFREVQqFdzc3NChQwds377doFegmJpjCQDDNFBPHxoWiUSQyWS4dOkS32E1OBkZGVVWCczMzODi4lJllcCQyxdu3boV1tbWDe7/fU0SgYqKCtjb2+v9ibxKpYKDgwMuX76s13lfZteuXXj99df5DqMKQ+3u+++KNC9KNOtLo9EgIiICTk5OvHW3vX37NpYtW4Y+ffrAxMQEr7/+OubPn1/npmWPHz/G4sWL4eTkhMDAwFe2THVjxhIAhmkEzp07B5lMBktLS3ZouJ7KysqqrBJ4eHjA2NgYAQEBkMvliIqKQmpqKt9hVjF37ly0aNECWVlZfIdSa8nJyVUSgX8fQBw5ciQiIiL0HtfcuXPh5uZW5ypGurB69Wq89dZbfIehZWjdfas7mPr333/rdM6SkhKMGDECrVq10mt1pqf38/v5+cHMzAz9+vXDTz/9VKtzNv9WWRXJ2dkZwcHB7Ma/EWMJAMM0Inl5eYiMjISvry+cnZ2hUChw7949vsNq8DIyMhAdHQ25XI6QkBCYmppqVwmUSiXi4+N5b+4zadIkBAQEGNShw9pISkqq0jyoMpnZvHkzXnvtNV5imjVrFtzd3Q2mW7dSqcTIkSP5DgOAYXX3rdyf3qZNG3h6eurtYOqjR4/Qo0cPBAUFITc3V+fzPX78uMp+fnt7e0gkEkRFRdW74llBQQGUSmWVqkhM48YSAIZphDQaDWJjYyGRSCAQCCCRSNihYQ49fvwY8fHxUKlUkEgkcHBwgImJSZVVAn3fNJaVlaFv374YOHAgZ3XM+fDvRODq1aswMTHhbXVj+vTp8Pb2NoimRgqFAlOnTuU7DIPp7ltQUKDdn96+fXtOa/i/zO3bt+Hj44MhQ4boNNlISUmpUp/f19cXCoWiTvv5q5Ofnw+lUgkbGxt24/+KYQkAwzRy6enpiIiIgIODA9q2bQuVSoXCwkK+w2p0/r1KYGJiAhcXF0gkEqhUKsTHx9eq6kZd5OXlwc/PD1OmTNHpPPpw6dIlbSLg4uKCpUuX8hKHRqOBTCaDn58fHjx4wEsMlSZMmMDLdqinGUJ335ycHERERFS5adXnw43ExEQ4Oztj2rRpnL8HT2/t8fHxgUAg0Nbn53I198GDB4iIiIC1tTVCQkJw+PBhzsZmGgaWADDMK6KkpER7aFgsFkMmkyEpKYnvsBqtwsJCxMfHQ6lUIiwsDPb29rC0tERISAjkcjmio6N1sm/63r17cHNzg1Kp5HxsPly8eBF+fn4wNjaGXC5Hdna23mNQq9V477330LFjR14rbr399tu8JUIA/919U1NTIZfLtfXrT506pfcYjh8/DisrK6hUKs7GfPDggbZqj7W1NZycnCCVShEdHc35w5rKBmiVTRPPnDnD6fhMw8ESAIZ5BSUmJkImk2mb4URHRzfobSMNRUpKCqKioiCXyxEQEIAmTZo8s0rAxeHt5ORk2NraYu3atRxEzb9z585BJBJh6NChEIlEUCgUeu+GW1FRAYlEgqCgIL2toJ08eRKnT5/G9evXcf/+ffTs2RObN2/Wy9z/xmd334sXL1Yp5clXdaajR4/CysoK69atq/dYlV14K7f2BAQEICIiAomJiTpZzcjJyYFCodAmT3/++SfnczANC0sAGOYVlpeXB5VKBS8vL7i4uEChUCAtLY3vsF4ZBQUFVVYJ7OzsIBQKq6wS3L9/v05jnz59GiKRCL/++ivHUfPDy8sL+/fvx4ULFyCRSLSJgD6fyJeWlmLAgAHo06ePXg6Zjh49Gk2aNIGxsTGICEQEgUAAZ2dn+Pj4oFevXnrpW8FXd9/KUp6V50H47IL+22+/QSwWIzo6uk6vz8/Px+7duyGTydC8eXNYWVlh2LBhiIqKqvPveE1kZ2dXufFPTEzU2VxMw8ISAIZhoFarqz00zOhfdasEXl5ekEqlUKlUSExMrPG+471790IkEuHkyZM6jlr3ZsyYgXHjxmn/+/Tp09o67/pMBIqKitCrVy/069dP55WffvnlFwiFQu3N/78/HB0dUVFRodMY9N3dV6PRICYmBl27doVYLK5SEYov27dvh1AorFUyrdFocO7cOcyfPx89e/aEiYkJ/P39MXPmTBw+fFjnZZrv3r0LuVwOCwsLhIWF4dy5czqdj2l4WALAMEwVKSkpUCgUsLe3h4+PD1QqVYMtLdkY5OfnIzY2FhEREQgLC4ONjQ1EIhFCQkKgUCgQExPzwsOpkZGRsLOzq3NDIEPxxx9/wMbG5pmD1KdOnaqSCOi67jvwTxWobt26YciQIdVunePqpvz+/fto0qRJtTf/FhYWWLhwISfzVEff3X3LysoQFRWFdu3awdnZGREREfUubcmFzZs3QywW49ChQy/93ocPHyI6OhoymQzNmjXTPnWPjIzU2+rF7du3IZfLtQ9yGvrvPaM7LAFgGKZalYeGK5/EyWQyg+uM+iqqqKhAcnIyoqKiIJPJ4OvrCyMjoxeuEnz99ddo1qxZg+4JodFo4OHhgYMHD1b79ZMnT+o1EcjLy0NAQABGjRpV5b1OS0vD66+/ztk5AW9v72oTAFNT03qfg3heoqLP7r4lJSWIiopCq1at0LJlS6hUKoPpvL1y5UqIxWIcOXKk2q+r1WokJiZq9/KbmJhoy3TGxsbqvOrX01JTUyGTybTnJK5fv663uZmGiSUADMO8VOWhYUtLS4SGhrJDwwYmLy+vyiqBtbU1xGIxQkNDERERgZiYGEyYMAF+fn68VrGpr2nTpmH8+PEv/J6EhASEhYXB1tZW50+Rc3Nz4evri7Fjx0Kj0SA1NRUuLi5o2rQpZ1ViZsyYAVNT0yo3/yYmJi99H2pi/vz5mD9/fpXP6au7b35+PlQqFVxcXNCpUydERUXpfDtTbXz//fews7N75rBsbm6u9im/i4sL7OzsIJFIEBkZWa8OvHV15coVjBw5EgKBABMnTsTt27f1HgPTMLEEgGGYGvv777+hUqng6ekJV1dXRERE6KUDJlM7lasEkZGRkEql8PX1RdOmTSEWi+Hk5IQVK1YgOTm5wTWGO336NOzs7Gq0fzo+Ph5vvPEG7OzsdJoIZGRkoGXLlpgyZQrs7e21B3adnJw4SZJjY2MhEAieSQCuXLlSr3EfPnwICwsLGBsbY/369QD00903OztbW4bSUBtPrVq1Cra2trhw4QIqKiq0T/lDQkJgbGyMgIAA7VN+vh6E3Lx5E1KpFBYWFpg8eXKDXt1j+MESAIZhao0dGm54srOzsXPnTri7u8PBwQECgeCZVQJ97J+vj8ptQP/73/9q/Jr4+Hj06dNHmwjo4sY2NjYW1tbWVar1WFpaYuvWrfUeu7i4uMoKgJGREbp161bvcT/++GNtYmFqaopVq1bBx8cHw4cP18nWlVu3blU5lPrHH39wPgcXdu3aBUtLS0yfPh1DhgyBlZUVXFxcMHbsWGzfvp33FbS7d+9W2eqTkpLCazxMw8USAIZh6uXmzZtQKBSws7ND586dERkZiSdPnvAdFvMceXl56NChAyZMmFDtKoGvry+kUikiIyMNcpXg008/xYQJE2r9uvj4ePTu3Rv29vacJgJnz56FSCRC06ZNn9mn7+Pjw8kc3bp1q1IG9Pfff6/XeGlpac9sKzIzM8OECRM472x7/vz5KjX8DfFQ6sOHD7Fjxw4MGDAATZo0gampKfr27YuFCxfiwoULBvE7UFnH38LCAhKJROdnM5jGjyUADMNwori4GFFRUejYsSOsrKwgk8nqvU2B0Y20tDS4u7tj7ty5VT6fmZmJmJgYREREIDQ0FObm5nByckJYWBgiIiIQGxvLe3J36tQp2Nvb13nrRXx8PHr16gV7e3solcp6/Tzx8fGwsLCAkZFRtQd1zczMcPTo0TqPX2nx4sWwsLAAEcHNza3eN+nDhw+HmZlZlViNjIxgYWGBCxcu1Dte4P9r+ItEIsjlcoPqL1JeXv7M4V1XV1eYmppi2rRpKCgo4DtErezsbHzyySewsLDAe++9xw73MpxhCQDDMJx7utNw5aFhQzrgx/xzeNDBwQHLli177vdU3iipVCpIpVJ4enrC2NgYvr6+kMlkiIqK0ksjqqdpNBq4u7vXe8tZfHw8evbsCQcHhxolAtXVwf/mm29gZmYGc3PzahOAJk2aoE+fPvWKEwAuXboEExMTCASCeh8urhyrunibNm0KOzs73Llzp05jq9VqxMTE4PXXX4eDgwMiIiL03rH5eVJSUhAZGQmJRAIrKys4OTlpD+/+/vvvEIvFBtU5u6SkBEqlEmKxGEOHDtX77xnT+LEEgGEYncnJyYFSqYSHhwe8vLygVCp12vWSqZ0LFy7A1tYWq1evrvFrMjIyEBMTA4VCgZCQEJibm8PFxaXKKoGuu+ROmzYNY8aM4WSs2NhYBAYGahOB6mI/fvw4mjRpUu2e/vv372PGjBnPTQRMTEzqffOm0WhgY2MDc3Pzem9d6tWrV5WzCtUlAcHBwVVWGX744Qdcu3btuWOWlpYiKioKPj4+8PDwgEql4n2l6OlqPe7u7tqHEUqlEomJidptPUlJSbCzs8P333/Pa7xPi4mJQcuWLREQEID4+Hi+w2EaKSMAIIZhGB3SaDR05MgRWrp0KcXFxVF4eDh9/PHHFBISwndor7zTp0/Tm2++SatXr6Zhw4bV+vXl5eV06dIlSkhIoLNnz9KJEycoIyODOnToQCEhIRQQEEDdu3cnT09PzmJOTk6moKAgysrKIpFIxMmYcXFxNHv2bLpz5w5NmzaN5HI5CQQCIiIKDg6mP/74g0xMTGjbtm00ZMiQZ15///59WrhwIS1fvpyIiEpLS4mIyMTEhN59913atGlTtfNqNBpKTU2l5ORkunr1KmVkZFBaWhrl5uZSSUkJFRQUEBFRbm4umZqakpeXFzk5OZGLiwu5ublR27ZtqV27dtSmTRsyNTV94c94/Phx6tevH5WVlVX5fNOmTcnY2JgcHR1p+vTp9MEHH5ClpSURES1YsIBmzZpFUqmUNm7cWOV1jx8/prVr19LixYvJ2tqaZs6cSSNGjCATE5OXvd2cKyoqovj4eIqLi6PY2Fi6fPkyBQQEUN++falv374UHBz8TFz37t2joKAgmjp1Kn3++ed6j/nfHj16RFKplC5evEjz5s0jqVRKRkZGepu/oqKCUlJSKCkpia5du0YZGRmUkZFBOTk5VFJSQo8fP9Z+r1AoJHNzc3J0dCQ3Nzfttejn50fe3t5kbGyst7gNRVpaGl25coWSkpIoLS2NMjIyKCsri0pKSigvL0/7fWZmZmRhYUH29vbk7OxMzZo1o1atWpGfnx/5+vqShYWFXuJlCQDDMHp18+ZNWrt2La1evZo8PT1JJpPRqFGj9PZHj3lWXFwcvf3227Rt2zYaMGBAvcfLzMyks2fP0smTJykhIYESExPJ1taWAgICqFu3bhQSEkJdunQhc3PzOs/RpUsXmjx5Mo0dO7be8T4tLi6OvvjiC7p79y5NmzaNunTpQgMGDNDeNJuZmdHOnTspLCys2tfn5ubSggULaMWKFdSkSRMqKSkhExMTunPnDrm6uhIR0dmzZyk2NpZOnDhBCQkJVFFRQT4+PuTj40MtWrQgZ2dncnFxISsrKzIyMiJra2v69ddfqVWrVmRjY0M5OTmUmZlJ6enpdPXqVbp8+TKVlJTQa6+9Rj179qQ+ffpQjx49qtyEAaCOHTtSUlISVf6zb25uTuXl5dS7d29SKBT0xhtvVLnhnD9/Pn3zzTdUUlJCxsbGlJKSQu7u7vTgwQNasWIFLV++nHx8fEihUFBYWJheb1bLy8vpzJkzdOTIETpy5AidOXOG3N3dKTQ0lPr27Uu9e/cmGxub576+oKCAunfvTsHBwbRy5Uq9xf08ly5doiFDhlDnzp1p3bp1nCW2L6LRaOiPP/7QXounT5+mJk2aaK9FDw8P7bUoEomoadOmJBaLqbCwkCoqKqiwsJCys7MpOzub7t27R1euXKErV66QWq2moKAg6tmzJ4WGhlJwcDBWk+UPAAAgAElEQVQ1adJE5z+Pvt26dYv2799PJ06coPj4ePr777/J29ub2rVrRy1atCA3NzdycnIiOzs77e9xSUmJ9uPBgweUlpZGWVlZdP36dUpKSqJHjx5R+/btqWfPntS7d2/q27ev9kEE5/hcfmAY5tVVeWi4Q4cOsLa2hlwuR2pqKt9hvbJ+/fVXiESiave619eTJ08QHx8PlUoFiUQCR0dHmJiYICAgAHK5HFFRUbXed75ixQp0796d81iBf7bc7NmzB506dYKHhweaNGnyTBfe/fv3v3CMu3fvYty4cTAxMYGRkRHGjBmDTz/9FB4eHrC2tsawYcOwYsUKJCUl1eh8zMsOPaempiIqKgrjxo1Ds2bNYGdnh7Fjx+LIkSPQaDTYvn27du9/ZQnYzz///LmHc7/77rsqB4XNzMwwduxYbSlPfdfwV6vV2qpV1e3jr83fjrKyMvTt2xcDBw40iLNJW7duhVgsRkREhF4qDv3555/46KOP4OrqCnt7e4wcORKrVq3CtWvX6j2/RqPB9evXsWbNGkilUjg6OsLZ2RmTJk3C6dOnOfoJ+HP37l3MmTMH/v7+MDc3x4ABA7Bw4UKcOXOGk/K5WVlZ2LFjB6ZOnYq2bdvCwsICQ4cOxa5duzjvOcESAIZheJeYmAipVAozMzN2aJhHGzduhLW1Nf766y+dz5WRkYHo6GjI5XKEhITA1NRUe5ZAqVQiPj4eJSUlz319Xl4eBALBC/em11dcXNwz1XIqP8zNzXHs2LEXvl6tVmPFihVwcHCAkZERpFIpDh48qJM6+0/TaDQ4c+YMZsyYAScnJ7Ru3RrW1tba0qQbN2584Xs7e/bsan9uY2NjDB06FBcvXtRp/JWePrhra2sLkUhU7T7+2po0aRI6d+6Mx48fcxxx7VRUVGhLKB88eFCnc5WVlWH9+vXo3LkzrK2t8dFHH+HIkSM6/ztbUVGB48ePY+rUqbCxsUGHDh2wZs2aGjXzMyRHjx5FeHg4zMzMMGTIEERHR6OwsFDn8165cgXffvstWrVqBRcXF0RERODBgwecjM0SAIZhDEZ2djaUSiXc3d3RsmVLKJVKzv7YMTWzYsUK2Nvb4/Lly3qd9/Hjx1VWCezt7bVPmuVyOaKjo5GTk1PlNe+++y5mz56ts5hee+21Z57+V34YGRlBIBA8t6FVdHQ0fH194enpiRUrVuD8+fO8lMUtLS2FTCaDnZ0d7O3toVKpXnjzNWvWrGc6D1d+CAQCzJkzR2exZmZmag/uNm/e/JmDu1z0KFiwYAE8PDyQmZnJQcR19+DBA4SGhqJ9+/Y6bealVquxbt06eHl5wdfXF2vXruXtgHZRURE2bNgAf39/tGjRAqtXrzb4Bz3Hjx9Hjx49YGdnh7lz5/J23Wg0Ghw6dAgDBw6ESCTCF198Ue/u5iwBYBjG4FSWE6ysRS+RSHDq1Cm+w3plfPvtt3B1deW9y2hKSgqioqIgl8sREBCAJk2awMXFBRKJBCqVCkuWLIGrq6tObiLi4uKqbe717yTA0tIS586d077u2rVr6NOnD5o1a2YwTzqLi4sBAL///js6deoEHx8fHDlypMr3aDQaTJky5bklTSs/hEKh9sl5fSt6ZWRkYMuWLfjggw/g5eUFU1NT9OjRA3PmzEF8fDzn793BgwdhZWWltxWM57l69So8PDwwYsQInd6Mnz17FoGBgWjVqhW2bt3KeZO3ulKr1dixYwfatm2LgIAA/Pnnn3yH9Izs7GxIpVJYWVlh3rx5BtUb4q+//kL//v3h7OyMzZs313kljCUADMMYtOvXr0OhUMDGxgYBAQGIjIzUeZlJBvjss8/g7u6Ou3fv8h2KVmFhIeLj46FUKhEWFgY7OzsYGRmhXbt22lWC3NzcOo8/efJk7RP9WbNmwcnJCU2aNEHTpk1haWkJoVD4TFJgZGQEsViMpKQkLFu2DEKhENOnT9fL9oC6qKiowPLly2FtbY2pU6eipKQEGo0GEyZMeO6T/6c/mjZtCqVSiREjRqBHjx61mjszMxO//PILZDIZ2rRpA2NjYwQGBkKhUOB///ufTm+Gr127Bmtra2zfvl1nc9TEX3/9BXt7e52upFRUVOCbb76BUCjE3LlztQmgoSktLcV3330HoVCIr776ivM97nW1d+9e2NvbY8SIEbyvFL3Irl270Lx5c/znP/+p00o5qwLEMEyDUFhYSFu3bqUff/yRsrOzaezYsTRp0iRq0aIF36E1SgDoo48+oiNHjtCJEyfIycmJ75CqNWXKFDp79iwFBgbSyZMn6fz58+Tk5KStNtStWzfq1KnTS6uQlJWVkVAoJAA0btw4WrhwIVlZWVFFRQWlpaVRamoq3b59m27fvk2XL1+mmzdvUnp6urZMp6mpKTk6OtKOHTsoKChIHz96vdy9e5fee+89KioqotatW1N0dDQR/fNzVFYPKisro4qKCjI3Nyc7OztycnIiY2NjSk5OJrVaTQDo8ePHzy37mZubS8ePH6eEhAQ6efIkXbx4UVsetlu3bhQaGvrCSj1cKSwspKCgIHrnnXdo7ty5Op/veY4fP06DBg2ir7/+mqZNm6aTOR49ekRDhw6l+/fv07Zt28jPz08n83Dp6tWrNHz4cLKysqLdu3eTvb09L3FoNBqaMWMGbdiwgVatWkXvvPMOL3HURn5+PslkMjp16hTt3r2bXnvttZq/mNt8hGEYRvf+fWg4JiZGL9UzXjVqtRrvvfceunTpgr///pvvcKqVmpoKc3Nz7ROwgoKCKqsEtra2EAqFCAkJgUKhQExMTLVPy86cOaPd/mJmZgZra2tERUXVaH4/Pz906tSpRt9vSMrKyhAUFARTU1P06dMHn3zyCZYsWYItW7bg+PHjuH79uvapfHZ2NsLCwqpsETI3N69S2SUnJ0d7sDsgIKBK1+jo6Gg8evRI7z+jWq3GwIEDMWjQIF63wMTExEAoFGLdunU6m+POnTvw8fHBu+++2+BWSYuLizFy5Ei0bt2al2pwxcXFGDp0KNq1a4fbt2/rff76WrZsGcRi8Uurkz2NJQAMwzRYWVlZUCqVaN68Oby9vaFUKvHw4UO+w2pUysrKMGTIEAQGBtb70Jmu9OjRAytWrKj2axUVFUhOTkZUVBRkMhl8fX1hZGQELy8vSKVSqFQqJCYm4vvvv4eFhcUzh16Dg4Nx/fr1asd+9OgR/Pz8IJVKDWKvf10plUrY29s/95BydHQ0RCIRTE1Nn6mE9MEHH+DDDz+Ej48PmjZtii5dumDGjBn47bff6t21mAvz5s1D69ateY1l06ZNEAqF2L17t87myMjIQIsWLfDxxx8bzF7/2tJoNJgxYwaaN2/+3PK0ulBeXo6BAweie/fuvCSpXNm1axdEIlGNK0qxBIBhmAavoqJCe2hYIBBAKpXiwoULfIfVaJSVlWHw4MEICAgwyH8gN2/eDF9f3xqvAt2/fx8xMTH44osv0KtXL1haWqJ58+bV7nk3NjaGqakpvvrqqyrlO0tKShASEoJ33nnH4CuZ1MQ333wDd3d3ZGRkaD93+/ZtdOvW7bmlUIkIrVq1wvTp07Fv3z6DSxCPHz8OoVBY5ZC2vq1fvx7W1tY4ceKEzubIz8+Hv78/Jk2apLM59GnKlCnw9fXV26rjmDFjEBAQYFAHfetqx44dEIvFSExMfOn3sgSAYZhG5dq1a5DL5RAKhdpDw4Z6CK4hKSsrw6BBgxAUFGQQT3afVlpaCmdnZxw+fLhOry8vL4etre0LD79aWFjA3d1dO8enn36KoKCgF9bTb2jGjh2LPn36oKKiAitXroS5ubm2edjzPgQCgUEmQLm5uXBzc8Pq1at5i2Hv3r06a673tJEjR2LAgAEG+f+hLtRqNcLDwzFs2DCdz7V27Vq4ubk9U2K4IVuyZAm8vLxe+neaJQAMwzRKBQUFiIyMhJ+fH5ycnKBQKGrdbZapqrS0FIMGDUJwcLDBJQFff/01Bg8eXKfXpqenP7fe/9MfRkZGICK8+eabEIvFvJdJ5dqTJ0/g5eUFNze3l74XlR+mpqa8PmGvjlqtRv/+/TF8+HDeYjh69CgsLS2xa9cunc6zfft2ODs7N6obWOCfPglubm7YuHGjzua4ffs2xGLxMyVxGzqNRoOBAwdizJgxL/w+lgAwDNPoxcfHQyKRwMzMDGFhYYiNjWWHhuuotLQU//nPfxAcHGxQS+aZmZkwNzev0wHCHTt2PFP/3sTEBEKhEEKhEMbGxiAiWFlZoW3bthCLxQgLCzOYsoVcKS8vx+zZs2Fubo7Ro0ejX79+8Pf3h729vfY9MDc3h0gkglgshkAggLGxMVQqFd+hVzF37ly0bt2at+vzwoULsLGx0emBX+Cfg6seHh7YunWrTufhy65du+Dm5qaz8rAjRozAhAkTdDI239LT0yEUCnH27Nnnfg9LABiGeWVkZmZCqVSiWbNmaN26NZRKpUHuaTd0lUlA165dDSoJGDZsGBQKRa1fN23aNJiamsLNzQ3BwcEYNWoU5syZg/Xr1yMuLg43btzQbiPbvn07WrVqVeU8gL4kJSWhV69eEAgEcHV1xVdffaWTbR+DBg2q9n3MycnBpUuXcODAAWzYsAHz5s3D5MmT8eOPP3IeQ10dPXoUQqGQtzNAN27cgIODA5YtW6bzuZYvX47XXnuNk4cZHh4eVRLg+fPnAwDmz59f5fNPH87Vx/XYtWtXLFmyhNMxgX9iFwqFnNT5r8t7B/zTbNDJyUln3cxnz56NsLCw536dJQAMw7xySktLER0djdDQUIhEIshkMt67gzY0paWlCA8PR0hIiMEkAQkJCbCxsal1E67KzrY1ERoaiu+//762odVbdnY27O3tMXnyZBQUFODs2bNwdnbGzJkzOZ/r6NGjcHJyanCVjbKzs+Hi4oL169fzMn9+fj58fX3rlITWRYcOHbB582bOxtu8eTOI6Jkb0uXLlz/zOX1dj9u3b0ebNm04X7H95JNPMH78eM7Gq817p1arERERgaCgIBgbG+ssAcjOzoaZmdlzmzmyBIBhmFfalStXIJfLYWlpiYCAAERFRTW4Gx++lJaWIiwszKCSgK5du+KHH37Qydg5OTkwMTHhZb+1QqGAvb19lWtz+fLlMDExqVf34+poNBp4enrWuJygIVCr1ejbty/Gjh3L2/zh4eHo37+/Xg7jXr58GSKRiPN6/++++y6aNm2KU6dOAfjnSXnPnj2f2e6mr+uxpKQENjY2nJ8zcXZ25vxwdk3fu5s3b2LRokWoqKiAmZmZzhIAAAgLC8PixYur/dqLWyMyDMM0cj4+PrR06VLKzMwkmUxGCxcuJHd3d5o1axalpaXxHZ5BMzU1pV27dpGNjQ0NGDCAHj9+zHdINHPmTFqyZAmVl5dzPvbJkyfJz8+PHB0dORnv8ePHZGRkpP2YNWsWERH997//1X5OpVIREdGBAweoW7duVbruhoaGUnl5OcXFxXESTyUjIyPq06cPJSQkcDquLs2fP5+ysrJoxYoVvMz/1Vdf0Y0bN2jbtm3UtGlTnc8XHx9PXbt2JYFAwOm4P//8Mzk7O5NUKqUHDx6QTCajtWvXartDV9LX9WhmZkbdunXj9Fq8efMm5efnc96xu6bvXatWrWjGjBl6uU7eeOON5753LAFgGIYhIrFYTDKZjJKSkmjTpk2UmppK3t7eFB4eTnFxcQSA7xANkqmpKe3cuZOsrKzorbfe4j0JGDRoEInFYtq+fTvnY587d44CAwM5G08oFFJJSQnZ2dnRnDlzSKlUEhHRl19+SZ999hlt2bKFPvnkEyL656bF3d29yusr//vWrVucxVTp9ddfp7Nnz3I+ri6cP3+elEolbd68mSwsLPQ+/+7du+nHH3+k3bt3k7W1tV7m5PparGRjY0MbNmyg1NRUat++PU2cOJFatmz5zPfp83rk+lo8d+4cde7cmUxNTTkbk6jm750+BQYG0rlz56r9GksAGIZhnmJkZEShoaEUHR1NqampFBAQQO+99x75+PjQggULKC8vj+8QDY6ZmRnt2LGDBAIBhYeH85oEGBkZ0aeffkoLFy7kPGnLzc0lZ2dnTsc0MzOjsWPH0urVq0mtVhMRUVlZGf322280dOhQIiICQMXFxWRubl7ltZX//eTJE05jIiJydnam3NxczsflWmlpKb3//vv05ZdfUocOHfQ+/5UrV2jMmDG0detW8vX11du8urgWK4WGhtLgwYMpKyuLbG1tn/m6vq9Hrq/F3NxccnFx4Wy8p73svdM3FxcXysnJqfZrLAFgGIZ5DldXV5ozZw6lp6fTt99+S3FxceTu7k4TJ06kpKQkvsMzKAKBgPbu3UsCgYD69u1Lf//9N2+xVC7BHzhwgNNxCwsLSSQScTomEdHEiRMpMzOT9u3bR0REu3btovDwcDIzMyOif5IagUBAJSUlVV5X+d+WlpacxyQWi6mgoIDzcbk2a9YsEolENGPGDL3P/eTJE5JIJDR9+nR666239Dr348ePdXItEhElJSVRXl4etWvXjj744AO6f/9+la/r+3rk+losLCwkoVDI2XhPe9l7p29isZhKS0uptLT0ma+xBIBhGOYlTE1NSSKRUGxsLB07doyIiIKDg6lbt260Y8cOnew3b4gEAgHt2bOHmjdvTt26daPMzExe4jAzM6OZM2fS119/zekqgJWVlU5uilu1akWhoaH0888/ExHRqlWraOLEiVW+x9vb+5kzKffu3dO+nmt5eXlkY2PD+bhcio+Pp7Vr19KGDRv0sp/63z788ENydnamL7/8Uu9zi8Viys/P53zckpISmjRpEq1Zs4Y2b95MeXl5NGHChGe+T5/XI9fXoq5+j2v63ulTXl4eCQQC7cOEp7EEgGEYphY6d+5MkZGRlJGRQaNHj6aIiAjy8PCgWbNmUXp6Ot/h8c7U1JS2bt1KQUFB1L17d7p9+zYvcXz44YeUlZVFBw8e5GxMJycnnSU1kyZNotjYWIqJiSGBQECenp5Vvv7mm29SQkJClWTz8OHDZGxsTKGhoZzHk5WVxdlhZ13Iz88nqVRKS5YsIW9vb73P/9NPP9HRo0f1duj335ycnCgrK4vzcWfOnEmTJk0iLy8v6tixI82dO5f27t1La9asqfJ9+rweub4WHR0ddfJ7XNP3Tp+ysrLIycmp+i/qrPYQwzDMK0CtViM2NhYSiQQCgQASiQSxsbF8h8U7jUaDTz/9FO7u7rhx4wYvMSxZsgQBAQGc1RCPiYmBn58fJ2P9W3l5OVxdXSESiRATE/PM17OysmBnZ4epU6eisLAQ58+fh4uLC6ZPn66TeEaPHo05c+boZGwuSKVS9O/fn5eO3ufPn4dIJEJ8fLze5660Zs0a9OnTh9Mx9+/fjxEjRlT5nFqtRvfu3SEUCnHr1i3t5/V5Pb755pv46aefOBsvNTUVZmZm2uZ+XKjNe/c0XZcBXbhwISQSSbVfYwkAwzAMR9LT0xEREQF7e3u0bdsWKpWq1k2pGhulUgknJydeOrMWFxfD1dUVBw4c4GS8Bw8ewMTEBBkZGZyM928RERFwd3d/bh35ixcvomfPnjA3N4ezszNmz56tk5rzarUazZs3x+HDhzkfmwt79uyBvb09srKy9D73o0eP4OnpyUszuKfduHEDFhYWnP19cXJy0nasHThwoPbzb7zxRpVutk83z9LH9VhUVASxWIykpCROx23WrBlnD2rq8t79u3swEeHnn3/mJJ6n9evXD0uXLq32aywBYBiG4VhJSQmio6MREhICsVgMmUyG5ORkvsPizYIFC2BjY6NtkKNPS5YsQWBgIGdPigcMGID58+dzMpahOnjwINzc3PTS0Kq2cnNz4eTkhO3bt+t9bo1Gg8GDB+M///kPLysP/9alSxesXbuW7zB0auPGjfD39+d83M8++wyjRo3ifFxDkpaWBjMzM2RmZlb7dXYGgGEYhmNmZmYkkUgoISGBjhw5QkT/1GOuPDRcUVHBc4T69dlnn9H8+fOpf//+nDetepmPPvqIsrOzaefOnZyMJ5PJ6Oeff36mAkpj8v3339MHH3zAy972l5k4cSKFhobSsGHD9D73ggULKDk5mTZu3EhGRkZ6n//fZDIZLV26VFs+trHRaDSkUqlo0qRJnI89fvx4+vXXX7UHlxujH374gcLCwp5f8lTPCQnDMMwr6e+//4ZKpYKnpydcXV0RERGBnJwcvsPSq19++QVCoRC//vqrXuddv349vLy8UFpaWu+xNBoNXn/9dSiVSg4iMzyHDh2CnZ0d8vLy+A7lGfv27YO9vT1yc3P1PvexY8cgFApx9uxZvc/9PGVlZfD29sa6dev4DkUnNm/eDE9PT5SUlOhk/HHjxjXaVYDbt2/DwsLihSvPLAFgGIbRo1f90PC+ffsgFAoRFRWltznVajU6duyIFStWcDLeyZMnIRKJcOXKFU7GMxR5eXnw9PTEjz/+yHcoz8jPz4ebmxu2bNmi97kzMjLg4uKC1atX633ul4mJiYGdnR3S0tL4DoVTmZmZcHR0xI4dO3Q2R3p6OqytrbFv3z6dzcEHtVqN3r1746OPPnrh9xkBrL89wzAMH1JSUmj16tW0Zs0a8vDwoIkTJ9LIkSN10tjJkMTFxdHQoUPp+++/p/Hjx+tlzt9//53Gjh1Lt27dIrFYXO/xvvrqK9q3bx+dPHmyUfz/AkDvvvsuFRcXU0xMjEFscXnapEmT6O7du5w3d3uZkpIS6tWrF/n7+9Pq1av1OndNTZw4ka5fv06HDh0iU1NTvsOpt/LycnrrrbfI3d2d1q1bp9O5tm/fTlOnTqWzZ89S8+bNdTqXvsydO5eio6Ppr7/+IgsLi+d/oz6yEYZhGOb5Kg8NBwcHw8rKCjKZDJcvX+Y7LJ06efIkbGxsMG/ePL3N2atXL85K7pWXl6N///548803UVZWxsmYfJo+fTpat26NBw8e8B3KM06fPg2hUIiUlBS9zz1u3Dh07dpVZ9tQuPDkyRMEBgbivffeg1qt5jucetFoNBg9ejS6dOmCx48f62VOuVwOHx8fg7z2a2vdunWwtbWt0b8fLAFgGIYxIImJiZDJZLC0tERoaCiio6NRXl7Od1g6cfnyZbi7u2P8+PF6+RkvXboECwsLXL9+nZPxCgsLERgYiEGDBqGoqIiTMfVNo9Fg1qxZcHV1xe3bt/kO5xmlpaXw9fXFDz/8oPe5Fy5cCHd39wZxVuf+/fto27YtxowZ02D/XpSXl2PChAnw9vbW63uuVqsxfPhwdO7cmZfSslxZv359rfpTsC1ADMMwBig3N5fWr19PK1eupPLycvrggw9o8uTJ5ODgwHdonMrMzKQBAwZQixYtaOvWrSQQCHQ63+TJkyktLY1iYmI4GS8/P5+GDBlCpaWltHv37ud33TRAJSUlJJPJ6PTp03Tw4EFq2bIl3yE9IyIigg4cOECnT5/Wa1WiQ4cO0dChQ+no0aPUpUsXvc1bHzk5OTRw4EBydHSkrVu3kpWVFd8h1VhBQQGNHDmS0tPTaf/+/c+vXKMjFRUVNHr0aDp58iQdOHCAfH199Tp/fQCgb7/9lpYuXUp79uyh7t271/iFDMMwjIGqPDQcFhYGc3PzRnlo+O+//0bPnj0RFBSE+/fv63SuR48ewcHBAb/99htnY5aUlOD999+Hs7MzDh06xNm4unT16lW0b98e3bt3N9gn3NeuXYOlpSXOnTun93ltbGx46TVQX4WFhRg0aBCaNWuG06dP8x1Ojfz1119o2bIlBg4ciPz8fL3Pf/PmTXz44YewsLDAG2+8AZFIhDVr1ug9jrrIyclB//794e3tXettoywBYBiGaSBu3rwJhUIBW1tbBAQEIDIyEk+ePOE7LE6UlJRg+PDh8PHxwZ07d3Q616pVq9CyZUsUFxdzOu7GjRshFovxwQcf6DyRqavi4mJ88803EAqF+PLLLw12u4harUa3bt3w+eef63XerKwseHp64uuvv9brvFy5ePEievfuDUdHRwiFQnz66ae83FTXRGFhIWbOnAlLS0ssWbJE783VEhMTIZVKIRAIIJVKcenSJQDA0aNH4ebmhgEDBuDWrVt6jamm1Go11qxZA3t7e4waNQoFBQW1HoMlAAzDMA1McXExoqKi0LFjR1hZWUEul/NyQJJrarUacrkcrq6uSExM1Ok8gYGBmDNnDudj37t3D0OHDoWdnR0WLVqkt4OML1NRUYFNmzahZcuWCAwM1On7y4WVK1fC29tbr2cr8vPz0alTJ4waNcogOv3WxsOHDyGXyyEQCCCXy5Gfn4/r16+jX79+cHV1xY8//mgwB5lLS0uxcuVKuLm5oXfv3notp6tWqxETE4PQ0FCIxWLI5fJqS6jm5eVBLpfD0tISM2bMMKhVsoMHDyIwMBBeXl71KmHKEgCGYZgGrPLQsIWFhfbQcEVFBd9h1UtkZCSsrKx02jDs4sWLEAqF2qd+XIuLi0NISAgcHR3xzTffICMjQyfzvExBQQF+/vlntG7dGq1atUJUVJTBV4p5+PAh7O3tsX//fr3NWVpair59+2LgwIEGuypSnfLyckRGRsLe3h5hYWHVPrHes2cPOnXqhObNm2PRokW8rU49fPgQS5Ysgbu7O/z9/bFz5069zZ2fn49ly5bB29sbHh4e+OGHH2r01PzixYsYPHgwhEIhPv74Y1y9elUP0T6rtLQU27ZtQ1BQEBwdHbFo0aJ6J8csAWAYhmkEsrOzoVQq4eHhgZYtW0KpVBrsNpSaOHjwIKysrBAREaGzOb7++mt06tRJp2U8Y2NjER4eDjMzMwwePBi//PKLzrdklJaW4sCBAxg/fjxEIhGCgoKwadOmBnNjK5PJMGTIEL3Np9FoIJVKERgYaDArNjURFxcHf39/tG3b9qXJkkajwd69exEaGgpzc3MMHz4cu3bt0vkWwqKiIuzZswejRtkDNacAACAASURBVI2CQCBA7969sWvXLr0loVeuXMHkyZMhEokQEhKCbdu21en34MKFC3j//fchEAjQs2dPrFy5EtnZ2TqI+P+p1WqcPHkS06ZNg5OTE9q2bQuVSsXZNcoSAIZhmEakukPDJ0+e5DusOrl06ZK2TKgubtJLS0vh7+8PpVLJ+dj/dvfuXXz77bfo2LEjzMzM0LdvX3z77beIj49HYWFhvcYuLS3F2bNnoVKp8Pbbb8Pa2hotWrTA9OnTcfHiRY5+Av1ITEyEUChEamqq3ub85JNP4O3tjdzcXL3NWR83b96ERCKBjY0NVCpVrW9ob9y4gdmzZ8PHxwcWFhYYOHAgFixYgNOnT9f7qXJxcTHOnDmDRYsWITw8HEKhEG3atMEXX3zBWfndl6nub+CpU6c4Gfvhw4dYtmwZevToARMTEwQHB+Pzzz/HgQMH8PDhw3qNXVFRgatXr2L16tWQSqVwcXGBg4MDJkyYgGPHjnG+LY2VAWUYhmmkbty4QevWraNVq1aRl5cXyWQykkqlOi+1yaXMzEwKDw8ne3t7io6O5ry04fnz56lHjx70xx9/6K30X2pqKh06dIhOnDhB8fHxlJGRQS1atCBfX19q0aIFOTs7k5ubGwmFQjI1NSVLS0sqLS2loqIiKi4upuzsbMrMzKT09HS6dOkS3blzh0QiEYWEhFCPHj2oT58+1LlzZ738LFzSaDQUEhJCAwYMoK+++kovcy5YsICWL19Op06dInd3d73MWVdPnjyhRYsW0eLFi2nkyJH03//+t95lga9evUqxsbEUHx9PCQkJdP/+ffL09KR27dqRu7s7OTs7k6urK1laWpK5uTkJBAIqLi6mkpISKioqoszMTMrKyqK0tDS6fPkypaamkp2dHXXr1o169OhBoaGh1K5dO47egRd78OABrVu3jn7++WciIvrwww9p/PjxZGdnp5P5cnJy6H//+x+dOHGCEhIS6MaNG+Ti4kLt2rWjFi1akKurK7m4uJCtrS0ZGxuTSCSiiooKKiwspLKyMsrNzaX09HTKysqi69ev05UrV6hp06b0+uuvU48ePah3794UEhKis/K3LAFgGIZp5B4/fky//PIL/fjjj3Tv3j0aPXo0ffrpp9SiRQu+Q6uRJ0+e0KhRo+jq1au0Z88eatu2Lafjz5o1iw4fPkwnT54kU1NTTseuiYcPH1JSUhJdvXqV0tPTKSMjg3JycujJkydUVFREarWamjZtqk0GnJycyNnZmZo1a0Y7duyg4OBgWrx4sd7j5tratWtp/vz5lJycTObm5jqfb+XKlTR79mw6duwY+fv763y+ugJAmzZtIoVCQT4+PqRSqah9+/Y6mSsnJ4eSk5O112JmZqb2WiwpKdFei+bm5tpr0cXFhZo1a0Y+Pj7k5+dHzs7OOomtOhqNhg4fPkxr1qyhmJgY6tatG02ePJnCw8P12jeCiKiwsJCuXLlCSUlJlJaWRhkZGZSdnU1PnjyhwsJC7XtXmQzY2tqSi4sLubm5kbe3N/n7+5Onpyc1adJEPwFzup7AMAzDGLTK0ndmZmYIDQ1FTExMg6h4otFooFQqYWVlhd27d3M6dllZGYKDgzFt2jROx9WHffv2oUWLFgZ/sPdlKvsz1KeqSW1ERUXVqmsqX/78808EBwejefPmiIqK4jscg5GZmQmlUgkvLy84OztDoVAYbMlOQ8USAIZhmFdQVlYWlEol3N3d0apVKyiVSjx48IDvsF7qt99+g7W1NRQKBac3vSkpKbC2tsaePXs4G1MfysvL4erqiri4OL5DqZfJkycjPDxcL3Pt2rULYrEYR48e1ct8dZGRkQGpVAqhUIiIiAjOe1Y0RJV7+yUSifYBRnR0tE4P8TdmbAsQwzDMK0ytVtP+/ftp2bJllJCQQBKJhKZNm0YdO3bkO7TnunHjBg0ePJi8vLxoy5YtnJ0L2LlzJ02cOJHOnz9v8PvBn/b555/TvXv3aMuWLXyHUidJSUkUHBxM58+fJ29vb53OdfDgQRo2bBht376d3nrrLZ3OVRfFxcW0bNkymjdvHvXv358WL17coK5FXUhPT6ctW7bQTz/9RBqNhkaOHEkffvgheXh48B1ag8YSAIZhGIaIiK5fv07r16+nyMhIatmypUEfGi4sLCSpVKo9F+Dj48PJuDKZjK5du0ZHjx7V+x7iukpNTSU/Pz/KyMggGxsbvsOpte7du1OfPn1o7ty5Op0nLi6O3n77bdqyZQuFh4frdK662LdvH3388cdkZ2dHS5cupa5du/IdEm9KS0spJiaGVq1aRQkJCRQeHk5SqZQGDBjQYH4vDR6/CxAMwzCMoSkoKEBkZCT8/f3h6OgIhUKBO3fu8B3WM9RqNT7//HPY2Njgt99+42TMoqIi+Pn5Yfbs2ZyMpy89evTAihUr+A6j1qKjo+Hm5qbz+vsJCQkQi8XYunWrTuepiytXrqB///5wdXVFZGRkgz/PUR9Xr16FQqGAg4MD2rRpA6VSaVBdeBsTlgAwDMMwzxUfH19lz60hHhresWMHrKysoFAoOGl2dfXqVVhbW2Pbtm0cRKcfGzZsQPv27fkOo1ZKS0u13Yl16fTp0xCLxVizZo1O56mthw8fQi6XQyAQQC6X16gzbWNUUFCANWvWIDg4GJaWlhgzZgwSEhL4DqvRY1uAGIZhmJfKysqijRs30ooVK0ggEND48eNpwoQJZGtry3doRER08+ZNGjZsGJmamtL27dvrXeL00KFDNHToUDp8+DAFBgZyE6QOFRcXk6urKx05coQ6derEdzg1snjxYtqyZQudPXtWZ6UPL168SH369KHvvvuOJk2apJM5aqu8vJzWr19Ps2fPpqCgIFq6dCl5eXnxHZZeaTQaOnbsGG3YsIF2795Nvr6+NG7cOHrvvfdILBbzHd6rge8MhGEYhmk4SktLER0djdDQUAgEAkilUly4cIHvsAD804VULpfD3t4e+/fvr/d4KpUKLi4uSEtL4yA63ZPJZJgyZQrfYdTIo0ePYGtrq9PqRdeuXYOTkxMWLVqkszlqKy4uDn5+fmjbti0OHDjAdzh6d/fuXW35ThcXF8jlcpw/f57vsF5JLAFgGIZh6uTq1auQy+UQCoUICAhAVFSUQZTk27lzJ6ytrSGXy+sdz6RJk9C5c2c8efKEo+h0548//oC1tTWKior4DuWlpk6disGDB+ts/Bs3bsDFxQXz5s3T2Ry1cePGDUgkEtja2kKlUnGyVa2hyM/PR1RUFEJDQ2Fubo6wsDBWvtMAsC1ADMMwTL0UFBTQtm3baNmyZfTw4UN6//336aOPPuK1fOGNGzdIIpGQtbU1bd26lVxdXes0Tnl5OfXr14+cnJxo69atZGRkxHGk3OrQoQPNmjWLRowYwXcoz3X9+nXq3LkznT9/nlq3bs35+Pfu3aMePXrQmDFjaM6cOZyPXxtPnjyhRYsW0eLFi2nkyJH03Xffkb29Pa8x6YNGo6FTp07Rpk2b6JdffqE2bdqQVCqlkSNHvhI/f0Ogp37DDMMwTGMlFotJJpNRcnIy7dixg1JTU6l169YUHh5OcXFxxMdzptatW9OZM2eoVatW1LlzZzp48GCdxjExMaGdO3fS2bNnaebMmRxHyb0xY8bQunXr+A7jhT777DOaMGGCTm7+MzIyqHfv3jRkyBBeb/41Gg1t3LiRWrZsSfHx8XT69GmKjIxs9De/9+7dowULFpC3tzcNGzaMzM3NKSEhgRITE+njjz9u9D9/Q8JWABiGYRjOZWZm0qZNm2j58uVkaWlJ48aNI5lMxkud+i1bttDkyZPp/fffpwULFpC5uXmtx0hLS6OQkBD66KOPaNasWTqIkhsPHz6kZs2aUXJyMrVs2ZLvcJ5x7NgxGjJkCN26dYvs7Ow4HTs3N5d69epFXbt2pdWrV/O2WvPnn3/Sxx9/TA8ePKB58+aRRCLhJQ59KSgooD179tCmTZsoISGBQkNDafTo0TRkyBAyNjbmOzzmOVgCwDAMw+hMWVkZ7d27l1atWkV//PEHjRgxgqZMmUL+/v56jePu3bsklUrp4cOHtGXLljp1Or58+TL17NmT5s+fTxMmTNBBlNyQSCTk6+ur88ZatQWAXnvtNRoxYgRNnz6d07Hz8vKoT58+1L59e1q3bp3Oqgq9SEZGBn3++ef066+/0vTp02nWrFl1SjYbgoqKCoqNjaUtW7bQ7t27qWPHjjR69GgaPnw4WVtb8x0eUwNsCxDDMAyjM6ampiSRSCg2NpaOHz9ORETBwcHUpUsX2rhxI5WXl+slDg8PDzp69CiNHj2aevToQQsWLCCNRlOrMdq1a0e///47zZgxg3bs2KGjSOtv3LhxtG7dOlKr1XyHUkV0dDQ9ePCApkyZwum4+fn51LdvX2rZsiWtWbNG7zf/xcXFtGDBAvL19SWif0rSzpkzp1He/P/111/0ySefULNmzejDDz8kDw8POnfuHJ06dYomTZrEbv4bELYCwDAMw+hVfn4+bd++nZYuXUqPHj2i999/nyZPnkzNmzfXy/x//vknjRw5kjw9PWnDhg21PiAcFxdHb7/9Nu3cuZP69eunoyjrTqPRkKenJ61atYr69+/PdzhERKRWq8nPz48UCgWNGTOGs3GLiorozTffJDs7O4qOjiYTExPOxq6Jffv2kVwuJwcHB1q6dCkFBwfrdX59SEtLo/9r796joqzTOIB/ZxiY4U6iYFwSQRBE8S4EKoUQKqCYgomISYquGpi7hltugtnG0a3AThp4SVBRQS3A0kRIEBGviInJpiByEUS5iAHDMPPbPzrOSuCFYS4gz+cczknmfZ/neT2e0++Z93dJTEzEd999h+rqavj4+CAoKAhTpkzp8YviydPRGwBCCCFKpa+vL100vGfPHhQXF2Po0KHw9/dXyqLhCRMm4NKlSzA2Nsbo0aPx448/dul+d3d3fPvtt/D390dOTo6CqpQdl8tFUFBQj1oMvHPnTjDGEBgYKLeYzc3N8Pb2hqamJg4cOKDUwX9+fj5cXV2xbNkyREZG4ty5cy/V4L+urg4JCQnw8PCAjY0NcnJy8NFHH6GiogIJCQlwd3enwX9vp5rdRwkhhJD/Ky8vZ+vXr2cDBgxgtra2LDo6mjU2Nio87759+5iBgQFbvHgxe/jwYZfuTUxMZHp6eiwzM1NB1cmuuLiYCQQCdu/ePVWXwpqbm5m5uTk7cuSIXGN6enoyNzc3pZ57cP/+fRYaGso0NTVZaGhol//N9GQtLS0sNTWV+fn5MQ0NDTZ27FgWHR3dI/4NEfmjNwCEEEJUztTUFBERESgrK8OGDRtw9OhRmJqaYunSpbh27ZrC8gYEBOC3335DdXU1RowYgczMzBe+d968eYiNjYWvry9++eUXhdUoi8GDB8PZ2RmJiYmqLgVbtmzBwIED4evrK5d4ra2tmDNnDv744w+kpqZCU1NTLnGfRSQSISYmBlZWViguLkZhYSFiYmKgq6ur8NyKJJFIkJOTg7CwMJiamiI8PBzDhg3D9evXpVt3DhgwQNVlEkVQdQdCCCGEdObixYssJCSEaWlpMRcXF5aUlKTQE1Tj4+OZnp4eCwkJYY8ePXrh+/bv38/09PRYRkaGwmqTxd69e9nw4cNVWkN9fT0zNDSU299Na2srmzlzJnNyclLat+/p6enM3t6e2dnZsePHjyslp6IVFBSwNWvWMDMzM/bqq6+y1atXs8uXL6u6LKJEtAiYEEJIj9bQ0IDdu3djy5YtaG5uRlBQEFauXAkzMzO557p9+zaCg4NRWlqK7777DpMnT36h+w4ePIiQkBB8//33cHNzk3tdsmhpaYGJiQlOnDiBcePGqaSGjz/+GBcvXsTPP//c7VhisRiBgYEoKipCRkaGws+U+P3337F69Wrk5ubik08+wcqVK6GmpqbQnIp069YtHDhwAPv378edO3cwa9YszJ8/H1OmTOnVz0VkpOoOhBBCCHkRYrGYpaenMz8/P6apqcn8/PxYenq63PNIJBIWGxvLdHR0WGhoKGtpaXmh+/bu3ct0dXXZjz/+KPeaZLV8+XK2bNkyleSuqqpiOjo67MKFC92O1dbWxgICAtioUaPYgwcP5FDd09XV1bHw8HCmra3NQkJCWE1NjULzKVJFRQWLjo5mLi4uTCAQMG9vbxYfH6+U9TWkZ6M3AIQQQnqdW7duYfv27di5cycGDBiApUuXYvHixdDW1pZbjhs3bmDhwoVobm7Gzp07MX78+Ofek5qaivnz52PLli1YtGiR3GqR1aVLlzBlyhRUVlZCS0tLqbk/+OADlJWV4dChQ92KI5FIsHDhQly9ehWZmZlyP0H4yTx79+7Fhx9+iOHDhyM6OhrDhw9XSC5FqqurQ1paGpKTk3HixAmMHz8efn5+mD9/Pvr376/q8khPoeoOhBBCCJFVS0sLS0pKYs7OztL5+4WFhXKLLxKJWHR0NNPV1X3htQGnTp1i+vr67D//+Y/c6uiOUaNGsYSEBKXmrKysZNra2qygoKBbcSQSCQsJCWG2trasqqpKTtV1lJeXxxwdHZm1tTVLSkpSWB5FaWpq6nQHn7t376q6NNJD0S5AhBBCei0+nw8/Pz+cOXNGuoPPhAkT4OHhgeTkZLS1tXUrPo/HQ1hYGPLz8/H777/DwcEBGRkZz7zH1dUVGRkZ2LRpE9auXdut/PLw+GRgZfr888/h5eUFBwcHmWMwxrB8+XKcOnUKmZmZMDY2lmOFfyovL0dQUBDc3d0xdepU/Prrr/Dz85N7HkUQCoVIS0tDUFAQjIyMsHbt2g47+AwcOFDVZZIeiqYAEUIIeanU19cjPj4eMTExEAqFWLJkCVasWNHt7QwZY9i+fTvWrFmDd955B5s3b4aent5Try8qKoKnpye8vLywZcsWlS20rK2thZmZGa5evYohQ4YoPN/du3dhbW2Ns2fPYsSIETLFYIzh/fffx4kTJ3Dq1Kkun9b8PE1NTfj666/x2WefwdfXF5s2beoVg2WJRILc3FwkJycjMTEROjo6mDt3Lt59913Y2tqqujzSi1ADQAgh5KUkkUiQmZmJuLg4HD16FN7e3ggJCYG7u3u34lZWVuJvf/sbLly4gK1btz5zf/vy8nJMnToVQ4YMwb59++S6RqEr5s2bB0tLS3z22WcKz7Vy5UrU1tZ26wyCDz/8EElJScjKysKgQYPkWB2QlpaG0NBQGBkZISYmBk5OTnKNL29isRjZ2dlITk7G4cOHIRAIMHfuXMybNw+jR49WdXmkl6IGgBBCyEvv5s2b2LFjB3bs2IFBgwZh6dKlCAwM7NbC2Pj4eHzwwQeYPn06vvzySxgZGXV6XWNjI9555x1UVFTg6NGjCtm+9HnS09OxaNEilJaWSt9EFBQUQF9fHxYWFjLHFQqF4PP50j9XVlbCxsYGFy5cgJ2dnUwx//nPfyIxMRFZWVndqu2vLl++jFWrVqG4uBj//ve/sWDBAnA4HLnFlyeJRILTp09LB/1cLhdz5syBv78/nJ2de2zdpPegNQCEEEJeekOGDEFUVBTKy8sRFhaGbdu2wcTEBEuXLsVvv/0mU8yFCxfi+vXrEIvFsLOzw86dO9HZd2q6urpISUmBs7MznJyccOXKle4+Tpe5u7uDz+fjyJEj2LZtG+zs7DBq1Cjk5+d3K669vT3Cw8Nx7949AMDGjRvh6+sr8+D/448/xr59+/DLL7/IbfD/4MEDhIWFwdXVFc7Ozrhx4waCgoJ63CD6yVN5zc3NMXv2bDx8+BBxcXEoLS1FTEwMXFxcelzdpJdS1epjQgghRJWePGnY3d2dJSUlsba2NpliZWZmMhsbGzZx4kR27dq1p14XHR3N9PT0lHpWgEQiYZmZmWzs2LGMx+MxHR0dBoBpa2uzxMREmeO2trYyDofDNDQ0mIaGBgsKCmKamprsxo0bMsX717/+xczMzNjNmzdlrumv9UVHRzN9fX3m7e3NiouL5RJX3q5du8bCw8OZqakpMzQ0ZAsWLGCpqakKPfWaEGoACCGE9GnV1dUsKiqKDRo0iFlaWrKoqCiZDn9qampi69evZ9ra2iw8PJw1Nzd3et2BAweYrq4u27ZtW3dLf6ba2lr26aefMhMTE6apqcl4PB4DIP3R0dFhu3btkjl+SUkJU1NTk8bj8/lMTU2NzZ07t8tNQGRkJDM2NmbXr1+XuZ4npaens2HDhrHRo0ezrKwsucSUp2vXrrH169czS0tL1q9fPxr0E6WjBoAQQghh/z9p2NvbmwkEAubn58dycnK6HKegoIA5OTkxKysrduLEiU6vyc3NZa+++ioLCQlhQqGwu6V3qqWlhY0bN67DwP/xj7a2Ntu6davM8XNycpiWllaHuOrq6ozH4zFfX9923+a3trYyLy8vdu/evXZxvvjiC2ZkZCSX8xuKioqYl5cXMzQ0ZNHR0TK/0VGEx4N+KyurdoP+1tZWVZdG+iBaA0AIIYQA4HK5cHd3R1paGq5evQpLS0vMmDED48aNQ1xcHJqaml4ojoODA86cOYNVq1bBz88PgYGBuHv3brtrXn/9dVy8eBEFBQV48803UVVVJffn4fP5OH78OAYOHAgej9fhc8YYWlpaZI5fVlYGLrfjMEIkEoHL5eL3339vd2rvvn378NNPP+H111+XrhmIjo5GVFQUMjIyMGzYMJlrqa+vx9q1azFmzBhYWVnh1q1bCAsLU9nWq48VFhYiIiIC1tbWmDRpEoqLi/HVV1+hqqoKCQkJ8PHxgbq6ukprJH0TNQCEEELIX1hbWyMqKgoVFRUIDQ3F1q1bYWpqirCwMJSUlDz3fi6Xi5UrV+LGjRvgcrmwsbFBREQEWltbpdeYmJggKysLtra2GDduHM6fPy/35zA0NMTJkyfb7dTzmEQi6XYD0NlBa3w+H5aWlsjOzoaBgYE0V2RkJBhjKCsrw4QJE7Bx40Zs2LABx44dw/Dhw2WqQSKRICEhAUOHDsWlS5dw7tw5xMTEQF9fX+bn6g6JRILs7GyEhobCzMwMb775JiorK/HNN9+gurqaBv2k51D1KwhCCCGkN7h48SJbsGAB4/P5XV40nJ2dzRwcHJiNjQ07fvx4h8+//PJLpqury/bs2fPUGGKxWObas7KymIaGRoepOuvXr5c55sqVKxmHw2kXk8/nM1tbW3b//v121yYnJzOBQCC9TkNDg+no6HT6d/GiTp06xUaOHMmsra1ZWlqazHG6SywWs9OnT7PQ0FBmYmLSbiEvTe8hPRU1AIQQQkgXVFVVsaioKPbaa68xKysrFhUV1WHA2xmRSMRiY2OZoaEh8/b2ZqWlpe0+T09PZwMGDGDLli1jLS0t7T6rq6tjrq6urKGhQea6d+7cyfh8vnQQzuVyWXh4uMzxpk6d+kKDf8YYGzZsWIdmQUNDg7322musoqKiS3nLysrYggULmI6ODlu/fn2HvytlaGtrkw76Bw4cyAYMGEALeUmvQlOACCGEkC4wNjZGeHg4SkpK8NVXX+HkyZMwMzODv78/zp49+9T7eDweQkJCUFRUBEtLS9jb2yMiIgJCoRDAn3v15+fn4+rVq3B2dkZxcbH03kWLFuH06dNYtmyZzHUHBwdj+fLl0NTUBPDndJUXXdfQmdLSUul/8/l8DB48GDk5Oe3m/QPAzz//jOLi4g5nJLS2tqKyshJOTk6orKwE8Oc6AbFY3Gm+pqYmREREYOjQoQD+PNwtIiKi0+lNiiAWi6X79JuYmCAgIAAAkJyc3G5Of2frLQjpcVTdgRBCCCG9XVFREQsPD2evvPIKGzt2LIuNjWVNTU3PvOfChQtswoQJzNramh07dkz6e5FIxMLDw5m+vj77+uuvmaurq3T6DJ/PZwcPHpS5TrFYzLy8vJimpiYDwIKDg2WOZWBgIP0m38bG5qlvQcaPH9/h2//HPzwej3E4HLZmzRqWmJjIOBwO+/bbb9vdL5FIWFJSEhs0aBBzdHRkeXl5MtfcVX/88Qc7cuQIW7BgATMwMGCWlpbsww8/ZOfOnWMSiURpdRAib9QAEEIIIXLy8OFDFhsbyxwcHJiRkRELDw9nJSUlT71eLBaz+Ph41r9/f+bt7c1u374t/SwpKandPvuPf7S1tdmdO3dkrvHRo0fM3t6ecTgcNnfuXJliiEQixuVyGY/He+bg/+zZsx3WHuCJrUI9PDzYpUuXWE5OjvQ6PT09VldXxxhj7NKlS2zixInM1NSUxcfHK2XQXVtbyxISEtisWbOYlpYWs7OzYx999BG7fPmywnMToiwcxjo5t5wQQggh3XLp0iXExMQgKSkJkyZNQmhoKLy9vcHhcDpcW1tbi8jISOzatQsLFy7E5s2bkZ2djWnTpnWYOqOhoQFHR0dkZWV1GutFVFZWYuTIkRg9ejTWrl2LiooK3L17F7W1tWhoaIBEIkFbWxt4PB4EAgF0dXVhbGwMU1NTmJmZQVdXF8OGDYONjQ1yc3M7TPt5bMqUKcjKypJO69HQ0IBEIoGXlxc+//xz2NnZoaSkBKNHj8bDhw/BGINAIEBgYCC0tLSwa9curFixAuvWrYOOjo5Mz/oiampqcOzYMSQnJ+PEiRMYMmQI/Pz84O/v363tSQnpqagBIIQQQhSoqqoK8fHx+OabbyAQCPDee+9hyZIl6NevX4dr8/Ly4ObmJt2//tGjR53G1NTURFRUFEJDQ1+4jsbGRuTk5CA7Oxu5ubm4cuUKGhsbMWTIEJiYmMDMzAzGxsZQV1eHrq4ueDwempqaIBQK8ccff6C8vBxVVVUoKSnBvXv3oK6uDnd3d3h4eGDy5MkYNWpUu3MBCgoKMH78eIhEIvD5fIjFYsyePRsbNmyAjY0NgD8bnzFjxqCyshIikUh6r5qaGt566y3ExsbC3Nz8hZ+xK0pKSpCamork5GScO3cOjo6O8PPzw+zZs2FmZqaQnIT0FNQAEEIIIUogzH+mewAAC8dJREFUFovx008/YcuWLThz5gzmzJmDv//97xg5cqT0mnPnzmHy5Mloa2uDRCJ5Zjw+n4/8/HzY2dk99Zq6ujocOnQI33//PTIyMjB48GBMnjwZkyZNgoODA/T09DB48OAuP0t2djbq6+tx8+ZNZGVlIScnB+rq6pg5cybmzJmDKVOmYPbs2UhJSYGamhreffddrFu3DoMGDZLGEIlEcHV1xeXLl6ULoR/j8XhwcnLC6dOnu1zbsxQWFiI5ORlHjx5FYWEh3N3d4ePjg5kzZ8LY2FiuuQjpyagBIIQQQpSsqKgIW7duxa5duzB06FCEhIQgKCgICxcuxJEjRzo9YOuv1NTUYGNjgytXrkBDQ6PdZ+fPn8fWrVuRlJSEMWPG4O2338bbb78NCwsLhTyPRCLBmTNncOTIERw6dAiMMdTU1GDBggXYsGEDTExM2l3PGENAQABSUlLQ3NzcaUyBQIBDhw7By8sLwJ8NQ1cP0JJIJMjNzcXRo0dx5MgR1NTUwMPDA97e3pg1axZ0dXVle2BCejlqAAghhBAVaWhoQEJCArZu3Yra2lrU1ta+0OD/MU1NTaxcuRKbNm0CAOTm5iIyMhLnzp1DcHAwQkJCYGtrq6jyOyUWixETE4Njx44hLy8Py5Ytw5o1a2BkZCS95pNPPsHmzZufexLxoEGD8N///hdZWVn4xz/+gfz8/HbTjDojFApx+vRppKWlISkpCW1tbZg2bRr8/Pzg6enZoVkipC+iBoAQQghRMcYYFi9ejMTExOcOigFIF//yeDy0tbXh8OHDOHjwII4fP45Vq1Zh1apVMDAwUHTZz3XlyhVs2LABmZmZiIyMxIoVK7B//34sXrwYra2t7a7l8/ngcDgQCoUwMzPD2LFj4eTkhPLycsTFxYExhtOnT8PR0bFDnqamJmRkZCA5ORkpKSno168fZsyYAR8fH7zxxhu0Nz8hf0ENACGEEKJiYrEYpqamqK6ufuZ1AoFA2iDweDyMGjUKNTU1KCsrw9y5c/HVV1/1yLns2dnZWLFiBYRCofSAs8ff5FtZWcHR0RFOTk4YNWoUHBwcoKWlhbKyMsyaNQuFhYVoaWmBhoYGVq9ejc8//xwA8ODBA/z4449ITk5Geno6rKys4OPjA29vb0ycOFFlz0pIb0ANACGEEKJiKSkpmDdv3lPnwz+mpqYGBwcH1NXV4c6dO5BIJNDV1cXu3bvh6OgIU1NTJVXcdeXl5fD09MTNmzexaNEivP/++xg6dGin386npKRg/vz5EAqF7aZEmZubY9WqVfjhhx+Ql5cHJycn+Pr6YtasWTItZiakr6IGgBBCCFExDw8PZGZmQiAQgMvlgsPhQCKRQCKRQCQSdVgXsHbtWmRkZIDL5SIwMBDz58/HK6+8oqLqu+aXX36Bv78/li9fjoiIiHZnGQiFQqxevRo7duzoMEUI+LMBcnFxQWBgIGbOnNluXQEh5MVRA0AIIYSoWFVVFerr69HQ0ICHDx+ivr5e+ueGhgbU1dXh/v37ePDgAWpqanDjxg2MGDECubm5Xd4ZpycoKirCtGnT4OPjg5iYGOnvfH19UVpa+tQ3IZqamti4cSNWr16tzHIJeelQA0AIIYT0Eo2NjXjzzTdhb2+P3bt3y3wScE9QVlYGFxcXLF26FObm5li6dCna2tqeuwvSuHHjcOHCBSVVScjLiRoAQgghpJcICAhAXV0dUlNTe+U3/3918eJFuLi4oLW1FTweDwKBQDr9CYD0QLQnTwnmcrm4d+8eDA0NVVU2Ib0e7YtFCCGE9AJ79uzBqVOnUFBQ8FIM/gHAwsICERER2LRpE+Li4sDlctHW1oaGhgYIhUI0NTXh0aNHEIlEqKmpgVAoxKNHj3Dr1i1qAAjpBnoDQAghhPRwjY2NsLa2xvbt2+Hj46PqcuRu7ty5MDAwQGxsrKpLIaRPoAaAEEII6eEiIyORnZ2NjIwMVZeiECUlJRg+fDjy8/NhY2Oj6nIIeelRA0AIIYT0YGKxGK+99hri4+Ph7u6u6nIU5t1334WhoSG++OILVZdCyEuPq+oCCCGEEPJ0J0+ehLq6Otzc3FRdikIFBwdj79690gXAhBDFoQaAEEII6cEyMzMxffp0cLmK+V/2o0ePwOFwpD9r164FAGzcuFH6u+joaIXkftLEiRPR2tqKwsJChecipK+jXYAIIYSQHuz8+fMIDg5WWHwdHR00NzfDxMQEYWFhWL9+PQBg3bp1aGxsxIgRIxAYGKiw/I9xuVw4OjoiLy8PI0aMUHg+QvoyegNACCGE9GDl5eWwsLBQaA6BQIDAwEBs374dYrEYANDa2oqUlBTMmTNHobmfZGFhgYqKCqXlI6SvogaAEEII6cHq6+thYGCg8DyLFy9GRUUF0tLSAACHDx+Gp6cnBAKBwnM/1q9fP9TW1iotHyF9FTUAhBBCSA8mEAjQ3Nys8DwODg4YP348tm3bBgD49ttvERISovC8T2pqaoKWlpZScxLSF1EDQAghhPRg/fv3R01NjVJyLV68GOnp6UhNTUVbWxvs7e2VkvexmpoaOuGXECWgBoAQQgjpwUaMGIHLly8rJde8efOgpaWFwMBALFmyRCk5n3Tx4kU4ODgoPS8hfQ01AIQQQkgP5uzsjKysLKXk0tXVhb+/P7hcLvz9/ZWS87GqqircunULTk5OSs1LSF9E24ASQgghPZivry9WrVqFsrIymJubKzyfh4cHtLS0lD4Xf8+ePXB3d4e+vr5S8xLSF9EbAEIIIaQHGzhwIN566y3ExcUpPFdLSwuSk5OVvvi3ra0N27dvx3vvvafUvIT0VdQAEEIIIT3cunXrEBMTg7t37yosPofDgZGREYYPH670efjbt2+HmpoaZs2apdS8hPRVHMYYU3URhBBCCHm2gIAAiMViHDx4UNWlyFVlZSVGjhyJ3bt3w8vLS9XlENInUANACCGE9AL379/HyJEj8emnnyI4OFjV5ciFWCyGp6cnLCwssGPHDlWXQ0ifQYuACSGEkF6gf//+2LdvH2bMmAELCwu4ubmpuqRuW7FiBWpqapCSkqLqUgjpU2gNACGEENJLvPHGG4iNjcXs2bNx5swZVZcjM8YYwsPD8fPPP+PYsWPQ1tZWdUmE9Cn0BoAQQgjpRebNm4eWlhZMmzYN8fHxvW7hrEgkwpIlS5CTk4OTJ0/CxMRE1SUR0udQA0AIIYT0MosWLYKxsTECAgKQl5eHjRs3Ql1dXdVlPdft27fxzjvvgDGG3NxcGBkZqbokQvokmgJECCGE9ELTp0/HhQsXkJGRARcXF+Tn56u6pKeSSCSIi4vD6NGj4eTkhOzsbBr8E6JC1AAQQgghvZS1tTVyc3Ph7e2NyZMnY8WKFaiqqlJ1We1kZ2fDyckJmzdvxoEDBxAdHQ0+n6/qsgjp06gBIIQQQnoxDQ0NfPLJJygoKEBNTQ2GDBmCDz74ALdv31ZZTYwxpKenw83NDb6+vvD19cWvv/4KT09PldVECPk/agAIIYSQl4ClpSWSkpKQl5eH6upq2NnZYfr06Th06BCampqUUsOdO3ewefNmDB06FEFBQXBzc8Pt27fx0UcfQSAQKKUGQsjz0UFghBBCyEuouroau3btwt69e3H79m1MnToVnp6emDRpEuzs7OSSo6WlBefPn0dWVhbS0tJQUFAAd3d3BAcHY8aMGb1iYTIhfRE1AIQQQshL7rfffsMPP/yAzMxMnD17FlpaWhg5ciTs7e1ha2sLMzMzmJiYYMCAAdDR0QGfz4eWlhbq6+shkUjw4MEDVFdXo7KyEqWlpbh+/TquXbuGa9euoX///nB1dYWnpyd8fHxgYGCg6sclhDwHNQCEEEJIHyISiVBQUICrV6+isLAQRUVFKC8vR1VVFe7du4fOhgV8Ph8DBw6EqakpzM3NYW9vD3t7e4wZMwYWFhbKfwhCSLdQA0AIIYQQKaFQ2G7NgIGBATgcjgorIoTIGzUAhBBCCCGE9CG0CxAhhBBCCCF9CA9AnaqLIIQQQgghhCjH/wDPJoeBialHRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.view_model()\n",
    "from IPython.display import Image, display\n",
    "display(Image(filename=\"causal_model.png\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                          \n",
      "─────(Expectation(y|Z1,X1,W1,Z0,W3,W2,W0,X0))\n",
      "d[v₀]                                        \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0,U) = P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0)\n",
      "\n",
      "### Estimand : 2\n",
      "Estimand name: iv\n",
      "Estimand expression:\n",
      "Expectation(Derivative(y, [Z0, Z1])*Derivative([v0], [Z0, Z1])**(-1))\n",
      "Estimand assumption 1, As-if-random: If U→→y then ¬(U →→{Z0,Z1})\n",
      "Estimand assumption 2, Exclusion: If we remove {Z0,Z1}→{v0}, then ¬({Z0,Z1}→y)\n",
      "\n",
      "### Estimand : 3\n",
      "Estimand name: frontdoor\n",
      "No such variable found!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "identified_estimand= model.identify_effect(proceed_when_unidentifiable=True)\n",
    "print(identified_estimand)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Model \n",
    "First, let us build some intuition using a linear model for estimating CATE. The effect modifiers (that lead to a heterogeneous treatment effect) can be modeled as interaction terms with the treatment. Thus, their value modulates the effect of treatment. \n",
    "\n",
    "Below the estimated effect of changing treatment from 0 to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                          \n",
      "─────(Expectation(y|Z1,X1,W1,Z0,W3,W2,W0,X0))\n",
      "d[v₀]                                        \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0,U) = P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+W1+Z0+W3+W2+W0+X0\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 10.690211742892949\n",
      "\n"
     ]
    }
   ],
   "source": [
    "linear_estimate = model.estimate_effect(identified_estimand, \n",
    "                                        method_name=\"backdoor.linear_regression\",\n",
    "                                       control_value=0,\n",
    "                                       treatment_value=1)\n",
    "print(linear_estimate) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EconML methods\n",
    "We now move to the more advanced methods from the EconML package for estimating CATE.\n",
    "\n",
    "First, let us look at the double machine learning estimator. Method_name corresponds to the fully qualified name of the class that we want to use. For double ML, it is \"econml.dml.DML\". \n",
    "\n",
    "Target units defines the units over which the causal estimate is to be computed. This can be a lambda function filter on the original dataframe, a new Pandas dataframe, or a string corresponding to the three main kinds of target units (\"ate\", \"att\" and \"atc\"). Below we show an example of a lambda function. \n",
    "\n",
    "Method_params are passed directly to EconML. For details on allowed parameters, refer to the EconML documentation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                          \n",
      "─────(Expectation(y|Z1,X1,W1,Z0,W3,W2,W0,X0))\n",
      "d[v₀]                                        \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0,U) = P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+W1+Z0+W3+W2+W0+X0 | \n",
      "Target units: Data subset defined by a function\n",
      "\n",
      "## Estimate\n",
      "Mean value: 12.12120239567344\n",
      "Effect estimates: [12.1212024 12.1212024 12.1212024 ... 12.1212024 12.1212024 12.1212024]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LassoCV\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "dml_estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     control_value = 0,\n",
    "                                     treatment_value = 1,\n",
    "                                 target_units = lambda df: df[\"X0\"]>1,  # condition used for CATE\n",
    "                                 confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\":LassoCV(fit_intercept=False), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=False)},\n",
    "                                               \"fit_params\":{}})\n",
    "print(dml_estimate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True causal estimate is 10.230462989973262\n"
     ]
    }
   ],
   "source": [
    "print(\"True causal estimate is\", data[\"ate\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                          \n",
      "─────(Expectation(y|Z1,X1,W1,Z0,W3,W2,W0,X0))\n",
      "d[v₀]                                        \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0,U) = P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+W1+Z0+W3+W2+W0+X0 | \n",
      "Target units: \n",
      "\n",
      "## Estimate\n",
      "Mean value: 11.64277316710624\n",
      "Effect estimates: [11.64277317 11.64277317 11.64277317 ... 11.64277317 11.64277317\n",
      " 11.64277317]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "dml_estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     control_value = 0,\n",
    "                                     treatment_value = 1,\n",
    "                                 target_units = 1,  # condition used for CATE\n",
    "                                 confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\":LassoCV(fit_intercept=False), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=True)},\n",
    "                                               \"fit_params\":{}})\n",
    "print(dml_estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CATE Object and Confidence Intervals\n",
    "EconML provides its own methods to compute confidence intervals. Using BootstrapInference in the example below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                          \n",
      "─────(Expectation(y|Z1,X1,W1,Z0,W3,W2,W0,X0))\n",
      "d[v₀]                                        \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0,U) = P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+W1+Z0+W3+W2+W0+X0 | \n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 11.884263413815686\n",
      "Effect estimates: [11.88426341 11.88426341 11.88426341 ... 11.88426341 11.88426341\n",
      " 11.88426341]\n",
      "95.0% confidence interval: (array([11.81925706, 11.81925706, 11.81925706, ..., 11.81925706,\n",
      "       11.81925706, 11.81925706]), array([12.69017104, 12.69017104, 12.69017104, ..., 12.69017104,\n",
      "       12.69017104, 12.69017104]))\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LassoCV\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from econml.inference import BootstrapInference\n",
    "dml_estimate = model.estimate_effect(identified_estimand, \n",
    "                                     method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     target_units = \"ate\",\n",
    "                                     confidence_intervals=True,\n",
    "                                     method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\": LassoCV(fit_intercept=False), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=True)},\n",
    "                                               \"fit_params\":{\n",
    "                                                               'inference': BootstrapInference(n_bootstrap_samples=100, n_jobs=-1),\n",
    "                                                            }\n",
    "                                              })\n",
    "print(dml_estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Can provide a new inputs as target units and estimate CATE on them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[11.87749573 11.87749573 11.87749573 ... 11.87749573 11.87749573\n",
      " 11.87749573]\n"
     ]
    }
   ],
   "source": [
    "test_cols= data['effect_modifier_names'] # only need effect modifiers' values\n",
    "test_arr = [np.random.uniform(0,1, 10) for _ in range(len(test_cols))] # all variables are sampled uniformly, sample of 10\n",
    "test_df = pd.DataFrame(np.array(test_arr).transpose(), columns=test_cols)\n",
    "dml_estimate = model.estimate_effect(identified_estimand, \n",
    "                                     method_name=\"backdoor.econml.dml.DML\",\n",
    "                                     target_units = test_df,\n",
    "                                     confidence_intervals=False,\n",
    "                                     method_params={\"init_params\":{'model_y':GradientBoostingRegressor(),\n",
    "                                                              'model_t': GradientBoostingRegressor(),\n",
    "                                                              \"model_final\":LassoCV(), \n",
    "                                                              'featurizer':PolynomialFeatures(degree=1, include_bias=True)},\n",
    "                                               \"fit_params\":{}\n",
    "                                              })\n",
    "print(dml_estimate.cate_estimates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Can also retrieve the raw EconML estimator object for any further operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<econml.dml.dml.DML object at 0x7f16ffc87100>\n"
     ]
    }
   ],
   "source": [
    "print(dml_estimate._estimator_object)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Works with any EconML method\n",
    "In addition to double machine learning, below we example analyses using orthogonal forests, DRLearner (bug to fix), and neural network-based instrumental variables. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binary treatment, Binary outcome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            X0        X1   Z0        Z1        W0        W1        W2  \\\n",
      "0    -0.187533  0.389821  1.0  0.839453  0.940559  1.415391 -0.380913   \n",
      "1    -0.748730  0.671668  0.0  0.475463  0.892004  1.294085  1.611890   \n",
      "2    -0.808333 -0.732631  0.0  0.868789  0.428347  0.675622  0.962348   \n",
      "3    -1.752145 -0.339531  1.0  0.914105  1.286270 -0.819968  1.832863   \n",
      "4    -0.108159 -1.170617  0.0  0.671865  0.987034  0.268190 -0.551639   \n",
      "...        ...       ...  ...       ...       ...       ...       ...   \n",
      "9995 -1.499062  0.679497  1.0  0.762626  2.221605 -0.350508 -0.256909   \n",
      "9996 -1.196394 -0.888150  1.0  0.562215  1.147892 -0.895030  2.146514   \n",
      "9997  0.246945  0.305430  1.0  0.942246 -0.310852  0.226676  1.843343   \n",
      "9998 -1.626697  1.613843  0.0  0.646475 -0.423104  2.829836  1.762817   \n",
      "9999 -2.709163  2.278675  1.0  0.743850  0.863928  0.361682  0.568979   \n",
      "\n",
      "            W3  v0  y  \n",
      "0     0.813579   1  1  \n",
      "1     0.351258   1  1  \n",
      "2     0.414451   1  1  \n",
      "3     0.184053   1  1  \n",
      "4    -0.152793   1  1  \n",
      "...        ...  .. ..  \n",
      "9995 -0.506154   1  1  \n",
      "9996 -0.534854   1  1  \n",
      "9997 -0.333092   1  1  \n",
      "9998  0.728452   1  1  \n",
      "9999  0.667374   1  1  \n",
      "\n",
      "[10000 rows x 10 columns]\n"
     ]
    }
   ],
   "source": [
    "data_binary = dowhy.datasets.linear_dataset(BETA, num_common_causes=4, num_samples=10000,\n",
    "                                    num_instruments=2, num_effect_modifiers=2,\n",
    "                                    treatment_is_binary=True, outcome_is_binary=True)\n",
    "# convert boolean values to {0,1} numeric\n",
    "data_binary['df'].v0 = data_binary['df'].v0.astype(int)\n",
    "data_binary['df'].y = data_binary['df'].y.astype(int)\n",
    "print(data_binary['df'])\n",
    "\n",
    "model_binary = CausalModel(data=data_binary[\"df\"], \n",
    "                    treatment=data_binary[\"treatment_name\"], outcome=data_binary[\"outcome_name\"], \n",
    "                    graph=data_binary[\"gml_graph\"])\n",
    "identified_estimand_binary = model_binary.identify_effect(proceed_when_unidentifiable=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using DRLearner estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                          \n",
      "─────(Expectation(y|Z1,X1,W1,Z0,W3,W2,W0,X0))\n",
      "d[v₀]                                        \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0,U) = P(y|v0,Z1,X1,W1,Z0,W3,W2,W0,X0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+W1+Z0+W3+W2+W0+X0 | \n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 0.6194410343673133\n",
      "Effect estimates: [0.61944103 0.61944103 0.61944103 ... 0.61944103 0.61944103 0.61944103]\n",
      "\n",
      "True causal estimate is 0.1553\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegressionCV\n",
    "#todo needs binary y\n",
    "drlearner_estimate = model_binary.estimate_effect(identified_estimand_binary, \n",
    "                                method_name=\"backdoor.econml.drlearner.LinearDRLearner\",\n",
    "                                confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{\n",
    "                                                    'model_propensity': LogisticRegressionCV(cv=3, solver='lbfgs', multi_class='auto')\n",
    "                                                    },\n",
    "                                               \"fit_params\":{}\n",
    "                                              })\n",
    "print(drlearner_estimate)\n",
    "print(\"True causal estimate is\", data_binary[\"ate\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instrumental Variable Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/25\n",
      "10000/10000 [==============================] - 1s 101us/step - loss: 11.6538\n",
      "Epoch 2/25\n",
      "10000/10000 [==============================] - 1s 59us/step - loss: 3.7950\n",
      "Epoch 3/25\n",
      "10000/10000 [==============================] - 1s 56us/step - loss: 3.0525\n",
      "Epoch 4/25\n",
      "10000/10000 [==============================] - 0s 48us/step - loss: 2.9154\n",
      "Epoch 5/25\n",
      "10000/10000 [==============================] - 1s 58us/step - loss: 2.8358\n",
      "Epoch 6/25\n",
      "10000/10000 [==============================] - 1s 53us/step - loss: 2.7859\n",
      "Epoch 7/25\n",
      "10000/10000 [==============================] - 1s 56us/step - loss: 2.7401\n",
      "Epoch 8/25\n",
      "10000/10000 [==============================] - 1s 58us/step - loss: 2.7221\n",
      "Epoch 9/25\n",
      "10000/10000 [==============================] - 1s 58us/step - loss: 2.7038\n",
      "Epoch 10/25\n",
      "10000/10000 [==============================] - 1s 52us/step - loss: 2.6872\n",
      "Epoch 11/25\n",
      "10000/10000 [==============================] - 1s 55us/step - loss: 2.6518\n",
      "Epoch 12/25\n",
      "10000/10000 [==============================] - 1s 53us/step - loss: 2.6039\n",
      "Epoch 13/25\n",
      "10000/10000 [==============================] - 0s 48us/step - loss: 2.5914\n",
      "Epoch 14/25\n",
      "10000/10000 [==============================] - 1s 52us/step - loss: 2.5683\n",
      "Epoch 15/25\n",
      "10000/10000 [==============================] - 1s 55us/step - loss: 2.5612\n",
      "Epoch 16/25\n",
      "10000/10000 [==============================] - 0s 47us/step - loss: 2.5512\n",
      "Epoch 17/25\n",
      "10000/10000 [==============================] - 0s 49us/step - loss: 2.5505\n",
      "Epoch 18/25\n",
      "10000/10000 [==============================] - 1s 53us/step - loss: 2.5420\n",
      "Epoch 19/25\n",
      "10000/10000 [==============================] - 0s 49us/step - loss: 2.5394\n",
      "Epoch 20/25\n",
      "10000/10000 [==============================] - 0s 50us/step - loss: 2.5293\n",
      "Epoch 21/25\n",
      "10000/10000 [==============================] - 0s 50us/step - loss: 2.5286\n",
      "Epoch 22/25\n",
      "10000/10000 [==============================] - 1s 51us/step - loss: 2.5297\n",
      "Epoch 23/25\n",
      "10000/10000 [==============================] - 1s 50us/step - loss: 2.5300\n",
      "Epoch 24/25\n",
      "10000/10000 [==============================] - 0s 49us/step - loss: 2.5244\n",
      "Epoch 25/25\n",
      "10000/10000 [==============================] - 1s 50us/step - loss: 2.5150\n",
      "Epoch 1/25\n",
      "10000/10000 [==============================] - 1s 122us/step - loss: 26619.1314\n",
      "Epoch 2/25\n",
      "10000/10000 [==============================] - 1s 66us/step - loss: 12579.7654\n",
      "Epoch 3/25\n",
      "10000/10000 [==============================] - 1s 69us/step - loss: 11734.7099\n",
      "Epoch 4/25\n",
      "10000/10000 [==============================] - 1s 80us/step - loss: 11379.2674\n",
      "Epoch 5/25\n",
      "10000/10000 [==============================] - 1s 70us/step - loss: 11188.0021\n",
      "Epoch 6/25\n",
      "10000/10000 [==============================] - 1s 69us/step - loss: 11216.7738\n",
      "Epoch 7/25\n",
      "10000/10000 [==============================] - 1s 69us/step - loss: 11229.4640\n",
      "Epoch 8/25\n",
      "10000/10000 [==============================] - 1s 68us/step - loss: 11196.7148\n",
      "Epoch 9/25\n",
      "10000/10000 [==============================] - 1s 68us/step - loss: 11116.8998\n",
      "Epoch 10/25\n",
      "10000/10000 [==============================] - 1s 68us/step - loss: 11170.7726\n",
      "Epoch 11/25\n",
      "10000/10000 [==============================] - 1s 67us/step - loss: 10935.9323\n",
      "Epoch 12/25\n",
      "10000/10000 [==============================] - 1s 65us/step - loss: 11147.3363\n",
      "Epoch 13/25\n",
      "10000/10000 [==============================] - 1s 70us/step - loss: 11035.0295\n",
      "Epoch 14/25\n",
      "10000/10000 [==============================] - 1s 65us/step - loss: 11127.9812\n",
      "Epoch 15/25\n",
      "10000/10000 [==============================] - 1s 66us/step - loss: 10992.1307\n",
      "Epoch 16/25\n",
      "10000/10000 [==============================] - 1s 67us/step - loss: 10733.7403\n",
      "Epoch 17/25\n",
      "10000/10000 [==============================] - 1s 75us/step - loss: 11034.8012\n",
      "Epoch 18/25\n",
      "10000/10000 [==============================] - 1s 75us/step - loss: 10997.0383\n",
      "Epoch 19/25\n",
      "10000/10000 [==============================] - 1s 77us/step - loss: 11014.5592\n",
      "Epoch 20/25\n",
      "10000/10000 [==============================] - 1s 80us/step - loss: 10885.7880\n",
      "Epoch 21/25\n",
      "10000/10000 [==============================] - 1s 63us/step - loss: 10905.1800\n",
      "Epoch 22/25\n",
      "10000/10000 [==============================] - 1s 64us/step - loss: 10852.0827\n",
      "Epoch 23/25\n",
      "10000/10000 [==============================] - 1s 63us/step - loss: 10641.3063\n",
      "Epoch 24/25\n",
      "10000/10000 [==============================] - 1s 67us/step - loss: 10856.0904\n",
      "Epoch 25/25\n",
      "10000/10000 [==============================] - 1s 72us/step - loss: 10879.3091\n",
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: iv\n",
      "Estimand expression:\n",
      "Expectation(Derivative(y, [Z0, Z1])*Derivative([v0], [Z0, Z1])**(-1))\n",
      "Estimand assumption 1, As-if-random: If U→→y then ¬(U →→{Z0,Z1})\n",
      "Estimand assumption 2, Exclusion: If we remove {Z0,Z1}→{v0}, then ¬({Z0,Z1}→y)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+W1+Z0+W3+W2+W0+X0 | X1,X0\n",
      "Target units: Data subset defined by a function\n",
      "\n",
      "## Estimate\n",
      "Mean value: 0.716738760471344\n",
      "Effect estimates: [ 0.84010315  2.239563   -0.25856018 ...  0.10548401  0.7332916\n",
      "  1.9136963 ]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import keras\n",
    "from econml.deepiv import DeepIVEstimator\n",
    "dims_zx = len(model.get_instruments())+len(model.get_effect_modifiers())\n",
    "dims_tx = len(model._treatment)+len(model.get_effect_modifiers())\n",
    "treatment_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_shape=(dims_zx,)), # sum of dims of Z and X \n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(64, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(32, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17)])                \n",
    "response_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_shape=(dims_tx,)), # sum of dims of T and X\n",
    "                                    keras.layers.Dropout(0.17), \n",
    "                                    keras.layers.Dense(64, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(32, activation='relu'),\n",
    "                                    keras.layers.Dropout(0.17),\n",
    "                                    keras.layers.Dense(1)])\n",
    "\n",
    "deepiv_estimate = model.estimate_effect(identified_estimand, \n",
    "                                        method_name=\"iv.econml.deepiv.DeepIV\",\n",
    "                                        target_units = lambda df: df[\"X0\"]>-1, \n",
    "                                        confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{'n_components': 10, # Number of gaussians in the mixture density networks\n",
    "                                                              'm': lambda z, x: treatment_model(keras.layers.concatenate([z, x])), # Treatment model,\n",
    "                                                              \"h\": lambda t, x: response_model(keras.layers.concatenate([t, x])), # Response model\n",
    "                                                              'n_samples': 1, # Number of samples used to estimate the response\n",
    "                                                              'first_stage_options': {'epochs':25},\n",
    "                                                              'second_stage_options': {'epochs':25}\n",
    "                                                             },\n",
    "                                               \"fit_params\":{}})\n",
    "print(deepiv_estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Metalearners"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            X0        X1        X2        X3        X4   Z0        Z1  \\\n",
      "0     1.889075 -0.148901 -3.321291  0.527714  0.408656  1.0  0.332018   \n",
      "1    -0.278462  0.502559  0.345945 -0.578400 -0.602006  0.0  0.857228   \n",
      "2    -0.868973 -2.328919 -1.114627  0.268724 -0.375727  1.0  0.172623   \n",
      "3     0.960801  0.437432 -0.439229 -0.478421 -0.482326  0.0  0.346306   \n",
      "4     1.747583  0.397312 -0.002248 -1.231510  0.872643  0.0  0.414153   \n",
      "...        ...       ...       ...       ...       ...  ...       ...   \n",
      "9995  0.647206 -0.211020 -0.175888  0.409196 -0.119107  0.0  0.802795   \n",
      "9996 -1.129477  0.329195 -0.385835 -1.099640 -1.171596  0.0  0.979090   \n",
      "9997  2.553457  1.340551 -0.524767  0.153859 -0.121225  0.0  0.605828   \n",
      "9998 -0.658872  1.386591 -2.662249  0.248188  0.063398  0.0  0.999788   \n",
      "9999 -0.476124  0.275236 -0.498938 -0.335998  1.497155  0.0  0.093024   \n",
      "\n",
      "            W0        W1        W2        W3        W4  v0          y  \n",
      "0     0.075309  2.149438  0.589466 -0.679515  0.481774   1  17.872428  \n",
      "1     0.332104  0.380258  0.176676 -0.066274 -0.130537   1  11.208160  \n",
      "2     1.347730  0.914570 -0.668152 -0.196757 -0.231634   1   3.795458  \n",
      "3     0.332007  0.577365  2.100847 -0.672990  0.513451   1  20.095006  \n",
      "4     0.716787  0.784172  0.792574 -1.049425 -0.882550   0   3.246151  \n",
      "...        ...       ...       ...       ...       ...  ..        ...  \n",
      "9995 -0.366199  0.152741  2.344528  0.286931  1.124170   1  20.967334  \n",
      "9996  1.710311  0.621336 -1.317947 -0.555407 -0.550753   1   3.074311  \n",
      "9997  0.575043  1.996095 -0.348195  0.046330  0.311212   1  31.600315  \n",
      "9998  0.749128  0.608419  1.121918 -0.629710  0.075430   1   7.718591  \n",
      "9999  0.526421  1.659487  0.514313 -0.880840 -1.048985   0   5.398943  \n",
      "\n",
      "[10000 rows x 14 columns]\n"
     ]
    }
   ],
   "source": [
    "data_experiment = dowhy.datasets.linear_dataset(BETA, num_common_causes=5, num_samples=10000,\n",
    "                                    num_instruments=2, num_effect_modifiers=5,\n",
    "                                    treatment_is_binary=True, outcome_is_binary=False)\n",
    "# convert boolean values to {0,1} numeric\n",
    "data_experiment['df'].v0 = data_experiment['df'].v0.astype(int)\n",
    "print(data_experiment['df'])\n",
    "model_experiment = CausalModel(data=data_experiment[\"df\"], \n",
    "                    treatment=data_experiment[\"treatment_name\"], outcome=data_experiment[\"outcome_name\"], \n",
    "                    graph=data_experiment[\"gml_graph\"])\n",
    "identified_estimand_experiment = model_experiment.identify_effect(proceed_when_unidentifiable=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                                      \n",
      "─────(Expectation(y|Z1,X1,X4,W1,Z0,X3,W3,X2,W4,W2,W0,X0))\n",
      "d[v₀]                                                    \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,X4,W1,Z0,X3,W3,X2,W4,W2,W0,X0,U) = P(y|v0,Z1,X1,X4,W1,Z0,X3,W3,X2,W4,W2,W0,X0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+X4+W1+Z0+X3+W3+X2+W4+W2+W0+X0\n",
      "Target units: ate\n",
      "\n",
      "## Estimate\n",
      "Mean value: 11.281848958239234\n",
      "Effect estimates: [ 8.24880761  8.85381156  1.98326246 ... 20.4788987   3.77948304\n",
      " 14.93714404]\n",
      "\n",
      "True causal estimate is 8.619365154692815\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "metalearner_estimate = model_experiment.estimate_effect(identified_estimand_experiment, \n",
    "                                method_name=\"backdoor.econml.metalearners.TLearner\",\n",
    "                                confidence_intervals=False,\n",
    "                                method_params={\"init_params\":{\n",
    "                                                    'models': RandomForestRegressor()\n",
    "                                                    },\n",
    "                                               \"fit_params\":{}\n",
    "                                              })\n",
    "print(metalearner_estimate)\n",
    "print(\"True causal estimate is\", data_experiment[\"ate\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Avoiding retraining the estimator "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once an estimator is fitted, it can be reused to estimate effect on different data points. In this case, you can pass `fit_estimator=False` to `estimate_effect`. This works for any EconML estimator. We show an example for the T-learner below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Causal Estimate ***\n",
      "\n",
      "## Identified estimand\n",
      "Estimand type: nonparametric-ate\n",
      "\n",
      "### Estimand : 1\n",
      "Estimand name: backdoor\n",
      "Estimand expression:\n",
      "  d                                                      \n",
      "─────(Expectation(y|Z1,X1,X4,W1,Z0,X3,W3,X2,W4,W2,W0,X0))\n",
      "d[v₀]                                                    \n",
      "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,Z1,X1,X4,W1,Z0,X3,W3,X2,W4,W2,W0,X0,U) = P(y|v0,Z1,X1,X4,W1,Z0,X3,W3,X2,W4,W2,W0,X0)\n",
      "\n",
      "## Realized estimand\n",
      "b: y~v0+Z1+X1+X4+W1+Z0+X3+W3+X2+W4+W2+W0+X0\n",
      "Target units: Data subset provided as a data frame\n",
      "\n",
      "## Estimate\n",
      "Mean value: 13.155958458684603\n",
      "Effect estimates: [14.58240352 11.70222833 15.95011529 10.52888709 13.01615806]\n",
      "\n",
      "True causal estimate is 8.619365154692815\n"
     ]
    }
   ],
   "source": [
    "# For metalearners, need to provide all the features (except treatmeant and outcome)\n",
    "metalearner_estimate = model_experiment.estimate_effect(identified_estimand_experiment, \n",
    "                                method_name=\"backdoor.econml.metalearners.TLearner\",\n",
    "                                confidence_intervals=False,\n",
    "                                fit_estimator=False,\n",
    "                                target_units=data_experiment[\"df\"].drop([\"v0\",\"y\"], axis=1)[9995:],                        \n",
    "                                method_params={})\n",
    "print(metalearner_estimate)\n",
    "print(\"True causal estimate is\", data_experiment[\"ate\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Refuting the estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding a random common cause variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Add a Random Common Cause\n",
      "Estimated effect:11.877495728144678\n",
      "New effect:12.27643164609447\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_random=model.refute_estimate(identified_estimand, dml_estimate, method_name=\"random_common_cause\")\n",
    "print(res_random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding an unobserved common cause variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Add an Unobserved Common Cause\n",
      "Estimated effect:11.877495728144678\n",
      "New effect:11.766787905625023\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_unobserved=model.refute_estimate(identified_estimand, dml_estimate, method_name=\"add_unobserved_common_cause\",\n",
    "                                     confounders_effect_on_treatment=\"linear\", confounders_effect_on_outcome=\"linear\",\n",
    "                                    effect_strength_on_treatment=0.01, effect_strength_on_outcome=0.02)\n",
    "print(res_unobserved)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Replacing treatment with a random (placebo) variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Use a Placebo Treatment\n",
      "Estimated effect:11.877495728144678\n",
      "New effect:-0.0003870075346435652\n",
      "p value:0.4908247078335508\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_placebo=model.refute_estimate(identified_estimand, dml_estimate,\n",
    "        method_name=\"placebo_treatment_refuter\", placebo_type=\"permute\",\n",
    "        num_simulations=10 # at least 100 is good, setting to 10 for speed \n",
    "        ) \n",
    "print(res_placebo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Removing a random subset of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refute: Use a subset of data\n",
      "Estimated effect:11.877495728144678\n",
      "New effect:11.924618691290394\n",
      "p value:0.4286479951555171\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res_subset=model.refute_estimate(identified_estimand, dml_estimate,\n",
    "        method_name=\"data_subset_refuter\", subset_fraction=0.8,\n",
    "        num_simulations=10)\n",
    "print(res_subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More refutation methods to come, especially specific to the CATE estimators."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}