{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Estimating the effect of a Member Rewards program\n", "An example on how DoWhy can be used to estimate the effect of a subscription or a rewards program for customers. \n", "\n", "Suppose that a website has a membership rewards program where customers receive additional benefits if they sign up. How do we know if the program is effective? Here the relevant causal question is:\n", "> What is the impact of offering the membership rewards program on total sales?\n", "\n", "And the equivalent counterfactual question is, \n", "> If the current members had not signed up for the program, how much less would they have spent on the website?\n", "\n", "In formal language, we are interested in the Average Treatment Effect on the Treated (ATT). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## I. Formulating the causal model\n", "Suppose that the rewards program was introduced in January 2019. The outcome variable is the total spends at the end of the year. \n", "We have data on all monthly transactions of every user and on the time of signup for those who chose to signup for the rewards program. Here's what the data looks like." ] }, { "cell_type": "code", "execution_count": 1, "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", "
user_idsignup_monthmonthspendtreatment
0011449True
1012583True
2013519True
3014581True
4015549True
..................
119995999908410False
119996999909409False
1199979999010398False
1199989999011401False
1199999999012394False
\n", "

120000 rows × 5 columns

\n", "
" ], "text/plain": [ " user_id signup_month month spend treatment\n", "0 0 1 1 449 True\n", "1 0 1 2 583 True\n", "2 0 1 3 519 True\n", "3 0 1 4 581 True\n", "4 0 1 5 549 True\n", "... ... ... ... ... ...\n", "119995 9999 0 8 410 False\n", "119996 9999 0 9 409 False\n", "119997 9999 0 10 398 False\n", "119998 9999 0 11 401 False\n", "119999 9999 0 12 394 False\n", "\n", "[120000 rows x 5 columns]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Creating some simulated data for our example\n", "import pandas as pd\n", "import numpy as np\n", "num_users = 10000\n", "num_months = 12\n", "\n", "signup_months = np.random.choice(np.arange(1, num_months), num_users) * np.random.randint(0,2, size=num_users) # signup_months == 0 means customer did not sign up\n", "df = pd.DataFrame({\n", " 'user_id': np.repeat(np.arange(num_users), num_months),\n", " 'signup_month': np.repeat(signup_months, num_months), # signup month == 0 means customer did not sign up\n", " 'month': np.tile(np.arange(1, num_months+1), num_users), # months are from 1 to 12\n", " 'spend': np.random.poisson(500, num_users*num_months) #np.random.beta(a=2, b=5, size=num_users * num_months)*1000 # centered at 500\n", "})\n", "# A customer is in the treatment group if and only if they signed up\n", "df[\"treatment\"] = df[\"signup_month\"]>0\n", "# Simulating an effect of month (monotonically decreasing--customers buy less later in the year)\n", "df[\"spend\"] = df[\"spend\"] - df[\"month\"]*10\n", "# Simulating a simple treatment effect of 100\n", "after_signup = (df[\"signup_month\"] < df[\"month\"]) & (df[\"treatment\"])\n", "df.loc[after_signup,\"spend\"] = df[after_signup][\"spend\"] + 100\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The importance of time\n", "Time plays a crucial role in modeling this problem. \n", "\n", "Rewards signup can affect the future transactions, but not those that happened before it. In fact, the transactions prior to the rewards signup can be assumed to cause the rewards signup decision. Therefore we split up the variables for each user:\n", "\n", "1. Activity prior to the treatment (assumed a cause of the treatment)\n", "2. Activity after the treatment (is the outcome of applying treatment)\n", "\n", "Of course, many important variables that affect signup and total spend are missing (e.g., the type of products bought, length of a user's account, geography, etc.). This is a critical assumption in the analysis, one that needs to be tested later using refutation tests. \n", "\n", "Below is the causal graph for a user who signed up in month `i=3`. The analysis will be similar for any `i`. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import dowhy\n", "\n", "# Setting the signup month (for ease of analysis)\n", "i = 3" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " user_id signup_month treatment pre_spends post_spends\n", "0 2 0 False 479.0 390.888889\n", "1 4 3 True 487.0 522.444444\n", "2 6 0 False 482.5 422.888889\n", "3 8 0 False 473.0 418.444444\n", "4 10 0 False 489.0 424.333333\n", "... ... ... ... ... ...\n", "5326 9987 0 False 480.0 414.000000\n", "5327 9990 0 False 495.5 421.777778\n", "5328 9992 0 False 473.0 405.666667\n", "5329 9996 0 False 490.0 415.555556\n", "5330 9999 0 False 482.0 420.111111\n", "\n", "[5331 rows x 5 columns]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwAAAAFpCAIAAACkh9hiAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydd1xT1///TwYzjLC3QABBlgMVFRwg1Za6R22t41O1trZa7bSbDj/VLlcdtdtVrXVVHFgFJ6IsKbL3HmGFJGQASX5/nF/vJ98QYgghNwnv5x953Nzc3Pu6yb3nvs77vM85FJlMhgAAAAAAAIYTVLIFAAAAAAAA6BowQAAAAAAADDvAAAEAAAAAMOwAAwQAAAAAwLADDBAAAAAAAMMOOtkCjBapVNrZ2SkUCkUiEbGys7NTKpXKb6DwLWtrazr9f3+KnZ0dsWxqaspgMBQ2AAB1EAqFbDZbIpFwuVzilU6n48sJvzo7O5ubm5OtFAB0R1dXV3d3N74diJW9vb08Hk9+M4ViHEOhUJhMpvwaXETLr8F3FpPJpFAo2tYOaAF4lD6e3t7etra29vb29vZ2hQWBQCAUCjs7O0UiUVdXF4/HE4lEPB4P31dDpAc/sRgMhrm5ua2trYWFhbm5uZ2dnbm5uYWFhb29vYODA/FKLJiamg6RHkB/6O3tLS0tzcvLy8vLq6qqavgXDoejztft7e3d3Nw8PDzc3NxYLFZISEhYWJifnx+NRhtq5QAwUHp7ezs7OzlydHZ24jW48snj8Xp7ezkcjlQq5XA42Nl0d3d3dXUp9TRDB3ZLuOjGPgmX2wwGw9TU1MbGxsTExNbW1s7OztbWlslkEq94QcFXAdqCAuMA9fb2stnsurq6xsbGurq6pqam2tra5ubm1tbWtra2trY2Lpc70H3a2tqamZlZWVnJB2zwtS6/DZWqpAmyo6NDXhtRF8ERI4FAIBaL5bdREysrK+yEnJycXF1d3d3d3d3d8aPO09PT1dXVxMRkoPsE9IGKiorbt2/funXr4cOHhYWFfZ03nU53cXFxcXExMTEhQj74IyIg1NPT09zc3NzcLF8Vxpibm48aNSoiImLatGnTp08fMWKELs4KGMb09PS0tLS0tLQ0NTWx2ezW1la80NLS0tHRgf0Nh8Pp6urS+BBEbIZKpcoHcmg0mo2NjfyWuGKpdCcymUyhaiEWiwUCAfEW2yylcaYBgaViM2RnZ+fk5OTk5OTo6Ojq6urs7Ozk5OTs7Ozq6mplZaXZ/octw8gAdXZ2lpeXV1RUlJeXV1dX19fX48qx0kKfgE6nK4RS5JdxGIbJZGI7z2Qyzc3NLS0tdXA6IpFIKBRyOByxWNzV1cXlckUiEZ/P7+jo6C9epSIoRaFQXFxccO3f3d3d29ubxWL5+fmxWCwHBwcdnA4wIDgczsWLF5OSkm7dulVXV0est7KyGjVqVFhYWEhIyMiRI93d3d3c3FxcXJRa7b5IJJLm5mZ8XxQXF+fn5+fl5RUWFsqX6T4+PtOmTYuPj4+PjyeMFAAMCLFYXFdXV1dXV1NTU1tbi8thbHHYbHZbW9tj94DLW9t/wc6AWMZVUDs7O2xoTExMrKyssJWxtLQ0MzPTwTkqhcPhSCSSzs7Onp4ePp+PDRNRrSXCV/KhLGLhsTu3sLDA9VvskNzd3T09Pb28vEaMGOHp6Wlvb6+DEzQsjNMA1dXVlZeXE3YHv/Z3Uzk4OBBREC8vL1dXV7zs7Oxsb2+vUBswaPh8fltbW2tra2NjI37I1dfXE3Gv5uZmpd9iMpmEGSJevb29oVVb97DZ7PPnz587dy4lJYWws0FBQdOmTZs2bdrkyZN9fX21/r9IpdLy8vK0tLRbt27dvn27rKwMrzc3N4+Li1u0aNG8efPAJQN9kclk9fX12OIQXqeurq62trapqam/b1GpVPz8xlENvIwLZEdHRwcHB8Lf6PJc9AHCIWGbiF+bmppwtKyxsbGlpUUoFPb3dUtLS29vb09PT09PzxEjRmBX5Onp6evr21+Iy+gxBgPEZrMfPXqUn5+PK6z5+fl9zTKNRvP09CSe3z4+Pp6enrgZCBI/Md3d3U1NTXV1dQ0NDZWVlYR3rK6u7u3tVdiYwWAEBweHhYXh15CQEA8PD1JkDwekUun169cPHTp04cIF/F9YW1vHx8cvWLAgJibGxcVFl2IaGhpSUlLOnz9/5coVHBkyNTVdtGjRSy+9NH36dLDFw5a6urqysrLS0lL51/7ybFxcXPDTFz+SPTw8iNYcJycnNQOWQF/4fD5uy25paWloaMB2s6amBsfbxGJx369QKBRPT8+AgAB/f3/51+FgMQ3PAPX09Pzzzz9ZWVnY6zx69Ki1tVV+AxMTk4CAgICAAPmghY+PD2QBa0Zvb29NTQ0RSMMUFRUp3Et2dnYhISE4bXbs2LFjx44dtrUKLcLhcH744YcffvihvLwcIWRra7t48eJFixbFxcWRXjwJhcKkpKRz586dO3eOz+cjhIKCgl5++eU1a9ZA05hxw+fzcW1T3u7It5NizM3N/f39vb29vby85NtivLy8SL96hyc4w1UhGldWVsZmsxW2pFKpXl5ehBkaNWpUaGio8eX/GYYBamhoSE1NvXv3blZWVnZ2tkKUz83NLSQkJDg4OCIiAi/Ao3eowa4oPz+/oKAAv+bl5clbIjqdPnLkyIiIiIiIiOjo6DFjxkBPogHB4/EOHDjw5Zdf4oT3iIiI9evXL1++XA/zHHk83okTJw4ePJiTk4MQcnBw2Lhx4+uvv25ra0u2NEALSCSS6upq4mbPysoqKioihvPAmJqaenp6slgsFosVHBwcEhKCq50QyzEIRCJReXl5QUFBhRyVlZUK9sDW1tbf35941I4ZM8bR0ZEszVpBTw0Ql8u9fft2WlpaRkZGRkaGfKY9g8EYO3bshAkTwsLCcBOMbpKOAdWIxeKCgoKCgoJHjx5lZGRkZmbK956ztraOiIiYMGHC5MmTp02bBikjKuByud98883evXs7OztpNNqzzz67ZcuW8ePHk63r8dy7d2/Xrl1nz56VSqUODg5vvPHGli1b4PY0OPh8flZWVmZmZl5eHs4uUGjJcnZ2Dg8PDw0NDQwMxEECLy8v8DpGBofDweG90tJS3NhSWlqqkA7h4+MTFhYWGho6ZsyYyMhIb29vstRqhh4ZIKFQmJWVlZqaev369du3bxM5njQaLTAwMOJfJk6cCI1ZBoF83C4rK0u+DGWxWHH/Ij/YI5CYmPjqq6/W1tZSqdTFixd/9tlnQUFBZIsaGPn5+V9++eXvv/8ukUg8PDy++OKLVatWkS0KeAwVFRXErZqRkSHfY9TU1NTf35+Ir48fP97NzY1EqQBZ9PT0lJSUEIHAgoIChSgRk8kcP358VFRUREREVFSU/vc7I9kA9fT0pKWl3bhxIyUl5f79+8RdZ2dnhzu2TJw4cdy4cVCJNHS6u7tzcnLS09Pv3r1748YNosmZTqdHRETExsbGxMRER0cP57bLoqKiV199NSUlBSG0cOHCHTt2jBw5kmxRmpOXl/f2228nJSUhhOLj47/77jsWi0W2KOB/NDc3379/Pz09/f79+xkZGfJjHzs4OERGRuIoe3h4uJ+fH0R3AKW0t7fn5ubm5eVlZWWlp6fLt41SqdRRo0ZNnDhx0qRJkZGRYWFhengVkWOABAJBcnLyn3/+mZiYSDRvWVpaTpkyJS4uLioqKjIyEsblM2IqKiquX79+/fr15OTk9vZ2vNLCwmLmzJlz586dN2+eq6sruQp1iUwm+/HHH7ds2SIUClks1t69e59++mmyRWmHxMTETZs2VVdXm5ub79ixY/PmzWQrGtbw+fz79+/jWy87O5so/HHGXnR0NK67BwcHQ28+QAN4PB7uopSamnr79m35oVWsrKwmTZqEo/7jxo3TkwtMpwaorq4uMTHx/PnzN2/exMEeOp0eHR09c+bMmJiYiRMngukZbkgkkpycnJSUlOTk5Fu3buFmMhqNFhUVNX/+/Pnz5/v5+ZGtcWhpbW1ds2ZNYmIinU5/77333nvvPSMLg3V1dX3yySc7d+6USqXLli07dOgQJEfrEoFAkJqaim+x7OxsYtBXDw+PadOmRUZGTpw4cezYsTAaCKB1KisrcZQxLS0tMzOTuPbc3NxmzpwZGxsbGxtLbtqQLgxQc3PzsWPHTp48mZWVhQ9nZWX15JNPzp8/Pz4+Xv+bCQHdwOfzk5KS/vrrr8uXLxNhodDQ0GeeeWb16tXG1wMTIXT//v0lS5bU19d7e3sfP348KiqKbEVDRXJy8sqVKxsbG1ks1tmzZ0ePHk22IiOnvLz83LlzFy9evH//PtE908HBYcaMGfjZExgYSK5CYFjB5XJv3bqVnJyckpKSl5dHGA8/P79Zs2YtXLhwxowZuo+ADKEBkkgkN27c+OGHH86fP9/T04MQcnJyevLJJ+fOnRsfHw+zuwH9IZFI0tLSLl68eO7cuZKSEoQQlUqdPHnyqlWrnn/+eaO5chITE5999lmBQLBgwYKff/7Z6GsCLS0tL7zwwqVLl6ysrE6fPj179myyFRkh+fn5Fy9eTExMTE1NxWsYDMbkyZNx08PYsWP1MA8DGG60tLTcvHnz7t27qampWVlZeCWTyXziiSfmzJmzcOFC3Q0kJhsCcnJyNm/e7OTkhA9ha2v74osvpqamSiSSoTgcYMRkZWVt3ryZGG3C1tZ2/fr19+7dI1vXYPnll1/wLLlbt24lW4vukEqlCQkJCCFTU9Pjx4+TLcdIkEgkd+7ceeONN3x9fYmC3cfH5/XXX79161ZPTw/ZAgGgX6qqqr777rvY2Fhi1nALC4sFCxYcPny4ra1tqI+uTQMkkUjOnz8/ZcoUfBpUKnXmzJlHjx4VCARaPAowDBGLxWfPnp07dy5xk4wePfrYsWMGWrh/++23CCE6nf7jjz+SrYUEdu7cSaVSqVTqTz/9RLYWw6a2tnbHjh3yWRTBwcFbt269c+eOVColWx0ADID29vZTp06tXLmSGOuVRqPFxcWdOnVq6Mp57RggsVj8888/jxo1Cuv29fX97LPPqqurtbJzACBoamr65ptvgoODiWrunj17+Hw+2boGwOHDhykUipmZ2YULF8jWQhonTpyg0+k0Gu3cuXNkazE8hELh8ePHY2Njia40kyZN+uqrr0pLS8mWBgCDhc/nnz59+vnnnydmIvfw8Pjwww/Ly8u1fqzBGiAul/v1118TE2GOGzfu5MmTvb29WhEHAEqRSqUXL16cNm0avuocHBwSEhJaWlrI1vV4Ll68SKfTqVTqH3/8QbYWkjl69CiFQjE1Nf3777/J1mIw/PPPP5s2bSLGDnVxcXn77beLiorI1gUA2ofP5//6669E1xAKhTJz5szff/9dJBJp6xCaG6Du7u5Dhw4RM1FHRUVduHAB4q6ALsnKylq5ciWeZYzBYGzdupXL5ZItql/y8/PxkJ4HDx4kW4tesGPHDoQQk8nE48kCKrhz586cOXNwyIdKpeKmge7ubrJ1AcCQU1xcvHXrVsJsODs7JyQkaCVDSEMDlJiY6O/vj03ZkiVLcP92ACCFoqKiF154Adsgd3f3I0eO6KERFwqFuO/3u+++S7YWPeLll19GCE2ZMsVA07mGGqlUevr0aWLUAE9Pz23bttXV1ZGtCwB0TXd39+nTp5944glcDbC2tt66dWtra+tg9jlgA9TY2Lhy5Uqi4fnOnTuDOTwAaIvCwsKlS5fiK3PatGmFhYVkK/o/vPbaawihCRMmQK1dHpFIhJ/uH3/8Mdla9I6///57woQJ+JIOCQn57bffxGIx2aIAgGRycnKef/553CHGxsbm008/1TjwPzADdPLkSSaTiRBydXU9efKkZocEDJE333xTfvSER48eka1IOUlJSXjOKQsLi/379+tJKOjGjRsUCsXGxmYo8vgMnby8PAsLCxqNlpmZSbYWfaG0tJSYDiU4OPjMmTN6ciUDwwqhUChf7M+ePZtsRf+jsrJy9erVOPDv5uZ29OhRDe4RdQ1QV1cXMaXziy++2NHRMdAjAUbAjRs39NwAyWQygUCwdetWfGPMmTNHB4NJqEYqlUZERCCEfvjhB90f/c8//0QDgZSY7jfffIMQiomJ0f2h9Q2xWPzRRx+ZmZnhBq/ffvtNz/uUtLe3R0REuLu7Z2RkkK0FGCqwHdcrA4TJz8+fP38+LrumTp060MC/Wgaovr4eR2KdnJzOnz+vkU7AGDAIA4RJS0vD84j5+fkVFBSQqOTo0aMIoVGjRpGS5oINUGJiYn8bCIVCYgCLgICArq4uXcrDiMVinFN48eJF3R9dfygsLMRe2cTE5LXXXtPnjH6C06dP44tny5YtZGsBhgq9NUCY5ORkPDaKhYXF7t271Q8FPd4AFRYWuru74/SF+vr6wekEDBsDMkAymay9vf2JJ55ACDGZzLS0NFI09PT0+Pj4IIQuX75MioDHGqDXX38d/6c0Go3EIbZPnDiBEAoNDR22bT3Hjh3D8+COHz+eXMs+INra2saOHevm5vbgwQOytQADprGxEd/+R48eVbGZnhsgmUwmFovfe+89HPhfsGABj8dT51uPmRemoKAgJiamoaFh6dKlt27dwk4IAAwCOzu7y5cvv/LKKxwOZ/bs2WlpabrXcOXKlaqqqujo6Keeekr3R38sd+7c2bNnD15+++23J0+eTJaSZcuWhYeH5+Xl3blzhywNZCGTyT744IOVK1eKxeIPPvjg3r17RExO/7G3t8/Ozm5oaJg4cSLZWoDhi6mp6RdffHHz5k1vb+/z589HRUVVV1c//msqzFFTU5OXlxdC6Nlnn4VOqoDM0CJABB988AFCyM7OrqSkRMeHXrBgAULoyJEjOj6uOvB4PJwwjhAKCwsjvYcRtmKrV68mV4aOkUqlr7zyCkKIwWCcPXuWbDnA8MJoIkAEra2t06dPRwh5eXk9ttNJvxGgnp6exYsX19bWLlq06Pjx48QcTABgcGzbtm3Dhg0dHR3z58/n8Xg6Oy6bzb506ZKNjc3ixYt1dlD1eeuttyoqKhBCJiYmR44cMTU1JVfPihUrzM3NT506xeFwyFWiS956660DBw44ODjcuXNn4cKFZMsBAMPGwcHh6tWr2L3ExcU1NDSo2ro/Z7R9+3aEUHh4uJptaVpHocFi69attbW1mzZtCgwMZDAYDAZj3Lhxu3btUhhVRWm3vezs7MWLF7u5ueEBlFxcXOS/UlVVtXnz5lGjRjEYDAsLCz8/v//85z+D6ZGbnJy8ePFiDw8PExMTS0vLkSNHrlix4vjx4/JDNml2dgPS3Pen6O3t3bNnz7hx46ysrCwtLUePHr1nz57+Ui4SExNnzZplb29vZmbm7e29Zs2akpISFREgdc6aRLq7u2NiYhBCmzZt0tlBf/31V4TQ2rVrdXZE9UlKSiKujc8//5xsOf+fJUuWIITOnDlDthAdgRPkmUymvo0lq87trE4faaIYMTU19fDweOaZZ1JTU4liBJOQkKB0b48tr65cuSL/FfkYxq5du+Q/amxsxOsHWfCqQKn+3bt3jx492sLCwsHBYfr06adPn8YbZ2ZmLl682MXFxdTU1NfX980331TxnK2trX3rrbfCwsKsra3Nzc1ZLNYLL7ygcMFo8OupSAlwcHBQ0EBEgAb0ECGR7u5u3DssOjpaRfuVcgNUVVVlYWFhYmKSk5MzZArVBU809sQTTzg5OW3fvr2+vl4oFN6/fx9POx8VFdXZ2dn3W8Qfdu3aNZxaSCBvgI4dO2ZmZubo6Pj7779zuVwul3vmzBlPT0+E0EcffaSB2s8++wwh9NRTT2VkZPD5fDabfebMGV9f3/4KCA3ObqCa8U8RFxc3Z86cL774orm5WSgUXr582cnJCSH02muv9f0KHrXPzc3tzz//5HK5HA7nxIkTPj4+n376Kf4NFQzQQM+aFGpqaqytralUqs4eNuvWrUMIHTt2TDeHU5+Ojg5i/r6JEyfqT0fr/fv3I4TeeOMNsoXogoqKCisrKwqFcunSJbK1/B8Gejv31z7Stxj5888/PT09iQFLlbakD7S8UtGIQ4wBQRggAs0eK+pA6J8/f/7u3bvb29ubm5u//fZbnJ+7c+fOc+fOxcbG5uTkCASCBw8ehIaGIoSmT58ukUj67u3kyZMWFhb29va//vprW1sbn8+/du1aSEgIhUJ5//33+zoPLf56j91tUlISnp5C6UOEdLq6usLCwpDKQVaVG6DNmzcjhF5//fUh0zYAiJJ6//798ut5PB7u57x06dK+38J/2JQpUzw8PN57772qqiqhUHj9+nVra2vCACUlJVGpVCqVqtBFqKCgADcH/PjjjwOSWlhYSKFQzM3NBQKBwnpTU1MVBkj9s9NAM/4pqFTqhx9+KL/++++/RwjRaDSFzn34IUSj0R4+fCi/Pisri2gJlS+5NDhrsvjqq68QQs8884xuDodTWfVwoqsVK1bg/9Hc3FyvhszOyclBCE2aNIlsIbrg+eef10O3p8HtrNQA9VeMECWVagOkfnk1GAM00MeKOhD6t23bJr8eD6RnbW0dEBAg764KCgqwkitXrijs6urVq1QqlUKhpKamyq9ns9nYeSgcQqbVX0+d3f74449Kd6sn4IvN0tKyv9ljlBggoVBobW1tYmJSW1s7xPLUAl+pLi4ufeup+B5DCPXt5EyMo6rg/jZs2IANUE9PD67TPP30030P+uyzzyKEnJ2dhUKh+lIPHjyIEHJ1de370aRJk1QYIDXPTjPNxLXLZrPl15eVleH9yw/qLRAI7O3tEUKLFi3qewhiyCn5kkuDsyYLLpfLZDJpNFpDQ8NQH6uzs5NCobi5uQ31gQbK2bNn0b/s2rWLbDn/h97eXhsbG1NTU6OfMKS0tJRKpdrb2+vboLIa3M59DZDqYgQ7v8caIHXKK9ngDNBAHyvqgPXTaDSFIVhxcydC6O2331b4Co7fK6xXXdp/8cUXCCE6na5Qv9Lir6fObsvLy5XuVn/AYcj+BqlSkgSdmprK4/FiY2Pxv6InTJ06FYcQ5Zk9ezZeIK4tBUxNTd944w35NQcOHGhqakIIXbt2rbKyEiE0Z86cvl+MjY1FCLHZbIXmatUwGAyEUFNT0y+//KLwUVpamnzWhQJqnt1gNAcEBOAoKMGIESPwAnEPIISuXLnS3t6OEJo1a1bfncycObPvSo3PWvdYW1svWrRIIpFcu3ZtqI+Fy1xcm9QfWlpa8PyjCKHp06fjWK/+QKPRvL29u7u7W1tbydYytODZLdasWYMnF9IftHI7qy5GcCcd1ahZXg0SzR4r6hAQEIAtIIGbmxteiIyMVNgYuzGFdF2itI+Pj++7f7yyt7f3559/Vnr0ofj1+u4W9xMf5G6HlDfeeINCoeC0wr6fKjFAODcKP031B6VuzNfXFzfKpKenK/1WcHCwra2t0o/u3r2LF3ATrALe3t544cGDB+qLnD17Ni7O1q5dO2PGjKNHj3Z0dKjzRTXPbjCaiduPwMTExMTEBCEkEomIlcSxAgMDlUrqu1LjsyYFfGHfu3dvqA+EH+EODg5DfaAB8fLLL7PZbISQlZXVr7/+irsF6BWOjo7o31/PiLl+/TpCiAip6g9auZ1VFyPqVK3VLK8GiWaPFXXoqx9Pb6L0I3Nzc9TnvIjSPigoqO/+iZVKx80aol9PN3+KdvH29h47dmxtbW1xcXHfT5UYIOxDlT7qSIS4euShUqnW1tYIoebmZqXfwoWpUgjHOnXqVEofiErAgIyts7Pz5cuX8dgqt27dWrVqlZOT0/Tp0/ft29fZ2anii2qe3WA093cIhTXEsZQaR1w7VEDjsyYFfGHroL7S1taG9MwAHTt2jGj/2rlzp77d45hhYoBqamoQQnj8fr1CK7ez6mJEoVeKUtQsrwaJZo8VdVAxqISa400QZZSdnV3fT83MzPDPqLSb9xD9err5U7QONov4jlNAiXSBQIDUu0Z1idL41WPpG97su0OFHD0FDh06NKAjTp48uaSk5Pz584sWLbKwsJBIJLdv3960aROLxUpMTHysGNUMkWalh1AaG+gvYKDZWZMC9nBdXV1DfSBcJcJ1O32goaEBN4cjhJ566qkXX3yRXD39gUsehW69xge+ApXWKEhn8Lez6mJEs8J8KNAfJX1RU5seBnH1DWxnlY4Ap8QA4Ua+lpaWoZY1ILAtU0AqleKzcnV1HegOiWk9tD7qGo1Gmz9//pkzZ1paWk6ePImHn2lvb1++fHltba3Sr6h5dkOnmYA4ltK4t4rHkgZnTQq4AcjZ2XmoD4QzAHAmhD6wdu1a/J/a2dn99NNPZMvpFxw5UxG7NQ5wa8JjRmkjj0HezqqLkSGqfkilUoU1j7XRWn+saBGitFf6G4rFYnx2fZulAAXq6+tRPz+UEgMUEBCAEMrKyhpqWQOiqqqq78qKiore3l6EkAbT0EydOhUvEL0QtQ6DwVi2bFlKSsqOHTsQQnw+Hzf890XNs9OBZiJBr6SkpO+n+EpSjfpnTQqZmZno34t8SMGNX/hxTjqHDh0i0lf37dunz5P66WfulNbBjV+3b98mW8hj0Ox2Vl2MKC3uNIZol+nrZh7rL7X+WNEiRGlfWFjY99OioiKFzQCl9PT03Lt3j06nK02lUmKA4uLiEEJXrlzpa6hJ5O7du93d3Qorr169ihdWrlw50B3GxcXhdu7+grrR0dFUKlXpxdcf27ZtMzU1FYvFCus3bNiAF/p+hFHz7IZCswJPPfUUfvYQR5fn1q1bfVdqfNakcOnSJfTvRT6k4BiGPkRSKysr33rrLby8ZMmS5cuX97elq6ur0g6GugT/YkZvgPCsF8ePHydbiCJauZ1VFyPJycmDlvk/bG1tLS0tEUJ9Q1Opqamqv6v1x4oWIUr7y5cv9/0Ur6TT6WvWrBnMUfQ/g2eQJCUltbe3T5s2TaFTHkbJybNYrPHjx1dUVOBHhZ7Q2dm5d+9e+TV8Ph+Pd/7MM8/07Vj4WOh0+sGDB6lUalJSUt9O0efPn09NTV2yZMlAp2Xu6ek5d+6cwkrs1ul0en/PXTXPbog0y2Nubo5HmLhw4QIeld3V4loAACAASURBVI6guLiYGFdDAc3OWvdkZmampqaOGDFCB9Oee3h42Nra5ufnKw2z6wyZTPbCCy/w+XyEkIuLCx7lRW/p6Oioqqpyd3fXt87hWufJJ5/09PS8du2a0koFuQz+dlZRjOTk5Fy8eFF7YhGVSp02bRrq46sKCwv//vtv1d/V+mNFixCl/ZUrVxR6rba0tOCZgz/66KNBdmVgMpm4yxsuIhBCycnJFAoFTyBt6Eil0o8++ggh1G/Ko9Is2t9//x0hNHr0aH2YBB6PkbB27drg4OBt27bV19eLRKIHDx5ERUUhNabCUL1zPNC4lZXVrl278J5LSkoSEhLMzc3j4uL4fP6ApH7++ecIIWtr62+++aa8vFwkEjU3N585c4bFYlGpVIUhRzU+u4FqVvFT4Ojx9u3bFda//vrrCCE3N7fTp0/zeDwul3v69GkvL69nnnkGXzbyI5hpcNZkgcvub7/9VjeHw4Og3L59WzeHU8rOnTuJm/38+fOqN3ZxcVE66prOwKFNjQfhNSzwKLoK4wKTjga3c38ljEIx0tnZefLkSV9f3/5m1FG9N1n/5VVaWhruWvXBBx+w2WyhUPj3339HREQkJCTgA/U3EOJAHyvq0J9+or96RkaGwkd4YKTFixf33Rsu7R0cHH777bf29vaurq5r167hYVDeeeed/qbCGNCvh4cUeuKJJ1pbW5ubm+fMmWNiYqIwX5AGu9UH8JSmo0ePVjrNiKy/qTB6e3snTJiAEPrvf/87lPLUAl+pW7du5XK5H374YWhoKIPBsLS07G/WOmKMc3m2bt3a3/7r6urefvvtsLAwPLubv7//3Llzz507p8EESQKB4NKlS6+++urEiRO9vLxwT8WRI0euW7euv35bAz27AWnu+1OEhITIZDKlWXXPP/+8/HcvXrw4a9YsOzs7U1NTd3f3RYsW3blzR+kshhqcNSn89ttvCCEWi9XV1aWbI37yyScIoR07dujmcEoZaC2WXAP0/vvvI4R27txJogad0dvbi5OLly1b1l8BrXvUv537SzGWrzYTxQieU/mFF14oKSlJSUnBWxYUFMjvcDDl1f379+fMmWNvb29iYuLl5fXiiy/W1NQoRKzlpwnSrOBVTX/6T5w40Vf/uXPnZMoGSvXz81PYbW1t7ZtvvhkSEsJgMMzNzX19fVevXp2enq7m0R/767HZ7P/85z9ubm50Oh1P2pqUlKSVP4VcUlJS6HS6mZmZiqnN+50N/uHDhyYmJnQ6/e+//x4aeepCXKnkyhgijPvs9Ifs7GxLS0sKhXL16lWdHRT7xWEys5VWwKnB+jY1+tDR2NiIuxq99NJLejil9hBBpDCSOIEUFLzGzd27d62srBBC3333nYrN+k2AGjNmzN69e3t7e5ctW/bw4cP+NgMA/ae8vHzevHkCgeDDDz9UOjb/EDFt2jRvb+/79+8PXa89Y+LevXsFBQUhISHjxo0jW4uOcHV1vXr1qoODw6FDh9asWdM3Idcoqa6uRggxmUxy+5kDxsqVK1eeeuopPp//3nvvbdy4UcWWqjLAX3755ddee62joyMuLk7fesUDgJqUlpbGxMTU1dUtW7aMSD7QDVQqdfXq1QihvtMqAX3BsxqtXbuWbCE6JTw8/Nq1a05OTr/99ltcXBwep8oIuHjxIoVCUTpo0M2bNxFC8fHxRt8FCdAxMpls165dc+fO5fF47777Ls7EV8Fjrr/du3dv2rSpvb19xowZffsFAICek5qaGh0dXVtbu3jx4qNHj+p+1NQXXniBSqX++uuv+jw/mj7Q2Nj4xx9/mJqarlixgmwtumbs2LHZ2dljx469c+dOWFjYhQsXyFakNfoOuZmbm3v+/HkLC4sPP/yQFEmAscJmsxcsWIBnP923bx/OgFbNYwwQhULZs2fPe++9x+fzlyxZ8umnn0okEi2pBYAhRCaT7du3LzY2ls1mr1mz5uTJk3jSPh3j4+Pz3HPPtbe3P7YuMsxJSEjo6upav369wnTTwwRPT89bt24tX74cF+Lr1q3TkyE0B8l///vfjz/+uKysTCwWt7S0HD58eNasWWZmZidPnhzMaB0AII9MJvv9999x5cHb2zslJeXVV19V95vqcPToUTyrUXR0dFVV1aBTlB4PnpReAZx/rmPU/A9wlyg10Z+zM0rYbDYe0I9Go5HepaiystLc3NzU1LSsrIxcJXpLYWEhnU63trZuamoiWwvJnDhxAk9+6eDg8P333+tP77CB0tvbe/PmzS1bthC9yRgMRmho6BtvvFFdXU2WKg0K3qEo/wEt8ujRoxkzZuB/YcWKFRwOR/3vqmuAZDJZXl5eeHg4QsjS0nLHjh36MEQQAPTl1KlTOIowYsSIW7dukS1HJpPJ8EDMc+fOHT49fdRHIpHg3uCff/452Vr0gsbGxvXr1+P8mFGjRp06dQouGwDoS3V19fr16/FAjn5+fpcuXRroHgZggGQymUAg2Lx5M55ifdy4ceSO8AYACmRnZ+MhxXBVoKOjg2xF/5+Ojo4RI0YghPbt20e2Fr0DNw4GBgbqbHwmg+DOnTvjx4/HF/P48ePPnj1ruNEgANAuZWVlL7/8Mh790sbG5vPPPxcKhRrsZ2AGCJOenj527Fh8Z86fP19hMCsA0D3V1dUrV67EleaAgIDr16+TrUiRO3fu0Gg0c3PznJwcsrXoEenp6SYmJmZmZtnZ2WRr0TukUumZM2fw2EgIoZEjR/7www+aFfQAYBxkZGQ888wzOApjaWn59ttvt7a2arw3TQyQTCbr7e09dOgQnl+eSqUuWbKk78CUAKAD8vPz//Of/+CqgJ2d3ddffy0SicgWpRw8MP/IkSPZbDbZWvSC2tpaLy8vNGyGftaM3t7ekydPRkREYBtkZ2e3ceNGsNHAsILD4Rw4cIC4C5hM5rvvvtt3hpOBoqEBwvD5/E8//ZSYZDUmJgbPIT9ITQCgDnfv3p07dy7u2W5lZfXWW2+1tbWRLUoVPT09eBjG8ePHc7lcsuWQTFtbW0hICEJo8eLFUGiow/Xr1+fMmYPrvvgqOnjw4IBSPgHAsJBKpbdv3161apWlpSW+7P38/L755httlZ+DMkAYHo+3c+dOXJNDCI0ePfrYsWMQpwWGCDxVNZ6zECHk5OT02Wef6bn1IeDxeHiWvbi4OLFYTLYc0ujq6poyZQpCaMaMGVBWDIja2trPP/+cmAPc0tJy9erVN27c0GDuQgDQW6qrq7/88svAwEB8nZubmy9fvjwlJUW7lSUtGCBMd3f34cOHcZUOIWRra7ty5cpr165B3Q7QFgUFBVu3biWGz/fx8dm9e7fBZc62tLQEBQXhiOnwrL63tbVFR0cjhMLCwvQnUd2wkEgkd+7cWb9+PVEztre3X7ly5alTp/h8PtnqAEBDysvLd+/eHRUVRQxaO2rUqB07drS0tAzF4bRmgDBSqfTChQvx8fFEnDY4OPirr75qaGjQ7oGA4UNbW9t3331H9IihUCgzZsw4ceKE4Q7EUF1djWs2Y8aMGW63RmVlJfZ/Y8aMgVF/Bk9bW9uePXuio6OJaSWsrKyWLl16/Pjx4WmvAYNDIpGkpaW98847I0eOJMZV8vT03Lhx41DnFlNkag/0NCAaGhqOHj3622+/FRUVIYTodPqsWbMWLVo0d+5cZ2fnoTgiYGRwOJzLly+fPXv24sWLYrEYIeTj47N69erVq1cT8X/Dpa2tbe7cuWlpaT4+PufOnRszZgzZinTB/fv3Fy1a1NjYGBcXd+bMGRsbG7IVGQ9sNvuvv/46f/58cnIyvl9MTU1jYmJmzZoVGxsbHh4OE28BekVLS0tKSkpycvKlS5caGhrwyqCgoAULFixcuHDChAm6mLloSO2VTCbLzMx87bXXHB0d8eGoVGpERERCQgJ0ngeUUlNTc+jQoTlz5uCOXQghc3PzpUuXXrhwwciyHEQi0dKlSxFCZmZmu3fvNu7GYqlUunv3bvyfLlmyBPJ+ho6urq4LFy6sXLnS1taWKOcdHR3nzJmzY8eOzMxMsgUCwxc+n3/t2rWtW7dGRETIO/Lg4OCEhATdX5xDFQFSQCgUJiUl/fXXXxcvXiTmuAkODp4/f/68efMmTJhANJkBwxCpVJqbm5uYmPjXX3/h8WAQQjY2Nk8++eSCBQvmzJljbW1NtsYhQSKRfPHFF3iKvQULFvz8889En0pjgs1mr1q16urVq6amptu3b3/99dd1PyvtMEQsFt++fTs5OTklJSU7O5uYxtHLy2vmzJmxsbHTpk3z9vYmVyRg9HR1daWnp6ekpKSkpKSnp/f29uL1bm5usbGxM2fOfOKJJzw9PUnRpiMDRIBb+y5evHj+/Pni4mK8ksFgTJ48OS4uLi4ubuzYsRCqHSZUVFRcv379+vXrN27caG1txSudnJyefPLJpUuX4kkTyVWoGx48ePDcc89VVlba29t//PHHmzZtMppbQCaTHT169K233mppafH29j5x4sTkyZPJFjUc4fP59+/fx7cbUcdACDGZzPHjx0dFRUVEREyZMsXBwYFcnYARIJFIioqKsv4lIyOju7sbf2RlZTVp0iT8rB83bhzpFSFdGyB5cnNz//rrr6SkJHlX6ODgMH369JiYmJiYGKJPGWA0lJeX37hx48aNGykpKU1NTXgljUYbM2bMrFmzFixYoKOmXz2jra3t1Vdf/eOPPxBC0dHR+/fvx/PuGTSZmZmvvPJKRkYGhUJZtWrV3r17IelHH2hsbMRhofv37xcXF0ulUryeSqUGBQVFRkZOnDhx0qRJoaGheJYlAHgsDQ0N6enpDx48ePDgQWZmJo/HIz6yt7ePjIycNm3azJkzx40bp1etPWQaIAKBQHDv3r3r16/fvXs3PT29p6cHr7exsQkLC4uOjo6Kipo4caKLiwu5OgEN4PF4//zzD64K3L17t7KykviIxWLhqsDMmTONsulnoFy7dm3Tpk3FxcV0On3lypUffvghi8UiW5QmFBUVbdu27cSJE1KpNCwsbP/+/VOnTiVbFKAE4vZMTU29fft2c3Mz8ZGJiUlAQEBISEhwcHBERERISIiBXo2A1unp6SkpKSkoKMjPz8/KyiooKKioqCA+pdPpI0eOxA/uiIiI4OBgva3T6oUBkofD4dy+fTslJeXGjRv5+flEuzVCyMfHJzIycsKECRMnThw3bhyDwSBRJ9AfYrH44cOHGRkZ6enpGRkZJSUlxDVGoVCCgoJweG/69Ol4znZAnu7u7m+++Wb79u18Pt8QbVBRUdHnn3/+xx9/SCQSJpOZkJCwceNGCCQYClVVVbgSn56e/vDhQ4FAIP+po6NjeHh4WFhYaGjo6NGjg4ODoRAeDkil0srKytzc3Ly8vEePHuXm5paVlck/mikUip+fH44dTpw4cezYsYaSvaB3BkgePp+fk5NDNCUWFBTIf+rm5obrJcHBwSEhISEhIebm5mRJHc40NDTgfyc/Px+/ikQi4lNXV9fx48dHRERERERMnjyZ6A8IqKC1tfXbb7/dt28ftkHz5s176aWX4uLi9DY3SCKRXLly5fvvv8eT4TCZzC1btmzevJnJZJItDdAc+Vs7KyuruLhY/rGHELKzs2OxWLgEZrFYLBYrKCgIXJFB09HRUVFRgQvzioqKioqK4uJiPp8vv42trW1oaCgRHRwzZoyVlRVZggeDXhsgBZqamnBQIT09PTMzs729Xf5TExOToKCg4ODgsLCw4ODgkSNHslgsCwsLstQaJWKxuLKysrS0tKCgIC8vLz8/v7CwUN7uIIRsbW3Hjh078V+IOVKAgdLa2rpz584DBw50dnYihPz8/F588cWVK1e6u7uTLe1/1NTUHD58+KeffqqpqUEI2dvbv/baa2B9jBKBQJCfn09EAvLy8uSbzDBUKtXLy8vf39/f3z8gIMDf39/Hx8fT0xPSq/UNkUhUV1dXW1tbXl5eVlZWWlqKX4VCocKW5ubmQUFBoaGhYWFhOATo4eFBimatY0gGSIGOjg4i5FBQUJCTk0P0JCLAFRQFfH199bZJUn/A9QAFqqurFaqApqam/v7+RBwuODhYn1t8DRGRSJSYmPjDDz9cv34drwkODl66dOlzzz1HTJSjeyorKy9cuPDnn3/eu3cPlyERERHr169fsWIFMTkDYPSIRKLy8nIiVICprKzs+1gxMzNzc3OztLQkYkVubm7u7u7+/v7y4xUBWgeX5A0NDY2NjfgPwstVVVVE/juBqampp6cn/oOIwJ63t7deZS5rEQM2QH2pqqrKz8/Py8srLCwsKysrLy8n+hnJY2Nj4+Pj4+Xl5ebm5uHh4S6Hi4vLsHp4s9nspqam2traxsbGhn+pr6+vqqpSCLBhnJyc/Pz8WCwW9jphYWG+vr562y5jZGRlZf38889nz54lqt3h4eE4myo6OloHCVXNzc23b9/GKXpEe7SHh8fixYvXrVsXFhY21AIAg6Czs5MIJ5SVldXU1FRXV9fU1CjUnQjs7e09PDxcXFxcXFycnJycnJxcXV2JBWdnZ7DU/SGRSFpaWlpaWthsdnNzM15ubGzEC/X19Y2NjUQPawVcXV29vLyw3SHCdV5eXsOqPDcqA9QXgUBQXl5eXl5eUVFBLFRVVREdzRQwMTFxcXHx9PR0dXV1dHR0+Bd7e3t7e3tiQf8zvHp6etrb29va2trb2xUW8K2C7w1ieAYF6HS6l5cX9jp+/8JisaAbM+lIJJLU1NSzZ8+eO3cOtzohhCgUSnBw8OTJk0P/RStdJhsaGvLz83FjR1paGp7WBuPr67t48eJFixZFRkYOqxITGBC5ubnffvvtyZMnu7u7aTRaTEzM0qVLRSJRTU0Nbn+prq5uamrqzxthGAwGYY/s7e1tbW1tbW2ZTKb8K8bW1tYIMu65XG5nZyeHw+ns7CQWiNeOjo4WOVQ/wZlMpqenp7e3N/Y6I0aM8PLywsv6/xTTAUZugJQikUgaGhrq6uoaGxvr6+sb5Kivr8f5FqqxsrLCfsjS0tLCwoLJZJqZmTEYDBsbG3NzcysrKysrK3NzcxsbGwaDQUzpYGZmJl+VsbGx6RtXlMlkHA6HeCsSiYgW2d7eXh6Px+fzRSIRl8sVCAQikYjD4QiFQrwgEom6urqw0eFyuY89CwaD4eXl5erq2tbW9ujRIwqFMm/evNdff93b29vDw8PExOSxewBIRCaTFRQU3Lp1C0dlGhsb5T91dHQMDAzErQxEdJNOp9va2tJoNMLIdnZ2SqVSDofT29tL2GL8WlxcrBAF9PT0nD59+tSpU6dNmzZq1CjdnSpggNy9e/fLL7+8dOmSTCZjMBjPP//8m2++KT/bJQG+9jD4od7U1IQXiJV4djN1YDAY2AlZWFjgMtbOzo5Kpdra2pqYmOCS2cLCwtLS0szMzNraGt8U8iZe4S3eXumxpFKpwvOCz+fL1667urq6u7u5XK5EIuno6MDb9/T04GJcKBQKBAKxWMzj8bq7uzs7Ozs6OvAtqebJOv2Li4uLs7Ozk5OTs7Mzjp85Ojp6eHgYaG6yzhiOBkg1AoGgrq6uqamp7f+iEEfpL3aiD9DpdBysUohdOTg44LAWjnLJzy9x5MiRV199lc/njx8//vjx40rLKUCfKS0tzc7OfvToEW4FrqioUL8YVQqNRvPz8wsLCwsJCQkNDR0/frwRzEELDDUSieTy5cvbtm1LT09HCDk7O2/YsGHTpk2DTILmcrnYFWGLIB8U6ejoUIiRqO+W9BA7OzuF4JbCq52dHfY3Tk5ORhDuIhcwQBrC5/Pb2tpEIhGfz+fxeGKxmMvldnV1icViIirT0dGBwzbEt3CFAC8rBHsw8mEhOp0u71FwMEk+2oRrMFZWVmZmZra2tpaWlubm5g4ODpo1VBUXFy9fvjw7O9vCwmL79u2bN2/WYCeAniAQCKqrq4lwzv79+2tqambOnGliYtLb24urpPgCw9ViMzMzXGV0dXXFTcA+Pj4wrgSgPjwe75dfftm5cydumfX399+4ceP69et13xUXR8QFAgERX+nt7e3s7JRIJBwORyKRcLnc7u7urq4uHIYRi8UKIx5xOBz5JyMuzxWOgtMJRowYQcxjhaNKxAY41IqDT/hew40AFhYWuKHAxMSEiFHRaDTsb4bsVwGUAAYI+B+9vb3btm37/PPPpVLpwoULf/zxR+i8ahx4eXnV1dVxOBwoYQGt09TU9P333+/du7ejowMhFBUVtXnz5kWLFhlr1yHM/v37N27c+PHHH3/66adkawE0BLIXgf9Bp9M/+eSTa9eueXh4nDt3LjQ0NCkpiWxRwGDh8/n19fXu7u7gfgDtkpub+9JLL/n6+n766aednZ1z5sy5d+/e3bt3ly5datzuByFkZ2eHEMKeDzBQwAABisTGxubl5T377LNNTU3x8fGbN2/W54Qn4LEUFRXJZLKgoCCyhQDGw927d+fOnTtmzJgffviBRqOtX7++sLAwMTFx8uTJZEvTEWCAjAAwQIASmEzmiRMnDh8+zGAw9u7dGxUVVVJSQrYoQEMKCwsRQtBvCxg8Uqk0MTExMjJy6tSpFy9edHJySkhIqK6uPnTo0HDrOQEGyAgAAwT0y6pVqzIzM8eNG5eZmTlmzJg9e/aQrQjQhOLiYoQQicNGA0YAj8fbs2ePr6/vvHnz0tPT/f39d+/eXVVV9cknnwzPTEFsgJQOGAsYCmCAAFUEBgY+ePAgISFBLBZv2bJl0aJFbW1tZIsCBgYewBCawADNaGpq+uSTT7y9vbds2VJTUxMVFXXq1KmioqLNmzcP58kW7e3tEUSADBzoBQaoRUpKyqpVq+rr611dXX/99dcnn3ySbEWAuoSGhubn59fU1MDEtMCAyM3N3b9//5EjR0QiEZVKjY+Pf//994dPlo9qenp6zMzMnJ2dlU64BBgEYIAAdeFwOBs2bDh58iSFQtm0adPXX39NDHIN6C29vb14xBEulzus5rkDBoP64zgPZ6ytrbu7uw163MVhDjSBAeoCmdGGSEVFhVgsDgoKAvcDPBbIcR4QdnZ23d3dCoMoAgYEGCBgYEBmtGEBCUCAOkCOswZARzBDBwwQMGAgM9qAwAYIuoAB/QE5zhoDHcEMHTBAgCbAmNGGAjZAMAgQ0JfhPI6zVoAIkKEDBgjQHBgzWv+BJjCgLzCOs1YAA2TogAECBgVkRus5xcXFNBrN39+fbCEA+UCOs3bBBojD4ZAtBNAQMECAFoDMaP2kubm5vb2dxWKZmZmRrQUgE8hxHgqsrKwQQnw+n2whgIaAAQK0A2RG6yHQ/gVAjvPQgX9A6AZvuIABArQGZEbrG2CAhjOQ4zzUWFpaIjBAhgwYIEDLQGa0/gAGaHgCOc66ARsgoVBIthBAQ8AAAdoHMqP1BDBAwwrIcdYxEAEydMAAAUMFZEaTDoyCOEyAHGdSYDAYCAyQIQMGCBhCIDOaRIRCYU1NjbOzMzwCjRjIcSYRiAAZOmCAgKEFMqPJori4WCqVwhjQxgrkOJMONkBdXV1kCwE0BAwQoAsgM1r3FBYWIkgAMkYgx1lPMDU1RQj19PSQLQTQEDBAgI6AzGgdU1xcjCAByIiAHGd9w8TEBIEBMmTAAAE6BTKjdQZ0ATMaIMdZPwEDZOiAAQJ0DWRG6wbcBAY5QAYN5DjrM2CADB2KTCYjWwMwTElJSVm1alV9fb2rq+uvv/765JNPkq3IeJBKpcRERVQq1HMMj9zc3P379x85ckQkElGp1Pj4+Pfffx+yfPSKyspKFos1YcKE9PR0srUAmgAlI0AakBk9dFRVVQmFwsDAQHA/BgfkOBsKdDodQQTIkIHCESATyIweIiAByOCAHGeDA5rADB0wQAD5QGa01gEDZEBAjrOBAgbI0AEDBOgFkBmtXcAAGQSQ42zQUCgUsiUAgwIMEKAvwJjRWgQMkJ7z6NEjGMcZAMgFDBCgX0BmtFYoKiqiUqkBAQFkCwEUwTnOo0ePhhxnQwf3oYY4kOECBgjQOyAzepC0t7e3tLR4e3vjuYoAfQBynI0PMECGDhggQE+BzGiNgVnA9ArIcTZWwAAZOmCAAP0FMqM1A8aA1hMgx9m4AQNk6IABAvQayIzWAJgGlXQgx3k4AAbI0AEDBBgAkBk9IKALGIlAjjMAGApggADDADKj1QcMkO6BHOdhiFQqRQjBbDOGC/xzgCEBmdGPRSwWV1ZW2tvbOzs7k61lWAA5zsMWsViMEDIzMyNbCKAhYIAAAwMyo1VTUlIikUggA1oHQI7zMAcMkKEDBggwPCAzWgXQ/qUDIMcZQGCADB8wQIChApnRSsEGCLqADRGQ4wwQYANkampKthBAQ8AAAQYMZEb3BRsgaALTLpDjDPQFIkCGDhggwOCBzGh5oAlMu0COM9AfYIAMHTBAgDEAmdEYmUxWUlJiamrq4+NDthaDB3KcAdWAATJ0wAABRsJwy4wuKSl57bXXDh48mJKS0tDQgFfW1tby+fyRI0fS6XRy5Rk0kOMMqAMYIEMHSknAqMCZ0Rs2bDh58mR8fPymTZu+/vpro8xSdHR0/O6774i3tra2gYGBTCYTIWRtbV1cXMxisUxMTMgTaJDcvXv3yy+/vHTpkkwmYzAY69evf/PNNyHLB1CKUChECEE40HABAwQYGzgz+qmnnnr11Vf37t17796948ePG98zzN7e3tPTs66uDr/t7OxMT0/Hy2lpaUFBQSYmJiwWKzg4ODAwMCgoaPLkycb3I2gLqVR66dKlbdu24d/Q2dl5w4YNmzZtgiwfQAVcLhchZGNjQ7YQQEPAAAHGyapVqyIjI5cvX44zo7dv375582ayRWmZ0aNHEwaoLz09PcXFxXhiVITQhQsXwAD1hcfj/fLLLzt37qypqUEI+fv7b9y4cf369VCtBx4LNkDW1tZkCwE0BHKAAKNloJnRYrF469ateIZng2D06NFqbhkdhsm2jwAAIABJREFUHT137twhFWNwQI4zMEh4PB6CCJAhAwYIMGYGlBn9/vvvf/XVV/v27dOlwsEQHh6u5pY7duwYUiWGBeQ4A1oBmsAMHTBAgPGjzpjR165d27VrF0LonXfeyc/PJ0PmgFEzAjR//vyoqKihFmMQ9B3HuaCgAMZxBjQDIkCGDhggYFigeszo1tbW1atX48YvkUi0YsUKg5hVIyAgwNLSUvU2NBrtv//9r2706IycnJwBba9iHGeYMwTQGIgAGTpggIBhRH9jRq9bt66xsZHYLCcn56OPPiJJ4wCg0WghISGqt1m9evVjtzEs9u/fP2XKlKamJnU2hnGcgaEDDJDBIwOAYUZPT09CQgKVSkUILVy48Ouvv1a4KSgUCpVKTU5OJlvp41m3bp2Ku9vMzKy6uppsjdrks88+w6f28ccfq96ysbExISHBzs4Ob49znHt7e3WjExgOREZGIoTKysrIFgJoCBggYJhy7do1d3d3hBB2Qn3x9PRsa2sjW+Zj2Lt3rwoD9Pbbb5MtUJts3boV21OEkKOjY1dXl9LNcnNz169fb25ujv9cnOOsY6nAcAC3n7a0tJAtBNAQaAIDhilxcXGZmZm2trZSqVTpBnV1devXr9exqoGiIg/a1tYWOwYjQCaTvf76619++SWFQpHJZAih1tbW3377TWEzyHEGdAmbzabRaPb29mQLATQEDBAwfPn66687OztVbHDmzJnjx4/rTI8GhIeH44hIX959913jSHORSCTr1q3bvXs34X4w33zzjUQiQZDjDJBBT08Ph8NxdHTsL4QM6D//p0ABgOHDtWvXZs+ejRBSfQvY2Nj8888/+jy5uo+PT3V1tcJKd3f30tLSx/YR038kEsmaNWuOHDmi4H4wx44da21thXGcAd3T2Njo7u4eGhr66NEjsrUAGgJTYQDDEfl+76rhcrkrV668efOm3g6RFx4e3tcAffLJJ0bgfrq7u5977rmzZ8+ifnzqhg0b8Fgs06dPf/PNN+fMmdNfPAwAtEtLSwtCyMnJiWwhgOZA7A4YdshksjVr1sj3e1fN3bt3v/rqqyGVNBj6pgGNHDnyhRdeIEWMFhEIBPPmzcPupz94PF5sbGx6evrNmzfnzp0L7gfQGWw2G4EBMnDAAAHDjoMHDyYmJqq/PYVCSUhIyM7OHjpJg6GvAdq+fTudbtjB3a6urrlz5169evWxW1pbW0+YMEEHkgBAHogAGQFggIBhB5PJXLRokfpzOMtksp6enhUrVggEgiEVphkKM4JNnDhx4cKFZInRChwO54knnkhJSVFn4wsXLhQWFg61JABQAAyQEQAGCBh2LF++/MyZMx0dHXfu3Nm6dWtERIQ63yosLHznnXeGWpsG+Pv7MxgM4i3uK06inkHCZrNjYmLS0tLU3F4mk+3evXtIJQFAX8AAGQFggIBhCo1Gi46O3rFjR2ZmZkVFxaFDh5YuXWplZdXf9hQK5cCBAxcvXtSlSHWgUqmhoaF4OT4+fsaMGaTKGRRNTU0zZ84c6FRfhw8fVnNmDADQFs3NzQgMkIEDBggAkK+v7/r160+dOtXc3HzhwoWXX355xIgRCtvgkUPXrl2La356BU4DolKp27ZtI1uL5lRXV0dHR+fl5Q30i2Kx+MCBA0MhCQD6A4+80LegAAwIGAcIAJBEIuFyufKvvb29ubm56enpDx48KCgowAPuYcLDw9euXSsUChFCUqmUGEoRfxchJBKJ+n7K4/F6e3v7E0B8VykCgUAsFqvQLxaLBQKBqampi4uLiYlJf5uZm5urGCDHwsICTx9Bo9GI+R1tbW3xOG8D/ZRKpeJPiVf5L/altLQ0Li4OP1Q0wMHBoaamxgh6/gOGQkhISEFBQUNDg5ubG9laAA0BAwQYJF1dXV1dXXw+v7OzUyAQdHV1dXZ28vn8np6ejo4OhBDxKpPJOByO/Cv2JfKvZJ/NMALbIDqdbm1tTbxSqdTs7GzVJg9DoVDs7OwcHR3d3Nzc3d1dXFxcXV1dXV2dnJymTp2qfmI7AAwSGxsbsVgsFAphJGjDBQwQQBoikaijo4PD4WDvwuFwsJXhcrk8Hg8vd3R04AUej8flcvEyh8PRlgYzMzNLS0scGpF/JUIaDAbD1NQUIWRqamppadnc3FxVVdXQ0DB//nxHR0f5T4lMZGIGcqWfKmBpaWlmZtafPBVflIfL5SYkJOzatQu/xVavv42xTVT6EXaNCKGurq7u7m6EUHd3d1dXl2afygfVFBxnf/OvDRQTExMrKytbW1sGg2FpaWlra2tlZcVgMBgMBpPJxCttbGysra3xsp2dHfNf9HZkS0D/6ejosLe3Z7FY5eXlZGsBNAcMEKA1uru7OX3AFkcpIpFIg6PgEIKNjY2lpaX8c87a2trW1tbU1NTa2ho/F+VfsZOQf8XWR+OTlUgkevUE5fP5KjK49RAiIJecnNzd3e3g4EClUplMptJwnUgkEggEHA6Hz+cLBAIc+evq6hIIBBoH8KytrZn/F3l71Pcj7Z4+YND8888/Y8aMmT59+s2bN8nWAmiOYY+WBugGPp/f0tLCZrNbW1vb2tpaW1vxW7xMWJwBDZPDYDA8PDyIBwxhZRgMhkKd3traWt7rqIiX6BK9cj8IIcNyPwghJpOJELKzs1u3bt0gd6UQLxQIBFwuVz5eSDSYKljw2tra2tpa9QUTZsjJycnJycnR0dHBwcHR0dHZ2dnJyQkv68n1CQwp+LKBDGhDBwzQsKanp4fwNGw2u6WlBS9ji9PS0oKX1UnOsLCwcHNzU1GBVlijIlcXAAaEtbW1Ztk/2B71pbOzs7/4pTpiCHtE4OzsTLglJycne3t7DdQC+gN0ATMOwAAZOSKRqLGxsaGhobm5ub6+Hr+y2ey6urrm5mY2m/3YNlAGg+Hq6uri4kKU4ESBjsHmBiq+gMGBI44eHh5qbo9jSO3t7bjC0CoH8batra2ioqKiokLFfnB/PU9PT2dnZ/zq4eHh6urq7u6O7zXIq9VzcATIy8uLbCHAoAADZPDw+fy6ujrC0zQ0NDQ1NeHXxsZGIjVVKTh0r+BpsMshHI+KjtMAMKywsrKysrLy9PRUvRmXy21ubiZiq33dUkNDg4rWNxqN5uzsjLu5Ea4Id3nz9PRUPdIBoBtwBAgMkKEDSdAGQ0dHR0VFRUNDQ2NjY8W/4LcqvmVnZ4dLUoVXOzs7b29vg0scAQDjQCgUdnR04Ohs39fm5mYVHeXs7OxYLBaLxcK3M+tfIFNbZ0yYMCEzM7OkpCQgIIBsLYDmgAHSL8RicW1tbU1NTU1NTVVVVXV1NV6ura1VmohDo9Hc3Ny8vLxwLN3FxUW+1gixdAAwREQiEY7jyrdc49eGhob+xiJnMpkjRozw9vb29vYe8S/e3t5ubm4GPT2cHsJkMnGKPUTjDBowQKTR2dlZVlZWWlqKX0tLS6uqqpqampT+I5aWlj4+PiPkwG89PDzodGjHBIBhhFAolK8d4cpSTU1NfX290tHGTU1Nvby8WCyWv79/wL/4+vriQaqAgdLY2Oju7h4YGFhUVES2FmBQgAHSBWKxuL6+Pj8/v6CggGi9qqys7PvjEy1WRFgbv/X19YU6HAAAqsEN5fLt4xUVFWVlZUpHS3JzcwsJCWHJERwcDDl/j+XWrVszZsyYN2/eX3/9RbYWYFBA8EDLSCSS8vLyR48elZSU4OhOSUkJnjdYHgqF4uXlhStkRLXM19cXSh8AADTGzs4uIiIiIiJCYX1bW1t5eTkONpeUlOCFxsZGhQxCOp3u4+NDlEijRo0KCwtzcXHR4RkYAMXFxQihwMBAsoUAgwUM0GDp6OjIz8/PysoqKCjIz8/PyckhJgcgwEmLwcHBRH0rMDAQEpABANANDg4ODg4OEydOlF8pHy6qqKjIz8/Pz88vKysrKytLSkoiNrOzs8NlV3BwcERExJgxY4Z52VVSUoIQGjlyJNlCgMECTWADg8/n5+fn5+bm5uXl5eXl/fPPP21tbfIbWFhYBAcHh4WFBQUFEdEdiOsAAGAQtLS04Lh1aWlpQUFBbm5uZWWlfJc0KpXq6+sbFhYWGhoaFhYWFhYWEBAwrDIR582bl5iYePv27alTp5KtBRgUYIAeQ0tLS2ZmZmZmZnZ2Ni4L5H8xGo3GYrHCw8NDQ0NDQ0PDw8P9/Pz0bZIEAAAAjenq6sJOCNf6cnNz2Wy2/AZmZmbBwcHh4eHjx48fP3786NGjjbvKFxgYWFJS0tTUBI2Dhg4YIEV4PN6DBw8yMjKw78EDXhG4urriSg+u/YSEhBj3rQ4AAKBAS0uLvB/Kz8+Xb/en0+mhoaERERHjx4+PjIwMDw83pjphT08PnqZQnXlRAD0HDBBCCFVVVaWmpqalpd29ezcvL08ikRAf+fj44GoNbvx2dHQkUScAAIC+IZVKKysrHz58iCuNWVlZ8ubAyspq0qRJU6ZMmTJlyuTJk21sbEiUOnhyc3NHjx49ZcqU1NRUsrUAg2UYNdwqUFZWduPGjRs3bty6dauhoYFY7+DgMGXKlMjISFyDAccDAACgAiqV6ufn5+fnt2TJEoSQTCYrKyvDZigtLS0rK+v69evXr1/HW4aGhsbExMTGxk6bNo3JZJKtfcBkZ2cjhMaNG0e2EEALDK8IUHNz89WrV1NSUm7cuEG0bVEolMDAwClTpkRFRU2ZMiUwMBBG3AEAANAKIpEoMzPz3r17OMpODGNNo9HGjh0bExMzc+bM6dOnm5ubk6tTTbZs2bJnz56ffvpp7dq1ZGsBBovxGyCpVJqZmXn58uVLly5lZ2cT3RmCgoJiYmJiYmJmzJjh5ORErkgAAIDhQEFBARF6b21txSstLS1nzpwZHx8fHx8/YsQIchWqZvr06bdv387KyoIgkBFgtAaop6cnJSXl9OnTFy5cIPosuLi4xMfHx8XFxcTEuLm5kasQAABg2CKVSh89enTjxo2kpKRbt26JRCK8PiwsbNGiRUuXLg0JCSFXYV9kMpmdnZ1AIOByuYYSsgJUYGwGSCKR/P3336dPnz5//nx7eztCiEqlRkREPP30008//fS4ceNgclAAAPQHkUgk35N09uzZ8oMQDhMEAkFycvLly5cvX75MJCcEBwcvWbJk2bJlwcHB5MojKC8v9/f3Hz16dE5ODtlaAG0gMxZKS0vfe+89d3d3fF40Gm3GjBn79+9vbGwkW5q6fPDBByr+KSsrq+Dg4FdeeSU/P59spUaOSCT6+eefn376aU9PT3NzcxMTE0dHx3Hjxq1YsWLv3r33798XiUQKX2lvb4+IiHB3d8/IyCBFM+nALzBInn76aYTQ7NmzB78rg/4vsrOz33///YCAAKLomzRp0o8//sjlcsmWJvvzzz8RQqtXryZbCKAdDN4ASSSSs2fPzpgxg8hcjo6OPnjwIJ5W3UDBJ7J48WJijUgkevjw4fLlyxFCdDp9165dJMozboqKikaOHEmn09etW5eent7a2ioWi6urq48ePUrURLdu3arwrdOnT+OPtmzZQops0oFfYJBo0QAZx3+Rk5PzwQcfEClBVlZWa9euzcvLI1HS+++/jxDavXs3iRoALWLABkgsFh8+fHjUqFH49rC3t1+/fn1ubi7ZurRAXwNEsHr1avzpn3/+qXthRg+fz8cF7oEDB/p+yuPxwsPDlRqgtra2sWPHurm5PXjwQCdK9Q74BfqDmHD06NGjKjbTogEypv9CIpFcu3Zt5cqVuKGQQqHExcVduHCBFDHx8fEIoZs3b5JydEDrGKQBkkqlR44c8fT0xMVKRETEyZMnu7u7ydalNVQYoPLycvxpUFCQ7oUZPfv27cN1zf4up7/++kupAQKA/tC9ATJK2Gz2p59+SvTYjYmJefjwoS4FSKVSR0dHKpXK4XB0eVxg6DC8jOD79+9HRkauWrWqrq5u1qxZycnJmZmZy5YtMzExIVuaLmCxWHgCmqKiIqJgBbTFnTt3EEKOjo79XU6xsbEwTBQA6B4nJ6ePP/64urp6//79vr6+N27ciIiIePHFFxUmJhs6CgsLW1tbw8LCbG1tdXNEYKgxJAMkkUg++eST6OjojIyMoKCgS5cuXb16NTY2lmxduoYYPpUYUgzQFp3/r737jmvyWh8AfkJCGCEshbBEhpMAooiD4URABEWv24KjrbVaodrealut2nrb2v5uFa+z1lVxVLAOEEFws0RREQFFZK8wQwgBQkJ+f5zb96YJIjNvCM/3j3zCedfznpDkyXnPe059PUKopKSkpaWl3RV0dHTa2tp+/PFHxcYFAEAIIS0trfXr17969Wrfvn1MJvO3336zt7e/ceOGAg6dkJCAEHJ3d1fAsYBi9JsEqLq6evr06bt27aLT6aGhoRkZGfhy7ABEzLODW4Obm5spUnx8fBBCT58+XbhwoZmZmZqaGoVCMTExkd5DSUnJP//5T0dHR11dXS0tLVtb2zVr1uAh3tsVFRXl7e09aNAgDQ0NCwuLJUuWJCUl3b17V/q4O3fu7FIkDQ0NERERQUFBeIY1vOeFCxfevn1b5ujyuxWLxaGhoU5OTtra2oMHD542bdqlS5fwymlpaQsXLjQxMdHQ0LCxsfn888/5fH7n6xa3rolEop9++qmTm7R71r1Yh2KxeP/+/c7Ozkwmk8FgODk57d+/XyI1ekVMTIz0JmFhYcSiffv2SS+qqKggFqWkpEgv2rp1a0lJSXBw8KhRo3R0dHR0dJydnfft29fa2tqTGujG6XTpQL31n9CZd0Q3zgVXMjHkWGBgILF5B3Ps9G4VdWarfkddXT0kJCQnJ2fZsmVVVVV+fn7ffPNNXx8UT/4FCZBKIfcKXCdxuVxnZ2eE0KhRo549e0Z2OH0OvzRd7QNE9CGIi4uTmaOexWIRq124cEFLS8vQ0PDkyZM1NTV8Pj8uLo7NZlMolK+++qqtrU1mt8HBwQghU1PT8PBwHo/H5XLDw8MtLCwWLVqEd56RkdGNSDZs2IAQMjMzO3/+fGVlpUAguHfvnqurK3p7Dxu8W09Pz3nz5u3bt6+2tpbD4fz73//Gc03/8ssvly9fnjFjxrNnzwQCwcOHD+3t7RFCU6dOFYvFnaz5M2fOEKGuXr369evXndxQ0mEfjm7Xoaenp5+f3/fff8/hcJqamqKjo3HWGxwcLLN+Bx1N8L27CKF2h4QwNzdHCM2aNcvIyOiHH34oLS1tampKSUnBr4Wbm1t9fX3Pa6Crp9OZA/XWf0JX3xG9+NK8c7cxMTE4Ke9eFQ2EHkWnT5/W0dFBfX/Lm42NDUKoqKioT48CFKkfJEAikcjNzQ0h5Orq2tDQQHY4itBBAtTBXWD4I8/V1dXc3PzLL78sKChoamqKj49nMplE2hEbG4tbYhITE6W3raysxJ+zu3fvli4/ePAgQohKpcr0N8zKyqLT6R1/eXccCU6AUlJSpDfk8/m4b3u797jh3aqpqckEGRQUhBBiMpnDhw+X/rbOysrCEd64cUN+b+1qaWkZO3YskQNRKBQPD4/9+/cXFxe/c9u3feX0pA7V1NS2bdsmXX7kyBG8t9LSUunyHiZACKGDBw9Klzc0NNja2iKEFi1a9M5zl465gwSo86fTmQP1yn9CV98R3TiXLiVA8rs9duxYt6toICRAEokkNTUV98v5+eef++gQ+EW0srLqo/0DUvSDBOjo0aMIIUdHx4HT914+AWppaXn27NmKFSsQQjQa7ZdffpHfCn/kIYS++eYb6fKPP/4Ypx2tra3W1tYIoTlz5shv/v333+Od5+fn4xKBQGBoaIgQWrBggfz6OJgOvrw7iEQikWzYsGHMmDHyu8XdayZOnPi2E6RSqTU1NdLlRLPNP//5T5lNcDolX96ByspKT09P9HcUCmXixIn//ve/q6ur37Zhu185PaxDNTW1yspK6fLc3Fy8yYULF6TLe5gAsVgskUgkswinbgih5OTkt521fMwdJ0CdPJ3OHKjn/wldfUd071y6mgDJ7JZo9O1GFQ2QBEgikSQmJtLpdG1t7T5qobl48SJC6L333uuLnQOyKHsfIKFQuH37dgqFcuTIkYHW9/7SpUvEtXwNDQ1XV9cnT56sX7/+2bNnmzZtettWdDp98+bN0iWHDh3CnT/i4uLy8/MRQu12n8KFIpHo+PHjuOTGjRt4OhEvLy/59adOndpx/B1EghA6cOBAu8PJ44GdHj9+LBAI2t3t8OHDcUpBIPpYTJw4UWZl/O1eVlbWcajSjIyM4uLioqOj58yZQ9wLJpFIHj58+Nlnn1laWv7444+STk8g08M6HD58uMxMvcS4cL17D6CHhwe+fiTN29sbP5G+MtgTvXs6Pf9P6Oo7QuboffHSyO92yJAhPd+tynN1dQ0JCREIBLt27eqL/eMOQPhaBFAZyp4ApaamVlZWzpgxY/LkyWTHomgyl8AaGxuzsrIOHjzY8RyBdnZ2b8sU8V0MCKFRo0bJLyUK8a3gCKHU1FT8ZOTIkfLrE+MwdSOSDujq6iKExGKxdI9dafKz2GpoaLxtEZ6wkJhnsfNmz54dFRVVUVFx8uRJPz8/4hACgeDLL7/EfXo6o4d1KH9G6urqOC3rxkl1oN1IrK2taTQakjqLHurd0+n5f0JX3xEdH71XXhqFveKqZ8uWLVQqNSoqqi92DreAqSRlT4DwL7Z2P55Auzq4u4T4BWlgYCC/VENDA3dYJn4lczgc/KTdPEamd3OXIsFu374dGBg4YsQIBoNBtHVNnz4dL21qamp3K6LfTJcWdY+hoeGqVasiIyOrqqpOnDhBJDEHDx58+fJlZ/bQwzokvtGl9cWEvm87EJPJRFJn0UdH6d7eev6f0NV3hMxS+cKevzQKe8VVz6BBg4yNjTkcTmduXeySurq69PR0Q0ND5ZmWFfQKZX9f4Q+mkpISsgPpN+QvZBA6eeGG8tdAf8T6lPaG/nvn3jqIBCG0bt26mTNnXrlyZe3atS9evCBGXr5z506XolUMJpO5evXq9PR0/BNQIpHEx8d3ZsMe1qHCKE8kitTVdwRQZgKBoLa2VldXt9cHxb1586ZIJPLy8oJMVMUo+8s5efJkbW3tmJgYoicg6DYzMzP8pK6uTn5pS0sLbnQhGuGJMXvaXb+xsbHbkfz555+4b/tPP/30+eefW1tb94uBvDU0NL7++mv8HPfseae+q8O3aWtrkyl5W1uatHZ7XLW1tTU0NCCps1AxXX1HAGV27NixlpaWvhgaFw+0OHv27F7fMyCXsidAOjo6n3/+eUtLy4YNG8RiMdnh9G8eHh74SXZ2tvxS4poOsRrRkzQnJ0d+/YKCgm5HEhsbi58o20Btu3btolKp7X4dYvimIfTXKJTv1Hd1KEO6l5LMos70AW83kry8PJFIhBCaMGFCT+NTSl19RwClVVhYuHPnTiqV2usjIkokkps3b1IolFmzZvXungHplD0BQght2bLF2to6NjZ21apV8r9uQed5enrisbyio6Pll+JCGo22Zs0aXDJ79uxBgwYhqXxF2q1bt7odCfE6yl+D6NIdW71OIpG0tbWlpKS8bYW8vDz8ZNq0aZ3ZYd/VoQw9PT1tbW2EUHFxscwifANLxxISEoRCoUwhEXNgYGBvxKh0uvqO6B64btLXSkpKZsyYweVy165dKz2OV694+vRpeXn5uHHjoCFQ9fSDd6a2tnZsbKypqWlYWNj8+fNramrIjqi/otFohw8fVlNTu3HjRlJSkvSiqqqq0NBQhND27duJRg5NTU08FMq1a9dkbll/9uxZT+62mDlzJn6CR9eQdvbs2W7vtrd888037V4Samlp+de//oUQWrx4Mb5d/536rg5lqKmpTZkyBcklVdnZ2Tdv3nzn5vX19fv375cu4fP5e/fuRQgtXrxY/q5y1dDVd0T36Ovr45vpiLk4bt26RaFQiMupoCdSUlLc3d3z8vL8/f3xS9a74PqXCusHCRBCaPjw4fHx8ZaWlteuXRs7duzdu3fJjqi/8vLyOnfunKam5ty5c0+fPl1XVycQCOLj42fMmMHhcL744ovt27dLr7927dpNmzaJxWJfX99Lly7x+Xwej/fHH38sWLBAZs0uWbJkyfz58xFCO3bs2L17d35+vkAgSE9PX7p0abstJQr2+PFjFxeX06dPFxQUCIXC1tbW4uLic+fOTZw4MTExcd68eSdPnuz83vqoDuXt2LGDTqc/fPhw27ZtVVVVzc3NcXFxgYGBW7dufee2a9asOXny5L/+9a+ysrKWlpbU1FQfH583b964ubnhkYhVVVffEd1Ap9PxKFB//vlnTU1NZWXlvn371NXV//GPf/TGGQxcIpFoz549U6ZMKSwsDAgICA8P74uuhJAAqTLFjLfYK7hc7pIlS3DYfn5+b968ITuiXva2X4TvHM6VmM1A2ttm1JJIJMXFxZ999hmbzWYwGJqamtbW1itXrkxNTX3b+lFRUV5eXgYGBhoaGkOHDl29enVOTg4xa2lWVlY3IhGLxcePH58yZYqenh6VStXV1XV2dt6+fTsxbDGG58SQ3y2bzZZIJOfPn5c/3OXLlyUSCZ7GQZqtre27XgFJc3NzUlLSL7/8smTJkvHjx1tZWTGZTBqNpq+vP2bMmA8//PD27dsym7yti3Fra2tv1SE+2XZ7Jq1YsUL6KCkpKX5+foaGhurq6kOGDPnwww+LiopkqvTYsWPSm+DDbdmyhcfjbdu2zd7ensFgaGtrjxs3bu/evcTdeR3ouAZ6cjry+ug/oZPviG6fS2Vl5apVq0xNTWk02qBBg6ZOnRoTE9Mrr3iXXgtVEhcX5+joiBCi0Wg7duzo/Hx/XVJbW0uj0QwMDFSvAoGkX0yFIePYsWO4/6mmpubmzZtLSkrIjmiAioyMxJ+tXZqiCEhThjokEiCyAgCgS+7fv0+MUT5hwoTHjx/33bEuXLiAEFqyZEnfHQKQqH9cApP2wQcf5OVw9QpGAAAgAElEQVTl7dixAyH0yy+/WFtbBwUFZWZmkh3XgFNYWIgQ0tfXV9V7pBUA6hCATmpra4uMjHR3d58yZUpsbKy5ufnRo0eTk5OdnZ377qDXrl1DcP1LdfW/BAghpKOjs3PnzpcvXwYHB9Pp9DNnzjg4OEyfPj0sLKwzQ56AzouKiqJQKPI3FiGEcE8sX19fuMmlY1CHAPRESUnJd999N2zYsLlz5yYmJtrY2Bw8ePD169dr167t0zeOQCC4du0anU739/fvu6MAEvXjj92hQ4eGhobi4R9YLNbdu3cDAwPNzMw2bNiQlJQkGZAj2/aR3377Tabk+fPnV65c0dLS2rZtGykh9TtQhwB0SWNj4x9//OHr62tlZfXNN9/k5+e7uLhcuHAhJydn/fr175xGpueioqL4fL6Pj4/MnLtAdZB9Da53iESiuLi4RYsWEXcBWFhYBAcHP3jwoI86xw0QuJMKlUrdvn3769evm5ubKysrT506xWKxGAzG1atXyQ6wH1DyOoQ+QECp4HaXwMBAHR0d/GGur6+/du3atLQ0BUeyYMEChNDZs2cVfFygMBSJarWUlJWVnT17Njw8/NGjR7jE3Nzc19fX19fX09OTeEeBThKLxQkJCVeuXElKSiovL6+srKTRaNbW1l5eXiEhIZaWlmQH2A8oZx2mpKRMnjxZppDNZr948YKUeMAAV1RUFB0dHR0dfevWLTwKl5aWlo+Pz5IlSwICAtqdI7ZPNTQ0sFgsiUTC4XB0dXUVfHSgGKqWABHy8/MjIiKkMyENDQ0PDw9fX99Zs2ax2WyY4BAAAEjU1NSUnJwcGxsbHR1NZN6ampo+Pj6LFy/28/NjMplkxRYWFhYYGLho0SL5wVqBylDZBIhQWVkZExMTFRV18+bN+vp6XGhkZDRx4kR3d3dPT89x48ZBMgQAAAogEonS09Pj4+Pj4+MTEhKam5txOYvF8vLy8vf39/b2VoYWF39//6ioqIiICBivUoWpfgJEaGlpefDgQUxMzO3bt9PT04npqMzMzNzd3V1dXV1dXceOHYsHrQcAANAruFxucnJycnJyYmJiSkoKMc8Mg8Fwd3efOXOmr68vm80mN0hpXC6XxWJpaGhwOBwF9LYGZBlACZC02trae/fu3blz5/bt23gQXlyura3t4uLi5uY2ceLE8ePHm5mZkRsnAAD0O2Kx+OXLl48ePUpOTk5KSsrKyiJ+cGpqak6ePHn69OkzZsyYMGFCX0xe0XMnT55cs2bNihUrwsLCyI4F9KEBmgBJ4/F4qampCQkJiYmJiYmJ0iMJmZqaOv9l4sSJxsbGJMYJAABKq6ysLO0viYmJ0vN4sFgsFxcXZ2dnd3d3d3d3TU1NEuPsDA8Pj4SEhMjISD8/P7JjAX0IEqC/aW1tTUtLS05OfvTo0ePHj3Nzc6Xrx8LCgs1mOzo62tvb29vbs9lsxd+bAAAApKupqcnIyHjx4kVGRkZGRkZmZiaPxyOWamhojBkzxsXFZcKECa6ursOGDSMx1K7Kysqyt7c3NzcvKCigUqlkhwP6EPR3+Rt1dfVJkyZNmjQJ/8nlch//JS0traCgoKSkhJiunEajDRs2zMHBwcHBwd7e3sHBwcbGBob0BQComKampqysLOmMp7y8XHoFOp3u5OQ0fvz48ePHu7i4ODg4KOe1rc44dOiQRCL58MMPIftRedAC1AX19fUvXrx48eLF8+fP8QeBzHTN2trabDYb50OjR48ePnz40KFDoVc1AKAfaWxsfP369evXrzMzM/EH3Zs3b8RiMbEChUKxsrLCDeG4RXzkyJH9N+OR1tjYaG5u3tjYWFBQgMcIBSoMEqAeqaury8zMTEtLy8rKyszMfPr0KXGDA6aurj5kyBCbv9jZ2bHZ7KFDh8JvCwAA6YRCYUlJSV5eXmZmZlZWVl5eXl5eXkFBAdFnGdPX12ez2Ww2G3+COTk5DR48mKyY+9SxY8fWrl37j3/8IyIiguxYQJ+DBKg3iUSi3Nxc3ESck5ODf0U1NDTIrKahoWFjYzNixIjhw4cPGzYMPw4ZMgSGIwIA9JGWlpa8vLycnJzc3Fz80ZSbm1tcXCz/FWBqajpixIhhw4aNGjUKN/AMnPthx48fn5aWdvPmzVmzZpEdC+hzkAD1OQ6HQyRD+KMnNze3sbFRZjVNTU0bGxtLS0tLS8shQ4YMHTp06NChlpaW5ubmqtG2DABQgMbGxsLCwsLCwqKiouLi4qKiooKCgqKiopKSEunLWJixsfFwKfj32ICdMig1NXXixIm2trY5OTnQm3MggO4pfY7FYrFYLA8PD+nCuro63NpMtDy/evUqKysrKytLfg8GBgb4CpqpqamZmRlxQc3AwEBRJwEAUC74M6SsrKy8vDzvL2VlZRUVFe3+rCU+RojL8Q4ODnp6eoqPXGkdOXIEIfTRRx9B9jNAQAuQspBIJBUVFYWFhfhHG/G7raioqLa2tt1N9PX1cYuRubm5iYkJTo9MTEzMzMxYLBZ0vgagXxMIBOV/V1ZWVlpaij8WWlpa5DehUqkmJiZWVlaWUqysrIYMGQK5TsdqamosLS3FYnFJSYmq9nACMuA7UllQKBRTU1NTU1PiJnxCc3NzWVmZ/K+9oqKi58+fP3/+vN0dGhgY4JTI1NTUwMAAP8GPlpaWJM4yCADAmpqacFqDH+vq6qT/fFtbDvqrH6FMq7Cpqam1tbW2traCz0I1HDp0SCAQrFy5ErKfgQNagPoxkUhUVlZWVFSEPzErKirwhyZ+rK6u7mBbQ0NDnG+xWKzBfzE2Nh48ePCgQYPwn3CrGgDd1tTUVFNTU11dzeFw8BOssrKSw+HgFEd63HkZ6urqxsbGuHGXaNnF6Y6lpSV8SfcugUBgZWVVXV39/Plze3t7ssMBCgItQP0YjUbDTdztLhUKhdXV1e3+siwvL8dX1jIzMzvYv6amJtF0ZPAXmT9NTU3h5jUw0DQ1NRHvqbq/yPzZQfsNpqGhQTTKSj/idxkMlqFIJ06cqKqq8vf3h+xnQIEWoAFKLBZzOJyKigoOh1NdXY1/oVZVVVVVVRE/VWtqamSGA5GnoaFhYGCgLwcXyi+CO9qAEmpsbOR2qK6ujnjkcrmdeV/gZlQjIyMjIyOiYRX/yWKxcKKjmLMDHROJRCNHjszLy3vw4IG7uzvZ4QDFgQQIvJVEIiEyIaL1Xjo9wm37fD6/8/tkMBjtZks6Ojp6enra2toMBkNPT09HR0dbW1tHR0dfX19bW1v5Z08ESoLH4wkEAoFAUFdX19jYKBAIGhoaeDwefi6dxEgnN62trZ3cP4VCwakMcaUYXzgmSoyNjY2MjAbsneT9EZ773c3NLSEhgexYgEJBAgR6SiQSyXyvvO051tzc3NVDqKmp6enpMZlMnCHp6+szGAxtbW1dXV0mk4mfGxgYaGtra2ho6Onp4fXV1NT09fUpFAr+qQ0/uJWZUChsbGzEjy0tLQKBAD82Nzc3NTXx+fzGxkbcTiMQCBobG3k8XkNDA35eV1eHkx7p+Tg7j8lkttta2W7Tpr6+fq+fOyCRWCwePXr069ev4+LiPD09yQ4HKBQkQEDRmpubZbKi+vr6+vp64vuMy+Xi3+v4Sw4/l5l2rXsoFApOifT19YkkSU9Pj0ql6urq4kcajcZkMvEjQoh4ghDC6yOEcKbVpaX9iEQi4XK5+HlDQ4NIJEII4XSk80vxk7a2tvr6evwoFot5PB5+FIlEDQ0Nra2tfD4fP/Y8bDqdjpNjnCXj5Bg/x8mxtra2np6e/MVZ6GozkB0/fvyDDz5wdXVNTEwkOxagaJAAgX4Dp0f4isbbLnM0NTXhxgOiOUEoFBLfsq2trcR3toIR6ZE8TU1NLS2tt23Y8VJ5OMnoxlIiiVEwnC/iRwaDgfMYOp2uo6Ojrq6OH9u9PKqnp4cb/3AhdC8DXdXa2jpq1Ki8vDxo/hmY4C4w0G/gH/FGRkY93xVOg4gGCbFYTDRUtLW1cblc3ICBfx7gzAn9dZkG/b0VRH4pQohoryKWIqn2EiWkpaWlqamJG65w/oEQwvkHXkFPTy81NbWxsXH06NFsNltmKW5XQ1KtXziPITIb3JELH0VLS4t4Ar27AIkOHz6cl5fn4eEB2c/ABC1AACiFjtOjziRPHecTPc82wsPDV65c2dTUtHr16iNHjuAkCYB+isvlDh8+vKamJiEhwdXVlexwAAkgAQIAdFZycnJAQEBlZeWMGTMuXboEPYJB//XFF1/8/PPPS5cuPX/+PNmxAHJAAgQA6IK8vDw/P7/s7Oxhw4Zdv359xIgRZEcEQJcVFBSMHj26ra0tKyvL1taW7HAAOWDOWwBAF9jY2CQmJk6fPj03N9fV1fX+/ftkRwRAl23ZsqW5ufmTTz6B7GcggxYgAECXiUSijRs3HjlyREND47fffnvvvffIjgiAzoqLi/Py8jI2Nn716hVcxh3IoAUIANBlNBrt8OHD+/bta21tDQoK2rlzJ/yUAv2CUCgMDg5GCP3888+Q/QxwkAABALopJCTk4sWLWlpau3btWrZsWTfG+AZAwf7v//7v5cuXbm5ugYGBZMcCSAaXwAAAPfLw4cN58+ZxOBxXV9crV670ykBNAPSF/Px8e3t7oVCYlpbm6OhIdjiAZNACBADokYkTJz5+/NjJySkpKWny5MkvX74kOyIA2iGRSD788EOBQBAcHAzZD0DQAgQA6BUNDQ3Lli27fv26gYHBn3/+OW3aNLIjAuBvfv31148++sja2vr58+fEIOZgIIMECADQO8RicUhIyMGDB+l0+q+//rpy5UqyIwLgv8rKythsNo/Hu3379tSpU8kOBygFuAQGAOgdVCr1wIED+/btE4lEq1atCgkJaWtrIzsoAP578YvL5a5btw6yH0CAFiAAQC+7cePGkiVLGhoaFi5c+Pvvv3dpNnsAet2RI0c+/vhjKyur58+fM5lMssMBygISIABA70tPT/f39y8uLp40adKVK1dYLBbZEYEBKjc3d+zYsQKBID4+fvr06WSHA5QIXAIDAPS+MWPGpKSkjBs3LiUlZfLkyVlZWWRHBAYikUj03nvv8fn8L774ArIfIAMSIABAnzAzM7t3756/v39+fv6kSZOio6PJjggMODt27Hj48OG4ceN27dpFdixA6cAlMABAHxKLxZs3b96/fz+NRgsNDV2/fj3ZEYGB4s6dO7NmzVJXV3/06JG9vT3Z4QClAy1AAIA+RKVSQ0NDjx49ihDasGED3BoGFKOiomL58uVisfjAgQOQ/YB2QQsQAEARYmNjFy9ezOPx5s+fHxYWpq2tTXZEQGWJRKIZM2Y8ePBg2bJl586dIzscoKQgAQIAKEhGRoafn19RUZGTk1NkZKSFhQXZEQHVtHXr1j179rDZ7IcPHzIYDLLDAUoKEiAAgOKUl5fPnTv38ePH5ubmkZGRY8eOJTsioGquXLmyYMECbW3t1NRUOzs7ssMBygv6AAEAFMfU1PTu3bsBAQGlpaVTpkyJiooiOyKgUl68eBEYGIgQ+u233yD7AR2DBAgAoFAMBuPPP//csWMHn88PCAj4z3/+Q3ZEQEXU1dXNnz8fj/qzdOlSssMByg4ugQEAyHH8+PGPP/64tbV17dq1Bw8epNFoZEcE+jGxWDxnzpzY2Fhvb+/r169TqVSyIwLKDhIgAABp4uPjFy1axOVyvb29L168qKurS3ZEoL8KCQnZv3+/lZXVo0ePBg8eTHY4oB+ABAgAQKbMzEw/P7+CggJHR8fIyEhLS0uyIwL9z88///zFF18YGBgkJCRA1x/QSZAAAQBIVl1dPX/+/ISEBFNT08jISGdnZ7IjAv1JeHj40qVLaTRadHT0zJkzyQ4H9BvQCRoAQLLBgwfHx8cvX768vLx86tSpV69eJTsi0G88ePAgKChIIpH89ttvkP2ALqHu3LmT7BgAAAMdjUZbsGABQiguLu6PP/6QSCTTpk0jOyig7LKzs728vBoaGvbs2QPTzIGugktgAAAlcurUqY8++kgoFH7wwQeHDh1SV1cnOyKgpKqqqlxdXXNzc9euXYsnmwOgSyABAgAol8TExICAgOrq6lmzZoWHh+vp6ZEdEVA6AoFg5syZKSkpc+bMuXLlCoyhALoBEiAAgNLJzc318/N79eqVvb19ZGSklZUV2REBJdLa2hoQEBAdHT1hwoQ7d+7AxLqge6ATNABA6QwbNiwpKWnq1KkvXrxwcXFJTEwkOyKgLEQi0fLly6Ojo21sbCIjIyH7Ad0GCRAAQBkZGhrGxsYGBgZWV1d7enqeP3+e7IgA+cRi8apVqyIiIiwsLOLi4oyNjcmOCPRjcBcYAEBJ0Wi0gIAALS2t2NjYP//8E24NG+AkEsm6detOnTrFYrFu3749fPhwsiMC/Rv0AQIAKLuLFy+uWrWqqalpzZo1R44cgVvDBiCJRLJ+/fojR44YGxvfuXMHhnsGPQcJEACgH0hOTg4ICKisrJw5c2ZERIS+vj7ZEQHFkUgkn3zyyaFDhwYPHnznzh17e3uyIwKqAPoAAQD6gcmTJycnJ48ePfrWrVsuLi45OTlkRwQUZ8uWLYcOHdLX14+JiYHsB/QWSIAAAP2DjY1NYmLi9OnTc3NzXV1d79+/T3ZEQBG++uqrn3/+WU9P7+bNmzBPHOhFkAABAPoNAwODmzdvrlu3rqamxsvLKywsjOyIQB+SSCSff/75Dz/8wGAwIiMjXVxcyI4IqBS4CwwA0J+oqan5+fkZGBjcuHHj8uXLEolk6tSpFAqF7LhALxOLxXg6FD09vRs3bri7u5MdEVA10AkaANAvXbp0KSgoSCAQLFmy5NSpU5qammRHBHqNUCh87733wsPDDQ0Nb9y4MWHCBLIjAioIEiAAQH/18OHDefPmcTgcV1fXK1euGBkZkR0R6AUCgeAf//hHTEyMqanpzZs3odcz6COQAAEA+rGSkhJ/f/9nz57Z2tpGRUWNGjWK7IhAj3C5XD8/v8TERGtr67i4OFtbW7IjAioLOkEDAPoxCwuL+/fvz5kz582bN25ubnfv3iU7ItB9lZWV06dPT0xMZLPZCQkJkP2APgUJEACgf2MymVevXt2wYUNtba23t/fp06fJjgh0R1FRkYeHx7Nnz1xcXO7du2dmZkZ2REDFQQIEAOj3qFTqgQMH9u3bJxKJVq1aFRIS0tbWRnZQoAsyMjLc3d1zcnKmT59+69atQYMGkR0RUH3QBwgAoDquXLmyYsUKgUCwcOHC33//XUtLi+yIwLtdv3596dKlfD5/8eLFv//+u4aGBtkRgQEBWoAAAKojICAgKSlpyJAhERERM2bM4HA4ZEcE3uHYsWMBAQF8Pj84OPj8+fOQ/QCFgQQIAKBSxowZk5KSMm7cuJSUlMmTJ2dlZZEdEWifRCLZuXPn2rVrEUKHDh0KDQ1VU4OvJKA4cAkMAKCC+Hz+8uXLIyMjmUzmH3/8MXv27LetyeFw+Hw+3HCkYM3NzatWrfrjjz+YTOaFCxd8fX3JjggMOJAAAQBUk1gs3rx58/79+2k0Wmho6Pr16+XXaWpqmjlzpqGhYVRUlOIjHLAqKirmzp376NEjc3PzqKgoJycnsiMCAxG0NwIAVBOVSg0NDT169ChCaMOGDfK3hrW1tQUFBSUnJ1+/fv3GjRskhalq+Hx+xytkZGRMmjTp0aNH48ePf/ToEWQ/gCyQAAEAVNnatWujoqJ0dXX379+/cOFCgUBALPr6668jIiLw802bNgmFQpJiVB1PnjyZMGFCQ0PD21a4du2au7t7YWGhj4/P7du3TU1NFRkeANIgAQIAqDhvb++EhARLS8vLly+7ubmVlJQghE6ePPnjjz8S08i/evXqP//5D6lh9nsCgWDFihXZ2dkhISHyS8Vi8VdffRUQENDQ0LB169aoqCgmk6n4IAEgQB8gAMCAUF5ePnfu3MePH5ubm+/atevjjz9ubW2VXkFHRycnJwfaJLpt48aNBw4cwM/Pnj27fPlyYhGPxwsKCrp69aqOjs6JEycWLVpEUowA/A8kQACAgaKxsXHFihVXr15VU1Nrd6jo1atXnzhxQvGBqYCbN2/6+PgghPB3iq6ubnp6upWVFULo1atXAQEBL1++tLW1vXz5soODA7mhAoBBAgQAGECqqqpGjhxZV1cnv4hCoVAolOTk5AkTJig+sH6Ny+U6OjoWFxdLF06YMCEhISEmJiYwMLC+vn7atGkXL140MjIiK0gAZEAfIADAQCEUChcvXtxu9oMQkkgkbW1tISEh8LOwq9atWyeT/SCEUlNTPT09582bx+Pxvvjii/j4eMh+gFKBFiAAwIAgkUiCgoLCwsLeuebvv/8eGBiogJBUQ1hYWLvVhTuYq6urHzt2LCgoSOFxAfAOkAABAAaEHTt2fPvttxTKuz/0WCxWTk6Orq6uYgLr10pLS+3t7blc7ttWYLFYmZmZMLs7UEJwCQwAoPrCwsK+++479FcX3Y5xOJwffvih74Pq9/BIkh1kPwghDofzwQcfKCwkADoPWoAAACquoKDAwcHhnSMUS6PT6ZmZmcOGDeu7qFTAL7/88tlnn3VmzcOHD69bt66v4wGgSyABAgCovtra2oiIiAMHDmRkZHRyk7lz5169erVPo+rXsrKynJ2dm5ubO7OyhoZGamqqo6NjX0cFQOdBAgQAGEDS0tJ+/fXXc+fOdaZBKDo6uoNp5AcyoVA4adKkp0+fdn4TBweH1NRUTU3NvosKgC6BPkAAgAHE2dn56NGjpaWlp0+f9vT07HjlzZs3y4wWDbCdO3d2KftBCGVkZOBuWAAoCWgBAgAMXNnZ2adPnz5x4kRVVVW7K/zyyy+bNm1ScFRKLikpacqUKW1tbZ38+rCzs/P39/f09Jw2bRqNRuvr8ADoJEiAAAADXUtLy5UrV06cOBEfHy8zRYauru7r16+NjY07sxOBQNDQ0CAUCuvr65ubm5uamng8nlAo5PF4TU1N0t1l+Hy+dNtSY2Oj/Fz02traGhoaxJ9aWlrS1490dHTodLq+vr6Ghoa2tjaTyaTT6Xp6epqamlpaWrq6ulQqtav10BmNjY1jx459/fp1x6tpamq6u7v7+fnNnz/f0tKyLyIBoIcgAQIAgP8qLCw8efLkyZMni4qKiMK5c+d+9NFHdXV1tbW1Mo/4SXNzc8e3gpOCSqXq6uoyGAxDQ0MDAwMDAwP8RP5x8ODBenp6ndztBx98cPz48bctNTIy8vHx8ff39/HxgcnegZKDBAgAMEC1tLRUVlaWlJRUVlaWlpZyOBz8WFZWVlBQUFdX18mPR3V1dR0dHQaDQafTDQwM6HQ6g8GQb6Ehrv7ItOXgNpt29yzTUCT9Z7stTPX19UKhsKGhQSAQtLS0vG3SD3kaGhomJiZmZmYsFsvCwsLY2Njc3JzFYpmZmZmYmBgbG+P2pGvXrs2bN09mWwqFMm7cOD8/Pz8/P2dnZzwANADKDxIgAICKEwgE+fn5+fn5BQUFxGNxcXF1dXUHWxkaGg4ePFgoFNbU1GhpaS1dulSm4YR4QqfTFXYuXSUQCGTaq+T/rKmpKS8vFwgEb9sJlUrFKdGLFy+IC3mampozZ86cN2/enDlzzMzMFHVCAPQaSIAAAKqjtrY2Ozs7KysrLy8PJzr5+fmVlZXya6qpqbFYLBaLZW5ubmxsLN/sId1Ik5KSMnToUFNTUwWeiqLx+XzpxrCysrKKiory8vLy8vKKioqamhr5TahUqrm5ubW1tZWVlbW19YgRI0aNGjVq1Ki3NWgBoFQgAQIA9Fd1dXWZmZk43cFP8vPz5T/TDAwMbORYWlrCHUmdl5CQcPny5eHDh6upqeXl5ZWVlZWXl+fl5bVb4aampmw2287Ojs1m29jYODo6dqYXOQAKBgkQAKB/EAqFGRkZT548SUtLe/78eXZ2tkzXYwqFYmVlNWrUKDabPWzYMNwsMXToUOl7qUDvamhowFcV8/LyXr169fLly8zMTPkxBczMzEaPHu3k5OTs7Ozs7Dxs2DA1NRiFDpAMEiAAgJJqbW3NyclJkyJ9JzmNRrO0tLSxscEtDXZ2dk5OTjo6OiQGDDAul/vmzRvcJocfCwoKpMcXYDKZjo6Ozn8ZPXo05ENA8SABAgAokdLS0vv37z948CAlJeXFixfS90AxmcyxY8eOGzfO2dnZyclpxIgRytz7GEhraGjIzs7GrXdpaWnyr6yTk5Orq+uUKVPc3Nw6f08+AD0BCRAAgGRlZWWJiYnx8fEJCQlZWVlEObQTqCqRSPTq1SuiYe/JkydNTU14EZVKHTlypLu7u5ub24wZMywsLMgNFagwSIAAACQoLi6Ojo6+e/fu/fv3y8rKiHJra+spU6bgloARI0bAoDIDgVAoTE9PT0hIuH//fkJCgvTwBKNHj/bw8Jg5c6aXl5e+vj6JQQLVAwkQAEBBxGLxs2fPIiMjo6Kinjx5Qnz42NjYuLm5ubu7z5o1y9ramtwgAeny8vISEhISExNv3rxZUFCAC6lUqpOTk5+fn7+//7hx4yAzBj0HCRAAoG9VV1fHxMRcv3795s2btbW1uHDIkCGzZ8/29vb28PAwMjIiN0KgtN68eXPv3r2YmJibN2/W19fjQisrK19f3zlz5kyfPh3GHALdBgkQAKBPNDU1RUVF/f777zExMSKRCMGPeNAD7TYfampqenp6BgUFBQQEqKurkx0j6GcgAQIA9CahUBgTE3P27NnIyEjcs9XQ0BD/Xvfy8jI0NCQ7QNDvFRQUREdHX79+/datWy0tLQghY2PjxYsXr1ixYtKkSWRHB/oNSIAAAL3j4cOHp06dCg8Px9MmaGlp+fv7r1ixwsfHB+5XB32hrq7u0qVLYWFhDx48wOMM2drarlixYtWqVdCZDLwTJEAAgB4Ri8XR0dF79uxJTExECKmpqU2ePDkoKGjp0qW6urpkRwcGhMz5Nq8AABTgSURBVNLS0oiIiPDwcOKf0NfXNyQkxNPTk+zQgPKCBAgA0E01NTXHjh07ePBgSUkJQsjGxmb9+vXLli2DucEBWTIyMk6fPn38+HE8Tcq4ceM+/fTTJUuWQBskkAcJEACgy/Ly8vbs2RMWFiYQCBBC06ZNCwkJmTt3LgxUCJQBn8///fffQ0NDc3JyEEKmpqbr168PDg6GJkkgDRIgAEAXlJaW/vTTT0ePHm1paaHT6fPmzfvss88mTpxIdlwAyJJIJLdu3QoNDb1+/bpEIjE0NNy4ceNnn33GZDLJDg0oBUiAAACdIhAI9uzZ8/PPPzc1NTEYjE8++WTTpk0sFovsuAB4hxcvXnz33XcRERFtbW1mZmZ79uxZsWIFjMIAIAECALzblStXPv3008LCQjqdvm7dui+//NLExITsoADogvT09G3btkVFRSGE3NzcDh065OjoSHZQgExwwR4A0BEul/vRRx/Nnz+/sLBwxowZT548CQ0NhewHyPj8888pUl68eEF2RLLGjBkTGRl569YtNpudmJjo7Oy8detWoVBIdlyANNACBAB4q3v37i1btqy8vNzY2PjAgQOLFi0iOyKg1O7evTt9+nSEUEZGhr29PdnhtK+1tfWnn3769ttvhUKhi4tLeHj40KFDyQ4KkABagAAA7du7d6+np2d5efmiRYsyMzMh+wGqQV1d/euvv378+LGTk9OjR4/Gjx8fHx9PdlCABJAAAQBkSSSSzZs3b968GSH0448/Xrx4cfDgwWQHRbKKigp8cScsLIzsWEAvcHBwSElJ+fDDD6urq318fE6dOkV2REDRaGQHAABQLhKJZM2aNadOnTIwMLh69aqHhwfZEQHQJzQ0NH799VdHR8eQkJD3339fKBSuXbuW7KCA4kALEADgb7788stTp06ZmJjcuXMHsh+g8j755JMzZ85QqdT169dfu3aN7HCA4kACBAD4nzNnzuzZs4fJZMbExIwZM4bscABQhOXLlx89erStrW3ZsmWZmZlkhwMUBBIgAMB/FRcXb9y4kUKhnD9/nsTsJyUlRfqG6q1bt5aUlAQHB48aNUpHR0dHR8fZ2Xnfvn2tra1v20NJSck///lPR0dHXV1dLS0tW1vbNWvWPHnyRH7N27dvL1y40MLCgk6nMxiMkSNHBgYGnjt3Dk9oLx2Pqakp/jMwMJCIrXtdozp/0G5XQmFh4aeffmpnZ6ejo6OtrT1s2LDVq1enpaVJr9Pc3Cx9CB8fH7FYvH//fmdnZyaTyWAwnJyc9u/f/7Y7haOiory9vQcNGqSpqWllZfX++++/fv265ydOltWrV2/ZskUgEKxatUokEpEdDlAICQAASCQSiWTJkiUIoY8//pjsQP7L3NwcITRr1iwjI6MffvihtLS0qakpJSXF1dUVIeTm5lZfXy+/1YULF7S0tAwNDU+ePFlTU8Pn8+Pi4thsNoVC+eqrr9ra2og1v/32W4TQ7NmzHz16xOfzKysrL126ZG1tjRDy9vaW2W15eTn+zDxz5kxPTqpLB+1eJYSFhWloaAwePPjcuXM8Ho/H4126dMnCwgIhtH37dvlDzJkzByHk6enp5+f3/fffczicpqam6OhoIyMjhFBwcLD8JsHBwQghU1PT8PBwHo/H5XLPnz9vZWW1a9cuXEsZGRk9PHHFEwqFY8eORQgdPHiQ7FiAIkACBACQSCSSnJwcKpVqYGDQblZBCvzdL/+F1NDQYGtrixBatGiRzCaxsbFqamoUCiUxMVG6vLKyEs/asXv3blySnZ1NoVA0NTUFAoH0mtnZ2XQ6vY8SoK4eVNL1SoiJiVFTU1NTU0tOTpYuz8rKwjOiHzt2TOYQOAFSU1Pbtm2bdPmRI0cQQlQqtbS0VLr84MGDuPzp06fS5WlpaTQard0EqBsnToqEhASE0NChQ4VCIdmxgD4HCRAAQCKRSL766iuE0FdffUV2IP+Dv/tZLJZIJJJZhL+DEULSX/Otra24RWHOnDnye/v+++8RQjQaLT8/XyKRHD58GCFkYmIiv+akSZP6KAHq6kElXayEjmtg6dKlCCFjY+OmpibpciIBqqyslC7Pzc3F+79w4QJRKBAIDA0NEUILFiyQP8S8efPaTYC6ceJkmTJlCkIoOjqa7EBAn4M+QAAAhBC6ceMGQmjx4sVkByLLw8ODSqXKFHp7e+MnZ86cIQrj4uLy8/MRQr6+vvL7wYUikej48eMIIQaDgRCqqKg4ceKEzJrJyckxMTG9dgJSun3QTlYCUQN+fn7yO5kxYwZCqLKy8s6dO/JLhw8fjq95ESwtLfETIvlDCN24caO2thYh5OXlJb+TmTNnths/KbXdPfgtgN8OQLVBAgQAQG1tbRkZGXp6eko4PSTuvCLD2toaX21JTU0lCvH1C4TQqFGj5DchCh88eIAQ8vb21tfXRwi9//7706ZNO3PmTF1dXW/HLqvbB+1kJRA10O40FMSEDw8fPpRfSvTyJqirq6urqyOEmpubiULiWCNHjmw3pHbjJ6W2uwe3AKWnp5MdCOhzkAABAFBNTY1IJDIzM6NQKGTHIktDQ0O+UE1NjclkIoQ4HA5RSDRUGBgYtLsfLS0thFBZWRlCyNjYODo62sbGBiF07969oKAgIyOjqVOnHjhwoL6+vg/OA/XkoJ2sBKIGPDw8KHKIFiPpFp13HkKmhDiWnp6e/Pq4pUceKbXdPfiaY7tVBFQMJEAAAIQbEpTz7l9Jpyds7uSaRJI3efLknJycK1euLFiwQEtLSywW379/f+PGjTY2NpGRkd0M9126d9BOnhqxmkz3ZBlHjx7tdvzEIdrNlTtIoEmp7W7A7wLc9AVUGyRAAACkr6+voaFRVlYmFovJjkWWQCCQL2xra2toaEAImZiYEIVmZmb4SbuXV1paWpqamtDfr/VQqdR58+ZdunSpqqrqwoULeCbz2tra5cuXFxcX9+p5/E83DtrJSiBqgMvl9knoUsdqt5JxDb8NKbXdVUVFRejv/1dAVUECBABAFArFxcWlsbGx3d4h5CooKJAvzMvLw7/UJ0yYQBQSE3dkZ2fLb/Ly5UuZ1aQxGIwlS5bcvn37xx9/RAjx+XwFzBDe+YN2shKIU8vKyuqLgBFCEydOxE9ycnLkl5aWlnZmJ6TUdifdunUL/f3/CqgqSIAAAAj9dZOUEk51npCQIBQKZQpjY2Pxk8DAQKLQ09MT9zKJjo6W3w8upNFoa9asQQjt3r2bTqe3tLTIrPbxxx/jJzKL5HvDdENXD0roZCUQNfC2i0ru7u5qamrtJoidNHv27EGDBkkfXdq9e/fa3arbJ65gEonk3Llz6C03EgIVAwkQAAAhhFavXq2pqXny5EncR1h51NfX79+/X7qEz+fv3bsXIbR48WKiQQIhRKPRDh8+rKamduPGjaSkJOlNqqqqQkNDEULbt28n7lRqbW29fPmyzOFwQxGNRvP09JQu19fXxz2l+Hw+Lrl16xaFQvn666+7dDpdOiihk5VA1EBMTExcXJzMTq5cuZKYmLhw4cLRo0d3KWZpmpqaeESla9euPXv2THrRq1evwsPD37Zh905cwa5evfr8+XNHR0d3d3eyYwF9TwFjDQEA+oWNGzcihObOnUt2IP+F78d5//337ezsdu/eXVpa2tzc/PDhQzc3N/SuqTAGDRp06tSp2traxsbGuLg4fFv4F198QUyF8d133yGEmEzm//3f/71586a5uZnD4Vy6dMnGxkZNTa3dyRBwq8CsWbOqq6s5HI6fn5+6unpaWlrnz6gbB+1GJeAa0NHR2bt3L14/Jydnx44dmpqanp6efD5fZn08EGK7oxHiW8N++OEHmfJNmzYhhExNTSMiIhoaGng8XkRExJAhQ4hxpGQGQuzGiStefX09Hing4sWLZMcCFAESIADAf3G53CFDhiClmQsJf/dv2bKFx+Nt27bN3t6ewWBoa2uPGzdu7969HUxWUFxc/Nlnn7HZbAaDoampaW1tvXLlytTUVOl1BALB9evXN2zYMGHChCFDhuCb5EeMGPHBBx+87RaqysrKVatWmZqa0mi0QYMGTZ06NSYmpktn1I2Ddq8S8HSwDg4OxGSo/v7+ly9flhlOmphng8BmsyUSSbsdnFesWCG9bVRUlJeXl4GBAZ1ONzMzW7BgwYMHD2SGWNyxY0e3T1zB2trali9fjt4yiDZQSRRJp28xBQCovFu3bvn4+CCErl69Sno3CAsLi9LS0i1btuCusgMTVIJibN++fffu3SYmJmlpacTNdEC1QR8gAMD/zJw588iRIyKRaMGCBfI9NgBQSXv27Nm9e7eWltbly5ch+xk4IAECAPzN+++///3337e0tCxZsuTUqVNkhwNAHxKJRBs3bty6dSuDwbh69eqkSZPIjggoDiRAAABZX3755YEDB8Ri8erVqzds2CB/AzYAKqCiomLmzJkHDhzQ19ePjY2dNWsW2REBhYIECADQjg0bNly7dk1fX//QoUOTJ09+8eIF2REB0JuuXbvm5OR0//59Ozu7lJQUfFcdGFAgAQIAtG/OnDmPHz92dnZ+8uTJ+PHjv//+e8U0BaWkpODJO/Gwwnv27KFQKO1Ob64k5KcdbdfOnTs7v89+Vwn9SHV19apVq+bNm8fhcJYtW/bw4cN2Z7YHKg/uAgMAdEQkEv373//+5ptvhELhsGHD9u7d6+fnR3ZQAHSHSCQ6ceLE119/XV1dbWBgEBoaKj2SOBhoIAECALzb8+fPN27ceP/+fYTQrFmzdu/eDZMlgX6kra0tIiJix44dL1++VFNTW7ly5Y8//mhsbEx2XIBMkAABADpFIpFcuHBhy5YteOLuuXPnfvnll3DXDFByYrH48uXLu3fvTk9PRwhNmjQpNDQU0neAIAECAHSJUCg8derUzp07y8vLEULOzs7BwcHLli1TV1cnOzQA/obH4508eXLfvn0FBQUIITabvWPHjoULF1IoFLJDA0oBEiAAQJc1NjYeO3bsP//5T15eHkJoyJAhGzZs+PDDDw0NDckODQCUm5u7f//+U6dONTQ0IITc3Nw2bdo0f/58NTW47wf8DyRAAIBuamtru3r1amho6L179xBCDAZj+fLl7733nru7O3zTAMVrbm6Oioo6ffp0dHR0W1uburr6okWLPv30UxcXF7JDA8oIEiAAQE89e/bs8OHDZ86caWpqQghZWFgsWLAgKCjI2dmZ7NCA6mtra0tKSgoPDz979mxNTQ1CSE9Pb+XKlZ9//jme3BeAdkECBADoHRwO5/z582fPnn38+DEucXBwWLFixbJlyywtLcmNDaikJ0+enD179sKFC2VlZQghGo3m5eW1fPnyBQsWaGlpkR0dUHaQAAEAellhYeGFCxdOnDiRk5ODS+zs7Pz9/f38/FxdXeHqGOiJ5ubmhISE+Pj4q1evvnz5Ehc6OzsHBgYuXbqUxWKRGx7oRyABAgD0lZSUlHPnzoWHh1dUVOASY2NjX19fX19fLy8vPT09csMD/UhJScn169ejo6Pj4+MFAgEutLOzW7Zs2fLly21sbMgND/RHkAABAPqWRCJJS0u7fv369evX09LS2traEELq6uru7u6zZs2aMmWKi4sLnU4nO0ygdHg8XkJCwv3792NiYvAoPgghTU3NqVOnzpkzZ86cOZD3gJ6ABAgAoDjV1dV37tyJjIyMjIzkcrm4UF1d3dHR0dPT083NbcqUKdAyNJBVVVWlpKQkJiYmJCSkpqa2trbicmNjY29vb39/f29vb11dXXKDBKoBEiAAAAlaW1sTEhLu3r17//79hw8f4tvHEELq6urOzs4eHh7u7u7Ozs7m5ubkxgn6Wltb2+vXrx89evTgwYMHDx5kZ2cTiwYNGuTm5jZ16tQZM2Y4OTmRGCRQSZAAAQBIJhKJ0tPTExISEhMTb926VVtbSyzS19dns9nOf2Gz2STGCXpLWVlZ2l+SkpKkX3ETE5Px48e7u7t7enqOHTsWusyDvgMJEABAiYjF4ufPn+NmobS0tNevX0t/RhkbG+NMyMnJic1m29rawhQcyq+xsfHly5eZmZlPnz5NS0t7+vQpn88nlqqrqzs4ODg7O7u6unp4eNja2pIYKhhQIAECACivhoaG9PR0orXg5cuXuA81RqPRLC0t7ezs2Gy2jY2NnZ2dk5OTjo4OiQEDLpf75s2bzMzMrKysvLy8zMxM+VdtxIgRzlJgzB5ACkiAAAD9RkNDA25FyMjIwN+sPB5PegU1NTUrK6vRo0cPHz7c+i9WVlaQFfWF6urq/Pz8goKC/Pz8/Pz8V69eZWVlcTgcmdVYLJadnd3o0aOdnJycnZ0dHByg3Q4oA0iAAAD9WElJSfZf8HWWqqoq+dWMjIyIZAg/sbS0NDc3ZzKZCg+5/6mqqqqoqMiXgpMePNWoNAqFMnTo0FGjRrHZ7FGjRuG8x8DAgJSwAegYJEAAAJVSU1OTnZ395s0b6S/ssrIy6aswBG1tbTMzMxMTE1NTU1NTUxMTEzMzMxaLZW5ubmxsPECGFW5tba2srCwrK6uoqKioqCgrK+NwOKWlpZWVlSUlJZWVlUKhUH4rOp0+dOhQIqG0trYePnz4yJEjGQyG4k8BgG6ABAgAoPqEQmFhYaH09Zri4mL8ld/c3Py2rdTU1AwMDAwNDdt9JJ5oa2vr6enR6XTlaU+qq6sTCoWNjY08Hq+urq62tpZ4xE+kS+rr6zvYlY6OjoWFBYvFks51rKyszM3N4RYt0K9BAgQAGNDq6urKy8ulWz6Ix8rKSuk7tDtDW1tbQ0NDX19fQ0ODwWDo6OjQ6XT8p7a2tsxqxJ9MJpNGo8nsqrm5mRgeCSHU2Ngo3RLD4/GEQiGPx2tqampubuZyuUKhkM/ny6z2TjQazdDQkMViva0lTDpsAFQJJEAAANARmfYS+cfGxkY+ny8UCrlcbktLCzFTFSl0dXXpdLqurq6Wlpampqaenh7RUjVo0CD5pizlabUCQMEgAQIAgF7W0NAgFArr6+txKw6PxxOLxfJNOGKxmPiTy+VKfxrTaDTp1KTdBiRcyGQy6XS6np6epqYm3E8OQOdBAgQAAACAAQe6sAEAAABgwIEECAAAAAADDiRAAAAAABhwIAECAAAAwIADCRAAAAAABpz/B4OhBdOPv25LAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "causal_graph = \"\"\"digraph {\n", "treatment[label=\"Program Signup in month i\"];\n", "pre_spends;\n", "post_spends;\n", "Z->treatment;\n", "pre_spends -> treatment;\n", "treatment->post_spends;\n", "signup_month->post_spends;\n", "signup_month->treatment;\n", "}\"\"\"\n", "\n", "# Post-process the data based on the graph and the month of the treatment (signup)\n", "# For each customer, determine their average monthly spend before and after month i\n", "df_i_signupmonth = (\n", " df[df.signup_month.isin([0, i])]\n", " .groupby([\"user_id\", \"signup_month\", \"treatment\"])\n", " .apply(\n", " lambda x: pd.Series(\n", " {\n", " \"pre_spends\": x.loc[x.month < i, \"spend\"].mean(),\n", " \"post_spends\": x.loc[x.month > i, \"spend\"].mean(),\n", " }\n", " )\n", " )\n", " .reset_index()\n", ")\n", "print(df_i_signupmonth)\n", "model = dowhy.CausalModel(data=df_i_signupmonth,\n", " graph=causal_graph.replace(\"\\n\", \" \"),\n", " treatment=\"treatment\",\n", " outcome=\"post_spends\")\n", "model.view_model()\n", "from IPython.display import Image, display\n", "display(Image(filename=\"causal_model.png\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More generally, we can include any activity data for the customer in the above graph. All prior- and post-activity data will occupy the same place (and have the same edges) as the Amount spent node (prior and post respectively). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## II. Identifying the causal effect\n", "For the sake of this example, let us assume that unobserved confounding does not play a big part. " ] }, { "cell_type": "code", "execution_count": 4, "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(post_spends|signup_month))\n", "d[treatment] \n", "Estimand assumption 1, Unconfoundedness: If U→{treatment} and U→post_spends then P(post_spends|treatment,signup_month,U) = P(post_spends|treatment,signup_month)\n", "\n", "### Estimand : 2\n", "Estimand name: iv\n", "Estimand expression:\n", "Expectation(Derivative(post_spends, [Z, pre_spends])*Derivative([treatment], [\n", "Z, pre_spends])**(-1))\n", "Estimand assumption 1, As-if-random: If U→→post_spends then ¬(U →→{Z,pre_spends})\n", "Estimand assumption 2, Exclusion: If we remove {Z,pre_spends}→{treatment}, then ¬({Z,pre_spends}→post_spends)\n", "\n", "### Estimand : 3\n", "Estimand name: frontdoor\n", "No such variable(s) found!\n", "\n" ] } ], "source": [ "identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)\n", "print(identified_estimand)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the graph, DoWhy determines that the signup month and amount spent in the pre-treatment months (`signup_month`, `pre_spend`) needs to be conditioned on." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## III. Estimating the effect\n", "We now estimate the effect based on the backdoor estimand, setting the target units to \"att\"." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "## Identified estimand\n", "Estimand type: nonparametric-ate\n", "\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "────────────(Expectation(post_spends|signup_month))\n", "d[treatment] \n", "Estimand assumption 1, Unconfoundedness: If U→{treatment} and U→post_spends then P(post_spends|treatment,signup_month,U) = P(post_spends|treatment,signup_month)\n", "\n", "## Realized estimand\n", "b: post_spends~treatment+signup_month\n", "Target units: att\n", "\n", "## Estimate\n", "Mean value: 97.57746107152285\n", "\n" ] } ], "source": [ "estimate = model.estimate_effect(identified_estimand, \n", " method_name=\"backdoor.propensity_score_matching\",\n", " target_units=\"att\")\n", "print(estimate)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The analysis tells us the Average Treatment Effect on the Treated (ATT). That is, the average effect on total spend for the customers that signed up for the Rewards Program in month `i=3` (compared to the case where they had not signed up). We can similarly calculate the effects for customers who signed up in any other month by changing the value of `i`(line 2 above) and then rerunning the analysis. \n", "\n", "Note that the estimation suffers from left and right-censoring. \n", "1. **Left-censoring**: If a customer signs up in the first month, we do not have enough transaction history to match them to similar customers who did not sign up (and thus apply the backdoor identified estimand). \n", "2. **Right-censoring**: If a customer signs up in the last month, we do not enough *future* (post-treatment) transactions to estimate the outcome after signup. \n", "\n", "Thus, even if the effect of signup was the same across all months, the *estimated effects* may be different by month of signup, due to lack of data (and thus high variance in estimated pre-treatment or post-treatment transactions activity)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## IV. Refuting the estimate\n", "We refute the estimate using the placebo treatment refuter. This refuter substitutes the treatment by an independent random variable and checks whether our estimate now goes to zero (it should!)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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", "/home/amit/py-envs/env3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:993: 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": [ "Refute: Use a Placebo Treatment\n", "Estimated effect:97.57746107152285\n", "New effect:1.251821060965955\n", "p value:0.430226053357455\n", "\n" ] } ], "source": [ "refutation = model.refute_estimate(identified_estimand, estimate, method_name=\"placebo_treatment_refuter\",\n", " placebo_type=\"permute\", num_simulations=20)\n", "print(refutation)" ] } ], "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 }