{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import sys\n",
"\n",
"sys.path.append(\"..\")\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import seaborn as sn\n",
"import utils\n",
"from sklearn.preprocessing import StandardScaler\n",
"\n",
"plt.style.use(\"ggplot\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tiltaksovervakingen: opsjon for kvalitetskontroll av analysedata\n",
"## Notebook 3: Outlier detection for whole water samples\n",
"\n",
"Exploring distributions for **single parameters** (as in notebook 2) is a reasonable starting point for quality assurance, but a more general approach is to look for \"outliers\" at the **water sample** level i.e. samples that are of questionable quality because *one or more* parameter values are unusual. If the suite of water quality parameters analysed is consistent (i.e. no data gaps), then each sample can be considered as a point in $n$ dimensional space, where $n$ is the number of parameters measured. Rather than looking for \"outliers\" along a single dimension (as in the distribution plots considered already), we can instead look for \"outliers\" in this higher dimensional space. A variety of algorithms are available to do this, some of which are explored below. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Terminology: \"outlier\" versus \"novelty\" detection\n",
"\n",
"\"**Outlier**\" and \"**novelty**\" detection are two different kinds of **anomaly** detection i.e. where we are interested in detecting abnormal or unusual observations. The difference between the two is important, but often overlooked:\n",
"\n",
" * **Outlier detection:** We have a **single dataset** that is believed to contain \"outliers\", which are observations that are \"far\" from the others (for some chosen definition of \"far\"). Outlier detection estimators thus try to fit the regions where the data is most concentrated, ignoring the deviant observations.\n",
"\n",
" * **Novelty detection:** We have access to a **reference dataset** that is *not* polluted by outliers, and we are interested in detecting whether a new observation (from a second dataset) is an outlier - a \"novelty\" - or not\n",
" \n",
"In the context of this project, we are primarily interested in **novelty detection**, because we have a reference historic dataset extracted from Vannmiljø and we would like to gauge whether observations in the \"new\" dataset are sufficiently unusual/unlikely to warrant further investigation and reanalysis. However, the [summary in the previous notebook](https://nbviewer.jupyter.org/github/NIVANorge/tiltaksovervakingen/blob/master/notebooks/02_distribution_plots.ipynb#3.-Summary) identified some possible issues in the Vannmiljø data, as well as in the \"new\" data. I am therefore reluctant to use the historic Vannmiljø dataset as a \"reference\" without additional cleaning. Instead, to begin with, at least, I will combine the \"new\" and \"historic\" results into a single dataset and then perform **outlier detection** (not novelty detection). This can be revised later if desired."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Read data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Choose dataset to process\n",
"lab = \"Eurofins\"\n",
"year = 2022\n",
"qtr = 3\n",
"version = 1"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"fold_path = f\"../../output/{lab.lower()}_{year}_q{qtr}_v{version}\""
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
vannmiljo_code
\n",
"
sample_date
\n",
"
lab
\n",
"
period
\n",
"
depth1
\n",
"
depth2
\n",
"
ALK_mmol/l
\n",
"
ANC_µekv/l
\n",
"
CA_mg/l
\n",
"
CL_mg/l
\n",
"
...
\n",
"
N-NO3_µg/l N
\n",
"
N-TOT_µg/l N
\n",
"
NA_mg/l
\n",
"
P-TOT_µg/l P
\n",
"
PH_<ubenevnt>
\n",
"
RAL_µg/l Al
\n",
"
SIO2_µg/l Si
\n",
"
SO4_mg/l
\n",
"
TEMP_°C
\n",
"
TOC_mg/l C
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
002-105961
\n",
"
2022-07-06
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
0.03
\n",
"
95.0
\n",
"
2.1
\n",
"
1.5
\n",
"
...
\n",
"
5.0
\n",
"
260.0
\n",
"
0.47
\n",
"
17.0
\n",
"
5.5
\n",
"
67.0
\n",
"
701.314913
\n",
"
0.40
\n",
"
NaN
\n",
"
19.0
\n",
"
\n",
"
\n",
"
1
\n",
"
002-105961
\n",
"
2022-07-19
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
NaN
\n",
"
NaN
\n",
"
2.5
\n",
"
NaN
\n",
"
...
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
5.9
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
15.0
\n",
"
NaN
\n",
"
\n",
"
\n",
"
2
\n",
"
002-105961
\n",
"
2022-08-03
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
0.06
\n",
"
130.0
\n",
"
2.7
\n",
"
1.4
\n",
"
...
\n",
"
5.0
\n",
"
250.0
\n",
"
0.35
\n",
"
17.0
\n",
"
6.2
\n",
"
48.0
\n",
"
607.806258
\n",
"
0.27
\n",
"
15.0
\n",
"
16.0
\n",
"
\n",
"
\n",
"
3
\n",
"
002-105961
\n",
"
2022-08-12
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
0.04
\n",
"
110.0
\n",
"
2.3
\n",
"
1.5
\n",
"
...
\n",
"
5.0
\n",
"
330.0
\n",
"
0.48
\n",
"
17.0
\n",
"
5.6
\n",
"
61.0
\n",
"
654.560586
\n",
"
0.32
\n",
"
NaN
\n",
"
21.0
\n",
"
\n",
"
\n",
"
4
\n",
"
002-105961
\n",
"
2022-08-30
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
NaN
\n",
"
NaN
\n",
"
2.6
\n",
"
NaN
\n",
"
...
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
5.7
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
13.0
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
"
5 rows × 25 columns
\n",
"
"
],
"text/plain": [
" vannmiljo_code sample_date lab period depth1 depth2 ALK_mmol/l \\\n",
"0 002-105961 2022-07-06 Eurofins new 0.0 0.0 0.03 \n",
"1 002-105961 2022-07-19 Eurofins new 0.0 0.0 NaN \n",
"2 002-105961 2022-08-03 Eurofins new 0.0 0.0 0.06 \n",
"3 002-105961 2022-08-12 Eurofins new 0.0 0.0 0.04 \n",
"4 002-105961 2022-08-30 Eurofins new 0.0 0.0 NaN \n",
"\n",
" ANC_µekv/l CA_mg/l CL_mg/l ... N-NO3_µg/l N N-TOT_µg/l N NA_mg/l \\\n",
"0 95.0 2.1 1.5 ... 5.0 260.0 0.47 \n",
"1 NaN 2.5 NaN ... NaN NaN NaN \n",
"2 130.0 2.7 1.4 ... 5.0 250.0 0.35 \n",
"3 110.0 2.3 1.5 ... 5.0 330.0 0.48 \n",
"4 NaN 2.6 NaN ... NaN NaN NaN \n",
"\n",
" P-TOT_µg/l P PH_ RAL_µg/l Al SIO2_µg/l Si SO4_mg/l TEMP_°C \\\n",
"0 17.0 5.5 67.0 701.314913 0.40 NaN \n",
"1 NaN 5.9 NaN NaN NaN 15.0 \n",
"2 17.0 6.2 48.0 607.806258 0.27 15.0 \n",
"3 17.0 5.6 61.0 654.560586 0.32 NaN \n",
"4 NaN 5.7 NaN NaN NaN 13.0 \n",
"\n",
" TOC_mg/l C \n",
"0 19.0 \n",
"1 NaN \n",
"2 16.0 \n",
"3 21.0 \n",
"4 NaN \n",
"\n",
"[5 rows x 25 columns]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Read from SQLite\n",
"stn_df, df = utils.read_data_from_sqlite(lab, year, qtr, version)\n",
"\n",
"# # Subset data to just the quarter of interest\n",
"# months_dict = {\n",
"# \"q1\": [1, 2, 3],\n",
"# \"q2\": [4, 5, 6],\n",
"# \"q3\": [7, 8, 9],\n",
"# \"q4\": [10, 11, 12],\n",
"# }\n",
"# months = months_dict[qtr]\n",
"# df = df[df[\"sample_date\"].dt.month.isin(months)]\n",
"\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Select parameters\n",
"\n",
"In order to perform outlier detection in a multi-dimensional space, it is necessary that all water samples have a **complete set of values for all parameters**. This is because outlier detection algorithms work by calculating distance metrics between samples, and this is not possible if a sample can't be located along one or more of the dimension axes due to missing values. It is therefore necessary to choose a set of parameters where the data are complete.\n",
"\n",
"The code below calculates the percentage of the time that each lab measures each parameter, where the percentage is calculated as:\n",
"\n",
" 100 * number_of_samples_for_par_X_from_lab_Y / total_number_of_samples_from_lab_Y"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plot\n",
"ax = pct_df.T.plot.bar(figsize=(15, 5))\n",
"ax.set_ylabel(\"Percent of samples\")\n",
"ax.set_title(\"Percentage of total samples analysed per parameter, split by lab\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Isolation Forests\n",
"\n",
"Isolation Forests are a type of random forest well suited to outlier detection. They have the advantage of making few assumptions about the underlying data distribution, which is useful in situations where - as here - most/all of the variables are strongly skewed.\n",
"\n",
"Isolation forests have a `contamination` parameter, which can be broadly interpreted as the \"expected proportion of outliers in the dataset\". In other words, setting `contamination=0.01` roughly translates to finding the most unusual 1% of data values. Without a strong theoretical basis for setting the `contamination` parameter, it must be found either by manual tuning or be fixed based on practical considerations (e.g. how many water samples can we realistically afford to reanalyse).\n",
"\n",
"### 4.1. `CA` and `PH` only\n",
"\n",
"The code below applies the isolation forest algorithm to `CA` and `PH`. Since these are almost always measured, this includes virtually all water samples (both new and historic) in the dataset."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/opt/conda/lib/python3.10/site-packages/sklearn/base.py:450: UserWarning: X does not have valid feature names, but IsolationForest was fitted with feature names\n",
" warnings.warn(\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"The total number of samples in the dataset is: 34852.\n",
"\n",
"The total number of outliers detected is 349:\n",
" 320 in the 'historic' period\n",
" 29 in the 'new' period\n",
"\n"
]
},
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
vannmiljo_code
\n",
"
sample_date
\n",
"
lab
\n",
"
period
\n",
"
depth1
\n",
"
depth2
\n",
"
CA_mg/l
\n",
"
PH_<ubenevnt>
\n",
"
pred
\n",
"
\n",
" \n",
" \n",
"
\n",
"
3930
\n",
"
021-28450
\n",
"
2022-07-08
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
7.7
\n",
"
6.7
\n",
"
outlier
\n",
"
\n",
"
\n",
"
3931
\n",
"
021-28450
\n",
"
2022-07-28
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
7.6
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
3933
\n",
"
021-28450
\n",
"
2022-08-17
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
12.0
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
3935
\n",
"
021-28450
\n",
"
2022-09-09
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
8.6
\n",
"
6.6
\n",
"
outlier
\n",
"
\n",
"
\n",
"
4871
\n",
"
021-46373
\n",
"
2022-07-05
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
6.1
\n",
"
7.2
\n",
"
outlier
\n",
"
\n",
"
\n",
"
4872
\n",
"
021-46373
\n",
"
2022-07-19
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
7.9
\n",
"
7.4
\n",
"
outlier
\n",
"
\n",
"
\n",
"
4873
\n",
"
021-46373
\n",
"
2022-08-11
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
7.4
\n",
"
7.4
\n",
"
outlier
\n",
"
\n",
"
\n",
"
5752
\n",
"
022-32018
\n",
"
2022-09-23
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
9.6
\n",
"
7.5
\n",
"
outlier
\n",
"
\n",
"
\n",
"
5847
\n",
"
022-32019
\n",
"
2022-08-01
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.4
\n",
"
7.0
\n",
"
outlier
\n",
"
\n",
"
\n",
"
5849
\n",
"
022-32019
\n",
"
2022-09-23
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
7.4
\n",
"
6.9
\n",
"
outlier
\n",
"
\n",
"
\n",
"
5940
\n",
"
022-32020
\n",
"
2022-07-04
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
7.4
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
5941
\n",
"
022-32020
\n",
"
2022-07-22
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
6.0
\n",
"
7.0
\n",
"
outlier
\n",
"
\n",
"
\n",
"
5942
\n",
"
022-32020
\n",
"
2022-07-26
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
6.2
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
5943
\n",
"
022-32020
\n",
"
2022-08-01
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
6.5
\n",
"
6.7
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6243
\n",
"
022-45769
\n",
"
2022-07-19
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
6.5
\n",
"
6.6
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6245
\n",
"
022-45769
\n",
"
2022-08-17
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.9
\n",
"
6.5
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6246
\n",
"
022-45769
\n",
"
2022-08-29
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.7
\n",
"
6.5
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6762
\n",
"
022-58904
\n",
"
2022-07-04
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
6.0
\n",
"
6.7
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6763
\n",
"
022-58904
\n",
"
2022-07-22
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.7
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6764
\n",
"
022-58904
\n",
"
2022-07-26
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.4
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6765
\n",
"
022-58904
\n",
"
2022-08-01
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.6
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
6770
\n",
"
022-58904
\n",
"
2022-09-29
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.8
\n",
"
6.8
\n",
"
outlier
\n",
"
\n",
"
\n",
"
14226
\n",
"
027-79278
\n",
"
2022-07-05
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
4.8
\n",
"
7.1
\n",
"
outlier
\n",
"
\n",
"
\n",
"
14227
\n",
"
027-79278
\n",
"
2022-08-02
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
4.9
\n",
"
7.2
\n",
"
outlier
\n",
"
\n",
"
\n",
"
17196
\n",
"
030-58838
\n",
"
2022-08-16
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
7.3
\n",
"
6.2
\n",
"
outlier
\n",
"
\n",
"
\n",
"
17197
\n",
"
030-58838
\n",
"
2022-08-30
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
5.6
\n",
"
6.4
\n",
"
outlier
\n",
"
\n",
"
\n",
"
19593
\n",
"
036-58751
\n",
"
2022-09-07
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
11.0
\n",
"
7.9
\n",
"
outlier
\n",
"
\n",
"
\n",
"
31621
\n",
"
079-58878
\n",
"
2022-09-13
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
6.6
\n",
"
6.2
\n",
"
outlier
\n",
"
\n",
"
\n",
"
33482
\n",
"
082-58870
\n",
"
2022-09-06
\n",
"
Eurofins
\n",
"
new
\n",
"
0.0
\n",
"
0.0
\n",
"
10.0
\n",
"
7.3
\n",
"
outlier
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" vannmiljo_code sample_date lab period depth1 depth2 CA_mg/l \\\n",
"3930 021-28450 2022-07-08 Eurofins new 0.0 0.0 7.7 \n",
"3931 021-28450 2022-07-28 Eurofins new 0.0 0.0 7.6 \n",
"3933 021-28450 2022-08-17 Eurofins new 0.0 0.0 12.0 \n",
"3935 021-28450 2022-09-09 Eurofins new 0.0 0.0 8.6 \n",
"4871 021-46373 2022-07-05 Eurofins new 0.0 0.0 6.1 \n",
"4872 021-46373 2022-07-19 Eurofins new 0.0 0.0 7.9 \n",
"4873 021-46373 2022-08-11 Eurofins new 0.0 0.0 7.4 \n",
"5752 022-32018 2022-09-23 Eurofins new 0.0 0.0 9.6 \n",
"5847 022-32019 2022-08-01 Eurofins new 0.0 0.0 5.4 \n",
"5849 022-32019 2022-09-23 Eurofins new 0.0 0.0 7.4 \n",
"5940 022-32020 2022-07-04 Eurofins new 0.0 0.0 7.4 \n",
"5941 022-32020 2022-07-22 Eurofins new 0.0 0.0 6.0 \n",
"5942 022-32020 2022-07-26 Eurofins new 0.0 0.0 6.2 \n",
"5943 022-32020 2022-08-01 Eurofins new 0.0 0.0 6.5 \n",
"6243 022-45769 2022-07-19 Eurofins new 0.0 0.0 6.5 \n",
"6245 022-45769 2022-08-17 Eurofins new 0.0 0.0 5.9 \n",
"6246 022-45769 2022-08-29 Eurofins new 0.0 0.0 5.7 \n",
"6762 022-58904 2022-07-04 Eurofins new 0.0 0.0 6.0 \n",
"6763 022-58904 2022-07-22 Eurofins new 0.0 0.0 5.7 \n",
"6764 022-58904 2022-07-26 Eurofins new 0.0 0.0 5.4 \n",
"6765 022-58904 2022-08-01 Eurofins new 0.0 0.0 5.6 \n",
"6770 022-58904 2022-09-29 Eurofins new 0.0 0.0 5.8 \n",
"14226 027-79278 2022-07-05 Eurofins new 0.0 0.0 4.8 \n",
"14227 027-79278 2022-08-02 Eurofins new 0.0 0.0 4.9 \n",
"17196 030-58838 2022-08-16 Eurofins new 0.0 0.0 7.3 \n",
"17197 030-58838 2022-08-30 Eurofins new 0.0 0.0 5.6 \n",
"19593 036-58751 2022-09-07 Eurofins new 0.0 0.0 11.0 \n",
"31621 079-58878 2022-09-13 Eurofins new 0.0 0.0 6.6 \n",
"33482 082-58870 2022-09-06 Eurofins new 0.0 0.0 10.0 \n",
"\n",
" PH_ pred \n",
"3930 6.7 outlier \n",
"3931 6.8 outlier \n",
"3933 6.8 outlier \n",
"3935 6.6 outlier \n",
"4871 7.2 outlier \n",
"4872 7.4 outlier \n",
"4873 7.4 outlier \n",
"5752 7.5 outlier \n",
"5847 7.0 outlier \n",
"5849 6.9 outlier \n",
"5940 6.8 outlier \n",
"5941 7.0 outlier \n",
"5942 6.8 outlier \n",
"5943 6.7 outlier \n",
"6243 6.6 outlier \n",
"6245 6.5 outlier \n",
"6246 6.5 outlier \n",
"6762 6.7 outlier \n",
"6763 6.8 outlier \n",
"6764 6.8 outlier \n",
"6765 6.8 outlier \n",
"6770 6.8 outlier \n",
"14226 7.1 outlier \n",
"14227 7.2 outlier \n",
"17196 6.2 outlier \n",
"17197 6.4 outlier \n",
"19593 7.9 outlier \n",
"31621 6.2 outlier \n",
"33482 7.3 outlier "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Columns of interest\n",
"key_cols = [\"vannmiljo_code\", \"sample_date\", \"lab\", \"period\", \"depth1\", \"depth2\"]\n",
"par_cols = [\"CA_mg/l\", \"PH_\"]\n",
"\n",
"# Run algorithm\n",
"data = df[key_cols + par_cols].dropna()\n",
"data = utils.isolation_forest(data, par_cols, contamination=0.01)\n",
"\n",
"# Summarise results\n",
"all_out = data.query(\"pred == 'outlier'\")\n",
"his_out = data.query(\"(pred == 'outlier') and (period == 'historic')\")\n",
"new_out = data.query(\"(pred == 'outlier') and (period == 'new')\")\n",
"\n",
"csv_path = os.path.join(fold_path, \"isoforest_ca_ph.csv\")\n",
"new_out.to_csv(csv_path, index=False)\n",
"\n",
"print(f\"The total number of samples in the dataset is: {len(data)}.\\n\")\n",
"print(\n",
" f\"The total number of outliers detected is {len(all_out)}:\\n\"\n",
" f\" {len(his_out)} in the 'historic' period\\n\"\n",
" f\" {len(new_out)} in the 'new' period\\n\"\n",
")\n",
"\n",
"new_out"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This initial approach identifies the **strangest 1% of the dataset overall**. \n",
"\n",
" * Most of the unusual values (300 out of 326) are actually in the historic dataset (i.e. they are already in Vannmiljø). \n",
" \n",
" * The 26 samples in the 'new' dataset that have been classified as outliers are predominantly those with unusually high concentrations of `CA` (greater than around 5 mg/l; see plots below)."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAImCAYAAABuJeE8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACa3ElEQVR4nOzdd3hUVf7H8fdt09MJJJTQi0hXActaUCyoIMIqioprWQsq9oauDUXUtay6NnRVdNWfClZUVLCArlhBQUFAitT0ZPrMvff3R2AgJkEIITNJvq/n8ZHcc2fmTM5k5jPnnqLYtm0jhBBCCJEEarIrIIQQQoiWS4KIEEIIIZJGgogQQgghkkaCiBBCCCGSRoKIEEIIIZJGgogQQgghkkaCiBBCCCGSRoKIEEIIIZJGgogQQgghkkaCiBCiQd16660oisInn3xS7biiKBx++OEN8hidOnWiU6dODXJfQojkkiAihNgld955J4qioCgKy5YtS3Z1GkRDhiMhRP1IEBFC/Cnbtnn66adRFAWA6dOnJ7lGQojmQoKIEOJPzZkzh99++43zzjuP1q1b89xzzxGNRpNdLSFEMyBBRAjxp5566ikAzjvvPMaPH09hYSFvvPHGXn1M27Z55JFH2HfffXG5XLRr146JEydSXl5e6/nl5eXce++9DBs2jPbt2+NwOMjNzWXkyJF88cUX1c599tlnE707n376aeKSk6Io3HrrrdXOGzNmDF26dMHtdpOens7BBx/M888/v9eetxAtjZ7sCgghUtvmzZt566232GeffRg8eDBut5sHHniAJ598klNOOWWvPe7ll1/Ov/71L/Lz8/n73/+OYRi8+eabLFy4kGg0isPhqHb+zz//zOTJkzn00EM5/vjjycrKYs2aNbz55pvMnj2bt956ixEjRgAwYMAAbrnlFm677TY6duzI2WefnbifHceMXHTRRfTu3ZtDDz2U/Px8ioqKePfdd5kwYQK//PILd9111157/kK0GLYQQuzE1KlTbcC+++67E8cGDhxoK4pir1y5ssb5t9xyiw3Y8+bNq3YcsA877LBdeswFCxbYgN21a1e7uLg4cTwUCtlDhw61Abtjx47VblNWVmYXFhbWuK/Vq1fbbdq0sXv27Fmj7M/qtGLFihrHwuGwffjhh9u6rtvr1q3bpecjhKibXJoRQtTJtm2mT5+OpmmceeaZieN/+9vfEmV7w3/+8x8AJk+eTHZ2duK4y+Vi6tSptd4mIyODVq1a1TjesWNH/vrXv7Js2TLWrl27W/Xo2rVrjWNOp5NLLrmEeDzO3Llzd+v+hBA1SRARQtRp7ty5rFy5kuHDh9O2bdvE8dNPPx2Hw8F//vMf4vF4gz/ud999B8Bhhx1Wo+wvf/kLul77VeUFCxZwyimn0KFDB5xOZ2LcxyOPPALA+vXrd6sea9euZeLEifTq1QuPx5O4v7Fjx9br/oQQNckYESFEnZ588kmAamMoAHJycjjxxBN5/fXXefvttxk9enSDPu62Aalt2rSpUaZpGjk5OTWOz5o1i7Fjx+JyuRg+fDhdu3bF6/WiqiqffPIJn376KZFIZJfrsGrVKgYPHkxpaSl/+ctfOProo8nIyEDTNFavXs1zzz23W/cnhKidBBEhRK12nBkzbtw4xo0bV+t5Tz75ZIMHkYyMDKBqoGyXLl2qlZmmSXFxMe3atat2/Oabb8bhcPDNN9+wzz77VCu74IIL+PTTT3erDvfffz/FxcX85z//qRHEXnrpJZ577rnduj8hRO0kiAgharVtrZD99tuPAQMG1HrOm2++yZw5c1izZg0dO3ZssMceNGgQ3333HZ9++mmNIPL555/XejloxYoV7LvvvjVCiGVZzJ8/v9bHUVUV0zRrLVuxYgUAY8aMqVG2u6FGCFE3CSJCiFptG4j673//m8GDB9d6Tm5uLnfffTdPP/00t99+e4M99tlnn8306dO58847GTVqVGLAajgc5oYbbqj1Np06deLXX39l/fr1id4S27a57bbbWLp0aa23ycnJYd26dXXeH8C8efMYOXJk4vgHH3wgK8sK0YBksKoQooZPPvmEZcuW0bdv3zpDCFQtcKYoCs8880ydPQv1cfDBB3PppZeycuVK+vTpw2WXXcZVV11Fnz59iMfj5Ofn17jNFVdcQWVlJYMGDeLiiy9m0qRJHHDAAdx7772ceOKJtT7OkUceyZo1axg1ahS33norU6ZM4bPPPgPg4osvxuFwcMoppzB+/HiuvfZaRowYwXHHHZcYrCqE2HMSRIQQNey4kurOdO3alcMPP5z169cze/bsBq3DQw89xMMPP0xGRgZPPPEEL730EscccwwfffRRjcXMoGocyH/+8x/y8/N57rnnePHFF+nQoQNfffUVgwYNqvMxTjvtNP73v/9xxx13cPPNNyem5Pbr14958+Zx0EEHMXv2bB577DEqKiqYOXMmF154YYM+VyFaMsW2bTvZlRBCCCFEyyQ9IkIIIYRIGgkiQgghhEgaCSJCCCGESBoJIkIIIYRIGgkiQgghhEgaCSJCCCGESBoJIkIIIYRIGgkiQgghhEga2WumFqWlpbVuqpXqcnNzKSwsTHY1RC2kbVKXtE3qkrZJXbW1ja7rZGVl7fZ9SRCpRTweJxaLJbsau0VRFKCq7rJYbmqRtkld0japS9omdTV028ilGSGEEEIkjQQRIYQQQiSNBBEhhBBCJI2MERFCCNEkxeNxgsFgsqvRong8HgzDaND7lCAihBCiyYnH4wQCAdLS0lBV6dxvDJZlUVlZidfrbdD7TYkg8n//93+89tpr1Y5lZGTw1FNPAWDbNq+++ioff/wxfr+f7t27c+6559KhQ4fE+bFYjBkzZrBgwQKi0Sh9+vThvPPOIycnp1GfixBCiL0vGAxKCGlkqqqSlpaG3+9v0PtNiSAC0KFDB26++ebEzzu+uN58803effddLr74YvLz85k5cyZTpkzhwQcfxO12A/Dss8/y7bffMmnSJNLS0nj++ee5++67mTZtmrxQhRCiGZL39sa3N37nKdOKqqqSmZmZ+C89PR2o6g2ZPXs2o0ePZsiQIRQUFDBx4kQikQjz588HqpLx3LlzOeuss+jXrx+dO3fm0ksvZe3atSxevDiZT0sIIYQQO5EyPSKbNm3iggsuQNd1unfvzmmnnUabNm3YsmULZWVl9O/fP3GuYRj07t2bZcuWMXz4cFatWoVpmvTr1y9xTnZ2NgUFBSxfvpwBAwbU+pixWKzawmWKoiR6WLYt2NJUbKtvU6t3SyBtk7qkbVKXtE3qa6i2SYkg0r17dyZOnEjbtm0pKytj5syZ3HTTTdx///2UlZUBVWNGdpSRkUFRUREAZWVl6LqOz+ercc6229dm1qxZ1camdO7cmWnTppGbm9swTywJ8vLykl0FUQdpm9QlbZO66mqbUCjU4LM3xK5xOBxAw/3dpEQQGThwYOLfBQUF9OjRg0svvZRPP/2U7t27AzWT164sK/tn54wePZoTTjgh8fO2xygsLGxye80oikJeXh6bNm2S5ZBTjLRN6pK2SV1/1jbRaLTJbcWxJ4YMGcJ5553H+eefn+yqEI1GAWq0ja7r9foinxJB5I9cLhcFBQVs3LiRAw44AKjq9dhxM52KiopEL0lmZibxeBy/31+tV6SiooKePXvW+TiGYdSZqJvqm5Jt20227s2dtE3qkrZJXdI2qauh2iZlBqvuKBaLsX79erKysmjdujWZmZnVBp3G43GWLl2aCBldunRB07Rq55SWlrJ27Vp69OjR6PUXQggh/mhbT4KoLiV6RJ5//nn2339/WrVqRXl5Oa+//jqhUIjDDjsMRVEYMWIEs2bNIj8/n7y8PGbNmoXT6eSQQw4BqlZ6GzZsGDNmzCAtLQ2fz8eMGTMoKCioNoBVCNHy6LqOLxZBCVRiWyakZRDQncRMM9lVE03c2LFjE1+IZ86ciaqqnHXWWVx77bUoisKQIUM47bTTWL16Ne+//z7HHHMMDz30EF9//TVTp05l0aJFZGVlcdxxx3HDDTfg8XgAKCoq4qqrrmL+/Pnk5uZy7bXXJvNp7nUpEURKSkp46KGHqKioID09ne7du3PnnXcmrjWNGjWKaDTK9OnTCQQCdOvWjcmTJydmuABMmDABTdN44IEHEguaXXfddTLPXIgWzKGqeDaspvS+m7FKiwFQvD4yLrwWfZ8BhJNcP9H0vfrqq4wbN463336bxYsXc+2119K+fXvGjx8PwOOPP87ll1/OpEmTAPj5558ZP34811xzDffddx/FxcXcdNNNTJ48mQceeACAK664gg0bNvDKK6/gcDi4+eabE5MzmiPFlotvNRQWFja5QVCKopCfn8/GjRvlemqKkbZJnuxYiMLLxkMtg89z7plORU4eeXl50jYp6M/+brZ9cU2msWPHUlRUxLx58xKTHe666y7mzJnDJ598wpAhQ+jTpw9PP/104jaXXXYZLpeLe+65J3Fs4cKFjBkzhl9//ZX169dz6KGH8vbbbzNo0CAAVqxYwWGHHcatt96aEoNVKyoq6NWrV422MQyj+QxWFUKIPeU0DILvvlJrCAHwvzwd98TJjVwr0dwMGjSo2qzO/fbbjyeeeAJz66W/Pw4P+PHHH1m9ejWzZs1KHLNtG8uyWLduHatWrULX9WprZ3Xr1q3GEhbNiQQRIUSzpFkm0VXL6iyPr/sNj9m0pumLpmfbuI9tLMvijDPO4Jxzzqlxbrt27Vi5ciXQshZykyAihGiWTFVD79SNyOJvai3X23bE1OQtUOyZ7777rsbPnTt3RtO0Ws/v27cvy5Yto3PnzrWWd+vWjXg8zqJFixJrbK1YsYLy8vKGrXgKkZGcQohmKRKL4T56FKi1fyD4xp1LWJG3QLFnNmzYwK233sqKFSt44403eOaZZzj33HPrPP/iiy/m22+/5cYbb+Snn35i1apVzJkzh5tuugmoCiJHHHEE11xzDd999x2LFy/mmmuuweVyNdZTanTydUAI0WwF3Wlk/+OflP3zFqzKqm+UistN+nlXEM3Nx7KsJNdQNHVjx44lHA5zwgknoGka55xzDmeccUad5/fu3ZvXX3+dadOmcfLJJ2PbNh07dmTkyJGJc+6//36uvvpqxo4dS6tWrbj22mvZsGFDYzydpJBZM7WQWTOiIUnbJJeua1XriFSWY1sWSnomAcNJ1LSkbVJYU5k107t3b26//fak1qOxyawZIYTYDfG4SZmiQ3rO9oOm9IQIkSrkAqkQQgghkkZ6RIQQQoh6eO2115JdhWZBekSEEEIIkTQSRIQQQgiRNBJEhBBCCJE0EkSEEEIIkTQSRIQQQgiRNDJrRgghdsKp67ijIYhGwOEk7HATrmNHXyHE7pMeESGEqEOGAtqcmZRcfibFl55OyWWnw5svkKnKKqwi9bRr1473338fgHXr1tGuXTt++umnJNfqz0kQEUKIWrg1lcjbL1P50nTsUBAAOxIh8MZ/Cc54DG/L2aW9RWhKS/z/85//ZPjw4Ts9p23btnz//ff06tWrkWpVf3JpRgghauGKhil69/9qLQt9+gHeU88hoDffHVFbAjscxJr1AixaCGYcNB36D0YdfQaKy5Ps6u0RTdNo3br1Ht1HNBrF4XA0UI3qJj0iQghRC9tfCXWNBbFtrNKSxq2QaFB2OIg19VqY9y4Ub4Gykqr/fzIba+q12OHgXnvsSCTCzTffTL9+/ejSpQsnnXQSP/zwAwCvvPIK++yzT7Xz33//fdq1a5cov//++1m6dCnt2rWjXbt2vPLKKzUeo7ZLM8uXL+fMM8+ke/fu9O/fn0svvZSSku2v47FjxzJ58mRuvfVW+vTpw2mnnbYXnn1NEkSEEKIWimvnvR2qx9tINRF7gzXrBdi4Dv54ScayYNPvVeV7yZ133sns2bN58MEHef/99+nUqRPjx4+ntLT0T287cuRILrjgAnr27Mn333/P999/z8iRI//0dps3b2bMmDH07t2b9957jxdffJGioiIuuOCCaue9+uqr6LrOG2+8wbRp0+r9HHeHXJoRQohaxF1eHD33JbpsSY0yvV1HTG9aEmolGsyihTVDyDaWVVV+2t8b/GGDwSDPP/88DzzwAMOGDQPg3nvvZejQobz88stkZ2fv9PZutxuv17vbl16ef/55+vbtyw033JA49s9//pMDDjiAlStX0rVrVwA6derETTfdVI9nVn/SIyKEELUIKCoZV9yG3rZDteNaqzZk3XgPAc1IUs3EnrJtu2pMyM6Y5l4ZwLp69WpisRgHHHBA4phhGAwYMIBff/21wR9vm8WLF/PFF1/QvXv3xH+HHXYYAGvWrEmc179//71Wh7pIj4gQollRFAXDMLBtm1gsVu/7sSyLMsNF+m0PoxRvJv77GrS89iit8ynXHZim2YC1Fo1JUZSqgak7o2lV5zWwbeHmj/dt2zaKoqCqao0AtCev4x3vf/jw4dx44401ytq0aZP4t9vt3uPH2l3SIyKEaDbSFJuMsiK02f+H8fGbZIX9eJT6f6u1LItyVMpbtSU06GAq8jpQpmgSQpqD/oNBreMjUFWryveCzp0743A4WLhwYeJYLBZj0aJFdO/enZycHPx+P8Hg9sGyS5ZUvzxoGAaWZe3W4/bp04dly5bRoUMHOnfuXO0/jye5M4QkiAghmoUMxSb89AMUX/03/K88Q+WMxym69HTsD9/CuwdhBKq+TZp7qateJIc6+gzIa18zjKgq5LWvKt8LPB4PZ555JlOmTGHevHksX76ca665hnA4zLhx4xg4cCBut5u7776b3377jVmzZvHqq69Wu48OHTqwdu1afvrpJ0pKSohEIn/6uGeffTZlZWVcfPHFfP/996xZs4ZPP/2UK6+8MunBWoKIEKLJ03Ud88dvCH/5SY0y/ytPY5QU7pVudtF0KS4P6g33wOEjIKc1ZOZU/f/wEag33LNX1xG58cYbGTFiBJdddhnHHnssq1ev5sUXXyQzM5OsrCwefvhhPv74Y4466ijeeOMNrrzyymq3HzFiBIcffjinnHIKffv25Y033vjTx8zLy+ONN97AsizGjx/PsGHD+Mc//kFaWhpqXT1DjUSxJeLXUFhY2CDX5BqToijk5+ezceNG+daWYqRt9r40LAJTriK2ZmWt5e4jjkWdMInwH/6upW1S15+1TUVFBenp6Q32eNvGaIg/V1FRQa9evWq0jWEY5Obm7vb9SY+IEKLJU2wby19RZ7lVVopi7941ddGySAhJHgkiQogmL2Y4cA4YUme5a+jhxBV5uxMiFclfphCiyQvFTbxjzkRxOmuUaTm5GAOHNrnLrUK0FBJEhBDNQoXTS869z+AcuLVnRNNwH34M2Xc9ToW+9zfuEkLUjyxoJoSoN1VVcWPjsOLYQEh3Eqlro7i9LG5ZlHszcF9yE754FFCIGA5KLRtb1v0QImVJEBFC1Iuhqvj8pVT+5xEqFi1EcbnxHHMSWcefQpmiJWUWimVZBADUrcuvmzJAVYhUJ5dmhBC7TVEUfP4yiq46h8gPX4FtY4eCBN74L+V3Xk2GLT0QQohdI0FECLHbPAr4Z/wb4jUHgMZWr8BauwpN05JQMyFEUyNBRAix2wwzTviHhXWWR76ch2HI7rRCiD8nQUQIsdtsQPV46yxXMzJlpVIhajF27Fj+8Y9/7NK569ato127dvz0008AfPHFF7Rr147y8vK9WcVGJ4NVhRC7LWQ48YwYi//lp2stdx12LKW7sBGXEKmisZZ4f+qpp+rdW7j//vvz/fffN+jS9qlAgogQYrdF43GyjhpJ5Lv/EVtefYvytLMvJeLxVXWbCJHCgjGTFxcVsvD3AHHLRlcVBrf3Mr5/Lh5j74xxysrKqvdtHQ4HrVu33qPHj0ajOBypta6OBBEhRL2UoZJ+zZ2w6Xci//sUxZeO65Ajibh9BG0SO3paVmpMoU21+ojkCsZMrvtgDevKo9Uy8+xlZSzeFGTaMR33ShgZO3YsvXv35vbbb2fIkCGMHz+e1atX884775CRkcGkSZM444wzar3tF198wV//+leWLl1KRkYGAF9//TVTp05l0aJFZGVlcdxxx3HDDTfg8VTtHjxkyBBOO+00Vq9ezfvvv88xxxzDQw891ODPa0/IGBEhRL3Ytk05KpVtO2Gfeh7miL9S6vJhKZBlRvGtWIJv1c9kmVFcavI2FHOqCllWDN9vy/D++hPZ8QjuJNZHpIYXFxXWCCEAFvB7eZQXFxU2Sj2eeOIJ+vXrxwcffMCECRO44YYbWLFixS7d9ueff2b8+PEcd9xxfPjhhzz22GMsXLiQyZMnVzvv8ccfp1evXrz33ntcfvnle+FZ7BnpERFC7BHLsohsHQ/iVWxYMJeiZx8Ba+taIppG+nlX4Bl8KMFG7oxwKwr6j99Q9Mhd26caqyq+U/6G76iRBJBA0lIt/D1Q59VDa2v5+fvv/XoMGzaMs88+G4CJEyfy1FNP8cUXX9CtW7c/ve1jjz3GSSedxPnnnw9Aly5duOOOOxgzZgxTp07F5XIBcPDBB3PhhRfuteewpySICCEahKqq6Jt/p+SZP3T7miYVT9xHTrd9CGW1brTZNIqi4PKXUfTgbdULLAv/y0+T1aMPaqcejVIXkVps2yZu7fx1aFp2owxg7d27d+LfiqKQm5tLcXHxLt32xx9/ZPXq1cyaNStxzLZtLMti3bp1dO/eHYB+/fo1bKUbmAQRIUSDcCkQnPl8neXBWS/iOu8qQo2074tD1wnOfq3O8sCr/8F99ZRGqYtILYqioP/J5TlNVRplFo2uV/8YVhRll8cxWZbFGWecwTnnnFOjrF27dol/bxsvkqokiAghGoRqxjELN9dZHt+yEcOKQyNdDtEsk+im9XWWm0WbUWUzvBZrcHsvs5eVUdtHvrq1PNX17duXZcuW0blz52RXZY/IYFUhRIOI6wZGz751ljt69SOuNd5qq3FVw9FnYJ3lRvfexHVZ/bWlGt8/l/YZjhofgirQPsPB+P65yajWbrn44ov59ttvufHGG/npp59YtWoVc+bM4aabbkp21XaLBBEhRIMIx028o8aBXktHq+HAM2IMkXi80eoTjcVwHXoMistds1DV8J1yDmFZ66TF8hga047pyIiembT2GuS4dVp7DUb0zNxrU3cbWu/evXn99df57bffOPnkkznmmGO455579nitkcam2LIOcw2FhYXEYjU380pliqKQn5/Pxo0bZWntFNOS2sahqXi2bKD84TuJr18DgN6hMxmX3UQgpw0xs3GnzeiaRlp5MeUPTyG2ajkAWl47MibeQLhtR6I2LaZtmpo/+7upqKho0BVGG2tl1eagoqKCXr161WgbwzDIzd39niQZIyKEaDBR08Js3Y60Wx9CDQdRANPtoVI1MHdxPIaqqriw0cw4pm4Q2jp7oT7ipklFRg7eG+9FD4ewLRPb7SWgO4jH4/LBIxLktZA8EkSEEA3KNE0qUMHlqzpgA7sYQjwKOLasx//KM8Q3/Y7RsRuZp/6NUEYO4T+Zbrmz+lSigHOHmQONeIlICLFzEkSEECnBqamo38yn+N93J46Zm9YTXvgZWTfeg9ltX2ISIIRodmSwqhAiJXhiESqmP1CzwLYpf+QuvDHZzVeI5kiCiBAiJdglRdjR2sOGVV4KgYpGrpFIZTK4OHka+ncvQUQIkRrUnb8dKYq8XYntdF0nEAhIIGlEtm0TCARqrAa7p2SMiBAiNWTloLi92KFAjSIttw2Wx5eESolU5fV6iUQiVFZWJrsqLYrT6UxsptdQJIgIIVJCQHeQefk/KL37etjxW66uk3nFrfgNp8x2EdU4nU6cTmeyqyH2kAQRIURKiJkWWrfetHroBYLvvkb8998wuvXGc8wo/C4vcQkhQjRLEkSEEA3K0DW8sQiKv6rL3Pal4TdcuxQkwpZN2J2G67S/4zTjmJpOSSwGjbwiqxCi8UgQEUI0GJcKjl+XUPrwnViV5QCoGVlkTvoH4Y7diOziomThbVssWE1rqwUhxO6TYehCiAahKAruyjJKp16XCCFQNfW2ZMpVeIIyqFAIUZMEESFEg3CpKoHXZ1QfaLqNZRF86yVceurvaCqEaFwSRIQQDUIzY8TXrqqzPPbbCjRTBpwKIaqTICKEaBCmbqC371hnuV7QGVOTYWlCiOokiAghGkTYtPCefFbthYqCd9TpRGT2ixDiDySICCEahG3bhDNzyLx6CorbmziueH1kXTeVkC9DluMWQtQg/aRCiAYTsmwcvQeS/eDzUFGOooCdnklAdxIzzWRXTwiRgiSICCEaVNQ0iaoGZLbaflBCiBCiDnJpRgghhBBJI0FECCGEEEkjQUQIIYQQSSNBRAghhBBJI0FECCGEEEkjQUQIIYQQSSNBRAghhBBJI+uICCGSQlEU3KqCMxrGjsWwnS6ChpNYXDbGE6IlSbkgMmvWLF566SVGjBjB2WefDVQtHf3qq6/y8ccf4/f76d69O+eeey4dOnRI3C4WizFjxgwWLFhANBqlT58+nHfeeeTk5CTpmQgh6qIoClm2if8/j1K5YC5YJmpWDulnTcTRdz8CtpLsKgohGklKXZpZsWIFH330ER07Vt/B88033+Tdd9/lnHPOYerUqWRmZjJlyhRCoVDinGeffZaFCxcyadIkbr/9dsLhMHfffTeWJZtsCZFq0m2TsjuvJvT5h2BVrbpqlRZT9tDtqD8vwjCMJNdQCNFYUiaIhMNhHn74YS644AK83u0bZtm2zezZsxk9ejRDhgyhoKCAiRMnEolEmD9/PgDBYJC5c+dy1lln0a9fPzp37syll17K2rVrWbx4cbKekhCiFoqiQNFmYqtX1Fpe+dyjeKLhRq6VECJZUubSzPTp0xk4cCD9+vVj5syZieNbtmyhrKyM/v37J44ZhkHv3r1ZtmwZw4cPZ9WqVZimSb9+/RLnZGdnU1BQwPLlyxkwYECtjxmLxYjFYomfFUXB7XYn/t2UbKtvU6t3SyBtU52u68RWLK2z3CzajBqLojg9e70u0japS9omdTV026REEFmwYAG//fYbU6dOrVFWVlYGQEZGRrXjGRkZFBUVJc7RdR2fz1fjnG23r82sWbN47bXXEj937tyZadOmkZubW89nknx5eXnJroKog7TNdsE2bessUxxODLeH/Db5jVYfaZvUJW2TuhqqbZIeRIqKinj22WeZPHkyDoejzvP+mLxs2/7T+/6zc0aPHs0JJ5xQ4zEKCwuJN7GR+4qikJeXx6ZNm3bpdyMaj7RNTZmduqO43NjhUI0yz/CRVKoa4Y0b93o9pG1Sl7RN6qqrbXRdr9cX+aQHkVWrVlFeXs7111+fOGZZFj///DPvv/8+Dz74IFDV65GVlZU4p6KiItFLkpmZSTwex+/3V+sVqaiooGfPnnU+tmEYdQ6Ka6ovfNu2m2zdmztpm+0qNQfZtz5EyW1XYIcCieOOPoNwjz6D0ljjfhGQtkkNTk3DEwtjV5ajGA7MEh1VVZvcF8OWoqH+bpIeRPr27ct9991X7dhjjz1G27ZtGTVqFG3atCEzM5PFixfTuXNnAOLxOEuXLmX8+PEAdOnSBU3TWLx4MQcddBAApaWlrF27NnGOECJ1xC0Lf+t2ZD/4PNb6NVilxRidexBPz6QMFZBQ0NKkKWDOe4fiV57BjkYA0NsWkHndXVSkZWGaZpJrKPaWpAcRt9tNQUFBtWNOp5O0tLTE8REjRjBr1izy8/PJy8tj1qxZOJ1ODjnkEAA8Hg/Dhg1jxowZpKWl4fP5mDFjBgUFBdUGsAohUkfcNClTDdSOPVA6KfhNc2v+kBDS0hiGAYu+onLGY9WOxzespXjyReT88zlK1aR/XIm9pEm07KhRo4hGo0yfPp1AIEC3bt2YPHlyYoYLwIQJE9A0jQceeCCxoNl1112HqqbMDGUhRC1krR/hjoWpePGJWstsfyXxZT+h9dlPekWaKcWWC6M1FBYWVpvW2xQoikJ+fj4bN26Ua90pRtomdUnbpIZMK0bxBWPqLPeOPA177NlEIpFGrJWoS11/N4Zh1GuwqnQXCCGESCpbVdFatamzXO/SQ3pDmjEJIkIIIZIqoDvxjTuv1jLF7cXYp5/MnGnGJIgIIYRIqng8jjpwCL6/ng369qGLWm4bcu58lErDlbzKib2uSQxWFUII0bxVWOA+diytjjoBu6y0aoXdnFaUmEhvSDMnQUQIIURKCFkWIdUB2W2qBkTmtMZshBV2RXLJpRkhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjQQRIYQQQiSNBBEhhBBCJI0EESGEEEIkjZ7sCgghhBDNlaqquLFRzThx3SBs2di2nexqpRQJIkIIIcRe4FNstN9XEXj1WeJFm3H06kPWmAn4PelELSvZ1UsZEkSEEEKIBuZWwfr4Hcr/+1TiWGjTekKff0TOlEex8gqIx+NJrGHqkDEiQgghRANzRSNUvvx0zQLTpPyRqXji0cavVIqSICKEEEI0IFVVia/+Feq4/BJfvwYtHGzkWqUuCSJCCCFEA1O0Pxn5oMjH7zbymxBCCCEakGVZaAVdQNNqLTe69CDucjdyrVKXBBEhhBCigQUNJxkXXlvjuOJ0kXHJZIKqzBXZRn4TQgghRAOLWDbqoANpdf+zBN7+P8wtG3DsOxD3sOOpcLgxTTPZVUwZEkSEEC2WoijJroJoxkIWhNNzcJ51CbplEtd0SmIxkBBSjQQRIUSLo6sqPjOKvX49wRVLyGzXkajLTcCWYCIalm3bhLetF2LFkluZFJUSQWTOnDnMmTOHwsJCANq3b8/YsWMZOHAgUNWQr776Kh9//DF+v5/u3btz7rnn0qFDh8R9xGIxZsyYwYIFC4hGo/Tp04fzzjuPnJycpDwnIURqMjQV78a1lEy5BjsUSBx3H3o0mRMuoUwWvBSiUaXEYNXs7GxOP/10pk6dytSpU+nTpw/33HMP69atA+DNN9/k3Xff5ZxzzmHq1KlkZmYyZcoUQqFQ4j6effZZFi5cyKRJk7j99tsJh8PcfffdWLKMrhBiB75omOJbL68WQgBCn80h+tkHOIyU+H4mRIuREn9x+++/f7WfTzvtNObMmcOvv/5K+/btmT17NqNHj2bIkCEATJw4kfPPP5/58+czfPhwgsEgc+fO5dJLL6Vfv34AXHrppVx00UUsXryYAQMG1Pq4sViMWGx7V5miKLjd7sS/m5Jt9W1q9W4JpG1Sh67rxH74EmK1r2oZmPUCmQcfRUxmNCSd/N2kroZum5T7a7Msiy+//JJIJEKPHj3YsmULZWVl9O/fP3GOYRj07t2bZcuWMXz4cFatWoVpmokQAlW9LAUFBSxfvrzOIDJr1ixee+21xM+dO3dm2rRp5Obm7rXnt7fl5eUluwqiDtI2qaFs0+91llkV5Th1jfw2+Y1YI7Ez8neTuhqqbVImiKxdu5bJkycTi8VwuVxcffXVtG/fnmXLlgGQkZFR7fyMjAyKiooAKCsrQ9d1fD5fjXPKysrqfMzRo0dzwgknJH7elu4KCwub3GZEiqKQl5fHpk2bZIvpFCNtkzo0TcPTq2+d5Xr7joRNi8qNGxuxVqI28neTuupqG13X6/VFPmWCSNu2bbn33nsJBAJ89dVXPProo9x2222J8j92Ae3KC/PPzjEMA8Mw6nXbVGXbdpOte3MnbZN88XgctWN3tJzWmMVbapSnTbiEoO7AbmJfRJoz+btJXQ3VNikxWBWqklReXh5du3bl9NNPp1OnTsyePZvMzEyAGj0bFRUViV6SzMxM4vE4fr+/xjnbbi+EEACVuoPsOx/F2Xe/xDE1I4vMy2/B7NyjyfWGCtHUpUyPyB/Ztk0sFqN169ZkZmayePFiOnfuDFR9q1m6dCnjx48HoEuXLmiaxuLFiznooIMAKC0tZe3atYlzhBACwDRNSnUXnkm3kBYNodk2ccNJwHASkxAiRKNLiSDy3//+l4EDB5KTk0M4HGbBggUsWbKEyZMnoygKI0aMYNasWeTn55OXl8esWbNwOp0ccsghAHg8HoYNG8aMGTNIS0vD5/MxY8YMCgoKqg1gFUIIqPqiEwCCDg/5+fkUb9wol2OESJKUCCLl5eU88sgjlJaW4vF46NixI5MnT06EiFGjRhGNRpk+fTqBQIBu3boxefLkxFRbgAkTJqBpGg888EBiQbPrrrsOVU2Zq09CCCGE+APFllFANRQWFlZbX6QpUBSF/Px8Nm7cKAO7Uoy0TeqStkld0japq662MQyjXrNmpLtACCGEEEkjQUQIIYQQSSNBRAghhBBJI0FECCGEEEkjQUQIIYQQSSNBRAghhBBJI0FECCGEEEkjQUQIIYQQSSNBRAghhBBJI0FECCGEEEmzy3vNTJw4EUVRdulcRVF4+OGH610pIYQQQrQMuxxEevfuvctBRAghhBBiV+xWj4gQQgghREPa5TEizz77LL/88sverIsQQgghWphd7hH5+eefee+998jIyGDw4MEMHTqU3r17o6oy3lUIIYQQ9bPLQWTatGls2bKF//3vf3z11Vd8+OGH+Hw+9t9/f4YOHUrfvn3R9V2+OyGEEEKIXQ8iAK1bt2bkyJGMHDmSkpKSRCiZNm0aLpeLQYMGMXToUAYMGIDD4dhbdRZCCCFEM1HvLozs7GxGjBjBiBEjKCsrY+HChXz11Vc88MADGIbBgAEDuPLKKxuyrkIIIYRoZhrkWkpmZiZHH300Rx99NH6/PxFKhBBCCCF2psEHdfh8PoYNG8awYcMa+q6FEEII0czUK4j8+9//rrNMVVU8Hg/dunVj8ODBMoBVCCGEEHWqV0pYsmQJwWCQYDCIqqqkpaVRWVmJZVl4PB4A3n33Xdq2bcstt9xCZmZmQ9ZZCCGEEM1EvYLIVVddxX333cf555/P0KFDUVUVy7L48ssvefHFF7nyyisxTZP77ruPl156iYsuuqih6y2EEEKIZqBeq5E9//zznHjiiRx00EGJBc1UVeXggw/mhBNO4LnnnqNnz56MGjWKH374oSHrK4QQQohmpF5BZOXKlbRv377Wsg4dOrB69WoAOnXqRGVlZb0rJ4QQQojmrV5BxO12s2TJklrLfvrpJ9xuNwDRaDTxbyGEEEKIP6rXGJFDDjmEN998E9u2OfDAA8nIyKC8vJwvvviCt99+mxEjRgCwatUq2rVr16AVFkIIIUTzUa8gcvrpp1NaWsobb7zBG2+8Ua3s4IMP5rTTTgOgR48eDBgwYE/rKIQQQohmql5BRNd1Jk2axJgxY1i6dCl+vx+fz0fv3r2rjR3p169fg1VUCCGEEM3PHq021r59+zoHrQohhBBC/Jk9CiLhcJiioiKi0WiNsi5duuzJXQshWghd11FVlXg8jmVZya6OEKKR1SuIVFRU8Pjjj/Ptt9/Wec4rr7xS70oJIZo/h6rgjUWIfjsfq2gz7r77obQtoELVJZAI0YLUK4g8+eSTLFmyhBEjRtCuXTvZT0YIsVscqorrt18omnodmGbVwdeeQ+/Qiax/PECpomPbdnIrKYRoFPVKED/99BNnnnkmRx11VEPXRwjRAnhjYYruvn57CNkqvm41/v8+ievMiYRMCSJCtAT1WtDM6XSSm5vb0HURQrQAmqYRX74E4vFay0Off4gzFmnkWgkhkqVeQeTQQw/lyy+/bOi6CCFaAEVRsMpL6z4hHq8zpAghmp96XZoZN24cjz32GPfeey+DBg3C5/PVOGfIkCF7XDkhRPNjmibGPnWvMaS37YDlcDZijYQQyVSvILJlyxZWrFjBxo0b+eabb2o9R2bNCCFqY9s2VmYOjj77Ef2p5sy79POuJKA7pFdEiBai3rNmgsEgEyZMoH379jJrRgixWypRybz8H4Tff53Au69jhwIYHbuSdu7lRNt2JC4hRIgWo14J4tdff+XCCy/kkEMOaej6CCFaANu2KUXBefw4so8+CWwbU9PxqzrmH2bSCCGat3oFkYyMDLxeb0PXRQjRwkTicSKKDgpgU2M6rxCi+avXrJmjjz6aDz/8sKHrIoQQQogWpl49IoqisHbtWq677joGDhxY66yZE044YY8rJ4QQQojmrV5B5MUXXwSgsLCQ1atX13qOBBEhhBBC/Jl6BZFHHnmkoeshhBBCiBaoXkFkd5Z3t22b119/naOOOorMzMz6PJwQQgghmql6DVbdHbZt8+qrr1JSUrK3H0oIIYQQTcxeDyJCCCGEEHWRICKEEEKIpJEgIoQQQoikkSAihBBCiKSRICKEEEKIpJEgIoQQQoikkSAihBBCiKRp0CBiWRbffPMN99xzz/YHUFUeeeQRCgoKGvKhhBBCCNEM1Gtl1T/asGED8+bN47PPPqOsrAyHw1GtfHdWYhVCCCFEy1HvIBKJRPjyyy+ZN28ev/zyCwDt27dn9OjRHHrooQ1WQSGEEEI0X7sdRJYvX868efP44osvCIfDuFwuDjnkEObPn8+5555L796990Y9hRBCCNEM7XIQeeedd5g7dy7r168HoFevXhxxxBEceOCBxONx5s+fv9cqKYQQQojmaZeDyIwZMwAYNGgQEyZMIC8vL1FmmmbD10wIIYQQzd4uB5FOnTqxevVqvvvuO8rLyzniiCM4+OCD8Xg8e7N+QgghhGjGdjmITJs2jdWrV/Pxxx+zYMECpk+fzvPPP8+QIUMYMmTI3qyjEEIIIZqp3Rqs2qlTJ84991zOOussvvrqK+bOncvnn3/O559/DsBXX31F27ZtyczM3Bt1FUIIIUQzU6/pu4ZhcMghh3DIIYewZcsW5s6dy6effsr777/Phx9+yH777cdVV13V0HUVQgghRDOzxwuatW7dmnHjxnHqqaeyaNEi5s6dy7ffftsQdRNCCCFEM7fLS7z7/X7uu+++OkOGoiiYpollWTz++OMNVkEhhBBCNF+7HETmzp3LmjVrGDBgQJ3nDBgwgHXr1vHBBx80RN2EEEII0cztchBZsGABRx55JJqm1XmOpmkceeSRfPPNNw1SOSGEEEI0b7scRDZu3EjXrl3/9LzOnTuzcePGPaqUEEIIIVqGXR6saprmTntDttE0jXg8vkeVEkIIIZLN0HU8sQgKYKoqQVWXlcT3gl0OIllZWfz+++9/uqnd77//LuuICCGEaNIyVIgv+JCymTOwSoowuvYk7exLieYXELKTXbvmZZcvzfTu3Zs5c+bstLcjHo8zZ84c9t133wapnBBCCNHYfAqEnv83FdMfwCopAiC2chklN1+CsfJnDMNIcg2bl10OIscffzzr16/nvvvuo6SkpEZ5SUkJ9957Lxs2bOCEE05o0EoKIYQQjcUI+Ql9+n6tZRVP3Y8nGm7kGjVvu3xppmPHjpx77rk8/fTTXHLJJXTp0oXWrVsDsGXLFlatWoVt25x33nkUFBTstQoLIYQQe4umacR/+7XOcrNoM0o4BG5fI9aqedutlVWPOuooCgoKmDlzJkuWLOHXX6say+FwMGDAAE466SR69Oix25WYNWsWCxcuZP369TgcDnr06MEZZ5xB27ZtE+fYts2rr77Kxx9/jN/vp3v37px77rl06NAhcU4sFmPGjBksWLCAaDRKnz59OO+888jJydntOgkhhGh5bNtG8e48ZCjGHi9KLnaw27/NHj16cP3112NZFpWVlQCkpaWhqrt8laeGpUuXcswxx9C1a1dM0+Tll19mypQp3H///bhcLgDefPNN3n33XS6++GLy8/OZOXMmU6ZM4cEHH8TtdgPw7LPP8u233zJp0iTS0tJ4/vnnufvuu5k2bdoe1U8IIUTLYFkWWvtOKE4ndiRSo9zZ/wCihgtkwGqDqfens6qqZGRkkJGRsccf8pMnT+bwww+nQ4cOdOrUiYsvvpiioiJWrVoFVCXU2bNnM3r0aIYMGUJBQQETJ04kEokwf/58AILBIHPnzuWss86iX79+dO7cmUsvvZS1a9eyePHiPaqfEEKIlsOvO8m6bir8YckKLSeX9AuvJYiSpJo1TynZvxQMBgHw+aq6x7Zs2UJZWRn9+/dPnGMYBr1792bZsmUMHz6cVatWYZom/fr1S5yTnZ1NQUEBy5cvr3Vp+lgsRiwWS/ysKEqid0VRmtYLbVt9m1q9WwJpm9QlbZO6ktk2cdsm3KkHuY+8THjh51gb12H03R+9e28qDCeYZot+zTR026RcELFtm+eee45evXolBr2WlZUBkJGRUe3cjIwMioqKEufoup4ILzues+32fzRr1ixee+21xM+dO3dm2rRp5ObmNtCzaXx5eXnJroKog7RN6pK2SV3JbhtXQedqP7uTVI9U1FBtk3JB5Omnn2bt2rXcfvvtNcr+mL5s+88v0u3snNGjR1ebarzt/gsLC5vc6rCKopCXl8emTZt26fciGo+0TeqStkld0japq6620XW9Xl/kUyqIPPPMM3z77bfcdttt1Wa6bFuptaysjKysrMTxioqKRC9JZmYm8Xgcv99frVekoqKCnj171vp4hmHUuTBNU33h27bdZOve3EnbpC5pm9QlbZO6GqptUmIqiW3bPP3003z11Vf84x//SKxPsk3r1q3JzMysNug0Ho+zdOnSRMjo0qULmqZVO6e0tJS1a9fWa0qxEEIIIfa+lOgRefrpp5k/fz7XXnstbrc7MabD4/HgcDhQFIURI0Ywa9Ys8vPzycvLY9asWTidTg455JDEucOGDWPGjBmkpaXh8/mYMWMGBQUF1QawCiGEECJ1pEQQmTNnDgC33nprteMXX3wxhx9+OACjRo0iGo0yffp0AoEA3bp1Y/LkyYlZLgATJkxA0zQeeOCBxIJm1113nawhIoQQQqQoxZaLbzUUFhZWm9bbFCiKQn5+Phs3bpTrqSlG2iZ1SdukLmmb1FVX2xiGUa/BqtJVIIQQQoikkSAihBBCiKSRICKEEEKIpJEgIoQQQoikSYlZM0IIIUQqUFUVr22iBwPY0QiKL42Qw0XYtJJdtWZLgogQQggBaJpGRthP2T2Tia1eUXVQ1/GNPI304/9KhWSRvUIuzQghhBBAWjxKyc2Xbg8hAPE4/pkzMOd/hMOQ7+57gwQRIYQQLZ6qqtgb1mIWb6m13P/6c3gi4UauVcsgQUQIIUSLp2kasbUr6yy3Ksoh3rQWumwqJIgIIYRo8UzTRG/fqc5yxZcGdezWLvaMBBEhhBAtnmVZqO07oWZm11ruG3U6QcPVyLVqGSSICCGEEECl7iTnjkfQ89tvP6iqeI4bgzHseKLxePIq14zJEGAhhBCCqssz5d4M0m57GC3oxw6HUDKyiDhclFuy8d7eIkFECCGE2Mo0TSoUDbwZVf8BSAjZqySICCGaDUVRcKsKjq2zG8KGg3DcTHKthBA7I0FECNEs6JpKeshP5YtPUPnVZyi6jnvY8WSffAZlqgPLkmUxhUhFEkSEEM1CeiRE0dXnYocCANimSfC9mUS+/ZKsO/9NKVqSayiEqI3MmhFCNHkuXSPwxouJELIjc8tGYj9+i67L9y4hUpEEESFEk+eIRoh880Wd5ZH5HyFLUQmRmiSICCGaPFtVUTzeOssVXzq2Km93QqQi+csUQjR5YcOJd+S4Oss9I8YQick+IUKkIgkiQogmLxaLoQ8cinPgkBpl3pHjiLfKw7ZlLQghUpGM3hJCNAtltkL6xTfgK95M+POPUJwuXH8ZTjQtA7+tJLt6Qog6SBARQjQbFagordriOPU8LNumLBaTnhAhUpwEESFEs2LbNpFIJNnVEELsIhkjIoQQQoikkSAihBBCiKSRICKEEEKIpJEgIoQQQoikkSAihBBCiKSRICKEEEKIpJEgIoRIeYqi4NZU0mwTr6qgbt03xqXrVccU0DQNAIeuk4aJDwvDkK3uRGoxDAOfLa/PHck6IkKIlOZUFTz+cgKvP0/416VouXmkn/I39DZt8b82A//3/0NNS8d70ngc3XsT/OgN/J9+gKIbuI8bg3fQgZQrmixsJpJKVVUyrBiReXMIzH0XRVVxH3UivgMPp4yW/fqUICKESFm6ruP6fRVFt0wCywQgvvF3Iou/Ie2Uv2EH/cQ3rAMgOu1G3AcPw+jcnfjvawCIPTYNR899ybj6TsqkA1gkUYYVp+wflxLf+HviWOzpB9E/epuMm+6jDC2JtUsu+csUQqQsbyxC+cN3JkLIjipfex7PYcdUOxZaMBejYzcUpzNxLLpsCdaqZYlLN0I0NkPXiX45t1oI2Sa+ZiXxxd+g6y23X0CCiBAiZanhIOaWjbUXWiZm8RbU9IxqhyNLfsDo2qvasdCHb+FQZeM7kRzOeJTQvPfqLA999A5OM96INUotEkSEEKlL+ZPwoGrYllX9JpoGfziGpgMSRESyKKDW/XGraBr2n73WmzEJIkKIlGW6POgdOtVeaDjQMrKw/ZXVDjv3HUB0xS/VjnmOO5mIWfPyjhCNIWI48Bwzus5y93EnE1Fa7sdxy33mQoiUF9AMMi67GcXhrFGWefYl+D94o9ox74gxRH76DuKxxDHn/gdjt+uI9cdeEiEaSSwWQx80FKPbPjXKHPsOQOvRB7MFB+WWOzpGCJHyTNMkkJNHq4deIPThW0R/WYye1x7PiaeipqVjxqJY/grU9Ey8J45Db51H+MdvcfbdD3Qdz/F/RenUnXK75XZ7i9RQjkbGdVMxl/1IaM6boKq4jx2N2nUfyhUNWvD0XcVuyZOX61BYWEgsFvvzE1OIoijk5+ezcePGFj0fPRVJ2zQMp66jWXEsVSdimti2jaHrGJaJralETBvLstB1HYdtgqISsdnpN01pm9TVXNum6vVZ1TsXQWmSPSF1tY1hGOTm5u72/UmPiBCiSYjEt84qsLbPLojF48Sg2vTeeDxO1RlN7w1eNH/bX59iGxkjIoQQQoikkSAihBBCiKSRICKEEEKIpJEgIoQQQoikkSAihBBCiKSRICKEEEKIpJHpu0KIJsswDHRsbEUlEos1q/UmhGgpJIgIIZocTdNIj0eJLfyCyMLPUTMyyThuDLHMHAKyiqoQTYoEESFEk6IoChmxMCU3XoRZvCVxPPjxu6SdfgGeYSMIShgRosmQMSJCiCbFpar4//tktRCyTeV/n8AZDiahVkKI+pIgIoRoUpyxCKH5H9VZHvl6AYZhNGKNhBB7QoKIEKKJsWEnG4XZkTCKIpdmhGgqJIgIIZqUmG7g7Ld/neXO/Q9ucrtnC9GSSRARQjQpIVTSzp2E4nDWKHMNPRQzM1um8QrRhMisGSFEk2JZFn5fFq0eeA7/q88RWbQQ1ZeO96TT0foPplxmzAjRpEgQEUI0OTHLotTpxXXWRDLi52MrKiHdQTweT3bVhBC7SYKIEKJJsm2bkGkTUra+jUkIEXuJoWl44xGoLEdRNSxvGgHDKcG3gUgQEUIIIergVkD/+XtKHpuG7a8EQM1uRdZVtxPMLyBqWkmuYdMng1WFEEKIWqiqirO0kLJ7b0qEEACrpIjiWy7DFw4ksXbNhwQRIYQQohYubPwvT6+9MB4n9NHbOGXxvD0mQUQIIYSohR6PEV+7qs7y2Iqf0ay6F9cTu0aCiBBCCFELUzfQ23ass1zv2A1T1RqxRs2TBBEhhBCiFiEUfOPOrb1Q1XAfcxIRWcV3j0kQEUIIIWphWRbR3HwyLrkRxelKHFd96WTfdB9Bty+JtWs+ZPquEEIIUYegDY5BB5Hzr/2xy0tQVBU7PZOA7iS2k80Xxa6TICKEEELsRNS0iKo6ZLXefrCBQoiiKLhUBT0ew9J0QihYVstam0SCiBBCCJEETlXFEygn8NrzhH9dgpabR9pfzyaeX0CgBe3bKEFECCGEaGS6ruPesJqimy+FrVOA4xt/J7L4G9LOuhjXYccRtlpGGpHBqkIIIUQj88SjlD18ZyKE7KjyhSdwxyJJqFVySBARQgghGpkaCmBuWl97oWVirl2FqraMj+iW8SyFEEKIFKIoys5P0HRsWy7NCCGEEGIvMF1e9IIutRfqBlr7jhJEhBBCCLF3BDSdzMtuRnE6a5RlXHwdQd2RhFolR0rMmlm6dClvvfUWv/32G6WlpVx99dUMHjw4UW7bNq+++ioff/wxfr+f7t27c+6559KhQ4fEObFYjBkzZrBgwQKi0Sh9+vThvPPOIycnJxlPSQghhKiTaZr4s3Jp9eALhD5+m+jSRWht2uE98VTC6VlEWsiMGUiRHpFIJEKnTp0455xzai1/8803effddznnnHOYOnUqmZmZTJkyhVAolDjn2WefZeHChUyaNInbb7+dcDjM3Xff3eIWhhFCCNE0xCyLEt2Jefw4nFfdgXrWJZSlZxNqORkESJEgMnDgQMaNG8eQIUNqlNm2zezZsxk9ejRDhgyhoKCAiRMnEolEmD9/PgDBYJC5c+dy1lln0a9fPzp37syll17K2rVrWbx4cWM/HSGEEGKXReNxghaE4/EWMy5kRylxaWZntmzZQllZGf37908cMwyD3r17s2zZMoYPH86qVaswTZN+/folzsnOzqagoIDly5czYMCAWu87FosR22HnREVRcLvdiX83Jdvq29Tq3RJI2+w+h2Gg2RZxRa32N9rQpG1Sl7RN6mrotkn5IFJWVgZARkZGteMZGRkUFRUlztF1HZ/PV+OcbbevzaxZs3jttdcSP3fu3Jlp06aRm5vbMJVPgry8vGRXQdRB2ubPWdEI5uaNBN57ldjqFTh69SHj8OPQW+ej6Hvv7UraJnVJ26SuhmqblA8i2/wxee1K99WfnTN69GhOOOGEGo9RWFhIPB6vRy2TR1EU8vLy2LRpU4vs2ktl0ja7xlBV3Gt+peSOqxKrTYa//YLKV58jZ8qjVOa2xWzg3U6lbVKXtE3qqqttdF2v1xf5lA8imZmZQFWvR1ZWVuJ4RUVFopckMzOTeDyO3++v1itSUVFBz54967xvwzAwDKPWsqb6wrdtu8nWvbmTttk5TzxCyX031Vjy2o5GKL3/FjLueIRytL3y2NI2qUvaJnU1VNukxGDVnWndujWZmZnVBp3G43GWLl2aCBldunRB07Rq55SWlrJ27Vp69OjR6HUWQtRDWQl2wF9rkblpPWqw9jIhRNOWEj0i4XCYTZs2JX7esmULq1evxufz0apVK0aMGMGsWbPIz88nLy+PWbNm4XQ6OeSQQwDweDwMGzaMGTNmkJaWhs/nY8aMGRQUFFQbwCqESF32n1wOtU2Zii9Ec5QSQWTlypXcdtttiZ+ff/55AA477DAmTpzIqFGjiEajTJ8+nUAgQLdu3Zg8eXJihgvAhAkT0DSNBx54ILGg2XXXXddiNg0SoqlTc3JBNyBec5aMmpYBvrQk1EoIsbcptlx8q6GwsHCvThncGxRFIT8/n40bN8r11BQjbbNrnJqCPv8jKp55qEZZ5jVTCPceRKyBB5FL26SuZLSNYRi4o2EUIKYbhCwZn1KbutrGMIzmOVhVCNEyREwb/aBhZHfsgv+l6cQ3rcco6Ipv/N+J5LRp8BAixI4yFZvY/DlUvP0KVmU5zv4HkHX636n0pBOTFbr3KgkiQoiUEbAVtA7dcF99J6oZx9QNKlBkqwaxV6UrFv5H7yLy/cLEsfCXnxD+egGt7nua8rRseQ3uRTKAQgiRUkzTxG9DhaoTsGz5ABB7laIoqMWF1UJIQjxG5dMP4UZeg3uTBBEhhBBNimEYZNgmmbEwmXYcxx6sumsYBpFvFtRZHvnxW4xaBlCLhiOXZoQQQjQZaYqN/c3nlP33SaySIhSvD9/I0/AceQJl9u7vfWLbNorHW/cJugGy381eJT0iQgghmgSnrmMt+JjyR+7CKqnaa8wO+Kl86SmCLzyOtx55IRaL4TzgkDrLPUccR1h31rfKYhdIEBFCCNEkuKMhKl+eXmtZ6JP3cERD9brfsMtL2oRLahzX8trhPeVswg28x5GoTi7NCCGEaBqCAexQsM5ic/NGlA5dd3vtj5ANnkOG02rgYEIfvoVVWozzoGHoPftQphqwiwOmFUWRdUfqQYKIEEKIJkExHDstV72+egeBIAohXxaOU89HBUKmWbUL+y6EEK8CjkgIq3gLqi8NKy0Dv+Zo8N2imysJIkKIZk++qTYPMZcHxz79if68qEaZmpUDmTl7dP+2bROJRHbrNpmqTeCZf1Gx4OPEMa1NW7Jv/iflnnQJI7tAxogIIZolRVFIV2yyAuWkrV5GVmUp6Viy/1QTFlRUMibdhJabV+244vGSffM/qdR33mPS0Fy6RvjNlwjvEEIAzM0bKP7HpaSZ0UatT1MlPSJCiGZHURSyrDhl99xAbMUvieNGx65k3ngPZbpTFkprgizLosxwk3nXY9jr1xJbuQytXQf0zj2oNFz16n1QVRWPbWJEQlU7QLs9BAzXLm0p4IqGKXl/Vu11LSnC3rQeJb+j9Mb9CQkiQohmx6vYVPzrjmohBCC2ZiVl99yI94Z7qETWhmiKLMuiDA21Q1e0Tj2IWFZVAKlHCNE0jYxQJeUP3UH58iVA1SWejPOvQu/Rh9Cf5YdoBDta96Ucc9MG1Had5fLMn5A+SiFEs2OEgkR++q7WstjKZehBfyPXSDQ0y7KIxWJ79CGfHo9QfMOFRLeGEACrtJjSe27EsWkdmqbt/A4crp0uhqZ36CQ9b7tAgogQotmxw3VP8QSwg4FGqolIVbquE1vyA1ZFWa3l/ucexW3t/PJMyOHCN+as2u+/XUfsnNZyWWYXSBARQjQ7ijcNdjIoVUnLaMTaiFSk6zqxH7+tszy6ahn6n/RmROJx9MOOxTd2AuwwtdjRZyBZ/7ifSlVGP+wK+S0JIZqdiNOF+/BjCc2dXaPMNfQwoi43sqFqy2ZZFnr7jnWWa7l5u/QSqbAVXCNOodVRI7GDfhSnk5jTQykKtlyW2SXSIyKEaHaCpo3n9AvwHDMatu3Mqmq4h43Ad96VBOTzocWLRqO4DjwC6hgH4vvr2YQcrl26r7BpUqoZlKVlUerw4LeRSzK7QXpEhBB7RNd1XPEoqm0R1w1CVmq8CZfZCu5x59Jq9HjsUBDF7SbicFNqSgoRVfwON9k3/ZPSaTdgh7fvU+MZMQa1/2BisVgSa9dySBARQtRbumJj//w9gVefxSwpwtGjD5mnn08wPYuIlfwwEopbhDQH+LZev5cQInYQtSzsjt3I+deLWJvWY4WC6B06EXG4qUj+y7fFkCAihKgXjwLRN14k8M7/JY6Fv/6c8LcLyL79Ycx2nav26hAihcVMi1JFR2nbCUVRqqbbSghpVDJGRAhRL85IqFoISbAsKh67B0989/bsECKZbNuWNT+SRIKIEGK3aZpGbMXPdZbH169BDYcbsUZCiKZKgogQYrfZto3i2PkGY4omby9CiD8n7xRCiN1mWRZ65551Tn109O5PzOFu5FoJIZoiCSJCiHoJGg4yL5lc47jqSyf9ousIKvL2IoT4czJrRghRLxHLRu27P63+9SKhD97A3LwBx4DBOA44hArDhSU7jgohdoEEESFEvYVsCLnTcJ5yLrptEwUCsVi9tmQXQrRMEkSEEHssEo0muwpCiCZKLuIKIYQQImkkiAghhBAiaSSICCGEECJpZIyIEEII8SecTieKbWPatuzK28AkiAghhBB1cKoq3kiA0IezMNeuwtF7AL4hh1JpuIjL7LAGIUFECCGEqIWhqbhW/0rhlKvBqgodof99ivLSdHLueoyKjBxMCSN7TMaICCGEELXwxqKU3js5EUK2sUMByh64Fa8ll2gaggQRIYQQ4g8URcEu3owdCtZaHl+7Cq2OMrF7JIgIIYQQtbAjkZ2Xy2WZBiFjRIQQQuwVmqZh6Do2EI1GsW072VXaZbZto+W3A1UFy6pRrmZmg8eXhJo1P9IjIoQQokEpikImFt5VP2M9MQ2e/RcZRRtJU5pOEAEIGy58fz271rKMC64mYDgbt0LNlPSICCGEaFCZmJRPu5HYiqWJY8G57+I5ehRpp5xDpa0ksXa7LmTZpB19Elmde+B/5WnimzdgdOlB2hkXEmmVTyweT3YVmwUJIkIIIRqMoetEPp1dLYRsE5zzJq4jRqDktm0yl2kqLdB69MV7432otklc1alQVKxaLteI+pFLM0IIIRqMOxYh+P6sOstD78/EaRiNWKM9Z5omlSiUKzoBGwkhDUyCiBBCiAZkY0fCdZeGgihNpDdENA4JIkIIIRpM1HDiGnJYneWuI44j2gA9Crquo+syuqA5kCAihBAtkKIo6LqOqjbsx0A4buI96XQUX1qNMqNTN7Ru++zRsuhuBbLjYVzfzsf5v3lkhf14m9hsHFGdxEkhhGhBFEUhzTZRS4qIrV6B1qo1avtOVOrOBts3pdzhptV9/yEw6wXCX8xDcTjwHHMSzmHHU4YG9bw041Vs7PkfUfjcI9Xuw3PcyaSPmUBFE5mNI6qTICKEaLE0TcOsLMdnx4npTsLNfDqmoihk2SZlU68jtmpZ4rjqSyf79n9RkZnbIGHENE1KNAeuceeTOeYsbKrW5AjG40D9QoiiKBilRRQ/+3CNsuB7M3EOGIzWva9sQtcEyaUZIUSLlK7YeJYvpuiWSfhvuQxefYbseARd05Jdtb3GrUDlfx6qFkIALH8FJbdOIi2+8yXNd1c4blKu6FQoOtE9DHlOXSP47v/VWR54fQZuq3kHyeZKekSEEC2OT7EJv/wUoY/eSRyLb1hH8ON3yLn3acq9Gc1yiqYzFqbyf5/WWmZVlGNv2YiSV5CSa3yolkW8uLDOcrOsBMWyQGm+QbK5kh4RIUSLY/grqoWQbexImMpnHsJtN78QAls3cdtJwLJKixt88GpDiSkqjkFD6yx39BlIXG9a65OIKqn5ihNCiL3EMAwi331ZS4ED1ZdOZPE3GPFo41esMbg8qL70Oov19p3qPcbCoetk2CaZtolPocEDTSwex3XQsNrrbzjwjTmLsJV6PTmpRtM00hWLrEiArGiQNMVGS/LlSLk0I4RoeXZYf0Jv15H0U/6G4nJh+SvQclqj6np9x1SmtKDDhe+086h46v4aZY6++xFPy9jt571tg7vIh29Q9s4rWJUVOPYdSPrZlxDKyiXSgOGg3HCRM+1JKp56gMgPX1XVu3tv0i+8hkqXDztJl9M8moozEsL2V6A4XZhuL35VT7nLew5NxVO4nvJH7iK+bjUARufuZFwyGX9WLrEkDfRV7FS8GJhkhYWFxGKxZFdjtyiKQn5+Phs3bkzJ67stmbRN6skKByi69DT0Dp3JPHcSpQ/fibnD+APXkEPxnX8Vpc1wOmiaAvY386n875NYFWWgG3iOPB7PX/9GGepuv0bTFYvAg3cQ/enb6gWqSs7Ux6nMbVevXpa6/m5UVcWNhSMaxbYtTIeLQCPt/eLUddzRMMSjoDsIOt144jFCrz5D8MO3Epe9jE7dyLzuLsqd3pSaxZMdCVJ4+Rnwh4HDitNJqwdfoETftd2E62obwzDIzc3d7XpJj4gQosWJuL34xk7A2bUnJfffglVRXq08/NVnaG3a4hp9FjZ21YdPJAIOB2Gni3A8tb7p7o5KG4yhh5M1aChEI2A4CDtclMZNdrc7RFEU1JKimiEEwLKoeOp+vNdPo5KGC3SWZREAArqj6oANNMKYngwVYh/OouSN/2IHA6i+dDIuuYHoymUEP3ij2rmx1SsouWUSmXf+G7/uRFVV4vF4UntIXLpK4JUXa4QQqBo7FPrwTZwjTycSa/yZRxJEhBAtTtCGtBFj4bflNUJI4pwPZpF7/F8JfvgWJW+/UrV/iuHAe/QoMkefQdku9JY4DR13JFT1ge9wEnK4iaTAWiWxWJwyRQfn1o+AeP2+teu6TvSHhTWPd+iE7/hT0LJy0KMRLKeboE3K9QiqqopD11Fsm6hl1dl74VEVwq/+h+D7MxPHLH8Fqm3jf/uVWm9jbtmIsmUjri2bsDasxbNPf/R2HYgZTvyojd5TYsTj+H9ZXGd5dMkPuEacQsNO4N41EkSEEC1SSNFwFm6qs9yORLCLt+B/7bntB2NRAu++ilVZhnfCpQR28gU3U7WJvPMKJe/8H3Y4hOL24B15GhnDR1GeWp/H9WbbNkpa9cGj7oOG4T54GBUvPkl8w1oAnP0PIOuCayh37fqlim2DXRVF2SsBJk2xUTeuJfT+TOxIGM+wEWjd9qFc0Ws8njMapmjOG7VVEjscqvMxYst+IvDRO4nfg57fnqzLbiLd7SOU3Qo9HAZsIg4nsb3cE2FqOlqrNsR/X1NrudY6H0vT6h1K94TMmhFCtEimaWJ07FZnuZqWgVlSiOL2ori91cpCn32IYyc7zHpUhdDLz+B/7bnEB5UdCuJ/5Wkis2bg1ur31utVICsWIrN4I1mhStIVq1Gn27p0jUyqZsa4NQ3TNHEMGAxKVe+Qmp6J98jjKfnnPxIfvgCRRV9TfOOFZJp/PhvJUFWyrBie5T9S+cZ/SS/cQAYWitJwl3fSFIvIi49TctPFhOZ/RPjr+ZRNu5GKu64lGxMfFk5j+/d0u6Ks9mnPloXictf5OFrrfMzyksTP8Y2/45/9Omo8ivH5HIJ3XkXgjisx5r5DFuZebcuwDd4xZ9VZ7hk5jrCZnEtH0iMihGiRbNvGbtUGvX3HWr8lpp3yN7Q27ci68GpAQU1Lxz/nTcJbFwSzykogt22t9+2Mhin6+O1aywLvzaTViacQ0nZtYOA2WYqN/6l/UvHVZ4ljRpeeZF53F2W6c6+OP9A0jYxYmMB/p1Py6QfYlon7oGFknXY+YaeHzMtuouyhO/AedSKVb71c64e2VVZCZOF8cg46nFJLqbW+uqbh3fI7xbdeXq2nwejSg8wb7qFM0fa4d0RVVdTNGwjNe69GWey3Xwl9+DZmyRbU9Cyyjx5FueFCcbpqva/g5x/iPXoU/rderlGmtc4Hy8IO+Ksd9xw0jNJHphJb/WviWOXz/yY4502y7niEUvbOVFrLsjDzO5A2/gIqX3pqexupGunnXU40Oxc7SdOfpUdECNFi+TWD3Nsfxtl/cOKY4nKTfsaFOHr2ofCGCyl54DZKHriVoqnX4xo4lKyJN4Cuo3p9dd6vXVle98Jhlontr9ytero1lcB/nyC8QwgBiK1aRumUq/FZcRRFwa2pZNhx0m0TQ2+475np8QjFN1xA8MO3sKMRiMcJfTaH4mvOxRWLEOt7AK0eeRnX4EOILl9S5/1El3xP4PUZZJi1z0pMVywCs15EcXv+8DyXE3jhMVzqnveKOHSd0A5jPf4oOG82jk7d8b/+PEVXTiAjVInpS0Nv26HGuaEvP8HZZyCeY0+GHXozjE7dyLrkRsqf/3e1842OXYkXbaoWQgC0vHboee2JL/m+Qdvtj/y2gn3E8bR69BUyr7mTzOum0urfr2AOPpxgEsdfS4+IEKLFsiwLPb8trstuJi0cxI5FUdwelFiUwktOq35yLErZY9NodcfD5N75GGZaBpi1f4Os6xt0otxRe7miKLgVcJhV4wVChpNIPI4rEqLokw9qvU183W9oleVkOZ34X36a0u++RPV48Z5wKr5DjqTM3v0puTsyDIPo53OxSopqlFn+CkIfvIF94mmUOj1kuj1oWa2Ih9bWck+gtc7DaN+J2NLvMQYcmFgmwaUoeKJBIp99gpbVisyzLwGg9Kn7E6EtNP8jWp12PiHNUe/nUsWuClPb6tSmLYqmE9+8HkwTOxZNrDNjB/xUPv9vnBdeR9aN91B88yVYpcXbb5vTGrV9Z/TufWg16jTsygoUlwvNcFB444XVzoWqsTLhrz6v9tiZ50zCLCsm9tuvWCXFpEeDVOguYnuphytkQ0h3ovbqD1T9DdR3N+SGIkFECNHiBWywHW5wuHHrKvFnHqr73Nmvo+bm4TnhFFBqfws1PT70gi7E166qUWZ06UHc460xU9ZQVXyBcvwzHqPy2y9RnC48w0eSdeKpEAmDVfcgQrt4M6VP3IdZtKXq8QN+Kv7zLxxffUbWFbdSUs/Ob1VVSTOjlH85r85zwgs/w3fcyUTQCGgG3rFnUf6vKTVPVBRcA4ZQ+fbLpP/1bCwFYlSFEOOnryl86I7EB2Lg/ZkYHbuSc/UdFE25umrKqWlW/b8eQUTTNDxmHC0aQrUNzCNGoGg63uEjia9fgx2NYHTuTuTH77BCASJLftj+/L6ej/ecSZR70smaNh1r/RrM9WvQC7qg5LfH73Bj2zYB3QlZVWtopCk2RpeeRL79onpFbDvRc6KmZ5I18QZKH7q92ho2FS89Rc4tDxLIL9hrYQRIqcXW5NKMEELsQIvHMTetr7M8vmUjqtNFaM6b1QY07siv6mRdPxUtt031+27TluxrphAxqo8PURSFtFAlxVefQ/jr+VWXb0IBAm+9RNltk9Bcrmqrwf6R4vJg/uHbN0B06Q9Ya1fi03b/koaqqmTGwoQ/ehvVU/dlKNXjw1KqPkpisRhavwNwH3VC9ZM0jcy/X0Xwk/dQvWlY/krcZUVkmVG8sRBlO4SQbWJrVhJaMBf3gUdUPU5WDio2abu51olDU0kr3kTgrqspnjiOwr+PQctuhaN3f4puv5Kyp+6n/LlHKbr1csySInwjT8MOB/9wLwqmaVKmaPgLuhE5eDix9p1RolGMBR9ifDKbrEA5PqWqbn5U0i+5gdypj5N58fU49h2A3r4jjj6D8BwxAgDf8WOpeHl6tRBS9cSjlNx5Dd5dGNjbXEiPiBBC7CCuGxg9ehNbs7LWcqNzd+Ibfye2ahm+o08iUstur5ZlUelJp9WUR4mt+IX4pvXobdqCaVI87QYyr7wd05eZmMrqVhX8L02vdskgUZ/f12DFYmRfdTvx1SsJ//A/osu2j8MwuvUivu63qh6DWoS/+gxv+04ENMduXaLx2SZld1+PWVFG5oSJhL9ZUOt53pHj0GwbVYGgqlNumqSNv5C0E04lsuR7FFVFz2+P//1ZuIccRmzFz5TcM7mqF6LnvrgHH1rnpYHgZx+Qc/3dhD7/kPSxEyh7dCpabj5pZ15E5S58j1YUBV/IT+ENF0K86jKQ4nBiByopr2WZ++C82Tj26YvnL8MJL5wPgGvoYUQNB2ydUWJZFmkqaCVbsLZswly9guCnH1D57MO4Dz2GzLMvQQ368b80nfAPX6F6fPhOPBX30MOIVZZjpGfi6LsfRvfeVLw0vdZ626EA9uYNKPkdU27tlb1BgogQIuXpuo4nHkE1TVBVgg430QbchkFRFHS9av2IcDxO9qjTCc6dXfPDXdfxHHoMRbdfidGhE7ZS94ehOxah+M5rMUsK0TKyMcuKEzMoSu+6lvQ7HqVi64epEY/i/+5/NevlSyP70slEv/+K4GdzwLZwH3g4aWMmUPrwnVVLiV86mcIbL6z7uTldRH74CseBRxKJ7PpyVVrQT2z1CgDMijI8hx1D8NPq41RcQw+FWJSi80/C6NKT9L9fRSg3n8q4RY7HS/jr+cTXr8Us3ITnyOOJrV+N3qkruVMfx6osR/WmEfqi7ss+diSClp1L9g1Tif36C9GfF8PPi3EfPQo1v+BPLy9k6SqBF59LhBAA56ChO33MwLuvkXnpjUDVFO60My+idOtsksTsoRefJfTZHGzLwj3kUFrddB+lT9xL5If/kVZyKltuuBBiVT0aZsBP+TMPEfp6Pr5J/6BM1cm67KaqAc07Yfkr99oaKqlGgogQIqV5FRv15++peP7fmJvWo6Zl4B11Op4jjqXM2vNZFGZZCeklm4n++C1qWgZGn0FE3Gnk3PYvyv41BXPLRqBqZkPmhIlUznoB4jE8J5xCxHBALIamabgsE9UyiRkOwnETLehPjBGJ/2GWjLlpPVqgErwZANgoqB4vZihQ7bzsiTdQ8dL0RCAAiK1cht6hM7n/fIawauDXddxDjyDwTu0rfLoGDiH45ScoBx1V5+/A43TgjIbBtok4XASjMexIhIyzL0FvWwCqgt6uI+7DjyG2ehV6q9boHTpj+Suo+O+TVfVatYziGy+k1X3PEMtoRRgdOxLGLNxUtTvuqNMwN28kOPddwt9+ifuAQ7AqynHsOwBmvVBrvRzdexP67kucPfpQssPCcqH3Xsfx96sJR+sOIl7Fxvx5EdFffqx2XPWlYZbWHHi7jVlajOJwkXb6BbgOP5oKw429dTXc9HiE4uv/Xm0QamjBx0QWf0P2NXcQWfQN5S88ngghO4r++C1K4SYys7KpeOkpfMeNIWfyvWCZRH76nsDH72AHt7e/XtAlpcZx7E0SRIQQKcuh6yjffE7Zo1MTx6zKcipfeAz3ulV4z5xI4A9fGFV1+wZouqbhi0dRAhVV3yx96QQMJ7Gtq0dmKjbF024ksvibHe+ArMtvQdt3ADk33I0dj6NoGmZZCeXPPkJ87SocfQah9z+AoGmSpoCy6mcCr8/AKi/B0Wc/skedhv1n+6tEI7hyXMRiMUyHi5zb/wXRSNWHd0kR4W8WEN+ysVoI2Sa+7jdC3/2PyJAjiMfi5I4aR+S7L6stIgbgG3Ua4cXf4DzwCPzRmh+OuqaRaUaIfvctFe+9hlVZgXPAYHJHjsNyOKj4/ENiK5eBouA58gTSx5yFXVFB5JcfsYIBnPsOwDNiLNFffkRNy8AOh6h87lHcl91MWNHwjj6D6M/XknXRtZQ/9yiRb79MPHb4f5/i6N2fzPOvxNFj35rTflWVtLETKJt+P1g26WdcCLZFfPMGrFAA5U86ChzhEJGVy1CzcmDj79V+d679Dyb8de2Xmhy9+qJ4fcSHj6IkHk/szWLoOpFP59SYCQNVr8nwt1/iPepEFLcbo1N3Qp9/iFm8ZfvTScsgvmYFRDrgO2ok/nf+j9CXn4Bt49rvIHKuv5vyZx4itnoF7sOPRVWqBr1WNsONF/9IgogQImV5omFKn3u01rLQpx/gPeVvBBwetG2Bo6wYq6wELb89qteHUllKfPMGcDiI/fwjgbnvkH7elWjd9sHSDCJzZlYPIVQtRBV491XSc1pTdPNEHN33wTfir+gdu5B9xa2YZhw1Mwc7HCQzHCL4yXsE3nk1cfv4hnWEPnmP3LufQGvTFnPzhpqVVzUURcHx3Rd4OnZFc7upeO05gp9/CPE4Wk4u2VffUWMdimrP/6O3yeq7H2ZZCabXR87N9xFZ/A2RRd+gery4hvyF6C8/Ef7+K7wnnoqqbp/Gq6oqLgXckSAVLzxOaP7H1eofnPceOdffjbmlagl8rVUe3iOPZ8t151Xbm0dxusi94xFa/+tFiEZQNB07HsNSVQJxE7r0xHfqOWCa1ULINtGli4gu/YHsa6bgf/MlAh+/ix0K4OjZh7SxZxH8+F3Mws1YZcUovjTCX8xD79CZ9JPPJGhXrSrrsS10y8RSIGRUBTtVVbE2rkNxOEkbfQbFSxdtf8xlS0gfdx5qekbNfYZUjfTTzyey9Ac8+wxkx+jmMOOEv/wExeHE6NoLgNjKX7CjERz79MM99DD8783Etm2cnbvjvvp2Yr/9iv+918k47e/YsQhWOITqchLftJ7Iom8Sl/7CX88nsuR7Wk2+l8jyJehZOWy55DScA4aQOfEGypp5x4hit4QLULupsLAwMb+9qZCt5lOXtE3dHLqOJxaGQCWKbhB3ewkoWqJHIzNYTvGkM7efv+8AfEePqlpWW1FRsnLw57YlzV9GyW2XYxZuRnF7yLluKpVv/JfID18lbuscMIS0UeMovmcyuXf+GzMtk9IrJ2BtvVbv6Lkv6aecQ3zzBsyiLTj26YtR0JXY6hXEN/2Omp5VtdFZWjpaTmuKbrucVjfdR9Gtl9f+3Hr2IX3cuRTddkWNMu+xo/EMO57wj9/gGXwoJQ/dTmzFL4lyo1svfMeeTPiHhYTmf1T9xrqBoqnoHbrgPeoEyh6/F9/oM/COGINixrEiYRTDQfjH71AMB64+A7H8FSguN2bhZrAttNy8qrknoSBbrj6n1vo7Bw6p6oWIRVHTsyiZdkOtA3jVzGxyrplC0e1XYEciGN16kXnu5dA6jwgaLkOn/N6biPxYyw69gKPHvqSfeRHBL+fh2qc/isNBbM0qAh+8kehRyJp4A1Y4iJ7XHmyb2G/LcR0xAjsWpfK5R7ACfhy9+uLouz+070JYVfFVlBD+ZkFVb1bRZvzvvlZV3/RM3Ecci++4MZQ/+yjh/31S9Wvt0JmM08/HikSwzCh2OEL04OGJzwKPpqL98D/0jCwiP31X9TvqM4jYujUYBZ0of/lpMs64kPDCz4ksXYSakUna6PGo3jSKbrsCq6Jse/t27ErGOZdRfPeN2DtcjvONGof7L0djrl+L4nKBpmOrKuGu+xCJps5nUl3vaYZhkJubu9v3Jz0iQoik8CqgfDufkmcfTixaZXTpQeb1d1ft7B6PoTvdOPvtT2TxN6SPOxfF46XsPw9XLa+u63gOO4aMU8+l6JZJiQ+ttJPGUznrRSKLqu8KG/nhK1DAN3IcFa88Q+Y5kxIhxOjak7SxEyi+96ZqS4unnXY+zl59iK1ZgR23UDQVZ9/9wLJofd9/IB4j57q7CH4xj9AXc6sNbo0u+wklvwPao7Owln6PNutZLH85aWMm4Nz/YCwzjl1ZgVm0mfiGdfhOOAX3Yceier1El/xA9LdfcQ8+BM9fhlP+wuN4hp+Ie8AQrPJSrEgYrVUb7HicvCdfJ/Ljt1T852HU9Azcg/9CdN1vuPruR6y4ELO8FNWXRtHtV2FuXFdVOV2n1W0PEV2yiFppGlp6JprHix1zgBmvusRRSxCxykqwgn7srZu2xVb8QuGNF5E79QmclokZDFQtElYHOxrBLNqMu/9giqdeV6Nc79gFR7/9sAo3Y5aXoue1Q8ttg3/mDBSnm4wzLiL87RdEfvoOq7wUzxEj8GVkoThdVLw8HeJxfCefQZuHZmBFItgBP8HPPqDypel4jzye9NPOw9yyCXPLBspnPI4dj5Ez+R5KH5pC5uBDscNB7HAIPSsH/8+LKfvwrUTdKme9SPZ1d1H5zitkjP87xVOvrxYsvEeMoOTB26uFEKiamux/62V8x55E5awXE8fD33yJ57BjKXv24apLQIqCa7+DyPj7VcQ0R7MdMyJBRAjR6DRNQ1v1M6WP3LX9oKqRdtJ4wnPewD/7dexgADU9k7S/TsB38pnEfRlUhqIYf7uC2AuPYhVuIvTFPNxDDqt2Ld7Rc18qXnqq1seNfP8VaSPHEfr8Q+KagbFPP2I/LyZ97ARKH5pSLYQ49z8Yre/+bMnI4+sDz+Sn4gg90hQOaqWRtXEl5tefE3jnVRSnE89hx9Lqxnsovmcy9g6b4UUs+K5CI5A3gKF3P0+6S6ciauHQFFzz38c1aCjxwi2k3fkEQW8WAUVlzYZiCtN60u3QfmRuWIEx+zVybrqP2PKf2HLd+dihrWtc6AYZZ10ENgR+X0d00KG4WufhcWpomTmYJUU48jtgKRAvLyXryltRsAktnI/n8ONQNAXsWj7YDAc510wh/M0CNk86EzsaQU3LwDfyVFz7HUj50zUXe7PDYRRdx45uDWK2TfkzD+Ea/BfCX3+O64BDagwa3cZ1wCEYXXpgx2Jk/v0qKl55Bqu8dGvZX8g8bxLF995MbMXP29um336knXY+qtNF4Y0XYfkrEmXBj94h8+Lr0Nt3SozvcPbsS2z9WoJzZ1ebhhz89AOcA4bgOewYyrZO51WzcsCyyL7yFkqmXkds1TK03Dakn/53AjuEkG0UTcfRpReVr8+oFkIUpxPFMGpdkRYg/N3/8I0YWy2IqL40Qv/7dPs4FNuuGiu0YS0Ztz1MWTNd+kuCiBBil/3ZdMLayl0uF7phYJkmwWDVh6jPNql4/t8YXXoQO+pkyG1LRoaX8EdvgScN46q7iMXiaOtXEVMdFHbsx//9sJ4t/jj7te/NAdNeJK1sM57yQso674tr+myUJV9j//gNFbkFxO5+HsPQUX9ZhLFwLigK6gGHYvYdSnlGJvYND7Iq7qDLFbcRXPkrofad0HoPgOVL0E46k/iAg4m3as1vwRiXvbaIwNYP2HeBf+sqDw/fh15t8gi88yp2JEJgzpvECzeRedG1BH7+CfWQo6FjdwKagx7uOGXBGLMWb8SybYb1aE22x6B00JG4DBVfd4XNpQHilsEVb/xIRXj7dvB98ltx11XTsNYvp+SB26rvXxOPUf7Mv8i6/RH+mzmYub8Wk7HZ4vQeHgYaXnwVFZglP1D+9IOJcKR4vGRM+gfFnmwimoPWBw1HC8eIDDqUuA3eUAVGTisKbQ0OPBaHBfFP38OqLKfixSdJO+VvuIYcWn3PG1VF7dAFxp6LNvfNxGJw0eVLSBtzFhUvPE7G3y5Da9sBtVN3okOHAwqOb+ZhLf8R1+BDCH46B2e//XEecAhZuW2qnqemo2bl1AghAJHF3+I5YgQVH79TLYRsU/b4fbR58HmgauaNWVaMoqi1roUS+eEr3EP+ktj80D30MALz3sPIa8e25W/dBx9JcN7sOl/3jh77UrnDrB4Axe1N9LjVyrYTs3G28Y0YS8UrT9c4Nb5hHfbGdajtOjfLXhEJIkI0I05Dx7l1Uayw4ST6hze6uuy4xbphGEBVr4XT6SQejxO2VCpjNmWhGG3S3Dh1Bcu20bFJdxvErar7CEbjaKqCU1OImhA1Lb7fUM5XqzczqF0mB3Zpha5AZUUl9t+vZ4neineWl1D8a5R7RhQQOeMKNldGePHrdawrC9Etdz/OHtKJ1m6Dc4d2JBizKA/H8EcsXLntiLctQAf8ER27/2F4Bw/D0FT0jBxsG5ztCzCOHomiQNS0sUwbp67iTPOSrUDMTMfXOh8bm/h1dxO3bAIRE8u2cRk6GS749yn9yUtz43KoxE2b0lAMBSjV8lFf/AxNVXDrChHTJqYoxPc/auvS3AqKZZHjMXDqCsfvm4dTV/liVRH75KfRMcuLadlETJvMDB+aovDQmP5g2yzeWI5tQZ+2GVRGLRzzP65zEz3/K0/T+6/Xstjn4OdNldxW6OfIHrlMGtST+IUjq51rBwOU3X09xrRnmboozOSje/CvrEOYP7eYgzrncPKAfZj+xW8s3VSJoSkc220kZx81Bm3qJKzyUvxvv0L2pJurBRHHiFP47zqT+XY/zrrwILpvXIL9xFQUlxs7HqtaQExRYMrTvLR4E+/9UoxtwzH9T2fc2fkoRWuI/PAV3mNGEXjzv/h3GPibe8/0GiFkG9WXTvSn72t/QVsm0d+Wo7UtwDloKFY4ROSbL2o/Fwh+Ngf30MMJfPAGvmNHs/mqv6Hn5uE74RTKVt2P6nJj1bFRYWzVMhy9B9SsQkU5Wqs2NW+wleL2VOuRcg09DK11HvEN62p/nF9+RCvoJkFE7B7d4SCCAxvQsNDi4Wb5ImrqNMNJWDGImTYOTSFNtzBNC2yLaC1THhuCruuoqkosFttpD4PD4cBWVJStdTEMA0XZuoW6qhGxNVTFxufQIRoj+Psa7M2/EyvoRqlLx3a4yfA4cegKbl2jImJiWhZuQ8Vl6JQEo4RjFooCXkPDoasoCkTiFssLA6woLKNLjpeebXykuxV8ToNwzERBI82pYwFl4ThxC3R163oYKIRMm3DMxKVr9GubzgEds7BtG9uGsGljOz2oXfZhoKYwpEcesbhNJG6RpilkOA3uHrkvmqoQiZvEq9Yww1BVMtwOfE6NecsLURUfW0wbr1Mjy+PAoaqJvTmcmoqqKDhUlYhpgWJj2wqqUhWOLMvGsmFFoZ/PVxZxbO880l1G1TYgNph2Vf3TXTo+p46mQkUojq4pZLl0wnGLmGWjKgoeQyWwdT2LomCEyojJ/333Owt+K8bQFI7vncf4AwpwG+DSFA7s0orPVhTx3FfryM9wc9w+bWjlNQibFoX+CFv8EfrkZ5DtdeBQFbR4FHP9mjpfI9am9bT3ahzYOYcLDu7CzEXr+Xh5IWf1Sifd6ap2qQgA20b78HUuHX0Rl838kQ3lYbI9BqP7t+WaWT9ibn09xkybt5cVs2iLmwcvvQ2mXFZ1WUjfGlRzcrFPmsAP7Qbw+GdVM4Ou3VjBuN5dOXXM3/AEKwgt+Bj3gYdT7Mriojd+YVPF9rq88uNmPl5VypNHtydtyGEQjRKYu0Ovg+HAriyr+4/oTwZ+W5UVpI85E3PLJhRVq/l72PGuwiGMbvvQ6sDDKb73JojHq5byz8oBILryF5z7DiS2anmN2/rfn0XuX4bj7LcfkcU7DMi1TKK//YpzwGAiPyyscbu0k8/E9FfiO/HUqnFHioL/o7p7XbQ2bYk0088PCSJ7gaIomA4fLy3exMxFGwjGTPbNS2PSYV3Jc1nY8Zazh0DKc6Xx8uJNvPbDesIxi3SXzvj9O9A5x8vyLZUM75mLx45gxRtmxLrucBJSXCzeXEmhP0TfthnkenTUSGW1QKJqOmHNzUcri/hhfQWdsj0c1zuPcNRE11Q0VeGN7zfwxeoSfA6N0wa2Zf88D+G2XViVXsCvWyrJsxWi/hCz5q3gxqP34cnvVvPe0s04dJV7T+rLW4s3MGfZFkzL5tCuOZx7UGd+WFnG2tIgHbO95KW7eH/pZtaUBvlL1xwuOqQLX60uoXvrNH4trCTL7aBnmzT+t7qEzRVhBnXIpEsrL4X+KD9vKqdzjo9VxX46ZXv5vSyEqij0apOGbdvc+t7POHWVMw8ooF2mm399soJR/drSJceDoWnM+7GQlUV+9s1P58geuURjNmtKgizdVEH7TA+DO2Xj0FTmryzmi9+KyXAbjO3fjgyPwZcrivho2RZ0VeWkfvl0z/WxotDP+z9vJm7ZHN49l33z0tlcGaZ1upPTDyjgmS9Xc2i3Vnz5Wwmzl24iErfI8To498BOdG3l5ffSEL3y0thYHiLd5SDDrfPa9+vpkO2hTZqT139Yjz8S55oje3DtGz8SjFVdyjEtm9cXbeDL1SXcP7ofFWGLS1/7odrll1e/+51bR+zDzxsrePn77Xvc7Jufzu0jejN9/homdt0HavkwA1A7dePTdQGeW1yIqsB1w3ti2TZf/F7JiZ26EV32U43bWOt+Q4tH2VBe9eE8sm9bXvx6XSKE7GhtaYgVWj49W+djFm5Cb9WG3EdeYnEFPLa4jJ8+qz49+eWlJZw0agQ5pesomnINGRdezQcbw9VCyDZFgSgfrA0w/sAjMMtLsMM7nBOLwh/25NmRWbgJvUPnquXta2G070Ro4ee4DzkKs2gzroFDa65VspX7wCOIrfuNkvtu2r5uSKduictM4W+/pNWtDxGYNzsxsDohHgeHk8xzLmfLDRdWGydS+coztJryKMG8dgQ+frdqBpIvHd+Jp6B6fYQ+eR9UlcC898i95QG0zMxa66c4nei9+ia2BGhumt303Q8++IC33nqLsrIy2rdvz9lnn80+++yzW/exp9N3LYeXye8tZ+mm6i9YBXj81AG0c8Ya/AUlU0TrweHm0S/XM+eXLTWKTt+/A4WVET5ZUcgDo/vS0WNhmbt2meOPtrVNYXEJv5ZbXPvmT0Ti27/Z9M1P544TeqMEy4CqSyLFlpOLXlmU+EAD0BSF24/vTY7PwVUzFyfGLQC09jm584Te3PTuUjZXbl/Gu12mm/tO6sv1b/7EmtKq8RnXD+/JzEXrWb6larnxnq19nD20E7fOXlqtXplug9uP780TC1Zx9pBOvL90E/sXZHHfx7+yf8csju7Vmrs+WJb48OqTn84pg9rz0jdrOX3/Ar7/vYy2GW4e+3xVtQ+4MQPaclLftpw5o2r9jhP2zWN0/7Ys21JJx2wvV7y+mOjWfT32aePjtuP35cqZi/m9bPtAUq9DY8oJ+/Lo5ytZUVj1xn/fSX15bP4qVhZVX530kC45HNg5h3s/3v5ttndeGv84bh9+XF/OzEUbOKx7Lt+tK2XhmtIa7XfBwZ1ZtL6cr1aXcM6BnejZxkdZMEbHLDdLN/t56JOqBceO650H2Ly3dHOtr4Mbhvfk+9/LeP/nmuWGpvD8mQdw9gvfVGuD4/fNo3cbH8e1sii/fHzNvWgUBeX2Jzh7gZ/yreFGUxXuO6kvv60vYvjrd9Z6acNz/F/5bMg47ppbNQvmzhP25eZ3l2DV8dYxcp9cLvz2GexYDM8RxxLsdzDHP1vHZRHgpqO6M2TGzcSW/YTjmru5/vdMftpYcywHQLdcL4+c0AN11rPE1qwg8v32KdfpZ15E5IeFtU799Z54Ku5BB1J0x1U1dib2HjMa1eulcuYLKB4vubc+iB2NUvLg7ZhF1X//anYrsi+7maJbJ1U7nn3NFMqf/3diHRijY1cyzr0c/zuvEv5m6z40g/9CxmnnE12/FnPLRly9+xOc/xGRn75DTc/Ee8RxmJWVVPY/mDRDwRmPoFg2le++SvCDWUDVCqoZEyaipmdgh0JUznyh2rRzxe0h++b7CeR1IJYiQaShp+82qyG4X3zxBc8++ywnn3wy06ZNY5999uGuu+6iqKju5XwbmqIobArEa4QQqBr29MAnK4lrdad80XjCllZrCAF4/Yf1DN+nDTHT5qZ3fyaqufb48QI4ueaN6iEE4MeNFTy3cC0Od9UOpzHNxW3vLasWQqDqcsGPG8r5v+9+rxZCAM4aUsCdc5ZVCyEAlmWzdFNFIoR4HRrpLj0RQgDOHNyRu+f8UqNeZaEYj3y2kov/0pXpX/zGyL5tuffj5Zi2zakD23P3h8uqBYxx+3Xgvo+Xc8qgDjw2fxUHdc7hkc9W1viW/foPG1hVHOToXq0BeGfJJmxgcMdsbnx7SSKEAFx/dC/u/Wh5tRACEIia3Pbez0wY3BGAQR2qPuj+GEIA5q8qxuPQaOXdvn380k2VfL6yiF55PpYX+umR66s1hAD895t1nNAnHxt4+svVuA2N177/Ha/ToH2mO3Fe77w0vq7jPgA+XVGIU6/9LTdm2vy8qZLR/dpWOz7nl80M7JDFrPUm6Xc8ipbXLlGmZmajXXEXz61TEiEEqnphVhUH+EvPfGK//VrzwTQN13FjKYpsfw0FYyYZbqPOurdyVK0MmnH6eeht2hGO7LxXV4sEybroGhz99sf6bTkuo+bGgNu4dA1r4zriG9aSNup0MLa3U+Vrz5E25kxc+x20/QaqivvQY3D1HYT/vddp9Y9/4hpyKGpWDkaXnmRfdRuOXn0Ss1HsYIDApx+g5bWj1c334R0xBjU9EzU9A+/xfyV3yqP43355+92nZ5B5xa3o3XujeryJ47F1vxH66lP0iZMpnTID/9QXiZ8xiS2rVuH/YSFG3/0on/ki8c0bcB/wF4wOnal47XkqC3py0/zNXDx3E4GMXPwfvYWjQydyJt9Dzo334DvmJMqnP0DxlGswi7aQefF1tH7webKuuIXsWx8k+8EX8LdpnzIhZG9oVpdm3nnnHYYNG8aRRx4JwNlnn82iRYuYM2cOp59+eqPUwTAMFv5a95vRL5sridkqjjrPEI1BUZQaH9o7isQtzK1fD8tCMcrDJjl7GNsXrS+v9iG7o9lLNjF+v/YYQDAOq4prfqAC9GuXwc3vLq1xPD/dxeqSP25dDt1zfXy7rmz7eRmuah/WhlY1bmLHD7IdLd/i3zpGQmHxhnIsG/LSXawvDxEztwcMBXBoChXhOF6nTu+8dD5eXljrfQK88t06LvpLl0QQXLSujL7tMykPVe+JdGgqX6+t/e+pLBRDURTchsYhXXKYuaiWFUy3mvdrIUM7ZfPOkk2JY+/+tImhnbLJchusLw/VedvKSBxd3T6Y95Vvf2fcfh2oCMcIx0ycukokbhGJW3idOkWB2j+k01wG4Xjd1/grIzH6t8/g5e+2L0ceM23ils3/1lXQfXABwb/fTWdHjGy3zi9BlSeXVPDjppq/n0jcYnVpkO6T7oDp0xLrWKjZrXBe8g8+KtXYNy8NTVUwLZsPf9nMqL5tefar2seiHLVvO+h+NtGfFlYttmbF6N82nUUbavZyqAr0yXWzRVWYO/wShrb1cKri4Zs62vHUAflYL9xG+slnUjnzeVrddC/+t/+P8A8LqxZmW7qYzPOvxD7zQuxYDNuyqr6FR0JElnxP5Kfv8Bx+LOnjzkVr0xa9TVuCH7+Lo3tvFK8Pz+HHYVeWU/rwnaSdci7uQ4/Gd+IpEI1ilpdUTRu/+Hp8wQC2GQePj4DhImBZeG+8Fz0UwI6EUXxpRJxu1kdUzp2z/RJax+xs2rc/kYpvwxz5lwkc385B+aqV4EuncOhIHvyhjF+29tpVWAptDjuGoqv+VuvvQsltQ4nmwPJmovUdjG3bVePBmunYkG2aTRCJx+OsWrWKk046qdrxfv36sWzZslpvE4vFql2CURQFt9ud+Hd92LZN5k6+WTh1FVWp//3XZdv9NfT9Nmdprp2//B3a9t9lzLKr1l2oh21tsrmy7sFykbhF3LJxKEqt1+m3seyqSzQmO4wn2Tq4tDaBaJx2mdt7cyrDcXJ26BlwaCqh2J9/0/I6dEqDVR+wHkOjMlI9uCgKxLf169s2GW6DdaU1g9E2Rf5otQ93j1MlXEs9wvGd160iHMNtqKiqsv3xaxE3bVS1evtFTQtVUQjFzJ3+zSpQra5FgShpLh1FodpuMnOXb+G43nk8Pn9VrfdzQp88Hv6k5r4x23TO8VLorx6OW3kdqIpCt1wv4ZjJPxZsJmpaXHlEd2YuWl9r+OT/27vz6Ciq/O/j76rqJd0J2cm+QggYlpBISAiMLOowioAow6Ayj8wjOCzqMzoOoh7nJ/wYefz9RgaPyzyOg+jAEUUIy09HR3BADaDIosCERZaAE5YEQtZOr1XPH02ahAQVBDvA93VOH+hKddVN307Vp++9dQvokxTOjJW7yIqN5Ne/e5E4nESFWjjgNGGK6cy8ZTvIT4nkP0fk8PTfd/PF4dPc2TeZ/JRItv27ptXv/tCQLLYerefGFBvWvgX4XG7QVB67sRv3L/2KhnM+C//nhi7UYGHK8t1E2EykpcTh1Q2KM6PZeKi61bpF6VH0To7Ec89DOLZ9TKex99K45n8wpWQQc9NI1E4RYLPhqzmFGhOH7mgEswVPiA0XKrF/+htUn8RwuzDFJeIu/5qqmb/GnJmN/cYRWPsWgsmE9+A+7DcMRw2xopgt6I0NKCF2jJQu1KomvF4v2DqdLdiZ1ocGRQN7uP8BoIPNoqFA4C/wcLWDw2fqYefROlKSe/NfB8NpdHtpdJ8NLJqioKDgiowhYtpMal/577MT4Kkqne6+H19SGoZhnB2QTsc8pl/q881VE0Tq6urQdZ2IiIhWyyMiIqipqWn3NStWrGDZsmWB55mZmTz77LMX1cfV0iDNwfz1B2jvsDiqdyIJkWFYTOE/aB/nk5CQcFm2ezXS6l0khIe0O4iuT3IEe074uy8smkpChI3EiOgftL/c5Mjz/iw50kaIWSOxcyJavZOYUAun2vlmvenQKX56XRzv7jr7zV43wGZufXBs9mVFLRMK0liy1f8t+0S9i85hVkLMKk6PTqPbR0yopd3XAoRZTfh0A4tJpVdiBCVfHaWitomszmGt1tMNMGtqoHXgRJ2T6+LDz9vd0TOxU6B7SFUgq3M4ETZz4Bt6M7vZf9XKuSe7ZgnhIdQ2edl2pIafdI3hnRaDPVsa1DWGN7e0vixycFZn/9VR+A+oETZzmxYZgKKM6FYn5z5JEXxT3URhZjSNLl8gBJYdr+eu61PJT41kW4tWKICf5yVzuNrBhP7p/Md7ZW3e61G9EzntcFN6oPUN1SYXZ7J6x1EGZcWyZs8JbsmJZ9XOYyz/qoJJAzLabR3LT40k2m6hT3IEJ+qdrD7qY0TPVOZt+YbCjDByVZVZI64jLcqOTzdY/L8KOHCygepGN0/8tDvVTR42HDxJiEmjb0qEv7VLgcX/OsG2o/V8XdmIzzDI6hzKgrvz+ee+Kr44cpr4MCs/z08m2m7hX8fq+b+jepEcGcKs93dz4GQjDw3O4mc5CZQeOIVhGIzs7R9IbDNrzC2HuPifkI+V6NGTsaoKoSaDet1ETLiN6PjWx2RbyyfJaYH/mpJSseb0BV1HsdnQwvzHWc0Wit5Qh+F0oIZHoNmS0ML8wcPe5h38dvVOL4O6xvLpgbZd/unR9sCVT+ca1r0zncNDsVvC0W8cgf36AbjLvwafjjmzG1pkDKrdzuU5M1wel+p8c9UEkWbtJbTzpbYxY8Zw2223tVmvqqrKn5Avktlk4fc/68HsD/a0OuBkxYZyd34yp6raH8z2QyiKQkJCAsePH5fBqt+Tpmn8cXRPHly+g9OOsyeg5Ehbq4P8lIEZqG4Hx47VXNR+musmOcJKVufQwODKlqYOysSmOzl2rBbNZOJ3w7KY+T9tTzJRdgvDsjvzeXk1VQ1ng8rmw9X8tEcc/zhnzItPN6isd/LYTdn819p9GMDfNh/mP27JYfb7u2ny+Nhw8BQjeyeyeuexNvubWJjOqh0VTCxM53idk4xoO+XVDr457aB/elSroPHurmP8qjCdd76s4LaeCcR3CmH5VxVtQoSmKtxbmMEDS/2DHScVZ2LWFHy6zj39Uvnb5rN3kN3+79P8qiidFz5uO7X4wC4xlB2rw2cYlB48yfw7c1m3r6pN10j3uDBCLaZW40yi7WZu65XAMx/uYc5tPVmw6RC//1kP/uPvu1uVNyPazl39Unl8tf/KkxCzypjcJPZX1uNwe9t06cz5xx4eHpLFfQMy+OfeKqxmlZu6x7Hx0Cn++NHXjOyVwF/vzuevG8spO1FP5zALt/dJwmrSSI60seGgP4jEhFqYXJxBcqSNuWv2sqeqgV8W+E+4Do+Pj/ZWUnrwFM+M7MnfNh9hz4l6wqwmRvZKYEyfZJxeH+PykokJtXKszsk31U08NCQLt9eL3aIRaQulyaPjOzMXTGZMKFaTiqZAWkQIkTnxmDUNRQGTprD+65Ms/ep4oNVJUxV+1iOBXcfq6J8Rxc9y4rGoCl6PF6vuIqdzCLVuA6fby4ODs5j74R7mrfuaSJuZW3LiuaNPEqG48NRX41UU/ndROve//SWLdrYOghMKUvl5VDjHjrX9bH6n+kb/o5lqBnsEeA2ob/A/LoKiKDwypAt1TjdfVZztmkqPsvPfo3uy+1jb8YEJ4SFMKU6nrrqK2sDxWUXN9F9Ioes61Nb6H1eA851vTCbTRX2Rv2qumvF6vUyYMIFHHnmE/v37B5YvXLiQ8vJyZs2a9b23dUluemey4MRC6aFqqhvdFGZEkRphRXM3XpagIFfNXBxN03BpNr6pdfHvGifp0XacHp0XPjlAqEVj0oB0siLN4Dn/GILv0lw3lZWVuMxhvFJ6kI/2VeHVDRLDQ5g6qAu5iaHgbHEAM1mpdCm8srGcfZUNJISHcG//dKLsZv60bj9TBmVy5JSDTw+dItyi8vOMEJKjQnlnXy3Lyk7i8urYzRrj8lMY2TuBmkYPJpP/clenx8dN3Ttj1lQOnGykpslD//Qotn5Tw+ufHaaywUVGtJ0pg7qQ1TnUPzGYpuA+M9nXW9u+4ZP9J3lyeA92H69n5Y6j1Lu8dOscxmM3Z3Oq0c2+ygb6p0Wio/DypwfYeWYsQWaMnd/dmM3Oo7UcPOVgdJ9EIkPMfLD7OIO6xBIXZmXrv2tYtPkIFbVNZMaEMndkTzYfqWHBxkOcbHRjN2uMyU3izr7J/L9PD7L5yGkibGbuLUwjNzmSlV8d5Z9fV2HRFEb1TmJIVixr951g2fajeHWDYdmdGdU7kf98fzfXJYRza04CNotGg8uL3WLiRJ2T4/VOusaGUe/y8uyavdQ0echJ6MSjN2YTZjWhYrC3qpGIEBMen8Fb2/5NjcNNUUY0w69L4NCpBjKjQ3HrOo+v3kVlg5shWbHcNyCDMKtGrdPLiToXdU4vNrNKZkwooVaNJo9Ok9uHSVM47XCREG7j66oGDlQ1kpscQWJECB6vgY5Bk8eHzazS6PJR6/Ti8elE281EhJjx6AZhVhMur45J9U8+F2JS8fgMdh6ro+FMfakKNLi8PPVeGV1jQ3lwcBY+n46BQVJECIq7CXwesIbh1JUzg54V0qPtuL06mmJgU32oHme7xx1FUdA0DbdqpdFroAB2k4LZ52x15aDJdOZy9b1VbDhUTaTNzC/yk8mOD6ep5lSHO6YpioJuttHgVThe7yTGbiHSqmLyNqGrZhoMEx+UVVLZ4GJQl2h6xodh9jqumnmkLvVVM1dNEAF44okn6NKlC5MmTQose/jhhykoKLigwaqX8u67zRNQeb3ey/ohlCDyw6iqGrhNulez4ENFM3RMuvsHX2rdsm4A1JAwHD4Fj+4/OYQpbpxNbYOOqqroJiteVDR0QlRwK2bcuoLT4yXMomFRdGw+D2r1CUBBj4nntBqCy+vDommYNQWXT8du0jDwT+KlGwaqAigKbq+OfmYyshCThsunY+DvZtFUcHp1zIriP6Eo/rESXsOgyaPj0XVCTSZ8hoHPMDCrCgZnx16Af5yN12fg8un4dAi1aoRZNBrcPv/4CsPApCr4DP84E5OqYlL9g2d1w//Nu5PV3+1U6/Ti9vn3E2Uz+SdiUxUcHh8KCjazisPtPzk7vTqKohBqUXF5DTpZNOrcXgwDTIpCg8fnL6PhHzsSZjGhGzpmVQ28Tz7dwKMbOD3+98RmVgkxqSgoeHXdX2D87yf4L0EMDzGhoPgnQHN68RoGbq9/H+FWE7qu4/aBpvrf4+ZuHYuq4lN0dJ+B0qJ7yuH2YTerdAox49V1zJpGg9OLzaLh0/1hxOH24vL6CLOY0BQF3dCx6C5cLn/3gKb5r1hp/hybzWY0kxmfrqOg4DA0DAPMio7mcwWm6W/vc6+q/hHbl+tYZjKb8SoaGoDPQ3x8fIc/pp3vtgdWqzVw7P8hLewdkdx991vcdtttvPDCC3Tp0oXs7GzWrl3LyZMnufnmm4NWpksVaMTlpev62cFhvqbAH8alvmDOMAx8TfVYASuAF843hFXXdXCfLYv/tNKEGTADOMGD/0FY5JmVnGg4/f3eHv+4D4t/NwEtx4O0HKLp4ez1/PqZh3bmX/C/F82dHtqZh8919nnzei236Tuzv8BwWS80nGktby6D55z1XbQ+MDnOZDRTi+Utu+Cb9+dx+v/vdZ1dr/kednVn9qmc2UfzGINzD6jn/rUqLdbF01wH7dOB6nZa+5vf07pzeuRaRs/zbdcEuJ1w6pzW/pafmcDnoens57Xl9s4NFOcO0m/5Xn/X5/1yf6P3evyfaB8dc5Bme84XkpqDoPhuV1UQKS4upr6+nuXLl3P69GlSU1N5/PHHf/DgUyGEEEJcHldVEAEYPnw4w4cPD3YxhBBCCPE9XFUzqwohhBDiyiJBRAghhBBBI0FECCGEEEEjQUQIIYQQQSNBRAghhBBBI0FECCGEEEEjQUQIIYQQQSNBRAghhBBBI0FECCGEEEEjQUQIIYQQQSNBRAghhBBBI0FECCGEEEEjQUQIIYQQQXPV3X33UjCZrty35Uou+9VO6qbjkrrpuKRuOq5z6+Zi60oxDMO4FAUSQgghhLhQ0jVzlWhqauKxxx6jqakp2EUR55C66bikbjouqZuO61LXjQSRq4RhGBw6dAhp4Op4pG46LqmbjkvqpuO61HUjQUQIIYQQQSNBRAghhBBBI0HkKmE2mxk7dixmsznYRRHnkLrpuKRuOi6pm47rUteNXDUjhBBCiKCRFhEhhBBCBI0EESGEEEIEjQQRIYQQQgSNBBEhhBBCBI1M4n+FW7FiBZs3b6aiogKLxUJ2djYTJkwgKSkp2EUTLaxYsYIlS5Zw6623MnHixGAXRwDV1dUsXryYL7/8ErfbTWJiIlOnTqVLly7BLto1zefz8c477/Dpp59SU1NDVFQUQ4YM4Y477kBV5bvzj6msrIzVq1dz6NAhTp8+zaOPPkr//v0DPzcMg3feeYePPvqIhoYGunXrxn333UdqauoF7UeCyBWurKyM4cOH07VrV3w+H2+99RZz5sxh3rx5hISEBLt4Ati/fz9r164lPT092EURZzQ0NPDUU0/Rs2dPnnjiCcLDwzlx4gR2uz3YRbvmrVq1ijVr1jB9+nRSUlI4ePAgL7/8Mna7nVtvvTXYxbumuFwuMjIyGDp0KM8991ybn69atYr33nuPadOmkZiYSElJCXPmzGH+/PnYbLbvvR8JIle4J598stXzadOmMWnSJA4ePEhOTk6QSiWaOZ1OXnjhBX79619TUlIS7OKIM1atWkVMTAzTpk0LLIuLiwtiiUSzffv20a9fP/Lz8wF/vZSWlnLgwIEgl+zak5eXR15eXrs/MwyDv//974wZM4bCwkIApk+fzuTJkyktLeXmm2/+3vuRdq6rjMPhACAsLCzIJREAf/3rX8nLy6NPnz7BLopoYcuWLXTp0oV58+YxadIkZsyYwdq1a4NdLAH06NGDXbt2cfToUQDKy8vZu3fveU+IIjgqKyupqakhNzc3sMxsNpOTk8PevXsvaFvSInIVMQyDN954gx49epCWlhbs4lzzNmzYwKFDh5g7d26wiyLOUVlZyZo1axgxYgRjxoxh//79LFy4ELPZzODBg4NdvGva6NGjcTgcPPzww6iqiq7rjB8/nkGDBgW7aKKFmpoaACIiIlotj4iI4OTJkxe0LQkiV5EFCxZw5MgRZs+eHeyiXPNOnjzJ66+/zpNPPonFYgl2ccQ5dF2na9eu3H333QBkZmbyzTff8OGHH0oQCbKNGzfy6aef8tBDD5Gamkp5eTmvv/56YNCq6FgURWn1/GIma5cgcpV47bXX2Lp1K7NmzSImJibYxbnmHTx4kNraWmbOnBlYpus6u3fv5oMPPuDNN9+UKwCCKCoqipSUlFbLUlJS+Pzzz4NUItFs8eLFjB49moEDBwKQlpZGVVUVK1eulCDSgURGRgIErmxqVldX16aV5LtIELnCGYbBa6+9xubNm3n66adlwF0H0bt3b/74xz+2WvbnP/+ZpKQkRo8eLSEkyLp37x4Yg9Ds6NGjdO7cOUglEs1cLlebvw9VVS/qm7a4fOLi4oiMjGTHjh1kZmYC4PV6KSsr45577rmgbUkQucItWLCA0tJSZsyYgc1mC/Tb2e126RIIIpvN1macjtVqpVOnTjJ+pwMYMWIETz31FCUlJRQXF7N//34++ugj7r///mAX7Zp3/fXXU1JSQmxsLCkpKZSXl/Puu+8ydOjQYBftmuN0Ojl+/HjgeWVlJeXl5YSFhREbG8utt97KihUrSExMJCEhgRUrVmC1Wi94PI/cffcKN27cuHaXT5s2TZoxO5inn36ajIwMmdCsg9i6dStvvvkmx48fJy4ujhEjRnDTTTcFu1jXvKamJt5++202b95MbW0t0dHRDBw4kLFjx2IyyXfnH9O//vUvZs2a1Wb54MGDmT59emBCs7Vr19LY2EhWVhb33XffBX/ZkiAihBBCiKCRjmohhBBCBI0EESGEEEIEjQQRIYQQQgSNBBEhhBBCBI0EESGEEEIEjQQRIYQQQgSNBBEhhBBCBI0EESGEEEIEjQQRITq49evXM27cuMBj/PjxTJkyhZdffpnq6mrAPwPiuHHj+Oyzz9rdxoIFC847C28wjRs3jgULFgS7GD+6kpISNm/eHOxiCNEhyHy5Qlwhpk2bRlJSEm63m927d7Ny5UrKysra3FxPdHwrVqygqKiI/v37B7soQgSdBBEhrhCpqal07doVgF69eqHrOsuXL+eLL74gOjo6aOU6ffo0VqsVu90etDJcLdxuN3V1dcTGxga7KEL8aCSICHGF6tatGwBVVVU/ehBpaGhg8+bNbNiwgV27dvHss8+SkZEBwEsvvURZWRkvvfRSq9csXbqUZcuWsXTp0jbbW7NmDe+++y5VVVXEx8czduxYBg4c2Gqdmpoali5dyrZt2wI3QxsyZAh33HEHmqYB/ruDPvDAA0yYMAFVVXn//fepq6sjLS2Ne++9l+zs7FbbPHDgAMuWLWPPnj243W6Sk5O5/fbbKS4uBqC8vJwZM2YwZcoUhg0b1uq127dvZ+7cucyYMYN+/foFfr/nnnuO5cuXs337diwWC3l5eUycODEQ1Jq7yD7++GM+/vhjAHJycnj66aepqanhwQcfpEePHgwcOJCioiLCw8MvpoqEuGJIEBHiCtV8e+6WJypd1/H5fG3WvRT3tnS5XGzdupXS0lK+/PJLVFUlLy+Phx9+mJSUlIve7pYtWwJjXKxWKx9++CHPP/88mqZRVFQE+EPI448/jqqqjB07lvj4ePbt20dJSQlVVVVMmzat1Tb/8Y9/kJycHLjT8dtvv83cuXN56aWXAoFg165dPPPMM3Tr1o3Jkydjt9vZuHEj8+fPx+12M2TIEDIyMsjMzGTdunVtgsj69euJiIggLy+v1fLnnnuO4uJihg0bxpEjR1iyZAlAoIxz5sxh9uzZ9OzZkzvvvBMgUKbY2FgeffRRNmzYwKJFi1i4cCF9+vRh4MCBFBQUYLPZLvp9FqKjkiAixBWiOWR4PB7KysooKSnBZrPRr18/KioqAJg/f/4l3afX62XHjh2UlpbyxRdf4PP56NOnD1OnTqWgoICQkJAfvI/6+nrmzp1LZGQkAPn5+fz2t7/lzTffDASRpUuX0tjYyLx58wLdFr1798ZisbBo0SJGjRrVKgzZbDZmzpyJqvrH40dFRfHEE0+wffv2QEvLggULSE1N5fe//32gRaVv377U1dWxZMkSbrjhBlRVZciQISxcuJCjR4+SlJQE+FuEtmzZwvDhwwOvbTZs2DBGjRoFQJ8+fTh+/Djr1q1j6tSpKIpCdnY2iqIQHh7epoVGVVUKCgooKCgIBL+NGzfyyiuv8Je//IX8/HwGDRpEXl4eZrP5B7/3QnQEEkSEuEI8+eSTrZ6npaUxadIkIiMjA0HknnvuoVevXm1eu3r1ajZt2nRB+zt8+DCzZs3C4XDQu3dvfvWrX1FYWEhoaOjF/xLt6NWrVyCEgP9kPGDAAJYtW8apU6eIiYlh27Zt9OzZk6ioqFYtPnl5eSxatIiysrJWQSQ/Pz8QQgDS09MBfzcW+FuTKioq+OUvfwnQapv5+fls27aNo0ePkpKSwk9+8hMWL17M+vXrufvuuwHYsGEDHo+HoUOHtvl9+vXr1+p5eno6Ho+H2traVr/nd7FarRQXF1NcXIzD4WDLli1s2LCBP/3pT1itVmbOnEmPHj2+9/aE6KgkiAhxhXjggQdITk5G0zQiIiKIiopqs058fHxgQGtLFzPOQNM07HY7DQ0NOBwOHA4HTqfzkgeR9k7Ozcvq6+uJiYmhtraWrVu3ctddd7W7jbq6ulbPw8LCWj1vbj1wu92Av6sHYNGiRSxatKjdbdbX1we2df311/PJJ58wfvx4VFVl/fr1ZGVlkZqa2uZ137Xvi+FyuWhsbMThcKDrOjabDZNJDt/i6iCfZCGuEMnJye2GjMslJSWFF198kX379lFaWsqqVatYtGgR3bt3p7i4mKKionZDhNlsxuPxtFnefGI/V3MoaG9Zp06dAv+mp6czfvz4drfRXij7Ns3B7Pbbb6ewsLDddZq7YQCGDh3KZ599xo4dO4iNjeXAgQNMmjTpgvZ5oerq6vj888/ZuHEjZWVlhIWFUVhYyF133cV1112HoiiXdf9C/FgkiAghvlV2djbZ2dlMnDiRnTt3UlpaypIlS1i4cCE5OTkUFxczePBgLBYLAHFxcdTW1lJTUxMIKl6vl6+++qrd7e/atavVurqus2nTJuLj44mJiQH83SXbt28nPj6+TYvDxUhKSiIxMZHDhw8Hulu+TW5uLtHR0axbt47Y2FjMZjODBg266P2bzeZ2W0i8Xi+ffPIJmzZtYteuXZhMJvr168eMGTPIzc2VVhBxVZJPtRDie1FVldzcXHJzc5k8eTLbtm2jtLSU119/nW7dugUu3y0uLubtt9/m+eefZ+TIkXg8Ht5//310XW93u506dWL27NnceeedgatmKioq+M1vfhNY5xe/+AU7d+7kqaee4pZbbglM7FZVVcX27duZPHlyILR8X5MnT2bu3Ln84Q9/YPDgwURHR9PQ0EBFRQWHDh3ikUceafW733DDDbz33nvYbDYKCwt/0LwpaWlplJWVsWXLFqKiorDZbCQlJVFdXc2rr75Kbm4u06dPp6CgAKvVetH7EeJKIEFECHHBLBYLRUVFFBUV4XA4Wg0MjYuLY8aMGSxZsoR58+YRFRXFiBEjqKurY9myZW221a9fP1JTU3nrrbc4efIkCQkJPPTQQ4G5PMDf9TJ37lyWL1/O6tWrOXXqFDabjbi4OPr27XtR41Z69erFM888Q0lJCW+88QYNDQ106tSJlJQUBgwY0Gb9oUOHsnLlSjweD0OGDLng/bU0ceJEFixYwPPPP4/L5QrMIxIVFcWrr756SVp9hLhSKMalmGBACCGEEOIiyE3vhBBCCBE00jUjxDVG1/XvnGn13Em6hBDicpEgIsQ1ZtmyZe2O1WjpxRdfJC4u7kcqkRDiWiZjRIS4xlRXV3P69OlvXSc9PV0uFRVC/CgkiAghhBAiaGSwqhBCCCGCRoKIEEIIIYJGgogQQgghgkaCiBBCCCGCRoKIEEIIIYJGgogQQgghgkaCiBBCCCGC5v8DRcM6mkEHhpUAAAAASUVORK5CYII=\n",
"text/plain": [
"