{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Causal Analysis of *How is the implementation of existing strategies affecting the rates of COVID-19 infection?*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*(Better displayed in [nbviewer](https://nbviewer.jupyter.org/) as red warnings in font tag may not be displayed on github)*\n",
"\n",
"We believe that a question such as *How is the implementation of existing strategies affecting the rates of COVID-19 infection?* requires a proper causal analysis of the data.\n",
"\n",
"Such question requires considering what the *causal effect* of certain strategies may be, and an evaluation of *counterfactual effects*.\n",
"\n",
"To answer this question we rely on the formalism of [*causal models*](http://bayes.cs.ucla.edu/BOOK-2K/) and the [*dowhy*](https://github.com/microsoft/dowhy) library.\n",
"\n",
"\n",
"### Disclaimer\n",
"In this notebook we setup a **preliminary template** for a causal analysis of this question. In particular, we setup a simple model that tries to estimate the **average causal effect** of a chosen policy.\n",
"\n",
"**This is a work in progress.** Publishing this notebook has the sole aim to share an initial analysis and draws suggestions and critiques for improvement.\n",
"\n",
"In no way, at the moment, these results should be taken to be significant. The data considered are limited, the model simplistic and debatable. **NO conclusions on real world policies should be drawn from this notebook**. \n",
"\n",
"Limitations will be highlighted in the notebook in RED (coloured fonts may not be displayed on github)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Setup and testing of dowhy\n",
"\n",
"The aim of this notebook is to (ii) gather some data from relevant sources, (ii) set up a mock causal model, and (iii) run a standard causal analysis (identification, estimation, refutation) using *dowhy*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Data Part\n",
"In this part we collect the data that will be used in the analysis.\n",
"\n",
"We will take into consideration four *variables* for different countries over different days:\n",
"- *deaths*: number of reported deaths from covid19\n",
"- *confirmed*: number of confirmed cases of infection from covid19\n",
"- *treatment*/*measure*/*policy*: the implementation of a specific policy\n",
"- *hospital beds*: number of hospital beds for 1000 citizens\n",
"\n",
"**The selection of these variables is arbitrary and should be enriched.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Importing libraries"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import dowhy as dy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Loading data\n",
"\n",
"We will use three data sets:\n",
"- **HDE**/HDE/acaps-covid-19-government-measures-dataset: this dataset provides information on the adoption of different measures to contain covid19.\n",
"- **johns_hopkins**/johns-hopkins-covid-19-daily-dashboard-cases-over-time: this dataset provides information on the effects of covid19.\n",
"- **world_bank**/: these datasets provide background information about a country."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"HDE_df = pd.read_csv('HDE/acaps-covid-19-government-measures-dataset.csv')\n",
"JHU_df = pd.read_csv('johns_hopkins/johns-hopkins-covid-19-daily-dashboard-cases-over-time.csv')\n",
"WB_df = pd.read_csv('world_bank/hospital-beds-per-1-000-people.csv')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Finding the common countries\n",
"\n",
"We perform a little bit of manual alignment on the countries contained in the two datasets.\n",
"\n",
"**This alignment is incomplete**"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def replace_HDE(original,replacement):\n",
" HDE_df.loc[HDE_df['country']==original, 'country'] = replacement\n",
"\n",
"replace_HDE('United States of America','US')\n",
"replace_HDE('Czech Republic','Czechia')\n",
"replace_HDE('Korea Republic of','Korea, South')\n",
"replace_HDE('North Macedonia Republic Of','North Macedonia') \n",
"replace_HDE('kenya','Kenya') \n",
"replace_HDE('Viet Nam','Vietnam') \n",
"replace_HDE('Russian Federation','Russia') "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"countriesHDE = set(HDE_df['country'].unique())\n",
"countriesJHU = set(JHU_df['country_region'].unique())\n",
"countries = list(countriesHDE.intersection(countriesJHU))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Finding the date range\n",
"\n",
"We convert the dates in the HDE dataset in datetimes, and we extract the range of dates from the JHU dataset.\n",
"\n",
"**Some dates in the HDE dataset are inconsistent** (hence the *errors='coerce'* option)\n",
"\n",
"**We rely on the JHU providing the same range of data for all countries**"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"HDE_df['date_implemented'] = pd.to_datetime(HDE_df['date_implemented'],format=\"%m-%d-%y\",errors='coerce')\n",
"\n",
"dates = JHU_df.loc[JHU_df['country_region']==countries[0],'last_update'].values"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Defining the measure/policy of interest\n",
"\n",
"We define a measure/policy of interest, which will be *treatment* of interest in the causal analysis. Use *list(HDE_df['measure'].unique())* to list the 32 measures tracked in the HDE dataset."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"treatment = 'Introduction of quarantine policies'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Creating a new dataframe\n",
"\n",
"We create a new dataframe, cast the dates in datetime, and add a new binary column that will track when the chosen policy of interest is activated."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"df = JHU_df.copy()\n",
"df['last_update'] = pd.to_datetime(df['last_update'])\n",
"df[treatment] = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set the new column/field as active (=1) from the day **after** the policy has been called for by a government."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Country Switzerland has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Djibouti has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Bhutan has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Trinidad and Tobago has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country US has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Rwanda has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Mongolia has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Libya has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Tanzania has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Pakistan has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Burkina Faso has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Turkey has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Zimbabwe has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Dominica has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Venezuela has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Morocco has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Luxembourg has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Maldives has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Singapore has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Ireland has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Lithuania has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Portugal has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Bolivia has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Kazakhstan has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Indonesia has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Antigua and Barbuda has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Papua New Guinea has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Slovenia has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Mexico has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Saudi Arabia has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Togo has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Timor-Leste has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Eswatini has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Eritrea has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Senegal has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Guinea-Bissau has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Botswana has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Zambia has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Denmark has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Nicaragua has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Sweden has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Japan has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Liechtenstein has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Kuwait has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Cambodia has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Belgium has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Sudan has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"Country Paraguay has not implemented policy Introduction of quarantine policies to date 2020-03-30\n",
"\n",
"113 countries have implemented policy Introduction of quarantine policies to date 2020-03-30\n"
]
}
],
"source": [
"ncountries_treatment = 0\n",
"for i in countries:\n",
" try:\n",
" activationdate = HDE_df.loc[(HDE_df['country']==i) & (HDE_df['measure']==treatment), 'date_implemented'].values[0]\n",
" df.loc[(df['country_region']==i) & (df['last_update'] > activationdate), treatment ] = 1\n",
" ncountries_treatment += 1\n",
" except IndexError:\n",
" print('Country {0} has not implemented policy {1} to date {2}'.format(i,treatment,dates[-1]))\n",
" \n",
"print('\\n{0} countries have implemented policy {1} to date {2}'.format(ncountries_treatment,treatment,dates[-1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Retrieving further data\n",
"We extract more background information, specifically considering what could be *confounders* between our treatment and effect.\n",
"In this notebook we take the last recorded value of *hospital beds per 1000 people* from the *world_bank* data set."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Country Egypt has no reported hospital beds between 1960-2019\n",
"Country Saint Lucia has no reported hospital beds between 1960-2019\n",
"Country US has no reported hospital beds between 1960-2019\n",
"Country Bahamas has no reported hospital beds between 1960-2019\n",
"Country Kyrgyzstan has no reported hospital beds between 1960-2019\n",
"Country Syria has no reported hospital beds between 1960-2019\n",
"Country Korea, South has no reported hospital beds between 1960-2019\n",
"Country Venezuela has no reported hospital beds between 1960-2019\n",
"Country Iran has no reported hospital beds between 1960-2019\n",
"Country Czechia has no reported hospital beds between 1960-2019\n",
"Country Saint Vincent and the Grenadines has no reported hospital beds between 1960-2019\n",
"Country Russia has no reported hospital beds between 1960-2019\n",
"Country Slovakia has no reported hospital beds between 1960-2019\n",
"Country Liechtenstein has no reported hospital beds between 1960-2019\n",
"Country Gambia has no reported hospital beds between 1960-2019\n",
"\n",
"146 countries have reported hospital beds between 1960-2019\n"
]
}
],
"source": [
"ncountries_confounder = 0\n",
"\n",
"confounder = 'hospital beds'\n",
"df[confounder] = np.nan\n",
"for i in countries:\n",
" df1 = WB_df.loc[WB_df['country_name']==i]\n",
" beds_nan = df1.iloc[:,4:].values\n",
" try:\n",
" beds = beds_nan[np.isfinite(beds_nan)][-1]\n",
" df.loc[(df['country_region']==i), confounder] = beds\n",
" ncountries_confounder += 1\n",
" except:\n",
" print('Country {0} has no reported {1} between 1960-2019'.format(i,confounder))\n",
" \n",
"print('\\n{0} countries have reported {1} between 1960-2019'.format(ncountries_confounder,confounder))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Conclusion\n",
"\n",
"In this part we collected data that we will be fed in our model.\n",
"\n",
"**Richer and more relevant data should be collected**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Model Part\n",
"In this section we rely on the *dowhy* library to build a causal model.\n",
"\n",
"**This sort of data implicitly requires a time-dependent analysis, but we will analyze data statically on a per-day basis.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Preparing the data for the causal model\n",
"We extract the data that is relevant for our analysis and we drop lines containing *nan*.\n",
"\n",
"**Missing values may be dealt better than just dropping**"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Dropped 2208 rows with nan values\n"
]
}
],
"source": [
"data = pd.DataFrame(data=\n",
" {'confirmed': df['confirmed'],\n",
" 'deaths': df['confirmed'],\n",
" 'hospital beds': df[confounder],\n",
" 'Introduction of quarantine policies': df[treatment].astype('bool'),\n",
" }\n",
" )\n",
"rows1 = data.shape[0]\n",
"data = data.dropna()\n",
"print('Dropped {0} rows with nan values'.format(rows1 - data.shape[0]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Setting up a Causal Model\n",
"\n",
"We now set up a causal model as illustrated by the graph below."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:dowhy.causal_graph:If this is observed data (not from a randomized experiment), there might always be missing confounders. Adding a node named \"Unobserved Confounders\" to reflect this.\n",
"INFO:dowhy.causal_model:Model to find the causal effect of treatment ['Introduction of quarantine policies'] on outcome ['deaths']\n"
]
}
],
"source": [
"model= dy.CausalModel(\n",
" data = data,\n",
" treatment = treatment,\n",
" outcome = 'deaths',\n",
" graph = \"\"\"graph[directed 1 node[id \"confirmed\" label \"confirmed\"]\n",
" node[id \"deaths\" label \"deaths\"]\n",
" node[id \"Introduction of quarantine policies\" label \"Introduction of quarantine policies\"]\n",
" node[id \"hospital beds\" label \"hospital beds\"]\n",
" edge[source \"confirmed\" target \"deaths\"]\n",
" edge[source \"hospital beds\" target \"deaths\"]\n",
" edge[source \"Introduction of quarantine policies\" target \"deaths\"]\n",
" edge[source \"hospital beds\" target \"Introduction of quarantine policies\"]\n",
" edge[source \"Introduction of quarantine policies\" target \"confirmed\"]]\n",
" \"\"\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAFbCAYAAADWYvcBAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeVxUZf8//tfAsO/ILnuAgKICyqLglgIqLpUV5ZKm4l2Z5rdutfJO7/Qu9VN3t3dlrpnmkpqpKe5gyaKQCKIsorLv+zoDDMz794e/ObejLAMCB/B6Ph7zgDNz5pz3nJnrel/nnOtcR0BEBIZhGIZhBhwlvgNgGIZhGKZnsCTPMAzDMAMUS/IMwzAMM0AJ+Q6gL6utrYVIJEJ9fT0qKysBAFVVVXi8G0NDQwPEYrHc+7S0tKCqqspNKysrQ1dXF0pKStDT04O2tja0tLSgpaXVOx+EYRjeyOoOsViMhoYGEBGqqqq41xsbGyESidp8v0gkQmNjY5uvy+qXtqioqEBbW5ub1tbWhoqKCoRCIXR0dAAAOjo6EApZOhiInptvtby8HLm5ucjPz0d5eTnKyspQWlqKkpISlJWVoaysDOXl5aipqYFIJEJ1dXWvxKWvrw8tLS3o6OjAyMiIe5iYmMhNW1lZwcrKCnp6er0SF8M8T6qqqlBZWcn9ra2thVgsRk1NDerq6iASiVBXV4eamhqIxWLU19ejqqoKYrEYYrGYS+TV1dWQSqUdJu6+6vEGgZ6eHpSUlKCjowMNDQ1oa2tDV1cXmpqa0NTUhL6+PjQ0NLj/NTU1oaGhAT09Pejp6cHAwAD6+vowMDCAsrIyz5/s+SUYKL3rq6urkZ6ejvT0dDx48AC5ubncIzs7W67Aqaurw8jICMbGxk8lU9mPWLbHrampCW1tbe4HL2sFyzzeGpapqalBS0sLNy0r8C0tLaipqUFtbS3q6+shEolQWVkJkUiE2tparrFRUlKC0tJSbvrxVryOjg6sra1hY2MDS0tLWFtbw8HBAY6OjnBycpJrsTPM86aurg7FxcVcGZI9Kisr5ZL4k9OtVYOyPWRtbW1oaGhAR0en3YQnEAi4PeLHk6Wuri6UlZWhpqYGTU1NAODmB8Ad4WvLk3viT2rtaOLj6uvr0dTUxE3L6qempibU19cD+F/j5PFlyRoujzdsqqurIRaLIRKJUFVVBZFIBLFYzP3f1hEHXV1dLuHLHk9Om5mZcfWxiYkJBg0a1OZnYhTX75J8WVkZEhISkJiYiPT0dNy/fx9paWkoLi4GAKiqqsLOzo7b87W2toa1tTU3bWVl1e8Ok9fU1CAvLw85OTlyDRfZ/1lZWWhubgYAWFpawsnJCY6OjnB2dsaIESPg7u4OfX19nj8Fw3RNY2Mj8vPzkZ+fj9zcXBQVFaGwsJBL4CUlJSguLkZpaSkaGhrk3qujowNjY2O5ZNJeopFN6+joyJ1yYxQjlUpRXV2N6urqpxpU7TWyKioqUFpaKrcsoVAIY2NjGBsbw9TUFCYmJtz04MGDuYeVlRXbuWlHn07yBQUFiI2NRWJiIhISEpCQkIC8vDwAgIWFBVxcXODk5AQnJycMGTIETk5OsLGxee7OLUkkEmRmZuLevXu4d+8e1/hJSUlBSUkJAMDOzg7u7u4YOXIk3N3d4e3tDWNjY54jZ553zc3NyMnJQVZWFvLy8pCbm4uCggLk5uYiLy8PBQUFXAMeeLRXa2pqyu31ySp9MzMzGBsbw8jICKampjA1NYWRkRHU1dV5/HRMZzQ3N8udRpU13EpLS1FUVCTXqCsoKJA7eqGnpwdLS0tYWlpyiV/WCLC1tYW9vf1z+1voU0k+IyMDUVFRiI6ORlRUFFJTU0FEMDc3h6enJ/cYPXo0zMzM+A63X6isrERycjLi4+O5R1paGqRSKezt7TF27Fj4+flh7NixcHV15Q4hMkx3qaysREZGRquPnJwc7iiUmpoaDA0NYWFhAXt7e5ibmz/1v7W19XPXiGdaJxaLUVhYiIyMDBQUFDz1v+yvjIGBAezt7Vt92NraQklpYF5sxmuSLyoqwvnz53HhwgWEh4ejvLwc2tra8Pb2hp+fH/z8/ODt7f3UOW/m2VRVVSEmJgbR0dGIjIzEzZs3IRaLYWZmhsmTJ2Pq1KkICAiAkZER36Ey/URLSwsyMzORkpKC1NRUpKSkICUlBQ8ePOB6kisrK8PKyqrNSpYdWWK6W11dHTIzM1ttYGZmZnJ9CNTV1eHg4ABnZ2e4uLjA1dUVLi4ucHZ2hpqaGs+f4tn0apInIsTGxuLMmTO4cOECEhISoKamhvHjxyMwMBD+/v4YOXIka6n3sqamJty8eRORkZG4ePEioqKiIJVKMWrUKEybNg3BwcHw8PDgO0ymDyAiPHz4EImJiUhNTUVycjLS0tKQlpbGVZjW1tZcRenk5MQlchsbG7lOqwzDJyJCfn4+l/Rl/buSk5Px8OFDNDc3Q1lZGXZ2dlzSd3V1xfDhwzF06NB+81vulSSfnJyM48eP49ChQ3jw4AHs7OwwZcoUTJ48GUFBQWxPvY8RiUSIiYnBmTNncOrUKeTk5MDW1hYzZ87EwoUL4e7uzneITC8pKCiQO9Vz48YNlJWVAQDMzc0xdOhQuLq6cn9HjBjByjPT70kkEuTm5iI5ORkpKSnc39TUVIhEIqioqMDR0VHuNLKHhwd39URf0mNJvqioCLt378bBgweRnp4OOzs7vP766wgJCcGIESN6YpVMD4mNjcXRo0dx7Ngx5OfnY8SIEZg/fz4WLVoEQ0NDvsNjuklhYSEiIyMRHR2Nmzdv4vbt26ivr4eqqiqGDRsGDw8PeHh4wN3dHSNGjICGhgbfITNMr2ppaUFaWhpu3brFPRITE1FTUwOhUAhXV1d4eHhwfZ2cnZ35Drn7k3x0dDS+++47/Pbbb9DV1cW8efMQEhICLy8v1qmrn5NKpYiMjMQvv/yCI0eOoKmpCW+++Sbee+89tnffDz148ABRUVG4du0aoqKicP/+fQiFQowcORKjR4/mkvqwYcPY5WQM0wYiwv3793Hr1i0kJCTgr7/+QmxsLEQiEUxMTODn54dx48bB398fI0aM6PWBgbolyRMRTp48iU2bNiEhIQGjRo3C8uXL8frrrz+3ly0MdHV1dfj555/x/fffIzk5GWPHjsVnn32GgIAAvkNj2lBeXo6LFy8iLCwMV69eRWFhITQ0NDB69GiMGzcOfn5+GDNmDDvczjDPSCKRID4+HlFRUYiMjERUVBQqKiqgo6MDf39/TJ06FdOnT4ednV2Px/LMST4sLAyfffYZEhMT8corr+DDDz+Et7d3d8XH9ANXr17FV199hXPnzsHf3x//+te/4O/vz3dYDICkpCSEhYUhLCwMN27cgJKSEvz9/TFlyhT4+/tj9OjRbC+dYXoYESE5ORmRkZGIiIjApUuXUFNTAxcXF0yfPh3Tp0+Hn59fz3Q6py66e/cujR07lgQCAc2cOZMSExO7uihmgIiOjqZJkyYRAJo6dSplZGTwHdJz6fbt27Rq1SqysrIiAGRqakqLFi2i48ePU3V1Nd/hMcxzr6mpicLDw+nDDz8kZ2dnAkB6eno0d+5cunjxIrW0tHTbujq9Jy+RSLB582Zs2rQJHh4e2LZtG7y8vLq/9cH0WxEREVi5ciUyMzPxxRdfYPny5QN2oIm+oqSkBIcPH8b+/fuRmJiIF154AXPnzkVwcDA8PT3Z9meYPuzhw4cICwvD0aNHERMTA0tLS8ydOxdvvfUWXFxcnm3hnWkRpKen0/Dhw0lTU5O+/vpram5u7rbWBjOwNDY20vr160lVVZXGjh1LeXl5fIc0IF25coVmzpxJKioqpKenR0uWLKHIyEiSSqV8h8YwTBfcu3eP1q1bR9bW1gSAvLy8aN++fdTY2Nil5Smc5P/8808aNGgQeXl50YMHD7q0MkXV1tb26PL7qpqamn61XEUkJSWRq6srWVpaUkJCAm9xDCRSqZSOHj1KI0eOJAA0ceJEOnz4MIlEIr5D4w2fv/GeMNA+T28YaNuspaWFIiIiaO7cuaSqqkrm5ua0detWqq+v79RyFEryBw8eJFVVVZozZ06PViQHDx6kF198kczMzHpsHYo6efIkWVpaUkpKSpeXIZFI6Nq1a/TJJ5/QhQsX2pzvu+++Iz8/P3J1de3yulqzY8cOGjduHA0ePLhbl9tZlZWVNHnyZNLW1qbz58/zGkt/d+XKFXJ3dyclJSUKCQmh+Ph4vkOiw4cPk5+fHwEgADRnzhy6du0a93plZSWtW7eONDU1SVtbm9atW0dVVVXdsu6eKjt86WyZraqqok8//ZT8/f1p6NChNH36dJoxYwatXr2aPv74Y/r222+7PUapVEr//ve/6csvvyQHBweaN28er0d1+0o915Py8vJo9erVpK2tTebm5rRjxw6Fz9t3mORPnDhBysrK9Pe//73HDwE2NzfThAkTyMjIqEfXo4hLly6Rh4fHU53HCgoKFF5GTEwMLVq0iADQnj172pxPIpGQm5sbOTs7dzne1jQ3N5Ofn1+faDQ1NTXRW2+9RVpaWmyPvgtKS0tp3rx5BICCg4MpKSmJ75Dk5OTkEIB2K9qVK1fS2rVru3W9PVV2+NKZMvv777+TmZkZjR07Vq6eqqiooPnz5xMA2rJlS7fHuGHDBlq2bBkREUVGRtKMGTNILBZ3+3oU1ZfquZ5WXFxMK1asIBUVFfLx8aE7d+50+J52e+Pcu3cPCxYsQGhoKLZu3drjg9koKyvD0tKyR9ehqClTpiA+Pl7uOsbKykrMmzdP4WX4+vri/fff73A+oVCIwYMHdynO9vSl7amiooI9e/bA19cXs2fP5m6By3Tsr7/+gqenJ/7880+cPn0aZ86cgZubG99hyZFdW9/efb1tbGy6/ffYU2WHL4qW2aioKLzyyiuwtrZGRESEXD1lYGCAAwcO4PXXX4dIJOr2GLdv3w5bW1sAgJ+fH37//Xdex0PpS/VcTzMxMcG2bdsQHx8PgUAALy8v/PTTT+2+p80kT0RYtGgRXF1dsW3btu6Otd8RiUQICQlBRkZGp97HrkH+H6FQiKNHj0IoFCI0NJTvcPqFy5cvY/z48XjhhRdw8+ZNzJw5k++QWiXbAWhvR0BVVZXdfKqbvP/++5BIJNi4cWObdcznn3/e7Um+oaEBJSUlbPRSnrm5uSEyMhKrV6/G4sWL8fnnn7c5b5tJ/uTJk4iNjcXOnTt5udtOUVERZs+eDUNDQ3h6eiI1NVXu9RMnTmD58uX46KOPMHXqVKxbt467CxYAJCYmYtGiRdiyZQtmzZqFKVOmAABSUlLw6aefwtXVFQUFBdw6vLy8cOPGDQCP9tj37t2LKVOm4NSpUwAebY/U1FSUlZVh6dKl+OqrrwAAxcXFWLp0KTZu3IilS5fipZdeQnl5eZc/9x9//IGgoCAYGhoiMDBQrlFBRNixYwfeeecdeHt7IyAgAPfv35d7/+nTpxEaGoo1a9bg/fffl7ufcnvbpbcYGhrihx9+wOnTpxETE9Or6+5vEhMTMWvWLMyZMweXL1+GiYkJ3yF1m8TERPz973+Hvb096uvrsWTJEhgZGcHLy+uphnRHZV2mvbLT3u++vXKVn5+PzZs3Y9iwYaioqEBgYCBsbGxw8uRJGBsbQyAQYN26ddyywsPDoauri/Xr13e4bJmOyuyT7t69i8TEROjr67c7wqSTkxPeffddhbajIt/H/v37sXTpUgDA8ePHsXTpUmzZskWh5R85cgS6urqwsrICAFRXV2Pjxo1QVlaGr6+vwjEous268p2Wl5fzXj92hrKyMjZs2IDt27djw4YN2LVrV+sztnUcf8qUKTRjxozuOpWgsHnz5pGWlhZ98MEHlJaWRklJSaSlpUXBwcHcPN988w2NGTOGmpqaiIiorKyMHB0dafz48Vy/AScnJ4qKiiIiIpFIRH5+fkREtHbtWtLX1ydlZWVatWoVXb16lU6cOEFGRkakqalJBQUFlJKSQqtWrSIA9Ouvv3LrDQ4OJltbW7l4J0yYQK+//jo3PWLECJo3bx43fffu3Q7PyRMRBQUF0aBBg+jtt9+m8+fP09dff02qqqpkYWHB9ab88ssv6aeffiKiR+ehXF1dyczMjHv90KFD5O3tzZ0fKy0tJSMjI7lzVW1tl97m5eUlt50Yec3NzTR06FCaOHFitw6M0VOqqqoIQLvnxr/77jvasWMHEREVFhbS5MmTCQC99957lJycTAkJCaSmpkYhISHcexQp64qUnfZ+9+2Vq/Pnz5OzszMpKyvT+vXradeuXeTl5UX5+fn01VdfEQD67bffuGVJJBLy9/fnYuuOMvukvXv3EgDy9PTs8HtRdDsq+n2UlZURANq0aVOnlk9EFBAQQJaWlnLvc3NzIx8fHyJS/DehyDbr6nfaV+rHzpJ1bG1tALJWk7xYLCY1NTU6cOBAjwf3pHnz5pGenh5JJBLuuYkTJ5K5uTkRPep4oKWl9VRs+/btIwD0888/U1NTEwkEAtq2bRv3+smTJ7n/33zzTVJRUeF+kEREx48fJwD02WefERHRH3/8oVCSnzhxIn3xxRfc9Ny5c2n48OHcdGeSvIWFhdxzX375JQGgbdu2UX5+PpmamspV+J999hkBoF9++YXq6+vJ3NycDh8+LLeMl156ifvxd7RdetNXX31FJiYm7HruNpw6dYqUlJTo/v37fIeikM4meSKijz/+mABQWVkZ95yfnx85OjoSkWJlnajjstPe776jckVEtHjxYgLw1HdRV1dHhoaG9Morr3DPnT17lr7//nuFlq1ImW3N1q1bCQAFBAS0Oc/jFN2OHX0fRK0neUWXP3v27KeSvI+PD5fkFYlBkW3W1e+0L9WPnSWRSMjOzo5WrVr11GutniDLzMxEY2MjRo4c+WzHE7pIRUVF7tydvb09rl+/DgC4ceMG6uvrYW1tLfee4OBgAI/GUZ83bx4CAwPxwQcf4O7du9i8eTNmz57NzaupqQllZWW50xCzZ8+Gmpoa7ty5AwAKnzuMiIgA8Ohc1aFDhxAXFwfq4u0AdHV15aYXLFiAjz/+GPHx8bCwsIBEIsGyZcvk5lmyZAk0NDQQGRmJwsLCpzpkqampcf+rqKi0u116k7u7O0pKSlBeXg4jIyNeYujLIiIi4O3tDQcHB75D6TGyu3E9XtYsLS3x4MEDAIqXdaD9stPe7z4mJqbdcgX8rz568rvQ0tLCggUL8P3336OsrAxGRkY4evQo14epo2UrUmZbIzvknZWV1e58Mopux46+j2ddviI6ikGRbdbV77Qv1Y+dJRQK8frrr+P8+fNPv9baG+rr6wG031O2Nz3eySM7OxsAUFFRITePkZERNDU1UVBQAODR+aGlS5di9+7dOHnyJI4dO4aJEye2uQ6hUAgLCws0Nzd3KraWlhZs3boVN2/exIoVK+Dt7c2d239WFhYW0NDQgFgsRmpqKrS0tLB79+5W55VVLB119Ovsdukpst9WXV0dS/KtKC8v71fn4DU1NaGkpASJRNLmPA0NDTA3N1d4mYqW9dY8XnaAtn/3HZWrjoSGhuI///kPDh48iIULF0JZWRkGBgYA0G1l9kmyYU4zMjLQ3Nzc4Q7Js2xHRfT08h+XlpYGoP1t9izfaV+pH7vC1NQUpaWlTz3fasc7Y2NjAI86v/U1sktF2url7uzsDOBR0j506BAOHToEoVCIoKCgpzrvPUkkEnHvV4RUKsW0adOQkpKCEydOYPz48Qq/V1ECgQDDhg2DpqYm8vLykJeX99Q8paWl3I9eVuDa0pXt0hOKioogEAi43xojz97eHnfv3u3yUaHepqKiAktLy3Y7nZaWlsLU1FThZSpa1tsiKztA27/7jspVR1xcXODv748ff/wRR48exdy5c7nXuqvMPmno0KEYMmQImpubERUV1eH8z7od+V7+4xTZZs/ynfaV+rErbt++3eqRv1aTvLW1NYyNjREZGdnjgXWWr68vdHV1uV7vMnl5eRCJRJg5cyYaGxu5noZvvvkmbty4ASLC1atX21xuYWEhSktLMWfOnDbnUVJSQl1dHTcdFxeHS5cuYcKECdxzEomk2yrmrKwsSCQSvPbaa3BzcwMRYc2aNXLzPHz4ENu3b8fw4cMBAEePHpV7XSqVoqWlBQC6tF16yrVr1+Dq6gotLa1eX3d/ILtc88nfeV82btw4VFVVIS4u7qnXpFIpYmJiMHbsWIWXp0hZb8vjZae9331H5UoRoaGhuHPnDg4cOIBJkyZxz3dHmW2NUCjkru75+OOP0dTU1Op8RUVF2L9//zNtxye1VrcpunyhUIi6ujq5z1ZXVwepVKrw+hXZZl39TvtS/dhZubm5OHbsGN54442nX2zrRP6yZcvIxcWl13v2Tp06lYRCodxg/DNnziSBQMD1SP3hhx9IIBDQlStXuHn+/ve/01tvvUVERA0NDeTu7s4NtdjU1ERGRkZ0/fp1IiJasmQJCQQCudvjvvvuu7Ro0SJuWtYR74cffuCe+9vf/kYA6ObNm3T16lWKiIggAOTv709JSUm0d+9eGjZsGGlra9Pt27epqKiIYmJiuA5A7Zk+fTqZmppSXV0dET0aOvLtt9/mOrhIpVIaPXo0AaCXX36Zfv75Z/r+++/pxRdfpNLSUiJ61AlQWVmZtm/fTvX19RQXF0cWFhYEgA4fPkwVFRXtbpfeIhaLydTUlOvkyLTurbfeIhMTE8rKyuI7FIVkZWWRjo4OWVtb0+3bt7nnMzIyaP78+fTf//5Xbv7333//qU5WkyZNIl1dXW66o7JO1HHZaa8+UKRczZs3jwQCAVVWVrb6ucViMRkYGND69evlnu+OMtveOOWbNm0igUBAvr6+FBcXxz1fWVlJR44cocmTJ1N+fr7C21GR7yMhIYEA0CeffCIXiyLL/+c//0kAaOPGjXTv3j3auHEjOTo6kp6eHt26dUvhGDraZnV1dV36TjvKG32VWCwmf39/cnV1bXXY+TaTfHJyMgmFQtq7d2+PBvi4AwcOkIGBAQGglStXUnV1Nf34449kaGjIPSdL/qdOnaKAgABavnw5/eMf/6Cvv/6a66nd0NBAo0ePpsDAQNq8eTOFhobS7t27ufUsWbKEVFVVadWqVfTqq6/S4sWLaePGjdz7w8PDady4cQSARo0aRZcuXSKiR/fptrS0JCcnJzp+/DgRPUr8Ojo65OPjQ1euXKFz586RkZERzZkzhyIiImjq1KkEgDw8PCgsLKzNz56UlEQhISEUGBhIoaGhtHLlSrme/URE5eXlNHfuXDIxMSFjY2NasGABV4iJiKqrq2nRokVkampK1tbWtGHDBgoNDaVFixbRlStXSCQStbtdesuXX35JWlpanRoi+HlUU1NDI0aMIFtb2x6/KVR3yc3NpQULFtCQIUPI0tKSXFxc6MUXX6TTp0/LzXflyhWytbUlAPTuu+9SSUkJHThwgLS1tQkAbdiwgats2yvrRB2XnY7qg/bK1a5du8jY2JgA0Pz587lk9KSNGzdSYWHhU88/a5ntaCcrMTGR3n77bbKxsSEjIyMaPXo0TZgwgX744Qe5K5Q62o6KfB9xcXH0xhtvEACys7OjQ4cOyd2DoKPvqbq6mmbMmEHa2trk4+NDf/31Fy1cuJDmzZtHv//+u8K/CUW2WVe+045+J31RbW0tBQUFkYGBQZtDXbd7P/kPPvgA+/fvR0JCAjeM4UCwdOlSHDx4kOuUw/SexMRE+Pr64rPPPsPHH3/Mdzh9Xnl5OYKCgvDgwQP89NNPmDVrFt8hMQzTB6SkpODVV19FaWkpzp8/D09Pz1bna3fs+s2bN8PGxgbBwcFP9ZxkmM7Kzc3FzJkz4evri9WrV/MdTr8waNAgREZGYs6cOZg9ezbeeOONDkdEYxhm4GpoaMD69evh4eEBfX193Lp1q80ED3SQ5NXV1XH27FnU1tZiwoQJyMnJ6faA+VBXV9etHeSYjqWkpGD8+PHQ09PDiRMnuOthmY6pq6tj9+7diIiIwK1bt+Do6Ii1a9eiqqqK79AYhuklUqkUx48fh6urK7766iv885//xLVr1zq8OU+7SR54NBDB9evXoaKiAm9vb/z111/dFjQffvjhB1y+fBktLS0IDQ1V6BIU5tlERERg7NixMDMzQ3h4OHcdMdM5EydORGJiItavX489e/bA1tYWq1evbvVSIYZhBob6+nr897//hYODA+bOncudvluzZo1iO0uKnuCvrq6moKAg0tTUpG+++YbrFMMwbWlsbKQNGzaQiooKvfHGG9TQ0MB3SANGVVUVbd68mSwsLEhFRYVmzZpFJ06ckLsqhWGY/kkqlVJUVBQtW7aM9PX1SUtLi5YvX97q2PQdUTjJEz0aH3fDhg2kqqpKPj4+lJyc3OkVMs+H2NhYcnNzIy0tLfrPf/7DxqjvIY2NjXTo0CEKCAggZWVlMjQ0pPfee49iY2P5Do1hmE7KzMykzz//nBwcHAgADR8+nL7++mu5Swo7q93e9W25e/cuFi9ejNu3b+Pdd9/Fxx9/zEYuYwAAOTk52LRpE3788UdMnDgRu3bt4kbEYnpWfn4+Dh48iAMHDiAlJQWOjo6YMWMGpk2bBn9//04Pn8owTM8iIty+fRthYWE4e/YsYmNjYWJigjfffBMLFizolvvHdCnJA4/GbN+xYwc2bdqEuro6rFixAh999BE73/qcKioqwpdffomdO3fC3Nwcn3/+OebNmyd33wGm9/z111/49ddfcfbsWaSkpEBXVxdTpkzBtGnTMH369E4NL8swTPepr6/HlStXcO7cOYSFhSE/Px/m5uaYPn06Zs+ejcDAQIVvkKaILid5GZFIhO+++w5bt25Fc3Mz3n77bbzzzjtwdHTsrhiZPuzOnTv4/vvv8fPPP8PAwACffvopFi9ezPYa+5DMzEycO3cOZ8+exR9//IGmpia4ublh3Lhx8PPzg7+/f6duHMMwjOJqa2sRExODqKgoXLt2DXFxcWhqasKoUaMQHByMadOmwcPDo8d2iJ45ycvU1NTghx9+wI4dO5CTk4OAgAC89957mDZtGpSUOuzEz/QjEokEJ0+exDDsWpEAACAASURBVPfff49r165hyJAhWL58ORYvXszdypHpm0QiEcLDwxEREYGoqCgkJiaiubkZDg4O8PPzw7hx4zB27Fg4OTnxHSrD9EvFxcWIiopCZGQkIiMjcfv2bbS0tMDJyQl+fn4YP348goKCeu0uk92W5GWkUikiIiKwbds2hIWFwdzcHHPmzMGrr76KsWPHssO3/ZTsBiPHjx/H0aNHUVpaikmTJmHFihUIDg5m32s/VV9fj4SEBERHR+PKlSuIjo6GWCyGrq4u3Nzc4OnpyT1cXFxYg51hHlNZWYnk5GTEx8dzj9TUVCgpKWHIkCHw8/PD2LFjMWHCBFhbW/MSY7cn+cfdv38fhw4dwi+//IJ79+7Bzs4OISEheOmll+Dp6ckqjD6upaUFN27cwIkTJ3Ds2DHk5+djxIgRCAkJwdy5c2FlZcV3iEw3a2pqQnx8PG7evIlbt27h1q1bSElJQXNzM3R1dTFy5Eh4eHjAw8MDrq6ucHZ2ZncSZAY8iUSChw8fIjk5Gbdv3+bKhmz0SXt7e65ceHp6wsfHB7q6ujxH/UiPJvnHJSQk4MiRIzh27Biys7NhbGyMwMBATJ06FQEBATAyMuqNMJgOFBUV4fz587hw4QIuX76MyspKODo6IiQkBCEhIXB1deU7RKaXNTQ0ICkpiavYbt26hTt37qCpqQkCgQA2NjZwdnbG0KFD4ezsDFdXV7i4uLBOuEy/IxaLkZaWhrS0NKSkpCA1NRWpqam4f/8+JBIJlJSU4OjoyCV02UNfX5/v0NvUa0n+cUlJSbhw4QIuXLiAqKgotLS0YNSoUVwnoDFjxvTa+YrnXUFBAaKiorhOIUlJSVBTU8O4ceMwdepUTJ06FUOGDOE7TKaPaW5uRkZGBpKTk+UqxLS0NNTX1wMAzMzM4OTkBHt7e9jb2+OFF17g/mflm+FLTU0NMjIynnrcv38fWVlZkEqlUFFRgYODA3e0StaAdXZ27nf9jnhJ8o+rqalBeHg4wsPDERkZibt370IqlcLZ2RljxozBmDFj4O7ujmHDhrEe289ILBbjzp07uHXrFmJiYhAdHY2MjAwIhUKMHDkSfn5+mDJlCiZMmABNTU2+w2X6ISJCdnY20tLSkJycjAcPHnCVaE5ODpqamgAA2traXMK3t7eHra0trK2tYWFhAUtLS5iamrLTeUyXlJeXo7CwEDk5OSgoKEBmZqZcMi8rKwMACAQCDB48mPsNOjg4cEeiHBwcoKKiwvMn6R68J/knVVdXIzo6GtHR0YiMjER8fDxEIhFUVFTg6uoKd3d3jBw5EiNGjICrqyvbI2hDfn4+0tLSkJCQgMTERCQmJiItLQ0tLS3Q0dGBl5cX/Pz84OfnBx8fH2hra/MdMjPAtbS0IC8vr9W9qKysLJSUlHDzqqiowMzMDFZWVhg8eDAGDx4MKysrrhFgbGwMMzMz6Onp8fiJmN4kEolQUlKCoqIiFBUVITc3F/n5+cjPz0dubi4KCgqQl5cndwtxXV1d2Nraws7OTq5RaW9vDzs7O6ipqfH4iXpHn0vyT2ppaUF6ejoSExORkJDAJS1Za0xfXx9OTk4YMmQInJ2d4eTkBAcHB1hZWWHQoEE8R9+zSkpKkJubiwcPHuDevXtIS0tDeno60tPTUVtbCwAwNzfHyJEjMXLkSLi7u8Pd3R0vvPAC6w3P9DmNjY2tVtp5eXkoKChATk4OioqK0NLSwr1HTU0NxsbGMDExgampKYyNjWFsbAxzc3MYGxvDyMgIxsbGMDQ0hL6+PgwMDNgRgj6iuroaVVVVqKysRHl5OYqLi1FaWorS0lIUFhZy/xcXF6OkpIQ7DSRjYmLCNfpkDcAnG4Vs56UfJPm2FBQUcEnt9u3b+PXXX6GiooKSkhKuEtDU1ISNjQ0sLS1haWkJa2trDB48GCYmJjAyMpKrAPoKIkJZWZnco7i4GAUFBcjOzkZeXh5yc3ORm5uLhoYGAIBQKIS9vT2GDBmCIUOGwMnJCU5OTnBxcWFHOpgBpaWlhUsGrSWCx5NEWVmZ3F6djJ6eHpfwDQwMuP9lf2UNBy0tLejp6UFDQwOampowMDCAhoYGNDQ0+nRHq55WV1cHsViM2tpa1NbWQiQSob6+HtXV1RCLxaivr0dlZSUqKyu5JN7a/1KpVG65QqGQa6SZmZnBxMSk1QacmZkZzMzMnou98O7Qb5O8zMOHDxEcHIza2lqcPXsWLi4uyMzM5JJhTk4OcnJyuOn8/HzU1NTILUMoFHJJX1tbG9ra2tDT04OmpiY0NTWhr68PLS0tqKqqQklJ6alDhPr6+tyesVQqRXV1tdzrVVVVICI0NjZCJBKhsrISIpEIIpEINTU1qKmpQV1dHZfUn/zxGxgYwMLCAtbW1rCysoKlpSUePnyII0eO4IsvvsCKFSsGzPkjhulOdXV1KC0tfSrRtJZ8ysrK8PDhQ1RVVUFDQwMikajdZWtqakJDQwN6enrQ1taGiooKVFVVuUsKZfWCuro611lLdsWBpqamXJJqrV55nI6OTptDnVZWVrb5vvr6eq4fBPCokSSr/2pra9Hc3AyJRIK6ujoAj/pItbS0cHUVEaGqqkquvmqPuro6tLS0nmo4tdaYevw1AwMDdoVVD+m+AXJ5EBMTg9mzZ8POzg5Xr16FmZkZAHC9INvS2NiIsrIylJeXc61/WYKtq6tDXV0dqqurUVlZifz8fFRXV7daIIDWk/rjSR/4XwGVVQCyRoOmpibs7e2hq6sLLS0trqFhamrK/W9kZNRqAm9paYGZmRk++ugj5OTk4N///rdi9xZmmOeIrNHe3k2SpFIpDh48iNWrV0MgEOCbb77B8uXLoayszO2dyhrnsv+rq6tRX18PkUjE7dE2NzdDLBajoaGBS44AUFFRgcbGRrm6oq6uDhKJhItBllRb8/iy2vqMbTXyZXVORUUFlJSUuOQK/K+hoayszF3TbWFhARUVFQiFQujo6AB4VJ/JGjSyuktDQwO6urrQ0dGBpqYmd9SDnQrpg7p8/zqeHT16lNTV1emll16i+vp6vsPhzeHDh0lDQ4OCgoKosrKS73AYpl+JjY0lb29vEgqFFBoaSqWlpXyH1CP+8Y9/kLGxMYlEIr5DYXpZv2t2ERE2bNiAkJAQvP/++/j111+f68u93njjDURHRyMlJQVeXl5ITU3lOySG6fPy8/OxYMEC+Pj4QEtLC7du3cLOnTsH7CHj5cuXo7a2FocPH+Y7FKaX9ask39LSgtDQUPzrX//Czp07sXXrVnZ4CIC7uzuuX78OQ0ND+Pj44MyZM3yHxDB9klgsxpYtW+Ds7IyYmBgcPXoU4eHhcHNz4zu0HmViYoLXX38d33zzDah/d8NiOovvQwmKamxspNdee43U1NToxIkTfIfTJzU0NNDChQtJWVmZNm/ezHc4DNOn/P7772RnZ0daWlq0fv16EovFfIfUq5KSkkggENClS5f4DoXpRf2id319fT1eeeUVREdH4+TJk5g8eTLfIfVpu3btwnvvvYdXX30Ve/fu7XfDMDJMd0pNTcWqVatw6dIlzJkzB1999RVvdwTj26RJk6Curo5z587xHQrTS/r8se7KykoEBAQgISEB165dYwleAaGhoTh79izOnz8PPz8/5Obm8h0Sw/S6iooKrFy5Em5ubigrK0NkZCSOHTv23CZ4AFi1ahUuXLiAlJQUvkNhekmfTvKFhYWYMGEC8vPzERkZCXd3d75D6jcCAwPx119/QSwWw8fHB3FxcXyHxDC9orm5Gbt27cKQIUNw/PhxbN++HXFxcRg7dizfofEuODgYQ4YMwbZt2/gOheklfTbJZ2VlYcyYMZBIJIiKioKTkxPfIfU7Dg4OiI2NxahRozBu3Djs37+f75AYpkdFRETAw8MDy5cvx5tvvom0tDSEhoayDrr/P4FAgPfffx8///wzNzQ4M7D1yV9+Tk4OJk2aBD09PVy7dg2WlpZ8h9Rv6ejo4LfffsMHH3yAhQsXYuXKlXJjfzPMQJCbm4sFCxbgxRdfhKmpKRITE7Ft2zZukBfmfxYuXAhNTU3s3LmT71CYXtDnknxubi4mTZoEHR0dXLlyZcBet9qblJWVsXnzZhw5cgS7d+9GcHBwuyNoMUx/IRKJsGHDBjg5OSE2NhZnz57F5cuX4erqyndofZampiYWL16MHTt2sAb/c6BP9a7Py8vDhAkToKKiIjdMLdN9EhISMHv2bKipqeH06dNwcXHhOySG6TQiwq+//oqPPvoINTU1WLt2LVatWgVVVVW+Q+sXMjIy4OjoiNOnTyM4OJjvcJge1Gf25IuLixEQEMASfA9jA+cw/V18fDz8/f0REhKC8ePHIy0tDWvWrGEJvhPs7e0xceJE7N69m+9QmB7WJ5J8cXExJk2aBKlUioiICJbge5iFhQX+/PNPvPzyy3jppZewZcsWvkNimA4VFhZi2bJl8PLygkQiQUxMDA4cOABTU1O+Q+uXli5dirCwMOTk5PAdCtODeE/yFRUVmDRpEogIf/zxB8zNzfkO6bmgpqaGffv2Yfv27Vi3bh3efPPNVu+9zTB8k0gk2LZtG5ydnREWFoZ9+/bhxo0b8Pb25ju0fu2ll17CoEGD2FU3Axyv5+TFYjECAgKQlZWFmJgYWFlZ8RXKc+3ixYsICQmBvb09Tp06xb4Hps+4cuUKVq5ciczMTKxYsQLr1q2DtrY232ENGKtXr8axY8eQkZHBLjMcoHj7VltaWjB//nykpKTg0qVLLLHw6MmBc2JjY/kOiXnOpaenIzg4GFOmTIG9vT2Sk5OxefNmluC72aJFi5CdnY1r167xHQrTQ3hL8qtWrcK5c+dYD+8+4vGBc8aPH4+ffvqJ75CY51BVVRXWrl2L4cOHIyMjAxcvXsSZM2dgZ2fHd2gDkouLCzw9PfHzzz/zHQrTQ3hJ8p9//jm2b9+OgwcPws/Pj48QmFY8PnDOokWL2MA5TK+RSqU4cOAAnJ2dsXv3bmzZsgV37txBQEAA36ENePPnz8evv/4KkUjEdyhMT+jt2979/PPPJBAI6Ntvv+3tVTOdcOTIEdLQ0KDAwECqrKzkOxxmAIuNjSVvb28SCoUUGhpKpaWlfIf0XCkuLiYVFRU6cuQI36EwPaBX9+QvX76MRYsW4R//+AeWL1/em6tmOikkJATR0dFITU2Fl5cXUlNT+Q6JGWDy8/OxYMEC+Pj4QEtLC7du3cLOnTvZKJe9zMTEBAEBAeyQ/QDVa0k+MzMTb7zxBl577TX885//7K3VMs+ADZzD9ASxWIwtW7bA2dkZMTExOHr0KMLDw+Hm5sZ3aM+t+fPn49KlSygtLeU7FKab9UqSF4vFmDNnDiwtLdkIS/0MGziH6U5nzpzB0KFDsXHjRnz44Ye4e/cuXn31Vb7Deu7NmDEDampqOHXqFN+hMN2sV5L822+/jaysLPz222/Q1NTsjVUy3YgNnMM8q9TUVAQFBWHWrFkYNWoUUlJSsGHDBqirq/MdGoNHN60JDAzEiRMn+A6F6WY9nuT/7//+D8ePH8fhw4dhb2/f06tjelBoaCjOnj2L8+fPw8/PD7m5uXyHxPRxFRUVWLlyJdzc3FBWVobIyEgcO3YM1tbWfIfGPOGVV15BREQEKioq+A6F6UY9muTDw8PxySefYOvWrQgMDOzJVTG9hA2cwyiiubkZu3btwpAhQ3D8+HFs374dcXFxGDt2LN+hMW2YMWMGlJWVcfr0ab5DYbpRjw1rm5WVhVGjRiEgIACHDx/uiVUwPKqtrcW8efNw8eJF7NixAwsXLuQ7JKaPiIiIwAcffIC0tDS888472LhxI3R1dfkOi1HArFmz0NzcjLCwML5DYbpJj+zJNzc3480334SFhQX27NnTE6tgeMYGzmGelJubiwULFuDFF1+Eqakpbt++jW3btrEE34+8/PLLCA8PR21tLd+hMN2kR5L8hg0bkJCQgIMHD7KOdgOYsrIyNm/ejCNHjmD37t2YPn06qqqq+A6L6WUikQgbNmyAk5MTYmNjcfbsWVy+fJkNV90PTZs2DRKJBOHh4XyHwnSTbj9cHxkZiYkTJ+L777/HsmXLunPRTB+WkJCA2bNnQ01Njd2P4DlBRPj111/x0UcfoaamBmvXrsWqVaugqqrKd2jMM/D29saIESOwa9cuvkNhukG37snX1NRg3rx5mDlzJkvwz5nHB87x9vZmA+cMcPHx8fD390dISAjGjx+PtLQ0rFmzhiX4AWDatGk4e/Yseqi7FtPLujXJf/jhhxCJRNixY0d3LpbpJ2QD57zyyits4JwBqrCwEMuWLYOXlxckEgliYmJw4MABmJqa8h0a002mT5+OwsJCJCUl8R0K0w26LcmHh4dj7969+OGHH2BiYtJdi2X6GTZwzsAkkUiwbds2ODs7IywsDPv27cONGzfg7e3Nd2hMN/P09IS5uTnrYT9AdMs5+draWgwbNgy+vr745ZdfuiMuZgC4ePEi3njjDdja2uLUqVNsAJR+6sqVK1i5ciUyMzOxYsUKrFu3Dtra2nyHxfSghQsXIjMzE3/++SffoTDPqFv25D/77DPU19fju+++647FMQNEYGAg4uLi0NjYCF9fXzZwTj+Tnp6O4OBgTJkyBfb29khOTsbmzZtZgn8OTJw4EbGxsewo3ADwzEk+KSkJ3333HbZu3cpuEck8xcHBATdu3MCoUaMwfvx4/PTTT23Ou2vXLtTU1PRecEyrqqqqsHbtWri5uSEjIwMXL17EmTNnYGdnx3doTC+ZOHEiGhsbWcN8AHimJC+VSvG3v/0NPj4+WLRoUXfFxAwwTw6cs2zZMjQ3N8vNc+bMGfztb3/D+vXreYqSkUqlOHDgAJydnbF7925s3boVd+7cQUBAAN+hMb3M2toatra27HD9QEDPYM+ePaSiokJ37tx5lsUwz5EjR46QpqYmBQYGUmVlJRERJScnk6amJgkEAlJSUqKEhASeoxxYtm3bRr/99lu788TGxpK3tzcJhUIKDQ2l0tLSXoqO6aveeustmjhxIt9hMM+oy0m+traWzM3Nafny5d0ZD/McuHXrFllbW5OjoyNdv36dbGxsSEVFhQCQUCgkd3d3amlp4TvMAeHMmTOkpKREVlZWJBaLn3o9Ly+P5s+fTwKBgCZNmkRJSUk8RMn0RT/++COpq6u3+rth+o8uH67funUr6urqsG7dum48rsA8D2QD5xgYGGDy5MkoKCiARCIB8Oi+B7dv38aPP/7Ic5T9X2JiIl577TUQEQoKCvDvf/+be00sFmPLli1wdnZGTEwMjh49ivDwcLi5ufEYMdOXTJgwAQ0NDbh58ybfoTDPoEtJXlZhrFu3jg2CwXSJhYUFRo0aBbFYzCV4GalUig8//BClpaU8Rdf/FRQUYOrUqZBIJCAitLS04PPPP0dOTg7OnDmDoUOHYuPGjfjwww9x9+5dvPrqq3yHzPQxdnZ2MDMzQ1xcHN+hMM+gS0n+yy+/hIGBAVasWNHd8TDPif3792P79u2QSqWtvi4Wi7FmzZpejmpgEIlEmD59OsrLy+U6OEqlUrz22muYNWsW/P39kZ6ejg0bNkBdXZ3HaJm+zNPTk+3J93OdTvKFhYXYu3cvPvnkE1Y5MF0SExODpUuXtjuPRCLBTz/9hOjo6F6KamCQSqV4/fXXcffu3aeOkEgkEsTFxWHnzp3Yv38/LCwseIqS6S9GjRrFknw/1+kk/69//QtGRkZ4++23eyIeZoArKSnBSy+99FQCao2ysjKWLl361OV2TNs++OADnD9/vs1tpqysjG+//bbNIygM8zhPT088ePAAFRUVfIfCdFGnknxhYSH27NmDTz75BGpqaj0VEzOAmZiYICIiAuvXr8fgwYMBACoqKq3O29zcjPT0dPz3v//tzRD7rV27duHbb79FS0tLm/M0NzcjOTkZ+/fv78XImP5q9OjRICIkJCTwHQrTRZ0au37dunXYvXs3srOz2aF65plJpVLuLmaHDh2CWCyGkpLSU0lKQ0MD6enpsLS05CnSvi8sLAwzZ85UaA9dIBDA0NAQGRkZ0NXV7YXomP7MysoKy5cvZ31k+imF9+TFYjF27tyJ9957jyV4plsoKSnBz88Pu3btQllZGY4ePYopU6ZAWVkZKioqEAgEAB7tfa5cuZLnaPuuxMTEDnvHKykpcfd6JyIoKyvj3LlzvREe08+5ubkhOTmZ7zCYLlJ4T37nzp1YuXIlsrOz2WVzTI/Kz8/H4cOH8eOPPyItLQ0CgQBEhAsXLiAwMBBVVVWQ/Wzr6uq48/sNDQ1P3VCjvr4eTU1NHa5TIpGgrq5OofiUlJSgp6en0Lw6OjoQCoVyzxkYGHD/a2trc6cr1NXVoaGhodByZQoKCuDh4YHS0lJuL162vubmZggEAgwePBje3t7w9PSEu7s73N3dWRlmFLZq1SrExMSwcez7KYWT/MiRI+Hh4cEGKWEA/C8pVlZWor6+HvX19airq+OSam1tLZqbm1FTU4OWlhZUV1dDKpVyCbqyshJSqRTV1dVoaWlBTU0NmpubUVtbCwBobGyESCRCS0sLGhoaFErUA4lQKISOjg43LWsYKCsrQ1dXF0KhEJqamoiPj+du6iMQCGBgYIDBgwfD0tIS1tbWeOGFF6Cnp8c1JnR0dKCqqgo9PT3o6OhAW1sbWlpaCjdamOfPjh07sHbtWlRVVfEdCtMFCiX5W7duwdPTE1FRURg7dmxvxMX0kIaGBlRVVaGysvKpv48nbVnCbmu6o6SrqakJNTU1ueQiFAqhq6sLZWVl6OnpQUlJCQYGBhAIBNDX15fbQ1ZRUZG7pamenh7u3LkDbW1tjBw5kttb1dDQ4E4fPfke4Olk2R5ZTB0Ri8VoaGjocD5Z4+VxsoaNjKzxA8gfdZA1cp5cjqxx1dTUhOjoaFRXV0NHRwcaGhpQUVFBQ0MDGhsb22xstUdHRwdaWlrQ1tbmGgGtTevo6EBfXx/6+vowMDDg/sr+l51mYQaGq1evYtKkSSgqKmJHgPohhZL8ihUrcOHCBdy7d48V4D5CJBKhrKwMRUVFKCsr4x4VFRWoqqpqM5G3lpxUVFRgYGAAbW1t7q+Wlha0tLRgYGDA/a+trQ19fX1uWkdHB3p6etxrsmTD+mz0XdXV1Vzyr66u5o7C1NTUoKamhpuurq5GbW0t17iTTdfX16O2thZVVVVyjZXHyRoATzYCHm8MGBsbw8TEBMbGxjAyMoKRkRGUlZV7eWswiigoKMDgwYPx559/Yty4cXyHw3SSsKMZmpqacOTIEfy///f/WILvQRKJBMXFxcjLy0NpaSnKyspQUlKCkpISLoGXlpaiuLgYZWVl3F6ejIaGBoyMjOT2qExMTODk5MRNt1Xpamlp8fSpmd4mO1JibGz8zMuSnXaRNSBba1TK/hYUFCAlJYWbr7S0FE/uXxgZGcklfVNTU7lpExMTmJmZwdzcHIMGDXrm+BnFWFhYQFdXF/fu3WNJvh/qMMn//vvvqKiowLx583ojngFJVskVFha2+Tc7O1vucKq6ujqXsC0sLGBubg4HB4ennpNNm5ubs0YY06tkl+IZGhp26f1isZgrA7Lk//h0Xl4e7ty5g8LCQuTm5soNoKSmpgZDQ0OuHLT1l5WL7uHg4ICHDx/yHQbTBR0erp8+fTqkUinOnz/fWzH1KzU1NcjKykJ2djYyMzORnZ3NPfLy8lBcXCy3x2JoaChXCQ0ePPipaRMTk073smaYga68vBxFRUVcQ+DJhnJ+fj4KCwvR2NjIvUddXR2DBw+GtbU1bGxsYGtrC1tbW9jY2MDGxgaWlpZtDsbE/M+MGTNgYGCAAwcO8B0K00nt7skXFxfj0qVLOHToUG/F0+c0NDTg/v37ePDgAbKzs5GVlcUl9ezsbFRWVnLzGhsbc5XH+PHjYWVlBQsLC7m9C3a+mmG6ZtCgQRg0aBCGDh3a7nzl5eVc0i8qKkJeXh5ycnKQlZWF69evIzs7m+uboqysjMGDBz/VALCzs4OjoyOsrKx646P1eebm5sjMzOQ7DKYL2k3yx44dg6amJmbOnNlb8fBGds4wIyMDGRkZSE5ORkpKitxhdAMDA9jb28Pe3h7jxo2DhYUFN+3o6MhGD2OYPkDWGBg2bFib81RWVnJlXXYkICMjA+fOncODBw+4ToWqqqqwtLSEq6srhg4dypV32eN5YW5ujpiYGL7DYLqg3SR/6tQpTJs2bcDsfUokEty7dw9JSUlISkpCWloa7t27h4yMDO7SJVlnNScnJ/j5+cHR0ZF7DJTtwDDPOwMDA3h6esLT07PV10tKSpCeno709HTcv38f6enpOHfuHO7fv88dBRg0aBAcHR0xZMgQDB06FCNGjMDw4cNhZmbWmx+lV5ibm6OwsJDvMJguaPOcfFVVFUxMTHDgwAGEhIT0dlzPrKSkBLdv3+YS+p07d5CcnIympiaoqqrCxcUFrq6uXEJ3dHSEk5MTGxSEYZg2SaVS5OTk4P79+1zyT0tLw927d5Gfnw/g0Wk7WcKXPVxdXfv1Tb1Onz6N2bNnQyQSsf5C/UybSf7YsWOYO3cuysrK+nziq6qqwo0bN7hHYmIiiouLATxqgQ4fPpwrdG5ubnBxcWGdbRiG6Vbl5eW4ffs27ty5g6SkJNy+fRvJycloaGiAUCjEkCFD4OnpCV9fX/j6+mLYsGH9ZmyAuLg4eHt7IysrCzY2NnyHw3RCm0n+nXfeQWJiIq5fv97bMbVLKpUiJSUF169fx/Xr13Hjxg2kpaWBiODg4AAfHx94eHhwid3IyIjvkBmGeU61tLTg/v37XNKPi4tDXFwcampqoK2tjdGjR2PMmDHw9vaGj49Pt4xf0BPS0tLg4uKCO3futNvXgel72kzyzs7OePnll/HFF1/0dkxPSU1NxYULF3Dp0iXExMSgpqYGWlpaGD16NHx9feHjBfa0hQAAIABJREFU4wMfHx+YmJjwHSrDMEy7WlpantpRuXfvHogIjo6OmDBhAoKCgvDiiy/2maOo2dnZsLW1RWxsLLy8vPgOh+mEVpN8WVkZjI2Ncf78eQQFBfV6UNXV1QgPD8fFixdx4cIF5OTkwNDQEJMnT8b48ePh6+sLNze3p+7uxTAM0x9VVFTgxo0buH79Oi5fvoybN29CIBDA19cXgYGBCAoKgru7u0L3VugJpaWlMDExwdWrVzFhwgReYmC6ptUkf/HiRQQFBaGkpKTXDh/l5eXh+PHjOHXqFGJiYiCVSjF69GjuB+7l5dVvzl8xDMM8i/Lycly+fBkXL17ExYsXUVhYCBMTEwQFBeHVV19FQEAAVFVVey2e+vp6aGtrIywsDNOmTeu19TLPrtVd4fj4eNjY2PR4gq+trcWRI0dw8OBBREdHQ0dHBzNnzsQ777yDKVOmsPGpu6C2tlbhu64xDNM3DRo0CCEhIQgJCQERISkpCRcvXsTvv/+OmTNnQl9fHy+//DLefvttjBkzpsfj0dDQgEAgeOqeGUzf1+qxn7S0tB7tXJGcnIzQ0FAMHjwYK1euhJWVFU6ePIni4mLukr2+luCvXLmCJUuWQCAQQCAQIDAwsE+NBLhz506MHz8eLi4uPb6uU6dOwcrKCqmpqT2+rs4gInzzzTfYvHkzHB0dMX/+/A5vrzqQ9NXvpS3Nzc2IjIzEp59+iosXL3LPd+Vz9LfP3hkCgQAjRozA6tWrERUVhezsbHz66aeIj4/H2LFjMWzYMHz33Xeor6/vsRiUlJSgpqbGknw/1GqSz8jI6JHRnOLj4zF79my4ubkhKioKGzduREFBAQ4dOoSZM2f26etIJ0+ejD179nBHN3788UfMnTu3U8voycEklixZAqlU2iNJ7cm4tbS0YGJi0ucGB/r8889x7949rF27Fvv27UN1dbXcTU0Gmv7yvbTlr7/+wr59+/DFF18gLy+Pe74rn6O/ffZnYWVlhQ8//BAJCQmIi4uDr68v1qxZA1tbW2zatAm1tbU9sl6hUIjm5uYeWTbTg6gVFhYW9PXXX7f2UpcUFxfT4sWLSUlJiby9venUqVPU0tLSbcvvTS+88AIBoNra2k69r6KigiZNmtRDUT0SEhJCZmZm3brM3oi7u5iYmNCXX37Jdxi9oj99L+25desWAaA9e/bwHUq/VlJSQp9++inp6+uThYUFHThwgKRSabeuQ0VFhQ4dOtSty2R6Xqt78qWlpTA1Ne2WRsTVq1cxcuRInD9/Hvv27cP169cxa9Ys3nqJPivZbSs7c/tKkUiEkJAQZGRk9FRYPaI/xd3Q0ICSkpLn4rai/el76Uhvdh4byIyNjbFp0yZkZGRgzpw5WLRoEWbNmoWKiopuWb5UKoVEIunTR1uZ1j2Vaevq6iCRSGBgYPDMC9+7dy8mT56MiRMnIi0tDQsWLBhwlXBiYiL+/ve/w97eHvX19ViyZAmMjIzg5eXFVcInT55EamoqysrKsHTpUnz11VfIz8/H5s2bMWzYMFRUVCAwMBA2NjYoLy8HAJw4cQLLly/HRx99hKlTp2LdunVyt9AEHg01GRoaijVr1uD999+XO3x75MgR6OrqcnfRqq6uxsaNG6GsrAxfX1+55Zw7dw7vvvsuVq5cCV9fX+zevbvNuCsrK7F3715MmTIFp06dkltOezErsp3a096y9+/fj6VLlwLA/8fencfVmP7/A3/VaaEdbRQpu6ZCVJIloRpLskZlz5ZlzJgha3YZxgczRjFDIZEsGZIkS8neprQi2hftqtNy/f7wPecna8s55z7ndD0fjx51Tt339TpL530v131d8Pf3h4uLCzw8PL65vuvXr8PFxQVr1qzBokWL8Pvvv2PcuHFNeu5yc3Ph4uKCbdu2wcXFBfb29tzX71uv77eWa+776UuvS2Ofc0IIjhw5giVLlsDU1BRjxoxBSkrKV5+7hIQErF+/Hn379kVWVhYmTpyI9u3bw8TEBA8ePGj06/Yl33p/fe19+rVlvve4oqOjMXfuXHh4eMDOzg6jR4/+ai5R0a5dOxw4cAB37txBVFQUzMzMeHKakPOa0SIvgj7dtc/OziYAyJ07d1p0iCAoKIiwWCyyefPmFq1H2HTv3p0AIOXl5YSQD8/XqFGjCADi6upK4uPjSVRUFJGVlSUODg7c5caNG0e6du3KvR0UFER69+7NfY68vLyIiYkJyczMJPv37yfm5uaEzWYTQggpKCggPXr0IMOHD+cegjt9+jQxNTUllZWVhBBC8vPziaqqaoPD9WPGjCHa2toN8hsYGBAzMzPubR8fH+Lg4MA9fbJjxw4CgISGhn4xd0JCAlm1ahUBQM6fP8+9/3uZG/s8fUljno+CggICgGzfvv2b6yKEEG9vb2JiYsJ9Devr60mfPn2IiopKk567ESNGkOnTp3NvGxkZEScnJ0LIt1/fby3X3PfTl16Xxq5r165d5MSJE4QQQmpra0nfvn2JpqYmqaio+OLzt3btWqKiokJYLBZZtWoVCQsLIwEBAURVVZXIycmRrKwsQkjjXrfnz583OFz/tffXt96nX1vme4+rZ8+eJDw8nBBCyPv374mFhcUXH6+oysnJIX369CGGhoakurq6ResqKioiAEhISAiP0lGC8lmRz8rKIgDIvXv3mr3Smpoa0rNnz+9+eIuiT4s8IYS4ubkRAKSgoIB7n4WFBenRowf39qcfyoQQMn/+fAKApKSkcO/Lzc0l8vLyxMfHp8HfHj9+nAAgJ0+eJBUVFaRjx47E19e3wd/Y29s3KPITJ078rFCZmZlxC1VeXh5RVlYmL1++5P4+Pz+fTJo0iSQkJHw19+3btxt8oDYmc2Ofp081dt2NLfLFxcVEVVWVBAQENLjfwcGhQZH/3nNHCCGWlpZk586d3NuOjo7E0NCQe/tLr29jlmvu++nT16Ux68rMzCQaGhoN+shs2rSJACB+fn7ka2bOnEmkpaW5BZwQQvz9/QkAsmnTpka/bp8W+S89jsa8Tz9d5nuPi81mEwkJCXLgwAHu7y9evPjVxyuqXr9+TeTl5VvcxyonJ4cAIHfv3uVRMkpQPrtOnnM4nXx5tNtGiYuLQ3JyMi5evNjsdYgSziA9H4/Ap62tjdTU1G8uJy0tDSkpKXTv3p1734MHD1BRUYEuXbo0+FvOoeSwsDCoqakhOzsbBgYGDf6mqYfSwsPDUV9fD11dXe59qqqqCAgI+OZyn4402JjMTk5OzXqeGrvuxrpx4wYKCgowYMCAbz6mxrh16xaAD/0BTp8+jUePHjX4v/nS69uY5Zr7fvrSY/jeuu7fv4+amhosWrSowXILFiz45mxjcnJyYLFYDSZ6mjhxImRlZREXF9ei1+3Tx9GY9+mny3zvcUlLS8Pa2ho//fQTnj9/jt27d2PixIlffbyiSkdHB3PmzEFAQAB+/vnnZq+nsrISAFrF1Qvi5rNPBc6LyHlRmyMrKwsA6GxFzZCeng4An3WYUVVVhZycHLKyspCYmAig5Z2Wnj9/jpqaGhBCWtRXojGZhWXdCQkJAMCT6TLr6uqwZ88ePHnyBCtWrICpqeln56R5uRw/vHjxAvLy8tzz2y0hJSWFTp06oba2lqevW3Pep415XAEBAXBxccHRo0dx8eJFnDt3DpaWlo3OJSp0dHTw33//tWgdnPP6mpqavIhECdBnHe+UlJTAYrFQVFTU7JVyBtK5fft2s9fRWnH2Vr7WGa13797c4s75IG0uJSUlVFVVcQvfx77VOepTjcncXLxeN2fP9lsdyxqjvr4eP/74IxISEhAQEIDhw4fzdTl+kZOTQ0ZGRoPr1Dny8/ObvL7379+jd+/ePH3dmvM+bczjkpKSwunTp3H69GlISUnBxsZGLAfTCQsL++yoX1NlZmZCUlKSFnkR9FmRl5SUhJKSUouKvI6ODiZNmoTVq1e3aD3CiHNYtamnMyQlJVFeXv7dvxs8eDCUlJQ+61mckZGB9+/fY8KECTA0NAQAnD17tsHffDoYjpSUFMrLyxvcV15ejvr6egDAoEGDAAAbNmzg3gd8GLTo6tWrjc7dmMzN1dh1N/b14IwIeObMmQb3l5aWNrj9vefu0aNHuHHjRoPJOjh7m9/S3OU+1dj30/cYGBiAEII1a9Y0uD8tLQ2HDx9u0rqys7ORn5+PKVOm8PQ90Zj36ae+97iqq6vh5eUFAJg5cyYePHgAQgjCwsIanUsU+Pv7IygoqEWH6oEPRV5dXb3B6RlKNHzxRKSWltYXt4Cb4sCBAxgyZAhsbW1x6dIlsdkC5BSDkpISKCgocH8G0GA0qLy8vAZDQHbq1AkFBQV4+vQpysrKYGJiwi0ixcXFUFFRAfBhzGoPDw8sXboUoaGhsLKyAgAcPHgQs2fP5h5OtLS0xIkTJ2BsbIzZs2cjPj4e4eHhyM/Px5kzZ2BnZwcDAwOcP38eu3btwrRp03Du3DlUV1fj7du3iIqKgrm5Off1sbKywpQpU5Ceno53797h2LFjX83NOXTH2SNqbObGPE+fauy6Oe/X7w27OWHCBHTt2hVeXl7o27cvRowYgcjISMTExDT4u+89d5zDxt7e3jAxMcHjx48RHx+P3NxcxMbGQkND44uvb2OWa+776dPXpTHP+ejRozFo0CD4+vqiqqoK9vb2KC0txYULF+Dn5/fN57K6uhoxMTEwMjICAGzfvh2zZ8/mTkXamNeN8//08ZCsnz6OxrxPP12mMY/r33//xZIlS8BisdCpUycoKyt/1ldDlPn7+8PJyQkrV65s8WmIrKwsaGlp8SgZJVBf6o03fvx44ujo2OJefcnJyaR79+5EU1OTXLt2rcXrY1JYWBhZunQpAUAAEFtbW+Ln50du3rxJunbtSgCQpUuXkry8POLj40MUFBQIAOLu7k5qa2tJTEwM0dbWJj179iT+/v7Ey8uLqKmpEQDE2dmZPHv2rEF7ly5dImPGjCHLli0jGzduJPv27WswglVJSQmZO3cu0dDQIF26dCHu7u5k4cKFZO7cueTmzZukrq6OlJSUkPHjxxMFBQViZmZGHj9+TObMmUOcnJxIYGAgIYSQiooKsmTJEqKlpUU0NDTIkiVLSHFxMbedT3OHhoaSYcOGEQBk4MCB5MaNG43K3Njn6Wu+te6nT5+SGTNmEABEV1eXnD59usFj+FRycjIZOnQoUVZWJkOHDiXXr18nTk5ODXrXN+a5W7x4MVFUVCRmZmbk5s2b5Nq1a0RVVZVMmTKF7N+//6uv77eWu3z5crPeT196XRr7nBcWFhJHR0eirq5O1NTUyKxZs0hmZubX/xkIIQsWLCAyMjJk1apVZOrUqWT+/Plk27Ztn42y9q3X7eHDh8TW1pYAIAMGDCBXr1796vvrW+/Try3zrcdVVVVFBg0aRKytrcnu3bvJwoULydGjR7/5mEVFWVkZWbFiBZGQkCDLly/nych3jo6OZPz48TxIRwnaF4v8Tz/9RAYNGsSTBkpKSoiDgwMBQOzs7EhSUhJP1ktRvPRpkae+bcGCBaRNmzZMx6A+UltbS3x8fIiWlhZRUVHh6RC0lpaWZNGiRTxbHyU4Xxxbtl+/foiNjeXJ5B5KSko4c+YMQkJCkJKSgr59+2L69Ol49uxZi9dNURTV2lVVVeHYsWPo06cP5s6di7FjxyI5ORkzZ87kWRuJiYmfXQpKiYYvFnljY2NUV1d/sTdrc40aNQqxsbE4ffo0kpOTYWxsDDMzM/zzzz98nSKRohrj/fv3YLPZLRofojXhDH9Nny/mJCQkYNWqVdDS0oKrqyuGDx+OFy9ewNPTkztbJi/k5eUhOzsb/fv359k6KcH5YpHv06cP5OXl8fTpU542xmKxuHvxYWFh0NPTg6urKzQ0NDBjxgxcvny5SZduUVRLZWVlwc3NDdevX8f79++/O7Y6Bfz9998ICQlBXV0dFi5ciPDwcKYjtRpv3rzB3r17YWxsDH19fVy+fBm//PILXr9+jaNHj6JHjx48b5Nz1JXTwZISLRLkK5vi5ubm6N+/P/766y++BigsLMS5c+fg5+eH8PBwyMvLw8rKCtbW1rC2tm4wyhVFUVRrUlNTg/v37yM4OBjBwcGIioqCiooKJk2ahBkzZsDS0pLvM3ru3r0bhw8fxps3b/jaDsUfXy3yy5cvx+PHjwU6EldGRgYuX76M4OBg3Lp1CxUVFejZsydsbGxgbW2NESNGQE5OTmB5KIqiBO3Vq1fcoh4aGoqysjJ0794d1tbWsLW1xejRowU6Re+MGTNQUVGBwMBAgbVJ8c5Xi7yvry/mzp2LwsJC7vXggsRmsxEeHs59s8fGxkJWVhaDBg2CmZkZzM3NYWZmJjbX31MU1frU1dUhLi4O9+/fx4MHD3D//n2kpaVBQUEBlpaWsLa2ho2NDbp168ZYxj59+mDq1KnYunUrYxmo5vtqkc/NzUXHjh1x9epV2NraCjrXZ7Kzs3Hjxg2Eh4fjwYMHSEhIQH19Pbp27cot+GZmZujXrx8dlYmiKKFUUFCAyMhIbkF/8uQJysvLoaSkBBMTEwwePBiWlpYYMmSIQPfWv5VXQ0MDFy5cgJ2dHdNxqGb4apEHPoz6ZW1tjb179woyU6OUlpbi4cOH3H+YBw8eoKioCG3btoWBgQGMjIxgYGAAQ0NDGBoaol27dkxHpiiqlaivr8fLly8RExOD2NhYxMXFISYmBi9fvoSEhAR69+4NMzMzDB48GGZmZtDX1+f7ufXm8PX1xZw5c5Cfnw9lZWWm41DN8M0iv2rVKoSFhSE6OlqQmZqFEILExEQ8fPgQ0dHR3H+qwsJCAEDnzp25BZ+zAdCzZ89mTTFKURTFUVRUhLi4OMTFxSE2NhYxMTF4/vw5KioqwGKx0L17d+5nD+fSYVHZ6Zg1axYyMjK40yNTouebRf6///7DhAkTkJ2dDQ0NDUHm4pmioiLEx8fj6dOnePr0KRISEhAfH4+qqipISUmhS5cu0NPTg56eHvr27Qt9fX3o6emha9euQrllTVGU4LHZbGRkZODly5eIj49HQkICXr58iZcvX+LVq1cghEBZWRk//PAD9PX10bdvXxgbG6N///6Ql5dnOn6z1NfXo1OnTvj555/x22+/MR2HaqZvFvmysjJ06NABPj4+cHBwEGQuvmKz2UhISEBiYiKSk5ORlJSElJQUJCcncyf0kJeXR8+ePdGjRw/ud11dXejo6KBTp070CABFiZny8nKkp6fj9evXePXqFZKSkpCcnIzk5GS8efMG9fX1kJSUROfOnbmfCb169UKvXr2gr68PbW1tph8CTz158gSDBg1CTEwMd+ZLSvR8s8gDwLBhw6Cnp4cTJ04IKBKz8vLyGvxzp6SkICkpCWlpadxBUqSkpKCtrQ0dHR3o6Oiga9eu3J91dHTQuXNnyMrKMvxIKIr62Lt375Cens79ev36dYPbnFN7AKCqqoqePXuiV69e3A19TmFv06YNg49CcLZt24a///4bmZmZ3NkTKdHz3SLv4eGBffv2IScnp1UfviaEIDs7+4sfDpz7ONN3SkhIoGPHjtDR0YGGhga0tbWhqanJ/a6lpQVNTU2oqqoy/KgoSvSx2Wzk5uYiIyOD+z0nJ6fB9zdv3qCsrIy7jKamZoMNc84X52gdE5cNC5uBAweif//+OHr0KNNRqBb4bpGPj4/HDz/8gMjISJiZmQkql0jKz89vUPjfvHmDnJwcZGZmIjs7G1lZWaisrOT+vaysLDp27AgtLS107NgRnTp1QseOHbkbAKqqqlBTU4O6ujoUFRUZfGQUJVh1dXUoKCjgfuXl5XG/Pi7eubm5yM3NbbCsmppag41pbW1tdOnSpUExby1748314sUL9O3bF7du3WrxXPQUs75b5AGge/fumDFjBrZt2yaITGKtqKiIW/A5X5zb2dnZyMzMRF5eHveoAEebNm24hV9DQ+OzjQA1NTWoqqqiXbt2UFFRQbt27ejogJTQKCoqQlFREYqLi1FYWIj8/Hzk5+ejoKAA+fn5yMvLa1DU8/PzGywvISEBVVVVqKurc4s3Z+P449uampr0VBkPrFu3Dj4+PkhPTweLxWI6DtUCjSryK1aswL179xAVFSWITBQ+zIrG2YP5+MOwoKAAubm5DT4Qc3NzuR0GPyYrK8st+B8X/6/9rKSkBGVlZcjLy0NeXh5KSkoMPHJK2NTW1qKsrAwlJSWoqKhAeXl5g6LdmJ8/JS0tzd1Q/XjDlbOx+vGGK+eLFhvBqK2tRdeuXTFr1izs3LmT6ThUCzWqyAcHB8PGxgavX7+Gjo6OIHJRTcRms1FYWPjND9yvfQh/fK7yU0pKSpCXl4ecnBxUVFSgoKDA3Qho164d92cFBQWoqKiAxWJBWVkZUlJSUFRUhLS0NBQUFCArKws5OTm0adMGbdu2hZycHN3j4qGSkhLU19ejuLi4wfeSkhJuka6pqUF5eTkqKyvx/v17lJSUoKysDBUVFaioqEBxcTHKy8u5t4uKilBRUQE2m/3FNqWkpLgbiI3ZmOTcbt++vchcJ94anT9/HtOnT0dqaiqdIEwMNKrIV1dXQ01NDXv27MHixYsFkYsSoLq6OhQVFaG0tJS7t1ZRUYHS0lKUlpZybzemKHAKSWPJy8tDRkYGCgoKkJaWhpKSEnePTVFRkXupYtu2bbnnUTkbDsCHw7gqKirc9X28/MeUlZW/23GUxWJ99+gFp2B+T2VlJaqqqj67/+OiWV1dzT0twynIHJyi/eky79+/R3V1NXc+97KyMtTW1n43z8ePT1ZWFvLy8o3aaOPcVlRU5B7lUVBQQLt27Wg/ETFlaWkJRUVFOiGNmGhUkQeASZMmobq6GlevXuV3JkoMsNlsVFRUcAsZp+h9rUiVlpairq4OxcXF4LwlP/6Z8/dAw+JYV1eH0tJSbrsfL8NRU1OD4uJi1NbWfrPDFSfb93xpg6GqqgrS0tLcDYyvbTBwjmQA4B7t4Ph47/bjDZyPl/n0SAhnI4nz95yNHE7Gdu3afbYhRFFfExUVBWNjYwQFBcHa2prpOBQPNLrIe3t7Y9GiRcjLy6PnaimRY2VlBTabjXv37vFl/YaGhrCwsMDhw4f5sn6KEoRp06YhOTkZUVFR9Np4MdHoC9/Hjx+Puro6XL9+nZ95KIrnbty4gVu3bmH79u18a8PJyQl+fn6NOhJAUcIoLS0NFy5cwLp162iBFyON3pMHPuwNaWhowNfXl5+ZKIpnCCEYPHgwVFVV8d9///GtnaysLHTp0gX+/v6wt7fnWzsUxS/z58/HnTt3kJSURK9kECNNGsLOzs4O165d+2pvW4oSNgEBAXj06BHfx3jo1KkTLC0tcfLkSb62Q1H8kJGRgVOnTsHNzY0WeDHTpCI/adIklJaWIiwsjF95KIpn6urqsGnTJjg4OKB///58b8/Z2RlXr15tMAY6RYmCvXv3Ql1dHc7OzkxHoXisSUVeW1sb/fv3x6VLl/iVh6J4xtvbGykpKdi8ebNA2ps8eTJkZGRw7tw5gbRHUbyQk5ODo0eP4pdffoGMjAzTcSgea/KMMxMnTsTFixe51/BSlDBis9nYvn075s+fj169egmkTXl5edjZ2dFD9pRI2bVrF5SVlbFw4UKmo1B80Kwin5ubi0ePHvEjD0XxxF9//YXs7Gxs2LBBoO06OzsjMjISycnJAm2XoprjzZs38PT0xMaNG+lcF2KqyUXewMAAPXr0wPnz5/mRh6JarLy8HLt378by5cuhra0t0LZHjx4NLS0tegUKJRK2bdsGTU1NzJ8/n+koFJ80a4L4adOm4dy5c/SQPSWU9u3bh+rqaqxZs0bgbUtKSsLBwQHe3t6fjbxHUcIkNTUV3t7ecHd3p+fixVizivz06dPx9u1bREZG8joPRbVIQUEB/vjjD6xevRodOnRgJIOTkxNev36NBw8eMNI+RTWGu7s7dHV14eTkxHQUio+aVeQNDAygr6+Ps2fP8joPRbXIrl27ICsri5UrVzKWoV+/fujVqxcuXLjAWAaK+paYmBicOXMG27Zt486RQImnZhV54MMhe39/f9TV1fEyD0U1W1ZWFv7++29s3LiR8RnS7O3tcf78eXrInhJKv/76KwYOHIipU6cyHYXis2YXeQcHB+Tk5ODOnTu8zENRzbZ582aoq6sLxaVAkydPxuvXrxEdHc10FIpq4MqVKwgJCcG+ffvoGPWtQLOLfM+ePdGvXz96yJ4SCikpKThx4gS2bNkCWVlZpuNg4MCB0NXVRUBAANNRKIqrtrYWbm5umDZtGiwsLJiOQwlAs4s88KED3oULF7jzfFMUUzZs2IBu3brB0dGR6ShcdnZ28Pf3ZzoGRXH9/fffSElJwY4dO5iOQglIi4q8g4MDCgsL6fSzFKNiYmJw/vx57Nq1S6g6EU2ePBnJycl48eIF01EoCsXFxdi6dStWrVqF7t27Mx2HEpAmTTX7JZaWlujQoQMdHIdijI2NDQoLC/Ho0SOhOsdYX18PbW1tLF26VOAj71HUp3766Sf4+voiJSUFysrKTMehBKRFe/IAMHv2bFy5cgUFBQW8yENRTXLv3j0EBwdj9+7dQlXggQ8D40ycOJGel6cY9/z5cxw+fBi7d++mBb6VafGefEVFBTp27IgdO3Zg+fLlvMpFUY0yZMgQSEtL4/bt20xH+aLQ0FCMGjUKqamp6NatG9NxqFaIEAJLS0u8f/8eDx48gKRki/ftKBHS4ldbXl4ekydPhre3Ny/yUFSjXblyBffv38fu3buZjvJVI0aMgKqqKi5evMh0FKqVOnHiBMLDw+Hp6UkLfCvU4j15ALh9+zYsLS0RExMDQ0NDXuSiqG+qr6+HsbEx9PT0hP5w+Lx58/DixQs6DDQlcEVFRejduzccHBxw4MABpuNQDODJZt3w4cPRrVs3+Pj48GJ1FPVdZ86cQVxcHLZt28Z0lO+ugE7xAAAgAElEQVSaNGkSHj58iOzsbKajUK3MunXrICEhgS1btjAdhWIIT4q8hIQEnJ2dcfLkSVRXV/NilRT1VTU1Ndi8eTOcnZ3Rt29fpuN8l5WVFdq0aYPg4GCmo1CtSGRkJLy8vLBv3z6oqKgwHYdiCM9O0CxcuBBFRUV08A+K744dO4a3b99i48aNTEdplLZt22LYsGEICgpiOgrVSlRXV2PBggWwtLTEzJkzmY5DMYhnRb5jx46ws7PD4cOHebVKivpMZWUldu7cicWLF0NPT4/pOI1mY2ODkJAQ1NbWMh2FagW2bt2K9PR0eHl5Cd2lpZRg8bSrpaurKyIjI/H06VNerpaiuA4ePIiioiKsW7eO6ShNYmtri6KiIjx8+JDpKJSYi42Nxe+//w4PDw+R2hCm+IMnves/ZmhoCBMTExw7doyXq6UoFBcXo1u3bli6dKlIdLj7VPfu3eHg4IDt27czHYUSU7W1tTAzM4OMjAzCw8PpJXMUb/fkAWDx4sXw9fXFu3fveL1qqpXbs2cPCCH4+eefmY7SLNbW1vS8PMVXe/fuRXx8PI4fP04LPAWAD0Xe2dkZ0tLSOHHiBK9XTbVieXl5+PPPP+Hm5oZ27doxHadZrKysEB0djcLCQqajUGIoLi4O7u7ucHd3R69evZiOQwkJnh+uB4CVK1fi8uXLSElJgbS0NK9XT7VCy5Ytw4ULF5Camgo5OTmm4zRLcXExVFVV4e/vD3t7e6bjUGKEzWbDxMQE8vLyuHv3LlgsFtORKCHBl+M5q1evRlZWFvz8/PixeqqVef36NY4dO4bNmzeLbIEHABUVFRgZGSEsLIzpKJSY2bRpE1JTU3HixAla4KkG+FLkO3fujOnTp2PXrl2or6/nRxOUGHrz5g2ioqI+u3/Tpk3o0qUL5s2bx0Aq3ho5ciRu3brFdAxKjNy/fx979+7FH3/8gR49ejAdhxIyfOuZsW7dOiQlJeHq1av8aoISMxERETA2Nsa0adOQkpICAIiPj4evry+2bt0qFqd+LC0tkZCQgJycHKajUGKgoqICc+bMgZWVFVxcXJiOQwkhvpyT5xg/fjzy8/Px4MEDfjVBiZEtW7Zg+/btkJCQQF1dHebNm4c3b94gOzsb0dHRYtFbuKysDB06dICPjw8cHByYjkOJuIULF+LChQuIi4tDx44dmY5DCSG+fmr+9ttvePjwIcLDw/nZDCUmEhMTQQhBTU0N6uvr4e3tjZs3b6Jz587Iz89nOh5PKCoqol+/foiIiGA6CiXi/P39cezYMXh5edECT30VX4v80KFDYWFhQWdAoholPj4edXV13NucYh8SEgJdXV2sXbsWpaWlDCbkDXNzczrtLNUib968weLFi7FkyRJMmjSJ6TiUEOPr4XoAuHfvHoYNG4abN2/CysqKn01RIk5eXh7v37//6u+lpKSgpKSEf/75BxMnThRgMt46e/YsnJycUFRUBAUFBabjUCKmtrYWw4YNQ1lZGR49eoS2bdsyHYkSYnw/yTl06FBYW1vDzc0NfN6eoERYdnb2Nws8h56eHoYPHy6ARPwzePBg1NbW0jkeqGbZsGEDYmNjce7cOVrgqe8SSE8mDw8PPH36FJcvXxZEc5QISkpK+ubvpaSkYG5ujlu3bonsiHccXbp0gba2Nj1kTzVZWFgYfv/9dxw8eBB9+vRhOg4lAgRS5I2MjDB58mRs2LChwTlXiuJITk6GlJTUF3/HYrFgY2OD4OBgKCoqCjgZf5iZmdEiTzVJZmYmHBwcMH36dLEYM4ISDIFdk7Rt2zYkJSXh1KlTgmqSEiHJyclfvEROUlISM2fOxMWLF9GmTRsGkvGHiYkJnjx5wnQMSkTU1NTAwcEBKioqOHLkCNNxKBEisCLfq1cvzJ8/H+vWrUNZWZmgmqVERGJiImpqahrcJyEhAVdXV3h7e391L19U9evXD1lZWcjLy2M6CiUCVq9ejejoaFy4cAFKSkpMx6FEiEBHF9m5cyfYbDadT5v6zPPnzz/rmLlp0yYcPHgQEhISDKXin379+gH4MHMYRX3LuXPncPDgQRw+fBj6+vpMx6FEjECLfPv27bFlyxb873//Q2JioiCbpoRYbW0tMjIyuLclJSXh6ekJd3d35kLxmZqaGjQ1NRETE8N0FEqIJScnw8XFBStXroSzszPTcSgRxPfr5D9VX18PU1NTqKioICQkRJBNU0IqOTkZvXr1goSEBFgsFnx9fTF16lSmY/GdjY0NNDQ04O3tzXQUSgiVlZXBzMwMysrKuHPnjljM3UAJnsBPdEpKSuLAgQOwsLDAxYsX6bzaIoQQguLiYpSUlKCkpAQ1NTUoLi7m/p7NZqOiooJ7W1ZWtsHUsAoKCpCWlka7du2gpKQEZWVlSEtLIzk5mfv3ly5dgrW1teAeFIOMjIwQHBzMdAxKCNXX18PZ2Rnv3r1DcHAwLfBUswl8T55j9uzZCAsLQ1xcHJSVlZmIQP2f/Px8pKWlISMjA1lZWcjNzUVmZib3e1FREUpLS/kypKycnBxYLBYqKipgZGSEvn37Ql1dHdra2tDQ0IC2tjZ0dXWhra0tFhPUfOz06dOYN28eKioqxK5jIdUy69evx++//46bN29i2LBhTMehRBhjRb6wsBD6+voYN24cjh07xkSEVoXNZuPFixeIjY1FYmIiUlNTkZqairS0NJSUlAD4cJRFXV0dGhoa0NLS4hbbdu3aQVlZGUpKSty9cCUlJcjIyEBZWZlbfFksVoOev5WVlaiqquLeLikpQV1dHYqLi1FcXMzdcDhz5gx0dXUhISGBzMxM5OXlISMjA3l5eWCz2QA+7OXr6emhe/fu6NatG3r27AlDQ0MYGBiIbG/jp0+fYuDAgUhOTqbzgFNcAQEBmDp1Kjw9Pen0sVSLMVbkASAwMBATJ07E1atXYWtry1QMsVNdXY0nT54gMjISMTExiI2NxYsXL1BTUwNZWVn06tWLWyy7d+/O/VlLS4uRPUo2mw0ZGZkv/i43NxcvX75ESkoKd6MkNTUVSUlJKCkpgYSEBHR1dWFkZARDQ0OYmJjA3NwcKioqAn4UTVdRUQFFRUVcvnwZ48ePZzoOJQSio6NhYWEBFxcX7N+/n+k4lBhgtMgDwIwZM3Dv3j3ExcWJ/HClTCkpKcGdO3cQERGBiIgIPHnyBNXV1dDU1ET//v1haGjILYK9evUSm0PDr169QmxsLPcrJiYGKSkpkJSURN++fWFhYYEhQ4ZgxIgR0NbWZjruF3Xu3BkrVqzAr7/+ynQUimG5ubkYNGgQevXqhaCgILH5P6WYxXiRp4ftmyc+Ph7//fcfbt68ibt374LNZkNPTw9DhgzhFre+ffuK5TXm31JaWopHjx4hPDwcERERCA8PR1VVFfT09DBq1CiMGzcOY8aMgaysLNNRAQCjR4+Gjo4Ofe+3cpWVlbC0tERhYSEePnyI9u3bMx2JEhOMF3ng/5+DunTpEiZMmMB0HKFECEFERAT8/PwQEBCAnJwcaGpqwsbGBjY2Nhg9ejT9YPiCqqoq3Lt3D0FBQbh+/TpevHgBBQUF/Pjjj3BwcICtrS2jw+UuW7YM0dHRCA8PZywDxSxCCBwdHREcHIz79++jV69eTEeixIhQFHkAmDdvHi5duoSoqCjo6OgwHUdoxMXF4eTJkzh79izevHkDfX19TJ8+HePGjUO/fv1a3Z56S71+/RpBQUE4d+4c7t69C0VFRdjb22PmzJmwsrISeA/+Q4cOwd3dHYWFhQJtlxIea9euxR9//IGgoCBYWVkxHYcSM0JT5CsqKmBiYgJFRUXcu3evVV8XymazcfnyZXh5eSE0NBSdO3fGxIkTMXXqVFhYWDAdT2wUFBTgwoUL8PHxwf3796GnpwcXFxcsWLAAHTp0EEiGkJAQjBkzBvn5+VBVVRVIm5Tw+Pfff7FgwQIcP34cs2fPZjoOJY6IEHn+/DmRk5Mj69evZzoKI4qKioi7uzvR1NQkUlJSxN7enoSEhJD6+nqmo4m9hIQE4urqSpSUlIicnBxZsGABSU1N5Xu7aWlpBAB5+PAh39uihMvt27eJjIwM2bhxI9NRKDEmVEWeEEI8PT2JpKQkCQ4OZjqKwJSWlpKtW7eSdu3akfbt25MNGzaQt2/fMh2rVSotLSWHDx8mPXr0INLS0mT+/Pnk1atXfGuvpqaGSElJkTNnzvCtDUr4xMbGEhUVFTJjxgy6EU/xldAVeUIImT59OlFXVyfp6elMR+Gr2tpacujQIdKhQweirKxM3N3dSXFxMdOxKPKh+B4/fpzo6ekRGRkZsmLFClJSUsKXtnR1dcmOHTv4sm5K+KSnpxMtLS0yYsQIUlVVxXQcSswJzTn5j5WXl8Pc3BwsFgvh4eGQl5dnOhLPRUVFYdGiRYiOjsbPP/+MNWvW0HEChFBNTQ2OHz+O9evXQ1ZWFgcPHsSkSZN42saoUaOgq6uLo0eP8nS9lPApLCzE0KFDwWKxcPfuXfo/T/GdUA4GrqCggMDAQGRkZGD27NmfzTMuyurq6rBx40aYmJhAVlYWUVFR2L17N/1nF1LS0tJYuHAhEhMTMXr0aEyZMgWTJk1qMDFPS+np6eHly5c8Wx8lnCorK2FnZ4eysjJcu3aN/s9TAiGURR4AunbtCj8/P1y+fBkeHh5Mx+GJd+/eYezYsdi7dy8OHjyIu3fvQl9fn+lYVCN06NABx48fx61bt/Do0SMMGjQIcXFxPFm3rq4uXr16xZN1UcKprq4Ojo6OSExMREhICDp37sx0JKqVENoiDwBWVlb4/fffsX79evz3339Mx2mR5ORkDBw4EC9evMC9e/ewZMkSeo27CBoxYgSePn2KTp06YfDgwQgMDGzxOrW0tJCVlSVWR6yo/48QgkWLFiE4OBhXrlxB7969mY5EtSJCXeQB4KeffsKsWbPg6OiI2NhYpuM0y8uXL2FlZQUNDQ08efIEAwcOZDpSq1JWVsbT9WloaCA0NBSOjo6YOnUqrl271qL1aWpqorq6mjsbICVefv31V/j4+ODcuXMYPHgw03GoVkboizwAeHl5wdTUFLa2tkhPT2c6TpNkZmZyC3xQUBDU1NSYjtQihBDs378fu3fvRo8ePeDs7IyAgAB07twZL168YDpeA56enhg+fDj69OnD83VLSUnhyJEjcHZ2xuTJkxEWFtbsdWloaAD4MEEJJV62bduG/fv3w8fHB2PHjmU6DtUKiUSRl5aWxvnz56GmpoYff/wRRUVFTEdqFEIIZs+ejbZt2yI4OFgkpj/9nq1btyIpKQlr167F8ePHUVJSAhkZGairqzM6BvyXLFiwAPX19airq+PL+iUkJODl5YXx48dj5syZzR6aVlNTEwCQk5PDy3gUww4fPozNmzfj8OHDcHBwYDoO1UqJRJEHACUlJVy7dg3l5eWwt7dHdXU105G+6/Dhw7hz5w68vb0FNkwqvx0+fBhdu3YFAFhYWCAwMBDjx4/H06dPoaury2y4T7BYLL5PMSspKYljx45BRkYGrq6uzVqHqqoqWCwW3ZMXI6dOncLy5cuxa9cuLFq0iOk4VCsmMkUeADp16oQrV64gOjoa8+bNE+qOSu/evYObmxvWrFmDQYMGMR2HJ6qqqpCXl0c7DH5CSUkJx44dw9mzZxEaGtrk5VksFtq1a4d3797xIR0laBcuXMDcuXPx22+/Yc2aNUzHoVo5kSryAGBoaIiAgACcP38eP/30E9NxvsrLywvS0tJYu3atQNq7du0ali5dipUrV2Lw4MGfDawSEBCAZcuWYfXq1bC1tcWGDRu4R0Oio6Px66+/Qk9PDxUVFViwYAFUVVVhYmLCvX7b29sbLi4uAAB/f3+4uLjAw8MDRUVF+OeffzB69GhcunQJwId+CLt378YPP/yAd+/ewdraGjo6Orh79y7WrVuHXr16ITMzE9u2bYOOjg709fURFhaGqqoqrFq1Ct26dUOXLl0QHBzc4DEQQnDkyBEsWbIEpqamGDNmDFJSUhr8zeXLl7Fw4UKsWbMGy5cvR3Z2Nl+e70+NHj0ao0aNwr59+5q1fJs2bVBZWcnjVJSgXbx4EQ4ODli8eDF27tzJdByKEq4JapriwoULREpKSmgnd+jXrx9ZunSpQNry8fEhDg4OpK6ujhBCyI4dOwgAEhoaSgghZP/+/cTc3Jyw2WxCCCEFBQWkR48eZPjw4aS+vp5kZ2eTUaNGEQDE1dWVxMfHk6ioKCIrK0scHBy47RQUFBAAZPv27dz7EhISyKpVqwgAcv78eUIIIUFBQaR3796ExWKRzZs3Ey8vL2JiYkKio6OJs7MzAUAWLlxInj59SkpLS4mpqSnR09Mjrq6uJCEhgZSVlRFzc3Oip6fX4HHu2rWLnDhxghDyYUjgvn37Ek1NTVJRUUEIIeT06dPE1NSUVFZWEkIIyc/PJ6qqqkRTU5MfT/tn/Pz8iJSUFCkoKGjysj169GjwvFKi57///iOysrLExcWFjkdPCQ2RLfKEEOLt7U0kJSWJh4cH01EaKC0tJZKSkuTixYt8bysvL48oKyuTly9fcu/Lz88nkyZNIgkJCSQ3N5fIy8sTHx+fBssdP36cACAnT54khBDi5uZGADQoUBYWFqRHjx7c218q8oR8mE3r4yJPCCHz588nAEhKSkqDv/3rr78IABIbG8u9b/PmzQQAiYqK4t63ceNGAoDk5eURQgjJzMwkGhoa3A0ZQgjZtGkTAUD8/PxIRUUF6dixI/H19W3Qnr29vcCKfGFhIZGQkCBXr15t8rKGhoatdvZFcUALPCWspJg4esArs2bNQmlpKVasWAFlZWWh6eDy+vVr1NfXo1evXnxvKzw8HPX19Q06vamqqiIgIAAAEBgYiIqKCnTp0qXBcuPGjQMAhIWFwcnJCSwWC8CHS8M4tLW1kZqa+t0MHy/DIS0tDSkpKXTv3r3B/Zx2JCX//5kiTuc4aWlp7n2cvAUFBVBTU8P9+/dRU1Pz2Wu8YMECtG3bFvfu3UN2djYMDAwa/F5WVva7+Xmlffv2UFdXR1paWpOXbdu2LT1cL6KuXbuGyZMnY9asWfD09KR9ViihItJFHgCWLVuGwsJCLF26FEpKSpgxYwbTkbgf1m3btuV7W8+fP0dNTQ0IIV/8cOGMK/Bppy5VVVXIyckhKyuL7xm/50u5OffV19cDAF68eAF5efmvTuJy4MABAICMjAyfUjaOnJxcs4q1rKws2Gw2HxJR/HTt2jVMmjQJzs7OOHLkCC3wlNARuY53X7J582asXLkSs2bNgr+/P9NxuBNPCKK3tJKSEqqqqpCQkPDZ76qrq7l7+F+bAEVUhtiUk5NDRkYGMjIyPvtdfn4+t7gzOVgSIQSFhYVo3759k5etrq4W6FEHquU+LvCenp4Njk5RlLAQm3flvn374OrqihkzZuDUqVOMZtHT04OcnByePXvG97Y4l+dt2LCBu9cLAE+fPsXVq1cxePBgKCkpcXu+c2RkZOD9+/eYMGFCo9siDF6yaGBgAELIZ5ckpaWl4fDhwzA0NAQAnD17tsHv+TkYzqdSUlJQWlqKH374ocnLVlVVCd1gQtTXBQUF0QJPiQSRP1zPISEhgf/9739o27Yt5syZAwBwcnJiJAuLxcLIkSNx6dIlLFiwgK9tmZubw9bWFpcuXYKVlRWmTJmC9PR0vHv3DseOHQMAeHh4YOnSpQgNDYWVlRUA4ODBg5g9ezYsLS0BgDtuem1tLXfdeXl5eP/+Pfc2Zy/64/sAcC9Ty8/P595XXl6Ouro6FBcXNxjpr7S09LN2OPcVFBRw7+OMN8+5zG/06NEYNGgQfH19UVVVBXt7e5SWluLChQvw8/ODqqoqLC0tceLECRgbG2P27NmIj49HeHg48vPzcebMGdjZ2UFOTq6Jz3DjXbp0CR06dICxsXGTl6V78qIjKCgI9vb2cHJyogWeEn7M9vvjDzc3N8JisT7rUS5I58+fJ5KSkiQxMZHvbVVUVJAlS5YQLS0toqGhQZYsWUKKi4sb/M2lS5fImDFjyLJly8jGjRvJvn37uL2Ab968Sbp27UoAkKVLl5K8vDzi4+NDFBQUCADi7u5OHj16RGbMmEEAEF1dXXL69GlSXFxMQkNDybBhwwgAMnDgQHLjxg3i5eVF1NTUCADi7OxMnj17RgghJDQ0lBgaGhIAxNHRkaSmppLbt2+T/v37EwDExsaGxMbGkvDwcDJgwAACgDg5OZG0tDRCyIfe646OjkRdXZ2oqamRWbNmkczMTO5jLCkpIXPnziUaGhqkS5cuxN3dnSxcuJDMnTuX3Lx5s0HPfF5js9lER0eHrFixolnLd+3alezZs4fHqSheu3btGpGVlSXz58/n6/uJonhFLIs8IYSsW7eOsFgs4u3tzUj7dXV1xMDAgIwcOZJeUtMKbN68mcjLy5P09PRmLa+urk4OHTrE41QUL50/f57IyMgQFxcXWuApkSG2x5l27NiBtWvXYt68efD29hZ4+5KSkvDx8cG9e/dw+PBhgbdPCU5UVBR27twJDw+Pzy5VbIz6+nq8e/cO6urqfEhH8cK///6L6dOnw8XFBUeOHKGH6CmRITbn5L9k+/btkJCQ4I5zzzlXLyj9+vXDunXrsHr1avTu3Zt7PpwSH2/fvsWkSZMwbNgwLF26tFnrKCgoQG1tLS3yQurQoUP46aef8Ouvv2L37t1Mx6GoJhHrIg98mM+ZU+irqqqwePFigba/adMmpKWlYcKECbh27RqGDx8u0PYp/snLy8OYMWOgqKiIs2fPNvsaac7sc7TICx8PDw+4ublhz549WL16NdNxKKrJxL7IAx/mQJeQkMDSpUtBCMGSJUsE1rakpCSOHz+OyspKjB8/HidPnoSdnZ3A2qf4IyEhAfb29pCUlERISEiLphLOy8sDAGhoaPAqHtVChBD8/PPPOHToELy8vPh+lQxF8UurObG0ZcsW7Nq1C66urgI/Ry4lJQVfX1/MmDED9vb22LBhg8Cu3aZ4z9/fH6amplBTU8OtW7daXJzz8vIgJSXFHUSJYlZdXR3mzZuHw4cPw8/PjxZ4SqS1ij15jjVr1qCqqgrLli0Di8US6Fj3MjIy8PT0hImJCZYtW4aIiAh4eXmhR48eAstAtUxxcTHWrFmDo0ePwtXVFX/88UeD8fabKz09HVpaWrQzlxCorq7GjBkzcOPGDVy5cgVjxoxhOhJFtUir+1TZvHkztm7diiVLlmD//v0Cb3/+/PkIDw9HUVERDA0NsXXrVu6AL5TwOnv2LPr06YPLly/Dz88Phw4d4kmBB4BXr15BT0+PJ+uimq+8vBzjxo3D7du3cePGDVrgKbHQ6oo88GEI2EOHDuGXX37B2rVrBd6+sbExnjx5gm3btmHPnj0wMDCAr69vg2FpKeEQEREBKysrzJgxA+PHj8eLFy8wbdo0nrbx8uVLWuQZ9u7dO4waNQrPnz/H7du3YW5uznQkiuKJVlnkAcDV1RXe3t7Yt28fli9fLvBx2aWkpLB69WrEx8fDzMwMs2bNgoGBAc6dO0eLvRB4+PAhbGxsYGFhgdraWoSHh8PLy4sv581pkWdWeno6hgwZgvz8fERERHDnQaAocdBqizwAODs74+TJk/D09MTixYsZKa46Ojrw8fFBSkoKLCwsMHPmTPTs2RMeHh4oLCwUeJ7WrK6uDleuXMHo0aMxePBglJaWIjAwEHfu3OHbnl1tbS3evn1LizxDnjx5AjMzM8jKyiI8PJy+DpTYadVFHgAcHBxw8eJFnDx5Eo6OjqipqWEkh66uLjw9PREXFwdra2vs2LEDXbp0gYuLCx4+fMhIptbi1atX2Lx5Mzp37gx7e3soKCggJCQE9+/fx/jx4/na9ps3b1BTU0OLCwOCg4MxcuRI/PDDD7hz5w46duzIdCSK4j1mR9UVHrdv3yaKiopk3LhxpLKykuk4pKSkhPz555+kT58+BADR09Mj69atI7GxsUxHEwuZmZlk//79xNTUlEhISBB1dXXi5ubW7LHnmysoKIgAIIWFhQJtt7Xz8vIiUlJSZN68eYTNZjMdh6L4RoIQBicJFzKPHz+Gra0tDAwMEBgYCEVFRaYjAfgwN7yfnx/Onj2Lt2/fom/fvhg7dixsbW0xZMgQyMjIMB1R6BFCEBMTg+vXr+PatWuIiIiAoqIi7Ozs4ODggFGjRvGst3xT7N+/Hx4eHsjJyRF4260RIQRbtmzB1q1bsWnTJri7uzMdiaL4ihb5T0RFRcHa2hpdu3bF9evX0b59e6YjcRFCEBERgYCAAAQFBSEpKQmKioqwsrLCmDFjMHToUPTt25deb/1/3r59i3v37iE0NBTXr19HVlYWNDQ0YG1tjYkTJ8LW1hZt2rRhNOPChQuRmpqKW7duMZqjNaiursbs2bNx6dIl/PPPP3B0dGQ6EkXxHS3yX5CYmIjRo0dDTU0NwcHBUFNTYzrSF7169QrXr19HUFAQbt++jbKyMqioqMDc3ByDBw/GkCFD0L9/f6ioqDAdle+qqqoQFxeHyMhI3L9/HxEREcjIyIC0tDRMTU1hY2MDGxsbDBgwoNljzPPD0KFDYWhoiL/++ovpKGKtsLAQEydORHx8PC5evEjnkKBaDVrkvyI9PR2jRo2ClJQUQkJCoK2tzXSkb6qrq0NsbCwiIiJw//59hIeH4+3btwA+9OA3MjKCoaEhjIyMoK+vDz09PcjKyjKcuunq6urw5s0bJCUlISYmBtHR0YiNjUVycjJqa2vRrl07mJubw9zcHEOGDMGgQYMgJyfHdOyvUlVVhbu7O5YtW8Z0FLGVmpqKsWPHoqamBteuXUPv3r2ZjkRRAkOL/Dfk5ORgzJgxKCsrw82bN9GtWzemIzVJVlYWtwh+XAzr6uogKSkJbW1tdOvWDd27d0e3bt2gra0NLS0taGhoQEtLC0pKSgLPXFlZiezsbGRnZ1+lKgkAACAASURBVCMnJwcZGRlIS0tDamoq0tLS8Pr1a7DZbACAtrY2jIyMYGRkhH79+sHIyAg9evQQqj31b8nLy4OGhgZu3rxJpyHmk8jISNjZ2UFPTw+BgYF0pj+q1aFF/jvy8/NhY2OD3NxchISEoE+fPkxHapHKykokJycjLS2N+8UpoFlZWdwCCgBt27ZFp06d0L59eygrK0NFRQXKyspQUlKCkpIS2rZtC1lZ2QZ7yioqKtwiW1ZWhtraWgAAm81GRUUFampqUFpaipKSEhQXF3N/fvfuHXJzc1FcXMxdl4SEBDQ0NKCnp9dgY6Rbt27o2bOnUPWXaI47d+5gxIgRyMrKopdv8cGpU6ewYMECjB07FqdOnULbtm2ZjkRRAkeLfCMUFxdj3LhxSEpKwtWrV2FiYsJ0JL7Jy8tDbm4uMjMzud+LiopQWlqKoqIilJSUcAtzdXU13r9/zx17v76+HiUlJdx1ycvLc3v+s1gsKCkpgcViNdhgUFZWRkxMDHJycrBx40Z06tQJHTt2RKdOnaCurg4pKfGdQ+nPP//Epk2b8O7dO6ajiJW6ujq4ublh7969+OWXX+Dh4UE7o1KtFi3yjfT+/XtMmzYNt2/fhr+/P2xtbZmOJDaePXuGgQMH4tKlS5gwYQLTcQRmzpw5yMrKwo0bN5iOIjbKysrg5OSE4OBgHDlyBHPmzGE6EkUxim7eNpKcnBwuX74MBwcHTJw4EWfOnGE6ktgYMGAA7O3tsWHDhlY1bv+TJ08wcOBApmOIjZSUFJiamuLJkye4c+cOLfAUBVrkm4TFYuHo0aNYvnw5nJyccODAAaYjiY2tW7ciPj4e58+fZzqKQJSXlyMxMRHGxsZMRxEL169fh4mJCdq2bYvIyEiYmpoyHYmihAIt8k0kISGBvXv3Yvfu3Vi1ahVWrVqFuro6pmOJPH19fUyfPh0bN27kdtYTZ1FRUairq8OgQYOYjiLyDhw4gHHjxsHW1hbh4eHo0qUL05EoSmjQIt9Mv/76K86dOwdPT0+MGzcOZWVlTEcSeVu2bMHLly9bxamQx48fQ01NjRakFqiursacOXPwyy+/YMeOHfD19aU96CnqE7TjXQs9ePAAdnZ20NTUxNWrV4V+0BxhN3fuXNy7dw8vXrxgZCx5QXF0dERxcTGuXr3KdBSRlJWVBXt7eyQmJuLUqVN8ny2QokQV3ZNvITMzM0RERKCqqgpDhgxBbGws05FEmru7O96+fQtvb2+mo/AV7XTXfJGRkRg4cCDKysrw5MkTWuAp6htokeeB7t27IzIyErq6uhgyZAguXLjAdCSRpaOjg7lz52Lbtm3c6+/FTUlJCVJTU2mnu2b4888/MWLECBgbG+PBgwfo0aMH05EoSqjRIs8j7du3x82bN+Hq6oopU6Zg7dq1repyMF7atGkT8vPzcezYMaaj8MWjR49QX19PO901QXl5OWbOnIkVK1Zg1apVuHz5MiPDLlOUqKHn5PnAy8sLy5cvx5gxY3D69Gn6YdQMK1euhL+/P1JTU4V6gpnmcHNzw8WLF5GYmMh0FJGQlJSEKVOmICcnB6dOnYK1tTXTkShKZNA9eT5YuHAhbty4gYcPH2Lo0KFITU1lOpLIcXNzQ0lJCY4cOcJ0FJ67desWRo4cyXQMkXDx4kWYmpqiTZs2ePz4MS3wFNVEtMjzyfDhw/Ho0SNISUlh0KBBCAwMZDqSSNHU1MSyZcuwa9cusbo8saSkBE+fPqVF/jtqa2uxdu1aTJo0CdOnT0dERAS6du3KdCyKEjm0yPNR165dcf/+fcyaNQsTJ07EypUrW8VAL7zy22+/gc1m488//2Q6Cs/cvn0bhBCMGDGC6ShCKyMjA8OHD8dff/2FM2fOwNPTkzvREUVRTUOLPJ/JysriwIEDOHr0KLy8vGBra4v8/HymY4mEDh06YOXKldizZw+KioqYjsMTt27dgpGREVRVVZmOIpTCwsIwcOBAFBYWIjIyEg4ODkxHoiiRRou8gMyfPx+RkZF4+fIlDAwMcP36daYjiYRffvkFkpKSYjNPAD0f/2X19fXYunUrRo0aBSsrKzx58gQ//PAD07EoSuTRIi9A/fr1Q1RUFEaOHIkff/wRK1euBJvNZjqWUFNWVsaqVauwb98+kT8CkpeXh/j4eFrkP5GVlYXRo0dj586dOHjwIE6fPg0FBQWmY1GUWKBFXsCUlJTg6+uLEydO4J9//oGFhQXtff8dP/30E+Tk5PDHH38wHaVFQkNDwWKxYGFhwXQUoRESEgJjY2O8fPkSt2/fhqurK9ORKEqs0CLPkFmzZuHJkyeoqamBsbExTp8+zXQkoaWgoIDVq1fj0KFDyM3NZTpOs4WGhsLExISOmwCgpqYG7u7usLGxwdChQxEdHQ0zMzOmY1GU2KFFnkG9e/fGgwcPMGfOHDg7O2PWrFkoLy9nOpZQWr58OVRUVODh4cF0lGapr6/H1atXMXbsWKajMO7169cYNmwY/vjjD5w4cQLnzp2DsrIy07EoSizRIs8wTu/769evcw9dPnv2jOlYQqdNmzb47bff8PfffyMjI4PpOE328OFD5OTkwM7OjukojPL390e/fv3AZrPx9OlTODs7Mx2JosQaLfJCYsyYMYiOjoaenh7MzMzg7u5Ox77/xOLFi6GpqYldu3YxHaXJAgMDoaenB319faajMKKsrAyzZs3C9OnTMXv2bNy/f59OLkNRAkCLvBDR0NDAtWvX8Pvvv2PXrl0YPnw47ZT3ERkZGbi5ueHo0aN49eoV03GaJDAwsNXuxT979gzGxsYICgrClStXcODAAcjKyjIdi6JaBVrkhYyEhARWrlyJ8PBwFBUVYcCAAfD09ASdR+iDuXPnokuXLtixYwfTURotLS0NCQkJmDBhAtNRBKq+vh579+7F4MGD0aVLF8TGxtI+CRQlYLTIC6lBgwYhKioK69evx/Lly+le/f+RlpbGxo0bceLECSQlJTEdp1EuX76M9u3bt6pL516/fo2RI0di/fr12LJlC27cuIGOHTsyHYuiWh1a5IWYtLQ01qxZg8ePH6OsrAxGRkbw8PBo9efqnZyc0LNnT2zfvp3pKI0SGBiIH3/8EVJSUkxHEQh/f38MGDAA+fn5iIyMxNq1ayEpST9qKIoJ9D9PBBgZGeHRo0f49ddfsX79etjY2ODNmzdMx2IMi8XCpk2b4Ovri9jYWKbjfNO7d+8QERHRKg7VFxcXw9HREdOnT8fUqVPx+PFjDBgwgOlYFNWq0SIvIqSlpeHu7o6IiAhkZGTAwMAAXl5erfZc/bRp02BgYIBt27YxHeWbrl69CklJSbGfBz04OBj6+voIDw9HaGgoPD09IScnx3Qsimr1aJEXMaampnj27BkWLVqEpUuXwsbGBq9fv2Y6lsBJSkpi8+bNCAgIQFRUFNNxvsrPzw/W1tZiO8pdZWUlVq5cCVtbWwwZMgRRUVGwtLRkOhZFUf9HgrTWXUEx8ODBA8ybNw/p6enYvHkzVq1aBWlpaaZjCZSpqSk0NDQQGBjIdJTPFBUVQVNTE//++y8cHR2ZjsNzDx8+xKxZs5CXl4c//1979x5VY7rHAfy7L23STXShISPRxSW17eHIpXRGOLFiGDI0WihzHDKYU4sxdTRGB+PouJ2YKVouo6aDSQpNDSVJRREt7eYohe4iu3a7vZ/zh2WPxt3s3dvO77NWa7H32/N83xb99vO+z/O8O3Z0ynMkRNfRSF6HjRo1Cvn5+fjmm2+wfv16iMViZGZmch2rXYWEhCAhIQFZWVlcR3lGXFwc+Hw+pk6dynUUjZLL5VizZg1cXV0xYMAAXL9+nQo8IR0UjeQ7iV9//RVLly7FqVOnMG/ePGzduhVmZmZcx2oX48aNg76+Pk6dOsV1lDY8PDzQs2dPxMbGch1FY7Kzs+Hn54fbt29j06ZNCAgIAI/H4zoWIeQFaCTfSdjY2CApKQnHjx9HWloa7Ozs3pmJeSEhITh9+jTOnj3LdRS1e/fu4ezZs5g9ezbXUTSiubkZwcHBGD16NHr27Im8vDwsWbKECjwhHRyN5DuhhoYGfPXVV9i5cyfGjBmD3bt3w8HBgetYWjVhwgQoFAqkp6e3e98PHjx4ZmLd9u3bsXbtWlRWVkJfX7/dM2nSk7kfFRUV2Lx5MxYvXkzFnRAdQSP5TsjExAQRERHIyMjA/fv34eLigvXr10Mul3MdTWvCwsKQkZGBlJSUNq9nZ2cjPDxcq317e3vD3d0dhw8fRnNzMwDgyJEj8Pb21ukC39TUhODgYIwZMwbW1ta4evUq/P39qcAToksY6dQUCgXbtm0bMzIyYn379mX79+/nOpLWeHp6MolEwlQqFbt8+TKbMmUKA8D69++v1X7FYjHj8XiMx+MxQ0ND5uvry3g8Hjtx4oRW+9Wm8+fPM3t7e2ZiYsIiIyOZSqXiOhIh5C1QkX9HVFRUsPnz5zMej8fc3d1ZQUEB15E0Ljs7mwmFQubh4cF4PB7T09NjAJienh5rbW3VWr8ODg4MgPrrSb+2trYsPDycVVVVaa1vTWtsbGSBgYGMz+ezKVOmsPLycq4jEUL+ALpc/46wsrJCTEwMLl68iKamJri4uCAgIAA1NTVcR9OIW7duYc+ePVCpVDh37hwYY1AoFAAAhUKB8vJyrfX9+9sgT/otKSnB2rVrYWVlhVmzZnX45Y3JyckYMmQIYmJiEBUVhcTERLz33ntcxyKE/AFU5N8xEokEmZmZ+P7773H8+HHY2dkhIiICSqWS62hvpaysDP7+/rC1tcX+/fuhUqnURfZp2nyC35P78L/HGINSqYRSqcSJEyc67DPU6+rqEBAQgMmTJ2PIkCEoKCjAp59+ynUsQogGUJF/B/F4PPj6+kIqlWLZsmUICgrC0KFDO9w681epq6vDqFGjsHfvXiiVyucWdwAQCoUoKSnRWo7XmdB46NAhiMVirWV4W3FxcbCzs8OJEycQHx+PhIQE9OnTh+tYhBANoSL/DjM0NERoaCiuXLmCfv36YdKkSfj44491Zi/8Hj164MyZM7CwsHjpdr4CgUCrI/mWlpYXvsfn8/HPf/4T06dP11r/b6OkpAQTJ07EnDlzMGPGDBQVFWHGjBlcxyKEaBgVeQJ7e3skJSUhISEB+fn5sLe3x6pVq1BXV8d1tFcaPHgwLl68iN69e7+w0Le0tKC4uFhrGV40khcKhZg/fz6++OILrfX9plpbWxEREQEnJydUVlYiMzMTkZGRMDIy4joaIUQLaDMc0oZCoUB0dDRCQkIgl8sRFBSE5cuXd/j13nfv3oWbmxv+97//Pfeyvb29PW7cuKGVvvl8/jM7C+rp6UEikSAtLQ0ikUgr/b6prKws+Pv7QyqVIiQkBKtWrYJQKOQ6FiFEi2gkT9rQ09NTF4KgoCBs2LABgwYNwp49ezr05LzevXsjPT0dgwYNem7hKi0t1coWvy0tLc+0KxQK0adPH5w4cULrBf7Ro0evPKa2thb+/v5wdXWFpaUlrl69iqCgICrwhLwDqMiT5zIwMEBQUBCKioowZcoULF26FE5OTjhx4gTX0V7IwsICGRkZcHZ2fubSfVNTEyorKzXe5+9n1gsEAnTr1g2nTp2Cqampxvt72rlz5+Dq6vrCCYeMMcTExMDR0REJCQmIjo7G6dOnMWDAAK3mIoR0HFTkyUtZWVkhMjISV69ehaOjI6ZOnYoPP/wQly9f5jrac3Xv3h1paWlwdXV9ZqSqjcl3z7sff+zYMQwcOFDjfT2tuLgY06ZNQ35+Pnbu3PnM+wUFBRg7diz8/Pzg7e2NoqIi+Pr60pa0hLxjqMiT12Jvb4/Y2FhkZmaiqakJI0aMwMcff4ybN2++8nsLCgqQmpraDikfMzAwQHJyMiZPngyBQADg8QhbG8vonh7J83g8REVFwd3dXeP9PK22thYTJ06ETCYDAKxbtw7V1dUAHl++Dw0NhUQiQXNzM7KyshAZGQkTExOtZiKEdEx0U468kT/96U9IT09HXFwcQkJCMHjwYCxYsADr1q2DtbX1c7/nH//4BxITE5GcnAw3N7d2ydmlSxf8+OOP8PHxwfHjx8EYe+5IvqGhAZWVlXj48CHu378Pxhjkcrm6gBoaGkJPTw8CgQDGxsbo3r07LCws1LPRny7ya9asga+vr1bPq6WlBdOnT0dFRYX6Mr1cLseaNWswbdo0LF26FI2Njdi0aROWLVsGPp8+xxPyTuNqP12i+5RKJYuNjWW2trZMJBIxf39/VlFR0eaYwsJC9cNbunbtys6ePduuGR88eKB+UM3gwYOZj48Pk0gkzMrKiolEojZ7zr/JV5cuXdh7773Hhg0bxgAwBwcHFhMTwy5dusSam5u1ci4qlYp98sknTCAQPJOHx+MxPp/PFi5cyGpqarTSPyFE99ASOvKHPVl2t379etTX12PRokVYs2YNLC0t4ePjg/j4eCgUCvD5fOjp6Wl1RF9YWIi0tDTk5OQgNzcXN27cgFKphEAggEgkwocffogBAwbA2toaZmZmMDc3h6WlJYyNjdUT5UQiEQwMDAAAjY2NUCgUYIzh/v37qK+vR1VVFWpqalBdXY3s7GwkJyfDwsICpaWlaGlpgZ6eHoYMGQKxWAyJRAJ3d3eN3KMPDQ3F+vXrn7tKQCgUwtHREVeuXKH77oQQNSryRGOampqwe/duhIeHQy6XY8GCBdixYwdUKpX6GD6fD5FIhOTkZIwfP/4P91lbW4uffvoJZ86cQWpqKiorK9G9e3dIJBKIxWL11/vvv48tW7ZofGMaqVQKQ0ND9OrVCyqVClKpFLm5ucjNzUVeXh4uXbqExsZG9O3bF+7u7vD09ISXlxeMjY3fqJ8jR47Ax8fnpcsAeTweDh8+jNmzZ//R0yKEdBJU5InGPXr0CDt27EBYWBhaWlqeWeL1pNCfOnUK48aNe+P26+rqEB8fjx9//BGpqakQCoUYN24cJkyYgAkTJsDFxUU94Y5rra2tuHjxIlJTU5Gamorz589DIBBg4sSJmDlzJry9vV+521xGRgYmTJiA1tbWlxZ5Pp8Pc3NzlJSUqK9EEELebVTkiVaUlZVhwIABaG1tfe77AoEAXbp0QWpqKkaOHPlabebm5mLPnj04cOAAVCoV/vznP2PWrFnw9vZ+45ExV+rr65GQkIC4uDicOXMGIpEIPj4+CAgIgIuLyzPHFxcXQyKRoLGx8bU3IwoJCUFoaKiGkxNCdBFNvSVasWHDhpfeG1YqlZDL5fDw8EB2dvYLj1OpVIiNjcXw4cMxYsQI5OTk4F//+heqq6uRkJAAX19fnSnwAGBqagpfX18kJCTg7t27CAsLQ3p6OsRiMUaPHo3ExET1sU8vlXtegefxeNDT01P/nLt164bRo0er5xAQQgiN5InGlZeXw8bG5oU7sT1NIBCga9euSEtLg0QiUb/OGMMPP/yAr7/+GkVFRZg9ezY+//zzNsd0JmfPnsWWLVuQmJgIFxcXrFmzBlu3bsX58+cBPL4ULxAI1D9TMzMzjBgxAmKxGM7OznB2doaNjQ2Xp0AI6YCoyBONCwwMxL///W/o6em9dqE3MDDAL7/8AmdnZ+Tm5mLZsmXIzs6Gj48P1q5dC3t7+3ZIzr28vDyEhYXh2LFjAB6P1vv27YsPPvgAYrEYw4cPh7OzMywtLTlOSgjRBVTkicbFx8fj5s2buH37Nm7dugWpVIqKigr1BjMA1BvMKBQK9aVoY2NjeHp6Ij4+Hq6urti+fTucnJy4Og3O5OfnY9++fUhMTMStW7ewcuVKrF+/vsM8zY4QojuoyJN2c//+fZSVlaGsrAylpaUoKyvD7du3IZVKIZVKUV9fDz6fj6+//hrBwcHv/HpvpVKJyMhIBAcHY+DAgTh8+DAGDRrEdSxCiA6hIk84FxkZicDAQIwePRobN26EoaEhBg8ezHWsDqO4uBhz585FUVERoqKiMGvWLK4jEUJ0BM2uJ5xhjOHLL7/EZ599huDgYKSkpGDkyJFU4H9n4MCBOH/+PBYsWIA5c+Zg27ZtXEcihOgIekAN4QRjDEuWLEF0dDSioqKwYMECriN1aCKRCNu3b0e/fv2wcuVK1NbWIiwsjOtYhJAOjoo84cRXX32F6Oho/Pe//4WXlxfXcXTG6tWrYW5uDj8/P1haWuJvf/sb15EIIR0YFXnS7qKiorBhwwZ8//33VODfwqeffop79+4hMDAQ1tbWmDZtGteRCCEdFE28I+1KKpVi+PDhWLZsGTZu3KjVvh4+fPjKfeH/qMbGRhgaGmq1jxdZuHAhfvrpJ1y7do3WzRNCnouKPGk3KpUK48aNg0wmQ1ZWltbWfUdGRuLQoUMoKSlBeXm5Vvo4ePAgoqOjUVhYiLt372qlj1dpbGzEsGHD4OTkhKNHj3KSgRDSsdHsetJujh49igsXLiA6OlqrG7ssWrQIKpXqtR/o8rqeLuZz5syBUql84QN42oOhoSH27t2LY8eOISMjg7MchJCOi4o8aTfh4eGYMWOG1nexEwgE6NOnj0bbrK+vx7x587Tax9vw8PDA2LFjtX7rgxCim6jIk3Zx4cIF5OTk4O9//zvXUd6YTCbDnDlz8Ouvv3Id5bmCgoKQlJQEqVTKdRRCSAdDs+tJuzh58iT69++vtafIHT9+HImJiTA1NYVMJnvmPjljDJGRkcjPz0deXh5MTEywc+dODBw4EABQWVmJL7/8EtbW1igrK0NNTQ2+++479OzZE0ePHsWNGzdQX1+PxYsXw87ODqtXr1a3fe/ePSxZsgTnzp1D//79ceDAATg4OAAArly5goiICNjb2yMzMxMymQxnzpzR6Ll7enrCxMQESUlJWLZsmUbbJoToOEZIO5BIJOyzzz7TStsHDx5kI0eOZE1NTYwxxqqrq5mZmRnr1auX+piNGzeyffv2McYYa21tZY6OjqxXr17s0aNHjDHG3Nzc2OzZs9XHOzk5sXnz5qn/7uXlxd5///02/c6bN48ZGBiwFStWsKKiIlZQUMAMDAyYl5eX+phBgwaxjIwMxhhjMpmMjRkzRsNn/9jMmTPZX/7yF620TQjRXXS5nrSLmzdvwtnZWePtymQyrF69GoGBgejatSuAx89aHzt2rPqYO3fuYNu2bZg/fz6Ax/fTZ86ciXv37iEhIQHA40e6Pj1XYMiQISgoKHhl/0KhEJs3b4adnR2GDh2KDz74ALm5uQAAhUKB4uJi9d/19fWxatUqzZz47zg5OaG4uFgrbRNCdBddrida19zcjIaGBq2s5U5PT8fdu3cxdOjQNq936dJF/efMzEwoFAoEBAS0OWbRokXQ19cHAKSmpqqzHjx4ENnZ2WCvsbpUT08PQuFv/41sbGxw4cIF9Xuenp5YsWIFrl27hvDwcHh7e7/dib5Cr169cO/ePa20TQjRXVTkidY9evQIALSyaUxRUREAvHRJ3o0bN2BgYIC9e/e+8BilUolNmzYhJycHy5cvx8iRI5GVlfXGeX7/eNz4+HgsXrwYe/fuxdGjRxEbGwt3d/c3bvdVjI2N0djYCMbYO/+IXkLIb+hyPdE6U1NTCIVCVFdXa7ztJ8W9tLT0hcd069YN5eXlz90Yp7q6GiqVClOmTMH169cRHx+P8ePHayyfUCjEwYMHcfDgQQiFQkyaNAk3btzQWPtPVFZWwtzcnAo8IaQNKvJE6/h8PszNzXHnzh2Ntz1s2DAAwJEjR9q8/vRmOEOHDgVjDEFBQW2OKSkpwa5du5CdnY3Tp0/Dzc1N/Z5CoWhzuZ7P56OxsfGNssnlcuzZswcAMHfuXGRlZYExhrS0tDdq53XcuXOHtrYlhDyDijxpFyNGjEB6errG23V1dYW7uzv27duH3bt3QyaT4dKlS8jIyEB1dTUOHz4MV1dXSCQSHDp0CB999BEOHDiAXbt2ISAgAEuXLlWPfvfv34+rV68iKioKhYWFqKysREFBASorK2FlZYWamhrk5ubil19+gUwmQ21tLe7fv4+WlhZ1nqqqKsjlcshkMgCPH8bz5MOGlZUVTExM4OLiovGfw7lz57S2PJEQosM4ndtP3hk7d+5kRkZGTC6Xa7zthoYG5ufnxywtLZm1tTULDQ1l/v7+zM/Pj6WkpDClUslqa2vZJ598wiwsLJi5uTnz9fVlFRUV6jaWLFnCjIyM2KhRo1hKSgo7efIkMzMzYzNnzmSNjY0sPz+f9enThw0aNIjFxcWxmJgYZmpqygCwwMBA1tDQwKKioliPHj3Urz148IBJJBLm6enJwsPDmb+/P9u7d6/Gz7+uro4JBAIWFxen8bYJIbqNHlBD2sXt27dhY2ODqKgo9VI2ohnffvstQkNDUV5eDhMTE67jEEI6ECrypN34+voiJycH165dA59Pd4o0QS6Xw8bGBnPnzsXmzZu5jkMI6WDoNy1pN8HBwbh58yaio6O5jtJpfPvtt6ivr8fnn3/OdRRCSAdERZ60G0dHR6xcuRIrV6586ZI38nquX7+OsLAwhIWFwcrKius4hJAOiC7Xk3bV3NwMsVgMU1NTpKSkqLeiJW+moaEBY8aMgYmJCc6dO0e3Pwghz0W/GUi76tq1K+Li4nD9+nXMmzcPKpWK60g6Ry6XY/r06aivr8ehQ4eowBNCXoh+O5B25+joqH407KJFi9Da2sp1JJ3R1NSE2bNnIy8vD0lJSbC2tuY6EiGkA6MiTzgxduxYxMfHIzY2Ft7e3ur97cmL1dXVYeLEiUhPT8fJkyefeSgPIYT8Ht2TJ5y6ePEipk6dir59++KHH37AwIEDuY7UIV25cgVz5sxBc3MzkpKS4ODgwHUkQogOoJE84dTIkSORmZkJHo8HFxcX7N+/n+tIHQpjDNu2bcOoUaPQu3dvZGZmUoEnhLw2KvKEc7a2tsjMzERAQAD8/Pzg5eUFqVTKdSzOFRQUwM3NDV988QXWrVuHlJQUWipHCHkjw6PFIAAAAu9JREFUVORJhyASibBlyxakpqbi1q1bGDp0KNatW4eHDx9yHa3d1dbWYsWKFRCLxZDL5cjKysLatWshEAi4jkYI0TF0T550OK2trdixYwdCQ0MhEAgQGBiI5cuXo3v37lxH06qqqips3boVu3btgr6+Pr755hv4+fnREjlCyFujIk86rPr6ekRERCAiIgKMMSxcuBD+/v6ws7PjOppG5efn4z//+Q9iYmJgZGSE1atXY8mSJTA0NOQ6GiFEx1GRJx1eQ0MDdu3ahT179qC0tBRubm5YuHAhpk2bBiMjI67jvZW6ujocPXoU3333HbKysmBvb4+//vWvWLhwIbp168Z1PEJIJ0FFnugMlUqF5ORkREZGIikpCQKBABMnTsRHH32EyZMnw9zcnOuIL3X37l0kJCQgPj4eaWlpEAgE8Pb2RkBAAMaPHw8ej8d1REJIJ0NFnuikuro6HD9+HPHx8Thz5gxaW1sxbNgweHh4YMKECRg1ahR69OjBacaqqipcuHABP//8M37++Wdcv34d3bp1w6RJkzBz5kx4eXnp7JUIQohuoCJPdN6DBw9w9uxZdTEtLCwEYwz9+/eHi4sLXFxc4ODggAEDBsDW1lbjl8MfPHiAkpISlJSUoLCwEHl5ecjLy0N5eTn4fD6GDx8ODw8PeHh4YOzYsXQ5nhDSbqjIk06nuroaOTk5yMvLQ25uLi5fvozS0lI8+afeu3dv9OvXD2ZmZjAzM4OFhQV69OgBAwMDiEQiCAQCGBsbA3g8H0ClUkEul0Mmk6G6uho1NTWoqalBVVUVSktLUVVVBQAQCARtPliIxWL1E/cIIYQLVOTJO6G5uRklJSWQSqUoKSlBeXl5m2JdX18PmUwGuVyO1tZW9fp8ExMT8Pl8dO3aFfr6+uoPBmZmZjA3N0ffvn1ha2sLW1tb9O/fHyKRiOMzJYSQ31CRJ4QQQjop2mWDEEII6aSoyBNCCCGdFBV5QgghpJMSAojjOgQhhBBCNO//5ExnTmGAZxEAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from IPython.display import Image, display\n",
"model.view_model()\n",
"display(Image(filename=\"causal_model.png\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This causal model is arbitrary, and proposed mainly as an example.\n",
"\n",
"We asssume that the number of deaths is influenced by the number of already confirmed cases, the availability of hospital beds (the number of beds is taken as a proxy measure of the availability and quality of healthcare), by adoption of the policy, and other unobserved confounder; the number of confirmed cases is taken to depend on the adoption of the chosen policy; the introduction of the chosen policy is affected again by the quality of healthcare expressed by the number of beds pro capite and by other unobserved confounders.\n",
"\n",
"We take as **outcome** variable the number of deaths, as **treatment** variable the specific policy under study, as **confounder** the number of hospital beds (plus potentially unobserved confounders). The number of confirmed cases constitutes a **mediator**.\n",
"\n",
"This model will allow us to answer (within this model) the question: *What is the effect of the chosen policy on the number of deaths caused by covid19?*\n",
"\n",
"Several questions may be asked about the meaningfulness of this model: does the choice of number of beds make sense as a confounder? Is the number of deaths affected only the mediator of confirmed cases? Is this model complete at all? These are valid questions, but we do not provide an answer as we lack domain knowledge. Instead we will study this model to provide an example of such analysis.\n",
"\n",
"**Proper modelling, using knowledge from epidemiology and sociology, should be done here.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Analysis\n",
"\n",
"We now go through the standard steps of causal analysis (after modelling): (i) *identification*, (ii) *estimation*, and (iii) *refutation*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Identification\n",
"\n",
"We evaluate the identifiability of the causal effect, that is whether we can evalute the causal effect of the policy under consideration on the outcome given the model we have specified."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:dowhy.causal_identifier:Common causes of treatment and outcome:['U', 'hospital beds']\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:[]\n"
]
}
],
"source": [
"identified_estimand = model.identify_effect()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that *dowhy* warns us about the presence of unobserved confounders that may affect our analysis. We go on with this analysis even if we are obviously aware of the fact that our simple model ignores many confounders.\n",
"\n",
"**Proper modelling, using knowledge from epidemiology and sociology, is necessary to identify confounders.**"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimand type: nonparametric-ate\n",
"### Estimand : 1\n",
"Estimand name: backdoor\n",
"Estimand expression:\n",
" d \n",
"──────────────────────────────────────(Expectation(deaths|hospital beds))\n",
"d[Introduction of quarantine policies] \n",
"Estimand assumption 1, Unconfoundedness: If U→{Introduction of quarantine policies} and U→deaths then P(deaths|Introduction of quarantine policies,hospital beds,U) = P(deaths|Introduction of quarantine policies,hospital beds)\n",
"### Estimand : 2\n",
"Estimand name: iv\n",
"No such variable found!\n",
"\n"
]
}
],
"source": [
"print(identified_estimand)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*dowhy* allowed us to discover an estimand based on the *backdoor* criterion to assess the causal effect of the treatment on the outcome; notice that *dowhy* reminds us that this estimand is based on an *unconfoundedness* assumption."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Estimation\n",
"\n",
"We run three methods for the estimation of the causal effect: a *propensity score matching* (by default, *dowhy* uses LogisticRegression for computing propensity scores and k-NN (k=1) to find a match), *propensity score stratification* (by default, *dowhy* sets up 50 clipped strata), and *inverse probability weighting*."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:dowhy.causal_estimator:INFO: Using Propensity Score Matching Estimator\n",
"INFO:dowhy.causal_estimator:b: deaths~Introduction of quarantine policies+hospital beds\n",
"/home/fmzennaro/miniconda2_1/envs/covid/lib/python3.7/site-packages/sklearn/utils/validation.py:760: 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(deaths|hospital beds))\n",
"d[Introduction of quarantine policies] \n",
"Estimand assumption 1, Unconfoundedness: If U→{Introduction of quarantine policies} and U→deaths then P(deaths|Introduction of quarantine policies,hospital beds,U) = P(deaths|Introduction of quarantine policies,hospital beds)\n",
"### Estimand : 2\n",
"Estimand name: iv\n",
"No such variable found!\n",
"\n",
"## Realized estimand\n",
"b: deaths~Introduction of quarantine policies+hospital beds\n",
"## Estimate\n",
"Value: 401.1329164185031\n",
"\n"
]
}
],
"source": [
"estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.propensity_score_matching\")\n",
"print(estimate)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:dowhy.causal_estimator:INFO: Using Propensity Score Stratification Estimator\n",
"INFO:dowhy.causal_estimator:b: deaths~Introduction of quarantine policies+hospital beds\n",
"/home/fmzennaro/miniconda2_1/envs/covid/lib/python3.7/site-packages/sklearn/utils/validation.py:760: 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(deaths|hospital beds))\n",
"d[Introduction of quarantine policies] \n",
"Estimand assumption 1, Unconfoundedness: If U→{Introduction of quarantine policies} and U→deaths then P(deaths|Introduction of quarantine policies,hospital beds,U) = P(deaths|Introduction of quarantine policies,hospital beds)\n",
"### Estimand : 2\n",
"Estimand name: iv\n",
"No such variable found!\n",
"\n",
"## Realized estimand\n",
"b: deaths~Introduction of quarantine policies+hospital beds\n",
"## Estimate\n",
"Value: 1667.7588405151612\n",
"\n"
]
}
],
"source": [
"estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.propensity_score_stratification\")\n",
"print(estimate)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:dowhy.causal_estimator:INFO: Using Propensity Score Weighting Estimator\n",
"INFO:dowhy.causal_estimator:b: deaths~Introduction of quarantine policies+hospital beds\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(deaths|hospital beds))\n",
"d[Introduction of quarantine policies] \n",
"Estimand assumption 1, Unconfoundedness: If U→{Introduction of quarantine policies} and U→deaths then P(deaths|Introduction of quarantine policies,hospital beds,U) = P(deaths|Introduction of quarantine policies,hospital beds)\n",
"### Estimand : 2\n",
"Estimand name: iv\n",
"No such variable found!\n",
"\n",
"## Realized estimand\n",
"b: deaths~Introduction of quarantine policies+hospital beds\n",
"## Estimate\n",
"Value: 1440.5259837012877\n",
"\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/fmzennaro/miniconda2_1/envs/covid/lib/python3.7/site-packages/sklearn/utils/validation.py:760: 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": [
"estimate = model.estimate_effect(identified_estimand, method_name=\"backdoor.propensity_score_weighting\", method_params={\"weighting_scheme\":\"ips_weight\"})\n",
"print(estimate)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All the models seems to agree on a strong effect of the treatment on the outcome, although in the opposite way we would expect. This outcome is reasonable given our simple model and limited data: the data cover only a limited span of time; very likely the number of deaths is following its trend and the introduction of a policy like a quarantine which has a delayed effect was not registered in the data. Again, this should remark the current **NON-relevance** of this preliminary results."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Refutation\n",
"\n",
"A last step in a rigorous causal analysis would be a refutation step. Our model is just a dummy model, so we are not really taking its conclusions in consideration. However we present the formal refutation procedure provided by *dowhy*.\n",
"\n",
"We run a refutation process based on shuffling the values of our treatment; if the treatment had a real causal effect the estimate average treatment effect is expected to drop to 0."
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:dowhy.causal_estimator:INFO: Using Propensity Score Stratification Estimator\n",
"INFO:dowhy.causal_estimator:b: deaths~placebo+hospital beds\n",
"/home/fmzennaro/miniconda2_1/envs/covid/lib/python3.7/site-packages/sklearn/utils/validation.py:760: 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:(1667.7588405151612,)\n",
"New effect:(26.190793359377157,)\n",
"\n"
]
}
],
"source": [
"from dowhy.causal_refuters.placebo_treatment_refuter import PlaceboTreatmentRefuter\n",
"refuter = PlaceboTreatmentRefuter(data=data, identified_estimand=identified_estimand, estimate=estimate, placebo_type='permute')\n",
"print(refuter.refute_estimate())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Further work\n",
"\n",
"This notebook has presented an outline for a static causal analysis of covid19 data. However this draft is far from being complete:\n",
"\n",
"- **Proper modelling and inclusion of meaningful and relevant variables is essential.**\n",
"- **Validation of standard causal hypothesis and assumptions is required.**\n",
"- **A time-dependent model would be in order to process this scenario.**"
]
}
],
"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.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}