{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"\n",
"[*NBBinder test on a collection of notebooks about some thermodynamic properperties of water*](https://github.com/rmsrosa/nbbinder)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"\n",
"[<- High-Dimensional Fittings](04.00-High_Dim_Fittings.ipynb) | [Water Contents](00.00-Water_Contents.ipynb) | [References](BA.00-References.ipynb) | [References ->](BA.00-References.ipynb)\n",
"\n",
"---\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Choosing the Best Fit with AIC\n",
"\n",
"In this final section, we use the *Akaike Information Criterion (AIC)* (see [Burnham \\& Anderson (2002); Bender (2000)](BA.00-References.ipynb)) to select the most efficient polynomial approximation."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"## Importing the libraries"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"import csv\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Akaike Information Criterion\n",
"\n",
"Beink $k$ the number of parameters in the model, $N$ the number of data points in the sample, and $E$ the mean quadratic error with the approximation with $k$ parameters, the **Akaike Information Criterion (AIC)** is given by\n",
"\n",
"$$\n",
" \\text{AIC} = N\\ln(E) + 2k.\n",
"$$\n",
"\n",
"The first term decreases with $E$, but the second term penalizes a high number of parameters.\n",
"\n",
"There are a number of other criteria, such as *Bayesian information criterion (BIC)*, *Mallow's $C_p$*, *AICc* (for a low number of samples), *adjusted $R^2$*; *ridge regression*; and *cross validation*."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### AIC for the approximating polynomials\n",
"\n",
"In our case, a polynomial model of degree $j$ has $k=j+1$ parameters. Denoting $E_j$ as the mean quadratic error associated with this polynomial, we find the $\\text{AIC}_j$ for each polynomial as\n",
"\n",
"$$\n",
"\\text{AIC}_j = N\\ln(E_j) + 2(j+1).\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Loading the data from file\n",
"\n",
"For the computation of the $\\text{AIC}_j$ values, we first load the data from file.\n",
"\n",
"This time, we use `csv` to read only the header of the file, contained in the first two rows, and then we use `numpy` to read the temperature and density values, which are in the remaining rows. This is a more efficient way to load the data."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Loading the header\n",
"\n",
"We start by reading the header. We read the first two lines of the file and make a dictionary from it."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{'temp': 'Temperature (C)',\n",
" 'density': 'Density (g/cm^3)',\n",
" 'viscosity': 'Viscosity (cm^2/s)'}"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"water_csv_reader = csv.reader(open('water.csv',\"r\"), delimiter=\",\")\n",
"line1, line2 = next(water_csv_reader), next(water_csv_reader)\n",
"header = dict([(line1[i], line2[i]) for i in range(3)])\n",
"header"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Loading the data\n",
"\n",
"Next we read the data values for temperature and density, which are in the first and second columns, so we skip the first two rows and ignore the last column."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"T, f = np.loadtxt(open('water.csv', \"r\"), delimiter=\",\", skiprows=2, usecols=(0,1), unpack=True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"We also define two auxiliary variables related to the number of rows."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"N = len(T)\n",
"N_half = int(N/2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Computing the mean quadratic error\n",
"\n",
"With the data in memory, we compute the $\\text{AIC}_j$ values."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"j=0: Error=1.75e-04\n",
"j=1: Error=9.22e-06\n",
"j=2: Error=1.33e-07\n",
"j=3: Error=3.16e-09\n",
"j=4: Error=3.27e-10\n",
"j=5: Error=2.64e-10\n",
"j=6: Error=2.64e-10\n"
]
}
],
"source": [
"A = list()\n",
"Err = list()\n",
"for j in range(N_half):\n",
" A.append(np.vstack([T**i for i in range(j+1)]).T)\n",
" Err.append(np.linalg.lstsq(A[j], f, rcond=None)[1][0]/N)\n",
" print(f'j={j}: Error={Err[j]:.2e}') "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Computing the AIC values"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"j=0: AIC=-125.74\n",
"j=1: AIC=-167.91\n",
"j=2: AIC=-229.52\n",
"j=3: AIC=-283.59\n",
"j=4: AIC=-315.63\n",
"j=5: AIC=-316.81\n",
"j=6: AIC=-314.83\n"
]
}
],
"source": [
"AIC = [len(T)*np.log(Err[j]) + 2*(j+2) for j in range(N_half)]\n",
"for j in range(len(AIC)):\n",
" print(f'j={j}: AIC={AIC[j]:.2f}')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Plotting the AIC value\n",
"\n",
"Finally, we plot the AIC values in terms of the degree of the polynomials."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAFRCAYAAAAirkrIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZyNdf/H8ddnxhqJogVJy20pRCZbiFJURCKksiSptO93d90q7X5td1q0WLKmktIiQkVJgwmRUqksRbpt2c3398f3GvdxzHJmzMx1Zub9fDzO45zzvbbPtZ7P+X6vxZxziIiIiEjBkBB2ACIiIiISOyVvIiIiIgWIkjcRERGRAkTJm4iIiEgBouRNREREpABR8iYiIiJSgCh5kzxlZrPM7Ll8mM4IM5uS19MpqMysv5n9amapZjYoG8MNMrMleRha5LRuM7OV+TGtMJlZLTP70sx2ZGd+zay6mTkzS8rD8AqksJeNmfU2s61hTDsn8uu4nBtycgwys5VmdltexRQPlLwVQGbWwMz2mtmcDLo7M+sSVVbfzCaY2e/Bj8aKIOGpmz9R57kbgcsOZgQF7QAcKzOrAAwFngCqAEMy6O+A7UbyxGBgG1ALOD29HvRnJNt+A44BUvJ6QhnsJxOAE/J62kXUEODMsIOIN0reCqargOeBOmZWO6uezaw98BVQFrgcqA10B9YCj+ZhnPnGObfJObcx7DjSmFnxsGOIcBxQDJjinFvrnCt0CWp2xMG6OQmY7Zxb6ZxbH3IsoTKzErkxHufcXufc7865PbkxvhxMf7tzbl0Y0y7snHNbnXMbwo4j7jjn9CpAL6A0sBGoB7wKDEmnHwd0CT4fAqwH3s1gfOUzKL8a+AMoFlU+FpgcfD4RmAz8DvwNLADaR/U/C3gu4vtK4LYs+ikBPAasCsb7NdA2i+UyAp+cRI7zeeBh4E9gHf4fXEIGw7cKllvka1As8UQMez4wD9gFtAcGAUuAXsF8bwWGB+O7Fl9bsAF4MjIuoDOwCNgO/AV8ChyVybxXAyYBW4LX20DVoFvvdOarejrjWBnVz8qgPG0eugM/BuN/B6gYNXwfYCmwA/geuDmjZR0xzB3BtrMVGBVMa2V2xgvUCJbPDmB5sA62Ar2D7tWD+ekBzAiW6cCgW7Ng2G3AauAFoFzEuC2I8cdguMXAZVnMUwJwb7BudwbDdIzaNw/YxqLGMSid/lpFzMvFwLQg7qXAOVHDnwy8H6yrdcA44Ogs4n40WH7bg23hcaBUVExLgH7Ar0F/+20HBPsg8C/8sSNtey8dtV++gN8X1wNfx7ANWzC/0wELysoCPxAcNyKWTVLUPnkeMD+I93OgKr4W55sgvinAERHxnQ58jD9mbAZmA01j2E96A1vTWV4Z7jf4P1RPAf8NXk8Fy2ZWJuspbb7a42sZdwTz1zCqv874bW8nflu8J23ZRR9zgfuAJelMaw7wbNS6vRG/r/w3WLeHRPRfEng6WPc7gLlA83Riz+46GRQZX1brKKPfmcL2Cj0AvbK5wnzN2TfB51b4g3PxqH4ik7eLgu/NsjmdCsGO3y6irAw+eekafD8VGADUxdcm3INPXGpFDLPvIBF8P2CnSqefMcGO3xLfFDEwGO+pmcQ7ggOTt03AA/gf+EuAPUCPDIYvERyY/gaODl5lY4kn4qC0GDg36KdScNDZiv8hqgO0Db5/iD/w1Q7Wz27g4mBcRwfjvhX/g1QH/4OZbvKG/2FbAHwRHNSSgliTg26lg+m6oPvRQGI646kU9NMv6KdSUJ42D5PwfxiaAr8AL0UMexW+FrcLcDzQAZ+UDcxkfV0SzOfVwfq5B38gXhnrePGJ0rfAJ0D9ILavguXZO+inejBfKyPGUxW/zW4NlvM/gMbAl8CbEdN/CJ/QtAuGuzTYPi7IZL5uDubj0mC+HgD2AvUj1u93+ORl3zYWNY6y+Ga4afxvWywRMS/fBcviH8BI/B+AtG31GPyP2mP47ase8B7+T0WGyTQ+4TwjmMb5+ATtwYjuadvBLKBB0O+3RPwpxO+DW4CJ/G97X02QAETsl1uA/8M3G9cmi204GK5yMF+3B99fC6ZfOmo9Rydv84AWwXJYgk9IPgnWdxLwM/CfiPjO4n+tE7WA5/CJSsUs9pPeHJi8ZbXf3BWM+2KgJvAM/pg1K5P1lDZf3wXLt06wvH8nSKSAhvht7n78NtgziOX69I65+P1hD9AoonvNYDqnRqzbTcDLwbI5F1+JcHfEMM/g99cLgn5eDqZ7zEGuk0Hsn7xluo4y+p0pbK/QA9ArmyvM1xTcFny2YCO9OKqfyOTtjuB7hRxMaxLwesT3y4IduFQmw8wF/hXxfd9BIvh+wE4VdSA5EUgFqkX18w7wfCbTHcGByduXUf1MA17JZBy9iTgAxxpPxEEpej0Mwv+7PCyi7E18jUOJDOb/tGBcx8W4js7BH6irR5SdEMTcJvieRAY1bhltN1HzsCNqHu4BVkR8/xW4PGq4m4ClmUzrC+DlqLLp7J+8ZTpe/I/XHqBKRPdmwXz0Dr5XD77fGjWeUcCrUWX1g36PxP9R2Q60iOrnaeCDTOZrNXBfOtv36IjvS0inxi2z7TlqXq6OKKsSlDUPvj8AfBI1XIWgn0aZTTNqmAFR63hQsJ1ViyhrHoz3HxExbyQiIcUfM3YCZSKWxaLsbsNBWadgXA8G76dGdEtbNtHJW2QN+cCg7LSo+Tqg1imiu+ETkssiytLbT3pzYPKW1X6zFrgralrfEVvy1jOirGyw3PsF38cAM9LZj1dFbZORx+UpwIsR3x8DkqO2x9+IaInBJ2fTg89l8H/GrojonoivdRx8MOskh+toJYU8edM5bwWImZ2E/8c7FsD5rXQM/l9ghoMdxCRHA53M7JDge098zcSOIJ4yZva4mS01s/8GJ/sn4ZtAcuq0IOalZrY17YX/N3diNse1KOr7GvwPc17Fk5zO8L865zZFfP8D+N45tyuqLC2ub/BJzBIze8vMrjGzSpnEVxtY45xbmVbgnPsJP68nZz17Mfklah72LccgtmOBl6KWz6Nkvr5q42u6Iu37HuN4a+HnfXXEOL7G/+hHi143DYHLosaddgHQifhlVwr4KKqfazKaLzMrh68hir6QaDa5ty5g/+16TfCetv00BFpGxfxb0C3D9WFmXcxsdnBB01Z8E170frzaOfdrxPev8Ms68rzbRW7/cyq/xNcaRk57ftR4Y9qGnXPv4I99/8L/Qfwmo/mJjCfi8x/B++Kosn3HBDM70sxeMrPvzWwTvpbwSHJ2TMtsvzkMX3M3L61jcDz/OsZx79tXguW9mP8tq9qkvw1WCbbR9LwMdDez0maWiK/ZejWqn6Vu/3MKI4+nJwLFI6frnNsbxBm97WdrnUTL5XVUYBULOwDJln74fzO/mu3LyQzAzI51zv2WzjDfB++18bUd2TEFX7PR0cw+Adrgq8vTDME3Kd2GP/9kG75GI7OTkFM5MKGMPIE8gf818e2O6m97NuOPHt6R/Yt0shPP3zHGkF5ZIvgDnpmdCzTBL+srgUfM7MwMfqwsGD49GZVnV2bLMe19ANnfvjITy3gzm/do0esmAXgFn6REW41v0gHfPPlrVPfo5REtvZhya13sN33nnAuOBZHr4338Phntj3TKMLMmwHh8M9vN+FqcC8ngquRcEL0uYtqGzawUfj/ciz9NIxaR68pXmzkXXRZ5TBgJHIVfDivxNXyfkPkxLZZppzetfTHlspwcE97HH78vxreulMefKxkps/mxiLKsppnddRItN9dRgaXkrYAws2L4E9/vxidVkV7Hn9j9QDqDpp3YeRf+gBw93vIug6s0nXM7zexNfI1bRfx5FZ9G9NIcGOWceysYVyn8P7Dvo8cVYT3+vJy06ZfC16AsDIoW4g8ERzvnZmYynrywiyCJipDv8QT/wL8EvjSzB/Dn9nTD18pFW4r/R109rebCzE7A1wAtzeakd3Pg/GcV6x9mtho40Tk3KhuDLsMnqK9FlDXJ5niX4ee9snMurQYqidgS9AXAKc65Fel1NLOl+B+F45xzM2IYH865zWa2Br9fRA7TnOyvi/S2xVgswJ9P+EvUD2JmzsDXqj2YVmBmx6XTX5WoP4mN8Mt6WUQ/dc2sjHMuLUFrgp+XHzOZfqzb8BP4k+LPAaaa2QfOuckxzmOsmgM3OOfeD+I4iojjVSDb+0k059wmM/sdvwxnBtMyfHL6ewyjaAL8FAxXBn/uW9p+sjSYj0jN8c2mWzKIZ4+ZjQD64pO3tzP6XcjACvx6bh4RVyL+XL+x2RhPLGJZR4WekreC4wJ8AvWyi7ps2szGA9eY2WDn3H5NRs65v82sHzDRzN7Hn7PzA3A4/mT504JxZ2Q0vhnveGBs1Pi/By4ys8n4A9q/8U1NmZkB9DWzd/GJ3D1E1Lw55743szHACDO7Ff9jdDj+fImfnHNvZzH+g7ESKGVm5+CTtm35HU9QC9IGmIqvKWmAbz7M6Md/Oj6pG2NmN+ATzf8EccaUdERYCZxtZp8CO51z/41xuEHAf8xsI/ABfn2ehj8X7ZEMhnkGGGVmX+PPv+mCP2H5r2yMdxr+goKR5m/IWRp/5e4esq7ReAyYa2YvAi/hm15qAR2cc1c757aY2RBgSPCj+hn+3KImQKpzblgG430CeMDMfsA3D16GPzm7YRbxRFsJnGdmNfEXJGzKvPd9huIv9JhgZo/h97ET8AndrRn8eH+PT5564v80tMVfnRttO35Z34Jf1i8C7zvnfojopxjwWvCnozK+mfvliGQuPVluw2bWDn9xSwvn3FfmbzT9ipl95ZyLJdmJ1ff45vSv8OdxPY5PSiKtJGf7SbRngDvM7Hv8/n01PglZG8Ow/zKz9fimy/uCGNOSpP8Dvg6W0Vh8Qngr8M8sxvkKcCe+deTcLPrdT/A78wLwqJn9ib/o4GZ8Ddnz2RlXDGJZR4WeznkrOK4EZkYnboGJ+Ht5tUlvwODfaVN8tfho/A/eRHxScEcW0/0M34x0cjBspFvwV7t+jr+Ccm7wOTOP4A/Ik/G1grPxB+lIffBXYz6OP4F3Cv5Kz1+yGPdBcc59gf9BGof/0UtbNvkZzyZ8TcgUfJL9f/ir/qKXfVrMDn8i93p8EjQT/8+9U9AtO24FWuPPkVqYRb+RMbyC/8d+Of5H+HOgP/4AntEwE/DJ2UPBtOriE6+Yxxv8kbgIXxszD9+c8hA+cduRRcyL8OuwOr42+Rv8thnZtHhvEONt+NrPafhmpQznC3gWn8A9jr8w4SL8hSzZvXnsy/garWT8uj0jloGCGsgz8D/AHwVxD8XXIu7MYJj3gpifxp+PdA4+IYi2Et+8+h5+H/4Jv29E+jSY5kz8BU8zyOIYk9U2HJz/OAJ/4vtXwWCPBtMZbhHnkOSCvvgkfT5+Xl/Dz3ekHO0n6RiCbzUZjj92gl9mmW67gbvwx4YF+KuO26clyM65BUBX/La6BL+sHsVflZmh4DzDT/GnCczK3qwAPvF7Az8/KfhTD9o552JJRrMjlnVU6Fn2j+8iIvHJzE7F/3AkOeeiT4yXHApqcbo45+pk0s8I/O0a2udXXIWNmS0A5jjnrs+geyt8clvJOfdnHkx/KTDGOfdQbo9bcpeaTUWkwDKzi/AnwP+Ar0V7El+LFl2bKxJXgvMK2+Jru4rha5VPDd7zO5Yj8U3l1fGnEUicU/ImIgXZofjz147F36hzFnBzDpqMRfJbKnAFvsk6AX/e23nOufRuOZTX/sBf2HZ1XtToSe5Ts6mIiIhIAaILFkREREQKECVvIiIiIgVIkTrnrWLFiq569ephhyEiIiKSpfnz5//pnDvgEYlFKnmrXr06yclhnAsqIiIikj1mlu79RNVsKiIiIlKAKHkTERERKUCUvImIiIgUIEXqnDcRERHxdu/ezapVq9ixI5bHqUpeKlWqFFWrVqV48eIx9a/kTUREpAhatWoVhx56KNWrV8fMwg6nyHLOsWHDBlatWsXxxx8f0zBqNhURESmCduzYwRFHHKHELWRmxhFHHJGtGtDQkzcz62pm35pZqpklRZSfY2bzzWxx8H5WRLdZZrbczFKC15HhRC8iIlJwKXGLD9ldD6Enb8ASoDPwWVT5n0AH51xdoBfwelT3ns65+sFrXT7Emal3Fq7mjEdncPxd73PGozN4Z+HqsEMSERGJe5MmTcLM+O677/aVrVy5kjp16uz7Pm/ePFq2bEnNmjWpVasW/fr1Y9u2bQc13VmzZtG+ffuDGkdYQk/enHPLnHPL0ylf6JxbE3z9FihlZiXzN7rYvLNwNXe/vZjVG7fjgNUbt3P324uVwImIiGRh3LhxNG/enPHjx6fb/Y8//qBr16489thjLF++nGXLltGuXTu2bNmSz5HGj9CTtxhdDCx0zu2MKBseNJnea5nUN5pZfzNLNrPk9evX50lwT0xdzvbde2ma8C3NEpYAsH33Xp6YekBOKiIiUiDlRQvT1q1bmTNnDq+++mqGydvQoUPp1asXTZs2BXwTY5cuXTjqqKP2669x48Z8++23+763atWK+fPnM2/ePJo1a0aDBg1o1qwZy5cf+Ns8aNAghgwZsu97nTp1WLlyJQCjR4+mUaNG1K9fn6uvvpq9e/eyd+9eevfuTZ06dahbty5PPfXUwS6KbMmX5M3MppvZknReHWMY9hTgMeDqiOKeQXNqi+B1eUbDO+eGOeeSnHNJlSod8HiwXLFm43YArk2czNgSD/N4sZcox9Z95SIiIgVZXrUwvfPOO7Rr144aNWpw+OGHs2DBggP6WbJkCQ0bNsxyXN27d+eNN94AYO3ataxZs4aGDRtSq1YtPvvsMxYuXMgDDzzAP//5z5jjW7ZsGRMmTGDOnDmkpKSQmJjImDFjSElJYfXq1SxZsoTFixfTp0+f2Gc6F+RL8uaca+Ocq5POa3Jmw5lZVWAScIVz7seI8a0O3rcAY4FGeRl/ViqXLw1Av9238fyeC+mc+DmflLydSw9NCTMsERGRXJHWwhQpN1qYxo0bR/fu3QGffI0bNy7H47rkkkuYOHEiAG+88QZdu3YFYNOmTXTt2pU6depw880371c7l5VPPvmE+fPnc/rpp1O/fn0++eQTfvrpJ0444QR++uknrr/+ej766CPKlSuX47hzIm6bTc2sPPA+cLdzbk5EeTEzqxh8Lg60x1/0EJrb29akdPFEdlKCx/d0p+OuwaynAg/tfhzG94Qtv4cZnoiIyEHJqCXpYFqYNmzYwIwZM+jXrx/Vq1fniSeeYMKECTjn9uvvlFNOYf78+VmOr0qVKhxxxBEsWrSICRMm7EsK7733Xlq3bs2SJUt477330r0lR7FixUhNTd33Pa0f5xy9evUiJSWFlJQUli9fzqBBg6hQoQLffPMNrVq1YujQofTr1y/HyyEnQk/ezOwiM1sFNAXeN7OpQaeBwEnAvVG3BCkJTDWzRUAKsBp4OYzY03RqUIVHOtelSvnSGLDxsNr8cOG70GYQrJgOzzWC+SMhaoMUEREpCNJamGItj8Wbb77JFVdcwS+//MLKlSv57bffOP7445k9e/Z+/Q0cOJCRI0fy1Vdf7SsbPXo0v/9+YMVI9+7defzxx9m0aRN169YFfM1blSpVABgxYkS6sVSvXn1fk+2CBQv4+eefATj77LN58803WbfO39Tir7/+4pdffuHPP/8kNTWViy++mAcffDDd5t68FPoTFpxzk/BNo9Hlg4HBGQyWdeN3PuvUoAqdGlSJKr0Zal8I794A790AiydCh2fgiBNDiVFERCQnbm9bk7vfXrxf02np4onc3rZmjsc5btw47rrrrv3KLr74YsaOHcudd965r+yoo45i/Pjx3Hbbbaxbt46EhARatmxJ586dDxhnly5duPHGG7n33nv3ld1xxx306tWLJ598krPOOuuAYdKmO2rUKOrXr8/pp59OjRo1ADj55JMZPHgw5557LqmpqRQvXpyhQ4dSunRp+vTps6+27pFHHsnxcsgJi66eLMySkpJccnJy/k84NRUWjIRp98He3dD6n9DkWkgMPXcWEZEiatmyZdSuXTvm/t9ZuJonpi5nzcbtVC5fmtvb1kyn0kJyKr31YWbznXNJ0f0qe8gPCQmQ1AdqtIX3b4Vp98KSt6Djc3B03bCjExERyVL6LUwShtDPeStSylWG7mOhy3DYvBqGtYJPHoDdsT/PTERERIo2JW/5zQzqdIbr5kHdS+Dz/4MXm8MvX4QdmYiIiBQASt7CcsjhcNELcNnbsHcnDD8PptwCOzaHHZmIiIjEMSVvYTvpbLjmS38BQ/Jr8HwT+H5q1sOJiIhIkaTkLR6ULAvtHoF+06FkORh7Cbx5Jfz9Z9iRiYiISJxR8hZPqibB1Z9Bq3/C0snw3OnwzQTd3FdERAqlsmXLZtnP008/zbZt2/I8lhEjRjBw4MBM+5k1axZffPG/c9RffPFFRo0aldehHUDJW7wpVgJa3QkDPvc3853UH8Z0hY2/hh2ZiIhIvstJ8rZ3796se8qB6ORtwIABXHHFFXkyrcwoeYtXR9aGvlOh3WP+StShTeCrl/wNf0VERAqRWbNm0apVK7p06UKtWrXo2bMnzjmeffZZ1qxZQ+vWrWndujUAH3/8MU2bNuW0006ja9eubN26FfCPuHrggQdo3rw5EydOpFWrVtx00000a9aMOnXqMG/ePMA/4qpTp07Uq1ePJk2asGjRogPiee+992jcuDENGjSgTZs2/PHHH6xcuZIXX3yRp556ivr16/P5558zaNAghgwZAkBKSgpNmjShXr16XHTRRfz3v/8FoFWrVtx55500atSIGjVq8Pnnnx/08lLyFs8SEqHJALj2S6jWBD68A15rC+u+CzsyERGRXLVw4UKefvppli5dyk8//cScOXO44YYbqFy5MjNnzmTmzJn8+eefDB48mOnTp7NgwQKSkpJ48skn942jVKlSzJ49e99D6f/++2+++OILnn/+efr27QvAv//9bxo0aMCiRYt4+OGH0605a968OXPnzmXhwoX7npdavXp1BgwYwM0330xKSgotWrTYb5grrriCxx57jEWLFlG3bl3uv//+fd327NnDvHnzePrpp/crzyk9YaEgqHAcXPYWLJoAH90FL7WAlrfDGTf5ZlYREZGD8eFd8Pvi3B3n0XXhvEdj7r1Ro0ZUrVoVgPr167Ny5UqaN2++Xz9z585l6dKlnHHGGQDs2rWLpk2b7uverVu3/frv0aMHAC1btmTz5s1s3LiR2bNn89ZbbwFw1llnsWHDBjZt2rTfcKtWraJbt26sXbuWXbt2cfzxx2ca+6ZNm9i4cSNnnnkmAL169aJr1677uqc9h7Vhw4asXLkypuWRGdW8FRRmcGp3uO5rqNUeZj4Ew86EVfPDjkxEROSglSxZct/nxMRE9uzZc0A/zjnOOeccUlJSSElJYenSpbz66qv7upcpU2a//s3sgO/pPdM9ur/rr7+egQMHsnjxYl566SV27Di4JyGlzVtG85VdqnkraMpWgq7Dod4l/qa+r7aBxtfAWfdAiTJZDy8iIhItGzVk+e3QQw9ly5YtVKxYkSZNmnDdddexYsUKTjrpJLZt28aqVauoUaNGusNOmDCB1q1bM3v2bA477DAOO+wwWrZsyZgxY7j33nuZNWsWFStWpFy5cvsNt2nTJqpU8c9xHTly5H6xbN584M30DzvsMCpUqMDnn39OixYteP311/fVwuUFJW8FVc3z4LhmMH0QzB0K302BDs/Aia3DjkxERCTX9O/fn/POO49jjjmGmTNnMmLECHr06MHOnTsBGDx4cIbJW4UKFWjWrBmbN2/mtddeA2DQoEH06dOHevXqccghh+yXnKUZNGgQXbt2pUqVKjRp0oSff/4ZgA4dOtClSxcmT57Mf/7zn/2GGTlyJAMGDGDbtm2ccMIJDB8+PDcXw34sverDwiopKcklJyeHHUbuWzkH3rsBNqyA+pdB28FQukLYUYmISBxbtmwZtWvXDjuMPNOqVSuGDBlCUlJS2KHEJL31YWbznXMHzIDOeSsMqp8BA+ZA81vgm3HwXCP49h3d3FdERKQQUvJWWBQvBW3+Df1nQbljYGIvmHAZbF4bdmQiIiL5btasWQWm1i27lLwVNsfUg34zoM39sGI6DG0M80eoFk5ERKSQUPJWGCUWg+Y3wTVf+GTuvRthZAfY8GPYkYmISBwpSue9x7Psrgclb4XZESfCFe/6q1DXfgMvNIPZT8Peg7/HjIiIFGylSpViw4YNSuBC5pxjw4YNlCpVKuZhdLVpUbF5LXxwm7+lyDGnwoXP+Vo5EREpknbv3s2qVasO+ga0cvBKlSpF1apVKV68+H7lGV1tquStKHEOlk72Sdy2v+CMG+HMO/3FDiIiIhJX4vZWIWbW1cy+NbNUM0uKKK9uZtvNLCV4vRjRraGZLTazFWb2rEU/10LSZwandILr5vlHbc1+El48A375IuzIREREJEahJ2/AEqAz8Fk63X50ztUPXgMiyl8A+gP/CF7t8j7MQuSQw6HT83D5JNi7C4afB1Nuhh0HPvJDRERE4kvoyZtzbplzbnms/ZvZMUA559yXzrf5jgI65VmAhdmJZ8G1c6HJdf52Is83geUfhR2ViIiIZCL05C0Lx5vZQjP71MxaBGVVgFUR/awKyiQnSpSBdg/DldOh1GEwrhu82Re2rg87MhEREUlHviRvZjbdzJak8+qYyWBrgWrOuQbALcBYMysHpHd+W4ZXXZhZfzNLNrPk9euVkGSoakPo/ym0+icsfReGng7fjNfNfUVEROJMsfyYiHOuTQ6G2QnsDD7PN7MfgRr4mraqEb1WBdZkMp5hwDDwV5tmN44ipVgJaHUnnNwR3r0eJl0Ni96ADk9D+WphRyciIiLEcbOpmVUys8Tg8wn4CxN+cs6tBbaYWZPgKtMrgMkhhlr4HFkL+n4E5z0Ov86FoU1g7ouQujfsyERERIq80JM3M7vIzFYBTYH3zWxq0KklsMjMvgHeBAY45/4Kul0DvAKsAH4EPsznsAu/hERofDVcNxeOawof3QmvtYV134UdmYiISJGmm/RK1pzzzacf3QW7tkKL2xXQDFsAACAASURBVKD5zb6ZVURERPJE3N6kVwoAMzi1m7+5b+0LYdbDMOxMWKVEWEREJL8peZPYla0EXV6FHhNgxyZ4pQ18dDfs+jvsyERERIoMJW+SfTXb+Zv7JvWFuc/7m/v+OCPsqERERIoEJW+SM6XKQfsnoc+HkFgCXr8IJl3jH3gvIiIieUbJmxyc45rBgDnQ4lZYNAGGNoJvJ+nmviIiInlEyZscvOKl4Oz7oP8sKFcZJvaG8T1hc4b3ThYREZEcUvImueeYetBvBpzzAPz4CQxtDMnDITU17MhEREQKDSVvkrsSi8EZN8I1X8Axp8KUm2DUhbDhx7AjExERKRSUvEneOOJE6PUedHgW1i6CF5rB7Kdh756wIxMRESnQlLxJ3jGDhr3guq/gpDYw/d/wylk+mRMREZEcUfImea/cMdB9DFwyCjavhWGtYPog2L097MhEREQKHCVvkn9O7ggD50H9HjD7KXjhDFg5J+yoREREChQlb5K/SleAjkPh8ncgdQ+MOB/eu8k/bktERESypORNwnFia7j2S2g6EBaMhKFN4LsPwo5KREQk7il5k/CUKANtH4Irp/saufE9YGIf2Lo+7MhERETiVrGwAxChakP/dIY5z8Bnj8NPM6HtI7yT2oInPv6eNRu3U7l8aW5vW5NODaqEHa2IiEioVPMm8aFYCTjzdhgwGyrWgHcGcMTkS7FNv+KA1Ru3c/fbi3ln4eqwIxUREQmVkjeJL5VqQp+PeLLYVTRgOVNL3EHbhK8B2L57L09MXR5ygCIiIuFS8ibxJyGB/2xtzbk7H2e5O5bniz9N54TPAFizUfeGExGRok3Jm8SlyuVLs4aKXLbrn3yZejJPlniRyxM/pnL50mGHJiIiEiolbxKXbm9bk9LFE9lGKa7cfTsf723Ig8VH8PLxs8C5sMMTEREJjZI3iUudGlThkc51qVK+NLsoweBD7uK3qh04edkz/hmpSuBERKSI0q1CJG51alBl/1uDpLaBD27ztxTZuQXO/z9I0P8PEREpWkL/5TOzrmb2rZmlmllSRHlPM0uJeKWaWf2g2ywzWx7R7cjw5kDyTUICXPB/cMZNkPwaTOoPe3eHHZWIiEi+ioeatyVAZ+ClyELn3BhgDICZ1QUmO+dSInrp6ZxLzrcoJT6YwTn3Q6ly8MkDsOtv6DIcipcKOzIREZF8EXrNm3NumXMuq5t39QDG5Uc8UkC0uBXOHwLLP4Cxl8DOrWFHJCIiki9CT95i1I0Dk7fhQZPpvWZmYQQlIWt0FVz0EqycDa93gu3/DTsiERGRPJcvyZuZTTezJem8OsYwbGNgm3NuSURxT+dcXaBF8Lo8k+H7m1mymSWvX68Hnhc6p3aHS0bC2m9gRHvYui7siERERPJUviRvzrk2zrk66bwmxzB4d6Jq3Zxzq4P3LcBYoFEm0x7mnEtyziVVqlTpYGZD4lXtDnDpBPjrJxh+Hmz8LeyIRERE8kxcN5uaWQLQFRgfUVbMzCoGn4sD7fEXPUhRduJZcPkk2LoeXmsHG34MOyIREZE8EXryZmYXmdkqoCnwvplNjejcEljlnPspoqwkMNXMFgEpwGrg5XwLWOJXtSbQ+z3Ys90ncL8rpxcRkcLHXBG6U31SUpJLTtbdRQq99d/DqI6w+2/o+RYce3rYEYmIiGSbmc13ziVFl4de8yaS6yrVgL4fQenDfRL306dhRyQiIpJrlLxJ4VThOJ/Ala8GY7rC8g/DjkhERCRXKHmTwuvQo6HPB3DUKTC+JyyaGHZEIiIiB03JmxRuhxwOvd6Fak3h7asgeXjYEYmIiBwUJW9S+JU8FC57E/5xLky5CeY8E3ZEIiIiOabkTYqG4qWh22g4pTNMuw8+eRCK0JXWIiJSeBQLOwCRfFOsBFz8CpQsC58PgZ1boN2jkKD/MCIiUnAoeZOiJSEROjwLJcvBl8/Brq3+e6J2BRERKRj0iyVFjxmcO9gncLMe9jVwF78CxUqGHZmIiEiW1F4kRZMZtLoT2j4Cy96FcT1g17awoxIREcmSkjcp2ppeCxc+Bz/NhNGdYcemsCMSERHJlJI3kdMuh4tfhVVfw8gO8PefYUckIiKSISVvIgB1OkP3cbB+OQw/HzavCTsiERGRdCl5E0lT41y47C2fuL3WDv76OeyIREREDqDkTSRS9eb+cVo7N/sEbt2ysCMSERHZj5I3kWhVToPeH/jPw8+H1QvCjUdERCSCkjeR9Bx1MvT90D+NYeSFsHJO2BGJiIgASt5EMnb4CdB3KpQ7xt9G5IdpYUckIiKi5E0kU+UqQ58PoWINfyPfbyeFHZGIiBRxSt5EslKmIvSeAlUawpt9YcHrYUckIiJFmJI3kViUOgwufxtOaAXvDoQvnw87IhERKaKUvInEqkQZ6DEeaneAqXfDrMfAubCjEhGRIkbJm0h2FCsJXUbAqZfCrIfh438pgRMRkXxVLOwARAqcxGLQcai/jciXz8HOLdD+KUhIDDsyEREpAuKi5s3MnjCz78xskZlNMrPyEd3uNrMVZrbczNpGlLcLylaY2V3hRC5FVkICnPc4tLgNFoyEt/rB3t1hRyUiIkVAXCRvwDSgjnOuHvA9cDeAmZ0MdAdOAdoBz5tZopklAkOB84CTgR5BvyL5xwzOvhfa3A/fvg3je8Lu7WFHJSIihVxcJG/OuY+dc3uCr3OBqsHnjsB459xO59zPwAqgUfBa4Zz7yTm3Cxgf9CuS/5rfBBc8CT98DGO6+mZUERGRPBIXyVuUvsCHwecqwG8R3VYFZRmVH8DM+ptZspklr1+/Pg/CFQFOvxI6vwy/fOEfp7Xtr7AjEhGRQirm5M3MiptZCzPrFnwvY2ZlsjH8dDNbks6rY0Q/9wB7gDFpRemMymVSfmChc8Occ0nOuaRKlSrFGq5I9tXrCt1Gwx/fwogLYMvvYUckIiKFUExXm5pZXeBdYCe+SXMCcCbQC+gWyzicc22ymEYvoD1wtnP77r2wCjg2oreqwJrgc0blIuGpdT70fAPGXQqvtYMrJkOF48KOSkRECpFYa95eAO5zztUC0i6p+xRonhtBmFk74E7gQufctohO7wLdzaykmR0P/AOYB3wN/MPMjjezEviLGt7NjVhEDtoJrXzStv0vGH4erP8+7IhERKQQiTV5OwUYHXx2AM65v4HSuRTHc8ChwDQzSzGzF4NpfAu8ASwFPgKuc87tDS5uGAhMBZYBbwT9isSHY0+H3u/D3l0+gVv7TdgRiYhIIWEuhrvDm9lC4CrnXLKZ/eWcO9zMGgHPOeca5XmUuSQpKcklJyeHHYYUJX+ugFEd/RWoPSdCtcZhRyQiIgWEmc13ziVFl8da83Yv8L6Z3Q+UMLO7gYnAv3IxRpHCp+JJ0PcjKFMRXu8EP84IOyIRESngYkrenHNT8DfErYQ/1+04oLNz7uM8jE2kcCh/rE/gDj8BxnaDZe+FHZGIiBRgMd8qxDm3wDl3rXPuAufcAOfc/LwMTKRQKXsk9J4Cx5wKb/SCb8aHHZGIiBRQsd4q5IGMujnn7su9cEQKsdIV4PJ3YHwPmHS1Pw+u0VVhRyUiIgVMTMkb+99TDeBo/H3eJuVuOCKFXMmycOlEmNgbPrjNJ3Atbgk7KhERKUBiSt6cc32iy4J7s/XI9YhECrvipaDb6/DONfDJ/bBjE7QZ5B90LyIikoVYa97S8zH+SQsikl2JxeGil6BEWZjztK+BO38IJMTj44ZFRCSexHrO2wlRRYcAl7L/w+FFJDsSEqH9U1CqHMx5BnZthY7PQ+LB/KcSEZHCLtZfiRXs/0D4bcBC/LNNRSSnzKDN/VCyHMx4EHZuhS6v+aZVERGRdMR6zpvackTyihm0vA1KHgof3gHjukH3sVCiTNiRiYhIHFJSJhIvGl8NnV6Anz+DUZ1g+8awIxIRkTiUYc2bmf1G8BD6zDjnquVqRCJFWf1LfY3bm1fCiPZw+SQoWynsqEREJI5k1mx6Wb5FISL/c3JHuLQMjL8MhreDKybDYVXDjkpEROJEhsmbc+7T/AxERCKc1MbXuo29BF4LErgjTgw7KhERiQMx35PAzOoDLYCK/O+qUz0eSySvHNcUer0HozsHCdw7cNQpYUclIiIhi+mCBTPrD8wBzgLuBOoCtwIn5V1oIkLl+tDnQ39PuOHnw6rksCMSEZGQxXq16R1AO+fcRcD24L0LsDvPIhMRr1JN6PsRlC4Pozr6q1FFRKTIijV5O9I593nwOdXMEpxzHwId8iguEYlUoTr0+chfuDC6Cyz/KOyIREQkJLEmb6vMrHrw+Xugo5m1AHblRVAiko5yx0DvD+Cok2FCT1j8ZtgRiYhICGJN3h4HagefHwBGAzOA+/MiKBHJQJkj4Ip34djG8FY/SB4edkQiIpLPYn081oiIzx+aWQWghHNua14FJiIZKFUOer4Jb1wBU27yD7Rvdn3YUYmISD6J9WrTp83s9LTvzrldStxEQlTiEP/805M7wcf/ghkPgcvygSgiIlIIxHqfNwMmm9nfwFhgrHNued6FJSJZKlYCurwG75WFzx6HnZuh7SOQoEcWi4gUZjEd5Z1zNwJVgWuBY4G5ZjbfzG452ADM7Akz+87MFpnZJDMrH5SfE0xjcfB+VsQws8xsuZmlBK8jDzYOkQIpIRE6/AeaXAtfvQjvXg+pe8OOSkRE8lDMf9Gdc6nOuWnOub5AHWAD8EQuxDANqOOcq4e/kvXuoPxPoINzri7QC3g9ariezrn6wWtdLsQhUjAlJEDbh+HMuyBlNLzZB/boQnARkcIq5uTNzMqa2WVm9j4+ydqDT6oOinPuY+fcnuDrXHwNH865hc65NUH5t0ApMyt5sNMTKZTMoPXdcO5DsHQyjO8Bu7aFHZWIiOSBWC9YmAj8DvQHpgDHOefOd86NzuV4+gIfplN+MbDQObczomx40GR6r5lZOsOIFD3NBkKHZ2HFJzD6YtixKeyIREQkl8V6wUIycKtz7tecTMTMpgNHp9PpHufc5KCfe/C1eWOihj0FeAw4N6K4p3NutZkdCrwFXA6MymDa/fFJJ9WqVctJ+CIFS8NeULIsvN0fRl4Il73t7w8nIiKFgrk4uL2AmfUCBgBnO+e2RZRXxd8MuI9zbk4Gw/YGkpxzA7OaTlJSkktO1oO9pYj4fqq/F1yF6nD5JChXOeyIREQkG8xsvnMuKbo89HsKmFk74E7gwqjErTzwPnB3ZOJmZsXMrGLwuTjQHliSv1GLFAA12vqb+W5aBa+1g79+DjsiERHJBaHXvJnZCqAk/upVgLnOuQFm9i/8lac/RPR+LvA38BlQHEgEpgO3OOeyvD+Cat6kSFo1H0Z3Zjsl6Jf6L77YXJHK5Utze9uadGpQJezoREQkAxnVvIWevOUnJW9SVH3y6SzqzehFInu4YtddLHEnULp4Io90rqsETkQkTuWo2dTMqplZnwy69Q7OSROROHffl6l02XUf2yjF+BKDaZMwn+279/LEVD0oRUSkoMnqnLf7gFIZdCsZdBeROLdm43Z+cUdz8c5B/OgqM6z4k1yT+C5rNupecCIiBU1WydtZQEb3chsDnJO74YhIXqhcvjQAf3A4l+y6jympTbiz+HhePOQl2L0j5OhERCQ7skreKuEvEEjPdqBi7oYjInnh9rY1KV08EYCdlOCG3QN5OrUbbVM/gxHnw5bfQ45QRERilVXythaon0G3U/FPXRCRONepQRUe6VyXKuVLY0CV8odQ/aJ/Q7fRsO47GNYa1iwMO0wREYlBplebmtkg/H3ULox4zihmVhmYBHzgnLs/r4PMLbraVCQdvy+GcT3g7z+h0/NQp3PYEYmICDm/Se9DwBrgBzObaWZjzWwm/t5ra4PuIlKQHV0XrpoJx5wKb/aBGQ9BamrYUYmISAYyTd6cc7udcxcCHYG5wNbg/ULnXCfn3J58iFFE8lrZStDrXah/GXz2OEy8AnZldLqriIiEKaYH0zvnpuOfZCAihVWxktDxOTiyNky7F15tCz3GQfljw45MREQiZJq8mdkDWY3AOad7vYkUFmbQbCBUqglv9oWXW0O3MVCtcdiRiYhIIKtz3o6N4SUihc0/zoF+06HkoTCyPSwcE3ZEIiISyLTmzTmX7qOx0phZVsmfiBRUlWpCv09gYm+YfC2sWwrnPAAJiWFHJiJSpOUo+TKzumb2BLAql+MRkXhyyOFw2Vtw+lXw5XMwrjvs2BR2VCIiRVrMyZuZVTKzG81sAZACNAJuzLPIRCQ+JBaHC4bABU/CjzPglXNgw49hRyUiUmRlmryZWXEzu9jM3gNWA1fjb867EejqnJuYDzGKSDw4/Uq4fBL8vQ5eORt+/izsiEREiqSsat7+AF4ClgNNnHMnO+ceBHbleWQiEn+ObwlXzYAyR8LrF8HXr4QdkYhIkZNV8rYIKA80Bk43swp5H5KIxLXDT/BXop54Frx/q3/t3R12VCIiRUZWT1hoBZwIfAzcBvweNKGWAYrneXQiEp9KlYMe46HZ9b72bXRn2PZX2FGJiBQJWV6w4Jz7xTn3oHPuH8DZ+GeapgLfmNnjeR2giMSphEQ4dzB0egF+nQsvnwXrl4cdlYhIoZetW4U452Y75/oDRwPXA3XzJCoRKTjqXwq9psCurfBKG/hhWtgRiYgUajm6z5tzbodzbpxz7rzcDkhECqBqjeGqmVDhOBh7CXzxHDgXdlQiIoWSnpAgIrmj/LHQdyrUag8f3wOTr4M9O8OOSkSk0FHyJiK5p0QZ6DoSzrwTUsbAyA6wdX3YUYmIFCqhJ29m9oSZfWdmi8xskpmVD8qrm9l2M0sJXi9GDNPQzBab2Qoze9bMLLw5EJH9JCRA639Cl+GwdhG83Bp+Xxx2VCIihUboyRswDajjnKsHfA/cHdHtR+dc/eA1IKL8BaA/8I/g1S7fohWR2NTpDH0/hNS98Oq5sOy9sCMSESkUQk/enHMfO+f2BF/nAlUz69/MjgHKOee+dM45YBTQKY/DFJGcqNwA+s+EI0+GCZfBp0/oQgYRkYMUevIWpS/wYcT3481soZl9amYtgrIqwKqIflYFZSISjw49Gnq/D/W6wczB8NaVsGtb2FGJiBRYxfJjImY2HX9vuGj3OOcmB/3cA+wBxgTd1gLVnHMbzKwh8I6ZnQKkd35bhn/lzaw/vomVatWq5XwmRCTnipeCi16CI2vD9Pthw4/QYxyUqxx2ZCIiBU6+JG/OuTaZdTezXkB74OygKRTn3E5gZ/B5vpn9CNTA17RFNq1WBdZkMu1hwDCApKQktdeIhMUMmt8MlWrBW/1gWGvoPhaqNgw7MhGRAiX0ZlMzawfcCVzonNsWUV7JzBKDzyfgL0z4yTm3FthiZk2Cq0yvACaHELqI5ETN8+DKaVCsBAw/DxZNDDsiEZECJfTkDXgOOBSYFnVLkJbAIjP7BngTGOCcS3vy9TXAK8AK4Ef2P09OROLdUSf7JzJUTYK3+/mm1NTUsKMSESkQzBWhK7+SkpJccnJy2GGISJo9u+CD22DBSKh5AXR+CUoeGnZUIiJxwczmO+eSosvjoeZNRIqqYiWgwzNw3uPw/Yfwalv47y9hRyUiEteUvIlIuMyg8dVw2VuweZV/IsMvX4QdlYhI3FLyJiLx4cSzoN8MKF0BRl4IC0aFHZGISFxS8iYi8aPiSdBvOhzfAt69Hj66G/buyXo4EZEiRMmbiMSX0hXg0onQ+BqY+zyMvQS2bww7KhGRuKHkTUTiT2IxOO9R6PAs/PwZvNIG/lwRdlQiInFByZuIxK+GveCKybD9L3jlLPhxRtgRiYiETsmbiMS36mfAVTOgXBUY3QW+GgZF6P6UIiLRlLyJSPyrUB2u/BhqtIUPb4cpN/kb/IqIFEFK3kSkYCh5KHQbA81vgfkj4PWL4O8NYUclIpLvlLyJSMGRkABt/g2dX4ZVX/sb+q5bFnZUIiL5SsmbiBQ89S6BPh/Cnh3+StTlH4YdkYhIvlHyJiIFU9WGcNVMOOIkGNcDZj+tCxlEpEhQ8iYiBddhVXwN3CkXwfR/w6QBsHtH2FGJiOSpYmEHICJyUEocAl1egyNPhpmDYcMK6D4WDj0q7MhERPKEat5EpOAzgzNvh0teh3VL/YUMa1LCjkpEJE8oeRORwuPkC6HvVMDgtXbw7aSwIxIRyXVK3kSkcDmmHvSf6d8n9oaZj0BqathRiYjkGiVvIlL4lD0Ser0H9XvCp4/Cm71h199hRyUikiuUvIlI4VSsJHQcCucOhmXv+WbUTavCjkpE5KApeRORwssMml0PPSbAf1fCsNbw27ywoxIROShK3kSk8KtxLlw5DUqUgREXQMq4sCMSEckxJW8iUjQcWQuumgHVmsA7A2DafZC6N+yoRESyTcmbiBQdhxwOl70Np/eDOc/4x2rt2Bx2VCIi2RJ68mZmT5jZd2a2yMwmmVn5oLynmaVEvFLNrH7QbZaZLY/odmS4cyEiBUZicbjg/+D8IbBiOrx6Dvz1c9hRiYjELPTkDZgG1HHO1QO+B+4GcM6Ncc7Vd87VBy4HVjrnIm+Z3jOtu3NuXf6HLSIFWqOr4PJJsOV3ePks+PnzsCMSEYlJ6Mmbc+5j59ye4OtcoGo6vfUAdIaxiOSuE87058GVqQSvd4Lk4WFHJCKSpdCTtyh9gQ/TKe/Ggcnb8KDJ9F4zs4xGaGb9zSzZzJLXr1+fm7GKSGFwxInQbxqc0Bqm3AQf3A5792Q9nIhISPIleTOz6Wa2JJ1Xx4h+7gH2AGOihm0MbHPOLYko7umcqwu0CF6XZzRt59ww51yScy6pUqVKuTpfIlJIlDoMLp0ATQfCvGEw5mLY/t+woxIRSVex/JiIc65NZt3NrBfQHjjbOeeiOncnqtbNObc6eN9iZmOBRsCo3ItYRIqchERo+xAcebKvgXv5bOgxHirVCDsyEZH9hN5sambtgDuBC51z26K6JQBdgfERZcXMrGLwuTg+6YuslRMRybkGPaHXFNi5GV5pAz9MDzsiEZH9hJ68Ac8BhwLTgnPYXozo1hJY5Zz7KaKsJDDVzBYBKcBq4OV8i1ZECr9qjeGqmVC+GoztCl8+Dwc0CoiIhCNfmk0z45w7KZNus4AmUWV/Aw3zOCwRKerKHwt9P/JPY5h6N6z7Fi540j/wXkQkRPFQ8yYiEp9KloWuo6DlHbBwNIzqCFt11bqIhEvJm4hIZhIS4Kx7oMtrsGahv6Hv7zrNVkTCE3qzqYhIgVDnYqhwPIy/FF49FzoP450dDXhi6nLWbNxO5fKlub1tTTo1qBJ2pCJSyKnmTUQkVlVO8xcyHFkLJvTk50kPsHrjNhyweuN27n57Me8sXB12lCJSyCl5ExHJjnLHQO/3mZrQkpsTxvNM8aGUZgcA23fv5Ympy0MOUEQKOyVvIiLZVbw0A7ZdzWO7u9Mh4Us+L3kTAxLfpQzbWbNxe9jRiUghp+RNRCQHKpc/hBf2XkiXXf/m29Tq3FV8PLNL3sjdZafAjk1hhycihZiSNxGRHLi9bU1KF09kgatBr9130XHnA6RQg/57xsJTdWHmw7Dtr7DDFJFCSMmbiEgOdGpQhUc616VK+dIY8OdhddnUaTRc/Rmc0BI+fQyergfT74e//ww7XBEpROzA58AXXklJSS45OTnsMESkKPjjW/hsCHw7CYqXhqS+0OwGOPSosCMTkQLCzOY755Kiy1XzJiKSF446BboOh+u+gtodYO7z8Ew9+PAu2Lwm7OhEpABT8iYikpcq1YTOw2BgMtTpAvOGwTOnwpRbYONvYUcnIgWQkjcRkfxwxInQaSjcsADqXwoLRsGzDeDd6+Gvn8OOTkQKECVvIiL5qUJ16PAM3LAQGvaGbybAfxrCpGvgzxVhRyciBYCSNxGRMJQ/Fi4YAjd+A42v9hc2DD0d3uoH674LOzoRiWNK3kREwlTuGGj3CNy0CJoOhO8+gOebwBu94PclYUcnInFIyZuISDwoeySc+yDctBha3AIrPoEXz4DxPWFNStjRiUgcUfImIhJPyhwBZ98HNy+GM++ClZ/DsDNhzCWwSvepFBElbyIi8al0BWh9t6+JO+tfsGoevHI2vH4R/PJl2NGJSIiUvImIxLNSh0HL230S1+Z+WLsIhreDEe3h58+gCD0lR0Q8JW8iIgVByUOh+U0+iWv7MPz5PYzsAMPP8+fHKYkTKTKUvImIFCQlDoGm1/lbjJz3BGz8FUZ3hlfawPdTlcSJFAFK3kRECqLipaFxf3+z3/ZPwdZ1MPYSf3HDsimQmhp2hCKSR+IieTOzB81skZmlmNnHZlY5KDcze9bMVgTdT4sYppeZ/RC8eoUX/f+3d+dBUtZ3Hsffnzk4RHQ8qMglrhZSZSkCjmgWo2iMSDSi4q5LeSRGRV3PTcV7jRuNbhJSW8Qj3hpjVAyCR9QoWK5RIR4DEtB4xHgBagFLEFAYZnq++0c/QDtpZgbo6ae75/Oq6prn+j3P9/kxRX/mOc3MUlTTHeq/n33t1ribYe1KeOgkuPWg7IN/HeLMKk5JhDdgUkQMjYhhwBPAj5LpY4HByWcicAuApB2Bq4EDgJHA1ZJ2KHrVZmaloroWhp8M5zXAcbdDZh1M/V72gb/zp0JLJu0KzaxASiK8RcTKnNFewPqLNsYBv4msl4E6SX2BMcDMiFgeEX8HZgJHFrVoM7NSVF0D+54I574CJ9wNqoLpZ8BN+8O8ByDTlHaFZraVSiK8AUi6TtJC4CQ2HnnrDyzMWWxRMm1T083MDKCqGvYeD+fMhn+9D2q3gUfPgRv3gzn3QvO6tCs0sy1UtPAm6VlJb+T5jAOIiCsjYiBwP3De+mZ5VhVtTM+33YmSGiQ1LF26tBC7YmZWPqqqYK9j4OwXYcIU2GYn+P0FcOMIePUOaFqbdoVmtpmKFt4i4vCI2DvP57FWiz4AjE+GFwEDc+YNAD5pY3q+7d4eX0BkBAAADMpJREFUEfURUd+nT5/C7IyZWbmRYMhYOPM5OGka9O4LT/0QbhgGL98CTWvSrtDMOqgkTptKGpwzegzwdjL8OHBqctfpgcDnEfEp8AxwhKQdkhsVjkimmZlZWyQYfDicPgNOfQx23AOevgwmD4VZN0Dj6rQrNLN21KRdQOKnkoYALcBHwNnJ9KeAbwPvAV8CpwFExHJJ1wKvJctdExHLi1uymVkZk2D30dnPh7PghZ/DzKtg1uTsQ4D3PxN6bJdujWaWl6ILPY27vr4+Ghoa0i7DzKw0LXwV/vhzeG8m9KiDA/8dDjgLetalXZlZlyRpTkTUt55eEqdNzcysBAwcCSc/nL0ubtA/w/PXw+R94LmfwJc+uWFWKhzezMzsq/rvBxMehLNezJ5WfWFSNsTNvBpW+659s7Q5vJmZWX59h8KJ98E5f4I9x8CsX2ZD3DNXwqrP0q7OrMtyeDMzs7Z9ba/s2xrOfRX2Ggcv/yp7d+pTl8Dni9OuzqzLcXgzM7OO6bMnHH9b9v2pQ/8FGu7KPifuif+AFR+nXZ1Zl+HwZmZmm2enPWDczXD+XBh2Esy9D24YDo+dC8vfT7s6s4rn8GZmZltmh0Hwnclw4Tyo/z7Mnwo31sP0s2DZX9OuzqxiObyZmdnW2X4AfHsSXDQfDjgb/vIY3LQ/PHw6LHkr7erMKo7Dm5mZFUbvXeDI6+GiBTDqAnjnD/CrA+GhU+DT+WlXZ1YxHN7MzKywtu0D37omG+K+8UN4/3m47Rvw4ARYPDft6szKnsObmZl1jl47wTevyoa40VfAR7PgjkPhtydkX8VlZlvE7zY1M7PiWLsSXrsDZt8Ea5Zn395w8CWw26i0KzPrkEdfX8ykZ97hkxVr6FfXk4vHDOHY4f07bXuberepw5uZmRVX42pouBtm3wBfLIVBB8Ehl/Doij2YNOPdon0xmm2OR19fzOXTF7CmKbNhWs/aav77+H067ffU4Q2HNzOzkrLuS5h7L7w0GVZ/xtzYk182HccfW4YC6vQvxq6g2EeKUhEBmXXQvBaa1/9szP7MNG4c3vCzvWVajSfz3/x4CVWZRrrRRHc1cVTj9XzOtvSv68msyw7rlF3bVHir6ZStmZmZtafbNnDgObDfafziZ1cxoWka93b7GR+39GEV29BENdWP18K8naGqJvuprs0/vGG8Fqqqc4ZroHr9MrVfHd4wrzanfXWy3OZsp2bjuJR2r27Q+kjR4hVruHz6AoDCBbjNCk75QlEboSonOLW7jq0mqOkBNd1zfnbPGe/BsuaeNLIdjdTSGN0Isv/Wn6xYU4Dtbx6HNzMzS1dtD25ePZrbOIjx1S9wcNV8ammmhgy1zRlA2S/qltXQ0gyZZmhpyhlOxjNN0JLZOK+lufj7sslg2IGQ2FY4zTvcdlB9/al3OTTTTHNVNdW00J0murc08bcnn4U1g9o5GtU6SLURzLZa+8GJHnWtlun21fHqbptYRzJc3Wq89TJVNe0G7yt++hyL8wS1fnU9C9AHm8fhzczMUtevrieLV6xhSuYwpmQ2noLqX9eTWadt4SmpiI0hLtPUznAS/DJN7QfD3PZbtO4822luTNads832giptX/b0Y4BueWY0AzPWj7QOTq1DUfdOCE6tQlQHglMpuHjMkLzXvF08ZkjRa3F4MzOz1HXKF6OUPSpVXQu1xT860ulaWtoIk01MuHUW/7fqC2rJ0EQ1jdSyLmrZcfvePPmDw8sqOJWC9aeaS+EaQoc3MzNLXSl9MZaNqiqoSo5k5XHi2O55A/GlR+4D3XsXq8qKcuzw/iXxO+nwZmZmJaFUvhgrhQNx5XJ4MzMzq1AOxJXJr8cyMzMzKyMOb2ZmZmZlJPXwJulaSfMlzZM0Q1K/ZPpJyfT5kmZL2jenzYeSFiRt/MoEMzMz6zJSD2/ApIgYGhHDgCeAHyXTPwAOiYihwLXA7a3aHRoRw/K9NsLMzMysUqV+w0JErMwZ7UXy1MGImJ0z/WVgQDHrMjMzMytFqYc3AEnXAacCnwOH5lnkdOAPOeMBzJAUwG0R0fqonJmZmVlFKsppU0nPSnojz2ccQERcGREDgfuB81q1PZRseLs0Z/KoiBgBjAXOlXRwG9ueKKlBUsPSpUsLvm9mZmZmxaSItt+NVkySBgFPRsTeyfhQ4BFgbES8u4k2/wWsjohftLf++vr6aGjw/Q1mZmZW+iTNyXdtf+qnTSUNjoi/JqPHAG8n03cFpgOn5AY3Sb2AqohYlQwfAVzTkW3NmTNnmaSPCroD/2hnYFknb6OrcZ8Wlvuz8NynheX+LDz3aWEVqz8H5ZuY+pE3SdOAIUAL8BFwdkQslnQnMD6ZBtAcEfWSdid7NA6y4fOBiLiu2HVviqQG3wFbWO7TwnJ/Fp77tLDcn4XnPi2stPsz9SNvETF+E9PPAM7IM/19YN9/bGFmZmZW+UrhOW9mZmZm1kEOb4Xnx5YUnvu0sNyfhec+LSz3Z+G5Twsr1f5M/Zo3MzMzM+s4H3kzMzMzKyMObwUk6UhJ70h6T9JladdT7iTdLWmJpDfSrqUSSBoo6X8lvSXpTUkXpl1TOZPUQ9Krkv6c9OeP066pUkiqlvS6pCfSrqXcSfpQ0gJJ8yT5QacFIKlO0sOS3k7+P/160WvwadPCkFQNvAt8C1gEvAZMiIi/pFpYGUvenLEa+M36BzfblpPUF+gbEXMl9QbmAMf6d3TLSBLQKyJWS6oFXgIujIiXUy6t7En6AVAPbBcRR6ddTzmT9CFQHxF+xluBSLoXeDEi7pTUDdgmIlYUswYfeSuckcB7EfF+RKwDpgDjUq6prEXEC8DytOuoFBHxaUTMTYZXAW8B/dOtqnxF1upktDb5+K/hrSRpAHAUcGfatZi1Jmk74GDgLoCIWFfs4AYOb4XUH1iYM74IfzFaiZK0GzAceCXdSspbcnpvHrAEmBkR7s+tNxm4hOyD223rBTBD0hxJE9MupgLsDiwF7klO7d+ZvO2pqBzeCkd5pvmvcCs5krYFpgEXRcTKtOspZxGRiYhhwABgpCSf3t8Kko4GlkTEnLRrqSCjImIEMBY4N7kcxbZcDTACuCUihgNfAEW/xt3hrXAWAQNzxgcAn6RUi1leybVZ04D7I2J62vVUiuS0yfPAkSmXUu5GAcck12lNAQ6T9Nt0SypvEfFJ8nMJ2VdLjky3orK3CFiUc5T9YbJhrqgc3grnNWCwpH9KLmD8N+DxlGsy2yC5wP4u4K2I+J+06yl3kvpIqkuGewKHA2+nW1V5i4jLI2JAROxG9v/Q5yLi5JTLKluSeiU3J5Gc2jsC8N37WyEiPgMWShqSTPomUPSbvlJ/t2mliIhmSecBzwDVwN0R8WbKZZU1SQ8Co4GdJS0Cro6Iu9KtqqyNAk4BFiTXaQFcERFPpVhTOesL3JvcaV4F/C4i/GgLKyVfAx7J/t1GDfBARDydbkkV4Xzg/uRAzfvAacUuwI8KMTMzMysjPm1qZmZmVkYc3szMzMzKiMObmZmZWRlxeDMzMzMrIw5vZmZmZmXE4c3MujxJv5b0k7TrMDPrCIc3MzMzszLi8GZm1gkk+SHoZtYpHN7MrMuRNFzSXEmrJD0E9MiZd7SkeZJWSJotaWjOvBGSXk/aTZX00PrTrZJGS1ok6VJJnwH3dGB9/SRNk7RU0geSLiheL5hZuXJ4M7MuJXmlzaPAfcCOwFRgfDJvBHA3cBawE3Ab8Lik7km7R4BfJ+0eBI5rtfpdknmDgIntrK8K+D3wZ6A/2XckXiRpTOfsuZlVCr8ey8y6FEkHA1OA/pH8ByhpNvAc2YC1LCKuyln+HWAiEGQD24Ccdi8Bz0fEf0oaDcwAtouItcn8W9pY31pgakTsmjPvcmDPiCj6uxLNrHz4mgwz62r6AYvjq3+5fpT8HAR8V9L5OfO6JW0iT7uFrda9dH1w68D6MkA/SSty5lUDL27uDplZ1+LwZmZdzadAf0nKCWK7An8jG8aui4jrWjeSdEiedgOTduu1PpXR1vq+DnwQEYO3bnfMrKvxNW9m1tX8CWgGLpBUI+l4YGQy7w7gbEkHKKuXpKMk9U7aZYDzknbjctptSlvrexVYmdzg0FNStaS9Je3fKXttZhXD4c3MupSIWAccD3wP+DtwIjA9mdcAnAnclMx7L1kut93pwArgZOAJoLGNbbW1vgzwHWAY8AGwDLgT2L5Q+2pmlck3LJiZbSFJrwC3RsQ9addiZl2Hj7yZmXWQpEMk7ZKcNv0uMBR4Ou26zKxr8Q0LZmYdNwT4HbAt2RsVToiIT9Mtycy6Gp82NTMzMysjPm1qZmZmVkYc3szMzMzKiMObmZmZWRlxeDMzMzMrIw5vZmZmZmXE4c3MzMysjPw/lSBpQDXiE/kAAAAASUVORK5CYII=\n",
"text/plain": [
"