{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import arviz as az\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import pymc3 as pm\n",
"import statsmodels.api as sm\n",
"\n",
"# R-like interface, alternatively you can import statsmodels as import statsmodels.api as sm\n",
"import statsmodels.formula.api as smf\n",
"import theano\n",
"\n",
"from scipy import stats\n",
"from scipy.special import logsumexp"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format = 'retina'\n",
"az.style.use('arviz-darkgrid')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.1"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" species | \n",
" brain | \n",
" mass | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" afarensis | \n",
" 438 | \n",
" 37.0 | \n",
"
\n",
" \n",
" | 1 | \n",
" africanus | \n",
" 452 | \n",
" 35.5 | \n",
"
\n",
" \n",
" | 2 | \n",
" habilis | \n",
" 612 | \n",
" 34.5 | \n",
"
\n",
" \n",
" | 3 | \n",
" boisei | \n",
" 521 | \n",
" 41.5 | \n",
"
\n",
" \n",
" | 4 | \n",
" rudolfensis | \n",
" 752 | \n",
" 55.5 | \n",
"
\n",
" \n",
" | 5 | \n",
" ergaster | \n",
" 871 | \n",
" 61.0 | \n",
"
\n",
" \n",
" | 6 | \n",
" sapiens | \n",
" 1350 | \n",
" 53.5 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {
"image/png": {
"height": 317,
"width": 716
}
},
"output_type": "display_data"
}
],
"source": [
"f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(8,3))\n",
"ax1.scatter(d.mass, d.brain, alpha=0.8)\n",
"ax2.scatter(d.mass, d.brain, alpha=0.8)\n",
"for i in range(len(d)):\n",
" d_new = d.drop(d.index[-i])\n",
" m0 = smf.ols('brain ~ mass', d_new).fit()\n",
" # need to calculate regression line\n",
" # need to add intercept term explicitly\n",
" x = sm.add_constant(d_new.mass) # add constant to new data frame with mass\n",
" x_pred = pd.DataFrame({'mass': np.linspace(x.mass.min() - 10, x.mass.max() + 10, 50)}) # create linspace dataframe\n",
" x_pred2 = sm.add_constant(x_pred) # add constant to newly created linspace dataframe\n",
" y_pred = m0.predict(x_pred2) # calculate predicted values\n",
" ax1.plot(x_pred, y_pred, 'gray', alpha=.5)\n",
" ax1.set_ylabel('body mass (kg)', fontsize=12);\n",
" ax1.set_xlabel('brain volume (cc)', fontsize=12)\n",
" ax1.set_title('Underfit model')\n",
" \n",
" # fifth order model\n",
" m1 = smf.ols('brain ~ mass + I(mass**2) + I(mass**3) + I(mass**4) + I(mass**5)', data=d_new).fit()\n",
" x = sm.add_constant(d_new.mass) # add constant to new data frame with mass\n",
" x_pred = pd.DataFrame({'mass': np.linspace(x.mass.min()-10, x.mass.max()+10, 200)}) # create linspace dataframe\n",
" x_pred2 = sm.add_constant(x_pred) # add constant to newly created linspace dataframe\n",
" y_pred = m1.predict(x_pred2) # calculate predicted values from fitted model\n",
" ax2.plot(x_pred, y_pred, 'gray', alpha=.5)\n",
" ax2.set_xlim(32,62)\n",
" ax2.set_ylim(-250, 2200)\n",
" ax2.set_ylabel('body mass (kg)', fontsize=12);\n",
" ax2.set_xlabel('brain volume (cc)', fontsize=12)\n",
" ax2.set_title('Overfit model')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.9"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6108643020548935"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = (0.3, 0.7)\n",
"-sum(p * np.log(p))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.10"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"94.92498968588757"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# fit model\n",
"m_6_1 = smf.ols('brain ~ mass', data=d).fit()\n",
"\n",
"#compute de deviance by cheating\n",
"-2 * m_6_1.llf"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.11"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [sigma, b, a]\n",
"Sampling 2 chains: 100%|██████████| 14000/14000 [00:10<00:00, 1284.28draws/s]\n"
]
}
],
"source": [
"# standarize the mass before fitting\n",
"d['mass_s'] = d['mass'] - np.mean(d['mass'] / np.std(d['mass']))\n",
"\n",
"with pm.Model() as m_6_8 :\n",
" a = pm.Normal('a', mu=np.mean(d['brain']), sd=10)\n",
" b = pm.Normal('b', mu=0, sd=10)\n",
" sigma = pm.Uniform('sigma', 0, np.std(d['brain']) * 10)\n",
" mu = pm.Deterministic('mu', a + b * d['mass_s'])\n",
" brain = pm.Normal('brain', mu = mu, sd = sigma, observed = d['brain'])\n",
" m_6_8 = pm.sample(2000, tune=5000) "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"theta = az.summary(m_6_8)['mean'][:3]"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"100.60410127191153"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#compute deviance\n",
"dev = - 2 * sum(stats.norm.logpdf(d['brain'], loc = theta[0] + theta[1] * d['mass_s'] , scale = theta[2]))\n",
"dev "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.12"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[This](https://github.com/rmcelreath/rethinking/blob/a309712d904d1db7af1e08a76c521ab994006fd5/R/sim_train_test.R) is the original function."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# This function only works with number of parameters >= 2 \n",
"\n",
"def sim_train_test(N=20, k=3, rho=[0.15, -0.4], b_sigma=100):\n",
" \n",
" n_dim = 1 + len(rho)\n",
" if n_dim < k:\n",
" n_dim = k\n",
" Rho = np.diag(np.ones(n_dim))\n",
" Rho[0, 1:3:1] = rho\n",
" i_lower = np.tril_indices(n_dim, -1)\n",
" Rho[i_lower] = Rho.T[i_lower]\n",
" \n",
" x_train = stats.multivariate_normal.rvs(cov=Rho, size=N)\n",
" x_test = stats.multivariate_normal.rvs(cov=Rho, size=N)\n",
" \n",
" mm_train = np.ones((N,1))\n",
" \n",
" np.concatenate([mm_train, x_train[:, 1:k]], axis=1)\n",
" \n",
" #Using pymc3\n",
" \n",
" with pm.Model() as m_sim:\n",
" vec_V = pm.MvNormal('vec_V', mu=0, cov=b_sigma * np.eye(n_dim), \n",
" shape=(1, n_dim), testval=np.random.randn(1, n_dim)*.01)\n",
" mu = pm.Deterministic('mu', 0 + pm.math.dot(x_train, vec_V.T))\n",
" y = pm.Normal('y', mu=mu, sd=1, observed=x_train[:, 0])\n",
"\n",
" with m_sim:\n",
" trace_m_sim = pm.sample()\n",
" \n",
" vec = pm.summary(trace_m_sim)['mean'][:n_dim]\n",
" vec = np.array([i for i in vec]).reshape(n_dim, -1)\n",
" \n",
" dev_train = - 2 * sum(stats.norm.logpdf(x_train, loc = np.matmul(x_train, vec), scale = 1)) \n",
" \n",
" mm_test = np.ones((N,1))\n",
" \n",
" mm_test = np.concatenate([mm_test, x_test[:, 1:k +1]], axis=1)\n",
" \n",
" dev_test = - 2 * sum(stats.norm.logpdf(x_test[:,0], loc = np.matmul(mm_test, vec), scale = 1)) \n",
" \n",
" return np.mean(dev_train), np.mean(dev_test)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 696.53draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 670.08draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 629.84draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 514.13draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 603.98draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 490.54draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 592.44draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 538.58draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 648.23draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 603.95draws/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"3\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 618.29draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 639.22draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 481.07draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 519.18draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 518.25draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 571.07draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 542.79draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 462.96draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:05<00:00, 390.67draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 684.51draws/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"4\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 804.55draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 611.91draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 458.24draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 689.49draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 556.63draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 720.40draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 570.87draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 737.06draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 754.59draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 613.29draws/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"5\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 723.48draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 803.98draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 718.26draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 503.21draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 704.96draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 667.57draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 527.39draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 734.26draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 548.40draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 667.49draws/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"6\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 721.66draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 819.60draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 480.44draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 683.64draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 727.10draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 769.77draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:01<00:00, 1027.99draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 759.98draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 677.01draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [vec_V]\n",
"Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 797.01draws/s]\n"
]
}
],
"source": [
"n = 20\n",
"tries = 10\n",
"param = 6\n",
"r = np.zeros(shape=(param - 1, 4))\n",
"\n",
"train = []\n",
"test = []\n",
"\n",
"for j in range(2, param + 1):\n",
" print(j)\n",
" for i in range(1, tries + 1):\n",
" tr, te = sim_train_test(N=n, k=param)\n",
" train.append(tr), test.append(te)\n",
" r[j -2, :] = np.mean(train), np.std(train, ddof=1), np.mean(test), np.std(test, ddof=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.14"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABscAAARBCAYAAAB5Fl5TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl0VeW5B+A3kAECQhgElEnUMrWiFLRapaioiLP0VhdesdpaqKi0DtepvbfWqS2rWkSp4qy9UocWEBVRi0WLFNEiFREQJ2Se5zGQ3D9cnuthykk4IYb9PGux1nl3vuE9YbkN+Z29d05paWlpAAAAAAAAQALUqOoGAAAAAAAAYG8RjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGLlV3QAAAABUpW3btsVnn30Ws2fPjiVLlsT69eujsLAwioqKomPHjnHooYdGTk5OVvbasmVL/Otf/4r58+fHihUroqioKJo1axZdu3aNwsLCrOwBAADsnnAMAADga2revHnRo0ePHY4//PDDcdxxx2W8zogRI+LGG29M1UOHDo2TTjopKz1WV8uWLYtXXnkl3njjjXj77bdj3bp1uxzbsGHDOO+88+KHP/xhNGzYsEL7rVu3LgYPHhyjR4+O1atX7/D1wsLC6NmzZ1xzzTWx//77V2gPAAAgM26rCAAAUM0MHjy4qluo1kaNGhXf+9734te//nX8/e9/320wFhGxYsWKuP/+++OMM86I8ePHl3u/6dOnx1lnnRV/+tOfdhqMRURs2LAhRo4cGWeddVZMnDix3HsAAACZc+UYAABANTNt2rT429/+lvirvypq1apVsW3btrRjhYWF0blz52jVqlU0aNAgNm7cGDNmzIi33347NXb58uVx+eWXx5AhQ3Z6Rd/OLFiwIPr16xfLli1LHatfv36ccMIJ0axZs1i6dGm88cYbsXTp0oj4Ioi7/PLL489//nO0b98+S+8YAAD4KuEYAABANXT33XfHiSeeGDVquCFIReXm5sYJJ5wQ559/fhx99NGRl5e3w5i5c+fGTTfdFJMnT46IiK1bt8a1114bL7/8cjRp0mS365eWlsbAgQPTgrEzzjgjbrnllqhTp07q2JYtW+LOO++Mxx57LCK+uIpswIAB8dJLL0VBQUEW3ikAAPBV/hUFAABQTeTm/v/nGz/88MN48cUXq7Cb6qtmzZpxxhlnxNixY+Pee++Nbt267TQYi4ho2bJlPPzww3HkkUemjm3YsCHuu+++MvcZO3ZsTJs2LVV/97vfjd///vdpwVhERH5+ftx4443Ru3fv1LH58+fHn//85/K+NQAAIAPCMQAAgGri1FNPjcLCwlR9zz33xNatW6uwo+qpT58+ceedd0bLli0zGp+fnx+33npr5OTkpI698sorUVpautt5Xw3QatSoETfffHPaGtu74YYbom7duqn6wQcf9PcLAACVQDgGAABQTTRs2DAuuuiiVD1nzpwYOXJkFXZUPX31CrxMtWnTJjp27Jiqly1bFgsWLNjl+M8++yxmzZqVqo877rho3br1bveoX79+nH766Wl7vP322+XuFQAA2D3hGAAAQDXy4x//OOrVq5eqhw4dGlu2bKnCjpJj+3Br+fLluxw7bty4tLpnz54Z7dGrV6/drgMAAOw54RgAAEA1Uq9evfjRj36UqhcuXBhPPfVUFXaUHOvXr0+rd3cF2jvvvJNWd+7cOaM9OnXqFDVr1tzlOgAAwJ4TjgEAAFQzF110UTRs2DBVDxs2LDZu3FiFHe37SktLY/bs2WnHmjVrtsvxn3zySep1rVq14uCDD85onzp16kSrVq1S9aefflrms80AAIDyEY4BAABUM3Xq1In+/fun6mXLlsWf/vSnKuxo3/f222+nPWOsbdu2aQHlVxUXF8e8efNSdbNmzSInJyfjvQ488MDU602bNsX8+fMr0DEAALAr5X8KMQAAAFWuT58+8cgjj8TixYsjIuKhhx6KPn36xH777bfXexk1alTcdNNNlbrHBx98UKnrl+Xee+9Nq0899dRdjl22bFls3bo1VR9wwAHl2qtp06Zp9aJFi6JFixblWgMAANg1V44BAABUQwUFBTFgwIBUvXr16njkkUeqpJeSkpLYtm1bpf6pSs8++2y89dZbqbqoqCguvPDCXY7fsGFDWl1YWFiu/erUqbPb9QAAgD0jHAMAAKimvv/970fLli1T9eOPPx4rVqyowo72PbNnz47bb7897dj1118f9evX3+Wc9evXp9UFBQXl2nP78cIxAADILrdVBAAAqKby8vLiiiuuiOuvvz4ivghlHnzwwVS9t/Tu3Tt69+69V/fcG1auXBmXX355bNy4MXWsZ8+eZb7XLVu2pNV5eXnl2jc/Pz+t3rRpU7nmAwAAu+fKMQAAgGrsrLPOikMPPTRVP/nkk6nnkFFxGzdujJ/+9KcxZ86c1LH27dvHHXfcUebc7cOt4uLicu29fbhWq1atcs0HAAB2TzgGAABQjdWoUSMGDhyYqjdv3hz33XdfFXZU/RUXF8eVV14ZU6dOTR1r0aJFPPjgg1G3bt0y52//zLDNmzeXa//tx5f3mWUAAMDuCccAAACquVNOOSW++c1vpuq//OUvMXfu3CrsqPoqKSmJ6667Lv7xj3+kju2///7x6KOPRpMmTTJaY/swa/tnkJVl+/HCMQAAyC7hGAAAQDWXk5MTP/vZz1J1cXFxDB06tAo7qr5uvvnmGDNmTKquX79+PPTQQ9GqVauM12jcuHHk5v7/I74XLVpUrh62H9+sWbNyzQcAAHYvt+whAAAAfN117949vv3tb8eUKVMiImL06NHxk5/8JA455JBK33vUqFFx0003VeoeH3zwQaWuHxExaNCgePrpp1N1YWFhPPDAA9G+fftyrZOXlxctWrSIzz77LCK+CLtKS0sjJycno/kLFy5Mva5du3Y0b968XPsDAAC758oxAACAfcRVV12Ver1t27YYMmTIXtm3pKQktm3bVql/Ktt9990XDz/8cKrOz8+PP/7xj3HEEUdUaL2vhpKbNm2KTz75JKN569atizlz5qTqgw46KONQDQAAyIxwDAAAYB9x1FFHxbHHHpuqX3755ZgxY0YVdlQ9PPnkkzF48OBUnZubG4MHD45jjjmmwmt26dIlrf7yir6yvPfee1FSUpKqu3btWuEeAACAnXNbRQAAgH3Iz3/+83jzzTcjIqK0tDQGDx4cPXv2rNQ9e/fuHb17967UPSrLc889F7feemuqrlGjRvzud7+LHj167NG6PXr0iEGDBqXql19+OX7wgx+UOW/s2LE7rAMAAGSXK8cAAAD2IZ06dUoLVMaPHx9Tp06two6+vsaNGxc33XRTlJaWpo7dcsstccYZZ+zx2gcddFDas8refPPNtNsl7szq1atjzJgxqbpx48Zx5JFH7nEvAABAOuEYAADAPuZnP/tZ1Kjx///c++tf/1qF3Xw9TZo0Ka666qrYunVr6tiNN96Y0dVdmbrssstSr0tKSuLmm29OC+K299vf/jbWrl2bqn/yk59Ebq4bvgAAQLYJxwAAAPYx7dq1i9NOOy1VfzUAIuL999+PAQMGxObNm1PHBg4cGBdffHFW9+nZs2ccdthhqXrixIlx7bXXxvr169PGbdmyJX7zm9/EiBEjUseaN28effr0yWo/AADAF3wEDQAAYB905ZVXxtixYwVjO/Hkk0/uEFANHTo0hg4dWq517rjjjjjnnHN2+fWcnJwYMmRI/Md//EcsX748IiJeeOGF+Mc//hEnnnhiNG3aNJYtWxavv/56LF26NDWvsLAw/vjHP0ZBQUG5+gEAADIjHAMAANgHHXTQQXHOOefEX/7yl6pu5WtnZ7c23LZtW7nXKSkpKXPMgQceGA888EAMHDgw5s+fHxFfPFts5MiROx3foEGDuPPOO9OeVwYAAGSX2yoCAADsoy6//PLIz8+v6jYS71vf+lY899xz0bdv36hfv/5Ox9SuXTvOOeeceP755+PYY4/dyx0CAECy5JTu7mnAAAAAQNZs2bIl3nnnnZg/f36sWLEi6tevHwcccEB07do16tSpU9XtAQBAIgjHAAAAAAAASAy3VQQAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGMIxAAAAAAAAEkM4BgAAAAAAQGIIxwAAAAAAAEiM3KpuAHZm5cqVVd3C105OTk4UFRVFRMSqVauitLS0ijsCqN6cVwGyy3kVILucVwGyy3m1+mvQoEHW1nLlGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxMit6gb2BQsXLoxp06bFggULYsOGDVFQUBCNGzeONm3aRPv27SM/P7+qWwQAAAAAACCEYxVWUlISo0ePjieeeCKmT5++y3F5eXnRuXPn6NevX3Tr1m2nY0488cSYP3/+Hvf0m9/8Jnr37r3H63ypb9++MXny5ArNffTRR+O73/1u1noBAAAAAADIBuFYBcybNy+uueaamDp1aplji4uLY/LkyXH44YfvMhzLlsLCwkpdHwAAAAAAoLoTjpXT7Nmz45JLLomlS5emjuXk5ESnTp2iffv20ahRo9i8eXPqVotz584tc82aNWtGzZo1y9VHaWlplJSUpOrCwsJKD9/K02NOTk4ldgIAAAAAAFAxwrFyWLFiRVx66aVpwdhJJ50UN954Y7Ro0WKncz788MMYMWJE7Lfffrtc99VXXy13L4899lj85je/SdWnnHJK1KlTp9zrZOqoo46KP/3pT5W2PgAAAAAAwN4gHCuH22+/PRYtWpSqf/rTn8ZVV1212zlt27aNG264Ieu9jBw5Mq0+99xzs74HAAAAAADAvqZGVTdQXUycODFeeOGFVN2jR48yg7HKMmPGjJg5c2aqbt68eXznO9+pkl4AAAAAAACqE+FYhh544IHU67y8vLjxxhurrJftrxo755xzPOMLAAAAAAAgA8KxDMydOzcmTZqUqrt37x4tW7askl62bt2adgVbTk6OWyoCAAAAAABkSDiWgTFjxkRpaWmqPv3006usl9dffz2WL1+eqrt06VJlQR0AAAAAAEB1IxzLwNSpU9Pqrl27VlEnEaNGjUqrXTUGAAAAAACQudyqbqA6eP/991OvmzZtGk2aNImIiDlz5sSIESNiwoQJsXDhwti4cWM0aNAg2rRpE8cee2ycffbZ0ahRo6z1sXLlyvj73/+eqmvXrh2nnnpq1tbfnQULFsQNN9wQ06ZNi6VLl0ZxcXE0aNAgmjZtGl27do3jjz8+unTpsld6AQAAAAAAqCjhWBnWrl0bS5YsSdUtW7aMrVu3xv333x/33XdfbN26NW38hg0bYv78+TFhwoS455574tJLL40BAwZETk7OHvfy4osvRnFxcao+5ZRTom7dunu8bibmzZsX8+bNSzv25XudMmVKPPDAA9G5c+f4xS9+EYcddthe6QkAAAAAAKC8hGNlWLVqVVrduHHjuOmmm+K5554rc+6GDRtiyJAhMXPmzPjDH/4Qubl79u0eOXJkWv11u6Xiu+++G3369Ilf/epX8YMf/GCP1spGmLiv+er3xPcHYM85rwJkl/MqQHY5rwJkl/MqXyUcK8PatWvT6jfffDN1rFatWnHxxRfHaaedFq1bt47i4uKYNWtWPPPMMzF69OgoLS2NiIhXXnkl7rrrrrjuuusq3MdHH32UdnvH5s2bx9FHH13h9TJVVFQUJ554YnTr1i3atWsXTZs2jVq1asWaNWti9uzZMX78+HjmmWdi3bp1ERFRXFwc//M//xNFRUVx8skn79G+7Fr9+vWrugWAfYrzKkB2Oa8CZJfzKkB2Oa+SU/plgsNOvfPOO/Gf//mfOxwvKiqKxx57LDp06LDTeaNHj47rr78+SkpKIuKLJHrUqFHRvn37CvUxaNCgePjhh1P1gAED4mc/+1mF1srU5MmT4/DDD4+CgoLdjlu2bFlcddVVMXny5NSxunXrxt/+9rdo0KBBpfYIAAAAAABQHsKxMrz33ns7vUXgkCFDomfPnrude8cdd8Tjjz+eqs8+++wYNGhQuXvYtm1bdO/ePZYuXZo69uqrr0arVq3KvVZl2bhxY1xwwQXxwQcfpI79+Mc/rvDVctvfzpIvAtYvP9GwevXq8J8uwJ5xXgXILudVgOxyXgXILufV6i+bd5xzW8Uy1KlTZ4djhxxySJnBWEREv379Yvjw4VFcXBwREePHj4+SkpKoUaNGuXqYMGFCWjDWtWvXr1UwFhFRu3bt+NWvfhXnn39+6tiYMWMqHI45Me1eaWmp7xFAFjmvAmSX8ypAdjmvAmSX8yrlS2kSaGfh2AknnJDR3MaNG0enTp1S9erVq+Ojjz4qdw+jRo1Kq88999xyr7E3HHHEEXHooYem6oULF8Znn31WdQ0BAAAAAABsRzhWhkaNGkVeXl7asW984xsZz2/btm1avXjx4nLtv3bt2hg3blyqrl27dpx66qnlWmNvOuKII9LqBQsWVFEnAAAAAAAAOxKOlSEvL2+HWxh+eV/STGw/dvXq1eXa/8UXX4zNmzen6pNPPjnq1q1brjX2pkaNGqXVK1eurKJOAAAAAAAAdiQcy8BXbxUYEbFly5aM524/tqCgoFx7b39Lxd69e5dr/t62cePGtLq87xcAAAAAAKAyCccycNRRR6XV5bk14qJFi9LqBg0aZDz3s88+i3fffTdVH3jggXH00UdnPL8qzJ07N61u2LBhFXUCAAAAAACwI+FYBk466aTIyclJ1VOmTMloXmlpaUydOjVV16xZM9q3b5/xviNHjkyrzz777LQ+vm42btwYb7/9dqrOzc0t1/PZAAAAAAAAKptwLAPNmjWLb3/726n6tddeixUrVpQ5b8KECbFgwYJU3alTp4yfF1ZaWhqjR49OO3buuedm2HHVePzxx2PdunWp+vDDD4/99tuvCjsCAAAAAABIJxzL0JVXXpl6vXnz5rj11lt3O379+vVx++23px3r27dvxvtNmjQpLVjr0qVLtG7dOuP525s3b160a9cu9efEE0/c7fhZs2aVa/0JEybEvffem3bs4osvLm+bAAAAAAAAlUo4lqFjjjkmjj/++FQ9ZsyYuP7662Pt2rU7jP3888/jkksuiU8//TR17LDDDotevXplvN+IESPS6r191diFF14Y/fr1i1dffTU2b968y3Fr1qyJIUOGRP/+/aO4uDh1/KijjopTTjllb7QKAAAAAACQsdyqbqA6GTRoUJx//vmp0GvUqFExbty46NatW7Rq1SqKi4vjww8/jEmTJqUFRQ0bNowhQ4ZEjRqZZZHr16+PV199NVXXrl27XMFaNpSWlsbrr78er7/+etSuXTvatWsXhx56aBQVFUWtWrVi7dq18dFHH8WUKVNi48aNaXPbtGkT99xzz17tFwAAAAAAIBPCsXKoX79+PPTQQzFw4MCYPn16RESsXbs2xowZs8s5Bx98cAwbNiwOPPDAjPcZO3ZsWuB08sknZ/ysssqwcePGmDp1akydOrXMsT179oxbbrklioqK9kJnAAAAAAAA5eO2iuXUokWLePrpp+Pqq6+O5s2b73JckyZN4r/+679i5MiR0apVq3LtMWrUqLS6d+/eFep1T/Tv3z++853vRGFhYZlj8/Pzo0ePHvHYY4/FkCFDBGMAAAAAAMDXVk5paWlpVTdRXZWWlsa0adPi008/jaVLl0ZOTk40bNgwOnToEO3bt6/q9rKipKQk5syZE3PmzIlFixbF2rVro7i4OAoLC6N+/frRpk2b6NixY+Tn52d135UrV2Z1vX1BTk5OKnhctWpV+E8XYM84rwJkl/MqQHY5rwJkl/Nq9degQYOsreW2insgJycnOnXqFJ06darqVipNjRo1ok2bNtGmTZuqbgUAAAAAAGCPua0iAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEyK3qBoDqZc2aNdG/f/+0Y8OGDYt69epVUUcAAAAAAJA5V44BAAAAAACQGMIxAAAAAAAAEkM4BgAAAAAAQGIIxwAAAAAAAEgM4RgAAAAAAACJIRwDAAAAAAAgMYRjAAAAAAAAJIZwDAAAAAAAgMQQjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGMIxAAAAAAAAEkM4BgAAAAAAQGIIxwAAAAAAAEgM4RgAAAAAAACJIRwDAAAAAAAgMYRjAAAAAAAAJIZwDAAAAAAAgMQQjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGMIxAAAAAAAAEkM4BgAAAAAAQGIIxwAAAAAAAEgM4RgAAAAAAACJIRwDAAAAAAAgMYRjAAAAAAAAJIZwDAAAAAAAgMQQjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEyK3qBgAAAACAr581a9ZE//79044NGzYs6tWrV0UdAUB2uHIMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGMIxAAAAAAAAEiO3qhvYFyxcuDCmTZsWCxYsiA0bNkRBQUE0btw42rRpE+3bt4/8/PyqbjFrFi9eHO+9914sXrw4Nm7cGE2bNo02bdrEYYcdVtWtAQAAAAAAlEk4VkElJSUxevToeOKJJ2L69Om7HJeXlxedO3eOfv36Rbdu3XY5rl27dhXuZfr06ZGbW7l/le+9914MHjw4/vnPf0ZJSckOX2/VqlX07ds3+vbtGzk5OZXaCwAAAAAAQEW5rWIFzJs3L/r06RPXX3/9boOxiIji4uKYPHlyvPXWW3upu+x78MEHo0+fPvHmm2/uNBiLiPj888/j9ttvj0suuSRWrVq1lzsEAAAAAADIjCvHymn27NlxySWXxNKlS1PHcnJyolOnTtG+ffto1KhRbN68OXWrxblz55Z7j5ycnKhRI/PcsjKv1Hrqqafi97//fdqxDh06RJcuXaJ27drx8ccfxxtvvBFbt26NiIh//vOfMXDgwHj44YcjLy+v0voCAAAAAACoCOFYOaxYsSIuvfTStGDspJNOihtvvDFatGix0zkffvhhjBgxIvbbb7+M97n88svjyiuv3ON+99SMGTPi1ltvTdX5+flxxx13xJlnnpk2bu7cuTFgwID48MMPIyLirbfeiiFDhsQ111yzV/sFAAAAAAAoi3CsHG6//fZYtGhRqv7pT38aV1111W7ntG3bNm644YbKbq1S/OEPf0hdERYRceutt+4QjEVEtGzZMh5//PE444wzYvny5RER8cQTT0Tfvn2jSZMme61fAKB62bx5c7z33nvx7rvvxpQpU2LKlCnx8ccfR2lpaUR88cGkbFq5cmU8+eSTMW7cuJgxY0asWrUqCgoKYv/994/9998/jjjiiDjuuOPimGOOiYYNG+4w/8wzz4w333xzh+O1a9eOevXqRVFRUXTo0CGOOOKIOO200+LQQw/Nav8AAABUb2vWrIn+/funHRs2bFjUq1evijpKLuFYhiZOnBgvvPBCqu7Ro0eZwVh19sEHH8Trr7+eqo888sg455xzdjm+YcOGcfXVV8cvfvGLiIjYtGlTPPLII9U2GAQAKt/VV18df/7zn/fKXq+88kpcfvnlqQ/yfGnLli2xdu3a+OSTT+Ktt96KYcOGRZcuXeLVV1/NeO2NGzfGxo0bY/HixTFr1qwYNWpU/PrXv45u3brFbbfdFt/61rey/XYAAACAPSAcy9ADDzyQep2Xlxc33nhjFXZT+caOHZtWX3DBBWXOOfPMM+N3v/tdrFmzJrWGcAwAyETdunXj8MMPj48++igWL16c1bUnTZoUffv2jeLi4qhZs2Z8//vfj549e0br1q2jZs2asWTJknjvvfdi3LhxMXny5IzWnDBhQur1tm3bYvXq1bFgwYJ455134vnnn4/FixfHG2+8ESeddFL89re/jYsvvjir7wkAAACoOOFYBubOnRuTJk1K1d27d4+WLVtWYUeV77XXXku9zsvLix49epQ5p6CgII4//vgYPXp0REQsXLgwPvjgg+jYsWOl9QkAVJ3hw4fHFVdcES1btox///vf5Z7fo0ePOPbYY6Nz587Rrl27qFGjRpx55plZD8d++ctfpoKxZ599No4//vgdxpx88slxzTXXxNy5c9Ount+VXf18c95558Wtt94aw4YNi9tuuy22bNkS1157bRxwwAHRs2fPPX0rAAAosHd9AAAgAElEQVQAQBbUqOoGqoMxY8aknn0REXH66adXYTeVb/Xq1TF79uxU3bFjxygoKMhobufOndPqt99+O6u9AQD7jt69e8cFF1wQHTp0iBo1KufH0oULF8aUKVMi4ouf4XYWjH1Vy5Yt48ILL9yjPQsKCmLgwIFx//33R0RESUlJ/PznP49Nmzbt0boAAABAdgjHMjB16tS0umvXrlXUyd7x8ccfp9UdOnTIeO72n6Lefi0AgL1p3rx5qdcHH3zwXt27d+/ecfbZZ0dExOLFi2P48OF7dX8AAABg59xWMQPvv/9+6nXTpk2jSZMmERExZ86cGDFiREyYMCEWLlwYGzdujAYNGkSbNm3i2GOPjbPPPjsaNWpU7v3eeuutmDlzZsycOTNWrFgRubm50aBBg2jdunUceeSRceqpp0arVq2y9v6298knn6TVBx54YMZzDzjggLT6008/zUpPAAAVkZ+fn3r94Ycf7vX9L7vssnjuueciIuKFF16IH/3oR3u9BwAAACCdcKwMa9eujSVLlqTqli1bxtatW+P++++P++67L7Zu3Zo2fsOGDTF//vyYMGFC3HPPPXHppZfGgAEDIicnJ+M9d3YrwjVr1sScOXPijTfeiMGDB0evXr3ipptuqlD4Vpbtn/PRrFmzjOfuv//+kZubm/q+LFq0KKu9AQCUR9u2baNWrVqxadOmeOmll+LZZ5+NH/zgB3tt/65du0bdunVj3bp18c4778TWrVsjN9eP4AAAAFCV/Mu8DKtWrUqrGzduHDfddFPqE8C7s2HDhhgyZEjMnDkz/vCHP2TtFyHbtm2LF154ISZPnhxDhw6NTp06ZWXdL23YsCGtrlOnTsZza9SoEbVq1Yp169btdK1MlSdMTIqvfk+q8vtTs2bNaNy48Q7H/J0B1c3X5bxanVX29zAbaxYWFkbfvn3jwQcfjNLS0ujfv3/cddddcfLJJ8eRRx4ZXbp0iebNm1dabzVr1oxvfetbMWnSpFi3bl0sXrw4WrRoUe79oDpwXgXIrq/DedXvAIB9ifMqXyUcK8PatWvT6jfffDN1rFatWnHxxRfHaaedFq1bt47i4uKYNWtWPPPMMzF69OgoLS2NiIhXXnkl7rrrrrjuuut2u1fr1q2jR48ecfTRR8ehhx4aDRs2jJo1a8aKFSvi/fffjzFjxsTYsWNj27ZtERGxZMmS6NevXzzzzDNZvc3i9oHWV29HlImCgoI9DseKiooqNC8p6tevX2V7FxUVxf/+7/9W2f4AlaEqz6vVWWFhYUR88eGYbP2/+6sfJsrWmnfffXfMmzcvXnrppYiImDVrVsyaNSv19VatWkWPHj3ioosuiuOPPz7rvTVt2jT1uqSkxM85JILzKkB2VdV51e8AgH2V8yo1qrqBr7vtw50vg7GioqJ46qmn4qqrrop27dpFrVq1Yr/99ouuXbvGoEGDYtCgQVGjxv9/ex955JGYOXPmLvd5+OGH4+WXX47rr78+unfvHs2bN4/atWtHfn5+NGvWLE466aS466674qmnnkp7rtfKlSvLDN3Ka/PmzWl1Xl5eueZ/NUzbtGlTVnoCAKio2rVrx4svvhjDhw+PY445Zoevf/755/Hoo4/GCSecEKeeemosXbo0q/vXrVs39Xr7D14BAAAAe58rx8qwq6umbrnllujQocMu55111lnx/vvvx+OPPx4REaWlpfHII4/EoEGDdjr+uOOOy6ifTp06xSOPPBLnnXde6pcr7777brz++uvRvXv3jNYoS0FBQVpdXFxcrvlbtmxJva5Vq1aFetj+dpZ8canvl59oWL16derKxL1t3bp1ccMNN6Qd++1vf5v2iz+A6uDrcF6tDufUBg0alDlmzpw5u70FxNChQ+OCCy7IaL+vPs812z8P9OrVK3r16hVLliyJSZMmxdSpU2Pq1KkxadKk2LhxY0REvPzyy/G9730vXn311R3+Hira24oVK9JqP+ewr/o6nFcB9iVfh/Nqdfh5FSBTzqvVXzbvxCIcK8POnrd1yCGHRM+ePcuc269fvxg+fHgqXBo/fnyUlJSkXVFWEQcffHD069cv7rzzztSxF198MWvh2Je3SPrSV8OuTHz1yrPt18qUf0jvXmlpaZV9j7Zt2xbLli3b4Zi/M6A6q6rzalLOqRX9/lbW92H//fePM888M84888yI+OIfJ4899ljccccdsWnTppg5c2bcd999ce2112alt+XLl6deFxUV7XN/v7AzVfnzKsC+yM+rANnlvIpwrAw7C8dOOOGEjOY2btw4OnXqFP/6178i4os0+qOPPoq2bdvucV/f//7346677kr9R/PPf/5zj9f80vaB1pfPD8tESUlJ6pPXO1sLAKheJkyYsMuvvfTSS3H77bfHAQccEM8+++wuxzVv3rwyWsuaunXrxhVXXBF169aNq6++OiIinnvuud2GY5kqKSmJDz74ICIi9ttvv7TnjwEAAABVQzhWhkaNGkVeXl7arQW/8Y1vZDy/bdu2qXAsImLx4sVZCccaNWoULVu2jM8//zwiIpYsWRLFxcXlfj7Yzmz/S5vFixdnPHfp0qWxbdu2VN2sWbM97gcAqDodO3bc5demTp0aERG5ubm7HVddXHDBBXHdddfF1q1b49NPP83KmpMnT0590OjII4+MmjVrZmVdAAAAoOL27P5+CZCXlxetWrVKO/blfUkzsf3Y1atXZ6WviC8Csq/K1vMrDjnkkLR6/vz5Gc9duHBhWn3wwQdnpScAgMqWn58fDRs2jIjY49tgf+n+++9PvT7jjDOysiYAAACwZ4RjGTj00EPT6vI8g2v7sQUFBVnpKSLSbl+YzbW3D8dmzJiR8dzp06en1cIxAKAqlee+7fPmzYulS5dGROzw4aiK+Otf/xqjR4+OiC+upu/Tp88erwkAAADsOeFYBo466qi0ujy3GVy0aFFa3aBBg6z0VFJSEgsWLEjVeXl5Ua9evaysXa9evbRbR86YMSM2bdqU0dx33303rT7yyCOz0hMAQEXMnDkzevfuHRMnTtztuE2bNsVVV12VCtN69epV4T03b94cQ4YMicsuuywivrgK7e67787qh6QAAACAivPMsQycdNJJcdttt6V+WTJlypS46KKLypxXWlqaehZHRETNmjWjffv2WelpypQpsWbNmlSdrXW/1KNHj5g9e3ZERBQXF8drr70Wp5122m7nbN68OV5//fVUfcABB+wTzx8BACrH4sWLY9y4cWnHlixZkno9fPjwtK8dffTR5b4qvbS0NMaPHx/jx4+Pgw8+OHr16hVdunSJAw88MAoLC2P58uUxZcqUeOKJJ2LOnDkREdGiRYu48sord7vuBx98kHpdUlISq1evjgULFsQ777wTo0ePTn2YqqCgIAYNGhQnn3xyufoGAAAAKo9wLAPNmjWLb3/72/Gvf/0rIiJee+21WLFiReqZFLsyYcKEtKu7OnXqFHXr1s1KT/fdd19afdxxx2Vl3S/17Nkz7RkZw4cPLzMce/7559MCu1NPPTWrPQEA+5bZs2fHFVdcscuvb/+1e++9t9zhWJ06daKoqChWrVoVn3zySQwdOnS34zt37hwPP/xwmVfkl/WzV05OTnzve9+L2267Lb75zW+Wq2cAAACgcrmtYoa++unhzZs3x6233rrb8evXr4/bb7897Vjfvn13GLdq1apy3aYxImLIkCExYcKEVF1QUFDmMyxGjBgR7dq1S/254YYbdju+Y8eOcfzxx6fqt99+O0aNGrXL8StWrIi77rorVdeqVSt+9KMflfFOAAAqV+vWrWPWrFkxcuTIuPrqq6N79+7RokWLqF27duTm5kb9+vWjY8eO0adPnxg+fHi8+uqrcdBBB5Vrj1q1akWTJk2ibdu2cc4558Svf/3rmDx5cowcOVIwBgAAAF9DrhzL0DHHHBPHH398jB8/PiIixowZE/n5+fHLX/4y9ttvv7Sxn3/+eVx77bXx6aefpo4ddthhO312xcKFC+O8886Lc889N84666zo0qVL5OTk7LSHzz//PO68884YO3Zs2vGf/OQn0bRp0z18hzv6+c9/HhMmTIitW7dGRMR///d/R25ubpxxxhlp4+bOnRsDBgyI5cuXp45ddNFF0aRJk6z3BAB8fVxwwQVxwQUXVHj+cccdFytWrMhiRzuXl5cX3bt3j+7du+/ROs8//3yWOgIAAACqknCsHAYNGhTnn39+KvQaNWpUjBs3Lrp16xatWrWK4uLi+PDDD2PSpElRXFycmtewYcMYMmRI1Kix8wv1tmzZEk8//XQ8/fTTUVRUFB07doyDDjoo6tWrFzVq1IiVK1fG9OnT4/3334+SkpK0ub169drt7Yj2RIcOHeKXv/xl3Hzzzak+r7nmmnjooYeia9euUatWrfj444/jjTfeSAVoERFHHXVUDBw4sFJ6AgAAAAAA2BPCsXKoX79+PPTQQzFw4MCYPn16RESsXbs2xowZs8s5Bx98cAwbNiwOPPDAjPZYtWpVTJw4MSZOnLjbcTVr1owBAwbEZZddtssrzbKhT58+sXbt2rj77rtTAdiMGTNixowZOx1/9NFHx9133x15eXmV1hMAAAAAAEBFCcfKqUWLFvH000/HI488Ek8//XTMnz9/p+OaNGkSP/zhD+PCCy+MWrVq7XK9Zs2axSWXXBKTJ0+OWbNmpV2BtTP169eP008/Pfr27VvuB9JXVL9+/eI73/lODB48OCZNmrTD1WsRES1btoy+ffvGRRddVKlhHQAAwK6sWbMm+vfvn3Zs2LBhUa9evSrqCAAA+DoSjlVAXl5e9O/fP/r16xfTpk2LTz/9NJYuXRo5OTnRsGHD6NChQ7Rv3z6jtRo0aBA33HBDRHxx28KPPvoo5s2bF0uWLIn169dHSUlJ1KtXL4qKiqJdu3ZxyCGHVCh86t27d/Tu3bvc8750+OGHx6OPPhqLFy+Of//737F48eLYtGlTNGnSJNq0aROdOnWq8NoAAAAAAAB7i3BsD+Tk5ESnTp2yFgzl5+dHx44do2PHjllZrzI0bdo0TjnllKpuAwAAAAAAoEJqVHUDAAAAAAAAsLcIxwAAAAAAAEgM4RgAAAAAAACJIRwDAAAAAAAgMYRjAAAAAAAAJIZwDAAAAAAAgMQQjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGMIxAAAAAAAAEkM4BgAAAAAAQGIIxwAAAAAAAEgM4RgAAAAAAACJIRwDAAAAAAAgMYRjAAAAAAAAJIZwDAAAAAAAgMQQjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGMIxAAAAAAAAEkM4BgAAAAAAQGIIxwAAAAAAAEgM4RgAAAAAAACJIRwDAAAAAAAgMYRjAAAAAAAAJIZwDAAAAAAAgMQQjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBiCMcAAAAAAABIDOEYAAAAAAAAiSEcAwAAAAAAIDGEYwAAAAAAACSGcAwAAAAAAIDEEI4BAAAAAACQGMIxAAAAAAAAEkM4BgAAAAAAQGIIxwAAAAAAAEgM4RgAAAAAAACJIRwDAAAAAAAgMYRjAAAAAAAAJIZwDAAAAAAAgMQQjgEAAAAAAJAYwjEAAAAAAAASQzgGAAAAAABAYgjHAAAAAAAASAzhGAAAAAAAAIkhHAMAAAAAACAxhGMAAAAAAAAkhnAMAAAAAACAxBCOAQAAAAAAkBjCMQAAAAAAABJDOAYAAAAAAEBi5FZ1A/uChQsXxrRp02LBggWxYcOGKCgoiMaNG0ebNm2iffv2kZ+fX671li5dGrNnz47PP/881qxZEzk5OVFUVBQtW7aMTp06RWFhYSW9EwAAAAAAgH2bcKyCSkpKYvTo0fHEE0/E9OnTdzkuLy8vOnfuHP369Ytu3brtdExxcXFMnDgxxo0bF5MmTYo5c+bscr3c3Nz/Y+/O43M68/+Pv082EUE2Sy1RWqlUxdRSVaMoHbWVdixjam1TVLW2X2uZblNTo6aL0o59nX6r2ilqG9TWjtpH1b4EQ4QQS8gikeQ+vz88nMlNIneSO064X8/Hw2Ou69zXdeVzos7g7VyXWrRooRdffFH16tUr9H3czlNPPaW4uLgCzV29erWqVavm5ooAAAAAAAAAAAAKh3CsAE6dOqXhw4dr165deY7NyMjQtm3bVLdu3RzDsZSUFLVo0UKXL1926WtnZmbqhx9+0Jo1a9SzZ0+9+eab8vX1zfc9AAAAAAAAAAAAeCLCsXw6cuSI+vbtq4SEBOuaYRiKiopSrVq1FBoaqvT0dGurxdjY2Nuul5WVdUswZhiGatasqUceeUShoaHy9vbWqVOn9PPPP+vSpUuSJNM0NW/ePF24cEEfffSRvLyK9vg4wzDy9TUMwyjCagAAAAAAAAAAAAqGcCwfLl68qOjoaKdgrFWrVho1apSqVKmS45zDhw9r4cKFKl26dJ7rR0REqEuXLmrfvr1CQkJu+Tw9PV3Tp0/XF198IYfDIUlavny56tWrpx49ehTwrlzTqVMnjRs3rki/BgAAAAAAAAAAQFEjHMuHDz74QPHx8VZ/wIABGjp06G3nREREaOTIkbcdU7NmTQ0ZMkStWrW67bgSJUpo0KBBCgsL07vvvmtdnzRpkrp27So/Pz8X7gIAAAAAAAAAAMBzFe1efPeQTZs2admyZVa/ZcuWeQZjrggMDNSSJUvyDMay+8Mf/qB69epZ/cTERG3evLnQtQAAAAAAAAAAANzrCMdcNG3aNKvt6+urUaNGuWVdLy+vAp0X1q5dO6f+nj173FIPAAAAAAAAAADAvYxwzAWxsbHasmWL1W/WrJmqVq1qY0VSeHi4U//8+fM2VQIAAAAAAAAAAHD3IBxzwYoVK2SaptW/+a0tO6SkpDj1fXw4Pg4AAAAAAAAAACAvhGMu2LVrl1O/QYMGNlXyP4cOHXLqV6xY0aZKAAAAAAAAAAAA7h68buSCvXv3Wu0KFSqofPnykqQTJ05o4cKF2rhxo86cOaOrV68qODhY1atXV5MmTdSxY0eFhoa6vR7TNLVs2TKna48//rjbv052Bw8e1NChQ7V//35rC8egoCBVqVJFDRs2VKtWrVSrVq0irQEAAAAAAAAAAKCwCMfykJSUpHPnzln9qlWrKjMzU1OmTNHkyZOVmZnpND41NVVxcXHauHGjJk2apOjoaA0cOFCGYbitpsWLFys2Ntbqh4eHq3bt2m5bPycHDhzQgQMHnK4lJyfr1KlT2rJliyZNmqQnn3xSb731lqpVq1aktQAAAAAAAAAAABQU4VgeEhMTnfphYWEaPXq0vv/++zznpqamauLEiTp48KA+/fRTt5wLlpCQoA8//NDp2quvvurW8K2gfvrpJ/3+97/XRx99pObNmxdqreJwP8VN9u+Jnd8fb29vhYWF3XKNnzMAd5vi8FzlmQrgXsJzFQDci+cqALgXz1VkZ5imadpdRHG2f/9+Pffcc1a/dOnSSkpKkiT5+/urT58+atu2rapVq6aMjAwdOnRI33zzjZYsWaLs39qXXnpJb775ZqFqycrK0ksvvaTNmzdb1xo1aqS5c+cWyS+ep556SpmZmWrZsqWeeOIJPfTQQwoNDZWfn58SExN14MABrVmzRosWLdK1a9eseSVKlNA//vEP1a1b1+01AQAAAAAAAAAAFAbhWB527NihF1544ZbrQUFBmjNnjiIjI3Oct2TJEo0YMUIOh0PS9SR68eLFhTqXa+zYsZo7d65TDYsWLVKlSpUKvObtbN68WQ0bNszzjbcTJ05o0KBBOnz4sHWtcuXKWrlypfz8/IqkNgAAAAAAAAAAgIJgW8U85BbuvP/++7kGY5L07LPPau/evVaYZZqmZs2apfHjxxeojjlz5jgFY76+vvr000+LLBiTpMaNG7s0rlq1apozZ446d+6s06dPS5Li4uL07bff5hgsuuLm7SxxPWAtW7asJOny5cuyK9dOTk7WyJEjna6NGzdOgYGBttQDAAVVHJ6rPFMB3Et4rgKAe/FcBQD34rl69wsKCnLbWoRjeShVqtQt1x544AG1bt06z7n9+vXTV199pYyMDEnShg0b5HA45OXlla8aFi9erHHjxll9wzA0btw4PfHEE/lapyiFhobq//2//6dhw4ZZ11asWFHgcIwXGm/PNE3bvkdZWVk6f/78Ldf4OQNwN7PrucozFcC9iucqALgXz1UAcC+eq8hfSuOBcgrHWrRo4dLcsLAwRUVFWf3Lly8rJiYmX19/7dq1+tOf/uT0i+Ptt99W+/bt87XOndC6dWunhHvXrl26evWqjRUBAAAAAAAAAAA4IxzLQ2hoqHx9fZ2u1axZ0+X5ERERTv2zZ8+6PHfz5s0aMmSIMjMzrWuDBw8u8NtYRc3Hx0d16tSx+pmZmTp37pyNFQEAAAAAAAAAADgjHMuDr6+vwsPDna7d2JfUFTePvXz5skvzdu/erYEDB+ratWvWtRdffFEDBw50+WvbITQ01Kl/6dIlmyoBAAAAAAAAAAC4FeGYCx588EGnfvbAKi83jy1RokSecw4dOqSXX35Zqamp1rWuXbtqxIgRLn9du9y8jaIr9wsAAAAAAAAAAHCnEI654LHHHnPq52drxPj4eKd+cHDwbcefOHFCL730khITE61r7dq105///GeXv6adYmNjnfohISE2VQIAAAAAAAAAAHArwjEXtGrVSoZhWP2dO3e6NM80Te3atcvqe3t7q1atWrmOP3v2rPr27auEhATr2lNPPaXx48fLy6v4/1SdOXNGR44csfqhoaEqX768jRUBAAAAAAAAAAA4K/6JSzFQsWJF1atXz+qvW7dOFy9ezHPexo0bdfr0aasfFRWlwMDAHMdeunRJffv2VVxcnHWtcePG+uyzz+Tj41OI6u+cv//97zJN0+o3adLEKVQEAAAAAAAAAACwG+GYi1577TWrnZ6erjFjxtx2fEpKij744AOnaz179sxxbHJysqKjo3X06FHr2qOPPqq///3v8vPzK0TV/7N161Y99NBD1o/capGun5N27NixfK3/3Xff6ZtvvrH6hmGoV69eBa4XAAAAAAAAAACgKBCOuahx48Zq3ry51V+xYoVGjBihpKSkW8aePHlSffv21fHjx61rderUUZs2bW4Ze+3aNQ0cOFB79+61rtWuXVvTp09XQECAe2/CRWlpaWrfvr2GDRumjRs3KjMzM9exCQkJ+vOf/6zRo0c7Xe/UqZPq1KlT1KUCAAAAAAAAAADky92xX18xMX78eHXr1s0KvRYvXqy1a9eqadOmCg8PV0ZGhg4fPqwtW7YoIyPDmhcSEqKJEyfmeG7Yrl27tHXrVqdrBw4cUKNGjfJVW8OGDTV37twC3FXOsrKytHz5ci1fvlyBgYGKjIxUjRo1VLZsWfn6+urKlSs6ePCgdu3a5XSvktSgQQO9//77bqsFAAAAAAAAAADAXQjH8qFs2bKaMWOGXn/9de3bt0+SlJSUpBUrVuQ6p0aNGpo6daoqVaqU4+fZz+i6weFw5Lu2gsxxVXJysrZv367t27fnObZ79+4aOXKk27aDBAAAAAAAAAAAcCfCsXyqUqWKFixYoFmzZmnBggWKi4vLcVz58uXVu3dv9ejRQ/7+/ne4ysLx9/fXgAEDtHXrVu3bt0/Xrl277fiAgAC1atVKvXr1YitFAAAAAAAAAABQrBGOFYCvr6/69++vfv36ac+ePTp+/LgSEhJkGIZCQkIUGRmpWrVqubRWo0aNdOjQoSKuOH9fx8/PT0OHDpUkZWZm6vjx4zp58qTi4+OVkpKizMxMlS5dWmXKlFHNmjX10EMPydvbuyjLBwAAAAAAAAAAcAvCsUIwDENRUVGKioqyu5Qi4+Pjo5o1a6pmzZp2lwIAAAAAAAAAAFBoXnYXAAAAAAAAAAAAANwphGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEY9uF59AAACAASURBVAAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8BiEYwAAAAAAAAAAAPAYhGMAAAAAAAAAAADwGIRjAAAAAAAAAAAA8Bg+dheQm5SUFB06dEiJiYlKSkqSaZrq1KmT3WXl6MyZM9qzZ49Onz6t1NRUlShRQmFhYapevbpq1aolPz+/Aq179uxZ7d69W2fPntXVq1dVoUIFVa9eXXXq1HHzHbguMTFRO3fuVHx8vJKTk1W+fHlVqVJFjz76qLy9vW2rCwAAAAAAAAAAwBXFKhzLzMzUkiVLNH/+fO3fv18Oh8Pp85zCsb1792rp0qWSJD8/Pw0fPvyO1OpwOLRkyRLNmzdP+/bty3Wcr6+vHn30UfXr109NmzZ1ae3du3drwoQJ2rx58y3fA0kKDw9Xz5491bNnTxmGUeB7yI/jx4/r448/1oYNG5SRkXHL5+XKlVO3bt3Uv3//AoeBAAAAAAAAAAAARa3YhGMxMTEaPHiwjh07JkkyTdPp89xCoPvvv18LFy5UcnKyJKlJkyZ6/PHHi7TWU6dOafjw4dq1a1eeYzMyMrRt2zbVrVvXpXBs+vTpmjBhgjIzM3Mdc/LkSX3wwQdat26dJkyYoKCgoHzVn1/ff/+93nvvPaWmpuY6JiEhQZ9//rnWrVunSZMmqUqVKkVaEwAAAAAAAAAAQEEUi3Bs79696tWrl65evZpjKHbztewCAwP13HPPad68eTIMQ8uXLy/ScOzIkSPq27evEhISnGqMiopSrVq1FBoaqvT0dGurxdjYWJfX/vrrr/XRRx85XYuMjFT9+vVVsmRJHT16VD/99JMVnG3evFmvv/66Zs6cKV9fX/fc4E1+/PFHjRw50ukNtvvvv1+NGjVSUFCQTp48qfXr1ystLU2StH//fvXv318LFixQYGBgkdQEAAAAAAAAAABQULaHY5cvX9Yrr7yi1NRUGYYhwzDUpk0btW3bVuHh4XruueeUlZV12zXat2+vefPmSZJ+/vnnIqv14sWLio6OdgrGWrVqpVGjRuX6ptThw4e1cOFClS5d+rZrHzhwQGPGjLH6fn5+Gjt2rDp06OA0LjY2VgMHDtThw4clSVu3btXEiROLZDvJc+fOaejQoVYwZhiGRowYod69e8vLy8sad/HiRQ0ePFjbtm2TdP0twHfeeUeffPKJ22sCAAAAAAAAAAAoDK+8hxSt6dOnKyEhQYZhyN/fX9OmTdMnn3yiVq1aKSIiwqU1oqKiFBQUJNM0debMGcXHxxdJrR988IHT2gMGDNAXX3xx2y0EIyIiNHLkSPXv3/+2a3/66adOWymOGTPmlmBMkqpWraq5c+cqNDTUujZv3jydO3cuP7fiksmTJyslJcXqv/baa+rbt69TMCZJISEhmjFjhh544AHr2ooVK7R//3631wQAAAAAAAAAAFAYtoZjDodDCxYssPqjR4926VyunNSuXdtqHz16tNC13WzTpk1atmyZ1W/ZsqWGDh3qlrX379+vH3/80eo3bNhQnTp1ynV8SEiIhg0bZvXT0tI0a9Yst9Ryw/nz5/Xtt99a/fDwcPXr1y/X8SVKlNDbb79t9U3T1OTJk91aEwAAAAAAAAAAQGHZGo79+uuvSkpKkmEYuu+++9S5c+cCr1WtWjWrferUKXeU52TatGlW29fXV6NGjXLb2itXrnTq//GPf8xzTocOHVSmTJlc1yistWvXKiMjw+p37do1z3PNGjdurBo1alj9H3/8UampqW6tCwAAAAAAAAAAoDBsDceOHz9utR9//HEZhlHgtbKf6ZWcnFyoum4WGxurLVu2WP1mzZqpatWqblt/3bp1VtvX11ctW7bMc06JEiXUvHlzq3/mzBm3bmOYvSZJeuaZZ1ya17p1a6udnp5epGfAAQAAAAAAAAAA5Jet4djFixetdsWKFQu1VvZzsLKysgq11s1WrFgh0zStfrt27dy29uXLl3XkyBGr//DDD6tEiRIuzX300Ued+tu3b3dbXf/5z3+sdlhYmMthYFHWBAAAAAAAAAAAUFi2hmM+Pj5Wu7CB1qVLl6x29u0G3WHXrl1O/QYNGrht7ZvPR4uMjHR57sMPP3zbtQrq7NmzSkpKKlY1AQAAAAAAAAAAuINP3kOKTnBwsNWOj48v1Fp79uyx2iEhIYVa62Z79+612hUqVFD58uUlSSdOnNDChQu1ceNGnTlzRlevXlVwcLCqV6+uJk2aqGPHjgoNDb3t2seOHXPqV6pUyeW67rvvPqd+9m0qC6MwNYWFhcnX19c6r8xdNQEAAAAAAAAAALiDreFY9erVrXb2bfzyKzY2VgcOHLD6UVFRhaoru6SkJJ07d87qV61aVZmZmZoyZYomT56szMxMp/GpqamKi4vTxo0bNWnSJEVHR2vgwIG5nqd29uxZp35+tpcsV66cfHx8rBoKGzC6oybDMFShQgWdOnUqx7UAAAAAAAAAAADsZGs4VqdOHQUFBSkxMVFxcXH66aef9OSTT+Z7nc8++0wOh0OSFB4eXujzy7JLTEx06oeFhWn06NH6/vvv85ybmpqqiRMn6uDBg/r000+dtpHMPia7UqVKuVybl5eX/P39lZycnONaBXXzOgEBAfman/0eMjMzlZ6e7vI5ajfkFiZ6suzfEzu/P97e3goLC7vlGj9nAO42xeG5yjMVwL2E5yoAuBfPVQBwL56ryM7WcMwwDLVv315ffvmlTNPU+++/r3/+858KCgpyeY0ZM2Zo2bJl1npdu3Z1a43Zz96SpJ9//tm65u/vrz59+qht27aqVq2aMjIydOjQIX3zzTdasmSJTNOUJK1evVqffPKJ3nzzzVvWvzmI8vPzy1d9JUqUKPJwLL/B1s3jU1JS8r1Gfv4b8ERly5a17WsHBQXpyy+/tO3rA0BRsOu5yjMVwL2K5yoAuBfPVQBwL56r8LK7gIEDByogIECGYejUqVPq0aOH0/lhuYmNjdXw4cP18ccfW6lqSEiIunfv7tb6bg6KbgRjQUFB+vrrrzV06FA99NBD8vf3V+nSpdWgQQONHz9e48ePl5fX/769s2bN0sGDB29ZPz093anv6+ubr/qyh2lpaWn5mpsbd9aU03oAAAAAAAAAAAB2sfXNMel6oPWXv/xFw4cPl2EYiomJUdeuXVW3bl3Vr1/fevtKkv75z38qNjZW27Zt0+7du+VwOKzPvb299dFHH+V7C8C85PYm1/vvv6/IyMhc5z377LPau3ev5s6dK0kyTVOzZs3S+PHjncbd/EZVRkZGvuq7du2a1fb398/X3Ny4s6ac1nPFzdtZ4vqbkTf+RcPly5edfm3cScnJyRo5cqTTtXHjxikwMNCWegCgoIrDc5VnKoB7Cc9VAHAvnqsA4F48V+9+7txxzvZwTJLatm2rixcv6q9//asVeP3666/69ddfrTGmaertt9926t94Y8zHx0fvvvuuGjdu7PbacjoD7IEHHlDr1q3znNuvXz999dVXVri0YcMGORwOpzfKbg7zbg6W8pL9rSx3BYM3r5PfN79uHp+fc9RusCv4uVuYpmnb9ygrK0vnz5+/5Ro/ZwDuZnY9V3mmArhX8VwFAPfiuQoA7sVzFbZvq3hDjx49NGfOHFWuXFmSczhiGEaOB9KZpqlKlSpp5syZ6tKlS5HUlVOw06JFC5fmhoWFKSoqyupfvnxZMTExTmNuDqJunB/mCofDoatXr+a6VkHdvE5KSkq+5mcf7+PjU6A3xwAAAAAAAAAAAIpCsQnHJKlhw4ZavXq1Pv74Yz355JMqVaqUleBm/+Hr66vHHntMY8aM0erVq9WoUaMiqyk0NPSWM7dq1qzp8vyIiAin/tmzZ536FSpUuO3nt5OQkKCsrCyrX7FiRZfn3k5hajJN02n8zWsBAAAAAAAAAADYqVhsq5idl5eX2rVrp3bt2snhcOjkyZO6dOmSrly5opIlSyo4OFjh4eF37G0kX19fhYeH6+jRo9a1G/uSuuLmsZcvX3bqP/DAA079uLg4l9c+c+aMU79GjRouz72dwtR0/vx5pzPK3FUTAAAAAAAAAACAOxS7cCw7Ly8v3X///br//vttrePBBx90Csfycy7YzWNvDvVuDqIOHDjg8tr79u1z6rsriCpfvrxKly6tpKSkfNe0f//+IqkJAAAAAAAAAADAHYrVtorF1WOPPebUz882g/Hx8U794OBgp36ZMmWctmk8cOCA0tLSXFr7l19+ceo3bNjQ5bryUr9+fat94cIFnTx50qV5O3fuLLKaAAAAAAAAAAAACotwzAWtWrWSYRhW/+YAKDemaWrXrl1W39vbW7Vq1bplXMuWLa12RkaG1q1bl+fa6enp+vHHH63+fffdp4cfftilulyRvSZJWrlypUvzVq1aZbVLlCihJk2auK0mAAAAAAAAAACAwioW4VhycrL1w475ealYsaLq1atn9detW6eLFy/mOW/jxo06ffq01Y+KilJgYOAt41q3bu3U/+qrr/Jce+nSpbpy5YrVf+aZZ/Kckx9PPfWUfH19rf63337rdJZYTjZv3qzjx49b/WbNmikgIMCtdQEAAAAofhITs9St2/92zejWLV6JiVk2VgQAAAAAubM9HFu4cKEaNmyohg0bqmfPngVao2fPntYay5cvd3OF17322mtWOz09XWPGjLnt+JSUFH3wwQe31JmThx9+WM2bN7f627dv1+LFi3Nd++LFi/rkk0+svr+/v1588cXb1rN161Y99NBD1o+8vtdhYWHq2rWr1T958qSmTZuW6/j09HT95S9/sfqGYeiVV1657dcAAAAAAAAAAAC402wPxxYtWiTTNCVJffr0KdAaffr0kWmaMk1T//znP91Y3f80btzYKcBasWKFRowYoaSkpFvGnjx5Un379nV6i6pOnTpq06ZNrusPGTJEPj4+Vv/tt9/WsmXLbhkXGxur3r1768KFC9a1Xr16qXz58vm9pTz1799fpUqVsvqTJk3S7Nmz5XA4nMZdvHhR0dHRiomJsa61bdvWrds8AgAAAAAAAAAAuINP3kOKTnJysn755RdJkp+fn1q1alWgdZ5++mm9++67SktL044dO3T16lWVLFnSnaVKksaPH69u3bpZodfixYu1du1aNW3aVOHh4crIyNDhw4e1ZcsWpy0IQ0JCNHHiRHl55Z5FRkZG6q233tJ7770nSbp27ZqGDx+uGTNmqEGDBvL399fRo0f1008/KTMz05r32GOP6fXXX3f7vUpShQoV9Mknn+iVV16Rw+GQaZoaN26cvv76az3++OMKCgrSiRMntH79eqWlpVnzHnzwQb3//vtFUhMAAAAAAAAAAEBh2BqO7d+/X5mZmTIMQ7Vr13Z6Syk/AgIC9Mgjj2jHjh3KzMzUvn371KBBAzdXK5UtW1YzZszQ66+/rn379kmSkpKStGLFilzn1KhRQ1OnTlWlSpXyXL979+5KSkrSZ599ZgVgBw4c0IEDB3Ic//jjj+uzzz5zOhvM3Zo3b66//vWveu+993T16lVJ0n//+1/997//zXF8ZGSkPv/88xzPVgMAoLhITMzSgAHntWpVkKTrZ+NMmRKmoCBvmysDAAAAAABAUbN1W8XsActDDz1UqLUiIiJyXNfdqlSpogULFmjYsGGqXLlyruPKly+vN954Q4sWLVJ4eLjL6/fr109fffWVnnjiiVzfNKtatapGjx6tOXPmKCgoKN/3kF+dOnXSwoUL1apVq1yDuHLlyunVV1/VN998oypVqhR5TQAAAAAAAAAAAAVh65tjiYmJVjs4OLhQa2UPiS5fvlyotfLi6+ur/v37q1+/ftqzZ4+OHz+uhIQEGYahkJAQRUZGqlatWgVev27dupo9e7bOnj2rX3/9VWfPnlVaWprKly+v6tWrKyoqKt9rNmrUSIcOHSpwTTVq1NAXX3yhS5cuaefOnYqPj1dKSorCwsJUtWpV1atXT97e/Gt7AAAAAAAAAABQvNkajhmGYbWzn9FVENnP4crKyirUWq4yDENRUVEFCqtcUaFCBf3ud78rkrULKjg4WC1btrS7DAAAAAAAAAAAgAKxdVvF7G+LnTt3rlBrnT171mrfia0GAQAAAAAAAAAAcPexNRyrWLGi1d62bVuh1tq6davVLleuXKHWAgAAAAAAAAAAwL3J1nCsXr168vG5vrNjfHy81qxZU6B11qxZo/j4eEmSl5eX6tev77YaAQAAAAAAAAAAcO+wNRwLCAhQgwYNZJqmTNPU2LFjdf78+XytkZCQoLFjx0q6fgbYb37zG5UpU6YoygUAAAAAAAAAAMBdztZwTJL69+8v6Xqwdfr0afXq1UsHDx50ae7BgwfVq1cvnTlzxrrWr1+/IqkTAAAAAAAAAAAAdz8fuwto3LixnnrqKa1bt06GYejYsWPq3Lmznn76abVr105169Z1OkPs3Llz2r17t5YuXaq1a9cqKytL0vVwrVmzZmrWrJldtwIAAAAAAAAAAIBizvZwTJL+9re/qVu3boqJiZFhGMrMzNTKlSu1cuVKSdeDL39/f6Wlpck0TWueaZoyDEOmaapmzZr6+OOP7boFAAAAAAAAAAAA3AVs31ZRkkqVKqWvvvpKTZs2tQIvSdZZZA6HQ6mpqXI4HNa1G0zTVPPmzfXll1+qVKlSdt0CAAAAAAAAAAAA7gLFIhyTpDJlymjatGkaN26cHnzwQacALCemaSoiIkIfffSRpkyZorJly96hSgEAAAAAAAAAAHC3KhbbKt5gGIY6deqkTp066eDBg9q2bZuOHDmixMREpaSkKDAwUEFBQYqIiNBjjz2miIgIu0sGAAAAAAAAAADAXaRYhWPZ1apVS7Vq1bK7DKBYSEzM0oAB57VqVZAkqVu3eE2ZEqagIG+bKwMAAAAAAAAA4O5SbLZVBAAAAAAAAAAAAIoa4RgAAAAAAAAAAAA8BuEYAAAAAAAAAAAAPAbhGAAAAAAAAAAAADyGj90FZJeamqrVq1frl19+0bFjx3T58mWlpqbKNE2X1zAMQ2vWrCnCKgEAAAAAAAAAAHC3KhbhmMPh0LRp0zRt2jRdvXrV6TNXgzHDMGSapgzDKIoSAQAAAAAAAAAAcA+wPRzLysrSa6+9pvXr11tB2I2gKz/yOx4AAAAAAAAAAACex/ZwbObMmVq3bp2k/4ViPj4+ioiIUOXKlRUQEMDbYAAAAAAAAAAAAHALW8OxjIwMTZ8+3QrFvL29NXDgQPXo0UNlypSxszQAAAAAAADAFomJWRow4LxWrQqSJHXrFq8pU8IUFORtc2UAANwbbA3Hdu7cqaSkJBmGIcMwNG7cOLVv397OkgAAAAAAAAAAAHAP87Lzix8/ftxqR0ZGEowBAAAAAAAAAACgSNkajl25csVqN2zY0MZKAAAAAAAAAAAA4AlsDceCg4OtdunSpW2sBAAAAAAAAAAAAJ7A1nCsUqVKVvvixYs2VgIAAAAAAAAAAABPYGs41rBhQ5UtW1aStHPnTjtLAQAAAAAAAAAAgAewNRzz8/NTly5dZJqmDh06pB07dthZDgAAAAAAAAAAuMckJmapW7d4q9+tW7wSE7NsrAh2szUck6TXX39dERERMk1To0aN0oULF+wuCQAAAAAAAAAAAPco28MxPz8/TZ06VREREYqNjdXzzz+vtWvX2l0WAAAAAAAAAAAA7kE+dhfw+eefS5KaNGmi48eP6+zZsxo0aJAqVqyo+vXrq0KFCipZsmS+1hw0aFBRlAoAAAAAAAAAAIC7XLEIxwzDsPqGYcg0TZ05c0bLly8v0JqEYwAAAAAAAAAAAMiJ7dsq5sQwDKfAzFWmaRZBNQAAAAAAAAAAALhX2P7mmESoBQAAAAAAAAAAgDvD9nDs4MGDdpcAAAAAAAAAAAAAD1Est1UEAAAAAAAAAAAAigLhGAAAAAAAAAAAADwG4RgAAAAAAAAAAAA8BuEYAAAAAAAAAAAAPAbhGAAAAAAAAAAAADwG4RgAAAAAAAAAAAA8ho/dBdzs559/1vr167V7927Fx8frypUrSk9Pd3m+YRjav39/EVYIAAAAAAAAAACAu1WxCcd2796t0aNH6+jRo9Y10zRtrAgAAAAAAAAAAAD3mmIRjq1Zs0ZDhw5VZmamUyBmGIbVzu36zZ8BKDopKSlavHixYmJilJWVJW9vb/n7++vy5csqU6aM3eUBAAAAAAAAAJAn28Ox2NhYjRgxQhkZGTIMQ97e3mratKkiIiI0ffp0maYpwzA0aNAgpaSk6MKFC9q7d6+OHTsm6XpQFhAQoO7du6tkyZI23w1wbzpy5IhmzZql+fPn68qVK7d83qhRI3Xu3FnR0dGKioqyoUIAAAAAAAAAAFxjezg2bdo0paSkSJICAwM1ffp0/eY3v5EkzZw5U1lZWZKkQYMGOc07fPiw/v73v2vlypW6evWqNmzYoBkzZui+++67szcA3OOmTJmit956Sw6HI9cxaWlp+vLLL/Xll19qyJAheuutt+Tl5XUHqwQAAAAAAAAAwDW2/u11RkaGli5dKsMwZBiG3nrrLSsYy0tERIQmTJigMWPGSJKOHTum6OhopaamFmXJgEf58MMPNXr0aCsY8/f3V9euXfXoo4+qYcOGqlevnsLDw1W6dGlrzoQJEzRs2DC2OwUAAAAAAAAAFEu2hmN79+5VWlqaTNNUSEiInn322Xyv0aVLF73yyisyTVPHjh3T1KlTi6BSwPMsWLBAH374odUfNGiQ9u3bp/Hjx6t8+fIKCQlRuXLlFBkZqS1btujdd9+13habN2+eJkyYYFfpAAAAAAAAAADkytZwLPu5YfXr15dhGLmOvbG9Yk4GDBigsmXLyjRNLVy4kDdWgEK6du2a3nvvPas/ZswYvf/++woODs5xfKlSpTR48GBNmzbNuva3v/1NiYmJRV0qAAAAAAAAAAD5Yms4duXKFasdHh5+y+fZzyy6du1aruv4+fnpt7/9rSTp/Pnz2rlzpxurBDzPsmXLdPbsWUlSmzZt9Oqrr7o07/nnn1d0dLSk6+eQ/d///V+R1QgAAAAAAAAAQEHYGo5lD7xKlix5y+elSpWy2hcvXrztWlWqVLHap06dckN1gOeaNWuW1e7fv3++5r788stWe/bs2dZ5ZQAAAAAAAAAAFAe2hmPZw6+0tLRbPi9TpozVjouLc3nd8+fPF64wwIMlJSVp06ZNkqQHH3xQTZs2zdf8mjVrWnOOHTumo0ePur1GAAAAAAAAAAAKytZwrGLFilY7p7OJqlevbrV37dp127ViYmKsdvbtGAHkT/a3NOvUqXPbswBzU6dOHat94cIFt9QFAAAAAAAAAIA72JoiPfDAA1b7+PHjt3z+8MMPW+2lS5fmuk58fLz+/e9/W/0KFSq4qULA82RkZFhtX1/fAq2Rfd7tzgsEAAAAAAAAAOBOszUcq169uoKDg2Wapg4ePHjL2UStW7e22jExMRo3bpxM03Qac/HiRQ0aNMj6C3jDMNSgQYOiLx64R5UtW9Zq52c70+xOnz5ttYOCggpdEwAAAAAAuLPS0tL03Xff6T//+Y82bdqkf//739q8ebOGDh2qjRs33vJ3dAAA3E187C7g8ccf17/+9S+lpKTo119/1aOPPmp9VqtWLTVs2FA7duyQJM2dO1fr169XkyZNVKZMGcXFxWn9+vVKSUmRdD0Ya9GihcqXL2/LvQD3grCwMFWpUkWnTp3Spk2bdOLECVWrVs3l+VeuXNHy5cslSYGBgU5viAIAAAAAgOItMTFREyZM0D/+8Q9dunTpls8XLVqkRYsWKSIiQgMGDFCvXr044gQAcNex/f+5brwdZppmjlsnvvvuuypZsqTVP3HihObPn6+pU6dq2bJlSk5Otj4rU6aMRo0aVfRFA/cwwzDUp08fSdd/Xc6ePTtf87/++mulpqZKkrp166ZSpUq5u0QAAAAAAFAETpw4odatW2vixIlOwZhhGPLxcf439ocPH9awYcMUHR2ttLS0O10qAACFYvubY82bN7cCrdKlS9/y+YMPPqjp06dr8ODBOn/+fI5rmKapihUr6osvvlCVKlWKtF7AKiJjCAAAIABJREFUE/To0UMffvihMjIyNHPmTD333HOqW7dunvNOnDihjz/+2Oq/+OKLRVkmAAAAAABwk/j4eHXo0EGnTp2SJPn5+aldu3Y6ffq0ypYtK8Mw5HA49Oyzz2r+/PnatGmTJGnx4sXKyMjQnDlz5O3tbectAADgMtvfHPP391fv3r3Vu3dvPf/88zmOqV+/vlatWqXhw4erfv36Cg4Olre3t8qUKaMGDRpo5MiR+te//qXatWvf4eqBe1P58uX1wgsvSJJSUlLUuXNnbdu27bZzDh8+rI4dOyohIUHS9bdCIyMji7xWAAAAAABQOKZpqnfv3lYwVrNmTW3dulWffvqpgoKCZBiGJMnLy0sdOnTQsmXLNH/+fAUEBEiSli9f7vSPZQEAKO5sf3PMVaVKldLLL7+sl19+2e5SAI8wduxYHThwQFu3btWFCxfUpk0b/e53v9Mf/vAHmaYpwzBkmqYSExP12muvaeXKlcrIyJB0/TfRkydPtvkOAAAAAACAKzZt2qTt27dLkqpWraolS5aoQoUKunLlSq5zWrdurXnz5qlbt27KysrSlClT9NprrzkdjwIAQHF114RjAO4sf39/ff311+rZs6c2btwo0zS1atUqrVq1Sl5eXvLx8VFWVpaysrKc5j3yyCNasGCBgoKCbKocAAAAAADkx8yZM632O++8owoVKrg076mnnlLnzp21YMECJSYmatGiRfrjH/9YVGUCAOA2tm+rCKD4Klu2rL799lu98847qly5snXd4XDo2rVrTsFYWFiYhg0bpuXLl+u+++6zo1wAAAAAAJBP586d07JlyyRJ5cqVU4cOHfI1P/t547NmzXJrbQAAFBXCMQC3VaJECQ0ZMkS//PKLvvzySz399NMKDAyUv7+/AgMDFRYWps8++0x79uzRW2+9pdKlS9tdMgAAAAAAcNGePXuUmZkpSfr9738vPz+/fM1v0KCBatasKUnatWuXHA6H22sEAMDd2FYRgEt8fHzUtm1b/fa3v1X//v2dPuvYsaNKlChhU2UAAAAAAKCgEhMTrXb2XWNcZRiGKleurCNHjsjhcCg5OVllypRxZ4kAALgd4RgAAAAA4J7hcDiUkZGhrKws65xcAEDufHz+99eDBX1m3njzTJK8vb0LXRMAAEWtSMOxyMhIq20Yhvbv33/bMe6Q29cBAAAAANybTNPUzp07NWXKFK1du9ZpS6969eqpe/fu6tu3r7XtFwDgf8LCwqz23r178z3/2rVrOnTokCTJ399fAQEBbqsNAICiUqRnjpmmaf3vjXZeY9zxAwAAAADgGfbs2aOnn35aTz/9tL777rtbzrq5fPmypkyZokaNGqlLly46deqUTZUCQPFUr149BQcHS5KWLFmi8+fP52v+8uXLlZCQIElq1aqVDMNwe40AALhbkYZjklwKqwi0AAAAAAD5tWHDBrVr1047d+60rvn4+KhcuXKqWLGiwsLCnM7GXbt2rVq3bs1uIwCQTcmSJfXCCy9Iuv4W2Ny5c12ea5qmpk+fbvWjo6PdXh8AAEWhSLdVHDRokFvGAAAAAACQ3S+//KIePXooNTVVklSrVi29+OKLWr16tdN5Nx9++KGWLl2qyZMn6/Tp0zpz5oy6dOmi1atXq3LlynaVDwDFSt++ffX5559Lkv72t7/pscceU9OmTW87xzRNjR07Vlu2bJEk1axZM885AAAUF4RjAAAAAIC7SlZWlqKjo61grG3btpo+fboyMjK0du1ap7FBQUF69dVX1a1bN3Xt2lW7du3SmTNnNGTIEH377bd2lA8AxU716tX18ssva/r06bp27Zq6du2q8ePHq02bNjmOv3DhgsaNG6eZM2dKkgzD0J///Ge2VAQA3DWKfFtFAAAAAADc6YcfftDx48clXT8rZ8aMGSpZsuRt54SFhembb75RpUqVJF3fYjEmJqbIawWAu8UHH3yg1q1bS5LS09M1ePBgPfHEEzpy5IguXLigxMREnTt3TsOGDdMjjzxiBWM35j7zzDN2lQ4AQL4RjgEAAAAA7iqzZs2y2m+++ab8/f1dmhcWFqYBAwZY/dmzZ7u9NgC4W/n4+Gju3Lnq3r27de38+fM6duyYduzYoa1bt+qXX37RwoULlZ6ebs2ZOHGi07MVAIC7AeEYAAAAAOCuER8fb22dGB4erpYtW+Zr/h//+EcrTJs/f76ysrLcXiMA3K38/Pz0+eef6/vvv9ezzz7rdIZjdkFBQRo4cKC2bNmiHj163OEqAQAovCI9c8wVH330kZ5//nnVqFHD7lIAAAAAAMXcf//7X5mmKUl65plncv2L29yEhIToiSee0Lp165SYmKjExESFhoYWRakAcFcyDENNmzZV06ZNdeTIEfXu3Vvp6enKysqSj4+PBg8erM6dO+e5nS0AAMWZ7eHYjBkzNHPmTNWtW1e///3v1aZNGwUGBtpdFgAAAACgGEpKSrLaZcuWLdAaQUFBTusRjgFAzipUqKBq1ao5XevYsSPBGADgrldstlX89ddf9c4776hp06YaMWKEtm7dandJAAAAuAu9+uqrCgkJUUhIiF599VW7ywHgZqVKlbLaqampBVojJSXFagcEBBS6JgAAAAB3F9vfHJMk0zRlGIZM09TVq1e1ZMkSLVmyRJUqVdLzzz+v5557TpUqVbK7TAAAAACAzSpXrmy1N2zYYP150lUpKSnavHmzJKlkyZIKDg52e40AAAAAijfb3xxbsWKFXnrpJZUrV87pummaiouL0+eff65WrVqpb9++WrZsma5du2ZTpQAAAAAAu1WrVk0NGjSQJO3bty/fu44sXLhQV65ckSR16tRJvr6+bq8RAAAAQPFmezhWo0YNvfHGG9qwYYOmTp2q1q1b3/KHE4fDoS1btuiNN95QkyZN9N5772n37t02VQwAAAAAsNNLL71ktSdMmCCHw+HSvKtXr2ry5Mk5rgMAAADAc9gejt3g5eWlZs2a6bPPPtPGjRv19ttvq3bt2jJN0xpjmqaSkpK0YMECdevWTe3bt9fs2bN14cIFGysHAADIP9M0tX37du3bt087duzQ1q1btWPHDo0aNUpbt251+j0QAMBZx44dFRYWJklavXq1/vSnP+UZkKWlpSk6OloHDx6UJNWrV0/16tUr8loBAACAG2JiYnT06FHt379fe/fu1cGDB/XNN98oOTnZ7tI8TrEJx7IrU6aMXnjhBX333XdaunSp+vTpo9DQUKcxpmkqJiZG48ePV7NmzTRw4ECtWbNGWVlZNlUNAACQt8zMTM2ZM0dNmzZVly5ddOrUKV24cEGJiYm6cOGC5s+frzZt2qhZs2b6xz/+Uax+b3Pt2jXNnz9fffr0Ub169RQeHq7KlSurbt266t69u2bNmuXyb+g7dOigkJAQhYSEaNy4cS7NGTdunDWnQ4cOt3x+47P58+db1+bPn29dv/nHV1995dqNAyh2/P39NWnSJHl5Xf8j7dSpU9W9e3dt3779ln9ckJWVpVWrVql9+/b617/+JUkKDAzUp59+esfrBgAAgOdxOBz6/vvv1bFjR7Vq1UoxMTGKjY1VXFycTpw4oTfffFO1a9fWyJEjFRMTY3e5HsPH7gLyUrNmTY0cOVJvvPGGfvrpJy1cuFAbNmxQRkaGNSYzM1Pr16/X+vXrFRISop9//tnGigEAAHKWnJysl156ST/88EOeY/fu3avBgwdr1apVmjZtmgICAu5Ahblbv369hg0bphMnTtzyWWxsrGJjY7Vq1SqNHz9eY8eO1fPPP29DlQA8SevWrTVhwgQNGTJEDodDP/zwg3744QcFBgYqNDRUPj4+ysjI0JNPPqm4uDhrXqlSpTRv3jzVqVPHxuoBAADgCZKSkvTiiy9q7dq1eY6bNm2aZs+erQkTJqh79+53qELPVezDsRu8vb3VokULtWjx/9n79zgtq3p//H/dM8NBEARBxTRRTFFMBZPSVDxhaVkCpn1L86y71Eza7TYatctzn9x9dmru+n0qzxqUQKbmNg9ppiImpgGKkodQVBCJ82GG+/fHbO4YDjIz3HOA+/l8PObxuNY111rX+wa9HOd1rbUOz3vvvZe77rorEyZMyLRp00rXFIvFzJ07tw2rBABYt2XLluULX/hCg5d4PvKRj2T58uXp3bt3ampqUltbm2OPPTa33357/vznPydJ7r333px00kkZM2ZMOnbs2Ca1//rXv865556b2tra0rnu3btnt912S4cOHfLyyy9nzpw5SZJ33nknZ599dmbNmpXzzjuvVes84ogjkiTTpk3LrFmzkiTbb7999txzz3Vev/3227dabUDLOPnkk9O7d+9ccMEFpefQwoUL1zuLdZdddsnPf/7zDBw4sDXLBACgAi1atCjDhw/PM888Uzq3yy67pGPHjunRo0eqq6uzfPny7LrrrrnrrruyZMmSrFixIuedd14WL15sf9wW1i6XVdyQnj175tRTT8348eMzduzY9O3bt61LAgB4XxdddFEpGOvevXvGjRuXO++8Mx/4wAfSsWPHVFVVpWPHjjn++OPz+9//PmPHjk23bt2SJI888ki+853vtEndU6dOzVe/+tVSMNa9e/dcc801mT59en7/+9/n3nvvzQsvvJCbb745ffr0SVL/wtJ3vvOdPPLII61a669//ev8+te/zmGHHVY6d9hhh5XOr/l1+OGHt2p9QMs4+uij89xzz+UnP/nJOvcQKxQKOeqoo/LLX/4yTz31lGAMAIAWVywW85WvfKUUjG211Va57bbb8uCDD2bnnXdOjx490q1bt/Tq1Svf//73M2XKlJxxxhml/v/+7/+ehx56qK3KrwibZDiW1C81dMkll+Tss8/O66+/3tblAACs15tvvplbbrklSbLFFlvkzjvvbBDgrMvQoUMzduzYdOrUKUlyww03ZPbs2S1d6lr+9V//NcuWLUuSdOnSJePGjcvJJ5/cYBZbVVVVjj322Nx9993p3bt3kvr/Ebjwwgvb1Z5pwOarc+fOOfHEEzNu3Lgccsgh+ehHP5r99tsvH/3oRzNx4sSMGTMmn/jEJ1JdXd3WpQIAUAH+/Oc/5+67706SdOvWLXfffXeOOeaY0p65a+rRo0d+8IMf5IILLkhSv0/ZpZdeutZ+upTPJhWOzZ07NzfccEM+85nP5IQTTsgdd9yRf/zjH21dFgDA+7rllltKIdG5556bj3zkI43q97GPfSz/8i//kiRZsWJFbr311harcV2effbZTJw4sdT++te/vs5ZGav069cvl1xySan92muv5Xe/+12L1giwpi5duqRnz57ZZptt0rNnz2y77bZtXRIAABXmF7/4Ren4kksuyV577bXBPoVCId/5zney7777Jkn+8pe/NFiSkfJq9+FYXV1dHnjggZx77rkZMmRI/s//+T956aWXGiSmxWIxPXv2zGmnnZbf/OY3bVgtAEBDdXV1ufnmm5PU76F62mmnNan/6aefnkKhkKR+9tjKlSvLXeJ6rXrLLamflXHWWWdtsM8JJ5xQWl5xzTEAAABgczd37tyMHz8+Sf1yiieccEKj+1ZVVTXYa2z1kI3yqmnrAtbnhRdeyPjx4/Pb3/427733XpL6EKxQKKRQKKRYLKa6ujqHHHJIjj/++Bx22GGpqWm3HwcAqFCvvfZaZs2alSQ5/PDDs8MOOzSpf9++fTNkyJA88sgjmTlzZt58883suOOOLVHqWiZNmlQ6PuCAA9K9e/cN9qmurs5RRx1VWkbyqaeearH6AAAAoL15/PHHS9sTfP7zn0+XLl2a1H/EiBH51re+lQULFth3rAW1qzRp3rx5ufvuuzNu3LhMmzYtSUozxFa9MV0sFrPrrrtmxIgROe6440r7WgAAtEerXvJJkp133rlZY/Tt27d0PG/evFYLx2bMmFE6/vCHP9zofqsvF/Haa6+lrq7OPj8AAABUhHfffbd0PGDAgCb379KlS3beeec8//zzmTt3bmnSEOXV5uHYypUr88c//jF33nlnHn744dTW1q4ViCXJlltumWOOOSbHH398ac1NAID2bvXNdpu7JOLq/da3eW9LWH1v16a8kNSrV6/ScbFYzPz589OzZ8+y1gYAAADt0er/D9/cF0VXrZK3av9yyq/Nw7FDDz00c+bMSbL2solJ/RI+I0aMyCc/+cl06tSpLUsFAGiy1YOiF198sVljTJ8+vXTcmiHTqmUgkqRDhw6N7rfmz2xLly4tW00AAADQnm211Val49dff73J/evq6vLGG2+UxjJrrGW03qvH6zF79uwG7WKxmA984AM5//zz8+CDD+bGG2/MZz/7WcEYALBJ+uAHP5h+/folSf70pz81CLoaY+rUqaV9u/r3758+ffqUvcb1Wf0H+oULFza634IFC9Y7zsbwxhwAAADt3f77718KtO64444m/7/s73//+7zzzjtJksGDB5e9Puq1+cyxpD4Q69y5c4466qgcf/zxOeCAA9q6pLUcccQRpbS2qe6///4Ge4WsbubMmTnyyCM3prSSBx98sKx7kPTv37/ZfadMmVKa+gkAlaxQKOT000/Pt7/97STJL37xi1x11VWN7v+zn/2sdHzmmWe26htjvXv3Lq2V/uqrrza63yuvvFI67tq16zo3H159JlptbW2jxl19mUcAAABoj3baaacceeSReeCBBzJz5sz8/ve/z9FHH93o/j//+c9Lx2eccUZLlEjawcyxfffdN5dcckkee+yx/OAHP2iXwdimoFAorPMXTwBA2/viF7+Yzp07J6kPu+67775G9bv77rtz0003JakPmU488cQWq3FdBg4cWDpeNXutMSZNmrTOMVbXrVu30vF7773XqHGnTp3aqOtWDxBXLdUNAAAAreXMM88sHV988cVrraC3PmPHjs2DDz6Y5J8hGy2jzcOxMWPG5MQTT8yWW27Z1qU0WqFQSHV1daO/3u8N76aOtb4x999//2y99dabxGcGgErTs2fPnH/++UnqN+Y99dRTc9NNN613xtSKFSvyi1/8ImeccUYp3Pna176W7t27t1rNSXLQQQeVjmfMmJGnn356g31mzpyZP/7xj+scY3Uf/OAHS8dTpkzZ4LhvvfVWo+6fpMHPlUuWLGlUHwAAACiXoUOHZp999klSvxLLZz/72bz22mvrvb5YLObmm2/OeeedVzp34YUXprq6usVrrVTWvWuGYcOGNWk5pPezww47NPot6NV95StfyUMPPVRqDx8+vCz1rM95552Xr371qy16DwDYnI0aNSp/+9vfMm7cuKxYsSIjR47M97///Wy55ZbZZpttUl1dndra2vzXf/1XxowZk1mzZpX6fv7zn8+//uu/tnrNw4cPz+jRozN//vwkyX/8x3/kt7/9baqq1v9+1Xe+852sXLkySVJdXZ2TTz55ndcNGjSodPz0009n5syZ77s89BVXXJHly5c3qu7tttuudDxjxoxG9QEAAIByqa6uzm233ZZPfOITmTVrVl588cUMHjw4n/zkJzN37tz06NEjhUIhK1asyM0335zbb789L7zwQqn/GWeckVNPPbUNP8HmTzi2CZo7d26DN7K7dOnSpDVLAYDWV1VVlZ/85CfZaqutcsMNNySpnw2VJC+//HLpuieffLJBv7PPPjtXXHFFm8zK7tq1a84///xcccUVSZInnngiF1xwQX74wx+mY8eODa5duXJlLrvsskyYMKF07ktf+tJ6A6+jjjoqW2yxRZYsWZKVK1dm5MiRuf322xvsRZbUvz139dVX59Zbb2103fvuu2/peMqUKXnwwQctRQEAAECr2mGHHXL33XfnxBNPzIwZM1JbW5t77rmn9P1CoZBisZiHH364Qb8vf/nLufTSS63O1sLaZTg2Y8aMPPnkk3nhhRcyb968LFiwIMVisbTnRqW76667smLFilL7E5/4RLp27dqGFQEAjVFTU5Orr746xx57bH72s5/lf/7nf0qzrFZXVVWVY445JmeddVYOPfTQNqj0n0aOHJkHH3wwEydOTJLcfvvteeqpp3LyySfnwx/+cKqrq/PSSy/ltttuy1/+8pdSv9122y2XXnrpesft3r17zjjjjPz4xz9Okjz44IM5/PDDc+aZZ2b33XfP8uXLM23atIwZMybPP/98ttxyyxx++OH57W9/u8GahwwZku233z6zZs1KsVjMCSeckN133z077bRTg1DvnHPOyZAhQ5r7RwMAAADva5dddsnvf//7XHfddbn55pszZ86c0vfW3CP7gAMOyFe+8pV85jOfae0yK1K7CscmTZqUa665Zq39JIrF4npT0ieeeCJf+9rXkiQdO3bMfffdt0ntX9Ycq7+RnbT8kooAQPkUCoUcfvjhOfzwwzN16tScfvrpWbJkSerq6lJdXZ0vfvGLOfnkk993icHWVF1dnTFjxuRLX/pSaeb6yy+/nO9+97vr7bP33ntn7NixG3x556KLLsqTTz6ZP//5z0mSqVOnrnP5yK5du+bnP/95/vznPzcqHKupqcn111+fk08+OYsWLUqSTJ8+PdOnT29w3ac+9akNjgUAAAAbo0ePHhk9enT+7d/+LWPHjs3ll1+e5cuXp66uLh06dMhnPvOZnHXWWdlrr73autSKsv4NI1rZddddl9NOOy1PP/10isVi6WtDDjzwwPTp0yfz58/Pu+++m9/97netUG3beeGFFzJt2rRSe4cddsjHPvaxNqwIAGiuHXfcMf369ctee+2VffbZJ3vttVfOPffcdhOMrdK9e/eMGzcuP/zhD7PTTjut97ptt9023/3ud3P//fc32Pdrfbp06ZJx48bl9NNPT03Nut/ZOvjgg/Pggw/mqKOOalLNhx56aP70pz9l5MiRGTx4cHr16rXWko0AAADQWjp16pTjjjsugwYNysc+9rF8/OMfz+DBg/O9731PMNYG2sXMsZ/85Ce57rrrGpzbaqutsuOOO2bKlCkb7H/CCSfk8ssvT5Lcf//9OeGEE1qkzvZg/PjxDdrHHXectUcBgBZXXV2d0047LaeddlqmTp2a5557LrNnz06xWEyvXr0yYMCADBw4sMk/l3Tr1i3/+Z//mdGjR+fRRx/NzJkzU1dXl+233z6DBw/OzjvvXLp21KhRGTVqVKPH3mmnnfLtb3+7SfUAAAAAm782D8eee+65/OhHPyr9ImWXXXbJRRddlIMPPjiFQiF77bVX6urq3neMo48+OldccUWKxWKefvrp1NbWrvft401ZbW1t7r777gbnLKkIALS2AQMGZMCAAWUds2fPnjnuuOPKOiYAAADAurR5gvSjH/2otKfYHnvskVtuuaXJe4Zts8026du3b1599dUsXbo0M2bMSP/+/Vuo4vqlDUeOHJmpU6eWNtDr0aNHdtxxxwwePDhDhw7NHnvsUfb7Pvroow027Nt///3fd2mjcpo4cWJeeOGFvPDCC5k7d25qamrSs2fP9O3bN4MHD87RRx/darUAAAAAAAA0V5uGY/Pnz8+TTz6ZpH5z+u9///tNDsZWGTBgQF599dUkySuvvNKi4di0adMa7PuVJAsXLszMmTPz5JNP5tprr82QIUMyevTo9O3bt2z3nTBhQoN2a84amzRp0lrn5s+fn9deey2PPvpo/uu//ivHHHNMLr744vTq1Wuj72epyIaqq6vSu/c//3Xt3bsm1dVVbfLnVF1dnd69e691zt8ZsClpL89Vz1Rgc+G5ClBenqsA5eW5ypoKxWKx2FY3f+ihh3LuueemUChk0KBBuf3229e6ZtWyioVCYa1AanVXXnllbrrpphQKhXzrW9/KySefXNZajzjiiLzxxhtN6tOtW7dcffXVOeywwzb6/vPmzcvBBx+cFStWJEm22GKLPPbYY80OExujOQHjtttumx//+MfZZ599WqAiAAAAAACAjdOmM8fefvvt0vHee++9UWN169atdLx48eKNGmt9tttuuxx55JH5+Mc/nv79+6dXr17p2LFj5s2bl2nTpuWBBx7I+PHjs3z58iTJggULcsEFF+SWW27Jvvvuu1H3vueee0rBWJJ84hOfaNFgbJW+ffvmyCOPzAEHHJAPfehD2XrrrVNdXZ25c+fmr3/9a+69997cd999pX3h3nnnnZxzzjkZO3asZRYBAAAAAIB2p82XVVxl9XCrOZYsWVI6rqkp/8e6/PLLM3jw4HWOvc0222SbbbbJkCFDcuaZZ+b888/P9OnTkyTLli3LyJEjc99996Vjx47Nvv/48eMbtFtjScWf//znOeigg9Y5pbNPnz7p06dPhg4dmtNOOy0XXHBBZs2alSR577338s1vfjO//OUvm33vefPmNbvv5mj+/JW56KK5ue22DyVJTjrp5Vx55dbp3r2q1WtZuHBhRo0a1eDcVVdd1SphLUC5tJfnqmcqsLnwXAUoL89VgPLyXN089OjRo2xjtWk4tnogtmjRoo0a66233iodl/MPaJUDDzywUdf17ds3N954Yz73uc/lzTffTJK88cYb+dWvfpWTTjqpWfeeMWNGnn/++VJ7hx12yAEHHNCssZri4IMPbtR1++yzT37xi1/kxBNPzIIFC5IkkydPziOPPJJDDz20Wfduw9U+26W6upWZM6e21J4zpzZ1dStTLLb+WrR1dXWZM2fOWuf8nQGbkvbyXPVMBTYXnqsA5eW5ClBenqusqfWnnaymV69epeOXX355o8Z6+umnS8fbb7/9Ro21sXr16pVvfOMbDc7de++9zR5v3LhxDdrHHXdcu9ugr1+/fjnnnHManLvnnnvaqBoAAAAAAIB1a9NwbODAgUnqZwn9+c9/zsKFC5s1zh/+8IfS/mU1NTUZNGhQ2Wpsrk9+8pMNpkI+++yzDZZ+bKyVK1fmrrvuanCuNZZUbI7jjz++QWj3xBNPtGE1AAAAAAAAa2vTcGy77bbLbrvtlqR+z7Abb7yxyWMsXLgwV155ZZKkUChk8OD1bmGYAAAgAElEQVTB6dy5cznLbJaamprsvffepXZtbW3eeeedJo/zpz/9qUG//fffPzvttFNZaiy3Xr165YMf/GCp/c4772TFihVtWBEAAAAAAEBDbRqOJckZZ5yRpH722E9/+tM8/PDDje777rvv5uyzz85rr7221njtwerLRibJe++91+Qxxo8f36A9bNiwjaqppa35mefNm9dGlQAAAAAAAKytzcOxYcOGZcCAASkUClmxYkXOP//8XHrppXn11VfX2+e1117Lddddl0996lN59tlnk9TPGvv4xz+egw8+uJUq37A1l1Hs1KlTk/ovWLAgDz74YKm9xRZb5JhjjilLbS1lYz8zAAAAAABAS6pp6wIKhUKuv/76nHDCCZkzZ07q6upy++235/bbb0/Pnj1TV1dXuvakk07K66+/njlz5iSpn222ao+r7bffPldffXWbfIb1+fvf/96gvfXWWzep/+9+97ssXbq01D7qqKMa7GPW3qxcuTJvvvlmqd2hQ4d07969DSsCAAAAAABoqM1njiVJnz59cvPNN5f2H0vqg6+5c+emUCiUArBnnnkms2fPTrFYTFIfrBWLxfTv3z+33HJLevbs2Sb1r8usWbPy0ksvldq9evXKtttu26Qxxo0b16A9YsSIstTWUp555pnMnz+/1N5jjz3asBoAAAAAAIC1tYtwLEl22WWX/PrXv855552XrbbaqnR+VRC2pmKxmC5duuRf/uVf8stf/jI77LBDa5XaKNdff32D2g866KBSyNcYr732WiZPnlxqf+ADH8gBBxxQ1hrL7b//+78btNvTEpcAAAAAAABJO1hWcXUdO3bMV7/61Zx99tl5+OGHM2nSpEydOjXvvfdeFixYkM6dO6dnz5750Ic+lAMPPDCHH354gyCtJSxfvjwzZ85Mv379Gt3nzjvvzNixY0vtQqGQU045pUn3HT9+fIP2cccd16RwbU3jxo3LRRddVGoPHz48V1111TqvnTdvXpYtW5btttuu0eNfc801eeyxx0rtTp065Qtf+EKz6wUAAAAAAGgJ7SocW6Vz58455phjcswxx7R1KVm6dGmOPfbYHH300RkxYkQOOOCA1NSs+49t9uzZuf7663P77bc3OD9s2LDsvffejb5nsVjMXXfd1eDc8OHDm158M82aNSsnnnhihg8fns9+9rP5yEc+st5g7vXXX89//ud/5r777mtw/uyzz25SuAYAAAAAANAa2mU41t7U1dXlnnvuyT333JMtt9wye+65Z/r165etttoqHTp0yPz58/PCCy/k2WefzYoVKxr03X///XPJJZc06X5PPvlk3njjjVL7Ix/5SPr27VuWz9JYy5cvz5gxYzJmzJj06NEjAwYMyM4775zu3bunqqoq7733XqZMmZK//vWvWblyZYO+xxxzTM4///xWrRcAAAAAAKAxhGNNtHDhwkyaNCmTJk3a4LVf+MIXMmrUqHTs2LFJ95gwYUKDdmvOGluXefPm5fHHH8/jjz/+vtdVV1fn3HPPzVe+8pWNWgISAAAAAACgpbSLcGzZsmWZPHlyXnnllfzjH//IP/7xjyxatChdu3ZNjx490r179+y6664ZOHBgk4OmjdW5c+d8+ctfzsSJEzNlypQsX778fa/v0qVLhg4dmlNOOaVJSymusmjRotx///2l9hZbbNHqy0v26dMnp59+ep566qm8+OKLqa2tfd/rt9pqq3z605/Ol770pSbtzQYAAAAAANDa2iwcW758eX7961/nrrvuypQpUzYYwCRJhw4d8uEPfzjDhw/PsGHD0qFDhxavs2PHjhk5cmSSpLa2Nq+88kpef/31vPXWW1m0aFFqa2vTrVu3dO/ePbvttlv69++f6urqZt+va9eumTx5crnKLxkxYkRGjBjRqGt79uyZUaNGJan/e3r55Zczc+bMvPPOO1m0aFFWrlyZ7t27p0ePHunfv3923XVXM8UAAAAAAIBNQpuEY2PHjs21116bOXPmJEmKxWKD768etKz+veXLl2fy5MmZPHlyrrvuulxwwQU5/vjjW6foJDU1Ndltt92y2267tdo921rHjh0zYMCADBgwoK1LAQAAAAAA2GitGo7V1dXle9/7Xn71q1+VQq9CoVAKw1adWzMsW3XdKsViMW+//XZGjx6dv/71r/n2t7+dqqqqVvgEAAAAAAAAbMpaNRz76le/mocffjjFYrFBIFZTU5Pdd989H/7wh7PtttumW7du6dKlSxYtWpSFCxfmnXfeyfPPP5+XXnoptbW1Dfr+8pe/zNy5c/OjH/2oNT8KAAAAAAAAm6BWC8duuOGGPPTQQ6WZYsViMf369cuXvvSlHHfccenSpcsGx1i8eHEmTJiQ2267LTNmzCiNc//99+emm27Kqaee2gqfBAAAAAAAgE1Vq6xF+O677+aaa64phVlJcsYZZ+S3v/1tvvCFLzQqGEuSLl265Itf/GLuuuuunHnmmUlSGvOaa67J3LlzW+wzAAAAAAAAsOlrlXDs1ltvzZIlS0rLKf77v/97vvnNb6a6urpZ41VXV+ff/u3fMmrUqNKYixcvzq233lrmygEAAAAAANictEo4NmHChNJyikOHDs1pp51WlnFPPfXUfOITn0ixWEyxWMyECRPKMi4AAAAAAACbpxYPx6ZMmZJZs2aVllO88MILyzr+1772tRQKhSTJrFmzMnXq1LKODwCbs8uvXJmDD6v/uvzKlW1dDgAAAAC0uBYPxyZNmpSkfm+wAw88MLvuumtZx991111z4IEHltpPPfVUWccHAAAAAABg89Hi4diLL75YOj7ggANa5B6rj7v6/QAAAGgdK1cW8vKM2jz3/IqsqK3O4sVtXREAAMC61bT0DV555ZXS8cCBA1vkHoMGDVrn/QAAAGg5xWIxk59Nfjk2efOdXjluxD/+9ztb5/87uZghh6zMiGGFDBqY0nL4AAAAba3Fw7E5c+aUjnfccccWucfq486ePbtF7gEAm6NvXVSVb13U1lUAsCl6cXoxl11RzCuvJsnawdfKlYX84ZHkD48Us8vOyeiLk/67C8gAAIC21+LLKs6bN6903L179xa5R7du3ZLUv7W4+v1gc2KZGgAA2otJTxdz/gWrgrENe+XV5PwLipn0dLElywIAAGiUFp85tmzZsiT1S2hsueWWLXKPLbfcMoVCIcVisXQ/2BxYpgag5ax66WDx4mLppYMePdq6KoD278XpxVw8upglS5vWb8nS5OLRxVx3jRlkAABA22rxcGzFihVJkqqqlp2kVlVVlbq6utTV1bXofaC1WKYGaA2XX7kyv/uf+uNjPlm/zOKazv/ayjz7l/rj009Nzjy9/prnni/m7nuL+etfk9mzk0JVsu22yf77JZ8bUciOO7a/Z5KXDgA2TrFY/zNqU4OxVZYsTS6/spibfuE5CwAAtJ0WD8eAppv0dNPexl21TM0VlyWD9/dLBqBlLV1azI+uLea396z9vVdfrf+acFcx3xiZHPvp9vNM8tIBwMab/GwavZTi+vztlfpx9htUlpIAAACarMX3HAOaZmOXqXlxun0cgJZTLCbf+e4/g7Hu3ZMP75UM3DfptfU/r6utTa76QTETn2ofzyR74wCUx7gJ5Xkujv+N5ysAANB2zByDdsQyNUB7N+E3ybx/JNv3Sb721UI+fmBSVVX/vCkWi3noD8mV3y9m6f8+x/7vNcXccUvbPpPsjQNQHosWFfPHP5ZnrEcfrR+va1fPVwAAoPW1SjhWKBRSV1eXU045pcXuYa8xNgeWqQHau1XB2E+vL2TrrRv+QrNQKOTIw5PaFcmlV9TPCJg5M3nu+WTffdqiWi8dAJTTO7OTupXlGatuZTJ7TtK1a3nGA9gcrVxZyMszarN4cTEraquzeHHSo0dbVwUAm4dWmzlWLBYzadKkFhu/UCikWLQ0B5u2ci5Ts98gv8QFWsY3v7F2MLa6o4YmP/n/1f/SM0n+8lzbhWNeOgAonyVLyjve4sXlHQ9gc1AsFjP52eSXY5M33+mV40b843+/s3X+v5OLGXLIyowYVsiggV7eAoCN0Wp7jvkPNry/llimBqDcdtwxGbz/+/83vaqqkH1WC8NeebXtnkf2xgEony22KO94XbqUdzyATd2L04s55fRiLhhZzONPFJI0/Ll75cpC/vBIcsHI+uvsOQ4Azdcq4VixWGyVL9iUtcQyNQDlts+HG3fdttv883jhwpapZUO8dABQXttuk1SX6f8gq6uTbXqXZyyAzcGkp4s5/4Jio1c9eOXV5PwLipn0tJ9RAaA5WnxZxQcffLClbwGbBcvUAJuCrbdu3HWdO//zeGkz9/vaWPbGASivrl0LOeSQYv7wyMaPNeSQ+vEAqJ8xdvHopu+Tu2RpcvHoYq67Jum/u2cqwIbYy5HVtXg4tsMOO7T0LWCzYJkaYFPQoUPT+7TV5G4vHQCU34hhhfzhkY1/sA8/zi9xAZL61ZYuu6LpwdgqS5Yml19ZzE2/sKUJwLrYy5H1afFwDGicVcvUlGOWg2VqALx0ANASBg1Mdtk5jV72a1367VI/DgDJ5Gc37pmaJH97pX6c/QaVpSSAzcaL0+tfQKh/zq4dfK3ay/EPjxSzy87J6IvNxK0krbLnGLBh9cvUlGcsy9QA2BsHoCUUCoWMvriQLTpv+Np12aJz8q2LCt7KBfhf4yaUZ5mF8b+x9xjA6uzlyIYIx6AdGTGsPL8ksEwNgJcOAFpK/90LueKypgdkW3ROrris4G1cgP+1aFExf/xjecZ69NH68QDY+L0cX5zueVoJhGPQjqxapmZjWKYG4J+8dADQMgbvX8h11xQa/bNrv12S664pZPD+nqcAq7wzuzxbKyT148yeU56xADZl5drLsdhWG6jTaoRj0I5YpgagvLx0ANBy+u9eyM03FHLN/y3koI8XkzT8BUJ1dTGHH5Zc838LuekXZowBrGnJkvKOt3hxeccD2BSVcy9HNm/CMWhnLFMDUD5eOgBoWYVCIfsNKuSibyYf2Pbd3DV+q9xxa/ds13tu7rglufS7VdlvkOcowLpssUV5x+vSpbzjAWyK7OVIYwnHoB2yTA1A+XjpAKB1VFUVs2u/muyzd4d0qKnzS1qADdh2m6S6TL+Zq65OtuldnrEANlX2cqQphGPQTlmmBqB8vHQAAEB707VrIYccUp6xhhxSPx5AJbOXI01R09YFAOtXv0xN/S9pzznn3fzkp7tm0aJiLrzw1Vx7Tc984APybYDGqn/poH7d8DG/Wpk/PZ4k//wFQnV1MUMOKWT4cYUMGhhLgAEA0OJGDCvkD49s/MyE4cf52RXAXo40hXAMNhGrlqlJYpkagGby0gEAAO3JoIHJLjsnr7za/DH67VI/DkCls5cjTSEcA4AK9q2LqvKti97/mut+1PTA6MzTq3Lm6c0sqpV46QAAgLZWKBQy+uLk/AuKWbK06f236Jx866KCVQ8A8s+9HMuxtKK9HDd/Xo8GAAAAgDbSf/dCrriskC06N63fFp2TKy6zBznAKvZypCmEYwAAAADQhgbvX8h11xSyy86Nu77fLsl11xQyeH+/uAVY3Yhh5Xku2stx82dZRQAAAABoY/13L+TmG5LJzyZjfrUyf3o8Sf75y9nq6mKGHFLI8OMKGTQwllIEWAd7OdJYZo4BAAAAQDtQKBSy36BCLvpm8oFt381d47fKHbd2z3a95+aOW5JLv1uV/QbZYwxgfer3cmz6UrWr2MuxcgjHAAAAAKCdqaoqZtd+Ndln7w7pUFOXLl3auiKATYO9HGkM4RgAAAAAALDZsJcjG2LPMQAAAAAAYLNiL0fej3AMAAAAAADY7NTv5Vg/M+ycc97NT366axYtKubCC1/Ntdf0zAc+YHG9SiUcAwAAAAAANmur9nJMYi9H7DkGAAAAAABA5RCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxahp6wLYtEyfPj0vvfRS3n777VRVVaVPnz4ZMGBAdtppp7YuDQAAAAAAYIOEY410xBFH5I033mhW3/vvvz99+/Zd7/cnTpyYU045pVlj77fffrnjjjua1bcpfvOb3+RnP/tZpk+fvs7vDxo0KOeee26GDBnS4rUAAAAAAAA0l2UVeV/Lli3LhRdemG9+85vrDcaSZPLkyTnnnHNy5ZVXplgstmKFAAAAAAAAjWfmWDMUCoVUVTU+VywUCk0av6qqqtF9ampa9q/w4osvzu9+97tSu1Ao5KCDDkr//v2zYsWKPP/885k8eXKSpFgs5sYbb8wWW2yRCy+8sEXrAgAAAAAAaA7hWDMMGzYsV111VYuNf+ONN+ZjH/tYi43fWLfeemvuvvvuUrtPnz75yU9+kj333LPBdX/605/yta99LQsWLEiS/Pd//3f2228/SywCAAAAAADtjmUVWafFixfn+uuvL7U7deqUG2+8ca1gLEkOOuigXH/99Q1mu/3whz+0vCIAAAAAANDuCMdYp7Fjx+bdd98ttc8666zssssu673+ox/9aD772c+W2tOmTcvDDz/cojUCAAAAAAA0lXCMdbrvvvtKx9XV1fn85z+/wT5f/OIX1zsGAAAAAABAeyAcYy3vvvtu/vKXv5TagwYNynbbbbfBfgMHDkyfPn1K7UceeSS1tbUtUiMAAAAAAEBzCMdYyzPPPJOVK1eW2oMGDWp039WvnTdvXl5++eWy1gYAAAAAALAxhGOsZcaMGQ3aAwYMaHTfNa9dcywAAAAAAIC2VNPWBWyKXnjhhYwcOTJTp07NnDlzkiQ9evTIjjvumMGDB2fo0KHZY489mj3+mDFj8uMf/zivvvpq/vGPf6RLly7p2bNn9txzzwwePDif+tSn0r1793J9nLX87W9/a9DefvvtG913zWvXHAsAAAAAAKAtCceaYdq0aZk2bVqDcwsXLszMmTPz5JNP5tprr82QIUMyevTo9O3bt8nj33PPPQ3aS5cuzdy5czNjxozcfffd+cEPfpBTTjkl5557bjp06LBRn2Vd3n777QbtpoRjq+85liRvvfVWWWoCAAAAAAAoB+FYC3n00Udz/PHH5+qrr85hhx1W1rEXLlyY66+/Pk888USuv/76bL311mUdf/HixQ3aXbp0aXTfrl27vu9YjVUoFJrVb3NVXV2V3r3/+a9r7941qa6uapM/p+rq6vTu3Xutc/7OgE1Je3mueqYCmwvPVYDy8lwFKC/PVdZUKBaLxbYuYlNwxBFHpLa2NkceeWQ+/vGPp3///unVq1c6duyYefPmZdq0aXnggQcyfvz4LF++vNSvU6dOueWWW7Lvvvuud+yJEyfmlFNOyaBBg3LEEUdk4MCB6devX7p3754VK1Zk9uzZeeaZZ/KrX/0qzzzzTIO+gwYNyk033ZROnTqV7bN++tOfzssvv1xqP/fcc40ef8aMGfnUpz5Vah922GH56U9/WrbaAAAAAAAANoZwrJGeeOKJDB48ODU17z/Z7rXXXsv555+f6dOnl87tsMMOue+++9KxY8d19nn77bczf/787LbbbhusY/z48fmP//iPLFu2rHTuy1/+ckaOHNnIT7JhQ4cOzd///vdSe9q0aamqqmpU37///e8ZOnRoqX3AAQfkpptuKlttAAAAAAAAG8Oyio104IEHNuq6vn375sYbb8znPve5vPnmm0mSN954I7/61a9y0kknrbPPdtttl+22265R4w8fPjwdO3bM17/+9dK5m2++OaeeemrZlldcc5bYihUrGj1zbPVZc0nSuXPnZtUwb968ZvXbXM2fvzIXXTQ3t932oSTJSSe9nCuv3DrduzcutCynhQsXZtSoUQ3OXXXVVdlyyy1bvRaA5movz1XPVGBz4bkKUF6eqwDl5bm6eejRo0fZxhKOtYBevXrlG9/4RoMA6957711vONZUn/70pzNu3Lg89thjSer39XrooYfyuc99rizjr7nH2LJlyxodjq0+o21dYzWWCY0N1dWtzJw5taX2nDm1qatbmWKx9deiraury5w5c9Y65+8M2JS0l+eqZyqwufBcBSgvz1WA8vJcZU2tP+2kQnzyk59skPY+++yzWbJkSdnGXzMIe+KJJ8o29pqB1qJFixrdd81rmxuOAQAAAAAAtAThWAupqanJ3nvvXWrX1tbmnXfeKdv4gwYNatBetYRjOay5xONbb73V6L6zZs1q0O7Tp09ZagIAAAAAACgH4VgL6tWrV4P2e++9V7ax19xfrJxj77rrrg3aTQne1gzS1hwLAAAAAACgLQnHWtCayyg2dt+uxli6dGmLjb1moDV16tRG950yZUqDdr9+/cpSEwAAAAAAQDkIx1rQ3//+9wbtNWd7lXPsNWepbYxBgwalquqf/2hMnjy50X2fffbZ0nGPHj3yoQ99qGx1AQAAAAAAbCzhWAuZNWtWXnrppVK7V69e2Xbbbcs2/iOPPNKgvccee5Rt7F69emXfffcttZ999tm8/fbbG+w3efLkBssqHnrooampqSlbXQAAAAAAABtLONZCrr/++hSLxVL7oIMOSqFQKMvYCxYsyC233NLg3CGHHFKWsVc5+uijS8d1dXUZM2bMBvvccccd6x0DAAAAAACgPRCObcDy5cvzt7/9rUl97rzzzowdO7bULhQKOeWUU9Z57euvv77W3mTvZ9myZbnwwgszd+7c0rk99tgjBxxwwPv2GzVqVPr371/6Gjdu3Ptef+KJJzZYqvFnP/tZXnnllfVe/9RTT+Wuu+4qtffcc88cfvjhG/o4AAAAAAAArUo4tgFLly7Nsccem69//et57LHHUltbu95rZ8+ene9973u5+OKLG5wfNmxY9t5773X2eeihh3LkkUfmuuuue9/wKUkmTZqUz3/+83nsscdK56qqqjJq1KiyzUpbpUuXLjn33HNL7WXLluX000/PtGnT1rr28ccfz7nnnttgptzXv/71stcEAAAAAACwsWwI1Qh1dXW55557cs8992TLLbfMnnvumX79+mWrrbZKhw4dMn/+/Lzwwgt59tlns2LFigZ9999//1xyySXvO/67776ba6+9Ntdee2122GGH7Lnnntl+++3TrVu31NXVZfbs2XnmmWfy6quvrtV39OjROfDAA8v5cUtOPvnkPPPMM7nnnnuS1O+jNnz48Bx00EHp379/amtr89xzz2Xy5MkN+n35y1/OkCFDWqQmAAAAAACAjSEca6KFCxdm0qRJmTRp0gav/cIXvpBRo0alY8eOjR7/jTfeyBtvvLHB63r06JHLLrssRx11VKPHbo4rr7wydXV1ue+++5IkxWIxjz32WIPZa6usWj7ywgsvbNGaAAAAAAAAmks4tgGdO3fOl7/85UycODFTpkzJ8uXL3/f6Ll26ZOjQoTnllFPWu5Ti6j72sY/l+OOPz6RJk/L6669v8Podd9wxJ554Yk488cT07Nmz0Z+juTp16pQf/ehHmTBhQn7+859n+vTp67xu4MCBOe+888wYAwAAAAAA2jXh2AZ07NgxI0eOTJLU1tbmlVdeyeuvv5633norixYtSm1tbbp165bu3btnt912S//+/VNdXd3o8ffcc89cccUVSZIFCxbkpZdeyptvvpk5c+ZkyZIlqa6uTvfu3dO7d+/svffe2W677Zr1Oa666qpcddVVzeqb1O+bNmzYsLz44ot56aWX8vbbb6e6ujrbbrtt9tprr/Tt27fZYwMAAAAAALQW4VgT1NTUZLfddstuu+3WIuN369Yt++23X/bbb78WGb8c+vfvn/79+7d1GQAAAAAAAM1S1dYFAAAAAAAAQGsRjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxhGMAAAAAAABUDOEYAAAAAAAAFUM4BgAAAAAAQMUQjgEAAAAAAFAxatq6AGDDevSozpgxfUrtMWP6pFgstmFFAAAAAACwaTJzDAAAAAAAgIohHAMAAAAAAKBiWFYRAAAAANoR2ysAQMsSjgEAAAAAAJstLx2wJssqAgAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAHhXIaAAACAASURBVAAAAABAxahp6wIAAADYtPXoUZ0xY/qU2mPG9EmxWGzDigAAANbPzDEAAAAAAAAqhnAMAAAAAACAiiEcAwAAAAAAoGLYcwwAqDj2xgEAAACoXGaOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDGEYwAAAAAAAFQM4RgAAAAAAAAVQzgGAAAAAABAxRCOAQAAAAAAUDFq2roAGioWi3n99dfz0ksv5a233srChQvTqVOn9OjRI7vvvnv22GOPVFdXt3WZAAAAAAAAmyThWCMdccQReeONN5rV9/7770/fvn3X+/0FCxbkgQceyB/+8IdMnDgx77333nqv3XLLLTN8+PCcdtpp2XHHHZtVT2PMnDkzRx55ZLP6brfddnn00UfLXBEAAAAAAMDGE461sSeffDJnnXVWVqxY0ajrFy5cmFtuuSV33nlnRo8eneOPP76FKwQAAAAAANh8CMeaoVAopKqq8du1FQqF9X5v4cKFawVjHTt2zL777pt+/fpl6623zvLly/PSSy9l4sSJWbZsWZJk8eLFufjii7N06dKcdNJJzfsgTVBVVfW+n2N1NTX+sQIAAAAAANonKUYzDBs2LFdddVVZxywUCjnwwAPz+c9/Pocffng6deq01jWzZ8/OpZdemv/5n/8pnbvsssvykY98JHvssUdZ61nT5ZdfnhEjRrToPQAAAAAAAFpa46c/0WKGDBmSCRMm5IYbbsjRRx+9zmAsSbbZZptcc801OfbYY0vnVq5cmR/+8IetVSoAAAAAAMAmTTjWxg477LD8v//3/5o08+vb3/52unbtWmo//vjjWbBgQUuUBwAAAAAAsFkRjrWx5uzP1aNHjxx88MGl9ooVKzJt2rRylgUAAAAAALBZEo5tonbaaacG7Tlz5rRRJQD/f/buPDyq8mD/+D3ZCdkghIAsKVGSgLKErQoIKrGIIrtVUCKWRaRgBakFK1IQBa1aK4tShVJ5CyIVEAVTgQrIDm8IIGRBoEBCEgIhkJA9md8f/Dwvk32ZZJLM93NdXNc8Z57znDvBnNK5c84BAAAAAAAAgPqDcqyeunnzpsXY2dnZRkkAAAAAAAAAAADqD8qxeio2NtZi7O/vb6MkAAAAAAAAAAAA9UflH3gFxcTEaPr06Tp16pRxO0MfHx+1bt1aPXv2VFhYmEJCQmrs+PHx8YqMjDTGXl5e6tixY40dT5L+/e9/a+vWrTpz5oyuXbsmV1dX+fj4qH379urZs6cGDRqk5s2b12gGAAAAAAAAAACA6qIcq4Lo6GhFR0dbbMvIyFB8fLwOHDigxYsXq1+/fnrttdcUEBBg9eMvXbpUZrPZGD/88MNycqrZv8qdO3dajLOyspSWlqb//ve/2rZtm959912NGjVKM2fOVOPGjWs0CwAAAAAAAAAAQFVRjtWQ3bt3a+TIkXr33Xf1wAMPWG3dffv2aePGjcbY2dlZzz//vNXWr6rc3FytWbNGBw4c0Mcff1ztUtBkMlkpWcNx+/fElt8fR0dHNWvWrNg2/s4A1Dd14bzKORVAQ8J5FQCsi/MqAFgX51XczmS+/RIklOqhhx5Sfn6+BgwYoN69eys4OFi+vr5ycXFRWlqaoqOjtX37dm3cuFG5ubnGfq6urlq9erW6dOlS7QyXL1/WiBEjlJKSYmybOnWqpk2bVu21i4qPj9eAAQMUHByssLAwde/eXXfddZd8fHxkNpt15coVHTt2TJs2bdLu3bst9g0ICNC6devUpEkTq+cCAAAAAAAAAACoDsqxCtq/f7969uxZ7u0Lz58/r6lTpyouLs7Y1qpVK0VERMjFxaXKx8/JydGzzz6ro0ePGttCQ0O1evVqOTs7V3nd0mRkZOj06dMKDQ0td+6uXbv0+9//XtevXze2DR48WO+9957VcwEAAAAAAAAAAFQH5VgNuHr1qkaNGqVLly4Z215//XU9/fTTVVqvoKBAv/vd77Rt2zZjW8uWLbVu3Tr5+/tXO681HDlyRM8++6zy8/Ml3bosdfPmzQoKCqrSemlpadaM1yCYTCZ5e3tLkq5fvy5b/ehmZGRo1qxZFtsWLVokDw8Pm+QBgKqqC+dVzqkAGhLOqwBgXZxXAcC6OK/Wfz4+PlZbi2eO1QBfX1/NnDlTM2bMMLZt3bq1yuXY3LlzLYoxHx8fffrpp3WmGJOkHj16aOTIkVq3bp0kyWw269tvv61yOUZnWzaz2Wyz71FBQYGuXLlSbBt/ZwDqM1udVzmnAmioOK8CgHVxXgUA6+K8CgdbB2ioBg4caNH2RkVFKSsrq9LrvPPOO1q/fr0xdnd319/+9jfdddddVslpTaNGjbIY79+/30ZJAAAAAAAAAAAASkY5VkOcnJzUqVMnY5yfn6/Lly9Xao2PPvpIK1asMMYuLi5aunSpunTpYrWc1nTPPfdYPP/s9ttKAgAAAAAAAAAA1AWUYzXI19fXYnzt2rUK7/vPf/5TH3zwgTF2dHTUX/7yF/Xu3dtq+azNwcHB4p6flfl6AQAAAAAAAAAAagPlWA0qehtFV1fXCu331Vdf6Y033jDGJpNJb731lsLCwqyaryZkZ2cbryv69QIAAAAAAAAAANQWyrEadPHiRYtx06ZNy91n+/btevXVVy0ewPf6669r2LBhVs9nbdeuXVN6eroxLnrlHAAAAAAAAAAAgK1RjtWQxMREnT592hj7+vqqefPmZe6zf/9+TZ8+Xfn5+ca2l19+WWPGjKmxnNa0e/dui3FISIiNkgAAAAAAAAAAAJSMcqyGLFu2zOLqrz59+shkMpU6//jx45oyZYpyc3ONbS+88IImTZpUozmtJS8vT5988onFtr59+9ooDQAAAAAAAAAAQMkox8qRm5urs2fPVmqfL7/8Ul988YUxNplMCg8PL3X+6dOnNXHiRGVmZhrbwsPD9dJLL1U+cCkWL16s4OBg48/ixYtLnZucnKy0tLQKr11YWKjXX3/d4ko5Pz8/PfbYY9XKDAAAAAAAAAAAYG2UY+XIzs7W4MGDNWPGDO3Zs8filodFpaSkaN68eXr11Vcttg8bNkydOnUqcZ/4+Hj95je/sSijRo0aVWyN2nTixAk99NBDevvtt3Xq1Kky58bExGjcuHHasGGDxfYZM2bI3d29JmMCAAAAAAAAAABUmpOtA9QHBQUF2rJli7Zs2SIPDw916NBBgYGB8vb2lrOzs27cuKGYmBhFRUUpLy/PYt8ePXpo/vz5pa69adMmXb582WLbxo0btXHjxkplnDJliqZOnVqpfcpy8+ZNrVy5UitXrpSfn586duyoNm3ayNPTU2azWampqTp27JhiY2OL7Tt+/HiNGDHCalkAAAAAAAAAAACshXKskjIyMnT48GEdPny43LmjR4/WrFmz5OLiUuqc259L9rOCgoJK5yppHWtJSUnRrl27yp3n5uamWbNmafTo0TWWBQAAAAAAAAAAoDoox8rh5uamyZMn6+DBgzp58qRyc3PLnO/u7q6wsDCFh4eXeivFui44OFhjxozRoUOHdPbsWRUWFpY538/PT8OHD9fTTz+tFi1a1FJKAAAAAAAAAACAyqMcK4eLi4umT58uScrPz9e5c+d04cIFJSUl6ebNm8rPz5enp6e8vLzUvn17BQcHy9HRscLrT5s2TdOmTaup+FU6Tps2bTR37lxJUmZmpn766SfFx8frypUryszMlMlkkqenp5o2baq7775bbdq0qcnoAAAAAAAAAAAAVkM5VglOTk5q37692rdvb+sotcbd3V2dO3dW586dbR0FAAAAAAAAAACg2hxsHQAAAAAAAAAAAACoLZRjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG5RjAAAAAAAAAAAAsBuUYwAAAAAAAAAAALAblGMAAAAAAAAAAACwG062DgAAAAAAgDV4eXlp7dq1to4BAAAAoI7jyjEAAAAAAAAAAADYDcoxAAAAAAAAAAAA2A3KMQAAAAAAAAAAANgNyjEAAAAAAAAAAADYDcoxAAAAAAAAAAAA2A3KMQAAAAAAAAAAANgNyjEAAAAAAAAAAADYDcoxAAAAAAAAAAAA2A3KMQAAAAAAAAAAANgNyjEAAAAAAAAAAADYDcoxAAAAAAAAAAAA2A3KMQAAAAAAAAAAANgNyjEAAAAAAAAAAADYDcoxAAAAAAAAAAAA2A3KMQAAAAAAAAAAANgNyjEAAAAAAAAAAADYDcoxAAAAAAAAAAAA2A0nWwcAAAAAAAAAUPd4eXlp7dq1to4BAIDVceUYAAAAAAAAAAAA7AblGAAAAAAAAAAAAOwG5RgAAAAAAAAAAADsBuUYAAAAAAAAAAAA7AblGAAAAAAAAAAAAOwG5RgAAAAAAAAAAADsBuUYAAAAAAAAAAAA7AblGAAAAAAAAAAAAOwG5RgAAAAAAAAAAADsBuUYAAAAAAAAAAAA7AblGAAAAAAAAAAAAOwG5RgAAAAAAAAAAADshpOtAwAAAMA+5eTk6Pjx4zp69KgiIyMVGRmpM2fOyGw2S5JSU1Oterxr167pn//8p3bs2KHo6GilpaXJ1dVVfn5+8vPzU9euXdW3b1/dd999atq0abH9H3/8ce3du7fY9kaNGsnLy0s+Pj7q0KGDunbtqkcffVR33XWXVfMDAAAAAADroBwDAACATcyYMUNr166tlWN99913+u1vf6urV69abM/NzVV6errOnj2rgwcPavny5erevbu2bdtW4bWzsrKUlZWl5ORkxcbGatOmTZo3b57uv/9+LViwQPfcc4+1vxwAAAAAAFANlGMAAACwOQ8PD3Xp0kU//fSTkpOTrbr2gQMHNHbsWOXl5cnR0VEjR47UwIEDFRAQIEdHR12+fFnHjx/Xjh07dOjQoQqtuWfPHuN1QUGBrl+/rkuXLunIkSP6+uuvlZycrN27dyssLEyLFi3SuHHjrPo1AQAAAACAqqMcAwAAQJWsWbNGU6dOVZs2bXTs2LFK7z9gwAD16dNHoaGhCg4OloODgx5//HGrl2OvvfaaUYytX79eDzzwQLE5Dz/8sF5++WVdvHhRu3btKnfNjh07lrj917/+td544w0tX75cCxYsUG5urmbOnKmWLVtq4MCB1f1SAAAAAACAFVCOAQAAwCZGjBhR48dITExUZGSkJOmxxx4rsRi7XZs2bfTMM89U65iurq568cUX1bp1a02YMEGFhYV66aWXdPToUbm5uVVrbQAAAAAAUH0Otg4AAAAA1JT4+HjjdWBgYK0ee8SIERo6dKgkKTk5WWvWrKnV4wMAAAAAgJJRjgEAAKDBcnFxMV7HxcXV+vFfeOEF4/U333xT68cHAAAAAADFUY4BAACgwQoKCjJuZfjtt99q/fr1tXr8Hj16yMPDQ5J05MgR5efn1+rxAQAAAABAcZRjAAAAaLAaNWqksWPHSpLMZrOef/553XvvvXr99de1efNmJSQk1OjxHRwcdM8990iSMjIylJSUVKPHAwAAAAAA5XOydQAAAACgJv3pT3/SuXPntH37dkm3bq94+y0WW7durX79+umpp55S3759rX78pk2bGq+vX7+u1q1bW/0YAAAAAACg4rhyDAAAAA1ao0aNtG7dOn3yySfq2bNnsffj4+O1Zs0aDRkyRKNGjdKVK1esevzGjRsbr9PT0626NgAAAAAAqDyuHAMAALARLy8vrV271tYxynT7VU+luXjxYpnzlixZojFjxlgzVqWZTCaNHDlSI0eO1OXLl3XgwAEdO3ZMR48e1cGDB5WVlSVJ+s9//qMhQ4bou+++M54VVl0ZGRnGay8vL6usCQAAAAAAqo4rxwAAAGBXmjdvriFDhmjOnDnasGGDYmNjNX/+fLm5uUmSYmJi9PHHH1vteFevXjVe+/j4WG1dAAAAAABQNVw5BgAAgFLt2bOn1Pe+/fZbvfnmm2rZsqXWr19f6rxWrVrVRDSr8fDw0NSpU+Xh4aEZM2ZIkr766ivNnDmz2msXFhbq1KlTkiRPT0/5+/tXe00AAAAAAFA9lGMAAAAoVceOHUt9LyoqSpLk5ORU5rz6YsyYMXrllVeUn5+vc+fOWWXNQ4cOGbdV7NmzpxwdHa2yLgAAAAAAqDpuqwgAAABIcnFxMZ6d5uBgnX8m3357xsGDB1tlTQAAAAAAUD2UYwAAAGiwzGZzhefGx8crJSVFktS2bdtqH/vLL7/U5s2bJUktWrTQ6NGjq70mAAAAAACoPsoxAAAANFgxMTEaMWKE9u3bV+a87OxsTZ8+3SjTBg0aVOVj5uTk6MMPP9QLL7wg6dZVaH/961/l6upa5TUBAAAAAID18MwxAAAA2ERycrJ27Nhhse3y5cvG6zVr1li8d++99yowMLBSxzCbzdq5c6d27typwMBADRo0SN27d9cdd9whd3d3Xb16VZGRkfrss890/vx5SVLr1q01bdq0Mtc9deqU8bqwsFDXr1/XpUuXdOTIEW3evFnJycmSJFdXV73zzjt6+OGHK5UbAAAAAADUHMoxAAAA2MTp06c1derUUt8v+t6SJUsqXY41btxYPj4+SktL09mzZ7V06dIy54eGhmrFihXy8vIqc17fvn3LfN9kMqlfv35asGCB7r777kplBgAAAAAANYtyDAAAAA1WQECAYmNjtW/fPv3www/63//9X505c0ZXr15VXl6eGjdurFatWqlLly56/PHH9atf/UoODpW787ibm5u8vLzk4+Ojjh07KjQ0VI8++qjuvPPOGvqqAAAAAABAdZjMlXlKOVBLrl27ZusIdY7JZJKPj48kKS0tTbb60b1x44aef/55i23Lly8v9zfsAaCuqSvnVQBoKDivAoB1cV4FAOvivFr/NWnSxGprVe7XYgEAAAAAAAAAAIB6jHIMAAAAAAAAAAAAdoNyDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2g3IMAAAAAAAAAAAAdsPJ1gFQuri4OJ0+fVrJyclycHBQixYt1LFjR7Vt29ZmmS5cuKBTp04pKSlJhYWF8vf3V/v27RUUFGSzTAAAAAAAAAAAABVFOVZBDz30kBISEqq073fffaeAgIAKz//qq6/06aefKi4ursT3Q0NDNWXKFPXr169Keapi165d+uijj3T06NES3w8ODtaECRM0ZMiQWssEAAAAAAAAAABQWdxWsQ7JycnRSy+9pFdeeaXUYkySjh49qkmTJmnhwoUym801mslsNuvNN9/UpEmTSi3GJCk2Nla///3vNWPGDOXm5tZoJgAAAAAAAAAAgKriyrEqMJlMcnCoeK9oMpkqNO/VV1/Vt99+a7Ffnz59FBwcrLy8PJ04ccIoqMxms1atWqVGjRrppZdeqtwXUAnvvfeePvvsM4tt3bp1U6dOneTo6KjY2Fjt27fPKOm2bNkiR0dH/fnPf66xTAAAAAAAAAAAAFVFOVYFw4YN06JFi6y65v/8z//om2++McYtWrTQxx9/rA4dOljM27t3r373u98pPT1dkvTRRx+pW7duNXKLxe+//16ffPKJMfby8tKHH36o++67z2LeqVOn9MILLygpKUmStHnzZnXv3l1PPfWU1TMBAAAAAAAAAABUB7dVrAMyMzO1bNkyY+zq6qpVq1YVK8YkqU+fPlq2bJnF1Wjvv/++1W+vaDab9f777xtjk8mkZcuWFSvGJKljx45atWqVXF1djW1LlixRdna2VTMBAAAAAAAAAABUF+VYHfDFF1/o6tWrxnjChAlq165dqfN79eqlIUOGGOPo6Gh9//33Vs20fft2i+eeDR06VD179ix1frt27TR+/HhjnJKSovXr11s1EwAAAAAAAAAAQHVRjtUBERERxmtHR0c9+eST5e4zZsyYUtewdiZJevrpp8vd56mnnpKjo2ONZQIAAAAAAAAAAKguyjEbu3r1qo4dO2aMQ0ND5e/vX+5+Xbt2VYsWLYzxrl27lJ+fb5VM+fn52r17tzFu2bKlOnfuXO5+/v7+6tq1qzE+evSoUlNTrZIJAAAAAAAAAADAGijHbCwyMlKFhYXGODQ0tML73j43LS1NP/30k1UyxcXF6caNG9XOVFBQoMjISKtkAgAAAAAAAAAAsAbKMRs7c+aMxbhjx44V3rfo3KJrWStThw4dqpzp7NmzVskEAAAAAAAAAABgDU62DlAfxcTEaPr06Tp16pSuXLkiSfLx8VHr1q3Vs2dPhYWFKSQkpEJrFS2PWrZsWeEcRedaq4gqus4dd9xh80wAAAAAAAAAAADWQDlWBdHR0YqOjrbYlpGRofj4eB04cECLFy9Wv3799NprrykgIKDMtZKTky3GlSnHbn/mmCQlJSVVeN/KZCp6HFtkQt3h5eWltWvX2joGAAAAAAAAAABVQjlWQ3bv3q2RI0fq3Xff1QMPPFDqvMzMTIuxu7t7hY/RuHHjMteqqqLrFD1ObWQymUxV2q8hu/17wvcHAKqP8yoAWBfnVQCwLs6rAGBdnFdxO8qxSvD399eAAQPUu3dvBQcHy9fXVy4uLkpLS1N0dLS2b9+ujRs3Kjc3V5KUnp6uF198UatXr1aXLl1KXLNoeeTq6lrhPEXn1lQ55uLiUuuZfHx8qrSfvfD29rZ1BABoUDivAoB1cV4FAOvivAoA1sV5FZRjFfTmm2+qZ8+ecnIq/i3z8/OTn5+f+vXrp/Hjx2vq1KmKi4uTJOXk5Gj69OmKiIgosWTKycmxGDs7O1c4U9H1srOzK7xvWYpmqkw5VlOZAAAAAAAAAAAArIFyrILuu+++Cs0LCAjQqlWrNGrUKF26dEmSlJCQoPXr1+vpp58uNr/olVZ5eXkVvnrs5yvUfubm5lah/cpT9PhFj1MbmdLS0qq0X0NmMpmM32i4fv26zGazjRMBQP3GeRUArIvzKgBYF+dVALAuzqv1nzXvOEc5VgN8fX01c+ZMzZgxw9i2devWEsuxos8Yy8nJqXA5VvQKr8o8r6wsRdepTDlmrUycmMpmNpv5HgGAFXFeBQDr4rwKANbFeRUArIvzKhxsHaChGjhwoDw8PIxxVFSUsrKyis0rWh7dvHmzwscoOremyrG6kAkAAAAAAAAAAMAaKMdqiJOTkzp16mSM8/Pzdfny5WLz/P39LcZJSUkVPkZiYqLFuEWLFpVMWbK6mAkAAAAAAAAAAMAaKMdqkK+vr8X42rVrxebceeedFuOfn1NWEUVLq6JrVZU1MwUGBlolEwAAAAAAAAAAgDVQjtWgordRLOlZYkWLqFOnTlV4/ZMnT1qMrVVE1cVMAAAAAAAAAAAA1kA5VoMuXrxoMW7atGmxOaGhoXJw+L+/hqNHj1Z4/aioKOO1j4+P7rrrriqkLK59+/by8vIq8TjluT2/o6OjunXrZpVMAAAAAAAAAAAA1kA5VkMSExN1+vRpY+zr66vmzZsXm+fr66suXboY46ioKCUnJ5e7/tGjRy1uYdi/f385OTlVM/UtTk5O6tevnzFOTEzUsWPHyt0vOTnZYl5oaGiJhSAAAAAAAAAAAICtUI7VkGXLlslsNhvjPn36yGQylTj3kUceMV4XFBRo3bp15a6/du3aUtewhkGDBlmM16xZU+4+n3/+uQoKCmosEwAAAAAAAAAAQHVRjpUjNzdXZ8+erdQ+X375pb744gtjbDKZFB4eXur8X//61/L19TXGn376qc6dO1fq/EOHDmnz5s3GuEOHDnrwwQfLzLR48WIFBwcbfxYvXlzm/AEDBigoKMgYf/XVVzp8+HCp88+dO6cVK1YYYz8/Pz3xxBNlHgMAAAAAAAAAAKC2UY6VIzs7W4MHD9aMGTO0Z88e5efnlzo3JSVF8+bN06uvvmqxfdiwYerUqVOp+7m7u2vKlCnGOCcnR88995yio6OLzd23b5+mTJlicVXajBkzSr0qrapMJpNmzJhhjM1ms6ZMmaL9+/cXm3vq1CmNGzdOOTk5xrapU6fKzc3NqpkAAAAAAAAAAACqyzoPqWrgCgoKtGXLFm3ZskUeHh7q0KGDAgMD5e3tLWdnZ924cUMxMTGKiopSXl6exb49evTQ/Pnzyz3GM888o8jISG3ZskXSred8DR8+XH369FFwcLDy8/N1/PhxHT161GK/yZMnWzwfzJoefPBBTZw4UZ988okk6caNGxo3bpy6deumdUZAigAAIABJREFUzp07y8HBQbGxsdq3b59FWTdkyBA99dRTNZIJAAAAAAAAAACgOijHKikjI0OHDx8u8xaDPxs9erRmzZolFxeXCq29cOFCFRQUKCIiQtKtq7X27NmjPXv2FJv7860aX3rppcp9AZU0Y8YMZWdna/Xq1ca2yMhIRUZGljj/0Ucf1YIFC2o0EwAAAAAAAAAAQFVRjpXDzc1NkydP1sGDB3Xy5Enl5uaWOd/d3V1hYWEKDw8v81aKJXF1ddVf//pXbdq0SStWrFBcXFyJ87p27arf/va3NXbF2O0cHBz02muv6f7779eyZcsUFRVV4rygoCBNmDBBQ4cOrfFMAAAAAAAAAAAAVWUy334/PJQpPz9f586d04ULF5SUlKSbN28qPz9fnp6e8vLyUvv27RUcHCxHR0erHC82NlanT59WcnKyHB0d1bx5c919990KCAiwyvpVcf78eZ08eVKXL19WQUGB/P39ja/bmq5du2bV9RoCk8kkHx8fSVJaWpr40QWA6uG8CgDWxXkVAKyL8yoAWBfn1fqvSZMmVluLcgx1EuVYcZy8AcC6OK8CgHVxXgUA6+K8CgDWxXm1/rNmOeZgtZUAAAAAAAAAAACAOo5yDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2g3IMAAAAAAAAAAAAdoNyDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2g3IMAAAAAAAAAAAAdoNyDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2g3IMAAAAAAAAAAAAdoNyDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2g3IMAAAAAAAAAAAAdoNyDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2g3IMAAAAAAAAAAAAdoNyDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2g3IMAAAAAAAAAAAAdoNyDAAAAAAAAAAAAHaDcgwAAAAAAAAAAAB2w2Q2m822DgEAAAAAAAAAAADUBq4cAwAAAAAAAAAAgN2gHAMAAAAAAAAAAIDdoBwDAAAAAAAAAACA3aAcAwAAAAAAAAAAgN2gHAMAAAAAAAAAAIDdoBwDAAAAAAAAAACA3aAcAwAAAAAAAAAAgN2gHAMAAAAAAAAAAIDdoBwDAAAAAAAAAACA3XCydQAAZUtLS1NcXJzOnz+vtLQ0mc1meXt764477lDXrl3l6elp64gAUG8UFhbqwoULOn/+vJKTk3Xjxg3l5ubK3d1dPj4+CgkJUfv27eXo6GjrqAAAAAAAwMoSExN14sQJXbp0SZmZmXJ1dVWzZs3Url07hYSEyMXFxdYRUUsox4A6prCwUEeOHNG2bdt04MABxcXFlTrXZDLpvvvu07hx49S/f/9aTAkA9UdqaqpWrFihyMhIRUdHKysrq8z53t7eGjJkiMaPH6+WLVvWUkoAaHgWLFig1atXW2wbPny4Fi1aZKNEAAAAsEeFhYXavHmzPvvsM508ebLUec7OzgoNDdWkSZN0//3312JC2ILJbDabbR0CwP/51a9+pfPnz1d6v8cee0zz58+Xh4dHDaQCgPrrxIkTGjVqVKX3c3d315w5czRixIgaSAUADVtUVJRGjx6twsJCi+2UYwBQsoceekgJCQlV2ve7775TQECAlRMBQMMQHx+vl19+WVFRURXeZ+LEiZo5c2YNpkJdwJVjQB2TmppabNsvfvELde7cWc2aNZOrq6uSkpK0f/9+JSUlGXO2bNmiy5cva8WKFXJ1da3NyABQrzRr1kxBQUEKCAiQt7e3HB0dlZaWpujoaEVFRRkf5GZmZmr27NnKzc3VU089ZePUAFB/5OXlac6cOcWKMQAAAKA2nT59Ws8995xSUlKMbSaTSZ07d1ZISIh8fX2Vk5Nj3Grx4sWLNkyL2kY5BtRRrVq10hNPPKHhw4erRYsWxd4vKCjQF198oYULFyonJ0eSdPjwYX3wwQf6wx/+UNtxAaDOcnR0VM+ePTVw4ED16dNHgYGBpc5NSEjQ/PnztXPnTmPbW2+9pd69e6tt27a1kBYA6r+//e1vxq3B/fz8LD6MAACUz2QyycHBoVLzAQCWUlNTNWHCBIt/i4aFhWn27Nlq3bp1ifvExcVpw4YN8vT0rK2YsCFuqwjUMUOGDNGzzz6rYcOGydHRsdz5u3bt0uTJk43fzHV2dtaOHTvk7+9f01EBoEEqKCjQxIkTtXfvXmPbuHHjNHv2bBumAoD64ezZsxo6dKhyc3PVqFEjvf766xbnT26rCAAlu/22ipwrAaD6Xn75ZX3zzTfGePLkyZo+fboNE6GuqfivoQCoFRs2bNDIkSMrVIxJUv/+/fXYY48Z47y8PO3YsaOm4gFAg+fo6KiXX37ZYtsPP/xgozQAUH+YzWbNmTNHubm5kqQpU6aoVatWNk4FAAAAe7Nv3z6LYmzAgAEUYyiGcgyoY5ycKn+309vLMUk6ceKEteIAgF26++675e7ubowTExNtmAYA6ofPP/9cR44ckSQFBQXpueees3EiAAAA2KO//e1vxmtnZ2fuBIMSUY4BDUDR5+BcuXLFRkkAoOFo3Lix8Zq7UANA2ZKTk/Xee+9JuvXsm3nz5snZ2dnGqQAAAGBvLl68qAMHDhjj/v37q02bNjZMhLqKcgxoAG7evGkxrsrVZwCA/5OVlaW0tDRjzD+kAaBsb7zxhtLT0yVJv/71r9WtWzcbJwIAAIA92rp1q8UvuBa94xbwM8oxoAGIjY21GLdo0cJGSQCgYYiIiFBeXp4xfuCBB2wXBgDquO+++07btm2TJPn6+hZ7biMAAABQW6KioizGPXr0sFES1HVcXgI0AJs3b7YY33vvvTZKAgD1X0xMjN5++21j7OPjo2effdaGiQCg7kpPT9cbb7xhjGfNmiVvb28bJgIAAIA9+/HHH43X/v7+at68uSTp/Pnz2rBhg/bs2aPExERlZWWpSZMmateunfr06aOhQ4fK19fXVrFhA5RjQD136NAhHTp0yBh7enqqb9++NkwEAPWL2WxWRkaG4uLiFBERoc8//1y5ubmSJHd3d3344Ydq1qyZjVMCQN30zjvv6PLly5Kk3r17a8iQITZOBAD1X0xMjKZPn65Tp04ZzxT38fFR69at1bNnT4WFhSkkJMTGKQGg7klPTzf+bSrdekRCfn6+Pv74Y3300UfKz8+3mJ+ZmamEhATt2bNHixcv1oQJEzRlyhSZTKbajg4boBwD6rGsrCzNmTPHYttzzz2nxo0b2ygRANR9Z8+e1eDBg41xYWGhxf3If9a/f3/Nnj1b7dq1q814AFBvHD58WOvXr5ckubi4aO7cuTZOBAANQ3R0tKKjoy22ZWRkKD4+XgcOHNDixYvVr18/vfbaawoICLBRSgCoe25/drgkNWvWTK+++qq++uqrcvfNzMzUhx9+qJiYGP3lL3+RkxPVSUPHM8eAemzevHn673//a4wDAwM1YcIE2wUCgHrAbDaroKDA+FO0GHNwcNDYsWM1f/58ijEAKEVubq7mzJljnEMnT56sX/ziF7YNBQB2ZPfu3Ro5cqR27txp6ygAUGekp6dbjPfu3WsUY25ubpo8ebI2b96sY8eO6ciRI/rnP/+poUOHWlwp9t133+n999+v1dywDepPoJ76xz/+oY0bNxpjFxcX/fnPf5arq6sNUwFA/VdYWKjVq1fr888/1zPPPKMZM2bIxcXF1rEAoE5ZunSpzp07J0lq166dJk6caONEAFD/+fv7a8CAAerdu7eCg4Pl6+srFxcXpaWlKTo6Wtu3b9fGjRuNW4Cnp6frxRdf1OrVq9WlSxcbpwcA28vMzLQY/1yW+fj4aNWqVerQoYPxnpubm3r06KEePXqob9+++sMf/qDCwkJJ0sqVKzVkyBBuYdvAUY4B9VBERIQWLVpksW3+/Pm65557bJQIAOqPO++8U7GxscY4NzfX+MAhIiJCX3/9tfLy8pSXl6e///3viouL08cff0xBBgD/X2xsrFasWGGM582bxzkSAKrpzTffVM+ePUu8jZefn5/8/PzUr18/jR8/XlOnTlVcXJwkKScnR9OnT1dERATnYgB2r7Tz4Pz58y2KsaKGDBmiH3/8Uf/4xz8k3brjzMqVK/XOO+/USE7UDdxWEahn9u/fr9///vfGbzJI0ssvv6zhw4fbMBUA1F8uLi5q3ry5+vfvr4ULF+pf//qXWrZsaby/d+9eLVmyxIYJAaDuKCws1Jw5c5SXlydJGj58uH75y1/aOBUA1H/33XdfhZ5vExAQoFWrVumOO+4wtiUkJBjPgAQAe9a4ceNi2+68804NHDiw3H0nTZokZ2dnY7xz506Lz1/R8FCOAfXI8ePHNWXKFOMWCpI0fvx4TZo0yYapAKBhCQkJ0SeffGLxj+J//OMfunbtmg1TAUDdsHr1ah07dkzSrdvTvPLKKzZOBAD2x9fXVzNnzrTYtnXrVhulAYC6o6Ry7MEHH6zQvs2aNVPnzp2N8fXr1/XTTz9ZLRvqHsoxoJ6Ii4vTxIkTLe6d+8QTT/CBBADUgPbt2+vRRx81xtnZ2dq1a5cNEwGA7WVnZ+uDDz4wxq+88oqaNm1qw0QAYL8GDhwoDw8PYxwVFaWsrCwbJgIA2/P19bX4RVfp1v+/r6igoCCLcXJyslVyoW7imWNAPXDhwgX95je/UVpamrHtkUce0fz5822YCgAatt69e+urr74yxrc/pwwA7FFubq7FL2rNmTNHc+bMKXMfs9lsMd60aZM2b95sjIcNG6a33nrLukEBwA44OTmpU6dO2r9/vyQpPz9fly9fVkBAgI2TAYDtODs7q23btjpz5oyxzdvbu8L7F517/fp1q2VD3cOVY0Adl5ycrHHjxiklJcXY1r9/f7377rtycOBHGABqSrNmzSzGGRkZNkoCAHVTQUFBuX+KPqfBbDaX+T4AoOJ8fX0txtwGHACku+66y2J8++NpylN0rqurq1UyoW7ik3WgDktNTdW4ceOUkJBgbOvVq5c+/PDDYpcIAwCsq2gZ5uXlZaMkAAAAQHFFb6PIh7gAcOuz09tV5taISUlJFuMmTZpYJRPqJm6rCNRRGRkZmjBhgs6ePWts69y5sz766CO5ubnZMBkA2IdTp05ZjFu2bGmjJABQN3h5eVX6FrMHDx5UeHi4MR4+fLgWLVpk7WgAYJcuXrxoMeY5kAAghYWFacGCBcbtvSMjIy3+PVoas9msqKgoY+zo6KiQkJAaywnb48oxoA7Kzs7W5MmTdfLkSWNbUFCQPv30U4sH7gIAakZ2dra+/vpri229e/e2URoAAADAUmJiok6fPm2MfX191bx5cxsmAoC6oUWLFurWrZsx/s9//qPU1NRy99uzZ48uXbpkjDt37sznsA0c5RhQx+Tn5+t3v/udDh8+bGz7xS9+ob///e+VeoAkAODW/cJjYmIqtU9hYaHmzp1r8Y/iLl26KDAw0NrxAAAAgCpZtmyZcVWEJPXp00cmk8mGiQCg7pg2bZrxOicnR2+88UaZ82/evKk333zTYtvYsWNrJBvqDsoxoA4xm82aNWuWdu7caWxr1aqVVq1apWbNmtkuGADUU9nZ2Ro2bJhefPFFff/99+U+iPfYsWMKDw/Xpk2bjG0ODg764x//WNNRAQAAYIdyc3MtHqdQEV9++aW++OILY2wymSp0yzAAsBf33XefHnjgAWO8detW/eEPf1B6enqxuRcuXNBzzz2nc+fOGds6deqkQYMG1UZU2BDPHAPqkEuXLhW7jdelS5c0YMCASq3TqlUrbdu2zZrRAKDeMpvN+ve//61///vfatSokUJCQnTXXXfJ29tbjRo10s2bN5WUlKQTJ04Ue26DyWTSggUL1KVLFxulBwAAQEOWnZ2twYMH65FHHtGIESN07733ysmp5I/rUlJStGzZMq1Zs8Zi+7Bhw9SpU6faiAsA9cY777yjJ5980ii9Nm3apB07duj+++9X27ZtlZeXp7i4OB04cEB5eXnGfk2bNtWHH34oBweuK2roKMeAOuT2WyLcvq2goKBS61R2PgDYi6ysLB09elRHjx4td66/v7/mzZunBx98sBaSAQAAwF4VFBRoy5Yt2rJlizw8PNShQwcFBgbK29tbzs7OunHjhmJiYhQVFWXxAa4k9ejRQ/Pnz7dRcgCou7y9vfXpp5/qxRdf1MmTJyVJ6enp2rp1a6n7BAYGavny5brjjjtqKyZsiHIMAAA0WI0bN9bbb7+tH374QYcPH1ZycnK5+3Ts2FHDhw/XiBEjePguAAAAalVGRoYOHz5s8Rzy0owePVqzZs2Si4tLLSQDgPqndevWWrdunVauXKl169YpISGhxHnNmzfXs88+q2eeeUZubm61nBK2YjKXdKkKAABAA3T58mWdOXNG8fHxun79unJycuTu7i4PDw+1bt1ad999t7y8vGwdEwAAAHYiNzdXS5cu1cGDB3Xy5Mlyn5Hr7u6usLAwhYeHcytFAKgEs9msEydO6Ny5c0pJSZHJZFLTpk3VoUMHhYSE2DoebIByDAAAAAAAALCx/Px8nTt3ThcuXFBSUpJu3ryp/Px8eXp6ysvLS+3bt1dwcLAcHR1tHRUAgHqPcgwAAAAAAAAAAAB2w8HWAQAAAAAAAAAAAIDaQjkGAAAAAAAAAAAAu0E5BgAAAAAAAAAAALtBOQYAAAAAAAAAAAC7QTkGAAAAAAAAAAAAu0E5BgAAAAAAAAAAALtBOQYAAAAAAAAAAAC7QTkGAAAAAAAAAAAAu0E5BgAAAAAAAAAAALtBOQYAAAAAAAAAAAC7QTkGAAAAAAAAAAAAu0E5BgAAAAAAAAAAALtBOQYAAAAAAAAAAAC7QTkGAAAAAAAAAAAAu0E5BgAAAAAAAAAAALtBOQYAAAAAAAAAAAC7QTkGAAAAAAAAAAAAu+Fk6wAAAAAAUJb4+HgNGDDAGC9cuFAjRoywYSIUdfDgQX355Zc6duyYUlJSdPPmTeO98PBw/fGPf7RhOgAAAACwRDkGAAAAVEHRwuZnK1asUN++fSu8zoYNGzR79mxjvHTpUoWFhVklI1DT8vPz9dprr2njxo22jgIAAAAAFUY5BgAAAFjRBx98UKlyDKjPPv744xKLMUdHxxJfA3XF4sWLjde9evXSL3/5SxumAQAAQG2jHAMAAACs6MSJE9q+fTtXf6HBy8zM1MqVK41xmzZt9Kc//Undu3dXo0aNbJgMKN+SJUuM11OnTqUcAwAAsDMOtg4AAAAANDR//etfVVhYaOsYQI3av3+/xbPF3nzzTfXt25diDAAAAECdRzkGAAAAWIGT0//dlCEuLk5btmyxYRqg5kVHRxuvXVxc1KNHDxumAQAAAICKoxwDAAAArOCRRx6Ru7u7MV68eLHy8/NtmAioWdeuXTNeN2nShGeLAQAAAKg3KMcAAAAAK2jatKnCw8ON8fnz57Vx40YbJgJqVmZmpvH69isnAQAAAKCu4//BAAAAAFYyfvx4rVmzRjdu3JAkLV26VEOHDpWLi4uNk1lHbm6uDh06pEuXLik1NVWenp7q1KmTOnfuXOZ+V65c0ZEjR5SQkKCCggL5+fmpV69eatWqlVVyXb16VZGRkUpOTlZWVpb8/PzUuXNnBQYGVnvt1NRURUZGKiUlRdevX5eHh4eaNWum7t27y8/Pzwrpbzl37pxiYmKUkpKizMxMtWrVSo8//rjV1v/Zf//7X/3444+6evWqcnJy1KRJE91xxx3q3r273NzcKrWW2Wy2er7KSkhI0IkTJ5SUlKSCggL5+/ure/fuatmyZZXXTE5O1unTp3Xx4kWlp6dLkry9vdWyZUuFhobK09PTWvFlNpt1/PhxnT9/XikpKSosLFSXLl3Uq1evUve5ePGizpw5o4SEBGVkZMjR0VHe3t5q06aNunTpYtVnvtW1n/m8vDwdO3ZMFy5cUGpqqkwmk5o2baqgoCB17NhRJpOpWutXR2FhoX788UedPXtW165dU35+vpo2bap27dqpS5cuVruyMicnR0eOHFFiYqKuXr0qV1dX9e/fX+3atStx/pUrV3Ty5EnFx8crIyNDZrNZjRo1kp+fn9q2baugoKAG879RAAAAlUE5BgAAAFiJl5eXfvOb3+iDDz6QJCUmJurzzz+3uKKsOmbNmmVcjdaqVSv95z//qdB+GzZs0OzZs43xjh071Lp16xLnPvTQQ0pISJAkDR8+XIsWLVJGRoaWLFmiDRs26Pr168X2CQ4O1rx58xQaGmqxPSEhQe+88462bdumgoKCYvs9+OCDmjt3bpWLjISEBC1cuFA7d+5UXl5esffvvvtuvfrqq1V6FtaOHTu0fPlyHT9+vNQSqGvXrpo2bZr69u1b7nqLFy/WkiVLjHFsbKwkafv27Vq2bJlOnjxpMd/T09Nq5VhBQYHWr1+vFStW6MKFCyXOcXNzU1hYmF566SW1adOm1LWCg4NL3J6QkFDiez//N1RVBw8etPj5+eyzz/TLX/5SMTExevvtt7V///5ifz8mk0m9evXSH//4x1Lz3q6wsFBHjhzR1q1btXfv3lK/R5Lk4OCg3r17a/LkyerZs2eFvoaxY8fq0KFDkqRevXpp9erVKigo0MqVK7VmzRpdunTJYv6AAQMsyrG8vDzt27dPERER2rt3r5KTk0s9lpOTk8LCwvTCCy8oJCSkQvnqw8/8xYsXtWTJEm3btk03b94scY6vr6+effZZhYeHl1oQ3n4Ovd2SJUssfj5vV95/w1evXtVHH32kr7/+WmlpaSXO8fT01JNPPqmJEyfKx8en1LV+dvt/t1OnTtW0adN0/fp1vffee9qyZYsyMjKK7VO0HNu/f7+WLVumw4cPl1lkOzs7KzQ0VGPGjNGgQYPKzQYAANBQcFtFAAAAwIrCw8PVtGlTY7x8+XJlZWXZMFH1JCcn64knntDf//73Ej8kl24VPeHh4frhhx+MbUeOHNGwYcMUERFR4ofkkvT9999r9OjRio+Pr3SuAwcOaMSIEdq2bVuJxZgknTx5Us8884yWL19e4XVv3Lih8ePHa8qUKTp27FiZHypHRUVp/Pjxmjt3bpWeL7dgwQL99re/LVaMWdOVK1f0xBNPaO7cuWWWPtnZ2frmm280aNAgbdiwocbyWMOWLVv05JNPat++fSX+/ZjNZh08eFAjR47Upk2byl0vNjZWY8eO1dq1a8v8Hkm3irQ9e/Zo7NixRgleWTdu3NDYsWP17rvvFivGSrJr1y5NmjRJGzZsKLMYk6T8/HxFRERo1KhRWrt2bZXy1bWf+VWrVmnQoEHatGlTqcWYdKukev/99zV8+PBy/x6tZevWrQoLC9Pq1atLLcYkKT09XZ9++qkGDx6sEydOVPo40dHRGjp0qNatW1diMVbUu+++q3HjxunQoUPlXuGZl5enQ4cO6V//+lelcwEAANRnXDkGAAAAWFHjxo31/PPPa+HChZJulROrV6/WpEmTbJys8nJzczVp0iSdPXtW0q0rsbp16yZPT09duXJFu3btMj6sz83N1cyZMxUREaHU1FRNmjRJN2/elLOzs3r06KGQkBC5ubnpwoUL+v77743nVSUmJmr27NlavXp1hXMlJibq7bffNj6MvueeexQaGioPDw8lJSVp9+7dunr1qqRbRcn777+vRo0alXsFX2pqqp555hmdOXPG2Obs7Kxu3bopKChIXl5eyszMVExMjA4fPmwUYp9//rmys7P19ttvV/hr+OSTT4yvuXHjxurdu7fatm0rR0dHJSQkKCoqqsJrlebatWt66qmndPHiRWObyWRSaGio7rnnHrm7uysxMVG7d+/WtWvXJN36oHz27NnKzs7WmDFjiq15+63hCgsLLT54L+m2cQ4O1v19zOjoaL333nvKzc2Vg4ODevTooQ4dOqhRo0a6ePGidu/ebdwK8eevpXHjxnr44YcrtL6zs7OCgoIUGBgoPz8/ubu7K/f/tXf3cVGV+fvArwEHlEEDEUQQI1ATWcEEIfWlppLtpiWmPaCCbruVmSgafktLpRI0caMVTV3xidUVTUzN1dxyE5eXD6gkqIhpCIIPIPIoKIMz8/uDF+c3h5lhziAPY1zvv849c58znzNnDu2ey/u+lUoUFRXh/PnzyM3NBVD3u1q3bh3s7OwwY8YMyfVrNBosWLAA586dAwB0794dQ4YMgbOzMx4+fIicnJxG12+ztrZGv3794OHhga5du8LGxgYPHz7ErVu3kJ6ejtu3bwvnHhUVBQcHB4wdO1ZyfeZ2z3/55ZfYvHmz6DUPDw8MHDgQTk5OAICCggKcPHlSuOevX7+OKVOmYO/evUKfehYWFsLvVDvAk8lkBn+rhl7fvn07li1bJroHXF1d4efnB2dnZ1haWqKwsBAnT54Ursvdu3cRFhaGXbt2oW/fvkbPH6i7j2fNmiUco2/fvvDz84O9vT3Kyspw4cIF0XSSycnJ2Lhxo+gYnp6e8PX1hZOTE+RyOe7fv487d+7gypUrwrUmIiIiam8YjhERERERNbOQkBBs3rxZeIickJCAkJCQZl2rqDX85z//QW1tLZydnREbG6uzDlJNTQ0+/fRTHDhwAABQVlaGzZs3IzU1FVVVVRg8eDBiYmLQq1cv0X6FhYV45513hKkF09LScOLECQwdOlRSXevXr4dSqTRYV21tLdatW4e1a9cKr8XGxmLYsGHw9PTUe8z60KI+GJPJZHjrrbcQHh4OBwcHnf75+fn45JNPcPr0aQDAvn37MGTIEAQHB0s6h7i4OADA1KlTMX/+fNja2oreVyqVko7TmMWLF4uCsd69e2PlypXw9vbW+az4+Hj84x//EF5bvnw5Bg0apDM1X1ZWlrDd1Gk+H0d9MNanTx+sWrVKp76qqiqsWLECu3fvBlAX4C1evBh+fn6iEZ3aLC0tERQUhIkTJ2LIkCFQKBQGP//s2bNYsmSJ8DtZtWoV/vjHP8LZ2VlS/enHs4hCAAAXkUlEQVTp6VCpVLCyssKiRYvw5ptv6oQvDa+9tbU1JkyYgAkTJsDf3x/W1tZ6j63RaPDzzz8jKipK+NuzePFiDB8+XPI6ZOZ0z3///feiYKxfv35YsmQJ/Pz8dPoqlUps3rwZq1evhkqlwt27dxEZGYnExERRv5iYGMTExAAQT134wQcfIDw8XMpXBAA4c+YMYmJihGDM1dUVS5YswciRI3XWPVOpVEhOTkZ0dDQePnyI6upqzJ07F/v375e01ldSUhJUKhVcXFywfPlyPP/883rPv966deuE7a5duyIuLk7vPvVu376NI0eOIC8vz2gtRERERL8nnFaRiIiIiKiZWVtbY9asWUK7vLxcZ/TDk6C2thZdunTB9u3bdR6SA3XnGR0dDVdXV+G1hIQEZGVlwdfXF5s3b9Z5SA7UjZaJi4sThQL1D9ulUCqV6Ny5M7Zs2aK3Lrlcjjlz5iAiIkK0T2xsrMFjfvfdd0hNTRXan3zyiTDyRh83Nzds2rRJtO5UfHy8wenkGlKpVHj77bexZMkSnWAMgKSH5o05deoUfvzxR6Hds2dPJCYm6gRj9Z/14YcfisIBpVIpjH40J0qlEq6urti2bZveNbUUCgW++OILvP7668JrpaWl+Oabbwwes2/fvli7di2CgoIaDcYAwN/fH//617+EMKy2thY7duyQXH/97+Pvf/87QkJC9I5Kanjthw8fjpUrV2LYsGEGgzGgLtAdPXo0duzYIZxHWVkZ9u/fL7k+c7nnKysrERUVJbT9/PyQlJSkNxgD6r6zmTNnYunSpcJrp0+fxsmTJxs936ZQqVRYuHChcC09PDywZ88evPDCCzrBGFAXvr7xxhuIj48X3s/JyZH8N0+lUsHOzg7bt283GHLV/2by8/NFgfjHH3/caDAGAD169MCMGTNE3x0RERFRe8BwjIiIiIioBUyaNAlubm5Ce9u2bSgpKWnDipomIiJCdB4NWVlZ4bXXXhPaarUaMpkM0dHRjQY8np6e8Pf3F9q//PKLSXWFh4fDw8Oj0T7vvvuuKEBJSUnRu8aTRqNBQkKC0B4xYgRCQ0ON1iCXy/HZZ58JD7wLCgpEazA1xtXVFfPmzZPUtykaTlm3dOlSg0Ffvffffx/9+/cX2qdOnRJG+piTRYsWGT2XhQsXwtHRUWjv27ev2db+s7OzE03ReezYMZP2Hz9+PEaPHt0stejj5uYmGsFoan3mcM/v3LlTWFvL2toasbGxkka/vfnmmxg4cKDoOM3thx9+EAVQK1asMDgqUduIESPwpz/9qUm1zZ8/XxRIGlJcXCxq/+EPf5D8GURERETtDcMxIiIiIqIWIJfLMXv2bKFdVVWlsw6MubO2thY9BDfEx8dH1Pb390efPn2M7ufr6yts5+Xloba2VlJdnTp1wuTJk432s7S0xLRp04S2Wq3G0aNHdfpduHBBtM6YlGCsnqenp2jtoPppFo2ZNGnSY48OM0SpVIpCOk9PT4wYMcLofpaWljrrZ7XGVImmcHV1xZgxY4z2UygUmDRpktCurKxs1lFE2tf82rVrqKqqkrxvSEhIs9VhiHZ9GRkZkvczl3t+3759wvaYMWMkBUP1tAOotLQ00ZpgzUG7Nl9fX9E5GfPyyy8L21lZWUIA2BgbGxu8+uqrko5vY2Mjal++fFlybURERETtDdccIyIiIiJqIa+++io2btyIa9euAQB27NiBGTNmoHv37m1cmTT9+/eXNFqj4XpLhqY+a2w/jUaD+/fvw97e3uh+gYGBRqe/q9dwhE5GRoZO+JWWliZsy2QyyfXXc3d3F0ZYSX0YHRgYaNJnmOLy5cuoqakR2lLCpHpBQUGwsLCAWq0GYPqIvpZmaOo6fcaMGYP169cL7YyMDKMjtgoLC3H48GGcP38eV69eRWlpKaqqqnRCHO3ARa1Wo6ioCM8884zRmqysrEwKUxrKy8vD4cOHceHCBVy7dg3l5eWoqqrSmc5Tu76SkhIolUpJYaw53PMlJSWisFp7tJkUTz/9tLBdWlqKwsJCyWvCGaNSqZCeni60Tf1boV2bWq1Gdna20fPz8fGRvGach4cHFAqFENYuW7YM9vb2GDZsmEl1EhEREbUHDMeIiIiIiFqIhYUF5syZgzlz5gAAampqsG7dOtFaOuZMaojX8MGtk5NTk/arrq6WFI7pW2vKEAcHBzg6OuLu3bsAIHroXk870NJoNKJ1xKSoD5KAujWepDA2JeTjuH79uqjt5eUleV+FQoGnn35aOEZOTk6z1va4TLn2ffv2FQV9+q59vbKyMqxatQrJycmi6ylVRUWFpH69evWCXC43+fi3b99GdHS0aB05U1RUVKBbt25G+5nDPZ+VlSVqL1u2DNHR0ZKOD0BnpFhZWVmzhWM3btwQjfbaunUrtm3b1uTjSfl7YcrfCrlcjmnTpmHDhg0A6sLBt99+G+7u7hg1ahQGDx6M5557TtI0kERERES/dwzHiIiIiIha0NixY+Ht7Y1Lly4BAPbs2YO//OUvja7pYy6kjlZoOJKnqftJnf7M2HpTDXXr1k0Ix8rLy3XeLy0tFbUbjsIxRWVlpaR+nTt3bvJnGNPwHLXX3pLC0dFRCMekhj6txZRr37FjR9ja2grnoO/aA3XrNIWFhTUanhmjPVKvMU257jk5OQgLCxN+w00htT5zuOcb3o9NCSu1Sb0npWju2qRMq9ilSxeTjhkeHo6cnBxRkJqbm4stW7Zgy5YtAOoCt8DAQLz00ksIDAyEhQVX3CAiIqL2h+EYEREREVELkslkmDt3Lt59910AQG1tLdauXYsVK1a0cWVPLqkP4vX1r66u1nm/OQMgqQFfU0YPSdVw/StTvy/tdYtMWUurNTRcU8mYTp06CddX37UHgEWLFomCsa5duyI4OBgBAQHo1asXnJycYG1tLZqW8PTp0wgLCxPaLXXdVSoVIiIiRMGYi4sLgoOD4efnBzc3Nzg4OMDKykpU3969e7Fw4UKT6zMHzR3IPm6Apa05gzZAWm0dOpj22EYulyM+Ph6HDh3Cpk2bhH+YoS0nJwc5OTnYuXMn3N3dERkZiRdffNGkzyEiIiJ60jEcIyIiIiJqYSNHjsSgQYOEtWoOHDiAd955B56enm1c2ZPpwYMHTe6vL1zRDo+6d++O48ePN704M9BwPTZTvy/tEEnq2m6txVDAZYixa5+ZmYmUlBSh7e/vj/Xr1xsd4WVqHU115MgRYT07AHjppZewatUqo+uHtVZ9LaFhmLtlyxYMHTq0jaoR69ixo6i9bNkyvP76621UjWEymQzjxo3DuHHjkJubi1OnTuHcuXP45ZdfkJ+fL+qbm5uL2bNnY/bs2QgPD2+jiomIiIhaH8fOExERERG1gnnz5gnbKpUKq1evNvkY2lOSmTIS5El+UK7PvXv3TOpfXFwsbD/11FM672uveVRUVCR5Cjpz1fAcTZ2OT7u/qVO6tbSSkhLJfWtqakTT1um79v/973+FbZlMhhUrVkia+lD7N9WSfv75Z2FboVAgJibGaDAGtF59LaHhGmQ3btxoo0p0mXNthri7u+Ott95CbGwsfvrpJ6SkpCAqKgoDBgwQ9VuzZg3Onj3bRlUSERERtT6GY0RERERErSAgIADDhg0T2keOHMHly5dNOob2iApTRgMVFRWZ9DnmLjs7W3LfkpISUdijb7Se9msajQZnzpx5vALb2DPPPCNqm/I7q6qqQl5entD28PBotrqagynncuXKFdG0dfquff3aakDduUpdC/DixYuS63gc2vX5+fnB1tZW0n6tVV9L6N27t6htTvdjr169ROGkOdUmlbOzM0JCQrBnzx7MmjVL9N63337bRlURERERtT6GY0RERERErSQiIkLY1mg0+Prrr03aX3sUT1lZmeT1oJ7EB7iNSUtLkzwa7ujRo6K2r6+vTp8hQ4aI2gcPHmx6cWagf//+sLa2FtoNv4PGHD16VBQoPffcc81a2+M6duyY5FGT2qPCAP3XXnsNKamj5Gpra036Th9HU+orKSnB6dOnW6qkFufm5oaePXsK7ZSUlGZf6wsQr+WlUqkk7dOxY0cMHDhQaGdkZDwRo8cMmTNnDrp27Sq0TfmHB0RERERPOoZjREREREStxMfHB2PGjBHax44dw/nz5yXvrz2iQqPRSJoCKzs7W1jr7Peiuroae/bsMdpPrVZjx44dQtvCwkL0/dcbNGgQXF1dhfb333+PX3/9tXmKbQNyuRwjRowQ2r/99pukddTUajW2bt0qek3f99WWbt68qRN66VNVVYXk5GSh3blzZ50QFBCvqXbz5k1JNXz77bcmT1XZVE2pb9OmTVAqlS1VUqsYP368sF1ZWYkNGzY0+2dof7cVFRWS93vllVeEbbVajbi4uGatqzXJZDJREFlbW9uG1RARERG1LoZjREREREStaO7cubCw+P//M1z7Ab4x2iMWACAxMbHR/g8ePMCiRYtMK/AJER8fj9zc3Eb7JCQkiKbhGzlyJFxcXHT6dejQAe+9957QfvToEWbPno3CwkKTajp79qzZhBKhoaGi9meffWZ0va4NGzbg0qVLQnvo0KHo27dvi9T3OGJiYoyey5dffimaTjQ4OFg0LWk97cC5qKgIKSkpjR738uXLiI2NNbHiptOuLyMjA1euXGm0///+9z9s2bKlpctqcTNmzBCFVwkJCTh06JBJxyguLm405NYOxC9cuCD5uBMnThTte+jQIWzcuNGk2iorK036TKny8/NNWpevrKwM165dE9ra50VERET0e8dwjIiIiIioFT377LN4+eWXhfajR48k79uzZ08MHjxYaKempiI2NlbvMa5evYrQ0FBcunQJcrn88Yo2M1ZWVqioqMCMGTP0Thn56NEjrFmzBl999ZVonwULFhg85uTJk0VrwuXl5eG1117DgQMHGp1yrby8HHv27MGUKVMwdepUPHz4sIln1bwCAwPx4osvCu2CggKEhYXpXbNLqVQiLi5ONM2nlZUVPvroo1ap1RRWVlYoKCjA9OnT9QZF1dXVWLp0KXbt2iW8Zmdnp7O2Ur2goCBR+6OPPjI4IvPf//43pk+fjurqatjY2DzGWUinXZ9arcbs2bP1Bj4qlQrbt2/HrFmzoFKpWq2+lmJvb48lS5YIbY1Gg/nz5yM6OhrFxcUG91OpVDh58iQ+/fRTjB49GidOnDDYV/sfG2RmZmLNmjWS1meUy+WIiYkRTcu4atUqfPjhh8jPz29034yMDMTExGDUqFE4cOCA0c8y1ZkzZ/DCCy9g4cKFSE1NbfS/LwUFBXj//fdFU9SOHTu22WsiIiIiMlcdjHchIiIiIqLmFB4ejh9++MGkYKxeREQEpk2bJqy7lJCQgMOHD2P48OFwcHBAZWUlLl68iPPnz0OtVsPR0RFTp041eX0zczZz5kwkJibi9u3bmDZtGnx8fODr6wtbW1sUFhbi+PHjOg/QIyMj4enpafCYlpaWiIuLQ2hoqBC6FBcXY8GCBYiJiUFAQABcXV1hY2ODqqoq3Lt3D9nZ2fjtt98kr1fU2r744gtkZ2cLD+yvXr2KiRMnws/PD97e3rCxscGdO3eQkpKiM9pk4cKF6NevX1uU3aj58+fjq6++wq+//org4GD4+/vD29sb1tbWKCgowPHjx0VT5FlYWGDZsmWidZW0eXt7IygoCD/99BMAoLS0FFOnToWfnx98fX3RsWNHFBcX48SJEygoKAAAdOrUCfPmzUN0dHSLn29QUBC8vb2FEX03btzAhAkTMHToUHh5ecHS0hKFhYVITU0Vpnrs1q0bQkNDn+jp/oC60X55eXn45ptvANQFZImJidi5cycGDBgALy8v2NnZoba2FuXl5bh+/TqysrJw//59ScefPHkykpKShDX24uPjER8fD7lcLvoHBa+88go+//xz0b7PP/88li5diqVLlwr7Hzx4EIcOHYKXlxcGDBgAe3t7aDQaVFRU4MaNG7h48SLKysqa46tpVE1NDfbu3Yu9e/dCoVDAy8sLnp6esLOzg1wuR3l5Oa5cuYL09HTRf4O8vb0xceLEFq+PiIiIyFwwHCMiIiIiamXu7u4IDg6WtG5WQ/7+/oiMjBRN7Xbz5k0kJSXp9O3evTvWrVtndCq2J02PHj3w9ddfY+7cuSgvL0dmZiYyMzP19pXJZIiIiMD06dONHvepp55CUlISFi9ejIMHDwqvl5aW4siRI0b379Chg2jKzLZmb2+PnTt34r333hPClfq16gyNjpLL5YiKisLkyZNbs1TJ+vfvj+XLl2PRokWoqalBWloa0tLS9Pbt0KEDoqOjRSPo9Fm+fDny8/NF98m5c+dw7tw5nb4KhQLx8fGiUUMtycLCAqtXr8a0adNw+/ZtAHUjyFJTU5GamqrT39HRERs3btQ7QvBJNHfuXLi7uyMqKkoY4VRbW4v09HRJaylaWVkZfM/b2xv/93//h5UrVwoBV/3xtdfeMjRV6htvvAFXV1dERkYK4bJarcalS5dE05M2pbbmUlVV1ej9Xs/X1xfr169vtd81ERERkTkwn//nRkRERETUjnzwwQdNfjj617/+FWvWrEGvXr30vm9tbY0JEyZg37598Pb2fpwyzdaQIUOQnJyMUaNGGXyg6+3tjX/+85+YOXOm5OPa2Njgb3/7G3bt2oWgoCB07Nix0f5yuRwBAQH4+OOPkZKSAltbW5POo6U5Ojpi9+7diIqKMvh7Aep+M+PGjcOhQ4fMNhirN378eCQlJSEgIAAymUznfZlMhoCAACQnJyM4ONjo8bp06YKkpCSEhYUZvN6dOnXC+PHjsX//ftH0m62hZ8+e2Lt3L4KDgw3+1rt06YKQkBDs378fXl5erVpfS5swYQKOHj2KmTNnonv37kb7u7u7IzQ0FLt378aUKVMa7fvnP/8Z3333HaZPnw4fHx9hdJVUw4YNw48//ojIyMhG7696Li4umDx5MrZu3YrIyEjJnyPV6NGjERUVhVGjRqFLly5G+z/77LP4/PPPsXPnToOjK4mIiIh+r2Sa+vlYiIiIiIjoiaLRaJCVlYWsrCyUlJRAoVDA2dkZAQEBkh6M/l4UFxcjPT0dhYWFePDgAZycnODj4wMPD4/HPrZSqURGRgby8/NRVlaGmpoa2NjYoGvXrnB3d0efPn2MBmjmJCcnB5cuXcK9e/dQU1MDe3t7uLi4wN/f3yzP4/Tp0wgLCxPaiYmJCAwMFNoFBQXIzMxEUVERVCoVnJyc4OfnBxcXlyZ93v3793H27Fnk5eXhwYMHcHBwgJOTE/z9/aFQKB77fB5XSUkJzp49i1u3bkGpVMLBwQEuLi7w8/NrlZFI5iAnJwfZ2dkoLS1FRUUFrKys0LlzZ7i5uaF3795wdHRss9pu3bqFzMxMlJSUoKKiApaWlrC1tUXPnj3h6enZ5N9lU2g0GuTk5CA3Nxd37tzB/fv3odFooFAo0KNHD/Tv379V6yEiIiIyNwzHiIiIiIiIyCwZC8eIiIiIiIiagtMqEhERERERERERERERUbvBcIyIiIiIiIiIiIiIiIjaDYZjRERERERERERERERE1G4wHCMiIiIiIiIiIiIiIqJ2g+EYERERERERERERERERtRsMx4iIiIiIiIiIiIiIiKjdkGk0Gk1bF0FERERERERERERERETUGjhyjIiIiIiIiIiIiIiIiNoNhmNERERERERERERERETUbjAcIyIiIiIiIiIiIiIionaD4RgRERERERERERERERG1GwzHiIiIiIiIiIiIiIiIqN1gOEZERERERERERERERETtBsMxIiIiIiIiIiIiIiIiajcYjhEREREREREREREREVG7wXCMiIiIiIiIiIiIiIiI2g2GY0RERERERERERERERNRuMBwjIiIiIiIiIiIiIiKidoPhGBEREREREREREREREbUbDMeIiIiIiIiIiIiIiIio3WA4RkRERERERERERERERO0GwzEiIiIiIiIiIiIiIiJqNxiOERERERERERERERERUbvBcIyIiIiIiIiIiIiIiIjaDYZjRERERERERERERERE1G78P99QJL4cJ93jAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"image/png": {
"height": 544,
"width": 867
}
},
"output_type": "display_data"
}
],
"source": [
"num_param = np.arange(2, param + 1)\n",
"\n",
"plt.figure(figsize=(10, 6))\n",
"plt.scatter(num_param, r[:, 0], color='C0')\n",
"plt.xticks(num_param)\n",
"\n",
"for j in range(param - 1):\n",
" plt.vlines(num_param[j], r[j,0] - r[j, 1], r[j,0] + r[j,1], color='mediumblue', \n",
" zorder=-1, alpha=0.80)\n",
"\n",
"plt.scatter(num_param + 0.1, r[:, 2], facecolors='none', edgecolors='k')\n",
"\n",
"for j in range(param - 1):\n",
" plt.vlines(num_param[j] + 0.1, r[j,2] - r[j, 3], r[j,2] + r[j,3], color='k', \n",
" zorder=-2, alpha=0.70) \n",
"\n",
"dist = 0.20\n",
"plt.text(num_param[1] - dist, r[1, 0] - dist, 'in', color='C0', fontsize=13)\n",
"plt.text(num_param[1] + dist, r[1, 2] - dist, 'out', color='k', fontsize=13)\n",
"plt.text(num_param[1] + dist, r[1, 2] + r[1,3] - dist, '+1 SD', color='k', fontsize=10)\n",
"plt.text(num_param[1] + dist, r[1, 2] - r[1,3] - dist, '+1 SD', color='k', fontsize=10)\n",
"plt.xlabel('Number of parameters', fontsize=14)\n",
"plt.ylabel('Deviance', fontsize=14)\n",
"plt.title(f'N = {n}', fontsize=14)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.15"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"data = pd.read_csv('Data/cars.csv', sep=',')"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [sigma, b, a]\n",
"Sampling 2 chains: 100%|██████████| 30000/30000 [00:31<00:00, 943.55draws/s] \n"
]
}
],
"source": [
"with pm.Model() as m_6_15 :\n",
" a = pm.Normal('a', mu=0, sd=100)\n",
" b = pm.Normal('b', mu=0, sd=10)\n",
" sigma = pm.Uniform('sigma', 0, 30)\n",
" mu = pm.Deterministic('mu', a + b * data['speed'])\n",
" dist = pm.Normal('dist', mu=mu, sd=sigma, observed = data['dist'])\n",
" m_6_15 = pm.sample(5000, tune=10000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.16"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"n_samples = 1000\n",
"n_cases = data.shape[0]\n",
"ll = np.zeros((n_cases, n_samples))"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"for s in range(0, n_samples):\n",
" mu = m_6_15['a'][s] + m_6_15['b'][s] * data['speed']\n",
" p_ = stats.norm.logpdf(data['dist'], loc=mu, scale=m_6_15['sigma'][s])\n",
" ll[:,s] = p_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.17"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"n_cases = data.shape[0]\n",
"lppd = np.zeros(n_cases)\n",
"for a in range(1, n_cases):\n",
" lppd[a,] = logsumexp(ll[a,]) - np.log(n_samples)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.18"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"pWAIC = np.zeros(n_cases)\n",
"for i in range(1, n_cases):\n",
" pWAIC[i,] = np.var(ll[i,])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.19"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"412.31568650728855"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"- 2 * (sum(lppd) - sum(pWAIC))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.20"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15.162740093813598"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"waic_vec = - 2 * (lppd - pWAIC)\n",
"(n_cases * np.var(waic_vec))**0.5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.21"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(17, 9)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d = pd.read_csv('Data/milk.csv', sep=';')\n",
"d['neocortex'] = d['neocortex.perc'] / 100\n",
"d.dropna(inplace=True)\n",
"d.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.22"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"a_start = d['kcal.per.g'].mean()\n",
"sigma_start = d['kcal.per.g'].std()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [sigma, alpha]\n",
"Sampling 2 chains: 100%|██████████| 4000/4000 [00:01<00:00, 2441.14draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [sigma, beta, alpha]\n",
"Sampling 2 chains: 100%|██████████| 40000/40000 [01:46<00:00, 376.87draws/s]\n",
"There were 14 divergences after tuning. Increase `target_accept` or reparameterize.\n",
"The acceptance probability does not match the target. It is 0.8907813033282677, but should be close to 0.8. Try to increase the number of tuning steps.\n",
"The number of effective samples is smaller than 25% for some parameters.\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [sigma, beta, alpha]\n",
"Sampling 2 chains: 100%|██████████| 4000/4000 [00:02<00:00, 1496.02draws/s]\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [sigma, beta, alpha]\n",
"Sampling 2 chains: 100%|██████████| 40000/40000 [03:02<00:00, 219.41draws/s]\n",
"There were 93 divergences after tuning. Increase `target_accept` or reparameterize.\n",
"There were 23 divergences after tuning. Increase `target_accept` or reparameterize.\n",
"The number of effective samples is smaller than 25% for some parameters.\n"
]
}
],
"source": [
"mass_shared = theano.shared(np.log(d['mass'].values))\n",
"neocortex_shared = theano.shared(d['neocortex'].values)\n",
"\n",
"with pm.Model() as m6_11:\n",
" alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)\n",
" mu = alpha + 0 * neocortex_shared\n",
" sigma = pm.HalfCauchy('sigma',beta=10, testval=sigma_start)\n",
" kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])\n",
" trace_m6_11 = pm.sample(1000, tune=1000) \n",
"\n",
"with pm.Model() as m6_12:\n",
" alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)\n",
" beta = pm.Normal('beta', mu=0, sd=10)\n",
" sigma = pm.HalfCauchy('sigma',beta=10, testval=sigma_start)\n",
" mu = alpha + beta * neocortex_shared\n",
" kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])\n",
" trace_m6_12 = pm.sample(5000, tune=15000)\n",
" \n",
"with pm.Model() as m6_13:\n",
" alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)\n",
" beta = pm.Normal('beta', mu=0, sd=10)\n",
" sigma = pm.HalfCauchy('sigma', beta=10, testval=sigma_start)\n",
" mu = alpha + beta * mass_shared\n",
" kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])\n",
" trace_m6_13 = pm.sample(1000, tune=1000)\n",
" \n",
"with pm.Model() as m6_14:\n",
" alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)\n",
" beta = pm.Normal('beta', mu=0, sd=10, shape=2)\n",
" sigma = pm.HalfCauchy('sigma', beta=10, testval=sigma_start)\n",
" mu = alpha + beta[0] * mass_shared + beta[1] * neocortex_shared\n",
" kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])\n",
" trace_m6_14 = pm.sample(5000, tune=15000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.23"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:993: UserWarning: For one or more samples the posterior variance of the log predictive\n",
" densities exceeds 0.4. This could be indication of WAIC starting to fail see\n",
" http://arxiv.org/abs/1507.04544 for details\n",
" \n",
" \"\"\"\n"
]
},
{
"data": {
"text/plain": [
"waic -16.9226\n",
"waic_se 4.87234\n",
"p_waic 3.00199\n",
"warning True\n",
"waic_i [-0.5820018909358816, -1.5105480636019253, -1....\n",
"waic_scale deviance\n",
"dtype: object"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"az.waic(trace_m6_14, m6_14)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.24"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:993: UserWarning: For one or more samples the posterior variance of the log predictive\n",
" densities exceeds 0.4. This could be indication of WAIC starting to fail see\n",
" http://arxiv.org/abs/1507.04544 for details\n",
" \n",
" \"\"\"\n",
"/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:993: UserWarning: For one or more samples the posterior variance of the log predictive\n",
" densities exceeds 0.4. This could be indication of WAIC starting to fail see\n",
" http://arxiv.org/abs/1507.04544 for details\n",
" \n",
" \"\"\"\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" waic | \n",
" p_waic | \n",
" d_waic | \n",
" weight | \n",
" se | \n",
" dse | \n",
" warning | \n",
" waic_scale | \n",
"
\n",
" \n",
" \n",
" \n",
" | m6_14 | \n",
" -16.9226 | \n",
" 3.00199 | \n",
" 0 | \n",
" 0.957972 | \n",
" 4.87234 | \n",
" 0 | \n",
" True | \n",
" deviance | \n",
"
\n",
" \n",
" | m6_13 | \n",
" -9.08342 | \n",
" 1.94869 | \n",
" 7.83922 | \n",
" 0.0190146 | \n",
" 4.05308 | \n",
" 3.33736 | \n",
" True | \n",
" deviance | \n",
"
\n",
" \n",
" | m6_11 | \n",
" -8.72321 | \n",
" 1.30887 | \n",
" 8.19943 | \n",
" 0.0158807 | \n",
" 3.55127 | \n",
" 4.60363 | \n",
" False | \n",
" deviance | \n",
"
\n",
" \n",
" | m6_12 | \n",
" -7.12231 | \n",
" 1.99539 | \n",
" 9.80033 | \n",
" 0.00713244 | \n",
" 3.23372 | \n",
" 4.75175 | \n",
" False | \n",
" deviance | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {
"image/png": {
"height": 411,
"width": 611
}
},
"output_type": "display_data"
}
],
"source": [
"az.plot_compare(compare_df);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.26"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.17745"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff = np.random.normal(loc=6.7, scale=7.26, size=100000)\n",
"sum(diff < 0) / 100000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Code 6.27"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare function already checks number of observations to be equal."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" m6_11 | \n",
" m6_12 | \n",
" m6_13 | \n",
" m6_14 | \n",
"
\n",
" \n",
" \n",
" \n",
" | alpha | \n",
" 0.657407 | \n",
" 0.364541 | \n",
" 0.701868 | \n",
" -1.059950 | \n",
"
\n",
" \n",
" | beta | \n",
" NaN | \n",
" 0.434416 | \n",
" -0.030158 | \n",
" NaN | \n",
"
\n",
" \n",
" | beta__0 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" -0.095699 | \n",
"
\n",
" \n",
" | beta__1 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" 2.753594 | \n",
"
\n",
" \n",
" | sigma | \n",
" 0.188607 | \n",
" 0.191718 | \n",
" 0.182601 | \n",
" 0.140068 | \n",
"
\n",
" \n",
"
\n",
"