"
],
"text/plain": [
" biomarker theta_mean theta_std phi_mean phi_std\n",
"0 GMI(FUS) 0.457107 0.044667 0.569150 0.035734\n",
"1 FCI(HIP) 2.794918 2.983066 -5.546734 2.887983\n",
"2 FCI(PCC) 11.583134 3.358987 1.781998 3.292063\n",
"3 FCI(Fusi) -19.767018 4.183856 -10.216922 3.021294\n",
"4 GMI(HIP) 0.341778 0.054571 0.482245 0.037704"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"theta_phi_kmeans"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Metropolis-Hastings Algorithm Implementation Soly Dependent on Theta_Phi_Kmeans"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def fill_up_pdata(pdata, k_j):\n",
" '''Fill up a single participant's data using k_j; basically add two columns: \n",
" k_j and affected\n",
" Note that this function assumes that pdata already has the S_n column\n",
" \n",
" Input:\n",
" - pdata: a dataframe of ten biomarker values for a specific participant \n",
" - k_j: a scalar\n",
" '''\n",
" data = pdata.copy()\n",
" data['k_j'] = k_j\n",
" data['affected'] = data.apply(lambda row: row.k_j >= row.S_n, axis = 1)\n",
" return data "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def compute_single_measurement_likelihood(theta_phi, biomarker, affected, measurement):\n",
" '''Computes the likelihood of the measurement value of a single biomarker\n",
" We know the normal distribution defined by either theta or phi\n",
" and we know the measurement. This will give us the probability\n",
" of this given measurement value. \n",
"\n",
" input:\n",
" - theta_phi: the dataframe containing theta and phi values for each biomarker\n",
" - biomarker: an integer between 0 and 9 \n",
" - affected: boolean \n",
" - measurement: the observed value for a biomarker in a specific participant\n",
"\n",
" output: a scalar\n",
" '''\n",
" biomarker_params = theta_phi[theta_phi.biomarker == biomarker].reset_index()\n",
" mu = biomarker_params['theta_mean'][0] if affected else biomarker_params['phi_mean'][0]\n",
" std = biomarker_params['theta_std'][0] if affected else biomarker_params['phi_std'][0]\n",
" var = std**2\n",
" likelihood = np.exp(-(measurement - mu)**2/(2*var))/np.sqrt(2*np.pi*var)\n",
" return likelihood"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def compute_likelihood(pdata, k_j, theta_phi):\n",
" '''This implementes the formula of https://ebm-book2.vercel.app/distributions.html#known-k-j\n",
" This function computes the likelihood of seeing this sequence of biomarker values for a specific participant\n",
" '''\n",
" data = fill_up_pdata(pdata, k_j)\n",
" likelihood = 1\n",
" for i, row in data.iterrows():\n",
" biomarker = row['biomarker']\n",
" measurement = row['measurement']\n",
" affected = row['affected']\n",
" likelihood *= compute_single_measurement_likelihood(\n",
" theta_phi, biomarker, affected, measurement)\n",
" return likelihood"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def metropolis_hastings_with_theta_phi_kmeans(\n",
" data_we_have, iterations, theta_phi, non_diseased_participants,\n",
" burn_in, thining, \n",
" ):\n",
" '''Implement the metropolis-hastings algorithm\n",
" Inputs: \n",
" - data: data_we_have\n",
" - iterations: number of iterations\n",
"\n",
" Outputs:\n",
" - best_order: a numpy array\n",
" - best_likelihood: a scalar \n",
" '''\n",
" num_participants = len(data_we_have.participant.unique())\n",
" num_biomarkers = len(data_we_have.biomarker.unique())\n",
" num_stages = num_biomarkers + 1\n",
"\n",
" all_dicts = []\n",
"\n",
" # initialize an ordering and likelihood\n",
" # note that it should be a random permutation of numbers 1-10\n",
" best_order = np.random.permutation(np.arange(1, num_stages))\n",
" biomarker_names = np.array(list(data_we_have.biomarker.unique()))\n",
" biomarker_best_order_dic = dict(zip(biomarker_names, best_order))\n",
"\n",
" best_likelihood = -np.inf\n",
"\n",
" # initialize participant_stages \n",
" # note that high should be num_stages; otherwise, no participants will be in the stage of 10\n",
" participant_stages = np.random.randint(low = 0, high = num_stages, size = num_participants)\n",
" participant_stages[non_diseased_participants] = 0\n",
"\n",
" for _ in range(iterations):\n",
" new_order = best_order.copy()\n",
" # randomly select two indices\n",
" a, b = np.random.choice(num_biomarkers, 2, replace=False)\n",
" # swapping the order\n",
" new_order[a], new_order[b] = new_order[b], new_order[a]\n",
" # biomarker - order dict\n",
" ordering_dic = dict(zip(biomarker_names, new_order))\n",
" # fill up S_n column using the ordering dict\n",
" # copy first in order not to change data_we_have\n",
" data = data_we_have.copy()\n",
" # now data_we_have has S_n column\n",
" data['S_n'] = data.apply(lambda row: ordering_dic[row['biomarker']], axis = 1)\n",
"\n",
" all_participant_ln_likelihood = 0 \n",
" for p in range(num_participants):\n",
" # copy participant_stages \n",
" participant_stages_copy = participant_stages.copy()\n",
" # this participant data\n",
" pdata = data[data.participant == p].reset_index(drop=True)\n",
"\n",
" \"\"\"If this participant is not diseased (i.e., if we know k_j is equal to 0)\n",
" We still need to compute the likelihood of this participant seeing this sequence of biomarker data\n",
" but we do not need to estimate k_j like below\n",
"\n",
" We still need to compute the likelihood because we need to add it to all_participant_ln_likelihood\n",
" \"\"\"\n",
" if p in non_diseased_participants:\n",
" # the following will update pdata's kj and affect columns\n",
" this_participant_likelihood = compute_likelihood(\n",
" pdata, k_j = 0, theta_phi = theta_phi)\n",
" this_participant_ln_likelihood = np.log(this_participant_likelihood)\n",
" else:\n",
" # initiaze stage_likelihood\n",
" stage_likelihood = np.zeros(num_stages)\n",
" for k_j in range(num_stages):\n",
" # even though data above has everything, it is filled up by random stages\n",
" # we don't like it and want to know the true k_j. All the following is to update participant_stages\n",
"\n",
" # likelihood for this participant to have this specific sequence of biomarker values\n",
" participant_likelihood = compute_likelihood(pdata, k_j, theta_phi)\n",
"\n",
" # update each stage likelihood for this participant\n",
" stage_likelihood[k_j] = participant_likelihood\n",
" likelihood_sum = np.sum(stage_likelihood)\n",
" normalized_stage_likelihood = [l/likelihood_sum for l in stage_likelihood]\n",
" sampled_stage = np.random.choice(np.arange(num_stages), p = normalized_stage_likelihood)\n",
" participant_stages_copy[p] = sampled_stage \n",
"\n",
" # if participant is at sample_stage, \n",
" # what is the likelihood of this participant having this sequence of biomarker data:\n",
" this_participant_likelihood = stage_likelihood[sampled_stage]\n",
"\n",
" # then, update all_participant_likelihood\n",
" if this_participant_likelihood == 0:\n",
" this_participant_ln_likelihood = np.log(this_participant_likelihood + 1e20)\n",
" else:\n",
" this_participant_ln_likelihood = np.log(this_participant_likelihood)\n",
" \"\"\"\n",
" All the codes in between are calculating this_participant_ln_likelihood. If we already know kj=0, then\n",
" it's very simple. If kj is unknown, we need to calculate the likelihood of seeing this sequence of biomarker\n",
" data at different stages, and get the relative likelihood before we get a sampled stage. Then we calculate\n",
" this_participant_ln_likelihood again. \n",
" \"\"\"\n",
" all_participant_ln_likelihood += this_participant_ln_likelihood\n",
" \n",
" \"\"\"\n",
" The key to both `metropolis_hastings_with_theta_phi_kmeans` and `metropolis_hastings` is to \n",
" compare best_likelihood and the likelihood of all participants having specific sequences of measurements\n",
" based on the assumed S_n. \n",
"\n",
" The difference lies in how to calculate all_participant_ln_likelihood. \n",
"\n",
" `metropolis_hastings_with_theta_phi_kmeans` tries to obtain k_j and calculate each likelihood exactly \n",
" whereas `metropolis_hastings` did not obtain exact k_j and calculate the average likelihood instead\n",
" \"\"\"\n",
" acceptance_ratio = np.exp(all_participant_ln_likelihood - best_likelihood)\n",
" random_number = np.random.rand()\n",
" if random_number < acceptance_ratio:\n",
" best_likelihood = all_participant_ln_likelihood\n",
" biomarker_best_order_dic = ordering_dic\n",
" participant_stages = participant_stages_copy\n",
" \n",
" if _ >= burn_in and _ % thining == 0:\n",
" all_dicts.append(ordering_dic)\n",
"\n",
" if (_+1) % 10 == 0:\n",
" print(f\"iteration {_ + 1} done\")\n",
" return biomarker_best_order_dic, participant_stages, all_dicts"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iteration 10 done\n",
"iteration 20 done\n",
"iteration 30 done\n",
"iteration 40 done\n",
"iteration 50 done\n",
"iteration 60 done\n",
"iteration 70 done\n",
"iteration 80 done\n",
"iteration 90 done\n",
"iteration 100 done\n",
"iteration 110 done\n",
"iteration 120 done\n",
"iteration 130 done\n",
"iteration 140 done\n",
"iteration 150 done\n",
"iteration 160 done\n",
"iteration 170 done\n",
"iteration 180 done\n",
"iteration 190 done\n",
"iteration 200 done\n"
]
}
],
"source": [
"biomarker_best_order_dic, participant_stages, all_dicts = metropolis_hastings_with_theta_phi_kmeans(\n",
" data_we_have, iterations, theta_phi = theta_phi_kmeans, \n",
" non_diseased_participants = non_diseased_participants,\n",
" burn_in = burn_in, thining = thining, )"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def get_biomarker_stage_probability(df, num_biomarkers):\n",
" # Create an empty list to hold dictionaries\n",
" dict_list = []\n",
"\n",
" # for each biomarker\n",
" for col in df.columns:\n",
" dic = {\"biomarker\": col}\n",
" # get the frequency of biomarkers\n",
" stage_counts = df[col].value_counts()\n",
" # for each stage\n",
" for i in range(1, num_biomarkers + 1):\n",
" # get stage:prabability\n",
" dic[i] = stage_counts.get(i, 0)/len(df)\n",
" dict_list.append(dic)\n",
"\n",
" dff = pd.DataFrame(dict_list)\n",
" dff.set_index(dff.columns[0], inplace=True)\n",
" return dff "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAHFCAYAAADcytJ5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACYyElEQVR4nOzdd1QT2dsH8G/oRaU3RSlKkyIoPwv2hopdLLgq2PvuYgMRFcVVVCxYQBc7KsWC6AqrsrZVsXdhsSOgSAcVlDrvHx7mNSYghMCAPJ9zcg65c+fmzmQSntw2PIZhGBBCCCGE1EESXFeAEEIIIaQ8FKgQQgghpM6iQIUQQgghdRYFKoQQQgipsyhQIYQQQkidRYEKIYQQQuosClQIIYQQUmdRoEIIIYSQOosCFUIIIYTUWQ0mUNm/fz94PB7u3LkjdPugQYOgr69fo3WIiYnBihUrkJOTU6OvU1ecP38etra2UFRUBI/HQ0REhNB8CQkJ4PF4fI8mTZqgTZs28PPzQ0lJCV/+Hj16oEePHjV/ADVk4sSJaNSoEdfVKNelS5fA4/Fw7NgxrqsCANDX18fEiRO5rkatys/Px4oVK3Dp0qVK5S/7DO3fv19sdWAYBqGhoejatSs0NTUhJycHXV1d9OvXD7t37xa5roRUlRTXFWhIYmJisHLlSkycOBHKyspcV6dGMQyD0aNHw9jYGKdOnYKioiJMTEwq3OfXX3/FL7/8AgDIycnBqVOnMG/ePCQlJWHjxo1svoCAgBqtOyFcy8/Px8qVKwGgUkG5jo4Orl+/jpYtW4qtDh4eHli3bh2mTZuGRYsWoXHjxnjz5g0uXLiAkydPYurUqSLVlZCqokCF1Ih3794hKysLw4cPR+/evSu1T4sWLdCxY0f2ef/+/fHkyROEhITwBSqtW7cWe31rQ35+PhQUFGr8dYqKisDj8SAlRR/vhkJWVpbvs1Ndnz9/hp+fH5ydnREYGMi3beLEiSgtLRXbaxHyIw2m60cUDMMgICAA1tbWkJeXh4qKCkaOHIlXr17x5YuOjsbQoUOhq6sLOTk5tGrVCjNmzEBGRgabZ8WKFVi0aBEAwMDAgO3iKGsu1dfXx6BBg3D69GnY2NhAXl4eZmZmOH36NICvXVdmZmZQVFRE+/btBbqw7ty5AycnJ+jr60NeXh76+voYO3Ys3rx5w5evrAssOjoakyZNgqqqKhQVFTF48GCB4yrP1atX0bt3bzRu3BgKCgqws7NDZGQk37Hq6uoCANzd3cHj8UTuVlNSUoK0tDRfmrCun6ysLMyePRvNmjWDjIwMDA0N4enpiYKCAr58PB4Pc+fOxb59+2BiYgJ5eXnY2trixo0bYBgGvr6+MDAwQKNGjdCrVy+8ePGCb//KvNdl54DH4+HevXsYOXIkVFRUKvy1e+3aNairq2PQoEHIy8sDADx//hy//PILNDU1ISsrCzMzM/j7+/PtV9ZNc/DgQSxYsADNmjWDrKwsXrx4gfz8fCxcuBAGBgaQk5ODqqoqbG1tERISUqlzL+x4YmNjMXbsWCgpKUFLSwuTJ09Gbm4um8/GxgZdu3YV2L+kpATNmjXDiBEj2LTKvmffSk9Ph4yMDJYtWyawLT4+HjweD1u3bmXT3r9/jxkzZkBXVxcyMjIwMDDAypUrUVxczOYp6zbx9fXFunXr2M9Qjx498OzZMxQVFWHx4sVo2rQplJSUMHz4cKSlpQm8flhYGDp16gRFRUU0atQI/fr1w/379/nylHX7vXjxAg4ODmjUqBGaN2+OBQsWsMedkJAADQ0NAMDKlSvZ74qKur+Edf1U9j0TJi8vDwUFBdDR0RG6XUJColJ1ffHiBSZNmgQjIyMoKCigWbNmGDx4MB4/fixQZmxsLOzt7aGgoAANDQ3MmTMHkZGRfN+TZf755x/07t0bTZo0gYKCAjp37ozz58/z5UlPT8f06dPRvHlzyMrKQkNDA507d8Y///xT4bGTOohpIPbt28cAYG7cuMEUFRUJPBwcHBg9PT2+faZNm8ZIS0szCxYsYM6cOcMEBwczpqamjJaWFvP+/Xs2344dOxgfHx/m1KlTzOXLl5kDBw4wbdq0YUxMTJjCwkKGYRgmKSmJ+fXXXxkATHh4OHP9+nXm+vXrTG5uLsMwDKOnp8fo6uoyFhYWTEhICBMVFcV06NCBkZaWZpYvX8507tyZCQ8PZ06cOMEYGxszWlpaTH5+PluHo0ePMsuXL2dOnDjBXL58mQkNDWW6d+/OaGhoMOnp6QLnoXnz5szkyZOZv//+mwkMDGQ0NTWZ5s2bM9nZ2RWex0uXLjHS0tJMu3btmLCwMCYiIoKxt7dneDweExoayh5reHg4A4D59ddfmevXrzP37t0rt8zXr18zAJh169ax70dGRgazZ88eRkpKivH09OTL3717d6Z79+7s88+fPzNWVlaMoqIis2HDBubcuXPMsmXLGCkpKcbBwYFvXwCMnp4eY2dnx3c+VVVVmXnz5jFDhw5lTp8+zRw+fJjR0tJirKysmNLS0iq91wzDMF5eXuxrubu7M9HR0UxERATDMAzj4uLCKCoqsnnDwsIYWVlZZtasWUxxcTHDMAwTGxvLKCkpMZaWlkxQUBBz7tw5ZsGCBYyEhASzYsUKdt+LFy8yAJhmzZoxI0eOZE6dOsWcPn2ayczMZGbMmMEoKCgwmzZtYi5evMicPn2aWbt2LbNt27YK3+OyMo8ePSpwPCYmJszy5cuZ6OhoZtOmTYysrCwzadIkNt+WLVsYAMyzZ8/4yoyKimIAMKdOnarye6anp8e4uLiwz4cPH840b96cKSkp4cvn5ubGyMjIMBkZGQzDMExKSgrTvHlzRk9Pj/nzzz+Zf/75h1m1ahUjKyvLTJw4kd2v7PrT09NjBg8ezJw+fZo5dOgQo6WlxRgbGzMTJkxgPys7d+5kGjVqxAwePJjvtVevXs3weDxm8uTJzOnTp5nw8HCmU6dOjKKiIhMbG8vmc3FxYWRkZBgzMzNmw4YNzD///MMsX76c4fF4zMqVKxmGYZgvX74wZ86cYQAwU6ZMYb8rXrx4Ue57VnYM+/btq/J7Vp5WrVoxjRs3ZjZu3Mj8999/fJ+DMj+q6+XLl5kFCxYwx44dYy5fvsycOHGCGTZsGCMvL8/Ex8ez5bx7945RU1NjWrRowezfv5+JiopiJkyYwOjr6zMAmIsXL7J5Dx48yPB4PGbYsGFMeHg489dffzGDBg1iJCUlmX/++YfN169fP0ZDQ4MJDAxkLl26xERERDDLly9nv6dI/dHgApWKHt8GKtevX2cAMBs3buQrJykpiZGXl2fc3NyEvk5paSlTVFTEvHnzhgHAnDx5kt3m6+vLAGBev34tsJ+enh4jLy/PJCcns2kPHjxgADA6OjpMXl4emx4REcH3pS9McXEx8+nTJ0ZRUZHZsmWLwHkYPnw4X/5r164xAJg//vij3DIZhmE6duzIaGpqMh8/fuR7LQsLC0ZXV5f9Miv74vT19a2wvG/zCntMnDiR/edd5vtAZefOnQwA5siRI3z51q1bxwBgzp07x6YBYLS1tZlPnz6xaWXn09ramu/L2M/PjwHAPHr0SGi9K3qvy/5JLF++XGC/bwOVtWvXMpKSksy6dev48vTr14/R1dVlA9kyc+fOZeTk5JisrCyGYf4/qOjWrZvA61hYWDDDhg0TWveKVBSorF+/ni/v7NmzGTk5Ofa8ZWRkMDIyMsySJUv48o0ePZrR0tJiioqKGIap2nv2faBy6tQpgTzFxcVM06ZNGUdHRzZtxowZTKNGjZg3b97wvcaGDRsYAGwAUXb9tWnThi/4KXv/hwwZwre/q6srA4B9bxITExkpKSnm119/5cv38eNHRltbmxk9ejSb5uLiIvS4HRwcGBMTE/Z5eno6A4Dx8vJiKqOiQOVH71l5bt26xbRo0YL9LDZu3JgZNGgQExQUxLdvVepaXFzMFBYWMkZGRsy8efPY9EWLFjE8Ho8vqGOYr5+DbwOVvLw8RlVVVSBQLCkpYdq0acO0b9+eTWvUqBHj6ur6wzqRuq/Bdf0EBQXh9u3bAo8uXbrw5Tt9+jR4PB7Gjx+P4uJi9qGtrY02bdrwNUWmpaVh5syZaN68OaSkpCAtLQ09PT0AwH///VfpullbW6NZs2bsczMzMwBfuzq+HdtQlv5tt86nT5/g7u6OVq1aQUpKClJSUmjUqBHy8vKE1mHcuHF8z+3s7KCnp4eLFy+WW7+8vDzcvHkTI0eO5Ju1IikpiQkTJiA5ORlPnz6t9PF+7/fff2ffj4sXL2LNmjU4cuQIxo4dW+F+Fy5cgKKiIkaOHMmXXtb8/H2TcM+ePaGoqMg+LzufAwYMAI/HE0j/9jxX9b12dHQUWmeGYTBjxgx4eXkhODgYbm5u7LYvX77g/PnzGD58OBQUFPiuPwcHB3z58gU3btz44eu0b98ef//9NxYvXoxLly7h8+fPQutSFUOGDOF7bmVlhS9fvrBdIWpqahg8eDAOHDjAjmPIzs7GyZMn4ezszI6bqep79q0BAwZAW1sb+/btY9POnj2Ld+/eYfLkyWza6dOn0bNnTzRt2pTvHA4YMAAAcPnyZb5yHRwc2C4N4P/f/4EDB/LlK0tPTExkX7u4uBjOzs58ryMnJ4fu3bsLdFvweDwMHjxY4Dx+300rLj96z8rzv//9Dy9evMCZM2ewZMkSdOrUCefPn4ezszOGDBkChmF++NrFxcVYs2YNWrduDRkZGUhJSUFGRgbPnz/n+7xcvnwZFhYWAuPPvv/sx8TEICsrCy4uLnznurS0FP3798ft27fZrtP27dtj//79+OOPP3Djxg0UFRX9sL6kbmpwo+3MzMxga2srkK6kpISkpCT2eWpqKhiGgZaWltByDA0NAQClpaWwt7fHu3fvsGzZMlhaWkJRURGlpaXo2LFjlf45qKqq8j2XkZGpMP3Lly9s2i+//ILz589j2bJl+N///ocmTZqAx+PBwcFBaB20tbWFpmVmZpZbv+zsbDAMI7TfumnTpgBQ4f4/oqury/fe9OjRAzweDx4eHjh79iz69esndL/MzExoa2vzBRkAoKmpCSkpKYE6iXqeRXmvy+vjLywsRFhYGMzNzdl/nN8eT3FxMbZt24Zt27YJ3f/7MTHCXmfr1q3Q1dVFWFgY1q1bBzk5OfTr1w++vr4wMjISWu6PqKmp8T2XlZUFAL5jnzx5Mo4fP47o6Gj069cPISEhKCgo4BtjUdX37FtSUlKYMGECtm3bhpycHCgrK2P//v3Q0dHhu0ZSU1Px119/CYxxKvP9ORT1ukhNTQXw9R+7MN8GPwCgoKAAOTk5vjRZWVm+z7M4VeY9K4+0tDT69evHntfMzEyMHDkSp0+fxt9//w0HB4cK958/fz78/f3h7u6O7t27Q0VFBRISEpg6dSrf62dmZsLAwEBg/++/f8vO9fcB7reysrKgqKiIsLAw/PHHH9i9ezeWLVuGRo0aYfjw4Vi/fr3Q7z9SdzW4QKWy1NXVwePxcOXKFfaD/a2ytCdPnuDhw4fYv38/XFxc2O3fD8KsSbm5uTh9+jS8vLywePFiNr2goABZWVlC93n//r3QtFatWpX7OmVfMikpKQLb3r17B+DreRMnKysrAMDDhw/LDVTU1NRw8+ZNMAzD948vLS0NxcXFYquTKO/19/+Iy8jKyuLixYvo168f+vTpgzNnzkBFRQXA1/Nc1ko1Z84coft//6Uu7HUUFRWxcuVKrFy5EqmpqWzryuDBgxEfH//D4xVVv3790LRpU+zbtw/9+vXDvn370KFDB75fy9V9zyZNmgRfX1+EhoZizJgxOHXqFFxdXSEpKcnmUVdXh5WVFVavXi20jLLgurrK6nrs2DG2de1npaamBldXV1y6dAlPnjz5YaBy6NAhODs7Y82aNXzpGRkZfEs0qKmpsUHIt77/nio719u2bSt3llNZcKOurg4/Pz/4+fkhMTERp06dwuLFi5GWloYzZ8788FhJ3UGBSjkGDRqEtWvX4u3btxg9enS5+cq+ZL8PZv7880+BvFX5JVMVPB4PDMMI1GH37t0Ci6WVOXz4MF93QUxMDN68ecOujSCMoqIiOnTogPDwcGzYsAHy8vIAvrY0HDp0CLq6ujA2NhbDEf2/Bw8eAPj6S7s8vXv3xpEjRxAREYHhw4ez6UFBQex2cajKe10ZNjY2uHz5Mvr06YMePXogOjoampqaUFBQQM+ePXH//n1YWVmxv+CrQ0tLCxMnTsTDhw/h5+dXo1Oly4IsPz8/XLlyBXfu3BE4R9V9z8zMzNChQwfs27cPJSUlKCgowKRJk/jyDBo0CFFRUWjZsiUbBNaEfv36QUpKCi9fviy3q6+qauq7orKKiorw4cMHgdYY4P+7OMsCvYrqyuPxBD4vkZGRePv2Ld+Pou7du2PDhg2Ii4vjC2hDQ0P59u3cuTOUlZURFxeHuXPnVvp4WrRogblz5+L8+fO4du1apfcjdQMFKuXo3Lkzpk+fjkmTJuHOnTvo1q0bFBUVkZKSgqtXr8LS0hKzZs2CqakpWrZsicWLF4NhGKiqquKvv/5CdHS0QJmWlpYAgC1btsDFxQXS0tIwMTFB48aNq1XXJk2aoFu3bvD19YW6ujr09fVx+fJl7Nmzp9yF5e7cuYOpU6di1KhRSEpKgqenJ5o1a4bZs2dX+Fo+Pj7o27cvevbsiYULF0JGRgYBAQHseifltSBURmJiIjv2Ii8vD9evX4ePjw/09PT4prV+z9nZGf7+/nBxcUFCQgIsLS1x9epVrFmzBg4ODujTp4/IdfpWVd7ryjIzM8OVK1fQp08fdOvWDf/88w90dXWxZcsWdOnSBV27dsWsWbOgr6+Pjx8/4sWLF/jrr79w4cKFH5bdoUMHDBo0CFZWVlBRUcF///2HgwcPolOnTjW+nsvkyZOxbt06/PLLL5CXl8eYMWP4tovjPZs8eTJmzJiBd+/ewc7OTmBBQW9vb0RHR8POzg6//fYbTExM8OXLFyQkJCAqKgo7d+5kp9FXh76+Pry9veHp6YlXr16hf//+UFFRQWpqKm7dusW2bFVF48aNoaenh5MnT6J3795QVVVlP9u1ITc3F/r6+hg1ahT69OmD5s2b49OnT7h06RK2bNkCMzMz9jNZUV0HDRqE/fv3w9TUFFZWVrh79y58fX0Fzrurqyv27t2LAQMGwNvbG1paWggODmZb/sq6zxo1aoRt27bBxcUFWVlZGDlyJDQ1NZGeno6HDx8iPT0dO3bsQG5uLnr27IlffvkFpqamaNy4MW7fvo0zZ85U+F1C6ijOhvHWsrLZLrdv3xa6feDAgQLTkxmGYfbu3ct06NCBUVRUZOTl5ZmWLVsyzs7OzJ07d9g8cXFxTN++fZnGjRszKioqzKhRo5jExEShI+E9PDyYpk2bMhISEnyj2fX09JiBAwcKvD4AZs6cOXxpwmbUJCcnM46OjoyKigrTuHFjpn///syTJ08EZkyUnYdz584xEyZMYJSVlRl5eXnGwcGBef78+Q/O4ldXrlxhevXqxZ6Tjh07Mn/99dcP61geYbN+5OTkGGNjY8bV1ZVJSUnhy//9rB+GYZjMzExm5syZjI6ODiMlJcXo6ekxHh4ezJcvX/jyVfZ8Mozw2S+Vfa/LZlx8OzW8zPfTkxnm6/tnamrK6OvrMy9fvmTrNXnyZKZZs2aMtLQ0o6GhwdjZ2fHNzBJWxzKLFy9mbG1tGRUVFUZWVpYxNDRk5s2bx07fLU9Fs36+P56y60nYTDY7OzsGADNu3Dihr1PZ9+z7a7hMbm4uIy8vzwBgdu3aJfQ10tPTmd9++40xMDBgpKWlGVVVVaZdu3aMp6cnO/OrKu//t8f8/XdJREQE07NnT6ZJkyaMrKwso6enx4wcOZJvyqyw955h/v/8fuuff/5hbGxsGFlZWQaA0HNQpqJZP1V5z8oUFBQwGzZsYAYMGMC0aNGCkZWVZeTk5BgzMzPGzc2NyczMrFRds7OzmSlTpjCampqMgoIC06VLF+bKlStCP8NPnjxh+vTpw8jJyTGqqqrMlClTmAMHDjAAmIcPH/LlvXz5MjNw4EBGVVWVkZaWZpo1a8YMHDiQfa++fPnCzJw5k7GysmKaNGnCyMvLMyYmJoyXlxffDEpSP/AYphJDt8lPY//+/Zg0aRJu374tdFAxIYTUFdOnT0dISAgyMzPF0gVK6ifq+iGEEMI5b29vNG3aFIaGhvj06RNOnz6N3bt3Y+nSpRSkNHAUqBBCCOGctLQ0fH19kZycjOLiYhgZGWHTpk34/fffua4a4Rh1/RBCCCGkzmpwK9MSQgghpP6gQIUQQgghdRYFKoQQQgipsyhQIYQQQkidRbN+CCGEkBpWnVW7v9UQ57/8tIFKTk4O11Wot5SVlen8iUhZWRkDBw7kuhr1UmRkZI3dQbghkJOTQ2FhIdfVqJdonZa67acNVAghhJC6QlwtKg0RBSqEEEJIDaNARXQUqBBCCCE1jAIV0dGsH0IIIYTUWdSiQgghhNQwCQlqFxAVBSqEEEJIDaOuH9FRiEcIIYSQOotaVAghhJAaRi0qoqNAhRBCCKlhFKiIjrp+CCGEEFJnUYsKIYQQUsOoRUV0dSZQSUpKQkJCAvLz86GhoQFzc3PIyspyXS1CCCGk2ihQER2ngcqbN2+wc+dOhISEICkpie+ukDIyMujatSumT58OR0dHmoNOCCGENECc/ff//fffYWlpiefPn8Pb2xuxsbHIzc1FYWEh3r9/j6ioKHTp0gXLli2DlZUVbt++zVVVCSGEkGqRkJAQy6Mh4qxFRUZGBi9fvoSGhobANk1NTfTq1Qu9evWCl5cXoqKi8ObNG/zvf//joKaEEEJI9VDXj+g4C1R8fX0rndfBwaEGa0IIIYTULApURMf5YNqbN2/i1KlTKCoqQp8+fWBvb891lQghhBBSR3AaqJw4cQKjRo2CnJwcpKSksHHjRmzcuBGurq5cVosQQggRK2pRER2nI3PWrFmDiRMnIicnBzk5OVi5ciX++OMPLqtECCGEiB2PxxPLoyHiNFB5+vQp3NzcICX1tWFn0aJFyMnJQUZGBpfVIoQQQkgdwWmg8unTJygrK7PPZWVlIS8vjw8fPnBXKUIIIUTMqEVFdJwPpj179iyUlJTY56WlpTh//jyePHnCpg0ZMoSLqhFCCCFi0VDXQBEHzgMVFxcXgbQZM2awf/N4PJSUlNRmlQghhBBSR3AaqJSWlnL58oQQQkitaKjdNuLAeYsKIYQQ8rOjQEV0nAYqp06dqlQ+GqNCCCGENEycBirDhg37YR4ao0IIIaS+oxYV0dEYFUIIIaSGUaAiOhqjQgghhNQwClRER2NUqunYsWM4dOgQMjMzYWBggHnz5sHGxkZo3oyMDGzZsgXx8fFISkrC6NGjMX/+/HLLPnfuHJYtW4Zu3bpV6W7T9QmdP9ENHDgQI0aMgKqqKhITExEYGIjY2FiheS0tLbF27VqB9BkzZiA5ORkA0KdPH8ybN08gz7Bhw1BUVCTeyteysLAw7N+/HxkZGWjZsiXc3NzQtm3bcvPfuXMHGzZswMuXL6GhoYGJEydi9OjR7PaioiLs2bMHf/31F9LS0qCvrw9XV1d07tyZzbNjxw7s3LmTr1w1NTVcuHBB/AdYg0JDQ7F//36kp6ejZcuWcHd3R7t27crNf/v2bfj6+rLnbvLkyXznDgAOHjyII0eOICUlBcrKyujbty9cXV0hKysrUN7u3buxZcsWjB8/Hu7u7mI/PlL31akxKjweDwzDCKTV1TEq0dHR2Lx5M9zc3GBlZYUTJ05g3rx5CA0Nhba2tkD+wsJCKCsrY9KkSQgJCamw7JSUFGzduhXW1tY1VHvu0fkTXdeuXTFt2jQEBATgv//+Q//+/bFy5UrMmjUL6enp5e43bdo0fP78mX2em5vLtz0vL49vHSMA9T5IOXPmDNavXw9PT09YW1vj2LFjmD17Nk6cOAEdHR2B/MnJyZgzZw4cHR2xZs0aPHjwAKtXr4aqqir69OkDANi+fTsiIyPh5eUFAwMDxMTEYN68eThw4ADMzMzYslq2bInAwED2eX1b9OvMmTNYt24dli5dChsbGxw9ehSzZs3CyZMnf3ju1q5di/v37+OPP/6AiooK+vbtCwA4ffo0/Pz84O3tDWtra7x58wZLly4FAIFA5MmTJzh27BiMjY1r/mBrWH177+sSTs9caWkp30NBQQEvXrzgS6urQQoAhISEYMiQIRg6dCgMDAwwf/58aGlp4fjx40LzN23aFAsWLICDgwMaNWpUbrklJSVYvnw5pk+fjmbNmtVU9TlH5090w4cPx7lz53Du3DkkJSVh165dyMjIgIODQ4X75ebmIjs7m318P06MYRi+7dnZ2TV5GLXi4MGDGD58OEaMGAFDQ0O4ublBW1sbR44cEZr/6NGj0NHRgZubGwwNDTFixAgMGzYMBw4cYPNERkZi6tSp6Nq1K3R1dTF69GjY2dkhKCiIrywpKSmoq6uzD1VV1Ro9VnELCgrCiBEj4OjoCENDQ7i7u0NbWxthYWFC8x85cgTa2tpwd3eHoaEhHB0dMXz4cOzfv5/N8/DhQ9jY2GDgwIFo1qwZ7OzsMGDAAMTFxfGVlZ+fj8WLF8PLywtNmjSpycOsFbSEvugoxBNRUVER4uPj0aFDB7709u3b4/Hjx9Uqe8+ePVBRUanTXV7VRedPdFJSUmjVqhXu37/Pl37v3j2+X/PCbN26FQcPHsTq1athZWUlsF1eXh779u3DgQMH4OXlBUNDQ7HWvbYVFRXhv//+Q6dOnfjSO3XqhIcPHwrd59GjRwL57ezsEBcXx7YuFRYWQkZGhi+PrKwsHjx4wJf25s0b9OnTBwMGDICbmxvbzVYfFBUVIS4uDnZ2dnzpdnZ2AsdZ5uHDhwL5O3fuzHfu2rZti7i4OPZznpSUhCtXrqBr1658+61evRpdu3YVeC9Iw0ODaUWUk5ODkpISgV9IampquHHjhsjlPnz4EKdOncKhQ4eqW8U6jc6f6Jo0aQJJSUnk5OTwpefk5EBFRUXoPllZWdi6dStevHgBaWlp9OrVC6tXr8bixYvZcS1JSUnYvHkzEhISoKCggCFDhsDX1xe//vor3r17V9OHVSOys7NRUlICNTU1vnQ1NbVy79KekZEhNH9xcTFycnKgoaEBOzs7HDx4EO3atUPz5s1x8+ZNXLp0ia8F2NLSEqtXr4aenh4yMzOxa9cuODs7Izw8nO9mrHVVRecuMzNT6D6ZmZk/PHcDBgxAVlYWnJ2dAQDFxcUYM2YMpk6dyu7z999/Iy4uDqGhoWI+Ku401NYQcaj3gUpBQQEKCgr40oQNyKop3198DMOIfEHm5eXBy8sLS5YsqRdfZOJA5090wsZzfZ9W5u3bt3j79i37PD4+Hurq6nB0dGQDladPn+Lp06dsnri4OGzduhWDBw/Gn3/+WQNHUHuqep0Jy/9tupubG7y9vTFs2DDweDzo6upi6NChOHnyJLtPly5d2L+NjIxgZWWFQYMG4dSpU+w/6fqovGuszI/O3e3bt7Fr1y4sXboUlpaWSEpKwtq1a6Guro6ZM2fi/fv3WLt2LQIDA2v1u7ymUaAiujoVqIjSB+fj44OVK1fypXl5ecHV1VWMNROkrKwMSUlJgV8WWVlZIvdDv337FikpKVi4cCGbVjaGwM7ODkeOHIGurq7ola5D6PyJ7sOHDygpKRFoPVFSUhJoZanI06dP0bNnz3K3MwyDZ8+eoWnTpqJWlXMqKiqQlJQUaD3JysoS+OVfRl1dXWh+KSkp9k7vqqqq8PPzQ0FBAXJycqCpqQk/P78Kz5WCggKMjIyQmJhYzaOqHWXnTthntLxzJ6yl6vtzt337dgwePBiOjo4AAGNjY+Tn58Pb2xvTp09HbGwssrKyMGbMGLaMkpIS3L17FyEhIbh79y4kJSXFeaikjuM0UFFRUeELTD59+gQbGxuB0dFZWVnlluHh4SEwRVVWVpZvZkNNkJaWhqmpKW7duoUePXqw6bdu3UK3bt1EKlNPTw/BwcF8aTt37kR+fj470PRnQedPdMXFxXjx4gVsbGxw/fp1Nt3GxqZK3WaGhoYVfrbK8iQkJIhaVc5JS0vDzMwMN27cQO/evdn0Gzdu8F1337KyssK///7Ll3b9+nW0bt0a0tLSfOmysrLQ0tJCUVERzp8/D3t7+3LrUlhYiFevXpU7/b6ukZaWRuvWrXH9+nW+c3f9+vVyA9w2bdrg8uXLfGkxMTF85+7z588CP0glJSXBMAwYhkHHjh0RHh7Ot33ZsmUwMDDA5MmT622QQi0qouM0UPHz86t2GbKyskKbB2s6UAGAsWPHYsWKFTA1NYWlpSUiIiKQmpqKESNGAAD8/f2Rnp6OFStWsPs8e/YMwNcR7Tk5OXj27BmkpKRgaGgIWVlZtGzZku81GjduDAAC6T8DOn+iO3HiBBYsWIDnz58jPj4e/fv3h4aGBqKiogAALi4uUFNTw6ZNmwAAQ4cORWpqKhITEyElJYWePXuiS5cuWL16NVvm2LFj8fTpU7x79w4KCgoYPHgwDA0NsWPHDk6OUVwmTJgAT09PtG7dGm3atMHx48eRkpKCUaNGAQC2bNmCtLQ09lyMGjUKoaGh8PX1haOjIx4+fIgTJ05g3bp1bJmPHj1CWloaTE1NkZaWhh07dqC0tBQTJ05k82zcuBHdu3eHtrY2srKysGvXLuTl5dWrQd7Ozs7w8PCAubk52rRpg6NHjyIlJYVdF8XPzw9paWlYs2YNAGD06NEIDQ3F+vXrMXLkSDx8+BDh4eFYv349W2aPHj0QFBQEMzMzWFpaIjExEdu3b0ePHj0gKSkJRUVFGBkZ8dVDXl4eysrKAun1CQUqouM0UHFxceHy5autb9++yM3Nxd69e5GRkQFDQ0Ns3ryZXV8gMzMTqampfPtMmDCB/Ts+Ph5nz56Fjo4OIiIiarPqdQKdP9FduXIFTZo0wdixY6Gqqoo3b97Ay8uLXUNFVVUVGhoabH4pKSlMmTIFampqKCwsZPPfuXOHzdOoUSP8+uuvUFFRQV5eHl6+fAl3d3c2OKyv+vfvj9zcXAQGBiI9PR2tWrWCv78/202TkZGB9+/fs/l1dXXh7+8PX19fhIWFQUNDA+7u7uwaKsDX1hF/f38kJydDQUGBDfq+nUabmpqKxYsXIzs7GyoqKrCyssLBgwfrVVda//79kZOTg507d7LnLiAggD2G9PR0pKSksPm/PXehoaHQ1NSEh4cHu4YKAEyfPh08Hg/btm1DWloaVFRU0L17d/z222+1fnykfuAxPxoZVUdUdZBlVfrqCT9lZWU6fyJSVlbGwIEDua5GvRQZGYkvX75wXY16S05ODoWFhVxXo176fqp5TRDXVP9Xr16JpZz6hLN1VMzMzBAcHPzDD9bz588xa9YsvmZXQgghpD6hBd9Ex1nXj7+/P9zd3TFnzhzY29vD1tYWTZs2hZycHLKzsxEXF4erV68iLi4Oc+fOxezZs7mqKiGEEFItDTXIEAfOApVevXrh9u3biImJQVhYGIKDg5GQkIDPnz9DXV0dNjY2cHZ2xvjx4xvEmhiEEEIIEcT5Oip2dnYCSy4TQgghPxNqUREd54EKIYQQ8rOjQEV0nAYqW7durVQ+mrZGCCGENEycBiqbN2/+YR4ej0eBCiGEkHrt+xXXSeVxGqi8fv2ay5cnhBBCagV1/YiOQjxCCCHkJxYQEAADAwPIycmhXbt2uHLlSoX5CwoK4OnpCT09PfbWJHv37q2l2griNFC5cOECWrdujQ8fPghsy83Nhbm5ucDNwQghhJD6hqsF38LCwuDq6gpPT0/cv38fXbt2xYABAyq8i/fo0aNx/vx57NmzB0+fPkVISAhMTU2rc/jVwvlNCadNm8Z3f4wySkpKmDFjBjZv3izy3XQJIYSQuoCrMSqbNm3ClClTMHXqVABf/++ePXsWO3bsgI+Pj0D+M2fO4PLly3j16hVUVVUBAPr6+rVZZQGctqg8fPgQ/fv3L3e7vb097t69W4s1IoQQQn4OhYWFuHv3Luzt7fnS7e3tERMTI3SfU6dOwdbWFuvXr0ezZs1gbGyMhQsX4vPnz7VRZaE4bVFJTU2FtLR0udulpKTYu8ESQggh9ZW4BtMWFBSgoKCAL01WVhaysrICeTMyMlBSUgItLS2+dC0tLb47hn/r1atXuHr1KuTk5HDixAlkZGRg9uzZyMrK4mycCqctKs2aNcPjx4/L3f7o0SPo6OjUYo0IIYQQ8ZOQkBDLw8fHB0pKSnwPYV043/o+SGIYptzAqbS0FDweD4cPH0b79u3h4OCATZs2Yf/+/Zy1qnAaqDg4OGD58uVCb+3++fNneHl5YdCgQRzUjBBCCBEfcQ2m9fDwQG5uLt/Dw8ND6Guqq6tDUlJSoPUkLS1NoJWljI6ODpo1awYlJSU2zczMDAzDIDk5WXwnpAo4DVSWLl2KrKwsGBsbY/369Th58iROnTqFdevWwcTEBFlZWfD09OSyioQQQkidISsriyZNmvA9hHX7AICMjAzatWuH6OhovvTo6Ohy77HXuXNnvHv3Dp8+fWLTnj17BgkJCejq6orvQKqA0zEqWlpaiImJwaxZs+Dh4QGGYQB8jTz79euHgICAcqM+QgghpL7gasG3+fPnY8KECbC1tUWnTp0QGBiIxMREzJw5EwDg4eGBt2/fIigoCADwyy+/YNWqVZg0aRJWrlyJjIwMLFq0CJMnT4a8vDwnx8BpoPLq1SsYGBggKioK2dnZePHiBRiGgZGREVRUVLisGiGEECI2XE1PHjNmDDIzM+Ht7Y2UlBRYWFggKioKenp6AICUlBS+NVUaNWqE6Oho/Prrr7C1tYWamhpGjx6NP/74g5P6AwCPKWvG4ICkpCRSUlKgqakJ4OsJ3bp1q1haUXJycqpdRkOlrKxM509EysrKGDhwINfVqJciIyOFjlcjlSMnJ4fCwkKuq1EvycjI1PhrtGvXTizlNMQlOzgdo/J9jBQVFYW8vDyOakMIIYTUDK5Wpv0ZcNr1QwghhDQEdPdk0XF65oRFiA01YiSEEEKIIE5bVBiGwcSJE9mpVV++fMHMmTOhqKjIly88PJyL6hFCCCFiQT/CRcdpoOLi4sL3fPz48RzVhBBCCKk51PUjOk4DlX379nH58oQQQgip42gwLSGEEFLDqOtHdBSoEEIIITWMAhXRUaBCCCGE1DAaoyI6OnOEEEIIqbOoRYUQQgipYdT1IzoKVAghhJAaRl0/oqMzRwghhJA666dtUVFWVua6CvUanT/RRUZGcl2FektOTo7rKtRrtXEXYCIa6voR3U8bqNDt4kVH/yyqh6490cjJydGXeTUwDIPS0lKuq1Ev1Ua3DF3boqOuH0IIIYTUWT9tiwohhBBSV9BgWtFRoEIIIYTUMOr6ER2FeIQQQgips6hFhRBCCKlh1PUjOgpUCCGEkBpGXT+io0CFEEIIqWEUqIiuzgQqSUlJSEhIQH5+PjQ0NGBubg5ZWVmuq0UIIYQQDnEaqLx58wY7d+5ESEgIkpKSwDAMu01GRgZdu3bF9OnT4ejoSP17hBBC6i36HyY6zs7c77//DktLSzx//hze3t6IjY1Fbm4uCgsL8f79e0RFRaFLly5YtmwZrKyscPv2ba6qSgghhFQLj8cTy6Mh4qxFRUZGBi9fvoSGhobANk1NTfTq1Qu9evWCl5cXoqKi8ObNG/zvf//joKaEEEII4QpngYqvr2+l8zo4ONRgTQghhJCaRV0/ouN8MO3Nmzdx6tQpFBUVoU+fPrC3t+e6SoQQQohYNdRuG3HgNFA5ceIERo0aBTk5OUhJSWHjxo3YuHEjXF1duawWIYQQQuoITtui1qxZg4kTJyInJwc5OTlYuXIl/vjjDy6rRAghhIidhISEWB4NEadH/fTpU7i5uUFK6mvDzqJFi5CTk4OMjAwuq0UIIYSIFc36ER2ngcqnT5+grKzMPpeVlYW8vDw+fPjAXaUIIYQQUmdwPpj27NmzUFJSYp+Xlpbi/PnzePLkCZs2ZMgQLqpGCCGEiEVDbQ0RB84DFRcXF4G0GTNmsH/zeDyUlJTUZpUIIYQQsaJARXScBiqlpaVcvjwhhBBSKyhQEV3DHEJMCCGEkHqB0xaVU6dOVSofjVEhhBBSn1GLiug4DVSGDRv2wzw0RoUQQkh9R4GK6GiMCiGEEELqLM5n/RBCCCE/O2pRER2NUSGEEEJqGAUqoqtTY1R4PB4YhhFIozEqhBBCSMNUp8aoNG7cGA8fPoShoSFHNSKEEELEr6HeUFAcaIwKIYQQUsOo60d0FOJVUVhYGAYMGID//e9/cHJywr179yrMf+fOHTg5OeF///sfHBwccOTIEb7tRUVF2LlzJwYOHIj//e9/GDVqFK5du8aXZ8eOHWjTpg3fo1evXmI/tpp2+PBh9OrVC5aWlhgxYgTu3LlTYf5bt25hxIgRsLS0RO/evRESEsK3PTw8HCYmJgKPgoICNk9wcDAGDx6Mtm3bom3bthgzZgwuX75cI8dX0+jaE59Zs2bh1atX+Pz5M+7cuYMuXbpUmF9GRgZ//PEHEhIS8OXLF7x48QKTJk1it0tJSWHZsmV48eIFPn/+jAcPHqBfv341fRi1Ijg4GH369EGbNm3g6OhYqc+to6Mj2rRpg759+yI0NFQgz4cPH+Dt7Y2uXbuiTZs2GDhwIN/nMiQkBEOHDoWtrS1sbW3h5OSEf//9V+zHRuoHalGpgjNnzmD9+vXw9PSEtbU1jh07htmzZ+PEiRPQ0dERyJ+cnIw5c+bA0dERa9aswYMHD7B69WqoqqqiT58+AIDt27cjMjISXl5eMDAwQExMDObNm4cDBw7AzMyMLatly5YIDAxkn9e3ZsSoqCj4+PjAy8sLbdu2RWhoKKZNm4bIyEg0bdpUIH9SUhKmT5+OUaNGwdfXF/fu3cPKlSuhqqrK9w+gUaNGOHPmDN++srKy7N/a2tpYuHAhWrRoAQCIiIjAnDlzcOLECRgZGdXQ0YofXXviM3r0aPj5+WH27Nm4du0aZsyYgb///hutW7dGUlKS0H2OHDkCLS0tTJkyBS9evICmpiakpP7/6/OPP/7A+PHjMW3aNMTHx6Nfv344ceIE7Ozs8ODBg1o6MvGLiorC2rVrsWzZMrRt2xZhYWGYMWMG/vrrL6Gf2+TkZMycORMjR47E+vXrce/ePaxatQqqqqqwt7cHABQWFmLKlClQVVXFli1boKWlhffv30NRUZEtR1tbG/Pnz2c/tydPnsTcuXNx/PjxevW5/Ra1qIiOx3w/epVDTZo0wcOHD2FgYFDtsr58+SKGGvEbN24czMzMsHTpUjZt2LBh6NmzJ37//XeB/Js3b8bly5cRERHBpq1atQrPnj3DwYMHAQB9+vTB1KlT4eTkxOZxdXWFvLw8fHx8AHz9VXvx4kWBX8Q1RU5OTuxljho1Cq1bt8bKlSvZtAEDBqBPnz5YsGCBQH5fX19cuHABf//9N5u2fPlyPH36FGFhYQC+tqisWbPmh7/wvte+fXssWrQIo0aNEvFoKkbXnmjk5ORq5cv8xo0buHfvHmbPns2mxcXFISIiAkuWLBHI369fP4SGhsLQ0BDZ2dlCy3z79i1Wr16NgIAANu3EiRP49OkTJkyYIP6DEIJhGLGvTTVmzBiYmZlhxYoVbNrAgQPRu3dvzJ8/XyD/hg0bcPHiRURGRrJpK1asQHx8PNuyEhoair179yIyMhLS0tKVrkvHjh2xcOFCjBw5UvQDKkdtBN/Ozs5iKScoKEgs5dQnnP40UlFRgaqqKvv49OkTbGxs+NJUVVW5rCKrqKgI//33Hzp16sSX3qlTJzx8+FDoPo8ePRLIb2dnh7i4OBQVFQH4+utCRkaGL4+srKzAr7A3b96gT58+GDBgANzc3JCcnFzNI6o9hYWFiI2NFWhe79y5M+7fvy90nwcPHqBz5858aV27dsWTJ0/YcwcA+fn56NmzJ7p164YZM2YgLi6u3HqUlJQgMjIS+fn5sLGxqcYR1S669sRHWloa7dq1w7lz5/jSz507Bzs7O6H7DBkyBHfu3GGP/enTp/D19eUL6GVlZQUC1M+fP/+wS6kuK/vcfv85rOrntnPnzoiNjWWvuwsXLsDa2hqrVq1Cly5dMHjwYPz555/lzu789nNrbW1d/QPjCI/HE8tDFAEBATAwMICcnBzatWuHK1eulJv30qVLQl83Pj5e1EOvNk67fvz8/Lh8+SrJzs5GSUkJ1NTU+NLV1NSQkZEhdJ+MjAyh+YuLi5GTkwMNDQ3Y2dnh4MGDaNeuHZo3b46bN2/i0qVLfB9aS0tLrF69Gnp6esjMzMSuXbvg7OyM8PBwKCsri/1Yxa28c6euro709HSh+2RkZEBdXZ0vrezcZWdnQ1NTE4aGhvDx8YGJiQk+ffqEoKAgjB07FidPnoS+vj6739OnT+Hk5ISCggIoKCjA398frVq1Evtx1hS69sRHXV0dUlJSSE1N5UtPTU2Ftra20H0MDQ3RpUsXfPnyBcOHD4e6ujoCAgKgqqqKKVOmAADOnj2L+fPn499//8XLly/Ru3dvDB06FJKSkjV+TDUlJycHJSUlQj+HVbnu1NXV+T63ycnJuHnzJgYNGoQ///wTCQkJWLVqFYqLizFnzhx2v2fPnmHs2LHs53bbtm316nNbV4SFhcHV1RUBAQHo3Lkz/vzzTwwYMABxcXFs15owT58+RZMmTdjnGhoatVFdoTgNVFxcXKpdRkFBAd/gSYB/jIK4fR/RMgxTYZQrLP+36W5ubvD29sawYcPA4/Ggq6uLoUOH4uTJk+w+3/4qMzIygpWVFQYNGoRTp06JrTmxNoj73FlbW/P9wmrbti2GDx+OQ4cO8XWRGBgYICIiAh8+fMC5c+fg7u6OQ4cO1bsvPbr2xEfYek3l9YJLSEiAYRiMGzcOHz58AADMnz8fx44dw5w5c/Dlyxf8/vvv2LVrF+Lj48EwDF6+fIl9+/bxDbj9WVT3uistLYWamhq8vb0hKSkJc3NzpKenY8+ePXyBir6+PsLDw/Hx40ecO3cOHh4eCAoKqnef2zJcjVHZtGkTpkyZgqlTpwL42kBw9uxZ7Nixg+3iFUZTU7PO/BipN6PiyvsS8fHxgZKSEt+jopMvKhUVFUhKSgr8ksjKyhL4BVFGXV1daH4pKSkoKSkBAFRVVeHn54cbN27g77//xsmTJyEvLy90oFoZBQUFGBkZITExsZpHVTvKO3eZmZkCv9bKCGttKTt35X14JCQkYGlpiYSEBL50GRkZ6OnpwdLSEgsWLICpqWm96uela098MjIyUFxcLNB6oqmpKdDKUiYlJQVv375lgxQA+O+//yAhIQFdXV223OHDh0NRURF6enowNTXFp0+f8Pr165o7mBqmrKwslusuMzOT73OroaEBPT09vtYmQ0NDZGRkoLCwkE0r+9xaWFhg/vz5MDExYcdX1Ufi6vopKCjAhw8f+B7f/1gvU1hYiLt377IDmcvY29sjJiamwvra2NhAR0cHvXv3xsWLF8V2HkTBWaBiZmaG4OBgvgtTmOfPn2PWrFlYt26d0O0eHh7Izc3le3h4eIi9vtLS0jAzM8ONGzf40m/cuIE2bdoI3cfKykog//Xr19G6dWuBQWSysrLQ0tJCcXExzp8/j549e5Zbl8LCQrx69arcf/J1jYyMDMzNzQWmvsbExJQ7VsTa2lrgg3T16lVYWFiUOwCPYRj8999/P2yiZBjmh9ddXULXnvgUFRXh7t276Nu3L1963759y/3ivnbtGpo2bco3K8XY2BglJSUC43UKCgrw7t07SElJwdHRka91qr4p+9x+f16q+rm9du0azM3N2euubdu2SExM5Bv4m5CQAA0NDYExU9+rT5/bmlKVH+cZGRkoKSmBlpYWX3rZTCthdHR0EBgYiOPHj7NLQPTu3ZvT6eGcdf34+/vD3d0dc+bMgb29PWxtbdG0aVPIyckhOzsbcXFxuHr1KuLi4jB37ly+EfrfkpWVFdrVUxMzLyZMmABPT0+0bt0abdq0wfHjx5GSksLOHtmyZQvS0tKwevVqAF9nuoSGhsLX1xeOjo54+PAhTpw4wRd0PXr0CGlpaTA1NUVaWhp27NiB0tJSTJw4kc2zceNGdO/eHdra2sjKysKuXbuQl5dXr+6BNGnSJLi5ucHCwgI2NjYICwtDSkoKO+Nk48aNSE1Nxfr16wEATk5OOHz4MHx8fDB69Gjcv38fx48fx8aNG9kyt2/fjjZt2kBfX58doxIfHw8vLy82z6ZNm9CtWzdoa2sjLy8PUVFRuHXrFnbv3l27J6Ca6NoTn02bNuHgwYO4c+cOrl+/junTp6NFixbYuXMnAGDNmjVo1qwZ2zUdHByMZcuWYd++ffDy8oK6ujp8fX2xd+9e9numffv2aNasGR48eIBmzZphxYoVkJCQYK/n+srFxQWLFy+GhYUFrK2tceTIEaSkpGDMmDEAvp7L1NRU9rpycnJCcHAw1q5di1GjRuHBgwcIDw/Hhg0b2DKdnJxw6NAhrFmzBuPGjcObN28QGBiI8ePHs3k2b96Mrl27QkdHh+9z++00+fpGXF0/Hh4eAjOufjTcoSrdxmXrUZXp1KkTkpKSsGHDBnTr1k3EWlcPZ4FKr169cPv2bcTExCAsLAzBwcFISEjA58+foa6uDhsbGzg7O2P8+PF1pp+sf//+yM3NRWBgINLT09GqVSv4+/uzTeUZGRl8Uaquri78/f3h6+uLsLAwaGhowN3dnV3HAvj6C8Hf3x/JyclQUFBAly5dsHr1ar5BTKmpqVi8eDGys7OhoqICKysrHDx4sMIm+rrGwcEB2dnZCAgIQFpaGoyNjREYGIhmzZoBANLT05GSksLmb968OQIDA+Hj44PDhw9DU1MTnp6efGuofPjwAcuXL0d6ejoaN26M1q1b49ChQ7CysmLzZGRkwM3NDWlpaWjcuDFMTEywe/dugZkJdR1de+Jz5MgRqKmpYfny5dDR0cGTJ0/g4ODAdmfp6OjwDTLMy8tD3759sW3bNty5cweZmZk4cuQI3zgoOTk5/PHHHzA0NMSnT58QFRWFCRMmIDc3t9aPT5wcHByQk5ODgIAApKenw8jICDt37iz3c6urq4udO3di7dq1CA4OhqamJpYsWcLX9aCjo4Pdu3dj7dq1GDZsGLS0tDBhwgR2DAXw9Xp2d3dnP9tl3xf17XP7LXEFKuX9OBdGXV0dkpKSAq0naWlpAq0sFenYsSMOHTpUpXqKU51aR0WcaqJFpaGoiXVUGhK69kRTW+uo/KxqYh2VhqI21lEpmyFWXXv27KlS/g4dOqBdu3Z8a/y0bt0aQ4cOrfR4zpEjRyIrKwsXLlyo0muLC61MSwghhNQwroLw+fPnY8KECbC1tUWnTp0QGBiIxMREzJw5E8DXrqS3b9+yEwz8/Pygr68Pc3NzFBYW4tChQzh+/DiOHz/OSf0BjgOVrVu3Virfb7/9VsM1IYQQQmoOV4HKmDFjkJmZCW9vb6SkpMDCwgJRUVHQ09MD8HVW27ez+AoLC7Fw4UK8ffsW8vLyMDc3R2RkJBwcHDipP8Bx109llsrn8Xh49epVlcum5nfRUddP9dC1Jxrq+qke6voRXW10/UybNk0s5ezatUss5dQnnLao1Oc1BgghhJDKoiBcdDRGhRBCCKlhFKiIjtNApbKrg9bnpboJIYQQClREx2mgMnHiRDRq1AhSUlLlLpHP4/EoUCGEEEIaKE4DFTMzM6SmpmL8+PGYPHky30JdhBBCyM+CWlREx+lNCWNjYxEZGYnPnz+jW7dusLW1xY4dO/hu/kUIIYTUd+K6KWFDxPndkzt06IA///wTKSkp+O2333DkyBHo6Ohg3Lhx5d4RkhBCCCENA+eBShl5eXk4Oztj5cqVaN++PUJDQ5Gfn891tQghhJBqoxYV0dWJQOXt27dYs2YNjIyM4OTkhP/973+IjY2FiooK11UjhBBCqo0CFdFxOpj2yJEj2LdvHy5fvox+/fph48aNGDhwICQlJbmsFiGEEELqCE4DFScnJ7Ro0QLz5s2DlpYWEhIS4O/vL5CP7vVDCCGkPmuorSHiwGmg0qJFC/B4PAQHB5ebh8fjUaBCCCGkXqNARXScBioJCQlcvjwhhBBC6jhOA5UvX77gn3/+waBBgwAAHh4efFOSpaSk4O3tTXfzJYQQUq9Ri4roOA1UDhw4gNOnT7OByvbt22Fubg55eXkAQHx8PLS1tTF//nwuq0kIIYRUCwUqouN0evLhw4cxefJkvrTg4GBcvHgRFy9ehK+vL44ePcpR7QghhBDxoOnJouM0UHn27BmMjY3Z53JycpCQ+P8qtW/fHnFxcVxUjRBCCCF1AKddP7m5uZCS+v8qpKen820vLS2lZfQJIYTUew21NUQcOG1R0dXVxZMnT8rd/ujRI+jq6tZijQghhBDxo64f0XEaqDg4OGD58uX48uWLwLbPnz9j5cqVGDhwIAc1I4QQQkhdwGnXz5IlS3DkyBGYmJhg7ty5MDY2Bo/HQ3x8PLZv347i4mIsWbKEyyoSQggh1dZQW0PEgdNARUtLCzExMZg1axYWL14MhmEAfH1D+/bti4CAAGhpaXFZRUIIIaTavp0oQqqG00AFAAwMDHDmzBlkZWXhxYsXAIBWrVpBVVWV45oRQgghhGucByplVFVV0b59e66rQQghhIgddf2Irs4EKoQQQsjPigIV0VGgQgghhNQwClRER6N7CCGEEFJnUYsKIYQQUsOoRUV0FKgQQgghNYwCFdH9tIGKnJwc11UgDRRde6IrW0uJiIbW6iA/o582UOnRowfXVai3Ll26xHUV6rWcnByuq1AvKSsrC72dBqkcOTk5ZGZmcl2NeklNTa3GX4NaVET30wYqhBBCSF1BgYroqJ2QEEIIIWKzf/9+5Ofni628KgcqRUVF6NmzJ549eya2ShBCCCE/Mx6PJ5ZHfeDh4QFtbW1MmTIFMTEx1S6vyoGKtLQ0njx5Um9OGCGEEMK1hhSoJCcn49ChQ8jOzkbPnj1hamqKdevW4f379yKVJ1LXj7OzM/bs2SPSCxJCCCHk5yUpKYkhQ4YgPDwcSUlJmD59Og4fPowWLVpgyJAhOHnyJEpLSytdnkiDaQsLC7F7925ER0fD1tYWioqKfNs3bdokSrGEEELIT6m+tIaIm6amJjp37oynT5/i2bNnePz4MSZOnAhlZWXs27evUjN0RQpUnjx5grZt2wKAwFiVhvpmEEIIIeVpaGvcpKam4uDBg9i3bx9evXqFYcOG4fTp0+jTpw8+f/6MpUuXwsXFBW/evPlhWSIFKhcvXhRlN0IIIaRBakg/4gcPHoyzZ8/C2NgY06ZNg7OzM1RVVdnt8vLyWLBgATZv3lyp8qq1jsqLFy/w8uVLdOvWDfLy8mAYpkG9GYQQQgjhp6mpicuXL6NTp07l5tHR0cHr168rVZ5IgUpmZiZGjx6Nixcvgsfj4fnz5zA0NMTUqVOhrKyMjRs3Vqm8goIC3Lp1CwkJCcjPz4eGhgZsbGxgYGAgSvUIIYSQOqUh/Yjv3r07OzzkW4WFhQgNDYWzszN4PB709PQqVZ5InWbz5s2DtLQ0EhMToaCgwKaPGTMGZ86cqXQ5MTExGDt2LJSVldGjRw+4urpi1apVGD9+PFq1agUjIyP4+vri48ePolSTEEIIqRMa0vTkSZMmITc3VyD948ePmDRpUpXLEylQOXfuHNatWwddXV2+dCMjo0oNjAGAoUOHYuTIkWjWrBnOnj2Ljx8/IjMzE8nJycjPz8fz58+xdOlSnD9/HsbGxoiOjhalqoQQQgipReUNA0lOToaSklKVyxOp6ycvL4+vJaVMRkYGZGVlK1WGvb09jh49ChkZGaHbDQ0NYWhoCBcXF8TGxuLdu3eiVJUQQgjhXH1pDakOGxsbtuWnd+/ekJL6/xCjpKQEr1+/Rv/+/atcrkiBSrdu3RAUFIRVq1YB+PoGlJaWwtfXFz179qxUGXPmzKn065mbm8Pc3FyUqhJCCCGcawiByrBhwwAADx48QL9+/dCoUSN2m4yMDPT19eHo6FjlckUKVHx9fdGjRw/cuXMHhYWFcHNzQ2xsLLKysnDt2jVRiiSEEEJIPebl5QUA0NfXx5gxYyAnJyeWckUao9K6dWs8evQI7du3R9++fZGXl4cRI0bg/v37aNmyZaXKUFVVRUZGBgBARUUFqqqq5T4IIYSQ+ozLwbQBAQEwMDCAnJwc2rVrhytXrlRqv2vXrkFKSgrW1tZVej0XFxexBSmAiC0q58+fR+/evbFy5UqBbdu3b8fcuXN/WMbmzZvRuHFj9u+G0CxGCCGkYeLqf1xYWBhcXV0REBCAzp07488//8SAAQMQFxeHFi1alLtfbm4unJ2d0bt3b6Smpv7wdVRVVfHs2TOoq6tDRUWlwuPNysqq0jHwGIZhqrQHAGVlZURHR+N///sfX7qfnx+WL1+ODx8+VLVIsavM/QOIcJcuXeK6CvVaTk4O11Wol5SVlfHlyxeuq1FvycnJITMzk+tq1Etqamo1/hrr1q0TSznu7u5Vyt+hQwe0bdsWO3bsYNPMzMwwbNgw+Pj4lLufk5MTjIyMICkpiYiICDx48KDC1zlw4ACcnJwgKyuL/fv3VxiouLi4VOkYRGpR2bx5MxwcHHD58mW0bt0aALBhwwasWrUKkZGRVS7v3r17kJaWhqWlJQDg5MmT2LdvH1q3bo0VK1aUOzOIEEIIqQ+4aFEpLCzE3bt3sXjxYr50e3t7xMTElLvfvn378PLlSxw6dAh//PFHpV7r2+Bj4sSJItW3PCIFKpMmTUJmZibs7e1x9epVhIWFYc2aNfj7779hZ2dX5fJmzJiBxYsXw9LSEq9evcKYMWMwYsQIHD16FPn5+fDz8xOlmoQQQkidIK5ApaCgAAUFBXxpsrKyQpcGycjIQElJCbS0tPjStbS08P79e6HlP3/+HIsXL8aVK1f4phf/SFV6Upo0aVLpvEA17vWzcOFCZGZmwtbWFiUlJTh37hw6dOggUlnPnj1jB+scPXoU3bt3R3BwMK5duwYnJycKVAghhNRr4rp7so+Pj8D4UC8vL6xYsaLcfb4PkspbkK2kpAS//PILVq5cCWNj4yrVS1lZ+YfBWNnrlpSUVKnsSgcqW7duFUjT0dGBgoICunXrhps3b+LmzZsAgN9++61KlWAYBqWlpQCAf/75B4MGDQIANG/enJ0ZRAghhDR0Hh4emD9/Pl9aeQutqqurQ1JSUqD1JC0tTaCVBfi6xP2dO3dw//59dlJMaWkpGIaBlJQUzp07h169egl9rYsXL4pyOJVS6UClvNsxS0pK4tq1a+z6KTwer8qBiq2tLf744w/06dMHly9fZgf9vH79WujJJIQQQuoTcXX9lNfNI4yMjAzatWuH6OhoDB8+nE2Pjo7G0KFDBfI3adIEjx8/5ksLCAjAhQsXcOzYsQpvFNy9e/dKHkHVVTpQqeztmEXh5+eHcePGISIiAp6enmjVqhUA4NixYyKNeSGEEELqEq6mJ8+fPx8TJkyAra0tOnXqhMDAQCQmJmLmzJkAvrbQvH37FkFBQZCQkICFhQXf/pqampCTkxNI/96jR49gYWEBCQkJPHr0qMK8VlZWVTqGKo9RKSoqgomJCU6fPs3O+KkuKysrgSgO+LoCrqSkpFhegxBCCGloxowZg8zMTHh7eyMlJQUWFhaIioqCnp4eACAlJQWJiYnVfh1ra2u8f/8empqasLa2Bo/Hg7DVT2p0jEoZaWlpFBQU1Ep0KM6V7QghhBCucLmo6ezZszF79myh2/bv31/hvitWrKhwoG6Z169fQ0NDg/1bnESa9fPrr79i3bp12L17d5WmL5VHQkKiwjexqtEXIYQQUpf87Kuvl7XQfP+3OIgUZdy8eRPnz5/HuXPnYGlpCUVFRb7t4eHhVSrvxIkTfM+Liopw//59HDhwQOgy/YQQQgipu54+fYpt27bhv//+A4/Hg6mpKX799VeYmJhUuSyRAhVlZWWRbtVcHmGjj0eOHAlzc3OEhYVhypQpYnstQgghpLb97C0q3zp27BjGjh3LDuAFgBs3bsDCwgLBwcEYNWpUlcoTKVDZt2+fKLtVWYcOHTBt2rRaeS1RDR06FE5OTlBTU8Pr16+xfft2oQODga+DjYQtXufs7MwOZvLz8xN6p8rr16/Dw8NDnFWvdYcPH8aePXuQnp4OIyMjLFmyBLa2tuXmv3XrFtauXYvnz59DU1MTU6dOxdixY9nt4eHhQs/Jo0eP2Ol7wcHBCAkJwdu3bwEARkZGmD17do1Opastx44dw6FDh5CZmQkDAwPMmzcPNjY2QvNmZGRgy5YtiI+PR1JSEkaPHi2wFsO3zp07h2XLlqFbt27w9fWtqUOoNWFhYdi/fz8yMjLQsmVLuLm5oW3btuXmv3PnDjZs2ICXL19CQ0MDEydOxOjRo9ntRUVF2LNnD/766y+kpaVBX18frq6u6Ny5M5tnx44d2LlzJ1+5ampquHDhgvgPsJYdP34cwcHB7LX3+++/l3uH3YyMDGzbtg1Pnz5FUlISRo0aBVdXV748kZGRWL16tcC+Fy9erPRU3LquIQUqbm5u8PDwgLe3N1+6l5cX3N3daydQqQ2fP3/Gtm3boKury3VVytWzZ0/MnTsXfn5+ePz4MYYMGYL169fDxcUFaWlp5e43fvx45Ofns8+/vYndsmXLIC0tzT5v0qQJ9uzZg8uXL9fIMdSWqKgo+Pj4wMvLC23btkVoaCimTZuGyMhING3aVCB/UlISpk+fjlGjRsHX1xf37t3DypUroaqqin79+rH5GjVqhDNnzvDt++0Xm7a2NhYuXMjeJTQiIgJz5szBiRMnYGRkVENHW/Oio6OxefNmuLm5wcrKCidOnMC8efMQGhoKbW1tgfyFhYVQVlbGpEmTEBISUmHZKSkp2Lp1a5Vv7V5XnTlzBuvXr4enpyesra1x7NgxzJ49GydOnICOjo5A/uTkZMyZMweOjo5Ys2YNHjx4gNWrV0NVVRV9+vQB8PUu8ZGRkfDy8oKBgQFiYmIwb948HDhwAGZmZmxZLVu2RGBgIPtcXKuTcumff/7Bli1bsHDhQlhZWSEiIgILFizA4cOHhV57RUVFUFZWhouLC0JDQ8stV1FRUWD7zxKkNDTv37+Hs7OzQPr48eNF+uEjcqBy7NgxHDlyBImJiSgsLOTbdu/evSqV9f0toRmGwcePH6GgoIBDhw6JWsUaN2rUKERFRbE3Yty+fTv+97//YejQodi1a1e5++Xk5ODTp09Ct338+JHvea9evfDly5d6f0fjffv2wdHRkY2kPT09cfXqVYSEhGDBggUC+UNDQ6GjowNPT08AX7/wHz9+jL179/IFKjwejx1pLsz3qyjOmzcPISEhePDgQb0OVEJCQjBkyBC223T+/Pm4efMmjh8/jjlz5gjkb9q0KXue//rrr3LLLSkpwfLlyzF9+nQ8ePBA4Hqsjw4ePIjhw4djxIgRAL7+2ouJicGRI0fw+++/C+Q/evQodHR04ObmBgAwNDREbGwsDhw4wAYqkZGRmDp1Krp27QoAGD16NGJiYhAUFMR3R1opKSmoq6vX9CHWqtDQUAwePBhDhgwBALi6uuLmzZs4ceIEZs2aJZBfR0cH8+bNAwCcPn263HJ5PF6t3MWYKw2pRaVHjx64cuUKuyZamatXr7KfmaoQKVDZunUrPD094eLigpMnT2LSpEl4+fIlbt++LfRL8ke+7w6RkJCAhoYGOnToABUVFVGqWOOkpKRgYmKC4OBgvvTbt2/D3Ny8wn137doFGRkZJCQk4ODBgxXePtvBwQEXLlzAly9fxFFtThQWFiI2NhbTp0/nS+/cuTPu378vdJ8HDx7wNaMDQNeuXXH8+HEUFRWxrU75+fno2bMnSkpKYGZmht9//73c9X1KSkpw5swZ5Ofnl9tFUh8UFRUhPj5e4BdL+/bty+12rKw9e/ZARUUFQ4YM+eFt3euDoqIi/Pfff5g8eTJfeqdOnfDw4UOh+zx69IjtVy9jZ2eHiIgI9torLCwUuKu7rKyswDl78+YN+vTpw94d/rfffqvTrcQ/UlRUhKdPn2LChAl86eK49j5//ozhw4ejtLQURkZGmDZtmkgDL+uqnz1QOXXqFPv3kCFD4O7ujrt376Jjx44Avo5ROXr0qEgTZEQKVAICAhAYGIixY8fiwIEDcHNzg6GhIZYvX46srKxKl7N3716MGzeO7/bQ9YWSkhIkJSWRnZ3Nl56dnQ1VVVWh+2RmZsLX1xfPnj2DtLQ07O3tsWnTJri6ugpdyc/U1BSGhoZYv359jRxDbcnOzkZJSYnAryV1dXWkp6cL3ScjI0Pgl6iamhqKi4uRnZ0NTU1NGBoawsfHByYmJvj06ROCgoIwduxYnDx5Evr6+ux+T58+hZOTEwoKCqCgoAB/f3+BSL8+ycnJQUlJicB1pqamhhs3bohc7sOHD3Hq1Kk63YpZVeVde2pqauXeRywjI0No/uLiYuTk5EBDQwN2dnY4ePAg2rVrh+bNm+PmzZu4dOkS31IKlpaWWL16NfT09JCZmYldu3bB2dkZ4eHhUFZWFvux1obyrj1VVdUqffd/T09PD56enmjZsiXy8vJw5MgRzJw5E0FBQWjevHl1q10n/OyByrBhwwTSAgICEBAQwJc2Z84cdlXcyhIpUElMTGSXtpeXl2ebhydMmICOHTti+/btlSpn2rRpGDRoEDQ1NQF8bZ6OiYnh+yfzI+Xd8rq2fL/yXkUXY1JSEpKSktjncXFx0NTUxJgxY4QGKg4ODnj16hXi4+PFV2EOVfYOnhXl/zbd2tqabxxF27ZtMXz4cBw6dAhLly5l0w0MDBAREYEPHz7g3LlzcHd3x6FDh+p1sAJU/XxWJC8vD15eXliyZEm9/SdaEXFfe25ubvD29sawYcPA4/Ggq6uLoUOH4uTJk+w+Xbp0Yf82MjKClZUVBg0ahFOnTgntv6/PhK1AWhUWFhZ8S7RbWVlh0qRJOHr0aIWDvkndUXZj4ZogUqCira2NzMxM6OnpQU9PDzdu3ECbNm3w+vXrKl2w3+f9+PFjlQ+2vFte17Tc3FyhvyyUlZWr9MsiLi4Offv2FUiXlZVFr169am2GVU1SUVGBpKSkwC/YzMzMcvvvhbW2ZGVlQUpKqtx/pBISErC0tERCQgJfuoyMDLsAkaWlJR4/foygoCCBEen1hbKyMiQlJZGZmcmXnpWVVW5r3o+8ffsWKSkpWLhwIZtW9lm0s7PDkSNH6mWXRXnXXlZWVrnjIdTV1YXml5KSgpKSEoCvLQh+fn4oKChATk4ONDU14efnJ3RgeBkFBQUYGRmJZblyrpRde99/x1XUkiwKCQkJmJqaIjk5WWxlcu1nGEjNFZEClV69euGvv/5C27ZtMWXKFMybNw/Hjh3DnTt32AFrtaW8W17X9ODT4uJiPH36FLa2trh69Sqbbmtry95JujKMjIwE/uEAX2cUycjIIDo6Wiz15ZKMjAzMzc1x7do1vqAsJiYGvXv3FrqPtbW1wG3Dr169CgsLC75ZUd9iGAb//fcfjI2NK6wPwzACA8DrE2lpaZiamuLWrVvo0aMHm37r1i1069ZNpDL19PQExlvt3LkT+fn5mD9/fr29i7m0tDTMzMxw48YNvmvtxo0bfOfuW1ZWVvj333/50q5fv47WrVsLXHuysrLQ0tJCUVERzp8/D3t7+3LrUlhYiFevXtXr8VHS0tIwMTHBrVu3+Kb43759W6RBkuVhGAbPnz9Hy5YtxVYm1372rp/v5eXl4fLly0In3Pz2229VKkukQCUwMJD9tTVz5kyoqqri6tWrGDx4cJX6nng8Ht+b9/3zyqjKLa/F7ejRo1iyZAmePn2K2NhYDB48GFpaWuygomnTpkFdXZ2dBTBy5Ei8f/8er1+/hrS0NPr27Yvu3btj2bJlAmU7ODjg6tWr+PDhQ60eU02ZNGkS3NzcYGFhARsbG4SFhSElJQVOTk4AgI0bNyI1NZUdj+Pk5ITDhw/Dx8cHo0ePxv3793H8+HFs3LiRLXP79u1o06YN9PX12TEq8fHxfC1qmzZtQrdu3aCtrY28vDxERUXh1q1b2L17d+2eADEbO3YsVqxYAVNTU1haWiIiIgKpqansDwV/f3+kp6fz3aPj2bNnAL4OQM7JycGzZ88gJSUFQ0NDyMrKCvxTaNy4MQDU+38WEyZMgKenJ1q3bo02bdrg+PHjSElJYWegbdmyBWlpaew6HqNGjUJoaCh8fX3h6OiIhw8f4sSJE1i3bh1b5qNHj5CWlgZTU1OkpaVhx44dKC0txcSJE9k8GzduRPfu3aGtrY2srCzs2rULeXl57GyZ+srJyQne3t4wMzODhYUFTp48idTUVHaMwo4dO5Ceno7ly5ez+5Rde58/f2avPWlpaRgYGAD4Oojb3NwczZs3R15eHo4ePYrnz5/ztfCR+uP+/ftwcHBAfn4+8vLyoKqqioyMDCgoKEBTU7N2AhUJCQm+ZqzRo0fzLYZUWQzDwNjYmA1OPn36BBsbG4EmsuoM0qpJFy9eRJMmTeDi4gJVVVW8fv0a7u7uSE1NBfB1AN63v0SlpKQwa9YsqKuro6CgAAkJCXB3d8fNmzf5ytXV1YWVlZXQabv1lYODA7KzsxEQEIC0tDQYGxsjMDAQzZo1AwCkp6cjJSWFzd+8eXMEBgbCx8cHhw8fhqamJjw9PfmmJn/48AHLly9Heno6GjdujNatW+PQoUN8txDPyMiAm5sb0tLS0LhxY5iYmGD37t0CM4rqm759+yI3Nxd79+5FRkYGDA0NsXnzZnZdkMzMTPY6LPPtTI34+HicPXsWOjo6iIiIqM2q17r+/fsjNzcXgYGBSE9PR6tWreDv789202RkZOD9+/dsfl1dXfj7+8PX1xdhYWHQ0NCAu7s7OzUZ+No64u/vj+TkZCgoKKBLly5YvXo1mjRpwuZJTU3F4sWLkZ2dDRUVFVhZWeHgwYMVdg/VB3369GGvvczMTBgaGmLDhg0VXnvfBnDx8fE4d+4ctLW12dutfPr0CevWrUNWVhYUFRVhbGyMgICAcmfw1UcNqUVl3rx5GDx4MHbs2AFlZWXcuHED0tLSGD9+vNAlAX6Ex4g4CurLly/sr4rvx5VU9hfDgQMHKpVPlFlB5TXrkh+r72u2cO3bBfxI5SkrK9frafhck5OTE9qNTH6sNtZvqWhtraqo66u1A18/yzdv3oSJiQmUlZVx/fp1mJmZ4ebNm3BxcanyBBGRWlTOnDkDZ2dnodP7eDxepe92XB+nJRNCCCGkfNLS0mwLkpaWFhITE2FmZgYlJSWRBpOLNAx57ty5GDVqFFJSUlBaWsr3qGyQUlXVnf5GCCGEcKVsyER1H/WBjY0N7ty5A+DrxJDly5fj8OHDcHV1haWlZZXLE+mo09LSqj0TwMzMDMHBwT+cffH8+XPMmjWLbyAbIYQQUp+UTRap7qM+WLNmDTtmadWqVVBTU8OsWbOQlpbGd++ryhKp62fkyJG4dOlStWYD+Pv7w93dHXPmzIG9vT1sbW3RtGlTyMnJITs7G3Fxcbh69Sri4uIwd+5czJ49W+TXIoQQQkjtsLW1Zf/W0NBAVFRUtcoTKVDZvn07Ro0ahStXrsDS0lJgbYHKTD3q1asXbt++jZiYGISFhSE4OBgJCQn4/Pkz1NXVYWNjA2dnZ4wfP/6nXCmTEEJIw1FfWkPEKS0tDU+fPgWPx4OJiUmFN5CtiEiBSnBwMM6ePQt5eXlcunRJYC2UqsyRtrOzY5fjJ4QQQn5GDSlQ+fDhA+bMmYPQ0FB23KqkpCTGjBkDf39/doXnyhJpjMrSpUvh7e2N3NxcJCQk4PXr1+zj1atXohRJCCGE/LQa0hiVqVOn4ubNmzh9+jRycnKQm5uL06dP486dOyJNrxapRaWwsBBjxoyp9gjkrVu3VipfVVexI4QQQgg3IiMjcfbsWb4bc/br1w+7du1C//79q1yeSIGKi4sLwsLCsGTJElF2Z23evPmHearalUQIIYTUNfVlarE4qKmpCe3eUVJSgoqKSpXLEylQKSkpwfr163H27FlYWVkJDKbdtGlTpcp5/fq1KC9PCCGE1Cv1pdtGHJYuXYr58+cjKCiInab8/v17LFq0SOi97X5EpEDl8ePH7B1Anzx5wretIb0ZhBBCCPm6yNu3//+fP38OPT09tGjRAgCQmJgIWVlZpKenY8aMGVUqW6RA5eLFi6LsJuDChQuYO3cubty4wXczLwDIzc2FnZ0dduzYIfKt6wkhhJC64Gf/EV929+yaIFKg8q3k5GTweDz2LrhV4efnh2nTpgkEKcDXvqwZM2Zg8+bNFKgQQgip1372QMXLy6vGyhZpdE9paSm8vb2hpKTENu0oKytj1apVAndSrsjDhw8rHAFsb2+Pu3fvilJFQgghhHDo7t27OHToEA4fPoz79++LXI5ILSqenp7Ys2cP1q5di86dO4NhGFy7dg0rVqzAly9fsHr16kqVk5qaKjAQl69yUlJIT08XpYqEEEJIndGQZv2kpaXByckJly5dgrKyMhiGQW5uLnr27InQ0NAqr1Ar0pk7cOAAdu/ejVmzZsHKygpt2rTB7NmzsWvXLuzfv7/S5TRr1gyPHz8ud/ujR4/YEcOEEEJIfdWQFnz79ddf8eHDB8TGxiIrKwvZ2dl48uQJPnz4INJyIyIFKllZWTA1NRVINzU1RVZWVqXLcXBwwPLly/HlyxeBbZ8/f4aXlxcGDRokShUJIYQQwoEzZ85gx44dMDMzY9Nat24Nf39//P3331UuT6RApU2bNti+fbtA+vbt29GmTZtKl7N06VJkZWXB2NgY69evx8mTJ3Hq1CmsW7cOJiYmyMrKgqenpyhVJIQQQuqMhtSiUlpaKnRYh7S0dJXGsZYRaYzK+vXrMXDgQPzzzz/o1KkTeDweYmJikJSUVKXbOWtpaSEmJgazZs2Ch4cHGIYB8PUN7devHwICAqClpSVKFQkhhJA6o74EGeLQq1cv/P777wgJCUHTpk0BAG/fvsW8efPQu3fvKpcnUqDSvXt3PHv2DP7+/oiPjwfDMBgxYgRmz57NVqoyXr16BQMDA0RFRSE7OxsvXrwAwzAwMjISaZldQgghpC5qSIHK9u3bMXToUOjr66N58+bg8XhITEyEpaUlDh06VOXyRF5HpWnTppWe3VMeIyMjpKSkQFNTEyoqKtiwYQO2bt1KQQohhBBSTzVv3hz37t1DdHQ025jRunVr9OnTR6TyKh2oPHr0CBYWFpCQkMCjR48qzGtlZVWpMsu6espERUXBx8enslUihBBC6oWGMj25uLgYcnJyePDgAfr27Yu+fftWu8xKByrW1tZ4//49NDU1YW1tDR6PJxBoAF+bt0pKSqpdMUIIIeRn0VC6fqSkpKCnpyfWOKDSgcrr16/ZRVrEdddjYaOYG8qbSQghhPyMli5dCg8PDxw6dAiqqqrVLq/SgYqenh77d6NGjaCmpgYASEpKwq5du/D582cMGTIEXbt2rfSLMwyDiRMnQlZWFgDw5csXzJw5E4qKinz5wsPDK10mIYQQUtc0pB/hW7duxYsXL9C0aVPo6ekJ/E+/d+9elcqr0mDax48fY/DgwUhKSoKRkRFCQ0PRv39/5OXlQUJCAps3b8axY8cqfRdFFxcXvufjx4+vSnUIIYSQeqEhBSrDhg0rd3iIKKoUqLi5ubHTiw4dOoRBgwbBwcEBu3fvBvB12dy1a9dWOlDZt29flStMCCGEkLonPz8fixYtQkREBIqKitC7d29s27YN6urq1Sq3SoHK7du3ceHCBVhZWcHa2hqBgYGYPXs2O5r5119/RceOHatVIUIIIeRn0xBm/Xh5eWH//v0YN24c5OXlERwcjFmzZuHo0aPVKrdKgUpWVha0tbUBfB2noqioyDdQRkVFBR8/fqxWhQghhJCfTUPo+gkPD8eePXvg5OQEABg3bhw6d+6MkpISSEpKilxulUM8mqVDCCGEkO8lJSXxTahp3749pKSk8O7du2qVW+WVaSuapVNQUFCtyhBCCCE/o4bwo76kpAQyMjJ8aVJSUiguLq5WuVUKVCozS8fZ2blaFSKEEEJ+Ng0hUPl+yRFA+LIjVV1yhMeIa/4QIYQQQoT6+++/xVLOgAEDxFJOTZg0aVKl8lV1xq/INyUkP7eWLVtyXYV66eXLl8jJyeG6GvWSsrIy11Wo9+jaE83Pfu0FBATA19cXKSkpMDc3h5+fX7mLs169ehXu7u6Ij49Hfn4+9PT0MGPGDMybN++Hr1NTS45QoEIIIYTUMK6mJ4eFhcHV1RUBAQHo3Lkz/vzzTwwYMABxcXFo0aKFQH5FRUXMnTsXVlZWUFRUxNWrVzFjxgwoKipi+vTpHBwBdf2QclCLimioRUV0P/uv2tpA155oauPaO3funFjKsbe3r1L+Dh06oG3bttixYwebZmZmhmHDhsHHx6dSZYwYMQKKioo4ePBglV5bXH7+FWgIIYSQn0RBQQE+fPjA9yhvxm1hYSHu3r0rENzY29sjJiamUq93//59xMTEoHv37tWuu6goUCGEEEJqGI/HE8vDx8cHSkpKfI/yWkYyMjJQUlICLS0tvnQtLS28f/++wvrq6upCVlYWtra2mDNnDqZOnSq2c1FVNEaFEEIIqWHimp7s4eGB+fPn86V9Ox24Mq/NMMwP63PlyhV8+vQJN27cwOLFi9GqVSuMHTtWtEpXEwUqhBBCSD0hKyv7w8CkjLq6OiQlJQVaT9LS0gRaWb5nYGAAALC0tERqaipWrFjBWaBCXT+EEEJIDZOQkBDLoypkZGTQrl07REdH86VHR0fDzs6u0uUwDMPpyvPUokIIIYTUMK5Wpp0/fz4mTJgAW1tbdOrUCYGBgUhMTMTMmTMBfO1Kevv2LYKCggAA/v7+aNGiBUxNTQF8XVdlw4YN+PXXXzmpP1BHApWCggLcunULCQkJyM/Ph4aGBmxsbNimJ0IIIYRU3ZgxY5CZmQlvb2+kpKTAwsICUVFR0NPTAwCkpKQgMTGRzV9aWgoPDw+8fv0aUlJSaNmyJdauXYsZM2ZwdQjcrqMSExODbdu2ISIiAoWFhVBWVoa8vDyysrJQUFAAQ0NDTJ8+HTNnzkTjxo25qmaDROuoiIbWUREdraNSfXTtiaY2rr2LFy+KpZyePXuKpZz6hLMxKkOHDsXIkSPRrFkznD17Fh8/fkRmZiaSk5ORn5+P58+fY+nSpTh//jyMjY0F+tgIIYSQ+kJc05MbIs66fuzt7XH06FGBW0KXMTQ0hKGhIVxcXBAbG4t3797Vcg0JIYQQ8WioQYY4cBaozJkzp9J5zc3NYW5uXoO1IYQQQkhdxOn05OzsbGzbtg0fPnwQ2Jabm1vuNkIIIaQ+oa4f0XEaqGzfvh3//vsvmjRpIrBNSUkJV65cwbZt2zioGSGEECI+FKiIjtNA5fjx4+xcbmFmzJiBY8eO1WKNCCGEEFKXcLqOysuXL2FkZFTudiMjI7x8+bIWa0QIIYSIX0NtDREHTltUJCUlK5zN8+7duyovGUwIIYTUNdT1IzpOowAbGxtERESUu/3EiROwsbGpvQoRQgghpE7htOtn7ty5cHJygq6uLmbNmgVJSUkAQElJCQICArB582YEBwdzWUVCCCGk2hpqa4g4cBqoODo6ws3NDb/99hs8PT1haGgIHo+Hly9f4tOnT1i0aBFGjhzJZRUJIYQQwiHOb0q4evVqDB06FIcPH8aLFy/AMAy6deuGX375Be3bt+e6eoQQQgjhEOeBCgC0b9+eghJCCCE/Ler6ER2ng2mfP3+OsWPHlrsy7S+//IJXr15xUDNCCCFEfGjWj+g4DVR8fX3RvHnzclembd68OXx9fTmoGSGEECI+FKiIjtNA5d9//8WoUaPK3T569GhcuHChFmtECCGEkLqE0zEqb968gaamZrnb1dXVkZSUVIs1IoQQQsSvobaGiAOnLSpKSkoVLpH/4sULod1ChBBCSH1CXT+i4zRQ6datW4V3R966dSu6du1aizUihBBCSF3CaaDi4eGBv//+GyNHjsStW7eQm5uL3Nxc3Lx5E46Ojjh79iw8PDy4rKKAw4cPo1evXrC0tMSIESNw586dCvPfunULI0aMgKWlJXr37o2QkBC+7eHh4TAxMRF4FBQUsHmCg4MxePBgtG3bFm3btsWYMWNw+fLlGjm+2jZu3DhcunQJcXFxOHnyJGxtbSvMLyMjgwULFuDff/9FXFwcLly4wLco4JgxYxAaGop79+7h3r17CAoKgpWVVU0fBieOHTuGYcOGoWvXrnB2dsb9+/fLzZuRkYFly5Zh1KhR6NixIzZt2lRh2efOnUOHDh2waNEicVebE/S5FS+69qqOWlREx/m9fo4dO4Z///0XnTp1gqqqKlRVVWFnZ4crV67gyJEjaNu2LZdV5BMVFQUfHx/MmjULERERaNeuHaZNm1bujRWTkpIwffp0tGvXDhEREZg5cyZWr16Ns2fP8uVr1KgRrl69yveQlZVlt2tra2PhwoU4fvw4jh8/jo4dO2LOnDl4/vx5jR5vTRs4cCCWLl2KgIAADB48GLdv38bevXuho6NT7j5bt25Fp06dsHjxYvTt2xeurq58U9g7dOiAv/76C+PGjcPIkSPx7t07HDhwAFpaWrVxSLUmOjoamzdvxqRJkxAUFARra2vMmzcP79+/F5q/sLAQysrKmDRpUoV3LAeAlJQUbN26FdbW1jVQ89pHn1vxomtPNBSoiI7HMAzDdSU+f/6MM2fOsCvTGhsbw97eHgoKClxXjc+oUaPQunVrrFy5kk0bMGAA+vTpgwULFgjk9/X1xYULF/D333+zacuXL8fTp08RFhYG4OsvszVr1vzwF9732rdvj0WLFlU4a6o6WrZsWSPlfuv48eOIjY3F8uXL2bSzZ88iOjoaGzZsEMjfrVs3bNmyBT169EBubm6lXkNCQgL37t3DypUrceLECbHVvTwvX75ETk5Ojb/O5MmTYWJiAnd3dzZtzJgx6NatG+bMmVPhvrNmzYKRkRHmz58vsK2kpAQzZ87E4MGD8eDBA3z8+LHWlghQVlaukXIb0ueWrj3R1NS19627d++KpZx27dqJpZz6hNMWFQBgGAbJyckwMTHBvHnz4ObmhmHDhtW5IKWwsBCxsbHo0qULX3rnzp3LbfZ88OABOnfuzJfWtWtXPHnyBEVFRWxafn4+evbsiW7dumHGjBmIi4srtx4lJSWIjIxEfn5+vb6ztLS0NCwsLHD16lW+9KtXr5bbita7d288fvwY06dPx7Vr1/DPP//Aw8OD71fs9+Tl5SEtLV0rX+C1paioCPHx8ejQoQNfevv27fH48eNqlb1nzx6oqKhgyJAh1SqnrqDPrXjRtUe4wOn05ISEBAwdOhRPnjwBADRv3hzh4eF1qrunTHZ2NkpKSqCmpsaXrq6ujvT0dKH7ZGRkQF1dnS9NTU0NxcXFyM7OhqamJgwNDeHj4wMTExN8+vQJQUFBGDt2LE6ePAl9fX12v6dPn8LJyQkFBQVQUFCAv78/WrVqJfbjrC0qKiqQkpJCRkYGX3pGRgY0NDSE7tOiRQvY2tqioKAAs2bNgqqqKlauXAklJSUsXrxY6D6LFi1Camoqrl27JvZj4EpOTg5KSkqgqqrKl66mpoYbN26IXO7Dhw9x6tQpHDp0qLpVrDPocytedO2JrqF224gDpy0q7u7u+PLlCw4ePIijR49CR0cHM2fOrFIZBQUF+PDhA9/j2wFt4vb9xcYwTIUXoLD836ZbW1tj6NChMDU1ha2tLfz8/KCvry/wgTUwMEBERATCwsIwduxYuLu748WLF+I4JE593/PI4/EE0r7fNm/ePDx69AiXLl3C6tWr4ejoKLRVZfr06Rg8eDBmzZqFwsLCGqk/l6p6LVYkLy8PXl5eWLJkSa00g9c2+tyKF117pDZx2qJy5coVhISEoHv37gC+Nh/q6enh8+fPkJeXr1QZPj4+fH3PAODl5YUVK1aIta4qKiqQlJQUaAHIzMwU+PVVRtivtqysLEhJSZX7gZSQkIClpSUSEhL40mVkZKCnpwcAsLS0xOPHjxEUFARvb2/RDohj2dnZKC4uFmg9UVNTEzjHZdLT05GamopPnz6xaS9fvoSEhAR0dHT4ztnUqVMxa9YsODs74+nTpzVyDFxRVlaGpKQkMjMz+dKzsrIEfulW1tu3b5GSkoKFCxeyaaWlpQAAOzs7HDlyBLq6uqJXmiP0uRUvuvZERy0qouM0UHn//j1MTU3Z57q6upCXl0dqaipf82lFPDw8BAZmVTRmQVQyMjIwNzfHtWvX0LdvXzY9JiYGvXv3FrqPtbU1Ll68yJd29epVWFhYQFpaWug+DMPgv//+g7GxcYX1YRimXrcSFBUV4cmTJ+jcuTPOnTvHpnfu3Bn//POP0H3u3r2LAQMGQEFBAfn5+QC+/mItKSlBSkoKm2/atGmYM2cOJk6cWO1+87pIWloapqamuHXrFnr06MGm37p1C926dROpTD09PQQHB/Ol7dy5E/n5+Zg/f369nTVFn1vxomtPdBSoiI7TQIXH40FCgr/3SUJCotymf2FkZWVrJDARZtKkSXBzc4OFhQVsbGwQFhaGlJQUODk5AQA2btyI1NRUrF+/HgDg5OSEw4cPw8fHB6NHj8b9+/dx/PhxbNy4kS1z+/btaNOmDfT19dm+7vj4eHh5ebF5Nm3ahG7dukFbWxt5eXmIiorCrVu3sHv37lo57pqyd+9ebNiwAY8fP8b9+/fh5OSEpk2bsl9aCxcuZKd4AsCpU6cwd+5crFu3Dlu2bIGKigoWL16MY8eOsd1906dPh6urK+bNm4fk5GT2V3N+fj4b3PwMxo4dixUrVsDU1BSWlpaIiIhAamoqRowYAQDw9/dHeno6X8vis2fPAHw9Fzk5OXj27BmkpKRgaGgIWVlZgZlejRs3BlA7M8BqEn1uxYuuPdFQoCI6TgOVsqnI376Bnz59go2NDV8Ak5WVxUX1BDg4OCA7OxsBAQFIS0uDsbExAgMD0axZMwBfuya+/WXfvHlzBAYGwsfHB4cPH4ampiY8PT3Rr18/Ns+HDx+wfPlypKeno3HjxmjdujUOHTrEt0hZRkYG3NzckJaWhsaNG8PExAS7d+8WmJlQ30RGRkJZWRm//vorNDQ08Pz5c0yZMoVd30JTU5NvTZX8/Hw4OzvDy8sLERERyMnJQWRkJN8CUuPGjYOsrCwCAgL4XmvLli3YunVr7RxYLejbty9yc3Oxd+9eZGRkwNDQEJs3b2bPV2ZmJlJTU/n2mTBhAvt3fHw8zp49Cx0dHURERNRm1WsdfW7Fi649Uts4XUflwIEDlcrn4uJSwzUh3/uZfsnUptpaR+VnRAMpq4+uPdHUxrX36NEjsZTzs660XRFOW1QoACGEENIQUNeP6Dhf8K2y6sACuoQQQgipZZwFKmZmZggODv7hCPjnz59j1qxZWLduXS3VjBBCCBEvuteP6Djr+vH394e7uzvmzJkDe3t72NraomnTppCTk0N2djbi4uJw9epVxMXFYe7cuZg9ezZXVSWEEEKqpaEGGeLAWaDSq1cv3L59GzExMQgLC0NwcDASEhLw+fNnqKurw8bGBs7Ozhg/fjwNsiOEEEIaKE4H0wJfVx60s7PjuhqEEEIIqYM4D1QIIYSQnx11/YiO00Clsgtw/fbbbzVcE0IIIYTURZwGKps3b/5hHh6PR4EKIYSQeo1aVETHaaDy+vVrLl+eEEIIqRUUqIiOxqgQQgghNYwCFdFxujLthQsX0Lp1a3z48EFgW25uLszNzfHvv/9yUDNCCCGE1AWcBip+fn6YNm0amjRpIrBNSUkJM2bMqNQ4FkIIIaQuo5VpRcdpoPLw4UP079+/3O329va4e/duLdaIEEIIET8KVETHaaCSmpoKaWnpcrdLSUkhPT29FmtECCGEkLqE00ClWbNmePz4cbnbHz16BB0dnVqsESGEEELqEk4DFQcHByxfvhxfvnwR2Pb582d4eXlh0KBBHNSMEEIIER/q+hEdp4HK0qVLkZWVBWNjY6xfvx4nT57EqVOnsG7dOpiYmCArKwuenp5cVpEQQgip1wICAmBgYAA5OTm0a9cOV65cKTdveHg4+vbtCw0NDTRp0gSdOnXC2bNna7G2gjgNVLS0tBATEwMLCwt4eHhg+PDhGDZsGJYsWQILCwtcu3YNWlpaXFaREEIIqTauWlTCwsLg6uoKT09P3L9/H127dsWAAQOQmJgoNP+///6Lvn37IioqCnfv3kXPnj0xePBg3L9/v7qnQGQ8hmEYrl781atXMDAwAI/HQ3Z2Nl68eAGGYWBkZAQVFRWuqkUAtGzZkusq1EsvX75ETk4O19Wol5SVlbmuQr1H155oauPae/nypVjKqep3c4cOHdC2bVvs2LGDTTMzM8OwYcPg4+NTqTLMzc0xZswYLF++vEqvLS6ctqgYGRmxs3pUVFSwYcMG6OnpUZBCCCGECFFQUIAPHz7wPQoKCoTmLSwsxN27d2Fvb8+Xbm9vj5iYmEq9XmlpKT5+/AhVVdVq111UnAYq3zfmREVFIS8vj6PaEEIIITVDXF0/Pj4+UFJS4nuU1zKSkZGBkpISgSEUWlpaeP/+faXqvXHjRuTl5WH06NHVPgeionv9EEIIITVMXDN2PDw8MH/+fL40WVnZKr02wzCVqk9ISAhWrFiBkydPQlNTs+qVFRNOAxVhg4Ma6vQrQgghPy9x/W+TlZX9YWBSRl1dHZKSkgKtJ2lpaT+cqBIWFoYpU6bg6NGj6NOnj8j1FQdOAxWGYTBx4kT2pH/58gUzZ86EoqIiX77w8HAuqkcIIYTUWzIyMmjXrh2io6MxfPhwNj06OhpDhw4td7+QkBBMnjwZISEhGDhwYG1UtUKcBiouLi58z8ePH89RTQghhJCfz/z58zFhwgTY2tqiU6dOCAwMRGJiImbOnAnga1fS27dvERQUBOBrkOLs7IwtW7agY8eObGuMvLw8lJSUODkGTgOVffv2cfnyhBBCSK3galjDmDFjkJmZCW9vb6SkpMDCwgJRUVHQ09MDAKSkpPCtqfLnn3+iuLgYc+bMwZw5c9h0FxcX7N+/v7arD4DjdVRI3UXrqIiG1lERHa2jUn107YmmNq69N2/eiKWcsgCjIaFZP4QQQkgNo4kiouN0HRVCCCGEkIpQoEIIIYSQOou6fgghhJAaRl0/oqPBtIQQQkgNS05OFks5urq6YimnPvlpW1QyMzO5rkK9paamRrMHRKSsrIzCwkKuq1EvycjIoLS0lOtq1FsSEhLo0aMH19Woly5dusR1FUgFftpAhRBCCKkrqOtHdDSYlhBCCCF1FrWoEEIIITWMWlRERy0qhBBCCKmzqEWFEEIIqWHUoiK6OhGoFBQU4NatW0hISEB+fj40NDRgY2MDAwMDrqtGCCGEVBsFKqLjNFCJiYnBtm3bEBERgcLCQigrK0NeXh5ZWVkoKCiAoaEhpk+fjpkzZ6Jx48ZcVpUQQgghHOBsjMrQoUMxcuRINGvWDGfPnsXHjx+RmZmJ5ORk5Ofn4/nz51i6dCnOnz8PY2NjREdHc1VVQgghhHCEsxYVe3t7HD16FDIyMkK3GxoawtDQEC4uLoiNjcW7d+9quYaEEEKIeFDXj+g4C1TmzJlT6bzm5uYwNzevwdoQQgghpC6qE4Npv/XlyxeEhYUhLy8Pffv2hZGREddVIoQQQqqFWlREx2mgsmjRIhQWFmLLli0AgMLCQnTq1AmxsbFQUFCAm5sboqOj0alTJy6rSQghhBCOcLrg299//43evXuzzw8fPow3b97g+fPnyM7OxqhRo/DHH39wWENCCCGEcInTQCUxMRGtW7dmn587dw4jR46Enp4eeDwefv/9d9y/f5/DGhJCCCHVx+PxxPJoiDgNVCQkJMAwDPv8xo0b6NixI/tcWVkZ2dnZXFSNEEIIERsKVETHaaBiamqKv/76CwAQGxuLxMRE9OzZk93+5s0baGlpcVU9QgghRCwoUBEd54Npx44di8jISMTGxsLBwYFv2fyoqCi0b9+ewxoSQgghhEuctqg4OjoiKioKVlZWmDdvHsLCwvi2KygoYPbs2RzVjhBCCCFc4zHfDhL5iWRmZnJdhXpLTU0NOTk5XFejXlJWVkZhYSHX1aiXZGRkUFpaynU16i0JCQn06NGD62rUS5cuXarx18jKyhJLOaqqqmIppz7htOvn33//FZqupKSEVq1aQVFRsZZrRAghhJC6hNNApaLoX1JSErNmzcLGjRshLS1de5UihBBCxKyhDoQVB04DlfKmHufk5ODWrVtYtGgRtLW1sWTJklquGSGEEELqAk4DFSUlpXLT9fT0ICMjgyVLllCgQgghhDRQde6mhN9q06YN3rx5w3U1CCGEkGqhrh/R1elA5d27d9DU1OS6GoQQQki1UKAiOk7XUalIWloali5dil69enFdFUIIIYRwhNMWFRsbG6FRZm5uLpKTk2FmZobQ0FAOakYIIYSQuoDTQGXYsGFC05s0aQJTU1PY29tDUlKyditFCCGEiBl1/YiO00DFy8uLy5cnhBBCSB3H6RiVvXv3oqCggMsqVNvx48fh6OiIHj16YNKkSXjw4EG5eTMyMuDl5QUnJyd07twZfn5+AnkiIyNhZ2cn8Kjv56k8x44dw7Bhw9C1a1c4Ozvj/v375ebNyMjAsmXLMGrUKHTs2BGbNm2qsOxz586hQ4cOWLRokbirzYnQ0FD0798f7dq1w+jRo3H37t0K89++fRujR49Gu3bt0L9/fxw5ckQgz8GDBzF48GDY2tqiT58+WLduXbnX2u7du2FpaYl169aJ5XhqU3BwMPr06YM2bdrA0dERd+7cqTD/rVu34OjoiDZt2qBv375Cu6A/fPgAb29vdO3aFW3atMHAgQNx+fJldntISAiGDh0KW1tb2NrawsnJqdzVuOuboUOHIiQkBOfOncOff/4JS0vLcvNaW1vj0qVLAo8WLVqwefz8/ITm8fHxqY3DIXUcpy0q06ZNw6BBg9iZPU2bNkVMTAz09fW5rFal/fPPP9iyZQsWLlwIKysrREREYMGCBTh8+DC0tbUF8hcVFUFZWRkuLi4Vjr1RVFQU2C4rKyv2+nMtOjoamzdvhpubG6ysrHDixAnMmzcPoaGhQs9fYWEhlJWVMWnSJISEhFRYdkpKCrZu3Qpra+saqn3tOnPmDNatW4elS5fCxsYGR48exaxZs3Dy5Eno6OgI5E9OTsacOXPg6OiItWvX4v79+/jjjz+goqKCvn37AgBOnz4NPz8/eHt7w9raGm/evMHSpUsBAO7u7nzlPXnyBMeOHYOxsXHNH6yYRUVFYe3atVi2bBnatm2LsLAwzJgxA3/99ReaNm0qkD85ORkzZ87EyJEjsX79ety7dw+rVq2Cqqoq7O3tAXy9FqdMmQJVVVVs2bIFWlpaeP/+Pd9tP7S1tTF//nz2H/LJkycxd+5cHD9+HEZGRrVz8DWgZ8+emDt3Lvz8/PD48WMMGTIE69evh4uLC9LS0srdb/z48cjPz2eff3s/sWXLlvGtQN6kSRPs2bOHL/Cr76jrR3Sctqh8fz/Ejx8/1qubkoWGhmLw4MEYMmQI9PX14erqCk1NTZw4cUJofh0dHcybNw8DBgxAo0aNyi2Xx+NBTU2N7/EzCgkJwZAhQzB06FAYGBhg/vz50NLSwvHjx4Xmb9q0KRYsWAAHB4cKz19JSQmWL1+O6dOno1mzZjVV/VoVFBSEESNGwNHREYaGhnB3d4e2trbAHcfLHDlyBNra2nB3d4ehoSEcHR0xfPhw7N+/n83z8OFD2NjYYODAgWjWrBns7OwwYMAAxMXF8ZWVn5+PxYsXw8vLC02aNKnJw6wRBw4cwIgRIzBq1Ci0bNkSS5Ysgba2drk/FkJDQ6Gjo4MlS5agZcuWGDVqFEaMGIG9e/eyecLDw5Gbm4vt27ejbdu2aNasGdq1awdTU1M2T8+ePdG9e3cYGBjAwMAArq6uUFBQwMOHD2v8mGvSqFGjEBUVhcjISCQmJmL79u1IS0vD0KFDK9wvJycHWVlZ7OPb7/qPHz/ybbO1tcWXL19q5WaBpO6rs9OT67qioiI8ffoU7du350tv3749Hj9+XK2yP3/+jOHDh2Po0KFYuHAhnj59Wq3y6qKioiLEx8ejQ4cOfOniOH979uyBiooKhgwZUq1y6oqioiLExcXBzs6OL93Ozq7crsaHDx8K5O/cuTPi4uJQVFQEAGjbti3i4uLY852UlIQrV66ga9eufPutXr0aXbt2RadOncR0RLWnsLAQsbGx6Ny5M196586dy+1mfPDggdD8sbGx7Lm7cOECrK2tsWrVKnTp0gWDBw/Gn3/+iZKSEqFllpSUIDIyEvn5+fW6lU9KSgomJia4ffs2X/rt27dhbm5e4b67du3C8ePHsXHjxh+eAwcHB1y4cAFfvnypbpXJT4DTrh8ej8fXHPb987osJycHJSUlArfcVlVVrdbtvPX09ODp6YmWLVsiLy8PR44cwcyZMxEUFITmzZtXt9p1RnnnT01NDTdu3BC53IcPH+LUqVM4dOhQdatYZ2RnZ6OkpESgZU1NTQ2ZmZlC98nMzBSav7i4GDk5OdDQ0MCAAQOQlZUFZ2dnAEBxcTHGjBmDqVOnsvv8/fffiIuLq7fLBJRdZ+rq6nzpampqyMjIELpPRkaGwLlTV1dHcXExsrOzoampieTkZNy8eRODBg3Cn3/+iYSEBKxatQrFxcWYM2cOu9+zZ88wduxYFBQUQEFBAdu2bUOrVq3Ef6C1RElJCZKSkgL3acvOzhb4LJfJzMyEr68vnj17Bmlpadjb22PTpk1wdXXFo0ePBPKbmprC0NAQ69evr5Fj4Ep9+d9WF3EaqDAMA2NjY/YN/PTpE2xsbCAhwd/QU9E//oKCAoHBf1yO5/i+O6uqLCwsYGFhwT63srLCpEmTcPToUcyfP7+61atzvv/wMgwj8gc6Ly8PXl5eWLJkCZSVlcVQu7rtR9easHP7bfrt27exa9cuLF26FJaWlkhKSsLatWuhrq6OmTNn4v3791i7di0CAwN/ujFSP7rOfnTuSktLoaamBm9vb0hKSsLc3Bzp6enYs2cPX6Cir6+P8PBwfPz4EefOnYOHhweCgoLqdbACCF57FZ3LpKQkJCUlsc/j4uKgqamJMWPGCA1UHBwc8OrVK8THx4uvwqRe4zRQ2bdvX7XL8PHxwcqVK/nSvLy88Ouvv1a77IooKytDUlJSIIiq6JeFKCQkJGBqaork5GSxlVkXlJ2/71sEsrKyRD5/b9++RUpKChYuXMimlfWD29nZ4ciRI9DV1RW90hxRUVEp91yVN35JWItBVlYWpKSk2JuBbt++HYMHD4ajoyMAwNjYGPn5+fD29sb06dMRGxuLrKwsjBkzhi2jpKQEd+/eRUhICO7evVvn1zkqu86EnYvyzp26urpA/szMTEhJSbEBsIaGBqSkpPiO39DQEBkZGSgsLISMjAwAQEZGBnp6egC+/gh5/PgxDh48KPCdVV/k5uYKbQlVVlauUktyXFwcO6j7W7KysujVq5dY/jeQnwengYqLi0u1y/Dw8BBoaZCVlcWnT5+qXXZFpKWlYWJiglu3bqF79+5s+u3btwX6+KuDYRg8f/4cLVu2FFuZdYG0tDRMTU1x69Yt9OjRg02/desWunXrJlKZenp6CA4O5kvbuXMn8vPz2YG69ZG0tDRat26N69evo3fv3mz69evX0bNnT6H7tGnTRmDGRExMDFq3bs3Orvj8+bPAL2FJSUkwDAOGYdCxY0eEh4fzbV+2bBkMDAwwefLkOh+kAF8DBXNzc8TExPD9Y4yJiSn39hxl02m/de3aNZibm7Pnrm3btjh9+jRKS0vZFuCEhARoaGiwQUp5CgsLq3FE3CouLsbTp09ha2uLq1evsum2tra4du1apcsxMjIS2m3Zs2dPyMjIIDo6Wiz1rUuo60d0dfqmhN8qr6lWVlZWaLN0TQcqAODk5ARvb2+YmZnBwsICJ0+eRGpqKrvi7o4dO5Ceno7ly5ez+zx79gzA138SOTk5bL+tgYEBgK8DQc3NzdG8eXPk5eXh6NGjeP78OV8rwc9i7NixWLFiBUxNTWFpaYmIiAikpqZixIgRAAB/f3+kp6djxYoV7D5l5y8/P589f1JSUjA0NISsrKxAQNe4cWMAqPeBnrOzMzw8PGBubo42bdrg6NGjSElJwejRowF8XYciLS0Na9asAQCMHj0aoaGhWL9+PUaOHImHDx8iPDycr9+/R48eCAoKgpmZGSwtLdkZHD169ICkpCQUFRUFptHKy8tDWVm5Xk2vdXFxweLFi2FhYQFra2scOXIEKSkpbEvRpk2bkJqayq4P4+TkhODgYKxduxajRo3CgwcPEB4ejg0bNrBlOjk54dChQ1izZg3GjRuHN2/eIDAwEOPHj2fzbN68GV27doWOjg7y8vIQFRWFW7duITAwsHZPgJgdPXoUS5YswdOnTxEbG4vBgwdDS0sLp06dAvB12Ql1dXV2DZSRI0fi/fv3eP36NaSlpdG3b190794dy5YtEyjbwcEBV69exYcPH2r1mEjdxlmgYmZmhmXLlmHkyJEV/gJ5/vw5Nm3aBD09PSxevLgWa/hjffr0QW5uLvbu3YvMzEwYGhpiw4YN7LoWmZmZSE1N5dtn4sSJ7N/x8fE4d+4ctLW12V+unz59wrp165CVlQVFRUUYGxsjICAArVu3rrXjqi19+/Zlz19GRgYMDQ2xefPmCs/fhAkT2L/j4+Nx9uxZ6OjoICIiojarXuv69++PnJwc7Ny5E+np6WjVqhUCAgLYdUDS09ORkpLC5tfV1YW/vz98fX0RGhoKTU1NeHh48LUqTJ8+HTweD9u2bUNaWhpUVFTQvXt3/Pbbb7V+fDXJwcEBOTk5CAgIQHp6OoyMjLBz50526rqwc7dz506sXbsWwcHB0NTUxJIlS9g1VICvSw3s3r0ba9euxbBhw6ClpYUJEybwDUTOyMiAu7s70tPT0bhxYxgbGyMwMFBgRlF9c/HiRTRp0gQuLi5QVVXF69ev4e7uzn5W1dTU+FovpaSkMGvWLKirq6OgoAAJCQlwd3fHzZs3+crV1dWFlZUVFixYUKvHQ+o+HlPd0Z8iunDhAtzd3fHixQvY29vD1tYWTZs2hZycHLKzsxEXF4erV68iLi4Oc+fOxZIlS6q0hkN5syHIj6mpqfEtxkQqT1lZuV437XNJRkamXq2jVNdISEjwdaOSyquN9Vq+XeyuOhQUFMRSTn3C2ToqvXr1wu3btxEZGQltbW0EBwdj7ty5GDduHFasWIHnz5/D2dkZycnJWLt2bb1caIoQQgjhWkBAAAwMDCAnJ4d27drhypUr5eZNSUnBL7/8AhMTE0hISMDV1bX2KloOzseolN3LhhBCCPlZcTWYNiwsDK6urggICEDnzp3x559/sitQf3u/pTIFBQXQ0NCAp6cnNm/ezEGNBdHKtIQQQshPatOmTZgyZQqmTp0KMzMz+Pn5oXnz5tixY4fQ/Pr6+tiyZQucnZ3ZpQy4xmmLytatWyuV72cb3EcIIYSIorxFToXNfi0sLMTdu3cFJqLY29sjJiamRuspTpwGKpVpVuLxeBSoEEIIqdfE1fVT3iKn3y7jUCYjIwMlJSUCa0iV3e27vuA0UHn9+jWXL08IIYTUK+UtcloRcd6qhAucD6YlhBBCSOWU180jjLq6OiQlJQVaT9LS0urVSt2cBipBQUGVyld2d1dCCCGEVI6MjAzatWuH6OhoDB8+nE2Pjo7G0KFDOaxZ1XAaqEycOBGNGjWClJRUuXeC5fF4FKgQQgghIpg/fz4mTJgAW1tbdOrUCYGBgUhMTMTMmTMBfO1Kevv2LV/DwYMHDwB8XSk9PT0dDx48gIyMDGcrpHMaqJiZmSE1NRXjx4/H5MmTYWVlxWV1CCGEkBrB1ZiQMWPGIDMzE97e3khJSYGFhQWioqLYu3qnpKQgMTGRbx8bGxv277t3/6+9+4+Jsn7gAP5+FORu3h0CC9DBAUnEQKgwoGPFkh/d8Lohsz8olpAlraHVXFm40mwU/XLUYjFsC3HJSDPINcOEBHJGAfHDSBQUVuoBlYHBgvjx9Ifzki9Zfg+6z/PE+7XdH889x7P3ffSP9z7P5/lcM8rKyhAQEIDe3l5nRrcTtoX+FV999RXee+89fPDBBwgODsbDDz+MjIyMWe9Eyy30Hcct9B3HLfQdxy30Z4db6DvOGVvo/+8jxY663vUp/yXCN3yLjY1FcXExbDYbHn/8cezbtw9Lly5FRkbGnP3DEhERkToJLypXaLVarFu3Djt27EBMTAzKy8vn7EeciIiIRJIkaU5e85Eiisr58+fx8ssv46abbkJ6ejqio6PR0dEBDw8P0dGIiIhIIKGLafft24eSkhLU1dXBbDZj586dsFgsWLhwochYREREc2q+zobMBaGLaRcsWACj0YiMjIy/3XzGkS30uZjWcVxM6zgupnUcF9PODhfTOs4Zi2nHx8fn5Dqurq5zch01ETqjYjQaIUkSysrKrvkZ/tYPERHR/CW0qIh6JpuIiMiZeOvHcUKLyujoKKqrq3HvvfcCuLxD3tWPJLu4uODFF1+ERqMRFZGIiIgEElpUSktL8cknn9iLSmFhIcLDw6HVagEAnZ2d8PX1nfFLkURERDQ/CH08ee/evVi/fv2098rKynD06FEcPXoUr7/+Ovbv3y8oHRER0dzgPiqOE1pUTp8+jZCQEPuxRqPBggV/RoqJicF3330nIhoREREpgNBbP0NDQ3Bx+TPCjz/+OO381NQUt9EnIiKax4TOqPj5+eHbb7+95vn29nb4+fk5MREREdHc460fxwktKqtXr8a2bdswOjo649xvv/2GHTt2wGKxCEhGRERESiB0Z9r+/n7ceuutWLRoETZu3IiQkBBIkoTOzk4UFhZiYmICLS0tf7tr7bVwZ1rHcWdax3FnWsdxZ9rZ4c60jnPGzrRz9X/76nWc84XQNSo+Pj44fvw4HnvsMTz77LO40pkkSUJycjLeeecdh0oKERGRkszX2zZzQWhRAYCgoCBUVVXh4sWL6O7uBgAEBwfD09NTcDIiIiISTXhRucLT0xMxMTGiYxAREZGCKKaoEBER/Vfx1o/j5t+qHCIiIlINFhUiIiJSLKGPJ89HY2NjyM/PR25uLtzc3ETHURWO3exw/BzHsXMcx45mi0XFyS5dugR3d3cMDQ3BYDCIjqMqHLvZ4fg5jmPnOI4dzRZv/RAREZFisagQERGRYrGoEBERkWKxqDiZm5sbtm/fzkVlDuDYzQ7Hz3EcO8dx7Gi2uJiWiIiIFIszKkRERKRYLCpERESkWCwqREREpFgsKkRERKRYLCpOUl9fD6vVimXLlkGSJFRWVoqOpBr5+fmIjo6GXq+Ht7c31qxZg1OnTomOpQpFRUWIjIyEwWCAwWCAyWTCp59+KjqWKuXn50OSJDz55JOio6jCCy+8AEmSpr18fX1FxyIVYlFxkpGREdxyyy0oLCwUHUV16urqkJOTg4aGBhw5cgQTExO45557MDIyIjqa4vn5+eGVV15BU1MTmpqakJCQgNTUVHR0dIiOpiqNjY3YtWsXIiMjRUdRlfDwcNhsNvvrxIkToiORCrmIDjBfpKSkICUlRXQMVaqqqpp2XFJSAm9vbzQ3NyM+Pl5QKnWwWq3Tjl966SUUFRWhoaEB4eHhglKpy/DwMDIyMvDuu+8iLy9PdBxVcXFx4SwKzRpnVEh1hoaGAACenp6Ck6jL5OQkysvLMTIyApPJJDqOauTk5MBisSApKUl0FNXp6urCsmXLEBQUhPT0dJw9e1Z0JFIhzqiQqsiyjM2bN+POO+/EihUrRMdRhRMnTsBkMmF0dBQ6nQ4VFRUICwsTHUsVysvL8c0336CxsVF0FNWJjY3Fnj17EBISgv7+fuTl5SEuLg4dHR3w8vISHY9UhEWFVGXjxo1ob2/HsWPHREdRjZtvvhmtra0YHBzEgQMHkJmZibq6OpaVf/DDDz/giSeewGeffQaNRiM6jupcfas7IiICJpMJy5cvR2lpKTZv3iwwGakNiwqpxqZNm3Dw4EHU19fDz89PdBzVWLRoEYKDgwEAt99+OxobG/HWW2+huLhYcDJla25uxsDAAFauXGl/b3JyEvX19SgsLMTY2BgWLlwoMKG6LF68GBEREejq6hIdhVSGRYUUT5ZlbNq0CRUVFaitrUVQUJDoSKomyzLGxsZEx1C8xMTEGU+pPPTQQwgNDcUzzzzDkvJ/Ghsbw8mTJ3HXXXeJjkIqw6LiJMPDw+ju7rYf9/T0oLW1FZ6enjAajQKTKV9OTg7Kysrw8ccfQ6/Xo6+vDwDg7u4OrVYrOJ2ybd26FSkpKfD398evv/6K8vJy1NbWzniSimbS6/Uz1kEtXrwYXl5eXB91HZ566ilYrVYYjUYMDAwgLy8Ply5dQmZmpuhopDIsKk7S1NSEVatW2Y+v3KPNzMzE7t27BaVSh6KiIgDA3XffPe39kpISZGVlOT+QivT39+PBBx+EzWaDu7s7IiMjUVVVheTkZNHR6D/u3LlzuP/++/HTTz/hhhtuwB133IGGhgYEBASIjkYqI8myLIsOQURERPRXuI8KERERKRaLChERESkWiwoREREpFosKERERKRaLChERESkWiwoREREpFosKERERKRaLChERESkWiwqRwg0MDODRRx+F0WiEm5sbfH19YTab8eWXXwIAJElCZWWl2JBERP8SbqFPpHBr167F+Pg4SktLceONN6K/vx81NTW4ePGi6GhERP86zqgQKdjg4CCOHTuGV199FatWrUJAQABiYmKQm5sLi8WCwMBAAEBaWhokSbIfnzlzBqmpqfDx8YFOp0N0dDSqq6unXdtms8FisUCr1SIoKAhlZWUIDAzEm2++af/M0NAQsrOz4e3tDYPBgISEBLS1tTnp2xMRsagQKZpOp4NOp0NlZSXGxsZmnG9sbARw+QcabTab/Xh4eBirV69GdXU1WlpaYDabYbVa8f3339v/dt26dbhw4QJqa2tx4MAB7Nq1CwMDA/bzsizDYrGgr68Phw4dQnNzM6KiopCYmMjZHCJyHpmIFO3DDz+UPTw8ZI1GI8fFxcm5ublyW1ub/TwAuaKi4h+vExYWJr/99tuyLMvyyZMnZQByY2Oj/XxXV5cMQC4oKJBlWZZrampkg8Egj46OTrvO8uXL5eLi4tl/MSKi68AZFSKFW7t2LS5cuICDBw/CbDajtrYWUVFR2L179zX/ZmRkBFu2bEFYWBiWLFkCnU6Hzs5O+4zKqVOn4OLigqioKPvfBAcHw8PDw37c3NyM4eFheHl52Wd2dDodenp6cObMmX/t+xIRXY2LaYlUQKPRIDk5GcnJydi2bRseeeQRbN++HVlZWX/5+aeffhqHDx/GG2+8geDgYGi1Wtx33334/fffAVy+rfNXrn5/amoKS5cuRW1t7YzPLVmyZLZfiYjourCoEKlQWFiY/ZFkV1dXTE5OTjv/xRdfICsrC2lpaQAur1np7e21nw8NDcXExARaWlqwcuVKAEB3dzcGBwftn4mKikJfXx9cXFzsi3SJiJyNt36IFOznn39GQkIC3n//fbS3t6Onpwf79+/Ha6+9htTUVABAYGAgampq0NfXh19++QXA5ds4H330EVpbW9HW1oYHHngAU1NT9uuGhoYiKSkJ2dnZ+Prrr9HS0oLs7GxotVpIkgQASEpKgslkwpo1a3D48GH09vbi+PHjeO6559DU1OT8wSCieYlFhUjBdDodYmNjUVBQgPj4eKxYsQLPP/88NmzYgMLCQgDAzp07ceTIEfj7++O2224DABQUFMDDwwNxcXGwWq0wm83T1qMAwJ49e+Dj44P4+HikpaVhw4YN0Ov10Gg0AC5vJHfo0CHEx8dj/fr1CAkJQXp6Onp7e+Hj4+PcgSCieUuSr3WzmojmlXPnzsHf3x/V1dVITEwUHYeICACLCtG89fnnn2N4eBgRERGw2WzYsmULzp8/j9OnT8PV1VV0PCIiAFxMSzRvjY+PY+vWrTh79iz0ej3i4uKwd+9elhQiUhTOqBAREZFicTEtERERKRaLChERESkWiwoREREpFosKERERKRaLChERESkWiwoREREpFosKERERKRaLChERESkWiwoREREp1h9qMeod/WM5iAAAAABJRU5ErkJggg==",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df = pd.DataFrame(all_dicts)\n",
"biomarker_stage_probability_df = get_biomarker_stage_probability(df, num_biomarkers)\n",
"sns.heatmap(biomarker_stage_probability_df, annot=True, cmap=\"Greys\", linewidths=.5, cbar_kws={'label': 'Probability'})\n",
"plt.xlabel('Stage')\n",
"plt.ylabel('Biomarker')\n",
"plt.title('Heatmap of Biomarkers Involvement in Stages')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Estimate Theta and Phi Based on Conjugate Priors"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def estimate_params_exact(m0, n0, s0_sq, v0, data):\n",
" '''This is to estimate means and vars based on conjugate priors\n",
" Inputs:\n",
" - data: a vector of measurements \n",
" - m0: prior estimate of $\\mu$.\n",
" - n0: how strongly is the prior belief in $m_0$ is held.\n",
" - s0_sq: prior estimate of $\\sigma^2$.\n",
" - v0: prior degress of freedome, influencing the certainty of $s_0^2$.\n",
" \n",
" Outputs:\n",
" - mu estiate, std estimate\n",
" '''\n",
" # Data summary\n",
" sample_mean = np.mean(data)\n",
" sample_size = len(data)\n",
" sample_var = np.var(data, ddof=1) # ddof=1 for unbiased estimator\n",
"\n",
" # Update hyperparameters for the Normal-Inverse Gamma posterior\n",
" updated_m0 = (n0 * m0 + sample_size * sample_mean) / (n0 + sample_size)\n",
" updated_n0 = n0 + sample_size\n",
" updated_v0 = v0 + sample_size \n",
" updated_s0_sq = (1 / updated_v0) * ((sample_size - 1) * sample_var + v0 * s0_sq + \n",
" (n0 * sample_size / updated_n0) * (sample_mean - m0)**2)\n",
" updated_alpha = updated_v0/2\n",
" updated_beta = updated_v0*updated_s0_sq/2\n",
"\n",
" # Posterior estimates\n",
" mu_posterior_mean = updated_m0\n",
" sigma_squared_posterior_mean = updated_beta/updated_alpha\n",
"\n",
" mu_estimation = mu_posterior_mean\n",
" std_estimation = np.sqrt(sigma_squared_posterior_mean)\n",
"\n",
" return mu_estimation, std_estimation"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def get_estimated_means_stds_df(biomarkers, data_we_have, theta_phi_kmeans):\n",
" '''To get estimated parameters, returns a Pandas DataFrame\n",
" Input:\n",
" - biomarkers: biomarkers \n",
" - data_we_have: participants data filled with initial or updated participant_stages\n",
"\n",
" Output: \n",
" - estimate_means_std_df, just like means_stds_df, containing the estimated mean and std_dev for \n",
" distribution of biomarker values when the biomarker is affected and not affected\n",
"\n",
" Note that, there is one bug we need to fix: Sometimes, data_full might have only one observation or no ob\n",
" '''\n",
" # empty list of dictionaries to store the estimates \n",
" means_stds_estimate_dict_list = []\n",
" \n",
" for biomarker in biomarkers: \n",
" dic = {'biomarker': biomarker} # Initialize dictionary outside the inner loop\n",
" for affected in [True, False]:\n",
" data_full = data_we_have[(data_we_have.biomarker == biomarker) & (\n",
" data_we_have.affected == affected)]\n",
" if len(data_full) > 1:\n",
" measurements = data_full.measurement\n",
" mu_estimate, std_estimate = estimate_params_exact(\n",
" m0 = 0, n0 = 1, s0_sq = 1, v0 = 1, data=measurements)\n",
" if affected:\n",
" dic['theta_mean'] = mu_estimate\n",
" dic['theta_std'] = std_estimate\n",
" else:\n",
" dic['phi_mean'] = mu_estimate\n",
" dic['phi_std'] = std_estimate\n",
" # If there is only one observation or not observation at all, resort to theta_phi_kmeans\n",
" # YES, IT IS POSSIBLE THAT DATA_FULL HERE IS NULL\n",
" # For example, if a biomarker indicates stage of (num_biomarkers), but all participants' stages\n",
" # are smaller than that stage; so that for all participants, this biomarker is not affected\n",
" else:\n",
" # DONT FORGTE RESET_INDEX; this because you are acessing [0]\n",
" theta_phi_kmeans_biomarker_row = theta_phi_kmeans[\n",
" theta_phi_kmeans.biomarker == biomarker].reset_index(drop=True)\n",
" if affected:\n",
" dic['theta_mean'] = theta_phi_kmeans_biomarker_row['theta_mean'][0]\n",
" dic['theta_std'] = theta_phi_kmeans_biomarker_row['theta_std'][0]\n",
" else:\n",
" dic['phi_mean'] = theta_phi_kmeans_biomarker_row['phi_mean'][0]\n",
" dic['phi_std'] = theta_phi_kmeans_biomarker_row['phi_std'][0]\n",
" # print(f\"biomarker {biomarker} done!\")\n",
" means_stds_estimate_dict_list.append(dic)\n",
" estimate_means_stds_df = pd.DataFrame(means_stds_estimate_dict_list)\n",
" return estimate_means_stds_df "
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def add_kj_and_affected(data_we_have, participant_stages, num_participants):\n",
" '''This is to fill up data_we_have. \n",
" Basically, add two columns: k_j, and affected, based on the initial or updated participant_stages\n",
" Note that we assume here we've already got S_n\n",
"\n",
" Inputs:\n",
" - data_we_have\n",
" - participant_stages: np array \n",
" - participants: 0-99\n",
" '''\n",
" participant_stage_dic = dict(zip(np.arange(0, num_participants), participant_stages))\n",
" data_we_have['k_j'] = data_we_have.apply(lambda row: participant_stage_dic[row.participant], axis = 1)\n",
" data_we_have['affected'] = data_we_have.apply(lambda row: row.k_j >= row.S_n, axis = 1)\n",
" return data_we_have "
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def metropolis_hastings_unknown_theta_phi(\n",
" data_we_have, iterations, non_diseased_participants, burn_in, thining, theta_phi_kmeans):\n",
" num_participants = len(data_we_have.participant.unique())\n",
" num_biomarkers = len(data_we_have.biomarker.unique())\n",
" biomarker_names = np.array(list(data_we_have.biomarker.unique()))\n",
"\n",
" all_dicts = []\n",
"\n",
" # initialize an ordering and likelihood\n",
" # note that it should be a random permutation of numbers 1-10\n",
" best_order = np.random.permutation(np.arange(1, num_biomarkers+1))\n",
" biomarker_best_order_dic = dict(zip(biomarker_names, best_order))\n",
" best_likelihood = -np.inf \n",
"\n",
" # initialize participant_stages \n",
" # note that high should be num_biomarkers + 1; otherwise, no participants will be in the stage of 10\n",
" participant_stages = np.random.randint(low = 0, high = num_biomarkers + 1, size = num_participants)\n",
" participant_stages[non_diseased_participants] = 0\n",
"\n",
" for _ in range(iterations):\n",
" participant_stages_copy = participant_stages.copy()\n",
" # when we update best_order below,\n",
" # in each iteration, new_order will also update\n",
" new_order = best_order.copy()\n",
" # randomly select two indices\n",
" a, b = np.random.choice(num_biomarkers, 2, replace=False)\n",
" # swaping the order\n",
" new_order[a], new_order[b] = new_order[b], new_order[a]\n",
"\n",
" # likelihood of seeing all participants' data \n",
" # biomarker - order dict\n",
" ordering_dic = dict(zip(biomarker_names, new_order))\n",
" # fill up S_n column using the ordering dict\n",
" # copy first in order not to change data_we_have\n",
" data = data_we_have.copy()\n",
" # now data_we_have has S_n column\n",
" data['S_n'] = data.apply(lambda row: ordering_dic[row['biomarker']], axis = 1)\n",
"\n",
" # add kj and affected based on the initial randomized participant_stages\n",
" data = add_kj_and_affected(data, participant_stages_copy, num_participants)\n",
" # print(data.head())\n",
"\n",
" # get estimated_theta_phi\n",
" estimated_theta_phi = get_estimated_means_stds_df(biomarker_names, data_we_have=data, theta_phi_kmeans=theta_phi_kmeans)\n",
"\n",
" all_participant_ln_likelihood = 0 \n",
" for p in range(num_participants):\n",
" # this participant data\n",
" pdata = data[data.participant == p].reset_index(drop=True)\n",
"\n",
" \"\"\"If this participant is not diseased (i.e., if we know k_j is equal to 0)\n",
" We still need to compute the likelihood of this participant seeing this sequence of biomarker data\n",
" but we do not need to estimate k_j like below\n",
"\n",
" We still need to compute the likelihood because we need to add it to all_participant_ln_likelihood\n",
" \"\"\"\n",
" if p in non_diseased_participants:\n",
" this_participant_likelihood = compute_likelihood(\n",
" pdata, k_j = 0, theta_phi = estimated_theta_phi)\n",
" this_participant_ln_likelihood = np.log(this_participant_likelihood)\n",
" else:\n",
" # initiaze stage_likelihood\n",
" stage_likelihood = np.zeros(num_biomarkers + 1)\n",
" for k_j in range(num_biomarkers +1):\n",
" # even though data above has everything, it is filled up by random stages\n",
" # we don't like it and want to know the true k_j. All the following is to update participant_stages\n",
"\n",
" # likelihood for this participant to have this specific sequence of biomarker values\n",
" participant_likelihood = compute_likelihood(pdata, k_j, estimated_theta_phi)\n",
"\n",
" # update each stage likelihood for this participant\n",
" stage_likelihood[k_j] = participant_likelihood\n",
" likelihood_sum = np.sum(stage_likelihood)\n",
" normalized_stage_likelihood = [l/likelihood_sum for l in stage_likelihood]\n",
" sampled_stage = np.random.choice(np.arange(num_biomarkers + 1), p = normalized_stage_likelihood)\n",
" participant_stages_copy[p] = sampled_stage \n",
"\n",
" # if participant is in sampled_stage, what is the likelihood of seeing this sequence of biomarker data:\n",
" this_participant_likelihood = stage_likelihood[sampled_stage]\n",
"\n",
" # then, update all_participant_likelihood\n",
" if this_participant_likelihood == 0:\n",
" this_participant_ln_likelihood = np.log(this_participant_likelihood + 1e20)\n",
" else:\n",
" this_participant_ln_likelihood = np.log(this_participant_likelihood)\n",
" all_participant_ln_likelihood += this_participant_ln_likelihood\n",
" \n",
" acceptance_ratio = np.exp(all_participant_ln_likelihood - best_likelihood)\n",
" random_number = np.random.rand()\n",
" if random_number < acceptance_ratio:\n",
" best_likelihood = all_participant_ln_likelihood\n",
" biomarker_best_order_dic = ordering_dic\n",
" participant_stages = participant_stages_copy\n",
" \n",
" if _ >= burn_in and _ % thining == 0:\n",
" all_dicts.append(ordering_dic)\n",
"\n",
" if (_+1) % 10 == 0:\n",
" print(f\"iteration {_ + 1} done\")\n",
" return biomarker_best_order_dic, participant_stages, all_dicts"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"