{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# DoWhy: Different estimation methods for causal inference\n", "This is a quick introduction to the DoWhy causal inference library.\n", "We will load in a sample dataset and use different methods for estimating the causal effect of a (pre-specified)treatment variable on a (pre-specified) outcome variable.\n", "\n", "First, let us add the required path for Python to find the DoWhy code and load all required packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os, sys\n", "sys.path.append(os.path.abspath(\"../../../\"))" ] }, { "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 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let us load a dataset. For simplicity, we simulate a dataset with linear relationships between common causes and treatment, and common causes and outcome. \n", "\n", "Beta is the true causal effect. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Z0Z1W0W1W2W3W4v0y
01.00.8568290.871424-0.792461-0.3363310.386621-0.068865True9.124501
11.00.4910770.197358-0.505399-0.4241400.3677620.168461True8.622930
21.00.6657950.945841-0.2889690.274395-1.3125872.382897True17.977266
31.00.9029051.268346-0.0595300.315513-0.932715-1.360252True8.367090
41.00.104740-1.342788-1.935350-0.649980-0.8524530.843568True-1.326686
..............................
99951.00.5773681.8469290.755214-2.9790111.525415-0.225743True9.653687
99961.00.1310651.880914-1.314365-0.538280-0.3034150.863559True11.305263
99971.00.739417-0.974042-0.707890-0.028049-1.3716080.100693True2.620035
99981.00.489953-0.363797-0.590689-1.905395-0.3743150.622429True1.844830
99991.00.4849421.118425-0.414818-1.1129580.6082691.865714True15.116874
\n", "

10000 rows × 9 columns

\n", "
" ], "text/plain": [ " Z0 Z1 W0 W1 W2 W3 W4 v0 \\\n", "0 1.0 0.856829 0.871424 -0.792461 -0.336331 0.386621 -0.068865 True \n", "1 1.0 0.491077 0.197358 -0.505399 -0.424140 0.367762 0.168461 True \n", "2 1.0 0.665795 0.945841 -0.288969 0.274395 -1.312587 2.382897 True \n", "3 1.0 0.902905 1.268346 -0.059530 0.315513 -0.932715 -1.360252 True \n", "4 1.0 0.104740 -1.342788 -1.935350 -0.649980 -0.852453 0.843568 True \n", "... ... ... ... ... ... ... ... ... \n", "9995 1.0 0.577368 1.846929 0.755214 -2.979011 1.525415 -0.225743 True \n", "9996 1.0 0.131065 1.880914 -1.314365 -0.538280 -0.303415 0.863559 True \n", "9997 1.0 0.739417 -0.974042 -0.707890 -0.028049 -1.371608 0.100693 True \n", "9998 1.0 0.489953 -0.363797 -0.590689 -1.905395 -0.374315 0.622429 True \n", "9999 1.0 0.484942 1.118425 -0.414818 -1.112958 0.608269 1.865714 True \n", "\n", " y \n", "0 9.124501 \n", "1 8.622930 \n", "2 17.977266 \n", "3 8.367090 \n", "4 -1.326686 \n", "... ... \n", "9995 9.653687 \n", "9996 11.305263 \n", "9997 2.620035 \n", "9998 1.844830 \n", "9999 15.116874 \n", "\n", "[10000 rows x 9 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = dowhy.datasets.linear_dataset(beta=10,\n", " num_common_causes=5, \n", " num_instruments = 2,\n", " num_treatments=1,\n", " num_samples=10000,\n", " treatment_is_binary=True,\n", " outcome_is_binary=False)\n", "df = data[\"df\"]\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we are using a pandas dataframe to load the data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Identifying the causal estimand" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now input a causal graph in the DOT graph format." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_model:Model to find the causal effect of treatment ['v0'] on outcome ['y']\n" ] } ], "source": [ "# With graph\n", "model=CausalModel(\n", " data = df,\n", " treatment=data[\"treatment_name\"],\n", " outcome=data[\"outcome_name\"],\n", " graph=data[\"gml_graph\"],\n", " instruments=data[\"instrument_names\"],\n", " logging_level = logging.INFO\n", " )" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "model.view_model()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7gAAAD7CAYAAABAHTgEAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydeVyUVfv/r1lhYNh3EAVcEJRcUxR3SS1Jc8EsGzWzSVvGLaUyoywfx1Ibt55QK1HTpEWjrL6iuWBqibkhKoqiCIrIvgzLzHx+f/hjHkZkEWfmHuC8Xy9eyj33fc7n/szNmbnOuc45PAAgBoPBYDAYDAaDwWAwmjff8blWwGAwGAwGg8FgMBgMhjFgAS6DwWAwGAwGg8FgMFoELMBlMBgMBoPBYDAYDEaLQMi1AIbpKCgoIACkVqupvLycAFBBQYH+9YqKCiorK6vz+rKyMqqoqKjzdT6fTw4ODnW+LhKJSCqV6n+3tbUlsVhMAoGA7O3tiYjIzs6OhEL2GDIsE61WS9nZ2ZSdnU0FBQWk1WqpuLiYNBoN2djYkJWVFUkkEnJ0dCRvb29ycnLiWnKLgvnPLcx/bmH+cwvzn1uY/9zS3P1nkQXHFBYWUn5+vv6npKSE1Go1FRYW6v9fXFxMxcXFVFZWRqWlpVRYWEhlZWWkVqspPz9fX45Op2swaLVUagbD9vb2JBAISCqVkkQiITs7O7K3tyeJREK2trbk4OBANjY2JJFIyMnJiSQSCdnY2JCDgwPZ29uTk5OT/kcgEHB8Z4zmgFqtpn/++YfOnz9PycnJlJycTNeuXaPs7GzS6XSNLkcikZCPjw8FBgZS165dqWvXrtSzZ08KCgoiHo9nwjto3jD/uYX5zy3Mf25h/nML859bWqr/PLaKsnEoKSmh7Oxsunv3Lt27d49ycnLo7t27VFBQYBDAPvj7w+yvHuGsGeDZ2dmRRCIhqVRaK9jj8/kklUpJJBI9NFAUi8Vka2tLRESOjo76B43H45Gjo2Od9yQUCsnOzq7O1x91BLi656eyspJKS0uJ6OGBeXXQXlRURGq1Wh/Uq9VqKisro4KCAlKr1foAv3qE+mHY2dnpg11HR0eD4NfR0ZGcnZ3Jw8OD3N3dyc3Njdzc3MjV1bXOe2K0DHQ6HZ04cYJ+++03Onz4MP3zzz9UWVlJzs7OFBISQsHBwRQYGEheXl7k7e1NHh4e5OzsTHw+X591UP18l5eXU15eHmVlZVFWVhZlZGRQSkoKpaSk0MWLF6myspLc3d1p4MCBNHToUBozZgz5+vpybQGnMP+5hfnPLcx/bmH+cwvzn1taif/fsQC3HiorKykzM5MyMzMpIyODbt++TXfu3NEHsXfv3qXs7GzKyckhtVptcK1UKiU3NzdydnauFVTV97udnR2JxWKO7rj5Up1+XVRUVG+HwoO/5+Xl0b179ww6GoRCoT7QfTD4bdOmDXl7e5OPjw+1bdvWIAWbYfkcPXqUvv32W/r555/p9u3b1L59exo6dCgNHjyYBg0aRG3btjVqfRqNhk6fPk1Hjhyhw4cP0+HDh6m4uJh69+5NEyZMoKlTp5KXl5dR67RkmP/cwvznFuY/tzD/uYX5zy2tzP/WG+BqNBrKyMig9PR0ysjIoFu3bul7H6r/f+fOHf35QqGQPD09ycPDgzw8PMjV1ZXc3NzI09NTH/y4u7uTh4cHubm5kbW1NYd3x3gUtFot5eTk6H/u3LlDOTk5dO/evVodGpmZmQadGfb29tSmTRtq06YN+fj4kK+vL/n4+JCPjw/5+flRQEAASSQSDu+OUVxcTFu3bqUvv/ySkpOTqVu3bjRhwgQaN24cde3a1axaKisr6cCBA7R7927avXs3FRYW0pgxY2jWrFkUHh5uVi3mgvnPLcx/bmH+cwvzn1uY/9zSiv3/jtCCycvLQ3JyMhISEhATE4OoqChERkaiV69ekEgkICIQEcRiMby8vNCrVy9ERkZCoVBAqVQiLi4OiYmJSEtLQ1VVFde3w7AQysrKkJaWhoSEBMTGxkKpVEIulyMiIgK9evWCl5eX/tkiIjg5OemfraioKMTExCAhIQFpaWnQaDRc306Lpbi4GCqVCh4eHrC2tkZkZCQSEhK4lqWnoqICcXFxCA8PB4/HQ7du3RAXFwedTse1NKPA/OcW5j+3MP+5hfnPLcx/bmH+Y2ezD3C1Wi2uXr2K+Ph4rFixAtOnT8eTTz4JJycnfYAhEAjQrl07DB06FK+88gqWLVuGnTt34u+//0Z2djbXt8BogZSUlOD8+fP4+eef8fnnn+Ott97C6NGjERQUBGtra/2zaW1tjeDgYEyYMAHvv/8+duzYgdOnT0OtVnN9C80WjUaD9evXw9nZGQ4ODliyZAny8vK4llUvJ0+exLPPPgsej4fQ0FAkJSVxLanJMP+5hfnPLcx/bmH+cwvzn1uY/3qaV4CblpaGH374AR9//DEmT56M7t27GwQLvr6+GDFiBObOnYsvvvgCf/zxB1JTU1FRUcG1dAZDj06nw61bt3DkyBFs2bIFixcvxoQJExAUFASRSAQiAp/PR/v27fHss89i0aJF+Oabb3Dq1Cn2LDfAyZMn0atXL4hEIixatMjiG/YHOXXqFAYNGgSBQIDXX38dBQUFXEt6JJj/3ML85xbmP7cw/7mF+c8tzH8DLDfAzczMRHx8PKKjoxEREQE3Nzd9IOvl5YXw8HAoFArExMQgMTERhYWFXEtmMB6bqqoqpKWlIT4+Xp/6HBYWBltbWxARhEIhgoODIZPJoFKpkJiYiNLSUq5lc45Wq8Xy5cshEokwePBgXLhwgWtJTUan02Hr1q3w8PCAn58fjh8/zrWkBmH+cwvzn1uY/9zC/OcW5j+3MP8fimUEuHfu3MH333+POXPmICwsDFKpFEQEkUiEHj16YMaMGVi/fj2OHTvGvswzWiUajQYXLlzAtm3bMG/ePAwePBgODg76oDckJATTpk1DTEwMUlJSWsw8ksaQm5uL8PBwiMVirFy5ssXc+927d/H0009DJBJh5cqVXMupE+Y/tzD/uYX5zy3Mf25h/nML879OuAlw09LSsGXLFrzyyisIDAzUz5Pt2bMnZs2ahY0bNyIpKYmlYzIY9aDT6XDlyhXs2rULUVFRGDZsmH6k183NDePGjcPq1atx8uTJFrtIWnp6OoKCgtCuXTucPHmSazlGR6fT4dNPP4VAIMCbb74JrVbLtSQDmP/cwvznFuY/tzD/uYX5zy3M/3oxT4Cbl5eHnTt3QiaTwcfHR7+4zsCBA7F48WL8/vvvLMWYwTACVVVV+Pvvv7Fq1SqMHTsWrq6uICJIpVKMGjUKa9aswdWrV7mWaRQuX74Mb29vdOvWDZmZmVzLMSk//vgjrK2tMWnSJItZeZv5zy3Mf25h/nML859bmP/cwvxvENMFuOfPn4dSqcSgQYMgFAohFAoxZMgQLFu2DImJiSgvLzdV1QwG4/+j0+lw4cIFfPnll5g0aRIcHR1BRAgMDMSCBQtw4MABVFZWci3zkcnMzISfnx9CQ0Ob3UIQTeXQoUOQSCR49dVXOU9DYv4z/80N859bmP/cwvznFuY/tzTBf+MGuOfPn8fbb7+Ndu3agYjg7u6OadOmYdeuXcjPzzdmVQwGowlUVVXh4MGDWLhwIYKDg0FEsLe3xwsvvIDff//dYnon66OkpAQhISEIDg7GvXv3uJZjVn755RcIhUIsXbqUMw3Mf+Y/VzD/uYX5zy3Mf25h/nPLI/r/+AFuTk4O1qxZg169eoGI4O/vj8WLF+PEiRMWl6/OYDAMuXbtGtatW4cBAwaAx+PB29sbCxcuRHJyMtfS6uTVV1+Fi4sLbty4wbUUTtiwYQMEAgEOHz7MSf3Mf+Y/lzD/uYX5zy3Mf25h/nPLI/jf9AD34MGDGDduHMRiMezs7PDyyy/j8OHDnA/dMxiMpnHlyhV88MEH8Pf3BxGhV69e2Lx5s0VNJ9i9ezd4PB52797NtRROGTduHHx9fc2ensT8vw/zn1uY/9zC/OcW5j+3MP+5pZH+P1qAq9Pp8OOPP6Jnz54gIgwcOBDbtm1jW/cwGC0InU6HQ4cOYerUqbCysoKnpyeWL1+OkpISTnWVl5fD398fU6dO5VSHJXDv3j24uLjgnXfeMVudzP//wfznFuY/tzD/uYX5zy3Mf25ppP+ND3APHTqEJ598Enw+HxMnTsQ///zz+CpbAC1t9efWMmHdmLRkz7KysvDuu+/Czs4OHh4eWL9+PWfzdFeuXAkbGxvcunWLk/otjc8//xwSiQQ3b940S33Mf0OY/9zC/OcW5j+3MP+5hfnPLY3wv+EANy8vD9OnTwePx8OoUaNw+vRp46p8DPbt24cXX3wRRAQiwtSpU5GSkqJ//ciRIxg7diyICIMGDcKePXuMUq9Go4FSqcSAAQMgFAqNUiaXlJeXY9myZejXrx8EAkGjromLi0NERAR69OiBESNGYMyYMXjjjTegVCrx9ttvm0zr2rVrERUVhaFDh2LgwIG4fPmyyeqqj6Z41pzJycnB/PnzIRaL0bt3b7O3AxqNBm3atMHChQuNWu7ixYthbW0NIsJTTz2FxMRE3Lp1C6+//rq+XRk/fjwOHjyov+bgwYPo27cv+Hw+3nrrLYM9hr/66itERkZi8eLFmDlzJnbs2GFUvTWpqKhAmzZtzNKL3Fz8z8zMxNdff41JkyahX79+RtX6IMz/2v5v3rwZ3bt3h1QqRbdu3fD1118bVW9NWoL/D5Keng4igoODA/r27YvRo0cjIiICERERGD16NIRCIYgI33zzjf4ac7Y5NWH+38ecbU5NWoL/xm5/Lly4gLFjx8LFxQWurq6YPHkysrKyjKq5GuZ/bf9rsnbtWhCZbifaRvhff4D777//wt/fH97e3vjxxx+Nr9AIlJeXg4jg6Oj40Pm/t2/fBhEZ/SFXq9VwdnY26RtoThp7Pzk5ORg6dCg6dOiAv//+W39cp9Nh+/btcHFxwSuvvGISjWvWrIFUKoVGo0FBQQHGjx/PaSZBS3sGGsOFCxcwYMAAWFtbY+PGjWar97fffgOPx8OVK1eMXvZ7770HIqrVMTNu3DgQEbZt21brmo0bN2L69OkGx5YuXQo/Pz/9ivH5+fnw8/PDmjVrjK65mvfffx+enp51fsgYi+bgfzU3b94EEaFz585G1/ogzP//8c477+Cll17Chg0bMGfOHEgkEhAR1q1bZ3TN1bQE/2ty+PBhDB069KFZQevWrQMR4bnnntMf46LNqUlr978ac7Y5NWkJ/hur/UlJScG4ceOwe/dunD59GjKZDESE4cOHG11zNcz/h3Py5EnY2NiY/LtxA/7XHeAeOXIEUqkUgwcPxp07d0yn0AjU17BotVoQkUlWdO7cuXOLCm4auh+dToewsDA4OzsjNzf3oeccOnQIkydPNpm+wMBAk5TdVFraM9AYdDodlEol+Hw+Fi1aZJY6p06digEDBpik7Ly8PNjY2KBNmzYG7cS///4LIkJ4eHita6ZPn44TJ07of7958yZEIhGWL19ucN6yZctgY2NjsuX809LSQETYv3+/ScqvxtL9fxBzfdlk/t8nIyMDU6ZMMTjn//7v/0BE6NChg0l0Ay3D/5ps3boVf/zxR63j586dg7W1Nby9vfVtCVdtTk1as/8PwkWA2xL8N1b7v2bNGpSVlel/r6qqgqOjI6RSqUl0A8z/h5Gfn4/FixcjMDDQ5N+NG/D/4QFucnIybG1tMXnyZJP3TBiDhhoWU5nc0oKbhu7nhx9+ABHh008/rbccU43229ramv0DpCFa2jPwKHz99dfg8/lYu3atyevy9/fHRx99ZLLyp0yZAiLCb7/9pj+m0+ng7OwMHo+HtLQ0/fHS0lL07NnT4Pr//Oc/IKJaGQXHjx8HEWHFihUm0x4QEIAPP/zQZOUDlu//g5jzyybzHzh69OhDO8Ld3Nxgb29vMt1Ay/C/muLi4lqZaGVlZQgODgafz8eBAwf0x7lsc2rSWv1/EC4CXKBl+G/s9h+4H+BKpVLMmTPHJJqrYf4bsmDBAhQWFprtu3E9/u/k0wMAoJdeeom6d+9O27ZtI6FQ+OApzZr4+Hh67bXXyNfXlwoKCmj69Onk6upKISEhdOrUKf15RUVFFBUVRe+++y4tWLCARo4cSQsWLKCCgoJaZV69epXGjBlDzs7O1KdPHzp06JD+taSkJAoNDSW5XE6LFi0ioVBIpaWlRERUXl5On376Kc2cOZOefPJJeuqppyg5OZl0Oh0dPnyY5s2bR/7+/pSVlUVDhgyhdu3a0ebNm8nFxYV4PB4tWbJEX89///tfEggEtGnTpnrLrkatVtOCBQvotddeoyVLltB7772n11UXP/30ExERDR8+vN7zxo8f32gfG/N+7N27l2bPnk2lpaV0584dmj17tv73hsrftGkT8fl84vF4RERUXFxMq1evNjjW2GeiMZ415T0tKCio9zmxVF5++WX66KOPaNGiRZSammqyerKzs+n69evUr18/k9Uxbdo0IiLavHmz/tjBgwfJ1taWABgc/+GHH2jMmDEG1x89epSIiNq0aWNw3NfXl4iIzp49axLdRET9+/en48ePm6z85uA/lzD/icLCwsjDw6NWuZWVlTRw4EATqb5PS/C/GqlUqv9cqmbevHmUkpJCixYtomHDhumPc9nm1KS1+m8ptAT/TdH+f/DBB6RSqUilUhlfcA2Y//9j3bp1NGnSJLK3tzeZ1gep1/8HQ97ff/8dPB4PFy5cMGXQbVToEUZwb926BalUCiLCsmXLcOPGDWzfvh1EhL59+wK434vXqVMng16Bu3fvolOnTggICNDPz6juoZg7dy4SEhIQExMDW1tbCAQCnDt3DgDQqVMnODs763sFIyMjcffuXQD3N2y+dOmSvo4RI0bAw8MD9+7dw7Fjx/Q57MuXL8f+/fsxc+ZMlJSU6OeC/P777/prb968iRdffFH/e11lFxUVQaPRoG/fvnj11Vf1r6elpekXUKiLJ598EkTU6JWjG+NjY96Pah58nxv7PrVv377WfdU81hgNjfWsqe9pfc+JJaPRaBAYGIjXX3/dZHX8/fffICKTbmyu1WrRpk0biEQi/UjUiy++iEOHDkEqlRrM8xg8eHCtuTDdu3cHEUGtVhscLysrAxGZdPGRpUuXmjR1vzn4/yANfSYYE+b/w/nrr78gkUjw77//mkw30DL8r4sff/wRRITevXujsrLS4DUu25yatFb/H8ScbU5NWoL/xmx/du/ejUGDBoGI4O/vj82bN5tMN8D8r+b48eNYvXq1/ndzjeDW43/tFOWFCxc2avjfkniUABfAQ3PDPTw8YGVlBeD+ymJEhNu3bxucs3XrVhCRft5h9RtYVFSkP2fNmjUgIkybNg3A/RQtIsKXX34J4P5cjsLCQv1D+7CfX3/91UBnXl6egY7Kykq0bdsWY8aM0R9bsmSJfmXbhspev349iAgXL140KLdTp071PpChoaEP9aUuGutjQ+9HNQ++z4/6PtXkwWMNaWiMZ4/zntb1nDQHoqOjTfrBvnfvXhCRyffhfeedd/SpfXl5eejVqxcA4JVXXgER4aeffsLVq1fRv3//WtdWf6CWl5cbHFer1SAifVmmYMOGDXB1dTVZ+c3B/wcx55dN5n9tNBoNBg8ejJ07d5pUM9By/H+QmzdvwsnJCVKpFKmpqbVe57LNqUlr9f9BuApwW4r/xmp/8vPzkZKSgvXr1+sHFLZs2WIy3cx/IDc3FzNmzDBI7zdXgFuP/7VTlHNzc8nd3f3BwxaNSCQinU730Ne0Wi2JRCKDYw+moBAROTk5UUVFBRER/fXXX0REZGdnZ3DOoEGDiIjo2LFjBsdrnvfcc88REVFKSgoR3U8dtrOzo1mzZlFYWBhVVFSQvb09nTx5krp27UoAav2MHj3aQKeTk1Ot+50zZw79+uuvdO3aNaqqqqLLly9T9+7diYgaLHvfvn1EROTn52dQLp9f63EwIDg4mIiILl68WO951TTWx4bej8ctvzE0pKExnj3Oe1rXc9Ic8PDwoJycHJOVr1ariYjI2traZHUQGabpbN++nSZPnkxERDNnziQioo0bN9KWLVtoypQpta7t3LkzEVGtKQz5+flEROTt7W0y3VKp1KTp7M3Bfy5h/tfmo48+ouHDh+vLMCUtxf+aaLVamjJlCuXn59O6deuoY8eOtc7hss2pSWv131JoKf4bq/1xdHSkoKAgeuONNygmJoaIiLZu3Woy3cx/otmzZ9NLL71EqampdPnyZbp8+bL+u/Ply5fp2rVrJtNdn/+1IpqAgAC6cOFCnQGjJeLn50eFhYUPfS0vL49cXV0fqbzqoCU9Pd3gePUcIwcHhzqvrT6nbdu2REQ0YcIEOnPmDI0cOZKOHTtG/fv3p23btlFubi5du3aNysrKapXRGO9nzpxJtra2tH79etqzZw9NnDhR/1pDZWdmZurPexQGDx5MREQnTpxo1PmP46MllF+Txnj2OO9pXc9Jc+Ds2bPUoUMHk5Vf3RnwsPnvxqRz587Up08funLlCn388cf6hjw0NJS6dOlC+/bto6+//pomTZpU69ouXboQEVFWVpbB8du3bxMR0YABA0ymOzc3l5ydnU1WfnPwn0uY/4b8+uuvZGtra7BGhClpKf7XZNmyZZSYmEiRkZE0ffr0Wq9nZ2dz2ubUpLX6bym0FP9N0f6PHTuWiIjEYrFJNBMx/4nur2MzbNgw6ty5s/7n+vXr+nJHjhxpMt31+V8rwH3++ecpMzOT4uLiTCbI2PTq1Yuys7NrBTpERIcOHXrkRS6qRwD37t1rcDwjI4OIiMLDw+u8tvqciIgIIiKKjo6mgIAA+uOPP2jnzp1UVVVF7733HnXu3JnKyspoxYoVBtdfvHiR1q9f36BGe3t7mjlzJn399de0a9cuGjdunP61hsqu7vl98P4a4qWXXqJevXrRmjVr9B+iD1JRUaHvLXscHxtDY8uvHjWtrKwkovsLqdXVIVIXjfHscd7Tup4TS+f27du0Y8cOeuGFF0xWh4uLCxGRSUeJq6nuxXzyySfJy8tLf/yVV14hnU5HvXv3fmiHmUwmI0dHRzp48KDB8T///JPEYjG9+OKLJtOck5Oj98gUNAf/uYT5/z8SEhLo1q1bFBUVZXDclIuwtCT/ie4vHrV06VJq27Ytbdy4sdbrAGju3Lmctjk1aa3+WwotyX9jt//V31OfeeYZ4wqtAfP//ijzg1mL1d+ZAdCVK1dMprle/x+WuCyXy+Hi4oKrV68aKUvatKSmpsLa2hq9e/dGRkYGgPvzVH/99Vd4enrq56ZW4+fnVys33MfHB0SEqqoqlJWVoWvXrmjTpo3B/M45c+YgLCxMP9k6KCio1nzK119/HWPHjtX/bmNjo9+EvaqqCg4ODujTpw/Ky8sREBAAIsKMGTPw7bff4v3338eIESP0c3qrddaVe3/9+nUIBAJ88sknBscbKvvMmTMQCoVwcXHBH3/8gbKyMvz555+wt7cHEeH69et1en3x4kW0a9cOAQEB+Omnn6DRaABAX8bw4cP1+2M11seG3g/g/l5dRISAgAD9OY0tv3rD6iVLluDKlSv4/PPP4ezsDCLCH3/8Aa1W26CGxnj2OO9pXc+JJVNRUYFhw4ahY8eOJp0fUlZWBrFYjB07dpisjmpyc3MhFouxa9cug+M5OTkQi8WIi4ur89oVK1agY8eOKC4uBgAUFRWhY8eOWLp0qUk1jxw5stYepMakufhfTfUiOx07djSVTAOY//fZv38/hg0bhvXr1+t/1q1bh3nz5uH99983meaW5H9+fj7atm0LgUCAI0eOPPSc9evX47nnngPAXZtTk9bsf02N5mxzatKS/H+c9mf16tX46quv9IuLlpeX47nnnsPzzz9fa+snY8L8fzjmmoNbj/8P3we3pKQEvXv3hq+vr8GKsJbM5cuXMXHiRAQEBMDf3x9+fn6YNGkSzp8/b3Dehg0b9Av/fPLJJygsLIRKpdIfe+edd6BWq1FcXIxFixZhxIgRWLBgARYtWoSlS5eioqJCX1ZCQgKeffZZDBkyBHK5HAqFAhs2bDDYLJmI0LNnTyiVSkyZMgURERH6ADI9PR1jxoyBs7MzPD09IZfLkZOTg9LSUixdulSvSS6X1wrSq5k7dy5yc3NrHa+r7GqOHDmCsLAw2NnZISAgAEqlEoMGDcKsWbNw4MABg3t4kOLiYqxYsQKjR4+Gv78/unbtiu7du2Px4sW1tDTkY2Pej5MnT2LWrFkgIvD5fHz00Uc4e/Zso8oH7neA9O3bF7a2thgxYgRSU1MxcOBAyGQyfPfdd/j8888b9Uw0xrOmvqf1PSeWSGlpKZ599lk4ODiYfJVUAOjTpw/efPNNk9cDADNnzqy1Mml9x2vy1VdfQSaTYfHixYiMjMTGjRtNJRPA/dUPHR0dsX79epPW01z8P3jwIORyOYgIIpEIn376Kc6cOWMqqcz//0/NFeIf/HlwH0Vj0tL8r96P0tnZGREREQY/Tz31FPz9/UFEWLBggf4ac7c5NWH+m7/NqUlL8x9oevv/4YcfokOHDnBycsLs2bMxZ84c7N+/35RSmf/1YI4AtwH/Hx7gAvdHzEJDQ2Fvb4/vv//edAoZDEaz4tKlSwgJCYGLiwuOHz9uljrfffddtGvXrt6Ol9bIgQMHQERISUkxaT3M/4fD/OcW5j+3MP+5hfnPLcx/bmnA/9qrKFfj5OREhw4doilTplBkZCRNnDhRv9AOg8FofVRUVNDHH39M3bt3J4lEQqdOnaLQ0FCz1D19+nS6efMmHThwwCz1NRe++uorCg0NpaCgIJPWw/x/OMx/bmH+cwvzn1uY/9zC/OeWBv1vTJR86NAhdO7cGTY2NlAoFMjOzjZuGM5gMCwWrVaLuLg4tG/fHhKJBNHR0QYp4OZiwIABePrpp81er6Vy48YNWFtbm3wj+2qY/4Yw/7mF+c8tzH9uYf5zC/OfWxrhf90pyg+iVquxevVquLm5wc7ODvPnz8eNGzeMo5TBYFgcpaWl2LBhA9q3bw+RSAS5XK5fxI0LDh48CCJCQkICZ6wSxpcAACAASURBVBosCZlMBn9/f5SXl5ulPua/Icx/bmH+cwvzn1uY/9zC/OeWRvjf+AC3mqKiIqxcuRK+vr4QCoWIiIjA999/b7Y3mcFgmJbjx49j9uzZcHJygkQiwezZsy1mRfVnnnkGISEhrb69+euvv8Dn8/Hdd9+ZtV7m/32Y/9zC/OcW5j+3MP+5hfnPLY30/9ED3GoqKyvx3XffYdSoURAIBHBycsKsWbNw7NixphbJYDA44saNG/jkk08QGBgIIkKXLl3w6aef4u7du1xLMyAtLQ329vaYO3cu11I4o6CgAP7+/nj66adNuv3Bw2D+M/+5hvnPLcx/bmH+cwvzn1sewf+mB7g1yczMxGeffYauXbuCiNC+fXvMnTsX+/bt42SuHoPBaJgzZ87gP//5D8LCwsDn8+Hm5gaFQoFTp05xLa1etm/fDh6P1ypXd9doNBg7diy8vLw463xg/jP/uYL5zy3Mf25h/nML859bHtF/4wS4NTl16hSioqIQEhICIoJUKsW4ceOwefNmZGVlGbs6BoPRSEpLSxEfH4/XXnsNvr6+ICJ4eHhgxowZiI+PR2VlJdcSG81bb70FKysr/Pnnn1xLMStyuRwSiQRHjx7lVAfzn/nPBcx/bmH+cwvzn1uY/9zyiP4bP8CtSXp6Or744guMHj0aEokEPB4PISEheOONN7Bz507cunXLlNUzGK2a4uJi7Nu3D0uWLMGQIUNgbW0NPp+P3r17Izo6GidPnmy2+6pptVpMnjwZ9vb2OHToENdyTI5Wq8XcuXMhFArx888/cy2H+W8Bepj/3Oph/nOrh/nPrR7mP7d6mP8NYtoAtyZlZWXYu3cv5s+fjz59+kAoFIKIEBAQgGnTpmHz5s24dOmSueQwGC2O7Oxs/PTTT5g7dy569+6t/xvr0KEDpk+fjm+++QZ37tzhWqbRqKiowKRJk2BlZYVdu3ZxLcdklJeXY/LkyRCLxdi5cyfXcvQw/7mF+c8tzH9uYf5zS0VFBcaOHQuxWMz854Caz7+5F7syJ4/hv/kC3AcpLS1FYmIilEolIiIi4ODgACKCvb09wsLCoFAoEBsbi+Tk5GY7ysRgmIq8vDwkJiZCpVJBJpMhODgYPB4PfD4fwcHBkMvliI2NRXp6OtdSTYpWq8WcOXPA5/Px7rvvoqqqimtJRuX69esIDQ2Fvb09bGxsIJPJLKojsLX47+DggAMHDnAtpxbMf25h/nML858bsrKyMH/+fNja2qJjx47Mfw7QarX4+eef4enpCR6Px/yvDXcB7oNUVlbi+PHjWL9+PWbMmIHu3btDJBKBiGBnZ4eBAwdizpw52LJlC06ePIni4mKuJTMYJqeqqgqXLl3Cjz/+iA8++AARERHw9vYGEYGI4O/vjwkTJmDZsmX4/fffUVBQwLVkTti0aRNsbGzQr18/pKWlcS3HKOzatQuOjo4ICQnBuXPnEBsbi86dO4PP5yMiIgJJSUlcS9TTUv13cHCAnZ0dvvnmG67l1EtL9b/6+U9JSeFaTr0w/7mF+W8ebty4AYVCAYlEAnd3dyiVSpSWljL/zYharcaXX36JDh06QCAQIDAwEAKBAFZWVsx/QywnwH0Y5eXlOHnyJGJiYiCXy9G7d29YWVmBiMDj8dCuXTuMHDkS8+fPx8aNG3H06FHk5uZyLZvBeGTUajVOnz6NnTt3YsmSJZg4cSK6dOkCsVisf947duyI559/HitWrEBCQgLy8vK4lm1RXLhwASEhIZBIJPj444+b7V5xV69exTPPPAMej4dZs2ahrKxM/5pWq0V8fDx69eoFIkJ4eLjFbM3WEv2fNm0annnmGRARXnzxRdy+fZtreXXSEv1/8Pm3ZJj/3ML8Nx3Xr1+HQqGAtbU12rZtC5VKVUvXuXPn4O3tzfw3EUVFRVCpVPDx8YFYLEZYWBhcXFzg7OwMlUqFc+fOseffEMsOcB+GRqNBamoq9uzZg+XLl2Pq1Kno3bs3pFKpflTL3d0dAwcOxLRp07B06VJs374dx44da1HzDxnNj6KiIpw5cwY//fQTVq5ciddffx2jRo1C+/btwePxQEQQCARo3749xo8fj/feew/bt2/HqVOnUFpayrX8ZkFlZSWUSqU+derbb7+FRqPhWlajyM7OxsKFC2FtbY2uXbs2uHhEQkIC+vXrByJCWFgY4uPjzaS0blqq//Hx8fDz84ODgwNUKpXF3lNL9b+5wPznjqqqKkRFRaFNmzbMfyNx9epVyOVyCIVC+Pv7Q6VSPTRwunr1Kvr37w8rKyu88sorzH8jkp2djejoaDg5OcHOzg7jx49HcHAwRCIR5HK5wXY5rP0xoPkFuPWRnp6OP/74A6tXr8bs2bMxcuRIdOzYUT8KRkSwtbVFSEgIxo4di3nz5mHt2rXYvXs3/v77b2RmZrL5vowmk5ubi/Pnz2Pv3r3YtGkTFi9ejBdeeAF9+/aFm5ub/hnk8Xjw8fHRd8JER0djyJAhCAgIAJ/PBxGhbdu2iIyMxMqVK3HkyBEW4D4iN2/exNSpUyEUCtG5c2fExsZCrVZzLeuh3LhxA2+//TZsbW3h4eEBlUr1SFs27d+/H8OGDQMRITQ0FD/88APnH2ot0f+ysjJER0fDysoK3bt3t5iR84fREv1vTjD/zcuNGzcwYMAAWFtbQ6VSMf8fk5MnT2LixIn6NT22bdtW52dKbGwspFIpunbtitOnTwNgz78xSElJgVwuh5WVFTw9PbFw4UI899xz+syt8+fP13kt8x9ASwtw60Kj0SA9PR1//vknNm/ebBB4eHh46AMPIoJQKESbNm3Qv39/TJw4EXPmzMHKlSuxY8cOJCYm4tKlS8jPz+f6lhhmpLS0FOnp6Thx4gT27NmDtWvX4p133oFMJsOQIUPQqVMnSCQSg+fIzs5O35Eyd+5crF27Fr/++itSUlLqbWiKior0i0dFRkbqn0+BQIDg4GDIZDKoVCokJSWxzphGkJqaiqlTp0IkEsHFxQXz58/HhQsXuJaFyspKxMfHIyIiAgKBAK6urli1atVjdWQcP34cY8eOBZ/PR0BAAFQqFYqKioyo+tFpLv57eXk12v/U1FSMGDECPB4PMpkMOTk5ZlDcNFqi/80J5r/p2b17N5ydnREUFISzZ88avMb8fzT+7//+D8OHDwcRoVevXvj+++/r/J5x584dPPvssxAIBIiKinroyC7z/9GoXjgqPDxcPy1tzZo1WLx4MaytrREYGIhffvml0eW1cv9bR4DbEOXl5bh27RoSExOxY8cOrFy5EnPmzMHEiRPRv39/+Pj46Ldcqf4Ri8Xw8fFBjx498PTTT2Pq1KmYP38+PvvsM2zZsgV79+7FP//8gytXruDevXssGLEgCgsLkZ6ejjNnzuDAgQP49ttvoVKpsHjxYrzyyisYM2YM+vXrh4CAANja2hq870QENzc3dOvWDREREZg1axaWLl2KLVu2YN++fbhw4YLRg4rMzEzEx8cjOjoa4eHhsLGxARFBKpUarDh+/fp1o9bbkrh9+zaWLVsGPz8/EBGCgoLw/vvvIykpyWyjnUVFRdizZw+mTp0KJycn8Hg8hIeHY9CgQXB1dTXaB8/Vq1ehUChga2sLOzs7KBQKzlfTtmT/v//++yb1GMfHx6NNmzb6OVCW3Ma3RP+bE8x/46NWq6FQKEBEkMlkKCkpqfNc5n/dVK/r0KdPH4PpLjqdrs5r4uLi4OLiAn9/fxw+fLjBOpj/9ZObm4uVK1ciICAAPB4Po0aNwq+//ootW7bAw8MDTk5OUCqVqKioaFL5rdT/nTwAIEaDaLVays7OppycHLpz5w7l5ORQTk4OZWdnU3Z2Nt27d4/u3r2rf02tVtcqw97enpycnMjR0ZGcnJwM/l/9r62tLUmlUnJwcCCJREI2Njbk6Oho8P/WSmlpKZWVlVFxcTEVFxeTWq2mkpISKioqorKyMiotLaWCggLKz8/X/1y7do0AUFFRkf6YTqczKFcoFJKbmxu5ubmRh4cHubu7k5ubG7m7u5Onpye5urqSm5sbeXp6kqenJ1lbW3PkwH00Gg1dvnyZTp06RadOnaK//vqLTp8+TTqdjry8vKhXr17Uq1cvGjBgAPXv359sbGw41WtJ6HQ6Onr0KP3000+0e/duunnzJjk4ONCAAQNo4MCB1LNnT+rSpQt5e3s/Vj0ajYZSU1MpOTmZTpw4QYmJiXTmzBnS6XTUv39/Gj9+PI0fP57atWtHZWVlNGrUKEpLS6OjR4+Sv7+/Ue61qKiIvvnmG1q1ahVlZmbSM888Q++99x7169fPKOU3BUv0/3EoLCykDz74gDZs2ED9+vWjL774gkJCQh6rTFPCtf9arZa6dOlCM2bMMIr/zQ2u/Tf2888Vly5dohdeeIGuXbtGX375Jb3wwguNuo75/z8qKipo165dtGzZMrp69So988wz9P7771Pfvn3rvKawsJAWLVpEmzZtoldffZVWrVpFUqm00XUy//8HADp8+DBt2rSJfvrpJxKLxTRt2jR68803KSsri+bNm0fJyck0Y8YM+uSTT8jNze2x62xl/n/HAlwTUVJSQjk5OZSfn68PumoGXw8eq/63tLSUSkpK6i3bxsaGJBIJOTg4kK2tLYnFYhKLxWRra0tERA4ODsTn88na2pokEgkRETk5ORERkUQiMQjQeDxevUGzVColkUj00NcKCgqorsentLSUKisr9b9rtVoqKirSe1NVVUUajYaKi4uJ6P6Xca1WS5WVlVRaWkoAqKCggMrKykitVlNhYWG9nlhZWZGtra1BZ4GdnR3t3buXrKysaNiwYRQeHk6urq61OheM0XBwTXFxMZ09e1Yf8B45coSys7NJIBBQYGCgQdDbo0cP4vP5XEvmHACUnJxMhw8fpsTERDp69ChlZWUREZGzszN16tSJvLy8qE2bNuTh4UH29vZkZWVFNjY2ZGVlRcXFxfpnuLCwkG7dukXZ2dl048YNunLlClVWVpJQKKTg4GAaPHgwDRo0iAYNGkTu7u61tBQWFtKwYcOosLCQjhw58tgfMDWprKykuLg4Wr16NZ0+fZoGDhxIc+fOpTFjxpBQKDRaPY+KOfzn8/nUoUMHGjlyZL3+Py5nzpyh2bNnU1JSEr3++uv0ySefkJ2dndHrMSZcPP9bt26lc+fOUXJy8iN9MW6JWFL705zYunUrvf766xQUFETfffcdtW/fvknltFb/c3Nz6b///S+tXbuWiouLafr06fT222836GNCQgLNmDGDNBoNbdq0iSIiIh5LR2v1/86dO7R9+3batGkTpaamUp8+fejVV1+lyZMnU15eHr3//vu0bds2Cg8Pp88//5y6du1qEh2twH8W4Foq1aOSZWVlBoFeQUEBlZaWklqtpqKiIv1DVl5erh81zs/PJyKisrIyqqioIJ1Opw8QHww8qwPKuqgu62FUB9cPo2bATWQYSFf/gfD5fHJwcDAoSygU6r8YPjhybWNjQzY2NmRvb09SqZRsbGz0o911BWw5OTm0atUqWrNmDbm7u9PixYtpxowZnH6xNxdZWVkGo7zHjh2jsrIykkql1K1bN33QO2jQIPLz8+NarkWQm5tL586do5SUFLp69SrduXOHbt26RXfv3qWioiKqqKjQ/w1Vd/7Y2dmRvb09+fj4kKenJ/n6+lLnzp2pS5cuFBQURFZWVo2qOycnhwYPHkwCgYAOHTpELi4uRr+/Q4cO0erVq2nv3r3k7e1NcrmcZs6cSV5eXkavqykY2/+pU6dS79696euvvza5dgC0bds2WrBgAYnFYlq+fDlNnTrV5PUak4b8LykpoYqKCqqqqmrS83/nzh0KDg6ml19+mVatWsXhnVomZ86codDQUFIoFFRVVWXW9sfSKS4uptmzZ9OOHTvorbfeos8++6zO7x9Nhcv239RcvnyZVCoVbd26laysrGj27NmkUCjIw8Oj3uvUajV99NFH9Nlnn9GECRPov//9r0k+m4harv9qtZr27NlD27Zto3379pFUKqWXXnqJXn31VerWrRuVlpbSZ599RitWrKC2bdvSqlWrHrsDoSm0MP+/Y3NwGa2CmzdvQqFQwMrKCn5+foiJieF8pVlzU1VVheTkZMTGxkKhUKBXr176VZu9vLwQERGB6OhoJCQkcL7YQmslIyMDfn5+6NOnj0kXiLp16xaio6Ph5uYGgUCAiIgIJCQk1DvvqjmyceNGWFlZmXX/2tzcXMjlcvB4PAwbNgwXL140W92mxtXVFV988cVjlbFp0ybw+XyLXoWaK5YtWwZnZ+cmz7VrqSQlJaFDhw5wc3PD3r17uZbTrEhMTERkZCQEAgECAgKgVCobvVDqiRMn0KlTJzg6OmLr1q0mVtqy0Gq1SExMhFwuh729PQQCAcLDwxEbG6ufL67VahEbGwtPT8/HnmfLqAVbZIrRukhPT9fv6xYUFITY2NhWF+jWhK3abHlcuXIFnp6eGDp0qMmX9q+oqEBcXBzCw8NBROjYsSOUSiXu3btn0nrNRXl5Odzd3fHhhx+ave4jR46ga9euEIvFiIqKsthtGh4FkUiEHTt2PFYZOp0Ow4cPxxNPPMH54i6WRpcuXfDGG29wLcNi0Ol0UKlUEIvFGDp0KDIzM7mW1CyobtdrLhwVFxfX6O86VVVViI6OhkAgwIgRI3Dr1i0TK245XLx4EdHR0QgICAARITg4GEqlslYn68GDB9G9e3cIhcJa+9kyjAILcBmtk2vXrkEul0MgEKBLly6Ii4trcaNXTYWt2sw9586dg7OzM5599lmzBQEXL15EVFQUHB0dYW1tjcjISBw9etQsdZuSxYsXw93dnZMAs6qqCiqVCnZ2dmjfvj1+++03s2swFsXFxSAio9xDamoqJBIJli9fbgRlLYOkpCQQEY4fP861FIsgJycHEREREAqFiI6ObtUd0Y3l7t27UCqVaNOmDUQiESIjIx/5eUpOTkbPnj0hkUigUqnY96JGkJeXh5iYGISFhYHH48Hb2xsKhQL//vtvrXNv3rwJmUwGHo+H4cOH49y5cxwobhWwAJfRurlw4QJkMhn4fD6eeOIJxMXFcS3J4mCpzdxw/PhxSKVSTJkyxayj50VFRYiJiUG3bt30+yHGxMTUuw2HJZOZmQmxWIwtW7ZwqkEmk4GIEBERgRs3bnCmpalkZGSAiIyWWrxs2TJYWVm1qBTux2Hu3Llo3749Cyhwf3TLx8cHvr6+SExM5FqOxZOamgqFQgEbGxs4ODhAoVDg5s2bj1RG9Wi5lZUVQkNDcfnyZROpbRmUl5cjPj4ekZGREIvFkEgkiIyMRHx8PKqqqmqdX1JSgujoaFhbW6NTp07su6bpYQEugwEA58+fR2RkJHg8HkJDQxEfH8+1JIuGpTabh/3798PKygqzZ8/mpP6kpCTI5XJIJBI4ODhALpcjOTmZEy2Pw5QpU9CjRw+uZeDAgQMIDAyEra0toqOjm9V8q/Pnz4OIkJKSYpTyqqqq0KNHDwwePLjVB3UajQaenp6cpNJbEhqNRp8aO3bsWOTm5nItyaJJTExEREQEeDweOnToAJVK1aSOyOvXr2PIkCEQiURstLweNBoN9u3bhxkzZsDR0VGfwr1t27Y6fX/YPNvy8nIzK2+VsACXwajJiRMnEBERoZ+38ueff3ItqdnAUptNw+7duyEUCvHee+9xpiE7OxtKpRL+/v4Gc7qayxzK6vTPQ4cOcS0FZWVl+p78J554otmkgScmJoKIkJWVZbQy//nnHwgEAmzevNloZTZHfv/9dxBRqx41y8jIwKBBg2Btbc1SY+uhoqICsbGxCAkJadL82geJjY2FnZ0dunTp8tCU2taOTqfDX3/9hTfffFPfkd+7d2+sWrWqwTnhJ06cQGhoqH6ebXZ2tplUM8ACXAbj4fz1118YPny4/gPk8OHDXEtqdrDUZuOxdetW8Pl8KJVKTnVotVokJCQgMjISQqEQXl5eiIqKahYpt2FhYRg3bhzXMvRcuXIFo0aNAo/Hg0wms/hFRn755RcQEcrKyoxa7rx58+Dg4NCqFxCaMmUK+vXrx7UMzvj555/h4uKCwMBAnD59mms5FklWVhY++OADuLu7w8rKCtOnT8eZM2eaXF52djbGjh0LHo8HhULBRhUfIDk5GdHR0ejQoQOICEFBQYiOjsalS5cavJbNs7UIWIDLYNRHYmIiBg8eDCJCeHg4kpKSuJbUrGGpzU1n3bp14PF4+PLLL7mWAuB+Wts777wDd3d3CIVCjBs3Dnv37rXY9La4uDjw+XxcvXqVaykGxMfHw9fXF05OTlCpVBb73G/btg1isdjo5ZaWliIgIACTJk0yetnNgZKSEkilUmzYsIFrKWanvLwcCoVC38lTXFzMtSSL48SJE3jxxRchFovh7u6OJUuWPPa2Zz/88ANcXV3h5+dnEVktlkJ6ejqUSiWCgoJARPD19YVCoWj0PPCa82w7duzI5tlyCwtwGYzGkJCQgCeffBI8Hg8RERGsl9mIsNTmxhMdHQ0+n//YW7UYk/LycuzYsQNDhw4Fj8eDr68vPvjgA6Snp3MtzQCNRgN/f3/MmzePaym1KCkpQVRUFIRCIXr16oV//vmHa0m1WLduHTw8PExSdnWK7p49e0xSviUTGxsLkUiEnJwcrqWYlevXryM0NBR2dnbYtm0b13IsiuptfsLCwkBE6NGjB2JiYh47e6KgoAByuRxEBJlMZtK91psLt27dgkql0q+A7OLiArlcjsTExEanyet0OsTGxsLLywuOjo5snq1lwAJcBuNRSEhIQM+ePcHn8xEZGdmodBXGo8FSm+vn7bffhkgkwi+//MK1lFrcvHkTSqUSbdu2BZ/P129sb+y01qby6aefws7ODgUFBVxLeShnzpxB//79wefzIZfLUVhYyLUkPR9//DECAwNNVv6UKVPg7e2N/Px8k9VhiTz11FMYM2YM1zLMSlxcHBwcHNCzZ0+kpqZyLcdiqF7rwNfXF3w+HxEREUhISDBK2QkJCfD19YWHhwd+/vlno5TZXMnMzMTatWsxcOBA8Pl8ODo6YsaMGUhISHjkDKTqebZ8Ph8ymYzNs7UcWIDLYDwqOp0O8fHx6Natmz7QvXLlCteyWjQstfl/6HQ6zJw5ExKJBAcPHuRazkPRaDQGc3WdnJwgl8tx9uxZTnXl5+dDKpVCpVJxqqM+qkcD3Nzc4OXlhdjYWK4lAbjfsdK3b1+TlX/v3j24u7vjjTfeMFkdlkZWVhYEAkGrSWUsKyvTpySzeZ//49SpUwar1SsUCqOta1BWVoaoqCjw+XxMnDix1WUKVJORkaEfqeXz+bCzs8MLL7yAPXv2NOk5zMjI0M+zHTZsGOefbYxasACXwWgqWq0WcXFx6NSpE0QiEWQyGdLS0riW1WpozanNGo0GkyZNgr29PU6ePMm1nHrJysqCUqlE+/btDfbV5Wq+3RtvvAF/f3+LnStcTV5eHhQKBfh8PoYMGWK07XmaysyZMzFy5EiT1lG9mFpr2ft05cqVsLe3t5gMB1OSnJyMrl27wtXV1SKzT8yNRqNBfHw8wsPDQUQIDAyESqUyalbS33//jcDAQDg4OCAmJsZo5TYXbt68aZB+bGNjg4iICMTGxjZ5X3c2z7bZwAJcBuNxqQ5027dvD7FYDLlc3qpXBOWK1pbaXFFRgaeffhqurq64cOEC13IaRKvVIjExUT9SYW9vr5/rZE5SU1PB5/ObzXzPkydPonfv3hCJRFAoFE3+Yva4TJw4Ec8//7zJ6xk5ciQCAwNbxehe9+7dMXPmTK5lmJzY2FjY2Nhg8ODBuHXrFtdyOCU/Px8qlUo/jaM6DdmY2yJVVVVBqVRCJBIhPDwcGRkZRivb0rlx44ZBUOvo6AiZTIb4+PjHalO0Wi2++uor/TzbVatWNat9zFshLMBlMIxFZWUlYmJi4OPjAysrK8jl8sde7ZDxeLT01ObS0lIMHDgQ3t7euHbtGtdyGk1+fj5iYmLQrVs3EBGCg4OhVCrNlj43evRoDB061Cx1GYOqqiqoVCrY29sjICAAv/76q9k1hIeH47XXXjN5Penp6ZBKpfjwww9NXheXXLhwwWL2ZjYVBQUFmDRpEgQCAaKjoy0+a8KUnD59Gi+//DKsra3h5OSEhQsXmmQhvuTkZPTs2RM2NjZYv359q9hPOD093SCodXJy0ge1xghC9+/fj27dukEoFGL27NkWv6UbAwALcBkM41NRUYGYmBh4eXnB1tYWCoWCLTxgQbS01OaCggL07NkT7du3b5aZA0lJSZDL5ZBKpbCyskJkZKTRRzQeZN++fSCiZrcaelZWFmQyGYgIERERZl2punfv3oiKijJLXStXroRYLEZycrJZ6uOCqKgotG3btll2rDWGv//+G/7+/vDw8MC+ffu4lsMJ1ashV6chd+rUCSqVyiRZGBqNBitWrICVlRVCQ0Nx+fJlo9dhSSQnJ+OTTz5Bz549QURwdXXFzJkz8ccff6CystIodVy+fBmRkZH6bSLZfrbNChbgMhimorS0FCqVCh4eHpBKpYiKimp1K4Q2B1pCavPdu3cRFBSErl274t69e1zLaRJFRUWIiYlBnz599F8Gly9fbrKUxieeeAIvv/yySco2NX/++SeCgoJgY2OD6Ohos6TKdejQAf/5z39MXg9wPx2wX79+6Nu3b4sMALVaLXx9ffHuu+9yLcXo6HQ6qFQqiEQiPPXUU7hz5w7XksxOZmYmPvjgA3h6ekIgEGD8+PE4cOCAyTrt0tLSMGjQIIhEohY7Uq7VanHs2DEsXLgQHTt21H8+v/baa0hISEBVVZXR6srNzUVUVBTEYjGCgoI4yZhhPDYswGUwTE1xcTGUSiWcnJzg7OyM6Ohoi9r+g1Gb5pjanJGRAT8/P/Tp06fZ72949uxZvPXWW3B2doZAIMCoUaOwc+dOqNVqWQ64ggAAIABJREFUo9WxceNGWFlZNdtpBJWVlVCpVLC1tUWnTp2Mtp1IXbi6uuKLL74waR01OXv2LEQiETZs2GC2Os3FgQMHQEQ4f/4811KMSnZ2NkaOHKkPtCylbTQXSUlJkMlkEIlEcHd3R1RUlEmzLHQ6HWJiYiCVShESEoJ///3XZHVxgUajQWJiIhQKBXx8fEBE8Pf3h0KhQGJiotGfr4qKCqhUKjg6OsLV1RUqlapFdha0EliAy2CYi6KiIiiVSjg4OMDV1RVKpdJiRwQZtWkOqc1XrlyBp6cnhg4datRgkCsqKioQHx+PyMhIiEQiODg4QCaTGSWFuby8HO7u7s1+rmdaWhqeeeYZEBEiIyNNNmImEomwY8cOk5RdF++88w7s7e1b3CI5L7/8Mnr27Mm1DKOSkJAAT09P+Pn54fjx41zLMRtqtRqxsbF44oknDFaJN/XK2Onp6Rg2bBiEQiGioqJazIJHpaWliI+Ph0wmg4ODg36NhujoaCQlJZms3vj4eAQEBMDGxgZRUVFsEKL5wwJcBsPc3Lt3D9HR0bC3t4ebmxuUSmWLCEZaG5aa2nzu3Dk4Ozvj2WefNdpcJEsgNzcXMTEx6NGjh35bjejo6MfqUFi8eDHc3d1bxN9ffHw82rVrB0dHR6OPPBQXF4OI8NtvvxmtzMZQXl6Ozp07Y/To0Wat15So1Wo4ODhg9erVXEsxClVVVYiOjtbvs9papuFcuXIFUVFRcHZ21q8d8Ndff5ml7ri4ODg5OSE4ONjit4lrDPfu3UNsbCwiIyNha2sLgUCAsLAwKJVKpKammrTuf/75BwMHDgSPx0NkZGSzWXuD0SAswGUwuCInJwdRUVGQSCTw9fWFSqVqFVtjtGQsJbX5+PHjkEqlmDJlSotME0xOTkZUVBTc3NzA5/MRHh7epL0NMzMzIRaLsWXLFhMpNS+lpaWIjo6GWCxGz549ceLECaOUm5GRASLCsWPHjFLeo3Do0CHweDx8//33Zq/bFHz33XcQCATNNjW+Junp6ejfvz8kEglUKhXXckyOVqtFQkICIiIiwOPx4OPjg+joaLOt/n779m2MGTMGfD4fcrm8WWeAJScnQ6lUIiwsDAKBABKJBGPHjsU333xjlnUkbt68CZlMBh6Ph9DQUE7aNoZJYQEug8E12dnZiIqKgrW1Ndq1a4eYmBijLpjA4BauUpv3798PKysrzJ4926jlWhLGSGGeMmUKevToYWKl5uXSpUsIDw8Hn8+HTCZ77C+M58+fBxEhJSXFSAofjRkzZsDT0xN5eXmc1G9MIiIiMGrUKK5lPDY//vijfhSxpc0lfpDqvWv9/f3B4/EQHh6OuLg4s35Ox8XFwcXFBf7+/s1ya6nKykrs378fc+bMQUBAAIgIbm5umD59On788UezBevFxcWIjo6GtbU1OnbsiLi4OLPUyzA7LMBlMCyFGzduQC6XQygUwt/fHzExMWyBgxaIOVObd+/eDaFQiPfee89I6i2XpqYwJyUltcj9SHU6HWJjY+Hu7g5PT0/ExsY2ed5yYmIiiAhZWVlGVtk4CgoK4O3tjZkzZ3JSv7G4e/cuRCIRtm/fzrWUJqNWq6FQKEBEkMlkzXoUsSFOnDiBadOmwdraGo6Ojpg3b57JU2Yf5O7du5gwYQJ4PB7kcjmKi4vNWv/jkJubi7i4OMhkMjg6OoKIEBAQAIVCYfSVjxuiqqoKMTEx8PDwgLOzM5RKJcuYa9mwAJfBsDSuX7+uD3SDg4MRGxvbItNMGf/DlKnNW7duBZ/Ph1KpNJF6y+NRU5jDwsIwbtw4M6s0D/n5+VAoFBAIBBg0aFCT9pb95ZdfQEQmXzinPnbt2gUej4f9+/dzpuFxWbduHWxtbZtVkFKTlJQUPPHEE3BwcMB3333HtRyTUFJSgk2bNun3V+3WrRtiYmJMsndtQ+zduxdeXl5o165ds3nu09LSoFKpEB4eDpFIZDCf9tKlS5xoSkhIQEhICEQiEeRyOe7evcuJDoZZYQEug2GpXLx4ETKZDAKBAF27dkVcXJzJ9tFjWB7GTG1et24deDwevvzyS9MLtyDKy8vx/fffIyIiAkKhEI6Ojnjttddw9OhRg7+luLg48Pl8XL16lUO1puXUqVPo06cPRCIRFArFIwVZ27Ztg1gsNqG6xjF27Fh07NiR00D7cejbty9kMhnXMppEbGwsbGxs0KdPH6SlpXEtx+hcvny51qJRpt56qy4KCgogl8v1K6Nbcmp+ZWUl/vzzT8yfP1+/P62LiwtkMhni4uI4XY045f+xd+dxUZbt28CPWdkF2VQEARERFELBFZcU0EzcRe0xMBfGMIU0bcgN9CkdrAzKzMElwBYbLQs0/Ym5hIopiuAuggsSCi4gCLLN+f7RyzySOzBzz8D1/Xz8Q5b7PoZF55zrus7z/HlVh/mAgADKzs7mLAujcazAZRhtd/bsWQoMDCQej0e9evWipKQkriMxHGjs1ua6TqeaHvWiLQoKCuizzz4jd3d31Va5ZcuW0aVLl6impoYcHR1p3rx5XMdUq9raWkpISCALCwtq3749JSQkvNTnffXVV9SmTRs1p3ux/Px8MjMzo48++ojrKK8sOzubeDwe/d///R/XUV5JSUkJvfXWW8Tj8SgsLKxZdWavqamhpKQk8vPzIx6PR05OTiSTyThd4duzZw/Z2tpSmzZt6Ndff+Usx/MUFhY+d+sx1z8jRUVFql0r3t7edOjQIU7zMJxgBS7D6IqsrCwKDAwkANS3b1+d2bLEqM+rbm1esGABiUQiSk5O5jo6p86ePUuRkZHk4OCgmrP45ptvkrGxMRUXF3MdT+0KCgpUHUQDAgIoNzf3uR//3//+l1xcXDSU7vnWrl1LQqGQTp06xXWUV7Js2TJq166dTvVVOH78ODk5OZG1tTXt3r2b6zhNJj8/n2QyGdnZ2amOMCgUCk6/Nw8ePCCJRKIaV6OJTsIvq7a2ltLT01Vdj3k8Hunr65Ofnx/JZDLOms/9W1lZGUVFRZGxsTF16NCBvv/+e7brreViBS7D6Jq0tDQKCAggAOTj40MHDhzgOhKjRV60tblr164kFovpxx9/5Doq52prayk1NZUkEgkZGxsTAOrcuXODRg7pooMHD1LXrl3JwMCAIiMjn9l0ZcGCBdS7d28Np3u62tpa6t+/P/Xs2VOnikVnZ2eaP38+1zFeilKppJiYGBKLxeTr68tZc7GmVDfiJzAwkIRCIbVp04akUqlWzD09fPgwderUiaysrGj79u1cxyGif2bTKhQKkkgk1LZtWwJADg4OJJFISKFQ0IMHD7iOqFJTU0NxcXHUrl07MjU1pVWrVunsMQamybACl2F01eHDh2nIkCEEgPz8/Oj48eNcR2K00L+3Nvfo0YN4PB4BIEtLyybt2qzLKioqaNiwYWRoaEgikYhatWpFQUFBlJSU1KzHdlVVVVFMTAwZGxuTs7PzU7fQzpw5k4YNG8ZBuqe7ePEi6enp0Zo1a7iO8lIOHz5MAHRi1bmwsJCGDx9OQqGQIiMjdb7B4f3790kul5ObmxsBIC8vL5LL5VRRUcF1NCovLyepVEp8Pp/efPNNys/P5zRP3WxaPz8/EgqF9RpEpaenc5rtWVJSUsjDw0PVQOrWrVtcR2K0AytwGUbXpaam0sCBA1WF7smTJ7mOxGi5O3fuUJ8+fcjIyIiGDRvWpF2bddnly5eJz+fTd999R3K5nHx8fAgAtW/fnsLCwpr171ZeXh5NmDBB1ZAlLy9P9b4JEybQpEmTOEz3pMjISDI0NNSJhkfvvvsuubm5cR3jhf744w+ysbEhe3t7OnLkCNdxGuXo0aMUHBxM+vr6ZGpqSnPmzKFz585xHUslLS2NXFxcyNTUlORyOScZioqK6IcffqCpU6eqVmnrxnH9/PPPWrVK+2/p6ek0ePBg1b9Xmh7fxGg9VuAyTHORkpJC3t7eqnN1p0+f5joSo8UePnxIAwYMIBsbG8rNzW3Srs26bMSIETR48GDV3y9cuECRkZHk5OSkOq/7MvN1dVVSUhI5ODiQqakpxcTEUE1NDfn5+dGsWbO4jlZPZWUlubm50RtvvMF1lOeqrKwkCwsLWrVqFddRnqm6ulrVhG7cuHFa3bX3eYqLi2nt2rWqRnLdu3enuLg4rTpuUFVVRZGRkSQQCGjo0KH1XkhSt+rqavrzzz9p8eLF1LNnT+Lz+SQSiWjQoEG0cuVKysjI0Pozqzdu3CCJREJ8Pp969erFGkgxz8IKXIZpTpRKJSUlJVH37t2Jz+dTYGAgXbp0ietYjJYqLi6mHj16kJOT0xPb4xrbtVlX7d27lwBQRkbGE+9LT0+nsLAwsrS0JD6fTz4+PiSXy7V6paMhysvLKTIykvT09MjT05O6dOlCUqmU61hPSEtLU624a6sdO3YQj8fT2hdEbty4Qf379yd9fX2KiYnhOk6DpKenk0QiISMjI9LX1+d0xM/zZGVlkaenJxkaGlJMTIxGismCggJKSEigwMBAVcdjR0dH1Vna+/fvqz1DU7h37x5JpVLS19cne3t7SkhI0PpinOEUK3AZpjmqra0lhUJBLi4uqkK3Oc/4ZBqusLCQXF1dqVu3bi/s3PmqXZt1lYeHB02bNu2Z76+oqKBt27bRqFGjSCwWk5GREb311luUlJRElZWVGkyqXpcvX6ahQ4eqVsOKioq4jvSE0NBQsrCwoNu3b3Md5anGjx9fb0eANtmxYweZm5uTq6srZWZmch3nlZSUlJBcLidPT08CQK6uriSTyeju3btcR3tCdXU1yWQyEovF1K9fP7XOYy0vL6eUlBSSSqXk5eVFAMjQ0FDV8Vhbz9I+S1VVFcnlcrK2tiZzc3OSyWTPbIbHMI9hBS7DNGd1ha6zs7OqCYMmt0QxuiEvL48cHByoV69er7wa2Ry3NsfFxZGenh4VFBS88GPv3LlD69atowEDBhCPxyNzc3OSSCR04MABnS7yH9eqVSsyMzMjc3NziomJ0arHVVJSQra2thQUFMRpjgcPHjyxolRcXEwGBga0adMmjlI9XUVFBYWFhREACgoK0qotvC9St1prbGys1au1dc6dO0fe3t5kYGBAMpmsyX93lEolZWRk0OrVq8nPz4/09fUJAHl4eNDChQtp3759OlsQJiUlkZOTE4nFYgoLC9OZ1WZGK7ACl2FagqqqKkpISKCOHTuSWCwmiUTSLEY/ME0nOzub2rZtS4MHD25Uh9HmsLX50aNHZG1tTVFRUa/0eXl5eRQTE/NEc6rU1FSd3k4nEolo48aNFBYWRgKBgPr3709ZWVlcx1LZuXMnAaA9e/ZwlkEul5OdnR0tWbKELl68qHqbvr6+Vj0xv3DhAnl6elKrVq3ohx9+4DrOS6lbre3evbvWr9bWqa2tpZiYGNLT06PevXvThQsXmuzaN27coM2bN9Nbb71F1tbWBICsrKxo8uTJ9O233+r8/+1paWnUv39/1UzgF83pZpinYAUuw7QklZWVJJfLqX379mRoaEhhYWGsrT6jkpWVRebm5jRy5Eiqqqpqsuvq4tbmxYsXk7W1dYOL/fPnz1NkZCR17tyZAJC9vT1JpdImfaKrCaWlpQSAfv/9dyIiysjIoD59+pBQKKSwsDCtOX88YcIEsre3p9LS0npvf/TokUZmi65bt454PB4JhUICQJ6enuTk5ESjR49W+71fVkJCAhkZGZG3t7dOHFnRtdXaOrm5uTRw4EASiUQUGRnZ6HnNZWVl9bYd83g8MjAwqLftWFv+3WyMS5cuUWBgIPF4PPL19W3WXesZtWMFLsO0RHWFbtu2bcnIyIikUqnOds5kmlZaWhoZGxvTlClT1PqkSdu3Nufn55NYLKb4+PhGX+vs2bMklUrJxsamXidmXRhxk5eXRwDo6NGjqrcplUpKSEggS0tLsrGxoYSEBA4T/qOgoIBat25N8+fPV73t8OHD5OzsTEZGRo0uMl7kq6++IpFIRAAIAPF4POLz+cTj8Wjw4MGUkJDwRPGtKQ8ePKApU6YQj8ejsLAwrT4n/uDBg3qrtV26dCGZTPbC/gDaQKlUklwuJ2NjY3J3d2/w3OOamhpKT09XzaQVi8UEgDp27KhqDqUtLyw1hTt37pBUKiWxWEyurq6kUCi4jsToPlbgMkxLVlZWRjExMWRtbU0mJiYklUqpuLiY61gMx/bt20d6enoUGhqqsXtq49bmKVOmUPfu3ZvserW1tZSamkphYWFkZWWl6sQcExOjtU2Szpw5QwDo/PnzT7zv7t27JJFIiMfj0ZAhQzhfnd6wYQPx+XxKSUmh0NBQ4vF4JBAICIDaV4NiY2NVhci//wgEAuLxeGRoaEgzZ85s0t0RL5Kenk6dOnUiKysr2rVrl8bu+6oeX63V09PTmdXaOteuXaMhQ4aQUCgkqVT6Si8iKJVKysrKotjYWBo1ahSZmJiodn3MmDGDtm7dSoWFhWpMz42HDx+STCYjU1NTsrGxIblcrvYXopgWgxW4DMP8sw1RJpORmZkZWVhYUGRkJJWUlHAdi+HQjh07SCgU0qJFizjLwPXW5vT0dAJABw8ebPJrP3r0iJKSkigoKIiMjY1JIBCQn58fJSQkaNXqTGpqKgF47rm+P//8k7p160b6+voUGRnZqDPcjaFUKmngwIFkbm6u2ioMgEQiEcXGxqr13l988cUzC9zH/6g7Rx2lUkkxMTEkFotp8ODBT4wB0wZ1q7U9evTQudXaxykUCmrdujW5ubnRiRMnXupzLly4QOvWraPAwECysrIiANS6dWsaM2YMff311816vF9d80t7e3vVDjJt+jePaRZYgcswzP/cvXuXIiMjydTUlCwtLUkmk1F5eTnXsRiOJCYmEp/PJ5lMxnUUFU1vbfbx8aGxY8c2ybWepbS0lL777jsaMWIEiUQiMjQ0pMmTJ9OOHTs4KxbrJCcnE4AX/jtQXV1NMTExZGJiQk5OTrR7924NJfzHrVu3aMqUKQRAtfL/+Arq+PHj1Xr/Tz/99LkFrlAopLffflutGeoUFRVRQEAACYXCJjkD2tSetVqra43Ybt26RaNGjSIej0cSieS5O0vy8/NJoVCQRCKhDh06EAAyMjJSnaNNTU3V6Mo+V3bv3k3u7u4kFArp3XffZT1AGHVhBS7DME+qOxNjaGhI1tbWJJPJOH+izXDjq6++Ih6PR+vXr+c6ylOpe2uzQqEgPp+vsaY8d+7coW+++YYGDRpEfD6fWrVqRUFBQZScnMzJ2cktW7aQWCx+6Y/Pz8+noKAgAkABAQF0/fp1Nab7Z6Vyw4YNZGxsXO8M7L//WFhYqDVHdHQ06enpPfXeIpGIvL29NTKu5eDBg9S+fXuys7Oj1NRUtd/vZf17tdbFxUUnV2vrbNmyhVq3bk2dOnWiw4cPP/H+v//+W1XQOjo6qubR+vj4kFQqpZSUFK0+C93UTpw4Qb6+vgSA/Pz8tKoLO9MssQKXYZhnKywsJKlUSvr6+mRnZ0dyuZyqq6u5jsVoWGRkJPH5fJ0ZK9KUW5tramrI0dGR5s2bp6H0/1NUVERyuZx8fHyIx+ORmZkZBQUFUVJSksZWe7766itq06bNK3/evn37yMXFhYyMjCgyMlJtT+blcvkLtwXX/cnOzlZLBiKilStXPrXAFQgEZGVlRTdv3mzU9cvLy597jZqaGoqMjCSBQECjR4/WmhE6zWW1tk5BQQGNHj2aeDwezZkzRzVD+NatW6qC1s3NTbVq7+XlpSpodXUebWNcvnyZJk6cSDwej/r27Ut//vkn15GYloEVuAzDvFheXh6FhYWRnp4eOTg4sEK3BVqwYAGJRCJKTk7mOkqDNGZr8+rVq8nExITTBmyPz9jl8Xhkbm6uKnbV+bv43//+l1xcXBr0ueXl5RQZGUn6+vrk4eFBR44caeJ0/0hMTCSxWFzv3O3TCs1vv/1WLfcnIvr444+fWuCKRCL666+/Gn39sLAw6tOnz1O/13l5eTRw4EDS19enmJgYzovHutVaLy+vequ1RUVFnOZqLIVCQZaWluTg4EA///wzJSUl1RvdIxAI6hW0Lfl4T1FREUmlUtLT0yMXFxdSKBSc/1wyLQorcBmGeXnXrl0jiURCQqGQunTpQgkJCVp3votRD6VSSTNnziQDAwM6cOAA13Ea7fGtzXWrLs/a2vz333+TsbExxcTEcB2biIiuX7+uKnbrtt/WFbtN/fu4YMEC6t27d6OukZ2dTcOGDSMej0dBQUEv1RH2VZuGnTp1imxtbZ+5TVkkEtGMGTMa+hBeaPny5U8tcDdv3tzoa+/fv594PB7xeDxavHhxvff99ttvZGFhQS4uLpSRkdHoezXGsWPHKCQkhExMTEhfX5+CgoK0apt0Qz1+1tbT05O6d+9OfD5fVdCGhYWRQqFgEwjon8kMMpmMWrVqRe3bt2cvhjNcYQUuwzCvLjc3lyQSCQkEAuratSt7dbaFqKmpoYkTJ1KrVq1euluoLnne1mZzc3MyNjamNWvWqK1rc0NcvXqVYmJiVKtl7du3p7CwMEpNTW2S38mZM2fSsGHDmiApUVJSEtnZ2VHr1q0pJibmuV/DuqLhVdy5c4eGDBmiGg307z+Ojo6NfQjPFBkZWa/AFQgEtGDBgkZft6SkhGxsbFSPicfjUUpKClVVVZFUKlW9aMDVjN07d+5QTEwMdevWjQBQt27dKCYmRmu2SDdUWVkZpaSk0OjRo+u9aPL4LFo2O/5/qqqqSC6XU9u2bcnY2Jh1Rma4xgpchmEa7vz58xQUFER8Pp/c3d1ZodsCVFZW0vDhw8nS0pLOnTvHdRy1q9vaPHfuXOLxeKoiRl1dmxvj7NmzFBkZSV26dCEAZGdn1+hid8KECTRp0qQmy1hWVkZSqVR1PvH48eNPfMypU6eIz+eTgYEBnTlz5pWuX1NTQ8uWLSMej/dEN2Uej6e2ecNLly5V/WwIhUIaPHhwk6ymBwcH19t6zefzydzcnNzd3cnExIS+++67Jkj/aurmOUskEjI0NCQDAwOdm1v7b7dv31ZtOfbx8alX1Lq6ulJCQoLONsRSJ6VSSQqFgjp16kRisZgkEonWzvRmWhRW4DIM03hnzpyhwMBA4vF41KdPH0pKSuI6EqNGDx8+pAEDBpCNjQ3l5uZyHUdjRowYQa+//vpLb21uSNfmplJX7Do7OxMAsre3b1Cx6+fnR7NmzWryfKdPn6Z+/fqRUCiksLAw1dzt2tpa8vLyIpFIREKhkDp06NCgwiIpKemJzso8Ho9+/vnnpn4oRES0aNEi1TlgR0dHun//fqOvmZSU9Mzt1s7OzhqflZqfn08ymYycnJwIAHl5eZFcLuds9bgxrl27RomJiRQSEkKurq6qFw88PT1p+PDhZGJiQra2trRv3z6uo2qtlJQU6tGjB/H5fAoMDGxR/xcwWo8VuAzDNJ3MzEwKDAwkAOTj40N//PEH15EYNSkuLqYePXqQk5MT5efncx1HI/bu3UsAnjjr2JRdm9Xhr7/+ooULF6rGlTg4ONAHH3xAR48efWGx6+3tTVKpVC25lEolJSQkkJWVFbVr144SEhJo/fr19VZeRSIR9e3bt0FdmC9fvkwuLi6qFVCxWPxS3bDv3r1LZ8+epdTUVEpJSaHk5GRSKBT0yy+/UEpKCu3fv58yMjIoPz9ftUobERGhGgVz9uzZV876b0VFRWRhYfHEKvTjK7mffvppo+/zIjU1NZSSkkKBgYEkFArJzMyMJBIJnT59Wu33bko5OTmqF6Xqfg/qdhHUbYe/dOkSjR8/XjXXlm2xfbp/j/zh+uw3wzwFK3AZhml6R48eVf0H6OPjQ4cOHeI6EqMGhYWF5OrqSt26dWsx2/c8PDxo2rRpL/y4xnRtVqd/b2O2tbUliUTyzG7MnTp1opUrV6o1U1FREU2fPp14PJ7q6/TvDsjvvfdeg65dWlpK48ePVxWKHh4eRPRPcX3mzBlKSEigDz/8kN58801ycnIifX39lx47VJfNxsaGbG1ticfjUWhoKP3555+Nnhs+fvz45871rbv30aNHG3WfZ7l06RJFRkaSnZ0d8fl88vPzo4SEBJ3oDFxTU0Nnz54luVxOgYGBZGVlRQDIyMio3hzaxx+LQqEgKysrsre31+mt1upU12SSz+dTr169mkWzQabZYgUuwzDqk5qaSq+//rrqld7m2JiopcvLyyMHBwfq1atXi1jxiIuLIz09PSooKHilz3uVrs2a2tpcV+zWze20tLR8Ys6upaUlrVu3TiN5RowY8dyibv369Q26rlKppE8//ZT4fD7x+XwaOXIkWVpaEgDS09MjT09PmjJlCn3yySeUkJBAKSkpdO7cOSooKKB79+6pvh9VVVV07949unPnDl2+fJlSU1Np69at9MUXX5CHhwc5ODhQ69atCQDp6+vTwIEDaenSpXTs2LFX2hb+3XffvXRxbWdn12TdeysqKkihUJCfnx/xeDyysbEhqVRKOTk5TXJ9damqqqL09HSSyWQUEBBAZmZmBIBatWpFfn5+JJPJKDU19am7AAoLC2nChAls1fY52MgfRgexApdhGPVLSUmhnj17qgrdU6dOcR2JaULZ2dnUtm1bGjx4cKNXrrTdo0ePyNramqKiohp9LW3a2pyTk/PU0UMCgYC2bNmi1nsTER0+fJh4PN4LC7pXXTXKy8ujjz/+mHr06KE6h9uvXz9as2YNnTx5sslGmDx+tvnq1auUkJBA06dPp44dO6pWyufMmUPHjh17YV4TE5MXfi3EYrFqR8Dj832Li4tfahv2486ePUtSqZQsLCxIIBCQn58fKRQKrR3vUlpaSikpKardEXWr7u3ataPAwMCX/p2pW7Vt166dzs73Vic28ofRYazAZRhGc/7dlOLixYtcR2KaSFZWFpmbm9PIkSNVq3/N1eLFi8na2lotxbw2bG3Ozc2lmJgY6tM2iwm9AAAgAElEQVSnj+redSu7DTkL+yLV1dXk5uZWr1vwswpcU1PTF64oKpVK2rNnD40ePZoEAgFZWVmRRCKhPXv2UG5uLu3cubPJH8PzZGZmUlRUFLm7uxMA8vT0pG+++YbKysqeyD106FBV8fq05lIAqG3btqpt5Y9/P3Jzc6lz584EgP7666/nZiouLia5XK4q/F1cXEgmk2llB9ybN2/STz/9RGFhYeTp6akameTq6koSiYQSExPp2rVrL329e/fukUQiIQAUGBjIxv38Cxv5wzQDrMBlGEazlEolJSUl0WuvvaYqdC9fvsx1LKYJpKWlkbGxMU2ZMkVr5sSqQ35+PonFYoqPj1f7vbjc2pyXl0cA6P333ycfHx/i8XhkZmZGQUFBpFAonijQXsXjjck+//zzlz7zKhKJqGvXrs+89+O7Req6/GrTudH09HSSSCRkZGRElpaWFBkZqeog/fXXX9drKiUQCFRFv7OzM0ml0md2wT5y5AiZm5uTSCQikUj0zM7Xj99fX19fNd5HW7acVldX04kTJyg2NpYmT55MdnZ29RpCvf/++/TLL79QYWFhg66fnJxMNjY21K5dO9bt/1/YyB+mGWEFLsMw3KitrSWFQkGdO3cmkUhEQUFBWn/Wi3mxffv2kZ6eHoWGhnIdRa2mTJlC3bt35+TemtrafObMGQJA58+fJyKiGzduqLYx8/l8MjQ0pICAAEpISHilUTF123B/+OEHIvqnA/HjK3N1RdqzilyhUEhjx46tV5SlpaWRl5cX8Xg8GjdunNZ3+S0qKqKIiAgyNjYmKysr+vjjj1VbbXk8HgmFQvLz86NvvvmGbt68+dxrKRQKEovFqq9fXUOluhc8CgoKKCYmhrp166Yq/GNiYuju3buaeKjPVVJSotpuHBAQQKampgSATExMyM/PjyIjIykpKanRY5fYqu3z/Xt3Ffu/mNFxrMBlGIZbdYXu468av+gJHaPdduzYQUKhkBYtWsR1FLVJT08nAHTw4EGuoxCRerY2p6amEgD6+++/n3hfXl4eyeVyCggIIKFQSAYGBqpi90XbGdeuXas6Yzpz5kzVVu+Kigo6evQoxcTE1Fu9A0AGBgb1Vjd5PB4tX76c7t27RyEhIcTn88nX15cyMzNf6evGtaKiInr//feJx+ORQCCgESNG0LZt215qS6hSqSSZTEY8Hu+JM7t8Pp8++OADGj16NAmFQjI3N6ewsDDOvz6Pj+txc3NT5e7YsaPazpxv376d2rRpQ3Z2drR79+4mu25zcOzYMdXEgxEjRlBWVhbXkRimKbACl2EY7VBVVUUJCQnk6OioKnSf9sSa0Q2JiYnE5/NJJpNxHUVtfHx8aOzYsVzHeKqm2NqcnJxMAF64xbeoqIgSEhIoICCARCIR6evrq4rduu23jxswYIAqi1AopM6dOz9zdmxhYSHt3LmTli1bRr6+vmRiYlKvkLO0tKR27dqpVoN1UU5ODm3atIl8fHxILBbT559//sItw48ePaL//Oc/z2xGVXdeuW6btqY6cz+urKyMUlNTVd2NLSws6o3rqZs/W1RUpJb73759mwIDA4nH49GMGTOarNt0c3DmzBkaM2aMqukaG+XHNDOswGUYRrtUVlaSXC4nGxsbMjIyorCwMLp16xbXsZgG+Oqrr4jH4zV4vIu2UygUxOfz6cqVK1xHeSmvurV5y5YtJBaLX+ked+/eVRW7YrGY9PT0KCAggORyOWVkZFB2dna9rbR1W5INDAxeqkhVKpV08eJF1WgXS0tLtc2C1bTa2lpatWoVCYVCevPNN59ZkBUVFVHfvn1f2JSLx+Np9GczPz+fFAoFhYWFqYr1x19Qed64nqZW1yHZxsaGnbV9zPXr10kikZBAIKCuXbuykT9Mc8UKXIZhtNPDhw8pJiamXifHxp7DYjQvMjKS+Hy+Tq+wPUtNTQ05Ojq+8lgWbfK8rc0dO3YkQ0PDBndtvnPnDm3atImGDx+uOiMqEAieOwLn7bfffu5qY3V1Nc2cOZOEQiHFxsY2yyfnaWlp1L59e/Lw8KjXjIvon5W39u3bP/eM8uMvHERGRlJycjJlZGQ0acbS0lI6ePAgRUdH09ixY6lt27aq8UV9+vShefPm0bZt2zR+3OTWrVs0fvx41Vzbp+0gaIkKCwtVs2wdHBxILpdTTU0N17EYRl1YgcswjHarm8XXunVrMjExIalUyraa6ZgFCxaQSCRqlrMmV69eTSYmJs3mZ/Lxrc3e3t4kFoubpGvzvXv3yNvb+4VFmVAopG7dulF2dvYT16iqqqLRo0eTkZFRs/xZetz169fJzc2NOnToQLm5uURE9H//939kZGT0wpXbx//o6ekRAPr0008bnKWmpobOnDlDGzdupJCQEPLw8FCtwrdr147GjBlD0dHRlJqayukcbIVCQZaWluTg4EApKSmc5dAmd+/epcjISDIxMSFbW1s2y5ZpKX7kERGBYRhGy5WWlmLdunWQyWQQCARYuHAh5s6dC0NDQ66jMS9ARJBIJPj++++xe/duDBo0iOtITaa4uBh2dnb4+OOPER4eznWcJrVw4UKkpqYiJSUFmZmZOHnyJI4cOYI///wTt2/fhkAggIuLC7y8vODl5YX+/fuje/fu4PP5T71emzZtUFhY+ML7ikQiCIVCbN68GZMnTwbwz8/Q1KlTsWPHDuzduxd9+/Zt0seqje7duwd/f3+UlpbinXfewdKlS6FUKl/5OjweDxMmTIBCoXipjy8oKEB6ejpOnjyp+p7fv38fIpEIHh4e8PHxUX3Pu3bt+sp5mtqtW7cwe/Zs/PrrrwgJCcFnn30GExMTrmNx6uHDh1i7dm29/y/DwsJgYGDAdTSG0YStrMBlGEan3Lt3D19++SW++OIL6Onp4YMPPmD/ceuA2tpa/Oc//8GePXvwxx9/wNvbm+tITWbOnDn4/fffkZ2dDYFAwHWcJhMSEoK8vDzs2bPniff9/fff9Qqgo0ePory8HMbGxnjttddUBdDAgQPh4OCAgoIC2NjYvHKGt99+G3FxcVixYgXWrFmDnTt3wt/fvykenk7Iz89H165dUVJSAgDQ09ODSCSCSCRSvbjXunVrAICJiQlKS0tx5coVVFZW4vGnd7a2tsjLy3vi+tXV1cjKysLhw4dV38/z588DADp27FivmO3Zsyf09PTU/ZBfybZt2xAaGopWrVph48aNGDJkCNeROFVVVYX4+HhERkairKwM7733Hj766COYmppyHY1hNIkVuAzD6KY7d+7gs88+w5dffgkLCwssWLAA7777rtY9AWP+p6qqCmPGjMGJEydw6NAhuLm5cR2pSWRnZ6NLly745ZdfMHr0aK7jNJnAwEAIBAJs3br1hR9bXV2NM2fO4NixYzh+/DiOHz+OixcvgojQoUMH2Nra4ujRow3K4ejoiGvXrmHTpk2YNm1ag66hq2pqapCZmQlfX19MnToVsbGxT/243NxczJkzB7t37wafz39ipZfH46GoqAglJSX1itn09HRUVlbCzMwM3t7eqoK2X79+sLCw0MRDbJCCggKEhoYiKSkJISEh+Pzzz2FsbMx1LM4olUr8/PPPiIiIwM2bN/HOO+9gxYoVaNOmDdfRGIYLrMBlGEa3FRYWYs2aNYiNjYW1tTUWL16M6dOnQygUch2NeYry8nK88cYbyM3NRWpqKhwdHbmO1CQCAgJQXl6O/fv3cx2lyfj7+8PJyQnr169v0OeXlJTg+PHj+Ouvv5CQkIArV6680ueLRCLo6enh4cOHMDMzw+7du9G7d+8GZdF1P/74I6ZMmYLff/8db7zxhurtFRUViI6OxqpVq0BEqK6ufuY1DA0NUV5eDgMDA/To0QO9evVC79690adPH9jb22viYTSJbdu24d1330Xr1q2xceNGvP7661xH4gwRYfv27Vi6dClyc3Mxbdo0LFu2DO3bt+c6GsNwaStrMsUwTLNw48YNCgsLY10idUBxcTH16NGDnJycms2s47179xKAJu9WyyVvb2+SSqVNcq1+/fo9t3sy/tVoys7Ojnbt2kXvv/8+tW3bttk08WqMiRMnkouLi6pJUFJSEtnZ2T0xdulpf8RiMY0aNYrS09OpqqqK40fSMNeuXSN/f3/i8/kkkUiorKyM60icSklJIS8vL+Lz+RQYGEiXL1/mOhLDaIsfn94JgmEYRsfY2dkhNjYWly5dwtChQ/Hee+/B3d0diYmJDWrMwqiPqakp9uzZA7FYjKFDh+Lu3btcR2o0f39/eHh44Msvv+Q6SpMpLi5ukrN7SqUSmZmZ9c6EPk9NTQ3y8vIwYsQIfPnllwgPD2dnCAHIZDJcu3YNK1euxPDhwzFq1Cjk5+ejtrb2hZ9bU1MDpVIJLy8viEQiDaRtOkSEuLg4uLu74+rVqzhw4ADkcjmMjIy4jsaJtLQ0DBkyBP7+/mjdujXS09OhUCjg7OzMdTSG0RqswGUYplmxt7eHXC7HmTNn4O3tjenTp+O1117Dtm3bXvoJNqN+VlZW2Lt3L8rKyvDmm2+itLSU60iNNmfOHPzwww+4desW11GaRHFxMczMzBp9nYsXL+Lhw4dPfR+fz4dIJKrXeVkoFMLOzg7t27eHsbEx9PT0nrv1VpM2b96MiRMnYsmSJQgJCcGPP/6osXu3bdsWHh4eiIqKwp49e57ZrfpplEol/vrrLzWmU4+rV6/C398f7733HmbPno0zZ85g4MCBXMfixNmzZzFx4kT069cPVVVVOHToEFJSUtC9e3euozGM1mFncBmGadbOnTuH5cuXY/v27XB3d8eSJUsQGBjIdSzm/7ty5QoGDBgAV1dX/P7779DX1+c6UoNVVlaiQ4cOmD17NiIjI7mO02hisRgJCQl46623GnWd+Pj4es2hBAIB2rRpAycnJzg7O8PBwQEODg5wdHSEg4MDbGxsUFNTg/bt22P+/Pn46KOPGvtQmsR///tfbN68GRkZGTAzM0NxcTG6d++OefPmISwsTO33VyqV+OOPPzB06FAsXLgQIpEIubm5uHz5MnJyclSdlnk8HvT09MDj8VBRUVHvGvn5+Q3qZq1pRIQNGzbggw8+gL29Pb799lv07NmT61icuHTpEj755BN8//338PLywtKlSzFy5EiuYzGMNmNncBmGaRmysrIoMDCQeDwe9e3bl1JSUriOxPx/WVlZZG5uTqNGjVKdL9RVixcvJmtra6qoqOA6SqOUlpYSAPr9998bfa39+/dTfHw8HTx4kK5du/ZSZ+N37txJPB6Pbt682ej7N4UbN26QSCSiVatW1Xv7J598QoaGhnTnzh2NZfH29iaJRPLE24uLi+nUqVO0fft2Wr16Nc2aNYsGDx5Mtra2qnO6v/32m8ZyNlROTg4NHjyYhEIhSaVSevToEdeROHHjxg2SSCQkFArJ1dWVFAoFKZVKrmMxjC74kRW4DMO0KGlpaRQQEEAAyMfHhw4cOMB1JIb++b4YGxvTlClTqLa2lus4DZafn09isZji4+O5jtIoeXl5BICOHj3Kyf2lUil17dqVk3s/zcqVKwkAHT9+vN7b09LSCABFR0drLEtERAS5ubm90udUV1dTTk4O5efnqylV49XU1NDnn39OhoaG9Nprr9GpU6e4jsSJgoICCg8PJz09PXJ0dKSEhATWMJFhXg1rMsUwTMvSp08fJCcn48iRI9DT08PgwYPh7++PEydOcB2tRevTpw9+/fVXbN++HXPnzuU6ToPZ2NggMDDwmfNKdUVxcTEANMkZ3IY4duwY+vbtq5F7bd++HRYWFuDxeFi6dKnq7d988w0EAgE2bNiAw4cPAwBsbW3rfa6dnR0AIDMzUyNZAaBfv364cOEC7t+//9KfIxQK0bFjR63dnnz27Fn069cPH330EaRSKU6cONHizpYWFRVhwYIFcHJywrZt2/D555/j4sWLCA4OhkAg4Doew+gUVuAyDNMi9evXD3/88QdSU1NRXV2NXr16wd/fH6dOneI6Wovl6+uLrVu3Ii4uDkuWLOE6ToPNmzcPGRkZOHToENdRGozrAvfatWtwcXHRyL0mTJiA5cuXAwB8fHxUbw8ICMDkyZMREhKCv//+GwDQunXrep9rbm4O4J9mSJrSuXNnEBHy8vI0dk91qa6uRnR0NLy9vVFdXY1jx45h2bJlOtfpuTHu3r2LqKgodOrUCQkJCVi2bBmuXLmC9957D2KxmOt4DKOTWIHLMEyL1r9/fxw8eBApKSkoLi6Gt7c3Ro4cidOnT3MdrUUaM2YMNm/ejFWrVmH16tVcx2kQLy8v+Pj46PQqLtcF7t27d2FhYaGx+82aNQsdOnTAN998o3rbhg0bsHDhQgBAq1atAPzTwOlxdX+vqqrSUFKovi537tzR2D3V4fTp0+jTpw+WL1+O5cuXt7hV29LSUkRHR8PJyQlr167FvHnzkJOTA6lUCgMDA67jMYxOYwUuwzAMAD8/P5w4cQJ79+5Ffn4+vLy8MHHiRFy+fJnraC1OUFAQYmNjERERAblcznWcBgkPD8dvv/2GnJwcrqM0SHFxMcRiMWdPtCsqKjR6b5FIhPDwcOzcuRO5ubmorq7GpUuX4OnpCQDo0qULgP8V/nXqtglrcutv3fzXZ41e0nYVFRWIiIiAt7c3DA0NkZGRAalU2mK24ZaVlSE6OhodOnTA6tWr8f777yMnJwdRUVGqF1IYhmkcVuAyDMM8xs/PDydPnsTWrVtx5swZuLq6YuLEibhy5QrX0VqUOXPmYNmyZZg9eza2bt3KdZxXNm7cONjb2+Prr7/mOkqDFBcXP7EdV5PMzMxe6YxpU5g5cyaMjIywdu1a/Prrr5gwYYLqfV27dgUA1VblOgUFBQD+2QmiKffu3QPwv+3RuuTw4cPo0aMHvvnmG3z++ec4dOiQxraic62usLW3t8cnn3yCWbNmqQpbU1NTruMxTLPCClyGYZh/4fF4CAwMxLlz57B161acPn0abm5uCA4O1uhZu5YuKioK8+fPR3BwMHbt2sV1nFciEAgQGhqKjRs3quaT6pLi4mLOticDgKWlJYqKijR6z1atWmHmzJnYvHkzfvrpJ4wdO1b1vqCgIJiZmeHAgQP1Pmf//v0Qi8X4z3/+o7GcdV8XS0tLjd2zscrLyxEREYFBgwahY8eOOHv2LMLDw8HnN/+noQ8fPkRsbCw6deqEjz/+GCEhIbh+/TpkMhmnv2MM05w1/39ZGIZhGojP56sK3Y0bN+Lo0aPo0qULZs2a9cRKDqMeq1evxtSpUxEYGKhzTZtCQkJARIiPj+c6yisrKSnh9Ml3165dkZGRofH7hoWFoaysDN27d4dQKFS9vXXr1vjoo4+wfv16lJWVAfjnDGVdQ7R/d1dWp1OnTsHAwACOjo4au2dj7NmzB66uroiLi8M333yDXbt2qbpPN2eVlZWIi4uDs7MzFi9ejEmTJiEnJwcymYzT3REM0xKwApdhGOYFRCIRgoODceHCBXz11VfYtWsXnJ2dER4ejlu3bnEdr1nj8XhYv349Ro4ciVGjRiE9PZ3rSC/NzMwMU6dORWxsLGpra7mO80q4XsHt168f0tLSNH5fBwcHzJ07F6GhoU+878MPP0RERARmz56NJUuWYMaMGVi4cGG90UKakJaWhp49e2p9h93i4mLMmjULw4cPR+/evXHp0iVIJBKuY6ldVVUV4uLi0LFjR8ybNw+BgYG4cuUKYmNjYW1tzXU8hmkReEREXIdgGIbRJVVVVYiPj0dUVBQePHiAOXPmQCqVslfl1aiqqgpjxozBiRMncOjQIbi5uXEd6aVkZ2ejS5cu+OWXXzB69Giu47y0wMBACAQCzs4/nz59Gt27d8fRo0c1Ng9XF9TU1MDOzg6hoaFYtmwZ13GeKTk5Ge+++y6USiW+/vprjBs3jutIalf3/8KKFStQVFSEd955B1FRUWjXrh3X0RimpdnKVnAZhmFekVgshkQiQXZ2Nj755BPEx8fD3t4eERERT3RZZZqGWCzG9u3b4erqiqFDh+rMWWhnZ2cMHz5c50YGcb2C6+npCU9PT2zevJmzDNpo586duH37Nt5++22uozzVrVu3MGHCBIwePRq+vr44d+5csy9uq6urkZiYCFdXV8ydOxcjRozA1atXIZfLWXHLMBxhBS7DMEwDGRkZITw8HFeuXMHixYshl8vh5OSkWtllmpahoSGSk5PRpk0b+Pv7qzrYarvw8HAcOHBAp2Yrc13gAsD06dPx448/smMAj4mJiYG/vz86duzIdZQnbNu2Dd26dcOpU6ewd+9eJCYm6mSn55elVCqxbds2uLm5YebMmfDz80NOTg7kcrlGx0YxDPMkVuAyDMM0krGxMaRSKW7cuIEPP/wQMTExcHJyQnR0NCoqKriO16yYmppiz549EIvFGDp0KO7evct1pBfy9/eHh4cHvvzyS66jvLTi4mLOR5eEhITAwsICUVFRnObQFjt37sShQ4e0bmvytWvXMGzYMEyePBnjx49HVlYW/Pz8uI6lNnWFraurK6ZMmYK+ffviwoULkMvlGm02xjDMs7ECl2EYpomYmJhAKpUiJycHM2bMwIoVK+Dg4IDo6Gg8evSI63jNhpWVFfbu3YuysjK8+eabKC0t5TrSC82ZMwc//PCDzqxGasMKrr6+PlasWIFNmzbh5MmTnGbhWnl5ORYuXIhx48bBx8eH6zgAACJCXFwcPDw8kJ+fjyNHjkAul8PY2JjraGrx+Irt5MmT8dprr+H8+fNITEyEk5MT1/EYhnkMK3AZhmGamIWFBWQyGa5du4Zp06Zh+fLl6Ny5M2JjY1FZWcl1vGbB1tYWKSkpuHHjBkaPHq31LyAEBwfD1NQUcrmc6ygvhesxQXWCgoLw+uuvY/LkyTrxQoa6zJs3D7du3cKaNWu4jgIAOH/+PPr37485c+Zg3rx5OHnyJPr06cN1LLVQKpXYunUr3N3d8dZbb6FXr164ePEiFAoFOnXqxHU8hmGeghW4DMMwamJlZQWZTIbLly9j7NixkEqlcHFxQVxcHGpqariOp/M6deqEvXv3IjMzE5MmTdLqr6menh5CQkKwbt06rS/Gy8rKUF1drRUFLp/PR0JCAoqLiyGRSKBUKrmOpHFbtmzBhg0bsGnTJtjb23OapbKyElFRUejevTuqq6uRnp6O5cuXQ09Pj9Nc6lBTU4PExES4ubnh7bffhqenJ86ePYvExEQ4OztzHY9hmOdgBS7DMIya2draIjY2FpcvX8awYcPw3nvvwd3dHYmJiTo3H1XbuLu7Y9euXdi/fz/eeecdrS6AZs+ejeLiYvz0009cR3muuk7g2lDgAoCNjQ22bt2KHTt24IMPPuA6jkbt2bMHM2bMwIcffsh5N+KjR4/Cy8sLn376KVasWIG0tDR4eHhwmkkd6roid+3aFdOmTYOHhwfOnj2L77//Hl26dOE6HsMwL4EVuAzDMBrSoUMHyOVyZGdnY+DAgZg+fTo8PDywbds2sJHkDdenTx/8+uuv2L59O+bOnct1nGeysbFBYGCg1o8M0rYCFwB8fX2RmJiIL7/8EosWLWoRvy979+7FhAkTMGXKFKxatYqzHA8fPkRERAQGDhwIOzs7nDt3DlKpFAKBgLNM6lBVVaUa9zNz5kz07t1btRWZFbYMo1tYgcswDKNhDg4OkMvlOHPmDLy8vFQNS1ih23C+vr7YunUr4uLisGTJEq7jPNO8efOQkZGBQ4cOcR3lmbSxwAWAiRMn4ttvv8Vnn32G6dOno7q6mutIarNlyxYEBARg7Nix2LBhA3g8Hic5du3aBTc3N8TFxWHdunXYvXs3HBwcOMmiLpWVlYiLi0PHjh0REhKCfv364cKFC2wrMsPoMFbgMgzDcMTV1RWJiYnIzMxEly5dMGnSJPTt2xfJyclcR9NJY8aMwebNm7Fq1SqsXr2a6zhP5eXlBR8fH61exdXWAhf4p1lXUlIStm/fDl9fX9y8eZPrSE2qqqoKCxYswNSpUzF//nwkJiZCKBRqPMft27cRHByMgIAA9O7dG5cuXYJEItF4DnUqKytDbGwsHB0dMX/+fIwfPx65ubmsKzLDNAOswGUYhuFYt27doFAocPr0aXTo0AGjRo2Cj48P9u/fz3U0nRMUFITY2FhERERobcfi8PBw/Pbbb8jJyeE6ylMVFxdDLBbDwMCA6yhP9cYbb+DIkSO4c+cOPD09sWPHDq4jNYns7Gz0798fcrkc8fHxkMlkGl+5JSLV+dPU1FTs2bMHCoUCVlZWGs2hTqWlpYiOjoa9vT2WLFmCwMBAZGdnIzY2Fu3bt+c6HsMwTYAVuAzDMFrCw8MDCoUCR48ehbm5OXx9fdG/f3+t3s6qjebMmYNly5Zh9uzZ2Lp1K9dxnjBu3DjY29vj66+/5jrKUxUXF6N169Zcx3guDw8PnDhxAmPGjMG4ceMwZswYXL9+netYDVJRUYHIyEh4eHigpqYG6enpCA4O1niO3NxcDB06FNOmTcP48eNx5swZDBs2TOM51OXu3buIioqCvb09PvnkE4SEhOD69euIjY1Fu3btuI7HMEwTYgUuwzCMlqnbppyamgqxWIzXX38d/v7+SE9P5zqazoiKisL8+fMRHByMXbt2cR2nHoFAgNDQUGzcuBElJSVcx3lCcXGxVm5P/jcjIyNs3LgR+/btw6VLl+Dm5obFixfj7t27XEd7KTU1NUhISEDXrl3xxRdf4JNPPsHx48fh4uKi8RyxsbHw8PBAYWEh0tLSIJfLYWxsrNEc6lJUVISoqCg4OTlh7dq1CAsLw40bNyCTyWBubs51PIZh1IAVuAzDMFqqf//+2L9/P1JSUvDgwQP07NkT/v7+yMjI4DqaTli9ejWmTp2KwMBArVsFDwkJAREhPj6e6yhPKCkp0YkCt46vry8yMzOxfPlybNiwAY6Ojli0aBHy8/O5jvZUFRUV2LRpE11Bf6gAACAASURBVFxdXRESEoIhQ4bg4sWLmD9/vsbP254+fRp9+/ZFREQEFixYgBMnTqBXr14azaAuhYWFiIiIgIODA9atW4f3338fOTk5iIqK0qmfb4ZhXh0rcBmGYbScn58f/vrrL6SkpOD+/fvw8vLCyJEjkZWVxXU0rcbj8bB+/XqMHDkSo0aN0qoVcDMzM0ydOhWxsbFaNwtZV1ZwHycWi7FgwQLk5uZiyZIl2Lx5MxwcHDBmzBj8/vvvqKmp4Toizpw5g/DwcNjY2GD27NkYOHAgLl68iI0bN8LGxkajWSoqKhAREQFvb28YGBggIyMDUVFREIvFGs2hDjdu3EB4eDgcHBwQHx+PZcuW4dq1a4iKioKpqSnX8RiG0QBW4DIMw+gIPz8/nDhxAr/99htu3ryJ7t27Y+LEicjOzuY6mtYSCATYsmULfHx8MHz4cJw/f57rSCrh4eG4fv06du7cyXWUenSxwK1jbGyMDz/8EDdu3MB3332HkpISBAQEoE2bNpg2bRqSk5NRVlamkSy1tbU4ceIEFi1ahC5dusDDwwO7du3CRx99hLy8PGzatAkdO3bUSJbHHTp0CJ6enpDL5fj8889x8ODBZjHn9dq1awgPD4eLiwt+/fVXrFq1ClevXoVUKoWhoSHX8RiG0SAesaGLDMMwOkepVOLnn3/G0qVLkZubi8mTJyMyMpKNt3iG8vJyvPHGG8jNzUVqaiocHR25jgQACAgIQHl5uVZ1zPb394eTkxPWr1/PdZRGO3/+PHr16oXJkyfj/Pnz+Ouvv8Dn8+Hl5YUBAwagd+/e6NatGzp16tTo7cH5+fk4d+4cTp48idTUVBw5cgQPHjyAo6Mjxo0bh3HjxqFPnz7g87lZW7h//z4iIiKwYcMGjBgxAuvXr28WXYNzc3MRHR2NzZs3w87ODuHh4Xj33Xehp6fHdTSGYbixlRW4DMMwOqyu0F20aBGuX7+OadOmYdmyZc3iiWtTKykpwZAhQ1BSUoLU1FSt6JyakpKCoUOHIiMjA56enlzHAQD07NkTvr6+kMlkXEdplIcPH6JXr14wMTHBn3/+CbFYjFu3buHQoUNITU3FoUOHcOHCBdTW1kJPTw+dO3eGnZ0d2rZtC1tbW5iamsLY2BgikQjGxsaorKxEeXk5KisrUVJSgtu3b+PmzZsoKCjA5cuXcf/+fQBA+/btMWDAAAwYMACDBg1C165dOf5KAN999x3mz58PPT09fP311xg1ahTXkRotKysLK1euxLZt2+Di4oJFixbhrbfegkAg4DoawzDcYgUuwzBMc1BdXY0ff/wRy5cvx82bN/HOO+8gKipKK4o4bVJUVIRBgwZBIBDg4MGDsLCw4DoSXnvtNXh5eWHz5s1cRwEAODs7Y/r06fjoo4+4jtIodR20T548CQcHh6d+zKNHj3D+/HmcO3cOly5dUhWs+fn5KC0tRWlpKaqrq1FWVgY9PT0YGhpCX18fJiYmaNOmDWxtbdG2bVt06tQJXbt2hbu7u1Z15r1y5QpCQ0Oxf/9+hIaGYuXKlWjVqhXXsRrlzz//RHR0NHbv3g13d3csXrwYEyZM4GxlnGEYrcMKXIZhmOakqqoK8fHxWL58OUpKSjBjxgwsWrQIbdq04Tqa1rh58yYGDBgAa2tr7Nu3DyYmJpzm2bBhA+bOnYtr166hbdu2nGYBACsrK6xYsQKhoaFcR2mwdevWYe7cudi1axfeeOMNruNoXHV1NdasWYOoqCh07twZ69evR9++fbmO1WBEhJ07d0Imk+Ho0aPw8fGBVCpFQEAAeDwe1/EYhtEuW9nLXQzDMM2IWCyGRCJBbm4u1qxZA4VCgU6dOiEiIkK1hbKls7W1RUpKCm7cuIHRo0fj0aNHnOYJDg6Gqakp5HI5pznq6NqYoH/LzMzEggULsGTJkhZZ3KampsLT0xMrVqyAVCrFiRMndLa4ra6uRmJiIrp164bRo0fD3NwcR44cweHDhzFy5EhW3DIM81SswGUYhmmG9PT0IJFIcOXKFSxZsgRxcXGwt7dHREQESkpKuI7HuU6dOmHv3r3IzMzEpEmTOB0jo6enh5CQEKxbt47zYrusrAzV1dU6W+Dev39f1cxp2bJlXMfRqPv372PWrFkYNGgQOnbsiHPnzuns6J+ysjLExsbCyckJM2fOhJeXF86ePYvk5GT069eP63gMw2g5VuAyDMM0Y0ZGRpBKpbhx4wYWL14MuVwOJycnREVFobS0lOt4nHJ3d8euXbuwf/9+vPPOO1AqlZxlmT17NoqLi/HTTz9xlgH4Z0QQAJ0scIkI06dPR3l5OX744YcW02yIiJCYmAgXFxckJycjPj4eycnJzzx3rM2KiooQFRUFe3t7LFmyBGPHjkVOTg4SExPh5ubGdTyGYXQEK3AZhmFaAGNjY0ilUuTk5GDOnDn44osv4OTkhOjoaFRUVHAdjzN9+vTBr7/+iu3bt2Pu3Lmc5bCxsUFgYCBiY2M5ywDodoEbHR2NnTt3QqFQaMVZZk3Izs6Gv78/pk2bhrFjx+LixYsIDg7mOtYru3r1KsLDw+Hg4KA6P339+nXExsbCzs6O63gMw+gYVuAyDMO0IObm5oiKikJOTg6mT5+O5cuXo3PnzoiNjUVlZSXX8Tjh6+uLrVu3Ii4uDkuWLOEsx7x585CRkYFDhw5xlkFXC9xDhw5h6dKlkMlkGDBgANdx1K6iogJRUVFwd3fHvXv3kJaWBrlcrnMdkjMzMxEcHIzOnTsjKSkJK1euxLVr1xAVFaVV3agZhtEtrIsywzBMC1ZYWIg1a9YgNjYW1tbWWLx4MaZPnw6hUMh1NI3bsmUL3nnnHaxatQoffvghJxn69+8Pa2tr/PLLL5zcf+fOnRg5ciTKy8thYGDASYZXdfv2bfTo0QNeXl747bffmn3joYMHDyI0NBR5eXlYunQpFixYoFPbsYkI+/btw2effYa9e/fC09MTUqkUgYGBOvU4GIbRWqyLMsMwTEtmbW0NmUyGy5cvY8yYMQgLC4OzszPi4uJQW1vLdTyNCgoKQmxsLCIiIjjraBweHo7ffvsNOTk5nNy/uLgYYrFYZ4pbpVKJoKAgGBgYIDExsVkXt7du3UJwcDAGDx6MTp064cKFC5BKpTpTFFZVVSEhIQGenp4YOnQoamtrsWfPHmRkZGDy5Mk68zgYhtF+rMBlGIZhYGdnh9jYWFy6dAlDhw7Fe++9B3d3dyQmJraoQnfOnDlYtmwZZs+eja1bt2r8/uPGjYO9vT2+/vprjd8b+KfAbd26NSf3boglS5YgNTUVCoVC57ZVv6y6JlLdunXD/v378fPPPyM5OVlnzqY+ePBA1RE5JCQELi4uOHbsGPbt24dhw4ZxHY9hmGaIFbgMwzCMir29PeRyOS5fvowBAwZg+vTpeO2117Bt2za0lBMtUVFRmD9/PoKDg7Fr1y6N3lsgECA0NBQbN27kZJxTcXGxzhSKv//+O6Kjo/HVV1+hR48eXMdRi6ysLPj4+GDGjBmYMmUKLl68iHHjxnEd66Xk5uYiPDwcNjY2WLZsGcaNG4ecnBwoFAr07t2b63gMwzRjrMBlGIZhnuDo6Ai5XI6srCy4ublh0qRJ8PT0xLZt27iOphGrV6/G1KlTERgYqPGmTyEhISAixMfHa/S+AFBSUqITBe6NGzcwdepUTJ48GTNnzuQ6TpMrLy9HVFQUevbsCT6fj4yMDMTGxsLY2JjraC908uRJBAcHq8YWLV26lHVEZhhGo1iByzAMwzyTm5sbFAoFMjMz4eLigkmTJqFv375ITk7mOppa8Xg8rF+/HiNHjsSoUaOQnp6usXubmZlh6tSpiI2N1fj2cF1Ywa2ursbkyZPRtm1bbNiwges4TS4pKQmurq748ssvsXbtWqSmpqJbt25cx3oupVKJ5ORk+Pj4wNvbG+fPn8emTZtw+fJlSKVSrf+ZYhimeWEFLsMwDPNC7u7uUCgUSEtLg6WlJUaNGoX+/fvjwIEDXEdTG4FAgC1btsDHxwfDhw/H+fPnNXbv8PBwXL9+HTt37tTYPQHdKHDnzZuHM2fOQKFQwNDQkOs4Tebq1asYNWoUxowZg4EDB+LixYsICQnR6sZZZWVliIuLg6urK8aMGQNzc3OkpKQgPT0dwcHBLbIbO8Mw3GMFLsMwDPPSevfujeTkZBw5cgT/r727j6v5bvwH/jqnTjeS0DC3WW7qJBbF6IaNQizCVXEhLlbR2FlZq+Gifd27fB86jK3MXSOtXOGbhNws2eamu6WrdEuhZHMpRffn/P7w02MmhNP5FK/n47E/ds7p/X6d9s9evd+f91tHRwejR4+Gg4MDLl26JHS0ZqGlpYWDBw9CKpVi7NixuHbtmlrm7devHxwdHSGXy9Uy32MtveBGRERg27Zt+PbbbyGVSoWOoxK1tbWQy+UYNGgQsrKycPz4cfzwww/o3Lmz0NGe6fbt2wgMDISRkRGWLFmCDz/8EP/5z38QHR0Ne3t7oeMR0VuOBZeIiF6atbU1Tp06hYSEBNTW1uKDDz6Ag4MDkpKShI6mcm3atEF0dDS6dOkCBwcHFBcXq2VemUyGs2fPIjU1VS3zAS274GZnZ8PDwwOLFy/GrFmzhI6jEmfPnoWFhQWWLl2KJUuWIC0tDWPHjhU61jNdvHgRM2bMQK9evbBjxw74+fmhsLAQwcHBMDU1FToeEREAFlwiInoNtra2+OmnnxAXF4eysjIMHToUTk5Oai1l6mBgYIDjx49DS0sLY8eOxd27d5t9TgcHBwwaNAhbtmxp9rkeKy0thYGBgdrma6rKykq4urrC1NQU//rXv4SO89qKi4vh7u6O0aNHw9jYGBkZGQgMDIS2trbQ0Z5SW1uL8PBwDB8+HMOHD0dOTg527tyJa9euISAgoFVdK0VEbwcWXCIiem329va4dOkSTp48iaKiIlhaWsLV1RVZWVlCR1OZTp064eTJk6ioqMCECRNQXl7e7HMuWrQIYWFhuH37drPPBbTcFdyFCxeisLAQ4eHhLbIENlVdXR3kcjlMTU3xyy+/ICYmBtHR0TAyMhI62lPKysogl8vRt29fzJw5E506dWp4vnb27NnQ0tISOiIRUaNYcImISGXs7e2RmJiIw4cPIzs7G2ZmZnB1dUVubq7Q0VSiR48eiIuLQ2FhISZPnoyqqqpmnc/d3R0GBgYIDg5u1nkea4nXBAUHByM0NBS7d+/Ge++9J3ScV5aQkIAhQ4YgICAAPj4+SE9Px4QJE4SO9ZTs7Own7q91dnZGfn4+n68lolaDBZeIiFRKJBLByckJycnJCA8PR2pqKszMzODu7o78/Hyh4722vn374uTJk/jtt9/g5uaGurq6ZptLW1sbHh4e2L59e7OX6YqKCtTW1raogpuWlgYfHx8sXboUkydPFjrOK7l79y68vLwwatQodO7cGSkpKQgMDISOjo7Q0RooFAqcOnUKTk5OMDU1RWxsLNauXYuioiLI5fIWucJMRPQsLLhERNQsxGIxXFxc8J///Afff/89fvnlF0ilUnh5eaGoqEjoeK9l4MCBiImJwZkzZzB37lwoFIpmm8vb2xulpaX48ccfm20O4NH2ZAAtpuCWlpZi6tSp+OCDD/D1118LHeelKRQKhIaGwsTEBDExMdizZw9OnTrVog5jKi8vR0hICMzNzTF27FhUVVXhyJEjyMrKgkwmg56entARiYheGgsuERE1K4lEAnd3d2RmZmLr1q2IiYmBsbExvLy81PZsaXMYPnw4Dh8+jIMHD2Lx4sXNNk+3bt3g4uLS7FcGtaSCq1QqMX/+fFRUVCAsLAwaGhpCR3opycnJsLa2xieffIKZM2ciMzMT7u7uQsdqkJ+fj4CAABgZGUEmk8HKygpXrlxBXFwcnJycWvTdu0REL8KCS0REaiGRSODp6Yn8/Hxs2bIF0dHR6Nu3LwICAvDf//5X6HivZMyYMQgPD0dISAiWL1/ebPP4+PggJSUF8fHxzTZHSyq4mzZtwv/93/8hIiICXbt2FTpOk927dw8ymQzDhg2Djo4OUlJSIJfLoa+vL3Q0KJVKnD59GpMnT0a/fv0QHh6OpUuXori4GKGhoRgwYIDQEYmIVIIFl4iI1EpLSwuenp7Izc3FmjVrsGfPHhgZGSEgIKChZLUmzs7O2LVrF9atW4eNGzc2yxyWlpawsbFp1lXcllJwL1y4gGXLlmHNmjUYOXKkoFmaSqlUNmxHjoyMxK5du3D27NkWURorKiqwfft2DBgwAPb29igtLUVERATy8vLwxRdfCP7fm4hI1VhwiYhIEG3atIFMJkNubi6WL1+OkJAQ9OnTB4GBgbh//77Q8V7K7NmzIZfLERAQ0GwnHstkMhw5cgR5eXnNMn5paSm0tLSgq6vbLOM3xZ07d/C3v/0N48aNg5+fn2A5XsZvv/0GW1tbzJ8/HzNmzMDVq1fh7u4u+Dbf3Nzchm3Ivr6+sLKyatgFMG3atFa37ZuIqKlYcImISFBt27aFv78/CgoK8OWXXyIoKAh9+vTBhg0bUFlZKXS8Jlu0aBFWrFgBb29vhIeHq3z8qVOnwsjICNu2bVP52MCjgtuhQ4dmGbspFAoFZs+eDU1NTezZs0fwgvgipaWlkMlksLS0hEgkQmJiIuRyOdq1aydYpsenIbu6usLU1BT79u3D4sWLcfPmTYSGhsLCwkKwbERE6sKCS0RELYK+vj78/f2Rl5eHTz/9FGvWrIGRkRE2bNjQ7FfkqEpgYCB8fX3h7u6OmJgYlY6toaGBhQsX4vvvv0dZWZlKxwYeFTYht6uuXLkS8fHx+Pe//w1DQ0PBcryIUqnEnj17YGpqigMHDmDHjh1ISEjA+++/L1im+/fvN5yG7ODggKKiIhw4cADXr19HYGAg3nnnHcGyERGpGwsuERG1KIaGhggMDEReXh7mzZuHr7/+Gv3794dcLkd1dbXQ8V5o48aNmDNnDlxcXFR+KJSHh0dDwVK1srIywQru6dOnsW7dOgQFBcHS0lKQDE2RkpICOzs7zJ8/H2PHjkVGRgb+8Y9/CLba/Pg6n27duuGLL76AnZ0drly5gvPnz8PFxQWampqC5CIiEhILLhERtUidOnXC+vXrkZ2djSlTpiAgIAAmJiYICQlBXV2d0PGeSSQS4bvvvoOTkxMmTZqExMRElY3dvn17zJkzB3K5HPX19SobFxBuBffGjRuYPn06XFxcsGDBArXP3xSPT0ceOnQoampqcOHCBYSGhgqyMqpQKBAdHQ0HBwdIpVLExsbin//8JwoKChAcHAxzc3O1ZyIiaklYcImIqEXr0aMH5HI5srKyMG7cOHz66afo378/QkJCVF7yVEVDQwM//PADbGxs4OjoiIyMDJWNLZPJUFBQgKNHj6psTECYgltbW4sZM2bA0NAQISEhap27KRQKxVOnI1+8eBFDhw5Ve5Y7d+5g7dq16N27N5ydndGmTRucOHECWVlZ8Pf3F/T5aSKiloQFl4iIWoVevXohODgYOTk5cHBwgLe3NwYNGoTIyEgolUqh4z1FS0sLBw8ehFQqxdixY3Ht2jWVjNuvXz84Ojqq/MogIQruF198gdTUVERFRbWIu2L/LDExETY2NoKfjpyUlAQvLy/07t0bGzduxJQpU5Cbm4sjR47AwcGhxR/GRUSkbiy4RETUqvTu3RvBwcG4cuUKLC0tMWPGjBZbdNu0aYPo6Gh06dIFDg4OKC4uVsm4MpkMZ8+eRWpqqkrGA9RfcCMjI7FlyxZ8++23MDMzU9u8L3L37l3IZDJ88MEH0NHRQUpKitpPR358aNT7778PKysrJCUlISgoCLdu3YJcLsd7772ntixERK0NCy4REbVKUqkUoaGh+O233yCVSuHm5oYRI0YgOjpa6GhPMDAwwPHjx6GlpYWxY8fi7t27rz2mg4MDBg0ahC1btqgg4SOlpaUwMDBQ2XjPk5OTAw8PD3h7e2P27NlqmfNF6urqEBISAhMTExw8eBC7d+/G2bNn1fpMa0ZGBmQyGbp37w6ZTAYTExP8/PPPSExMhKenJ/T09NSWhYiotWLBJSKiVm3AgAGIiIhAamoqevXqhUmTJsHGxganT58WOlqDTp064eTJk6ioqMCECRNQXl7+2mMuWrQIYWFhuH379kv/bGlp6VMHdalrBbeqqgpubm7o3bs3Nm3a1OzzNcW5c+dgaWmJRYsWYebMmQ3bkdWhuroakZGRcHBwwIABA3DixAksX74ct27dQkREBKytrdWSg4joTSFStrT9XERERK/hwoULWLNmDY4ePQobGxusWbMGo0aNEjoWACA3Nxd2dnaQSqU4duwYdHR0Xnms6upq9OrVC97e3li5ciWAR1tbd+/ejczMTHz33XfP/Nnx48fjxIkT0NHRgb6+PgwMDFBZWQkjIyO89957aN++Pdq3bw8DAwPMnTsXnTp1euWcfzVv3jwcPnwYiYmJMDY2Vtm4r6K4uBj+/v7Yt28fRo8eja1bt0Iqlapl7pycHOzcuRM7d+7E/fv3MXnyZHh6emLMmDF8rpaI6NWFQ0lERPQGOn/+vPKjjz5SAlDa29srL126JHQkpVKpVKalpSk7duyonDRpkrK2tva1xlq2bJmyc+fOytTUVOWCBQuUOjo6SgDK999//7k/t3nzZqVIJFICeOofsVislEgkSrFYrHz33XdfO+Offf/990qRSKQ8dOiQysZ8FTU1NcqgoCBlu3btlD169FDu3btXLfPW1dUp4+LilB9//LFSJBIpu3fvrvT391feuHFDLfMTEb0FDnCLMhERvZFsbGxw5swZJCQkoKamBsOGDYODgwOSk5MFzTVw4EDExMTgzJkzmDt3LhQKxSuPZWVlhbKyMlhYWGDXrl2oqqoCAPzxxx/P/blx48Y980AuhUKB2tpaaGho4NNPP4WmpuZLZUpLS8OGDRueGv/KlSv47LPP4O/vD2dn55caU5XOnj2LwYMH46uvvoKPjw9ycnKafTtycXExNmzYAGNjY4wbNw5VVVX48ccfcf36daxfvx49evRo1vmJiN4mLLhERPRGs7W1RXx8POLi4nDv3j1YWVnByckJv/32m2CZhg8fjsOHD+PgwYNYvHjxMz/XWPmtrq5GaGgopFIppkyZ0nAXcE1NTcNn7t2799z5pVIpunfv/tzPKJVKeHh4PPczjYmMjERAQAAmTpzYkKO8vByurq4YOnQoVq1a9dJjqkJBQQGmTZuG0aNHo1+/fsjIyEBgYOBrbRN/HoVCgVOnTsHV1RW9evXC5s2bMWPGDOTm5iIuLg4uLi4v/ccDIiJ6MRZcIiJ6K9jb2+Py5cs4cuQIbt68iSFDhsDV1RXZ2dmC5BkzZgzCw8MREhKC5cuXP/V+cXExPvzwQ1RUVDS8dvbsWbz77ruYN28esrKyAOCpw6IAoLKyErW1tc+d/+OPP4aWllaj70kkEri4uKBLly4v85UAAAcOHAAAxMXFYeDAgbh8+TLmzZuHe/fuISwsTO2lrrKyEhs2bIC5uTnS0tIQExODQ4cOoXfv3s0yX2lpKUJCQjBw4EA4ODggPz8f27Zta1it5RU/RETNiwWXiIjeGiKRCE5OTkhKSkJ4eDjS0tIglUrh6uqKvLw8tedxdnbGrl27sG7dOmzcuLHh9dzcXAwdOhQJCQnYu3dvw+sjR46Evb09RCLRc+/8VSqVL7yOaNy4cc8swbW1tc9dWX6W9PT0ht9jXV0dSkpKYG1tjYsXL+LAgQPo1q3bS4/5OqKjozFgwACsWrUKS5YsQXp6OiZMmNAscyUlJcHLywvdu3eHn58fbG1tkZaW1nDFT3OtFBMR0ZNYcImI6K0jFovh4uKCjIwMhIeHIyUlBVKpFF5eXrh165Zas8yePRtyuRwBAQEIDg7GlStXMGLECNy5cwcA8K9//athq7KGhgbCwsIwZswYSCSS5477ooJrb28Psfjp/w0Qi8UwNzfHiBEjXvq7/Pvf/34iV11dHerq6nDz5k3I5XKUlZW99JivIisrC46Ojpg8eTKsrKyQmZmJwMBAaGtrq3Se8vJyhISEYPDgwbCyskJSUhI2b96MoqIiBAcHY+DAgSqdj4iIXowFl4iI3lp/Lrrff/89Tp06BWNjY3h5eaG4uFhtORYtWoQVK1bA29sb1tbWKCsra1hdLSwsxNGjRxs+K5FIcOjQIVhZWT235L7ooCl9fX0MHz680StpfH19X+l7HDhwoNFVYaVSiWPHjmHIkCG4cuXKK43dFKWlpQgICMCgQYNw584dnDt3DhEREejZs6dK57l69SoCAgJgZGSEzz77DP369UNcXFzDaq2enp5K5yMioqZjwSUioreeRCKBu7s7MjMzsXXrVhw9ehR9+/aFTCZDSUmJWjJ8+OGHkEgkqKqqeqIkisXiJ7YvA4Curi5iY2NhYmLyzJL7ohVcAJg4ceJTz8Tq6enBzc3tpfNnZ2c3PBfcmNraWhQUFGDYsGENz+mqilKpRGhoKExNTbFjxw5s3LgRly5dgq2trcrmqK6uRmRkJBwcHCCVShEVFQV/f3/cunULERERsLe3V9lcRET06lhwiYiI/j8tLS14enoiPz8fmzdvRkREBPr27YuAgIAXnkz8Oo4ePYpx48Y1bOn9s/r6evz888+4fPnyE68bGBjgzJkz6Nmz51MlV1NTs0kFd/z48U+UaYlEAm9vb7Rp0+alv8Nftyc3RqlUwszMDIMGDXrp8Z8lMTER1tbWmD9/PiZPnoysrCzIZDJoaGioZPzMzEwsWbIEPXr0wMyZM9GhQwecPn0aWVlZ8Pf3h6GhoUrmISIi1WDBJSIi+gttbW14enoiNzcXq1evxu7du2FkZISAgACVP0f6ww8/wNnZGXV1dQ1X/vyVRCKBXC5/6vVOnTrhp59+wjvvvPPESqyGhsYLtygDgIWFBd55552Gf6+rq4OXl9crfItnVy8spAAAErZJREFUb08GHuWXSCRYu3YtLl26hAEDBjx3rMLCwhfOV1xcDC8vL3zwwQfQ1dVFcnIygoODn/g+r6qyshKhoaGws7ODmZkZDh06hM8//xwFBQWIiIjA6NGjG93aTUREwmPBJSIiegY9PT3IZDLk5eVh2bJlCA4ORp8+fRAYGIj79++/9vgZGRlYsGABFApFo3fePlZbW4sff/wRRUVFT73Xs2dP/PTTT9DX129YtWzKKcrAo1OlJ0yYAIlEAk1NTUycOPGVrrG5du0a0tPTG31PLBbDwsICaWlp8Pf3f+HK6jfffIPRo0c/ca/vn9XW1kIul8PU1BTHjh3D7t27cebMGZUc6JSZmYmAgAD07NkTHh4eaN++PSIiIpCTk4Nly5aha9eurz0HERE1LxZcIiKiF2jbti38/f2Rl5eHRYsWYfPmzejTpw82bNiAysrKF/78pk2bUFpa+tTrZmZmuH79Or788suGVc5nEYlE2L59e6Pv9e/fHydOnIC2tjbEYjHq6+ubVHABwNHRsWFr9KtcDQQABw8efOpZXolEAi0tLaxduxa//vorTE1NXzhOTEwMZDIZ8vPzG12xPn36NCwsLPDVV19h4cKFyMzMhLu7+ytlfqyqqqrh2VozMzNERUXBz88PN27cQHR0NFxcXFS23ZmIiJqfSPm8i/SIiIjoKX/88Qc2bdqErVu3om3btvD19YVMJmv0rtOsrCyYmZnB0tISZ86cQdu2bRsds7CwEKtWrcLu3bshFosb3e7brl07FBUVPfOU3tOnT8PR0RG1tbVwdHTEsWPHGt6rr69HSUkJSkpKUFpaivr6epSXl+PevXv45JNP0KVLF0RGRqJDhw7o1q0bOnTo0OTfh6WlJVJSUhru5hWLxRg2bBj27t2L/v37N2mM9PR0DB8+HJWVlVAoFNDV1UVeXh66du2K3NxcLF26FJGRkfj444+xZcuWV1pp/rOkpCSEhITgwIEDqKmpwaRJk+Dp6YkxY8Zw+zERUesVzoJLRET0in7//Xf87//+L+RyOTp37oxly5Zh3rx5T6xmTp8+HVFRUVAqlbC2tsaJEycaLcKPXb16FatXr0ZYWBg0NTWfKLoaGhrYtm3bc5+T/fHHHzFjxgz07NkTjo6OSE9PR35+PkpKSp67DfqvdHV10b17d5iYmMDc3Bzm5uYYMmQIpFLpEwXw5s2b6NWrF5RKJSQSCcRiMb7++mv4+fk1es9uY4qLizFkyBD88ccfDYdsSSQSuLm5NayUGxkZQS6XY9y4cU3+Dn91//59hIeHIzg4GMnJyTA1NcXcuXMxf/58lTy7S0REgmPBJSIiel03btzApk2bEBwcjK5du+Krr77C/PnzcfXqVQwcOLBhZVNTUxOjRo1CTEwMtLW1nztmeno6VqxYgUOHDkFTUxN1dXUQiUTo3bs3cnNzG8qjQqHAhQsXcOzYMcTHx+PSpUuoqamBhoYGbG1tYWZmBhMTE3Tt2hXdunVDly5d0LFjR4jFYujr60NTUxPr1q3D3//+d2hpaeG///0vioqKUFRUhBs3biAjIwMZGRnIzMxETU0NOnfuDDs7O3z00UeYNGkSoqKi8PnnnwMA7O3tsXPnTvTq1avJv7uHDx/C1tYW6enpT61ai0Qi6OvrY/Xq1Vi4cOFT26Cb6vFq7f79+1FfXw8nJyeu1hIRvZlYcImIiFQlPz8f//M//4N9+/ahf//+aNeuHZKTk59ahf34448bfW61MefPn0dAQAB+/vlniMViKBQKxMTEoF27dti/fz+OHDmC4uJi9OnTBx999BFGjRqFkSNH4syZM5g7d26TctfW1r7wip+6ujqkpKTg3LlziI+PR3x8PMrLy6Gnp4fa2lqsX7++oeg2lUKhwOTJk3H8+PGnrkcCHv1BYNCgQUhMTMS5c+dw6dIl+Pn5NWns0tJSREREYNu2bUhLS4OZmRnc3d3h4eGBjh07vlROIiJqNVhwiYiIVC0rKwuff/45jh8/3uj7Ghoa+Pvf/449e/Y0eRvviRMn4OfnhytXrkBPTw8PHjzA+++/j2nTpmHKlCkwNzdX5Vd4oZqaGhw8eBBfffUVKioqUF5ejkmTJmHBggWwt7dv0hifffYZtm/f/szrkYBHq7geHh7YtWsX2rdvj+Li4uf+YeDxau2+ffugoaGBGTNmYPbs2bC1tX3p70hERK1OOE9RJiIiUjETExOIRKJnrorW19dj//79TT61uKKiAlevXsWdO3cgkUigoaGBkJAQpKam4p///Kfayy0AaGlpYcqUKSgoKEBxcTH279+PsrIyjB07FhYWFoiMjMTz/oYeEhKCrVu3PrfcAo+uPAoJCWk4GTo2Nvapz5SUlEAul8Pc3BxWVlZISkrC5s2bUVRUhODgYJZbIqK3CAsuERGRiiUmJuL48eONnoT8mEKhwHfffQdfX99nfqa+vh7btm2DkZERVq5cCU9PT5SUlODevXtwdHRsjugvRVdXF8Cjsuvi4oK4uDhcunQJvXr1gpubG6ytrZGUlPTUz8XExGDhwoUvNZdSqYSGhgZ27NgB4NHv79SpU3B1dUXPnj2xcuVK2NjYICUlBYmJifD09HzmidVERPTm4hZlIiIiFbO3t8e5c+eeW3AfE4lEWLVqFZYtW/bE64mJiViwYAHS0tLg4+ODgICAl7q6R2jJycnw8fHBzz//DC8vL6xduxYGBgZITU2FtbU1qqurGz3VWVNTE/X19c9c/dXQ0ICvry8iIiJQUFAAOzs7eHh44G9/+1tD4SYiorcWn8ElIiJSpYSEBIwcORKampoQi8Woqalp0s8FBQVBJpNBoVBg48aNWLFiBaytrbF9+3aYmZk1c+rmoVQqsW/fPvj5+UFXVxdbtmyBh4cHSkpKIBKJnrgGycDAAIMHD4axsTEOHz6M+/fvP/PgKRMTE4wfPx6ffPIJTE1N1f21iIio5WLBJSIiUqXbt28jJSUFeXl5yM3NRU5ODrKzs1FQUNBQ5jQ1NSGRSFBTU9PwDKpIJEJQUBCio6Nx7tw5rF27Fr6+vm/ENTa///47Zs2ahbi4OCiVSnTv3h1WVlawtLSEhYUFBg8ejB49eiAuLg5TpkxBTU3Nc1e/33vvPeTl5b0RvxsiIlIpFlwiIiJ1UCgUuHnzJvLy8hr+yc3NRWZmJq5du4aHDx8CAN555x3ExsbCyspK4MSqlZycjODgYOzcuRMLFy6EXC5/4gTpHTt2YOHChVAqlY1uXf6r8+fPw8bGpjkjExFR68OCS0REJLTs7GyMGjUK+vr6WLJkCebMmQMdHR2hYzWLqKgozJw5E5MmTUJYWBgAYNmyZdiwYUOTx5BIJJg1axZ27drVXDGJiKh1YsElIiISUlFREWxsbPDuu+/i+PHjMDAwEDpSs4uPj4ejoyPc3NxQUlLS6NU/jdHQ0IBYLIZCoYCWlhZKSkqgr6/fzGmJiKgVYcElIiISyoMHDzBixAjU19fj3LlzMDQ0FDqS2uzZswfz5s2DUqlsuM5HX18fYrEY+vr60NTURLt27SCRSKCvrw9dXV3o6upCR0cHOjo60NXVxZw5cyCVSgX+JkRE1IKEawqdgIiI6G3l4+ODoqIiJCcnv1XlFgDc3Nzw8OFDfPbZZ4iJicHIkSOFjkRERG8AruASEREJ4PDhw5g6dSqioqLg7OwsdBzBTJ06FYmJibhy5cpbsT2biIiaFbcoExERqVt1dTWkUins7Oywd+9eoeMI6u7duzAxMYGHhwfWrVsndBwiImrdwsUv/gwRERGp0jfffIOSkhKsXbtW6CiCMzQ0xPLlyyGXy3Hjxg2h4xARUSvHgktERKRG9fX1CAoKwqefforu3bsLHadF8Pb2hqGhIbZv3y50FCIiauVYcImIiNTo5MmTuHXrFjw9PYWO0mJoaWlh7ty52LNnD+rq6oSOQ0RErRgLLhERkRqFh4fDxsYGffv2FTpKi/KPf/wDt2/fRnx8vNBRiIioFWPBJSIiUqOEhAQ4ODgIHaPFMTY2hrGxMc6fPy90FCIiasVYcImIiNSkpKQE165dw4gRIwTLsH//fujp6UEkEmHDhg2or68HAISFhUFbW1vQU52tra3x66+/CjY/ERG1fiy4REREalJQUAAAMDExESzDzJkz4evrCwBwcnKChoYGAMDOzg4TJ07EnDlzBMvWv39/XL9+XbD5iYio9WPBJSIiUpM//vgDwKOrcYTk4+MDfX19BAUFNby2f/9+zJ8/X8BUj34vd+/eFTQDERG1biy4REREalJZWQkA0NHRETRHx44dsXjxYuzduxdFRUUAgNOnT2P8+PGC5mrbti0ePHggaAYiImrdWHCJiIjUpEOHDgCA0tJSgZMAvr6+0NLSQlBQEJKSkjBs2LCG7cpCuXv3Ljp27ChoBiIiat00hQ5ARET0tni8Nfn3338XfJuyoaEhFi5ciO+++w4lJSVYsWKFoHmAlvF7ISKi1o0ruERERGrSv39/aGlpISUlRegoAIAlS5agpqYGhYWF6NOnj9BxkJycjIEDBwodg4iIWjEWXCIiIjXR1dWFhYUFfvnlF6GjAAC6dOkCBwcHwQ+XAgCFQoGLFy8KeoUSERG1fiy4REREajRmzBhER0dDoVAIHQUPHz7E1atXMW3aNKGj4KeffkJpaSlGjx4tdBQiImrFWHCJiIjUaO7cuSgsLMTp06eFjoJt27Zh8eLF0NXVFToKdu7cieHDh0MqlQodhYiIWjGRUqlUCh2CiIjobWJnZwd9fX0cO3ZM7XNfvHgRnp6eePjwIerr63H16lVoaWmpPcefFRYWwsTEBN98802L2C5NREStVjhXcImIiNRs1apViI2NxalTp9Q+t56eHu7fvw+xWIywsDDByy0ALF++HF27dsWsWbOEjkJERK0cV3CJiIgEMHHiRNy4cQOXL1+Gtra20HEE88svv8DOzg5hYWFwc3MTOg4REbVu4Sy4REREAsjPz8fgwYMxb948bN68Weg4gigrK8PgwYNhamqKmJgYiEQioSMREVHrxi3KREREQjA2Nsb27dshl8tx8OBBoeOoXX19PebMmYOqqirs3buX5ZaIiFRCU+gAREREb6uZM2fi4sWLmDVrFgwNDfHRRx8JHUltvL29cfLkScTFxaFTp05CxyEiojcEV3CJiIgEFBQUhClTpsDZ2Rnx8fFCx2l2CoUCPj4+2LVrF8LDw2FjYyN0JCIieoOw4BIREQlILBZj7969GD9+PMaNG4eIiAihIzWb6upqzJw5E9u3b8cPP/yASZMmCR2JiIjeMCy4REREAtPS0sKBAwewYMECzJgxA0uXLkVdXZ3QsVTq+vXr+PDDDxEbG4vY2FhMnz5d6EhERPQGYsElIiJqAcRiMYKCghAcHAy5XI6RI0ciPz9f6FgqERERgcGDB+PBgwf49ddfMXr0aKEjERHRG4oFl4iIqAX55JNPcPnyZVRUVMDc3ByrV69GdXW10LFeSV5eHiZOnIjp06dj+vTpuHjxIqRSqdCxiIjoDcaCS0RE1MKYmZkhKSkJK1euxPr16zFw4ECEhYWhvr5e6GhNcufOHXz55ZcwNzdHYWEhzp49i2+//Ra6urpCRyMiojccCy4REVELJJFI4O/vj8zMTIwYMQJz5syBubk5QkNDUVVVJXS8RhUWFsLPzw/GxsYIDQ3F+vXrkZycjFGjRgkdjYiI3hIipVKpFDoEERERPV9OTg5Wr16NAwcOoF27dpgzZw7mz58PMzMzQXPV1tbi+PHjCAkJQWxsLDp37owvvvgCCxYsQJs2bQTNRkREb51wFlwiIqJW5Pbt29i1axd27NiB69evQyqVYtq0aXB2doaFhQU0NDSaPUN5eTnOnDmDqKgoREdHo7S0FGPGjIGXlxcmT54MiUTS7BmIiIgawYJLRETUGikUCpw/fx5RUVE4dOgQCgsLYWBgAFtbW9jZ2WHIkCEYMGAAunXr9lrz1NXVITs7G+np6bhw4QISEhKQmpoKhUIBa2trTJ06FVOnToWRkZGKvhkREdErY8ElIiJq7ZRKJdLT0xEfH4+EhAScP38eRUVFAICOHTuif//+6Nq1K3r06IEuXbqgXbt20NbWRps2baCtrY3y8nLU1dWhvLwcZWVluHnzJkpKSlBQUICcnBzU1NRAU1MTZmZmGDVqFEaOHImRI0eic+fOAn9zIiKiJ7DgEhERvYnu3r2LtLQ0ZGRkIDc3F7dv38bNmzdx584d3L9/H9XV1Xjw4AFqamrQtm1bSCQS6Ovro127dujevTveffdd9OzZE6amphgwYACkUim0tbWF/lpERETPw4JLREREREREb4RwXhNEREREREREbwQWXCIiIiIiInojsOASERERERHRG0ETQKTQIYiIiIiIiIhe08X/B4KC3GYdhE6FAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from IPython.display import Image, display\n", "display(Image(filename=\"causal_model.png\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We get a causal graph. Now identification and estimation is done. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_identifier:Common causes of treatment and outcome:['W1', 'W4', 'W0', 'W2', 'Unobserved Confounders', 'W3']\n", "WARNING:dowhy.causal_identifier:If this is observed data (not from a randomized experiment), there might always be missing confounders. Causal effect cannot be identified perfectly.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "WARN: Do you want to continue by ignoring any unobserved confounders? (use proceed_when_unidentifiable=True to disable this prompt) [y/n] y\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_identifier:Instrumental variables for treatment and outcome:['Z0', 'Z1']\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Estimand type: nonparametric-ate\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "─────(Expectation(y|W1,W4,W0,W2,W3))\n", "d[v₀] \n", "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W1,W4,W0,W2,W3,U) = P(y|v0,W1,W4,W0,W2,W3)\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" ] } ], "source": [ "identified_estimand = model.identify_effect()\n", "print(identified_estimand)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 1: Regression\n", "\n", "Use linear regression." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_estimator:INFO: Using Linear Regression Estimator\n", "INFO:dowhy.causal_estimator:b: y~v0+W1+W4+W0+W2+W3\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "*** Causal Estimate ***\n", "\n", "## Target estimand\n", "Estimand type: nonparametric-ate\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "─────(Expectation(y|W1,W4,W0,W2,W3))\n", "d[v₀] \n", "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W1,W4,W0,W2,W3,U) = P(y|v0,W1,W4,W0,W2,W3)\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", "## Realized estimand\n", "b: y~v0+W1+W4+W0+W2+W3\n", "## Estimate\n", "Value: 10.000000000000021\n", "\n", "## Statistical Significance\n", "p-value: <0.001\n", "\n", "Causal Estimate is 10.000000000000021\n" ] } ], "source": [ "causal_estimate_reg = model.estimate_effect(identified_estimand,\n", " method_name=\"backdoor.linear_regression\",\n", " test_significance=True)\n", "print(causal_estimate_reg)\n", "print(\"Causal Estimate is \" + str(causal_estimate_reg.value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 2: Stratification\n", "\n", "We will be using propensity scores to stratify units in the data." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_estimator:INFO: Using Propensity Score Stratification Estimator\n", "INFO:dowhy.causal_estimator:b: y~v0+W1+W4+W0+W2+W3\n", "/home/amshar/python-environments/vpy36/lib/python3.6/site-packages/sklearn/utils/validation.py:744: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n", " y = column_or_1d(y, warn=True)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "*** Causal Estimate ***\n", "\n", "## Target estimand\n", "Estimand type: nonparametric-ate\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "─────(Expectation(y|W1,W4,W0,W2,W3))\n", "d[v₀] \n", "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W1,W4,W0,W2,W3,U) = P(y|v0,W1,W4,W0,W2,W3)\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", "## Realized estimand\n", "b: y~v0+W1+W4+W0+W2+W3\n", "## Estimate\n", "Value: 10.173499320316472\n", "\n", "Causal Estimate is 10.173499320316472\n" ] } ], "source": [ "causal_estimate_strat = model.estimate_effect(identified_estimand,\n", " method_name=\"backdoor.propensity_score_stratification\",\n", " target_units=\"att\")\n", "print(causal_estimate_strat)\n", "print(\"Causal Estimate is \" + str(causal_estimate_strat.value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 3: Matching\n", "\n", "We will be using propensity scores to match units in the data." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_estimator:INFO: Using Propensity Score Matching Estimator\n", "INFO:dowhy.causal_estimator:b: y~v0+W1+W4+W0+W2+W3\n", "/home/amshar/python-environments/vpy36/lib/python3.6/site-packages/sklearn/utils/validation.py:744: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n", " y = column_or_1d(y, warn=True)\n", "/mnt/c/Users/amshar/code/dowhy/dowhy/causal_estimators/propensity_score_matching_estimator.py:62: FutureWarning: `item` has been deprecated and will be removed in a future version\n", " control_outcome = control.iloc[indices[i]][self._outcome_name].item()\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "*** Causal Estimate ***\n", "\n", "## Target estimand\n", "Estimand type: nonparametric-ate\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "─────(Expectation(y|W1,W4,W0,W2,W3))\n", "d[v₀] \n", "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W1,W4,W0,W2,W3,U) = P(y|v0,W1,W4,W0,W2,W3)\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", "## Realized estimand\n", "b: y~v0+W1+W4+W0+W2+W3\n", "## Estimate\n", "Value: 10.036816324727294\n", "\n", "Causal Estimate is 10.036816324727294\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/mnt/c/Users/amshar/code/dowhy/dowhy/causal_estimators/propensity_score_matching_estimator.py:77: FutureWarning: `item` has been deprecated and will be removed in a future version\n", " treated_outcome = treated.iloc[indices[i]][self._outcome_name].item()\n" ] } ], "source": [ "causal_estimate_match = model.estimate_effect(identified_estimand,\n", " method_name=\"backdoor.propensity_score_matching\",\n", " target_units=\"atc\")\n", "print(causal_estimate_match)\n", "print(\"Causal Estimate is \" + str(causal_estimate_match.value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 4: Weighting\n", "\n", "We will be using (inverse) propensity scores to assign weights to units in the data. DoWhy supports a few different weighting schemes:\n", "1. Vanilla Inverse Propensity Score weighting (IPS) (weighting_scheme=\"ips_weight\")\n", "2. Self-normalized IPS weighting (also known as the Hajek estimator) (weighting_scheme=\"ips_normalized_weight\")\n", "3. Stabilized IPS weighting (weighting_scheme = \"ips_stabilized_weight\")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_estimator:INFO: Using Propensity Score Weighting Estimator\n", "INFO:dowhy.causal_estimator:b: y~v0+W1+W4+W0+W2+W3\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "*** Causal Estimate ***\n", "\n", "## Target estimand\n", "Estimand type: nonparametric-ate\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "─────(Expectation(y|W1,W4,W0,W2,W3))\n", "d[v₀] \n", "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W1,W4,W0,W2,W3,U) = P(y|v0,W1,W4,W0,W2,W3)\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", "## Realized estimand\n", "b: y~v0+W1+W4+W0+W2+W3\n", "## Estimate\n", "Value: 10.722320441623154\n", "\n", "Causal Estimate is 10.722320441623154\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/amshar/python-environments/vpy36/lib/python3.6/site-packages/sklearn/utils/validation.py:744: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n", " y = column_or_1d(y, warn=True)\n" ] } ], "source": [ "causal_estimate_ipw = model.estimate_effect(identified_estimand,\n", " method_name=\"backdoor.propensity_score_weighting\",\n", " target_units = \"ate\",\n", " method_params={\"weighting_scheme\":\"ips_weight\"})\n", "print(causal_estimate_ipw)\n", "print(\"Causal Estimate is \" + str(causal_estimate_ipw.value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 5: Instrumental Variable\n", "\n", "We will be using the Wald estimator for the provided instrumental variable." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_estimator:INFO: Using Instrumental Variable Estimator\n", "INFO:dowhy.causal_estimator:Realized estimand: Wald Estimator\n", "Realized estimand type: nonparametric-ate\n", "Estimand expression:\n", " -1\n", "Expectation(Derivative(y, Z0))⋅Expectation(Derivative(v0, Z0)) \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", "Estimand assumption 3, treatment_effect_homogeneity: Each unit's treatment ['v0'] is affected in the same way by common causes of ['v0'] and y\n", "Estimand assumption 4, outcome_effect_homogeneity: Each unit's outcome y is affected in the same way by common causes of ['v0'] and y\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "*** Causal Estimate ***\n", "\n", "## Target estimand\n", "Estimand type: nonparametric-ate\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "─────(Expectation(y|W1,W4,W0,W2,W3))\n", "d[v₀] \n", "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W1,W4,W0,W2,W3,U) = P(y|v0,W1,W4,W0,W2,W3)\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", "## Realized estimand\n", "Realized estimand: Wald Estimator\n", "Realized estimand type: nonparametric-ate\n", "Estimand expression:\n", " -1\n", "Expectation(Derivative(y, Z0))⋅Expectation(Derivative(v0, Z0)) \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", "Estimand assumption 3, treatment_effect_homogeneity: Each unit's treatment ['v0'] is affected in the same way by common causes of ['v0'] and y\n", "Estimand assumption 4, outcome_effect_homogeneity: Each unit's outcome y is affected in the same way by common causes of ['v0'] and y\n", "\n", "## Estimate\n", "Value: 6.7777521025251435\n", "\n", "Causal Estimate is 6.7777521025251435\n" ] } ], "source": [ "causal_estimate_iv = model.estimate_effect(identified_estimand,\n", " method_name=\"iv.instrumental_variable\", method_params = {'iv_instrument_name': 'Z0'})\n", "print(causal_estimate_iv)\n", "print(\"Causal Estimate is \" + str(causal_estimate_iv.value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 6: Regression Discontinuity\n", "\n", "We will be internally converting this to an equivalent instrumental variables problem." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:dowhy.causal_estimator:Using Regression Discontinuity Estimator\n", "INFO:dowhy.causal_estimator:\n", "INFO:dowhy.causal_estimator:INFO: Using Instrumental Variable Estimator\n", "INFO:dowhy.causal_estimator:Realized estimand: Wald Estimator\n", "Realized estimand type: nonparametric-ate\n", "Estimand expression:\n", " \n", "Expectation(Derivative(y, local_rd_variable))⋅Expectation(Derivative(v0, local\n", "\n", " -1\n", "_rd_variable)) \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", "Estimand assumption 3, treatment_effect_homogeneity: Each unit's treatment ['local_treatment'] is affected in the same way by common causes of ['local_treatment'] and local_outcome\n", "Estimand assumption 4, outcome_effect_homogeneity: Each unit's outcome local_outcome is affected in the same way by common causes of ['local_treatment'] and local_outcome\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " local_rd_variable local_treatment local_outcome\n", "1 0.491077 True 8.622930\n", "17 0.526249 True 1.107572\n", "18 0.557455 True 4.576484\n", "24 0.416279 True -5.730869\n", "25 0.554845 True 1.196812\n", "... ... ... ...\n", "9974 0.531890 True 1.377569\n", "9982 0.575699 True 14.511282\n", "9995 0.577368 True 9.653687\n", "9998 0.489953 True 1.844830\n", "9999 0.484942 True 15.116874\n", "\n", "[1924 rows x 3 columns]\n", "*** Causal Estimate ***\n", "\n", "## Target estimand\n", "Estimand type: nonparametric-ate\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "─────(Expectation(y|W1,W4,W0,W2,W3))\n", "d[v₀] \n", "Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W1,W4,W0,W2,W3,U) = P(y|v0,W1,W4,W0,W2,W3)\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", "## Realized estimand\n", "Realized estimand: Wald Estimator\n", "Realized estimand type: nonparametric-ate\n", "Estimand expression:\n", " \n", "Expectation(Derivative(y, local_rd_variable))⋅Expectation(Derivative(v0, local\n", "\n", " -1\n", "_rd_variable)) \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", "Estimand assumption 3, treatment_effect_homogeneity: Each unit's treatment ['local_treatment'] is affected in the same way by common causes of ['local_treatment'] and local_outcome\n", "Estimand assumption 4, outcome_effect_homogeneity: Each unit's outcome local_outcome is affected in the same way by common causes of ['local_treatment'] and local_outcome\n", "\n", "## Estimate\n", "Value: 22.999383345332262\n", "\n", "Causal Estimate is 22.999383345332262\n" ] } ], "source": [ "causal_estimate_regdist = model.estimate_effect(identified_estimand,\n", " method_name=\"iv.regression_discontinuity\", \n", " method_params={'rd_variable_name':'Z1',\n", " 'rd_threshold_value':0.5,\n", " 'rd_bandwidth': 0.1})\n", "print(causal_estimate_regdist)\n", "print(\"Causal Estimate is \" + str(causal_estimate_regdist.value))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.9" } }, "nbformat": 4, "nbformat_minor": 2 }