diff --git a/convergence_test.ipynb b/convergence_test.ipynb new file mode 100644 index 00000000000..adaf3a1910f --- /dev/null +++ b/convergence_test.ipynb @@ -0,0 +1,1121 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setting up and running the TARDIS example simulation " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from tardis.base import run_tardis\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import numpy as np\n", + "import logging\n", + "import time\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d8df4fb6a2b14601ae1641d20ab5c96e", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Tab(children=(Output(layout=Layout(height='300px', overflow_y='auto')), Output(layout=Layout(height='300px', o…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d8eb68745d68421592fe03b32b9b1431", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Iterations: 0/? [00:00" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sim.iterations_t_rad[-2:]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iteration 19 values: [10952.63231877 11065.29655507 11146.8359289 11168.51126187\n", + " 11162.76266113 11173.84014227 11125.28131485 11103.5044254\n", + " 11057.89578328 10982.17284381 10916.74656888 10853.86934121\n", + " 10797.7863424 10741.66583036 10657.1863048 10579.63315852\n", + " 10512.26505999 10441.9998283 10359.22941609 10257.68110852] K\n", + "Iteration 20 values: [10974.413526 11079.82475989 11155.4144222 11184.4515512\n", + " 11180.5101573 11186.08014564 11145.73248128 11121.47797056\n", + " 11075.74310553 10993.74039084 10924.91062725 10853.81618167\n", + " 10798.87369014 10731.34443456 10658.81913311 10581.85516831\n", + " 10513.56980191 10441.74819424 10359.80427222 10261.24669165] K\n" + ] + } + ], + "source": [ + "last_two_shells = sim.iterations_t_rad[-2:]\n", + "\n", + "for i, iteration in enumerate(last_two_shells, start=19):\n", + " print(f\"Iteration {i} values: {iteration}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Percentage differences between last two shells\n", + "T_rad 0: 0.20%\n", + "T_rad 1: 0.13%\n", + "T_rad 2: 0.08%\n", + "T_rad 3: 0.14%\n", + "T_rad 4: 0.16%\n", + "T_rad 5: 0.11%\n", + "T_rad 6: 0.18%\n", + "T_rad 7: 0.16%\n", + "T_rad 8: 0.16%\n", + "T_rad 9: 0.11%\n", + "T_rad 10: 0.07%\n", + "T_rad 11: -0.00%\n", + "T_rad 12: 0.01%\n", + "T_rad 13: -0.10%\n", + "T_rad 14: 0.02%\n", + "T_rad 15: 0.02%\n", + "T_rad 16: 0.01%\n", + "T_rad 17: -0.00%\n", + "T_rad 18: 0.01%\n", + "T_rad 19: 0.03%\n" + ] + } + ], + "source": [ + "last_two_shells = sim.iterations_t_rad[-2:]\n", + "\n", + "differences = last_two_shells[1] - last_two_shells[0]\n", + "\n", + "percentage_differences = (differences / last_two_shells[0]) * 100\n", + "\n", + "print(\"Percentage differences between last two shells\")\n", + "for i, percentage_difference in enumerate(percentage_differences):\n", + " print(f\"T_rad {i}: {percentage_difference:.2f}%\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Analysis of damped convergence iterations in TARDIS (0.4)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "c949d2c71b8e4519a6a81aeb98b5c267", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Tab(children=(Output(layout=Layout(height='300px', overflow_y='auto')), Output(layout=Layout(height='300px', o…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "958762d98764479e94b8c21455bb17b5", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Iterations: 0/? [00:00" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sim.iterations_t_rad[-2:]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iteration 19 values: [10952.63231877 11065.29655507 11146.8359289 11168.51126187\n", + " 11162.76266113 11173.84014227 11125.28131485 11103.5044254\n", + " 11057.89578328 10982.17284381 10916.74656888 10853.86934121\n", + " 10797.7863424 10741.66583036 10657.1863048 10579.63315852\n", + " 10512.26505999 10441.9998283 10359.22941609 10257.68110852] K\n", + "Iteration 20 values: [10974.413526 11079.82475989 11155.4144222 11184.4515512\n", + " 11180.5101573 11186.08014564 11145.73248128 11121.47797056\n", + " 11075.74310553 10993.74039084 10924.91062725 10853.81618167\n", + " 10798.87369014 10731.34443456 10658.81913311 10581.85516831\n", + " 10513.56980191 10441.74819424 10359.80427222 10261.24669165] K\n" + ] + } + ], + "source": [ + "last_two_shells = sim.iterations_t_rad[-2:]\n", + "\n", + "for i, iteration in enumerate(last_two_shells, start=19):\n", + " print(f\"Iteration {i} values: {iteration}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Percentage differences between last two shells\n", + "T_rad 0: 0.20%\n", + "T_rad 1: 0.13%\n", + "T_rad 2: 0.08%\n", + "T_rad 3: 0.14%\n", + "T_rad 4: 0.16%\n", + "T_rad 5: 0.11%\n", + "T_rad 6: 0.18%\n", + "T_rad 7: 0.16%\n", + "T_rad 8: 0.16%\n", + "T_rad 9: 0.11%\n", + "T_rad 10: 0.07%\n", + "T_rad 11: -0.00%\n", + "T_rad 12: 0.01%\n", + "T_rad 13: -0.10%\n", + "T_rad 14: 0.02%\n", + "T_rad 15: 0.02%\n", + "T_rad 16: 0.01%\n", + "T_rad 17: -0.00%\n", + "T_rad 18: 0.01%\n", + "T_rad 19: 0.03%\n" + ] + } + ], + "source": [ + "differences = last_two_shells[1] - last_two_shells[0]\n", + "\n", + "percentage_differences = (differences / last_two_shells[0]) * 100\n", + "\n", + "print(\"Percentage differences between last two shells\")\n", + "for i, percentage_difference in enumerate(percentage_differences):\n", + " print(f\"T_rad {i}: {percentage_difference:.2f}%\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Analysis of damped convergence iterations in TARDIS (0.4)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "c949d2c71b8e4519a6a81aeb98b5c267", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Tab(children=(Output(layout=Layout(height='300px', overflow_y='auto')), Output(layout=Layout(height='300px', o…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "958762d98764479e94b8c21455bb17b5", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Iterations: 0/? [00:00estimated_value (now it works)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "492c4b3cd91b437083d655b001dd4dde", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Tab(children=(Output(layout=Layout(height='300px', overflow_y='auto')), Output(layout=Layout(height='300px', o…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "e3c3ae3aa2194169b5d6806c18084dd2", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Iterations: 0/? [00:00" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "logging.basicConfig(level=logging.INFO)\n", + "logger = logging.getLogger(__name__)\n", + "\n", + "# Defining the control function (cos(x))\n", + "def control_function(x):\n", + " return np.cos(x)\n", + "\n", + "# Defining the damped convergence function\n", + "def damped_converge(value, estimated_value, damping_factor):\n", + " return value + damping_factor * (estimated_value - value)\n", + "\n", + "# Checking the convergence status\n", + "def check_convergence_status(x, estimated_x, threshold):\n", + " convergence_value = abs(x - estimated_x) / abs(estimated_x)\n", + " return convergence_value < threshold\n", + "\n", + "# Example of running the convergence check on cos(x)\n", + "def run_convergence(x0, damping_factor=0.5, threshold=1e-6, max_iterations=100):\n", + " start_time = time.perf_counter() # Starting timing\n", + " iterations = 0\n", + " x_values = [x0] # Storing initial guess values\n", + " while iterations < max_iterations:\n", + " estimated_x = control_function(x0)\n", + " if check_convergence_status(x0, estimated_x, threshold):\n", + " end_time = time.perf_counter() # End timing\n", + " logger.info(f\"\\n\\tConverged in {iterations+1} iterations \"\n", + " f\"\\n\\tSimulation took {(end_time - start_time):.6f} s\\n\")\n", + " return estimated_x, x_values, iterations\n", + " x0 = damped_converge(x0, estimated_x, damping_factor)\n", + " x_values.append(x0) # Storing each new value of x\n", + " iterations += 1\n", + " end_time = time.perf_counter() # Ending timing\n", + " logger.info(f\"\\n\\tDid not converge in {iterations} iterations \"\n", + " f\"\\n\\tSimulation took {(end_time - start_time):.6f} s\\n\")\n", + " return estimated_x, x_values, iterations\n", + "\n", + "# Initial guess\n", + "x0 = 1.0\n", + "\n", + "# Running the convergence check and printing results \n", + "converged_value, x_values, iterations = run_convergence(x0)\n", + "print(f'Converged value: {converged_value} after {iterations} iterations')\n", + "\n", + "# Plotting the convergence process\n", + "iterations_range = np.arange(len(x_values))\n", + "plt.plot(iterations_range, x_values, marker='o', label='Damped Fixed Point Iteration')\n", + "plt.axhline(y=converged_value, color='r', linestyle='--', label='Converged Value')\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel('Value of x')\n", + "plt.title('Convergence of Damped Fixed-Point Iteration')\n", + "plt.legend()\n", + "plt.grid(True)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Undamped convergence for cos(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged in 34 iterations.\n", + "Convergence took 0.006329 seconds.\n", + "Fixed point: 0.7390848229131413\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHFCAYAAAAaD0bAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABoNklEQVR4nO3dd3yT1f4H8M+TNKMbaEsXpS3IpuyN7CVDEa4/qigKggq4EPRecbFUxIEgCjjAXrwqqAwRUamCyFLZGwUZZbRAC7SlM03O74+Qp02TtklI84T28369+pI+OTk5z0nafj3ne86RhBACRERERNWISukGEBEREXkaAyAiIiKqdhgAERERUbXDAIiIiIiqHQZAREREVO0wACIiIqJqhwEQERERVTsMgIiIiKjaYQBERERE1Q4DoGriwIEDGDNmDOLj46HX6xEQEIA2bdrgzTffxJUrV5RuHrnZihUr0KxZM/j6+kKSJOzbt89uuaSkJEiShF27dtl9fMiQIYiLi3Nr20aPHu32Or1Rz5490bNnzwrLxcXFQZIku1/Xr19XrL8sn43Tp09XWHb06NFW7dbpdGjUqBGmTZuG/Px8p19bkiRMnz7d+UYDWLhwIZKSkhwuHxcXhyFDhsjf5+bmYvr06fj1119den13Ka8dzrw3VDYfpRtAle/jjz/GxIkT0ahRIzz33HNo2rQpDAYDdu3ahcWLF2PHjh1YvXq10s0kN7l8+TJGjRqFO+64AwsXLoROp0PDhg2VbhaVo2vXrnj77bdtrvv5+eHll1/G008/rUCrnOPr64uNGzcCAK5evYovv/wSM2fOxLFjx7BixQqn6tqxYwfq1KnjUjsWLlyI0NBQjB492qXn5+bmYsaMGQDgUABbWcprx+DBg7Fjxw5ERkYq0LKqgwFQFbdjxw5MmDAB/fr1w5o1a6DT6eTH+vXrhylTpuDHH39UsIU3z2g0oqioyOreqrO///4bBoMBDzzwAHr06KF0c8gBNWrUQKdOnew+Vr9+fQ+3xjUqlcrqHgYOHIjTp0/jq6++wty5cxEdHe1wXWX1xa3MYDBAkiT4+Nz8n92wsDCEhYW5oVXVG6fAqrjXX38dkiTho48+shsgaLVa3HXXXfL3JpMJb775Jho3bgydTofatWvjwQcfxLlz56ye17NnTzRv3hw7d+5Et27d4Ofnh3r16uGNN96AyWQCYB6J0Gq1ePnll21e99ixY5AkCe+99558LS0tDY899hjq1KkDrVaL+Ph4zJgxA0VFRXKZ06dPQ5IkvPnmm3j11VcRHx8PnU6HTZs2AQC+/fZbtGjRAjqdDvXq1cP8+fMxffp0SJJk9fpCCCxcuBCtWrWCr68vatasiXvuuQcnT550+j4trl27hilTpqBevXpy3w0aNAjHjh2TyxQWFuLVV1+V+zcsLAxjxozB5cuX7b+BpaxduxadO3eGn58fAgMD0a9fP+zYsUN+fPTo0bj99tsBAImJiZAkya3/F2vp/7fffhtz585FfHw8AgIC0LlzZ/z+++825ZOSktCoUSPodDo0adIEy5Yts1vvjBkz0LFjR9SqVQtBQUFo06YNlixZgtJnNVumK9atW4fWrVvD19cXTZo0wbp16+TXa9KkCfz9/dGhQwebqb3Ro0cjICAAhw8fRp8+feDv74+wsDA88cQTyM3NtSrr6GdECIE333wTsbGx0Ov1aNOmDX744Qen+7YspafAli9fDkmS8P7771uVmzZtGtRqNZKTk+Vru3btwl133YVatWpBr9ejdevW+Oqrr2xe4/fff0fXrl2h1+sRFRWFqVOnwmAw3HTbLYHMmTNnAAApKSl44IEHULt2bfkz8c4779j8LJWeArNM+WzatAkTJkxAaGgoQkJCMHz4cFy4cEEuFxcXh8OHD2Pz5s3ydJwz04enT5+WA4sZM2bIdZQcTTp+/DhGjhxpdQ8ffPCBVT2//vorJEnCZ599hilTpiA6Oho6nQ4nTpzA5cuXMXHiRDRt2hQBAQGoXbs2evfujS1btjjcjrKmwJYuXYqWLVtCr9ejVq1aGDZsGI4ePWpVxvIzcOLECQwaNAgBAQGIiYnBlClTUFBQ4HBfVQmCqqyioiLh5+cnOnbs6PBzHn30UQFAPPHEE+LHH38UixcvFmFhYSImJkZcvnxZLtejRw8REhIiGjRoIBYvXiySk5PFxIkTBQDx3//+Vy43bNgwERMTI4xGo9Xr/Pvf/xZarVakp6cLIYRITU0VMTExIjY2Vnz44Yfi559/FrNmzRI6nU6MHj1aft6pU6cEABEdHS169eolvvnmG7FhwwZx6tQp8cMPPwiVSiV69uwpVq9eLb7++mvRsWNHERcXJ0p/1B955BGh0WjElClTxI8//ii++OIL0bhxYxEeHi7S0tKcvs+srCzRrFkz4e/vL2bOnCl++uknsXLlSvH000+LjRs3CiGEMBqN4o477hD+/v5ixowZIjk5WXzyySciOjpaNG3aVOTm5pb73nz++ecCgOjfv79Ys2aNWLFihWjbtq3QarViy5YtQgghTpw4IT744AMBQLz++utix44d4vDhw2XW+emnnwoAYufOnXYfHzx4sIiNjbXp/7i4OHHHHXeINWvWiDVr1oiEhARRs2ZNce3aNZu6hw4dKr777jvxv//9T9x2223y+1zS6NGjxZIlS0RycrJITk4Ws2bNEr6+vmLGjBlW5WJjY0WdOnVE8+bNxZdffinWr18vOnbsKDQajXjllVdE165dxapVq8Tq1atFw4YNRXh4uFW/PvTQQ0Kr1Yq6deuK1157TWzYsEFMnz5d+Pj4iCFDhli9lqOfkWnTpgkAYuzYseKHH34QH330kYiOjhYRERGiR48eZfZ9yXsaNGiQMBgMVl+Wn5mHHnrIpr/Gjx8vtFqt/L798ssvQqVSiZdeekkus3HjRqHVakW3bt3EihUrxI8//ihGjx4tAIhPP/1ULnf48GHh5+cnmjZtKr788kvx7bffigEDBoi6desKAOLUqVMV3sNDDz0k/P39ba4PGzZMABB///23uHTpkoiOjhZhYWFi8eLF4scffxRPPPGEACAmTJhg9TwAYtq0afL3ls9SvXr1xJNPPil++ukn8cknn4iaNWuKXr16yeX27Nkj6tWrJ1q3bi127NghduzYIfbs2VNu22NjY8XgwYOFEELk5+eLH3/8UX4/LXWcOHFC7qvg4GCRkJAgli1bJjZs2CCmTJkiVCqVmD59ulznpk2b5N9T99xzj1i7dq1Yt26dyMjIEMeOHRMTJkwQy5cvF7/++qtYt26dGDt2rFCpVGLTpk0OtcPSHyXfm9dff10AEPfdd5/4/vvvxbJly0S9evVEcHCw+Pvvv63eK61WK5o0aSLefvtt8fPPP4tXXnlFSJJk8/NW1TEAqsLS0tIEAHHvvfc6VP7o0aMCgJg4caLV9T/++EMAEC+88IJ8rUePHgKA+OOPP6zKNm3aVAwYMED+fu3atQKA2LBhg3ytqKhIREVFiX/961/ytccee0wEBASIM2fOWNX39ttvCwDyH3HLH+D69euLwsJCq7Lt27cXMTExoqCgQL6WnZ0tQkJCrAKgHTt2CADinXfesXr+2bNnha+vr/j3v//t9H3OnDlTABDJycmiLF9++aUAIFauXGl1fefOnQKAWLhwYZnPNRqNIioqSiQkJFgFk9nZ2aJ27dqiS5cu8jXLL9+vv/66zPosXA2AEhISRFFRkXz9zz//FADEl19+adXeNm3aCJPJJJc7ffq00Gg0Nn/QS9+rwWAQM2fOFCEhIVbPj42NFb6+vuLcuXPytX379gkAIjIyUuTk5MjX16xZIwCItWvXytceeughAUDMnz/f6jVfe+01AUBs3bpVCOH4Z+Tq1atCr9eLYcOGWZXbtm2bAOBwAATA5uvFF1+U21y6v/Lz80Xr1q1FfHy8OHLkiAgPDxc9evSwek8aN24sWrduLQwGg9VzhwwZIiIjI+XPUWJiovD19bUK6oqKikTjxo2dDoAswdvly5fF/PnzhSRJon379kIIIZ5//nm7P0sTJkwQkiSJv/76S75WVgBU+nfTm2++KQCI1NRU+VqzZs0c6neLkgGQEEJcvnzZ5vUtBgwYIOrUqSMyMzOtrj/xxBNCr9eLK1euCCGKfwa7d+9e4esXFRUJg8Eg+vTpY/U5Kq8dpQOgq1evCl9fXzFo0CCrcikpKUKn04mRI0fK1yw/A1999ZVV2UGDBolGjRpV2N6qhFNgJLNMI5VOHuzQoQOaNGmCX375xep6REQEOnToYHWtRYsW8nA3YM4DiIiIwKeffipf++mnn3DhwgU8/PDD8rV169ahV69eiIqKQlFRkfw1cOBAAMDmzZutXueuu+6CRqORv8/JycGuXbtw9913Q6vVytcDAgJw5513Wj133bp1kCQJDzzwgNVrRUREoGXLljarLhy5zx9++AENGzZE3759UZZ169ahRo0auPPOO61et1WrVoiIiCh31clff/2FCxcuYNSoUVCpin9sAwIC8K9//Qu///67zRROZRo8eDDUarX8fYsWLQAUT3VY2jty5Eir6cfY2Fh06dLFpr6NGzeib9++CA4OhlqthkajwSuvvIKMjAxcunTJqmyrVq2s8kmaNGkCwDxd6efnZ3O95Ptkcf/991t9P3LkSADFPwOOfkZ27NiB/Px8m/q6dOmC2NhYm9cty+23346dO3dafU2cOLHM8jqdDl999RUyMjLQpk0bCCHw5Zdfyu/JiRMncOzYMbldJe9h0KBBSE1NxV9//SXfc58+fRAeHi7Xr1arkZiYaPWaJpPJqh6j0Wj1eE5ODjQaDTQaDcLCwjBp0iQMHDhQXmCxceNGNG3a1OZnafTo0RBCyAnU5Sk5XQ/Yfu4qU35+Pn755RcMGzYMfn5+Nn2an59vMw38r3/9y25dixcvRps2baDX6+Hj4wONRoNffvnFZrrKUTt27EBeXp7N7+6YmBj07t3b5ne3JEk2vxdL/06rDpgEXYWFhobCz88Pp06dcqh8RkYGANhdWRAVFWXzwxESEmJTTqfTIS8vT/7ex8cHo0aNwoIFC3Dt2jXUqFEDSUlJiIyMxIABA+RyFy9exHfffWcV1JSUnp5u9X3pNl69ehVCCKtf4halr128eLHMsgBQr149p+/z8uXLqFu3rt36Sr7utWvXrAK0kkrfY0kVvTcmkwlXr161CgAcYUnILP3HzKKoqMjue1K6Tyz5ZZY+sbQ3IiLC5rkRERFWuQt//vkn+vfvj549e+Ljjz+Wc8DWrFmD1157zaqfAaBWrVpW31v6s6zrpZdh+/j42LTf0k5Lux39jFR0n44KDg5Gu3btHC4PALfddhu6deuG77//HhMmTLD6bFy8eBEA8Oyzz+LZZ5+1+3zL5y0jI8Oh9s+cOVNelQSYg9mS76Ovry9+++03AObPQ2xsLIKCguTHMzIy7ObjREVFyY9XpKLPXWXKyMhAUVERFixYgAULFtgtU9HvKQCYO3cupkyZgvHjx2PWrFkIDQ2FWq3Gyy+/7HIAVNHvh5J5YYB5daFer7e6ptPpXNqy4FbGAKgKU6vV6NOnD3744QecO3euwmWlll8uqampNmUvXLiA0NBQl9oxZswYvPXWW1i+fDkSExOxdu1aTJo0yWoEITQ0FC1atMBrr71mtw7LL0mL0knNNWvWhCRJ8i/+ktLS0qy+Dw0NhSRJ2LJli93EcFdWk4WFhdkkipdmSdwsa9VdYGBgmc8t+d6UduHCBahUKtSsWdOJFptZ/sCfP3/e7uPnz58vMwgoj6W9pfve3rXly5dDo9Fg3bp1Vr+U16xZ4/TrOqKoqAgZGRlWf0wtbbJcc/QzUtF9Vub+PZ988gm+//57dOjQAe+//z4SExPRsWNHuf0AMHXqVAwfPtzu8xs1agTAfA+OvE+PPvqo1X45pftFpVKVG8SFhISU+fkt2WZvVbNmTajVaowaNQqPP/643TLx8fFW35f+PQUA//vf/9CzZ08sWrTI6np2drbLbavo94O3961SOAVWxU2dOhVCCDzyyCMoLCy0edxgMOC7774DAPTu3RuA+Qe0pJ07d+Lo0aPo06ePS21o0qQJOnbsiE8//RRffPEFCgoKMGbMGKsyQ4YMwaFDh1C/fn20a9fO5qt0AFSav78/2rVrhzVr1ljd5/Xr1+UVQiVfSwiB8+fP232thIQEp+9x4MCB+Pvvv8sdxh8yZAgyMjJgNBrtvq7lD5I9jRo1QnR0NL744gurlVE5OTlYuXKlvDLMWZ06dUJAQIDdfVqOHDmCw4cPlzutV157IyMj8eWXX1q198yZM9i+fbtVWcvS4JIBcV5eHj777DOnX9dRn3/+udX3X3zxBYDi/VYc/Yx06tQJer3epr7t27dX6nTCwYMH8dRTT+HBBx/Eli1b0KJFCyQmJuLq1asAzP3foEED7N+/327727VrJwfcvXr1wi+//GL1Pw9Go9HmMxEVFXVTPyd9+vTBkSNHsGfPHqvry5YtgyRJ6NWrlytdYaP06KwrzwdsR5X8/PzQq1cv7N27Fy1atLDbp/ZGi0uzbBRZ0oEDB6xWc5bXDns6d+4MX19fm9/d586dw8aNG13+3V3VcQSoiuvcuTMWLVqEiRMnom3btpgwYQKaNWsGg8GAvXv34qOPPkLz5s1x5513olGjRnj00UexYMECqFQqeR+Pl19+GTExMXjmmWdcbsfDDz+Mxx57DBcuXECXLl1s/tjPnDkTycnJ6NKlC5566ik0atQI+fn5OH36NNavX4/FixdXOII1c+ZMDB48GAMGDMDTTz8No9GIt956CwEBAVa7XXft2hWPPvooxowZg127dqF79+7w9/dHamoqtm7dioSEBEyYMMGp+5s0aRJWrFiBoUOH4vnnn0eHDh2Ql5eHzZs3Y8iQIejVqxfuvfdefP755xg0aBCefvppdOjQARqNBufOncOmTZswdOhQDBs2zG79KpUKb775Ju6//34MGTIEjz32GAoKCvDWW2/h2rVreOONN5xqr0VgYCBmzJiBKVOmwGQyITExETVr1sTBgwfx+uuvIzY2Fk899ZTT9apUKsyaNQvjxo3DsGHD8Mgjj+DatWuYPn26zdTK4MGDMXfuXIwcORKPPvooMjIy8Pbbb1favk5arRbvvPMOrl+/jvbt22P79u149dVXMXDgQHkLAUc/IzVr1sSzzz6LV199FePGjcP//d//4ezZs3bv011ycnIwYsQIxMfHY+HChdBqtfjqq6/Qpk0bjBkzRh45+/DDDzFw4EAMGDAAo0ePRnR0NK5cuYKjR49iz549+PrrrwEAL730EtauXYvevXvjlVdegZ+fHz744APk5OS4td3PPPMMli1bhsGDB2PmzJmIjY3F999/j4ULF2LChAlu26wzISEBy5cvx4oVK1CvXj3o9XqngrXAwEDExsbi22+/RZ8+fVCrVi2EhoYiLi4O8+fPx+23345u3bphwoQJiIuLQ3Z2Nk6cOIHvvvvOoTymIUOGYNasWZg2bRp69OiBv/76CzNnzkR8fLzVlh/ltaO0GjVq4OWXX8YLL7yABx98EPfddx8yMjIwY8YM6PV6TJs2zeH7r1aUyr4mz9q3b5946KGHRN26dYVWqxX+/v6idevW4pVXXhGXLl2SyxmNRjFnzhzRsGFDodFoRGhoqHjggQfE2bNnrerr0aOHaNasmc3r2FuxIoQQmZmZwtfXVwAQH3/8sd02Xr58WTz11FMiPj5eaDQaUatWLdG2bVvx4osviuvXrwshilchvfXWW3brWL16tUhISJCXOr/xxhviqaeeEjVr1rQpu3TpUtGxY0fh7+8vfH19Rf369cWDDz4odu3a5dJ9Xr16VTz99NOibt26QqPRiNq1a4vBgweLY8eOyWUMBoN4++23RcuWLYVerxcBAQGicePG4rHHHhPHjx+3e08lrVmzRnTs2FHo9Xrh7+8v+vTpI7Zt22ZVxplVYBZfffWVuP3220VgYKDw8fERdevWFRMmTLBaGSRE+f0POytWPvnkE9GgQQOh1WpFw4YNxdKlS+323dKlS0WjRo2ETqcT9erVE7NnzxZLliyxWYVUesVOydd+/PHHK2yrZbXSgQMHRM+ePYWvr6+oVauWmDBhgvwZK92uij4jJpNJzJ49W8TExAitVitatGghvvvuO9GjRw+HV4HZu6eSbS7ZXw888IDw8/Oz2d7g66+/FgDEu+++K1/bv3+/GDFihKhdu7bQaDQiIiJC9O7dWyxevNjqudu2bROdOnUSOp1OREREiOeee0589NFHN70MvrQzZ86IkSNHipCQEKHRaESjRo3EW2+9ZbNNRunPUlmrFS2fdcvycSHMKw379+8vAgMDBYByVxwKYb//f/75Z9G6dWuh0+kEAPHQQw/Jj506dUo8/PDDIjo6Wmg0GhEWFia6dOkiXn31VZt22fsZLCgoEM8++6yIjo4Wer1etGnTRqxZs8buz0VZ7bC3DF4I889bixYthFarFcHBwWLo0KE2n5Oy3ivLdg7ViSREqZ3GiKoQg8EgrxrasGGD0s0hhY0ePRrffPMNrl+/rnRTiEhhnAKjKmXs2LHo168fIiMjkZaWhsWLF+Po0aOYP3++0k0jIiIvwgCIqpTs7Gw8++yzuHz5MjQaDdq0aYP169e7lMhLRERVF6fAiIiIqNrhMngiIiKqdhgAERERUbXDAIiIiIiqHSZB22EymXDhwgUEBgba3cqciIiIvI8QAtnZ2YiKirI6ONoeBkB2XLhwATExMUo3g4iIiFxw9uzZCk8PYABkh+WMnLNnz1qdZuwOBoMBGzZsQP/+/cs8+Zzch/3tWexvz2J/exb727Nc6e+srCzExMSUe7i0BQMgOyzTXkFBQZUSAPn5+SEoKIg/QB7A/vYs9rdnsb89i/3tWTfT346krzAJmoiIiKodBkBERERU7TAAIiIiomqHARARERFVOwyAiIiIqNphAERERETVDgMgIiIiqnYYABEREVG1wwCIiIiIqh0GQERERFTtMAAiIiKiaocBkIelZubjeKaE1Mx8pZtCRERUbTEA8qDPfz+Dnu/8hvePqNHznd+wYmeK0k0iIiKqlhgAeUhqZh5e+vYQTML8vUkAL6w6hNTMPGUbRkREVA0xAPKQU+k5EML6mlEInE7PVaZBRERE1RgDIA+JD/WHJFlfU0sS4kL9lGkQERFRNcYAyEMig30xuV9D+XuVBLw+vDkig30VbBUREVH1xADIg+5uFQ0AUEPg1yndkdi+rsItIiIiqp4YAHmQzsfc3SYAEUE6ZRtDRERUjTEA8iDtjQBIQEKRSVRQmoiIiCoLAyAPsgRAAFBYZFKwJURERNUbAyAP0qpLBEBGBkBERERKYQDkQT5qFVQ3lsJzBIiIiEg5DIA8zDINxhEgIiIi5TAA8jDLNFhhEZOgiYiIlMIAyMPkESBOgRERESmGAZCHySNAnAIjIiJSDAMgD+MIEBERkfIYAHkYR4CIiIiUxwDIwzgCREREpDwGQB7GAIiIiEh5DIA8TKs274TIKTAiIiLlMADyMI4AERERKY8BkIcxCZqIiEh5DIA8jCNAREREymMA5GEcASIiIlIeAyAPKx4B4llgRERESmEA5GGcAiMiIlIeAyAP4xQYERGR8hgAeRhHgIiIiJTHAMjDOAJERESkPAZAHsYRICIiIuUxAPIwBkBERETKYwDkYTwLjIiISHkMgDyMI0BERETKYwDkYUyCJiIiUh4DIA/jCBAREZHyGAB5mBwAcQSIiIhIMYoHQAsXLkR8fDz0ej3atm2LLVu2lFv+888/R8uWLeHn54fIyEiMGTMGGRkZ8uNJSUmQJMnmKz8/v7JvxSE8C4yIiEh5igZAK1aswKRJk/Diiy9i79696NatGwYOHIiUlBS75bdu3YoHH3wQY8eOxeHDh/H1119j586dGDdunFW5oKAgpKamWn3p9XpP3FKF5BwgToEREREpRtEAaO7cuRg7dizGjRuHJk2aYN68eYiJicGiRYvslv/9998RFxeHp556CvHx8bj99tvx2GOPYdeuXVblJElCRESE1Ze34BQYERGR8hQLgAoLC7F7927079/f6nr//v2xfft2u8/p0qULzp07h/Xr10MIgYsXL+Kbb77B4MGDrcpdv34dsbGxqFOnDoYMGYK9e/dW2n04yzICZOAIEBERkWJ8lHrh9PR0GI1GhIeHW10PDw9HWlqa3ed06dIFn3/+ORITE5Gfn4+ioiLcddddWLBggVymcePGSEpKQkJCArKysjB//nx07doV+/fvR4MGDezWW1BQgIKCAvn7rKwsAIDBYIDBYLjZW7WigjnwKSgyub1usmXpY/a1Z7C/PYv97Vnsb89ypb+dKSsJIRTJxr1w4QKio6Oxfft2dO7cWb7+2muv4bPPPsOxY8dsnnPkyBH07dsXzzzzDAYMGIDU1FQ899xzaN++PZYsWWL3dUwmE9q0aYPu3bvjvffes1tm+vTpmDFjhs31L774An5+fi7eoX3p+cCsvT7QqQTe7Gh0a91ERETVWW5uLkaOHInMzEwEBQWVW1axEaDQ0FCo1Wqb0Z5Lly7ZjApZzJ49G127dsVzzz0HAGjRogX8/f3RrVs3vPrqq4iMjLR5jkqlQvv27XH8+PEy2zJ16lRMnjxZ/j4rKwsxMTHo379/hR3orLMZ2Zi1dweMUGHQoAFurZtsGQwGJCcno1+/ftBoNEo3p8pjf3sW+9uz2N+e5Up/W2ZwHKFYAKTVatG2bVskJydj2LBh8vXk5GQMHTrU7nNyc3Ph42PdZLVaDQAoayBLCIF9+/YhISGhzLbodDrodDqb6xqNxu0fcj+9+XWKTAJqtQ9UKsmt9ZN9lfFeUtnY357F/vYs9rdnOdPfzrwvigVAADB58mSMGjUK7dq1Q+fOnfHRRx8hJSUF48ePB2AemTl//jyWLVsGALjzzjvxyCOPYNGiRfIU2KRJk9ChQwdERUUBAGbMmIFOnTqhQYMGyMrKwnvvvYd9+/bhgw8+UOw+S7IkQQPmlWB6lVrB1hAREVVPigZAiYmJyMjIwMyZM5GamormzZtj/fr1iI2NBQCkpqZa7Qk0evRoZGdn4/3338eUKVNQo0YN9O7dG3PmzJHLXLt2DY8++ijS0tIQHByM1q1b47fffkOHDh08fn/2WJbBAzcCIA0DICIiIk9TNAACgIkTJ2LixIl2H0tKSrK59uSTT+LJJ58ss753330X7777rrua53ZadfGUFzdDJCIiUobiR2FUN5IkQS2Z85UYABERESmDAZACLLNgDICIiIiUwQBIAT43ZsF4HAYREZEyGAApgCNAREREymIApADLCFABAyAiIiJFMABSAEeAiIiIlMUASAHMASIiIlIWAyAFcASIiIhIWQyAFCCPADEAIiIiUgQDIAX4qG5shGg0KtwSIiKi6okBkAI4BUZERKQsBkAK4BQYERGRshgAKcAyAsR9gIiIiJTBAEgBXAZPRESkLAZACmAOEBERkbIYACmAOUBERETKYgCkAI4AERERKYsBkAKYA0RERKQsBkAKkDdC5AgQERGRIhgAKYBTYERERMpiAKQAyxRYAafAiIiIFMEASAEcASIiIlIWAyAFcBk8ERGRshgAKYAjQERERMpiAKQANZfBExERKYoBkAI4AkRERKQsBkAKYABERESkLAZACvCRbmyEyCkwIiIiRTAAUgBHgIiIiJTFAEgB8kaIDICIiIgUwQBIAcUjQEZlG0JERFRNMQBSAE+DJyIiUhYDIAUwB4iIiEhZDIAUYBkBMgmgiKNAREREHscASAE+JXqd02BERESexwBIAVYBEKfBiIiIPI4BkALUEqDiifBERESKYQCkEO2NYSDuBUREROR5DIAUolWbu545QERERJ7HAEghlhEgToERERF5HgMghcgjQAyAiIiIPI4BkELkESBOgREREXkcAyCFcASIiIhIOQyAFMIcICIiIuUwAFJI8TJ4nghPRETkaQyAFKJVm3dC5D5AREREnscASCGcAiMiIlIOAyCFcCNEIiIi5TAAUghHgIiIiJTDAEghXAZPRESkHAZACuEIEBERkXIYACmEO0ETEREphwGQQjgFRkREpBwGQAop3giRARAREZGnMQBSCJfBExERKYcBkEKYBE1ERKQcBkAKYQBERESkHAZACrGcBeZqAJSamYft/6QjNTPPnc0iIiKqFnyUbkB1dTPL4FfsTMHUVQdhEoBKAmYPT0Bi+7rubiIREVGVxREghbi6DD41M08OfgDAJIAXVh3iSBAREZETFA+AFi5ciPj4eOj1erRt2xZbtmwpt/znn3+Oli1bws/PD5GRkRgzZgwyMjKsyqxcuRJNmzaFTqdD06ZNsXr16sq8BZe4mgN0Kj1HDn4sjELgdHquu5pGRERU5SkaAK1YsQKTJk3Ciy++iL1796Jbt24YOHAgUlJS7JbfunUrHnzwQYwdOxaHDx/G119/jZ07d2LcuHFymR07diAxMRGjRo3C/v37MWrUKIwYMQJ//PGHp27LIZYRoAInp8DiQ/2hkqyvqSUJcaF+7moaERFRladoADR37lyMHTsW48aNQ5MmTTBv3jzExMRg0aJFdsv//vvviIuLw1NPPYX4+HjcfvvteOyxx7Br1y65zLx589CvXz9MnToVjRs3xtSpU9GnTx/MmzfPQ3flGFdHgCKDfTF7eIL8vSQBrw9vjshgX7e2j4iIqCpTLAm6sLAQu3fvxvPPP291vX///ti+fbvd53Tp0gUvvvgi1q9fj4EDB+LSpUv45ptvMHjwYLnMjh078Mwzz1g9b8CAAeUGQAUFBSgoKJC/z8rKAgAYDAYYDAZnb61clvpUMAc+BQaj068xvFUkZnx3BLmFRozvHofhrSLd3s6qwtIv7B/PYH97Fvvbs9jfnuVKfztTVrEAKD09HUajEeHh4VbXw8PDkZaWZvc5Xbp0weeff47ExETk5+ejqKgId911FxYsWCCXSUtLc6pOAJg9ezZmzJhhc33Dhg3w86ucqaV9e3YB8EFm9nWsX7/e6efnG9QAJPxz4iTWF55we/uqmuTkZKWbUK2wvz2L/e1Z7G/Pcqa/c3Mdz4dVfBm8JFkntAghbK5ZHDlyBE899RReeeUVDBgwAKmpqXjuuecwfvx4LFmyxKU6AWDq1KmYPHmy/H1WVhZiYmLQv39/BAUFuXJbZTIYDEhOTsbtnTth3qFd8NHqMWhQD+fqMJpg2vEzAKB2nVgMGtTErW2sSiz93a9fP2g0GqWbU+Wxvz2L/e1Z7G/PcqW/LTM4jlAsAAoNDYVarbYZmbl06ZLNCI7F7Nmz0bVrVzz33HMAgBYtWsDf3x/dunXDq6++isjISERERDhVJwDodDrodDqb6xqNptI+5H46LQDzPkDOvkaesXiIL6fQ+edXR5X5XpIt9rdnsb89i/3tWc70tzPvi2JJ0FqtFm3btrUZ2kpOTkaXLl3sPic3NxcqlXWT1Wo1APMoDwB07tzZps4NGzaUWadSbuYojLxCo/zvrDzORRMRETlL0SmwyZMnY9SoUWjXrh06d+6Mjz76CCkpKRg/fjwA89TU+fPnsWzZMgDAnXfeiUceeQSLFi2Sp8AmTZqEDh06ICoqCgDw9NNPo3v37pgzZw6GDh2Kb7/9Fj///DO2bt2q2H3aczM7QZcMgLLzi9zWJiIioupC0QAoMTERGRkZmDlzJlJTU9G8eXOsX78esbGxAIDU1FSrPYFGjx6N7OxsvP/++5gyZQpq1KiB3r17Y86cOXKZLl26YPny5XjppZfw8ssvo379+lixYgU6duzo8fsrj+UsMINRwGQSUJXe3KcceYYSI0D5HAEiIiJyluJJ0BMnTsTEiRPtPpaUlGRz7cknn8STTz5Zbp333HMP7rnnHnc0r9JYRoAA8yiQXqV2+LklAyCOABERETlP8aMwqivLTtCA89Ng+cwBIiIiuikMgBSiKRkAOZkInVsiALpeWART6cPBiIiIqFwMgBSiUknQ3MgDcjYAKjkFJgSQXcBpMCIiImcwAFKQZRrsZgIgAMhmIjQREZFTGAApyNWl8PmlAqCsPI4AEREROYMBkIJc3QyxZA4QwBEgIiIiZzEAUpDOx7z0vcDZKbBSAVAWl8ITERE5hQGQglwdASo9BcYRICIiIucwAFKQnATtZA5Q6SRo7gVERETkHAZACnJ1BKj0FBh3gyYiInIOAyAFuZwEfWMEyPJ8ngdGRETkHAZACtLJy+CNFZS0ZjkKIzxIB4AjQERERM5iAKSgm90IsXagHgBHgIiIiJzFAEhBLucAGTgCREREdDMYACnIEgC5ug+QPALEVWBEREROYQCkoJtdBh8eZA6AOAJERETkHAZACrrZZfC1A81TYMwBIiIicg4DIAXdfA6QJQmaI0BERETOYACkoJs9CsOSBF1YZLI5HoOIiIjKxgBIQToXcoAMRhMMRgEACA3QQZLM15kHRERE5DgGQApyZQSo5Dlgfjo1AnQ+AJgHRERE5AwGQApyJQCy7AKtksyryIL0GgAcASIiInIGAyAFWZbBFzgxBWYZAfLVqCFJEgL1N0aAuBcQERGRwxgAKUjrowbg2hSYr9Yc+AT5cgSIiIjIWQyAFOTKFFhuoSUAMj83SM8cICIiImcxAFLQzeQA+WrMo0fFOUAMgIiIiBzFAEhBrhyFUTIHCECJHCBOgRERETmKAZCCdDexDF5vGQHy5QgQERGRsxgAKehmcoD8tKVGgJgETURE5DAGQAqSAyAnpsDy5VVgzAEiIiJyFQMgBck5QM5MgRVaT4EF3giAmANERETkOAZACrKMABW4sg+QnAPEZfBERETOYgCkoOIcIMdPcrcEQMU5QNwIkYiIyFk+rj7RYDAgLS0Nubm5CAsLQ61atdzZrmrBpWXwNvsAcQSIiIjIWU6NAF2/fh0ffvghevbsieDgYMTFxaFp06YICwtDbGwsHnnkEezcubOy2lrluLQM3pIDVGoE6HpBEUwm4eYWEhERVU0OB0Dvvvsu4uLi8PHHH6N3795YtWoV9u3bh7/++gs7duzAtGnTUFRUhH79+uGOO+7A8ePHK7PdVYJlCswkgCIHR4HK2ghRCOB6oWvTYKmZedj+TzpSM/Ncej4REdGtxuEpsO3bt2PTpk1ISEiw+3iHDh3w8MMPY9GiRVi6dCk2b96MBg0auK2hVZElAALM02A+6orj0fxSAZBeo4bWR4XCIhOy8gzysnhHrdiZgqmrDsIkAJUEzB6egMT2dZ2qg4iI6FbjcAD09ddfO1ROr9dj4sSJLjeoOtGWCHgKi0zw01b8nOLDUNXytSC9BunXC8xL4Ws6/vqpmXl4ftVBiBszZyYBvLDqELo3DENksK/jFREREd1iXFoFdvHixTIfO3DggMuNqW581CqoJPO/Hc0DKj0FBhQnQju7GeKp9Bw5+LEwCoHT6blO1UNERHSrcSkASkhIwNq1a22uv/322+jYseNNN6o6cXYvoDw7I0CBN84Dc/Y4jPhQf0ilrqklCXGhfk7VQ0REdKtxKQD6z3/+g8TERIwfPx55eXk4f/48evfujbfeegsrVqxwdxurNGeXwpfOAQJcHwGKDPZFy5hg+XuVBLw+vDmnv4iIqMpzaR+gKVOmoG/fvnjggQfQokULXLlyBZ06dcKBAwcQHh7u7jZWaVofNYAi56fASuUAAUBWnvN7AQXoipOmv3ykEzrWC3G6DiIioluNyztB16tXD82aNcPp06eRlZWFESNGMPhxgbN7AeUW2o4ABcojQM4vg0+/XiD/2xJcERERVXUuBUDbtm1DixYtcOLECRw4cACLFi3Ck08+iREjRuDq1avubmOV5uyJ8KVPgweAIDkHyPkRoJIBUMoVJj8TEVH14FIA1Lt3byQmJmLHjh1o0qQJxo0bh7179+LcuXNl7hNE9jlzIrzBaILBaF62ZTUCpHNtBMhoEriSUyh/n5LBAIiIiKoHl3KANmzYgB49elhdq1+/PrZu3YrXXnvNLQ2rLrROTIHll5ii0mtufgToam4hSp6ewREgIiKqLlwaASod/MiVqVR4+eWXb6pB1Y0zy+AtS+BVUnHuEOB6DlDJ6S+AARAREVUfDgdAy5cvd7jSs2fPYtu2bS41qLpxZhl8yU0QJal4Bx9XV4GlZ5unvyzB1NkruRCld0YkIiKqghwOgBYtWoTGjRtjzpw5OHr0qM3jmZmZWL9+PUaOHIm2bdviypUrbm1oVeXMFJi9JfDAzY8ANY8OhiQBOYVGq5wgIiKiqsrhHKDNmzdj3bp1WLBgAV544QX4+/sjPDwcer0eV69eRVpaGsLCwjBmzBgcOnQItWvXrsx2VxlOBUA3psBK5v8ArucAWQKgqBq+SL2WhwuZ+Ui5kouQAJ1T9RAREd1qnEqCHjJkCIYMGYKMjAxs3boVp0+fRl5eHkJDQ9G6dWu0bt0aKpXLWwtVS8UBUMV78FhGgPzKGAFy9iiM9Ovm0Z4Qfy1iavnJAVDruk6cqEpERHQLcmkVWEhICIYOHerutlRLOmdygOxsgggUjwAVFpmQbzDajBCVxTICFBaoQ06BH/44dQVnmQhNRETVgEvDNWfPnsW5c+fk7//8809MmjQJH330kdsaVl24kgNUOsAJ0PrAkhPtTB6QJQAKDdCibi3zAahnuBcQERFVAy4FQCNHjsSmTZsAAGlpaejbty/+/PNPvPDCC5g5c6ZbG1jVuZIDVDoJWqWSEKCzTIM5ngeUcWMKLDRAh7oh5gCIS+GJiKg6cCkAOnToEDp06AAA+Oqrr5CQkIDt27fjiy++QFJSkjvbV+VZlsEXODAFZu8keAvLUnjXRoB0iLkxAsQpMCIiqg5cygEyGAzQ6cwrhX7++WfcddddAIDGjRsjNTXVfa2rBpwZAcotYwQIKJEI7eBeQEIIeQQoJEArT6ulZuWjoMgInY9jeURERES3IpdGgJo1a4bFixdjy5YtSE5Oxh133AEAuHDhAkJCQtzawKrOpX2A3DAClJVXJCdehwboEOKvhZ9WDSGA81fzHKqDiIjoVuVSADRnzhx8+OGH6NmzJ+677z60bNkSALB27Vp5aowc47YAyNe5HKDLN6a/AnU+0N/YWdqSCO1qHlBqZh62/5OO1EwGUERE5N1cmgLr2bMn0tPTkZWVhZo1i/eMefTRR+Hn5+e2xlUHzhyFkV/uFJhlBMixACjDkv8TWLzpYd1afjiWlu1SHtCKnSmYuuogTMJ8Vtns4QlIbF/X6XqIiIg8weVdC9VqNYqKirB161Zs27YNly9fRlxcnNM7QC9cuBDx8fHQ6/Vo27YttmzZUmbZ0aNHQ5Ikm69mzZrJZZKSkuyWyc/Pd/VWK5XODUdhAECQnAPk2BRYyU0QLVwdAUrNzJODHwAwCeCFVYc4EkRERF7LpQAoJycHDz/8MCIjI9G9e3d069YNUVFRGDt2LHJzHf/juWLFCkyaNAkvvvgi9u7di27dumHgwIFISUmxW37+/PlITU2Vv86ePYtatWrh//7v/6zKBQUFWZVLTU2FXq935VYrnUtJ0HamwJwdASq5AszC1aXwp9Jz5ODHwigETqdzRRkREXknlwKgyZMnY/Pmzfjuu+9w7do1XLt2Dd9++y02b96MKVOmOFzP3LlzMXbsWIwbNw5NmjTBvHnzEBMTg0WLFtktHxwcjIiICPlr165duHr1KsaMGWNVTpIkq3IRERGu3KZHyAHQzS6D93XuOAw5AAosHgGKcXEzxPhQf5Q4nB4AoJYkxIVyOpSIiLyTSwHQypUrsWTJEgwcOBBBQUEICgrCoEGD8PHHH+Obb75xqI7CwkLs3r0b/fv3t7rev39/bN++3aE6lixZgr59+yI2Ntbq+vXr1xEbG4s6depgyJAh2Lt3r2M3pgCt2hzMFNzkFJjzI0DFmyBa1C2xF5AQwu7z7IkM9sUdzayDzNeHN0dksK/DdRAREXmSS0nQubm5CA8Pt7leu3Zth6fA0tPTYTQabeoJDw9HWlpahc9PTU3FDz/8gC+++MLqeuPGjZGUlISEhARkZWVh/vz56Nq1K/bv348GDRrYraugoAAFBQXy91lZWQDM+x0ZDM6dsF4RS32W/6phDnwKDMYKXyu3wDy6o5FgU9ZfYx6CuZZb6FCbL2eZ83Nq+vrI5cP9zUdq5BQacTEz1yo/qCISRIl/A70ahri971xRur+pcrG/PYv97Vnsb89ypb+dKetSANS5c2dMmzYNy5Ytk3Nr8vLyMGPGDHTu3NmpuqRScydCCJtr9iQlJaFGjRq4++67ra536tQJnTp1kr/v2rUr2rRpgwULFuC9996zW9fs2bMxY8YMm+sbNmyotFVtycnJAICj1yQAaqRfuYb169eX+5yLGWoAEg7u2w3DaesRmmM36rlw+WqF9QDAiXPmus78dQjr0w/K14M1alwrlLBi3c+IC3T8fv48Ya5PLQkYhYQFX/+CdmGOjyJVNkt/k2ewvz2L/e1Z7G/Pcqa/nclDdikAmjdvHgYOHIg6deqgZcuWkCQJ+/btg06nw4YNGxyqIzQ0FGq12ma059KlS3ZHl0oSQmDp0qUYNWoUtNryRylUKhXat2+P48ePl1lm6tSpmDx5svx9VlYWYmJi0L9/fwQFBTlwN44zGAxITk5Gv379oNFoEHLqChYf3QW9fwAGDepa7nPn/rUVyM1Fj66d0Da2ptVj0ecysejoH4DGF4MGda+wHW8d2wIgD/27W9f1v9Sd2Hn6KmKatMagFpEO3VN2vgHpO8xnw93TNgYrdp3DVd9oDBrUwqHnV6bS/U2Vi/3tWexvz2J/e5Yr/W2ZwXGESwFQQkICjh8/jv/97384duwYhBC49957cf/998PX17G8D61Wi7Zt2yI5ORnDhg2TrycnJ2Po0KHlPnfz5s04ceIExo4dW+HrCCGwb98+JCQklFlGp9PJR3uUpNFoKu1DbqnbV2cO4AqNpgpfK7/InAMU4KuzKVszwDwSl11Q5FCbLcdgRNTwtyofG+KPnaev4kJmgcP3/leK+QMXXcMX97QzB0BbTmRAUqnho3Z5pwW3qsz3kmyxvz2L/e1Z7G/Pcqa/nXlfXAqAZs+ejfDwcDzyyCNW15cuXYrLly/jP//5j0P1TJ48GaNGjUK7du3QuXNnfPTRR0hJScH48eMBmEdmzp8/j2XLllk9b8mSJejYsSOaN29uU+eMGTPQqVMnNGjQAFlZWXjvvfewb98+fPDBB67caqVzah+gcjZCtByFcb2gCCaTgEpV9jRibmGRnFBdciNEwLW9gA5fyAQAJEQHo3VMDQT7apCZZ8CelGvoEF/L4XosUjPzcCo9B/Gh/kykJiKiSuFSAPThhx/aJB8D5jPC7r33XocDoMTERGRkZGDmzJlITU1F8+bNsX79enlVV2pqqs2eQJmZmVi5ciXmz59vt85r167h0UcfRVpaGoKDg9G6dWv89ttvXntEhzP7AOUbzGXs7wNkfiuFAK4XFskBkT3p2ebRH71GBf9SwVSsC3sBHTxvDoCaRwfBR61Cz0Zh+HbfBWw8dsnpAIg7ShMRkSe4FAClpaUhMtI2PyQsLMzp0+AnTpyIiRMn2n0sKSnJ5lpwcHC5SU7vvvsu3n33XafaoCT5KIwKAqAio0neK8jPzgiQXqOG1keFwiITsvIM5QZAlnPAQvx1NgnnMfJSeMd3cT4kB0DBAIDejWvj230XsOnYJTw/sLHD9ZS1o3T3hmEcCSIiIrdyKUEjJiYG27Zts7m+bds2REVF3XSjqhNHN0K0TFkB5mDHHstxGBWdCJ9u5xwwC8sU2IXMPBQUGW0eLy2noAgn03MAAM2izAFQj4ZhUEnAXxezcf6a44EUd5QmIiJPcSkAGjduHCZNmoRPP/0UZ86cwZkzZ7B06VI888wzNnlBVD5LAGQwCphK//UvwRIASVJx3lBpllGfrLzy90GwJECHBdiuoAvx18JPq4YQwPmrFQcvR1KzIAQQEaRH2I2AqoafFm3qmleWbTx2qcI6LOJD/VE6dUmSwB2liYjI7VyaAvv3v/+NK1euYOLEiSgsvJFPotfjP//5D6ZOnerWBlZ12hLBTKHRBL3K/uhOfmFx/k9Z+yQFOjsCFGA7AiRJknwqfMqVXNQLCyi3rtLTXxa9GtfGrjNXsenYJYzqFGvvqTYig33RMb4Wdpy8Il/z06ih97HfJ45IzczH8UwJqZn5qBvKVRtERGTm0giQJEmYM2cOLl++jN9//x379+/HlStX8Morr7i7fVWeVm0dAJUl12AOauzl/1gE+d4YAargOAxLABRiZwQIKJkHVPHUU8kE6JJ6N64NANj+T7p8hllFMvMM2H/OXN9/BjRCfKgfcgqNmLXuiEPPL23FzhT0fOc3vH9EjZ7v/IYVO+0fsktERNXPTW3SEhAQgPbt26N58+Z299GhilkFQOUkQluWwJeV/wO4ZwQIcG4p/OHz5j2AEkqNADWOCERksB75BhN2/JNRYT0A8PWus8gtNKJheADG96yPd0a0giQBq/aex69/OT6VBpSdUJ2a6XhOUun6tv+T7vLziYjIu3jHLnXVmEolQaM2T2mVGwCVcxK8haM5QJZl8DcbAOUVGnH8UjYA2ykwSZLQ68YokCN5QEaTQNL20wCAMV3jIUkS2tStiTFd4gEAL64+hOsFjp10DwB7zlx1W0L1ip0p6PrGRoz8+A90fWPjTY8kMZgiIlIeAyAv4MhS+PxyToK3kEeAKggU0nMcHQEq/w/00bQsmIS5ntp2VpT1blQcAFV0uvzPRy/i3NU81PDT4O5W0fL1Zwc0RJ2avjh/LQ9v//RXuXVYpF8vwJs/HrP72Fe7UuTRNEfYG0mauuqgy8GLO4MpRwIpd5VxvK7inKtbq93eVcbxutjfnm2T5/rbve2+Nfu7srmUBE3upfVRIafQWG4OUN6NJOjypsAcHwEyB0BhgfZzgOqGFOcAlXc47eHzlh2gg+yW6XJbCLQ+Kpy/loe/L15Ho4iyT1dduvUUAOC+DnWtgjw/rQ9mD0/AqCV/4r87TuPOlpFoG1v25orXcgsxasmfOHMlD8G+PsjOL4JJmE+oFwBW772A/WczMe/eVmhRp0aZ9QDA8YvZmLnuiM1IkkkA9yzajuFt6mBAswg0izLff1k7WBcWmXA0NQub/76Mucl/W9Xz/MqDEAJoE1sTdWv5ye+vvbpMJoEruYW4lFWAr3efRdK20xAw39vY2+MxtFU0gnx9EOyrQaBeg292n7W7qaTRJMy7gRca8dXuc3hnw18QwrziblLfBri7VTR0PuZ9pXQ3vlbuOWdV12vDEjCsdTSKTELeo2r1nvOY8+MxmIQaC4/+hucHNsGw1tFQqyT5a83e83jl20NyPTOHNsfwNtEwmgRMwnx0zcrd5/Da+qNymRcHNcGwNnWgkgAJEiQVsHrPOcz47ohcZtqdzTC8TTQEzJuBCiEgBLB673m8+n1xuZcGN8XwNtGQIJk7DsDqvecws0Rdr9zZFMPb1IEcswtg1d5zmLWuuMzLQ5pieOs6ch0O1QNg1R479bSpY/PZc6RccRlzf5dfxpF6bq6Mp19PuTZVfn97Y19WZruV2vBWEhX9r3k1lJWVheDgYGRmZlbKYajr16/HoEGD5DNLOr7+My5mFWDdk7fbTCVZfL3rLJ775gB6NgpD0hj7u1onbTuF6d8dweCESHxwfxu7ZQqKjGj00o8AgL0v90NNf9sgKN9gRJNXfoQQwJ6X+6GWnTIA8O9v9uOrXefwZO/bMKV/I7tlHlr6Jzb/fRn/uaMxJvSsb7fM4QuZGPzeVqhVErb8uxeiathuevjs1/vxze5zuK12AL5/6nbo7KwMy8434IElf2L/2WsIDdDhq8c6wUcS+Gr9JowY1AspVwsw5et9uJhVAB+VhEl9G2BCz9ugLrX2ft/Za1i46QQ2HLlot72lRdfwRXyoP7b9kw5x4wd6RLsY+GrV2Hf2Gg5fyHJop29JAqKCfeGrVeHEpRz5ep2avjCaBC5nF6ConK0SHKFVSyg08keeiLyHWpKw9fleNhve2vt7WRFn/n5zCswLWJbCFzgyBVbeCJADq8AsewD5qCQE+9r/QOk1akQEmQ9XPZORY7cMABy6kQBdVtAGAH2amKfBNpWTB5S07TQA4I7mEXaDHwB4aXAThAbocOLSdXyw8YTN47mFRRibtAv7z15DTT8NPh/XEfXCAhAZrEeDYIHIYD1ubxCKnyZ1x6CECBSZBN7e8DcSP9yB3WeuYPuJdKzdfx73f/I77v5gmxz8DGgWjid63Qb1jREutSRh2p1NMS+xFe5oFgG9xjzCtfVEuvx/+iYBLN95Fp9uO429KddQWGRCDT8NOtWrBalUuyWYE8YD9T7mvZeu5VkFPwBw7moeUjPz5eDHsuFlabX8NeWuEgTgUPCj91HJeWlERJVNqQ1vOQXmBRzJAXIkCTrQMgVWziowSwAUEqAt98DUmFp+SM3MR8qVXLS+salhSfkGI/6+aD8BuqRejWoDOIzdKVeRmWtAsJ910JVxvQDf7r8AAHi4a3yZ9dTw02Lm0GaY+PkeLPz1HwxMiESTyCC5LY99tht/nr6CQL0PPhvbsczpthp+Wnwwsg1W7TmPaWsPY9eZq/jXoh1WZdQqCUNbRWFCj/poEG6u5/5OdXE6PRdxoX7y/6Xc3ToaeYVGfLzlpNXUlsWAZuEY2DwSrWJqIDbED5IkYcXOFLyw6hCMQkAtSXh9eHMktq8LIQSu5BTiuwMXMH2t7bL/V+9ujj5NaiM0QIf06wXo+sZGq6k5tSTh+6e6ITLYF4VFJhy/lI07F2y1KqOSgK/Hd0Z8aAD8tGpcySnA7XM22dSz6bmeiAw2jzoVFplw9kou7pj/m01dP03qjphafvBRSbiUnW+3ri3/6YnwIHNd56/moc/cX23q+WVyT0TW0EMlmevp/mbpeoDf/t0L4UF6CAAXruWh19u29Wyc0hMRweZ6JAlIy8xDj7dsy21+ricign0hhLlMTzt1/fqsuYwkmacje9qp59cb/QQAqdfcUw/g4OvdTJkbbQIqvn9Hy7izrqrcJrbbfhm1JCmy4S1HgLyA9sZ0jkM5QOXtA2RJgi4nB6iiJfAWdSvYC+jvi9koMgnU9NMgKlhfZj0xtfzQoHYAjCaBzccv2zz+xR8pKCwyoWWdYLSpW6PcNg1sHoH+TcNRZBJ4fuUBGE0CBqMJT3yxB1uOp8NPq0bSmA7lBmSAeYXav9rWwX8fbm/7GICvH+uEuSNaycEPYN6ksXP9EJshWl+tGv/Xro7NDtZqScL0u5rh7tbRiAv1l3OkEtvXxdbne+HLRzph6/O95HlvSZIQEqDDgGYRduvq06Q2IoN9oVGrEBnsi9nDE6xGpV4f3lxum9ZHhWZRwTZlZg9PQNvYWqjlr4Veo0ZUDb9y61GrJPhq1WgYEWi3rgbhgdBr1PBRq+S6LG1XScDrw5sjqoYf1CoJWh8V4sP87dYTH+Yvn2VXp6a9NiUguqYffNQqaNQqxIbYrycutLgejVqFmFr2y8XU8odGrYLWR4W6ZdRVN8RfrqduGfXUvVGPRu2+ehx+vRtlSvZ3WWXKapMj9+9oGXfW5c1tsunvW6Td3t7fJX/veBJzgOzwdA7Q0A+2Yf/Za/j4wXbo1zTc7vNm/3AUH24+iXG3x+OlIU3tljlyIQuD3tuC0AAddr3U126Zr3adxb+/OYDuDcOw7GH7uUQA8N4vxzE3+W+MaFcHb97T0ubxL/5IwQurD6Jbg1B8NrZjufc8e/1RfPjbSQxrHY13E1vJ1wuLTOg6ZyMuZxdgXmIr3N06uuxKbriYlY++czcjO78Ij/esj90pV/H7ySvQ+ajw6Zj26FI/1Kp8eXPI2/9Jx8iP/7B5jS8f6YTO9UMqbEtJZY3suMLRulIz82xGpSqrjKPlUtKz5ZyruqH2R+G8sd3eVsbRcuxvz7bJk/3tzrpu1f6u7BwgToF5AZ0jy+ALK14GH+RrfjvLywEqHgGyn9hsUdFeQAfLOALDnl6Na+PD307i178uwWgSctLx+oOpuJxdgNqBOgxKiKywHgAID9LjxUFN8Pyqg/jg13/k6yM71rUJfipiOXvMHUOxie3ronvDMId+Obirrshg3wpfx11lHK+rOOeqstvk3nZ7VxnH62J/e7ZNnutvd9Z1q/Z3ZeMUmBcoPhG+7P1pLDlA5e8EbY6QC4tMZR4/YdkEMayCKbDi4zDs79Fw+MKNACiq4gCobWxNBOp9cDXXgH1nrwEwL1P+dJt56fsDnWKtzkSrSPeGtoHOsu1nnN5PoqKpJGeVNU2mdF1ERGSLI0BeQA6Ayk2CLj4MtSyBOh9IknkflOz8IrvBUkYFmyBaWEaALmTmobDIZH1oa5EJx1LNCdClj8CwR6NWoXvDMHx/IBWbjl1C29ia2JNyDfvPZULro8LIjs5NFZ3OsB2VsqwicDZgcOfIDRER3To4AuQFHFoF5sAUmEolIUBrOQ/M/jSYPAVWxiaIFqEBWvhp1fLS7JKOX8pGodGEIL0PYmo5FjBYdoX+5cZyeMvoz9CWURUGY6VZpq5KuplVBBxtISKqfhgAeQFH9gHKc+A0eKDkXkD2l8JbpsBC/MsPOiRJKjMP6HCJ/X/K2iW6tJ6NwiBJwNHULOxJuYofDqUBMJ/75Sx3T10REVH1wykwL1CcA3Rzp8EDJU+Er2AEyIFRl5hafjiWlo2UjBwAYfJ1ZxKgLUICdGhZpwb2nb2Gp5fvhdEk0DG+FppGubbKjlNXRER0MxgAeQF35QABJc8Dsx0BMt44SwqoeAoMKHsl2KEbCdDNnAxe+jSujX1nr8mJ1a6M/pTkDasIiIjo1sQpMC/grtPggfJHgK7kFMqHXtbycy0AKjKaD/YEHEuALslQaoTr6o1gjIiIyNMYAHkBnQMjQLmF5hGdCkeAyjkPzDL9VctPCx91xW99cQBUnAT9z+Uc5BtMCND5IC7Ev8I6LFIz8/D+JuszvF5afcjppetERETuwADICziTA1TRCJB8HIadJGhLABRSwSaIFjEljsOwbBhuyf9pGhVU7llipZ1Kz7HacBBQ7gA8IiIiBkBewLEpMMdygOQDUe2cB+ZMAjQA1Klpzq+5XlCEq7nm+g6dd3wDxJLcvXSdiIjoZjAA8gIVJUEXGU3y6FDFU2BljwBZToJ3NADSa9SICDJv+W7JA5J3gI52LgGaS9eJiMibcBWYF5D3ASpjCiy/RGBUcRJ02TlAl50cAQLMeUBpWflIuZKLhOhgHL7gWgI0wKXrRETkPRgAeYGKRoAsCdCSVJwwXRZ5Gby9HCDLJogO5gAB5jygP09fwdkruTiVnoPcQiN8NWrUCwtwuI6SuHSdiIi8AafAvIDOxzyqU1YAlF9YPP1V0c7LlmXw5eUAVXQQakmxIeYcnTMZOXL+T9OoIPlEdyIiolsRR4C8QEUjQJaT4CvK/wGKl8HbzQHKcewcsJJK7gUUpLckQLu2ezMREZG3YADkBeRVYGXkAFkCoIqOwQBKjADZ2wco27kkaKDkUvji/XqauZD/Q0RE5E0YAHmBijZCtOQAVXQQKlCcA3S9oAgmk5D36hFCFI8AOZkEDQAXMvOQeWNazZUEaCIiIm/CHCAvUNEUmKPHYADFI0BCANcLi6fBMvMMMBjNOxHW8nd8Ciw0QAtfjdpcX0ERtD4q3FbbtQRoIiIib8EAyAtUtBN03o0kaEemwPQatVxfyTwgSwJ0oN7HoXosJEmSR4EAoElkEDQOHKNBRETkzfiXzAtUtBO0M0nQQPFxGCVXgqXf2ATRmRVgFjElAiAmQBMRUVXAAMgLyBshui0Asl0J5uwxGCWVHAGKqcmjK4iI6NbHAMgLFOcAGe0+nudEEjRgfy+g9GznDkItKf16vvzvOT8dw4qdKU7XQURE5E0YAHmBCpfBW3KAHAyA5L2ACmynwJwdAUrNzMN3B1Ll74UAXlh1CKmZeeU8i4iIyLsxAPICFS2Dd3YKrHgEqHgKzJUl8ABwKj0HQlhfMwqB0+m5TtVDRETkTRgAeQHLFJhJmE9+Ly3fxRygklNgly2bIDqxCzQAxIf6o/SpF2pJQlwoc4GIiOjWxQDIC2hLHHBqbxrMshGiI/sAAcUjQNkFN58EHRnsi9nDE6C+cQaZWpLw+vDmPNCUiIhuadwJ2gtoS+yrU1hkgl+pQZo8Q/FhqI6wNwJUHAA5nwSd2L4uujcMw+n0XMSF+jH4ISKiWx4DIC/go1ZBJZmnwOzlAeUVOr4TNFBiBOjGMnghxE0tgwfMI0EMfIiIqKrgFJiXKG8vIKdzgG6sArMciJpbaET+jVEkVwMgIiKiqoQBkJcobym8M6fBA0CgZQrsxgiQZfTHV6OGv46DfkRERAyAvITWxxzc2JsCy70xBeboRoiWozCyb+QAydNfTq4AIyIiqqoYAHmJ8vYCcuY0eMB2BMiyBD7En9NfREREAAMgr1HeifByErTDOUA3NkK8kQPk6iaIREREVRUDIC9R3onwruYAFRaZkG8wIv3GCFAYp8CIiIgAMADyGtpypsDynJ0C0/ngxr6FyM4vuukl8ERERFUNAyAvUdYyeKNJyEGRn4MjQCqVhACtZS8ggxwAhfhzBIiIiAhgAOQ1yloGbxn9ARwfAQJK7gVUhAzLSfCBHAEiIiICGAB5jbKmwCwJ0EDxSjFHFO8GbeAUGBERUSkMgLxEWQFQyV2gJUmyeV5Zis8DK8JlBkBERERWGAB5ieIAyGh13dlNEC0sI0Dp1wvkM8HCGAAREREBYADkNXQV5AA5ugTewpIDdCo9BwCgUUvy/kBERETVHQMgL1FRDpAzCdBA8QiQJQAK8dc5NYVGRERUlTEA8hKO5AA5w5IDZAmAeA4YERFRMQZAXsKyDL6gjCkwZwMgywjQuau5AJgATUREVBIDIC9R1ghQrotTYJYcIJMwf8+DUImIiIoxAPISZeYA3eQIkAWnwIiIiIopHgAtXLgQ8fHx0Ov1aNu2LbZs2VJm2dGjR0OSJJuvZs2aWZVbuXIlmjZtCp1Oh6ZNm2L16tWVfRs3rcwcIFdHgG7kAFlwCTwREVExRQOgFStWYNKkSXjxxRexd+9edOvWDQMHDkRKSord8vPnz0dqaqr8dfbsWdSqVQv/93//J5fZsWMHEhMTMWrUKOzfvx+jRo3CiBEj8Mcff3jqtlxS0VEYzi6DtxkBYgBEREQkUzQAmjt3LsaOHYtx48ahSZMmmDdvHmJiYrBo0SK75YODgxERESF/7dq1C1evXsWYMWPkMvPmzUO/fv0wdepUNG7cGFOnTkWfPn0wb948D92Va3QV5AA5uxGiJQfIggEQERFRMcV2xissLMTu3bvx/PPPW13v378/tm/f7lAdS5YsQd++fREbGytf27FjB5555hmrcgMGDCg3ACooKEBBQYH8fVZWFgDAYDDAYDA41BZHWeorXa9aMmcr5xuKrB7LKTD/W6u2fU55fEvFS8F6ldvv5VZQVn9T5WB/exb727PY357lSn87U1axACg9PR1GoxHh4eFW18PDw5GWllbh81NTU/HDDz/giy++sLqelpbmdJ2zZ8/GjBkzbK5v2LABfn5+FbbFFcnJyVbfH7ksAVDjQtolrF+/Xr5+/KQKgAopJ09g/frjDtdvMAEl3959v2/ByWqcB126v6lysb89i/3tWexvz3Kmv3Nzcx0uq/jZCKV3JxZCOLRjcVJSEmrUqIG77777puucOnUqJk+eLH+flZWFmJgY9O/fH0FBQRW2xRkGgwHJycno168fNJriaSrpUBr+d+IAgmqGYNCg9vL15K8OAJfS0LJ5UwzqEmuvyjJN3ZUMg1FAJQH/d9dAqFXVbyfosvqbKgf727PY357F/vYsV/rbMoPjCMUCoNDQUKjVapuRmUuXLtmM4JQmhMDSpUsxatQoaLXWwxoRERFO16nT6aDT2ebIaDSaSvuQl67bV2e+D4NRWF0vMJqnxgL0WqfbEuyrQfr1QtTy10Kvq8bDP6jc95Jssb89i/3tWexvz3Kmv515XxRLgtZqtWjbtq3N0FZycjK6dOlS7nM3b96MEydOYOzYsTaPde7c2abODRs2VFin0io6C8zZJGgACLyxFJ4J0ERERNYUnQKbPHkyRo0ahXbt2qFz58746KOPkJKSgvHjxwMwT02dP38ey5Yts3rekiVL0LFjRzRv3tymzqeffhrdu3fHnDlzMHToUHz77bf4+eefsXXrVo/ck6vcvQweAPQ3gqoAvfPPJSIiqsoUDYASExORkZGBmTNnIjU1Fc2bN8f69evlVV2pqak2ewJlZmZi5cqVmD9/vt06u3TpguXLl+Oll17Cyy+/jPr162PFihXo2LFjpd/PzXD3afArdqbgaFo2AGDX6WtYsTMFie3ruqGlREREtz7Fk6AnTpyIiRMn2n0sKSnJ5lpwcHCFWd733HMP7rnnHnc0z2PK2gfIldPgUzPzMHXVQatrL6w6hO4NwxAZ7HuTLSUiIrr1KX4UBpnJI0DGm98I8VR6jnwIqoVRCJxOd3x5IBERUVXGAMhLyDlAZRyG6kwOUHyoP0qveFdLEuJCK2dPIyIiolsNAyAvUeFp8E6MAEUG+2L28ASob+x9pJYkvD68Oae/iIiIblA8B4jMSk6BWTZuNJqEHBA5kwMEAInt66J7wzCcTs9FXKgfgx8iIqISGAB5CUsABJiDIJ2PWk6ABpwPgADzSBADHyIiIlucAvMSlhwgoHgazJIADQB6Dd8qIiIid+FfVS9hLwAquQTekfPRiIiIyDEMgLyESiVBozYHOZal8K4kQBMREVHFGAB5kdJL4eVdoF3I/yEiIqKyMQDyIqWXwltygJj/Q0RE5F78y+pFLAFQQakcID8tF+sRERG5EwMgL1L6OIw8F84BIyIioooxAPIiZeUA6ZkETURE5FYMgLyI1scc6MgBkDwCxLeJiIjInfiX1YuUToLOK2QOEBERUWVgAORFdGr7OUDOnARPREREFWMA5EVsRoCYBE1ERFQpGAB5kbKmwHy1fJuIiIjciX9ZvYhlFVhBkTnw4U7QRERElYMBkBcpvRFi8VlgTIImIiJyJwZAXoQbIRIREXkGAyAvUjoHKN/AHCAiIqLKwL+sXoSnwRMREXkGAyAvoivjNHjmABEREbkXAyAvUjoHKJ85QERERJWCAZAXsZkCYwBERERUKRgAeZEyd4JmEjQREZFb8S+rF5H3ATJa5wDxLDAiIiL3YgDkRUqOABlNQh4J4mnwRERE7sUAyIuUzAGyJEADzAEiIiJyNwZAXqTkCFBeiQDIsjyeiIiI3IN/Wb2IrsQy+Dw5/0cFlUpSsllERERVDgMgL2JvBIj5P0RERO7HAMiLaNXmXJ/CIhOPwSAiIqpEDIC8SMmdoC0jQHoN3yIiIiJ3419XL2JvCsxXyxEgIiIid2MA5EUsy+ALOAVGRERUqRgAeZHiESBjcQDEJGgiIiK3YwDkRXR2coB8mQNERETkdvzr6kVK5gDl8yR4IiKiSsMAyItYcoBMAsjOLwLAJGgiIqLKwADIi2hLHHmRmWcAAPhqmANERETkbgyAvIjdAEjLt4iIiMjd+NfVi/ioJEg3jv0qHgHiFBgREZG7MQDyIpIkyXlAlgBIzwCIiIjI7RgAeRnLNFjxFBgDICIiIndjAORldD7mgOdarjkA8mMARERE5HYMgLyMZTPELOYAERERVRoGQF6m5InwAHOAiIiIKgMDIC9jSYK24AgQERGR+zEA8jIl9wICmARNRERUGRgAeZnSARCToImIiNyPAZCXKT0FxhwgIiIi92MA5GVspsAYABEREbkdAyAvwxwgIiKiyscAyMuUDoD0PgyAiIiI3I0BkJfRlcgB0mtUUKkkBVtDRERUNTEA8jIlR4CY/0NERFQ5GAB5GQZARERElU/xAGjhwoWIj4+HXq9H27ZtsWXLlnLLFxQU4MUXX0RsbCx0Oh3q16+PpUuXyo8nJSVBkiSbr/z8/Mq+FbcouQxezwRoIiKiSuGj5IuvWLECkyZNwsKFC9G1a1d8+OGHGDhwII4cOYK6devafc6IESNw8eJFLFmyBLfddhsuXbqEoqIiqzJBQUH466+/rK7p9fpKuw93KjkCxE0QiYiIKoeiAdDcuXMxduxYjBs3DgAwb948/PTTT1i0aBFmz55tU/7HH3/E5s2bcfLkSdSqVQsAEBcXZ1NOkiRERERUatsrC6fAiIiIKp9iAVBhYSF2796N559/3up6//79sX37drvPWbt2Ldq1a4c333wTn332Gfz9/XHXXXdh1qxZ8PX1lctdv34dsbGxMBqNaNWqFWbNmoXWrVuX2ZaCggIUFBTI32dlZQEADAYDDAbDzdymDUt9ZdXrU2LRl85H5fbXr24q6m9yL/a3Z7G/PYv97Vmu9LczZRULgNLT02E0GhEeHm51PTw8HGlpaXafc/LkSWzduhV6vR6rV69Geno6Jk6ciCtXrsh5QI0bN0ZSUhISEhKQlZWF+fPno2vXrti/fz8aNGhgt97Zs2djxowZNtc3bNgAPz+/m7xT+5KTk+1eP3FBAmAe+cnMuIz169dXyutXN2X1N1UO9rdnsb89i/3tWc70d25ursNlFZ0CA8zTVSUJIWyuWZhMJkiShM8//xzBwcEAzNNo99xzDz744AP4+vqiU6dO6NSpk/ycrl27ok2bNliwYAHee+89u/VOnToVkydPlr/PyspCTEwM+vfvj6CgoJu9RSsGgwHJycno168fNBqNzeNX/0jBmjPHAADxMdEYNCjBra9f3VTU3+Re7G/PYn97Fvvbs1zpb8sMjiMUC4BCQ0OhVqttRnsuXbpkMypkERkZiejoaDn4AYAmTZpACIFz587ZHeFRqVRo3749jh8/XmZbdDoddDqdzXWNRlNpH/Ky6vbVFV/z1/vwh8xNKvO9JFvsb89if3sW+9uznOlvZ94XxZbBa7VatG3b1mZoKzk5GV26dLH7nK5du+LChQu4fv26fO3vv/+GSqVCnTp17D5HCIF9+/YhMjLSfY2vRCWToHkSPBERUeVQdB+gyZMn45NPPsHSpUtx9OhRPPPMM0hJScH48eMBmKemHnzwQbn8yJEjERISgjFjxuDIkSP47bff8Nxzz+Hhhx+Wk6BnzJiBn376CSdPnsS+ffswduxY7Nu3T67T22nVxUEPV4ERERFVDkVzgBITE5GRkYGZM2ciNTUVzZs3x/r16xEbGwsASE1NRUpKilw+ICAAycnJePLJJ9GuXTuEhIRgxIgRePXVV+Uy165dw6OPPoq0tDQEBwejdevW+O2339ChQweP358ruAyeiIio8imeBD1x4kRMnDjR7mNJSUk21xo3blxuRvi7776Ld999113N8zirAIgbIRIREVUKxY/CIGslj8JgAERERFQ5GAB5GU6BERERVT4GQF5GxwCIiIio0jEA8jLMASIiIqp8DIC8jFUOEEeAiIiIKgUDIC/DESAiIqLKxwDIy5QMgK7nFynYEiIioqqLAZCXWXfggvzvB5b8gRU7U8opTURERK5gAORFUjPzMPO7I/L3JgG8sOoQUjPzFGwVERFR1cMAyIucSs+BSVhfMwqB0+m5yjSIiIioimIA5EXiQ/2hkqyvqSUJcaF+yjSIiIioimIA5EUig30xe3gC1JI5ClJLEl4f3hyRwb4Kt4yIiKhqUfwwVLKW2L4uujcMw+n0XMSF+jH4ISIiqgQMgLxQZLAvAx8iIqJKxCkwIiIiqnYYABEREVG1wwCIiIiIqh0GQERERFTtMAAiIiKiaocBEBEREVU7DICIiIio2mEARERERNUOAyAiIiKqdhgAERERUbXDAIiIiIiqHZ4FZocQAgCQlZXl9roNBgNyc3ORlZUFjUbj9vrJGvvbs9jfnsX+9iz2t2e50t+Wv9uWv+PlYQBkR3Z2NgAgJiZG4ZYQERGRs7KzsxEcHFxuGUk4EiZVMyaTCRcuXEBgYCAkSXJr3VlZWYiJicHZs2cRFBTk1rrJFvvbs9jfnsX+9iz2t2e50t9CCGRnZyMqKgoqVflZPhwBskOlUqFOnTqV+hpBQUH8AfIg9rdnsb89i/3tWexvz3K2vysa+bFgEjQRERFVOwyAiIiIqNphAORhOp0O06ZNg06nU7op1QL727PY357F/vYs9rdnVXZ/MwmaiIiIqh2OABEREVG1wwCIiIiIqh0GQERERFTtMAAiIiKiaocBkActXLgQ8fHx0Ov1aNu2LbZs2aJ0k6qM3377DXfeeSeioqIgSRLWrFlj9bgQAtOnT0dUVBR8fX3Rs2dPHD58WJnG3uJmz56N9u3bIzAwELVr18bdd9+Nv/76y6oM+9t9Fi1ahBYtWsibwXXu3Bk//PCD/Dj7unLNnj0bkiRh0qRJ8jX2uftMnz4dkiRZfUVERMiPV2ZfMwDykBUrVmDSpEl48cUXsXfvXnTr1g0DBw5ESkqK0k2rEnJyctCyZUu8//77dh9/8803MXfuXLz//vvYuXMnIiIi0K9fP/ncN3Lc5s2b8fjjj+P3339HcnIyioqK0L9/f+Tk5Mhl2N/uU6dOHbzxxhvYtWsXdu3ahd69e2Po0KHyHwH2deXZuXMnPvroI7Ro0cLqOvvcvZo1a4bU1FT56+DBg/JjldrXgjyiQ4cOYvz48VbXGjduLJ5//nmFWlR1ARCrV6+WvzeZTCIiIkK88cYb8rX8/HwRHBwsFi9erEALq5ZLly4JAGLz5s1CCPa3J9SsWVN88skn7OtKlJ2dLRo0aCCSk5NFjx49xNNPPy2E4Ofb3aZNmyZatmxp97HK7muOAHlAYWEhdu/ejf79+1td79+/P7Zv365Qq6qPU6dOIS0tzar/dTodevTowf53g8zMTABArVq1ALC/K5PRaMTy5cuRk5ODzp07s68r0eOPP47Bgwejb9++VtfZ5+53/PhxREVFIT4+Hvfeey9OnjwJoPL7moehekB6ejqMRiPCw8OtroeHhyMtLU2hVlUflj621/9nzpxRoklVhhACkydPxu23347mzZsDYH9XhoMHD6Jz587Iz89HQEAAVq9ejaZNm8p/BNjX7rV8+XLs2bMHO3futHmMn2/36tixI5YtW4aGDRvi4sWLePXVV9GlSxccPny40vuaAZAHSZJk9b0QwuYaVR72v/s98cQTOHDgALZu3WrzGPvbfRo1aoR9+/bh2rVrWLlyJR566CFs3rxZfpx97T5nz57F008/jQ0bNkCv15dZjn3uHgMHDpT/nZCQgM6dO6N+/fr473//i06dOgGovL7mFJgHhIaGQq1W24z2XLp0ySayJfezrChg/7vXk08+ibVr12LTpk2oU6eOfJ397X5arRa33XYb2rVrh9mzZ6Nly5aYP38++7oS7N69G5cuXULbtm3h4+MDHx8fbN68Ge+99x58fHzkfmWfVw5/f38kJCTg+PHjlf75ZgDkAVqtFm3btkVycrLV9eTkZHTp0kWhVlUf8fHxiIiIsOr/wsJCbN68mf3vAiEEnnjiCaxatQobN25EfHy81ePs78onhEBBQQH7uhL06dMHBw8exL59++Svdu3a4f7778e+fftQr1499nklKigowNGjRxEZGVn5n++bTqMmhyxfvlxoNBqxZMkSceTIETFp0iTh7+8vTp8+rXTTqoTs7Gyxd+9esXfvXgFAzJ07V+zdu1ecOXNGCCHEG2+8IYKDg8WqVavEwYMHxX333SciIyNFVlaWwi2/9UyYMEEEBweLX3/9VaSmpspfubm5chn2t/tMnTpV/Pbbb+LUqVPiwIED4oUXXhAqlUps2LBBCMG+9oSSq8CEYJ+705QpU8Svv/4qTp48KX7//XcxZMgQERgYKP9trMy+ZgDkQR988IGIjY0VWq1WtGnTRl42TDdv06ZNAoDN10MPPSSEMC+nnDZtmoiIiBA6nU50795dHDx4UNlG36Ls9TMA8emnn8pl2N/u8/DDD8u/N8LCwkSfPn3k4EcI9rUnlA6A2Ofuk5iYKCIjI4VGoxFRUVFi+PDh4vDhw/LjldnXkhBC3Pw4EhEREdGtgzlAREREVO0wACIiIqJqhwEQERERVTsMgIiIiKjaYQBERERE1Q4DICIiIqp2GAARERFRtcMAiIjIjri4OMybN0/pZhBRJWEARESKGz16NO6++24AQM+ePTFp0iSPvXZSUhJq1Khhc33nzp149NFHPdYOIvIsH6UbQERUGQoLC6HVal1+flhYmBtbQ0TehiNAROQ1Ro8ejc2bN2P+/PmQJAmSJOH06dMAgCNHjmDQoEEICAhAeHg4Ro0ahfT0dPm5PXv2xBNPPIHJkycjNDQU/fr1AwDMnTsXCQkJ8Pf3R0xMDCZOnIjr168DAH799VeMGTMGmZmZ8utNnz4dgO0UWEpKCoYOHYqAgAAEBQVhxIgRuHjxovz49OnT0apVK3z22WeIi4tDcHAw7r33XmRnZ1dupxGRSxgAEZHXmD9/Pjp37oxHHnkEqampSE1NRUxMDFJTU9GjRw+0atUKu3btwo8//oiLFy9ixIgRVs//73//Cx8fH2zbtg0ffvghAEClUuG9997DoUOH8N///hcbN27Ev//9bwBAly5dMG/ePAQFBcmv9+yzz9q0SwiBu+++G1euXMHmzZuRnJyMf/75B4mJiVbl/vnnH6xZswbr1q3DunXrsHnzZrzxxhuV1FtEdDM4BUZEXiM4OBharRZ+fn6IiIiQry9atAht2rTB66+/Ll9bunQpYmJi8Pfff6Nhw4YAgNtuuw1vvvmmVZ0l84ni4+Mxa9YsTJgwAQsXLoRWq0VwcDAkSbJ6vdJ+/vlnHDhwAKdOnUJMTAwA4LPPPkOzZs2wc+dOtG/fHgBgMpmQlJSEwMBAAMCoUaPwyy+/4LXXXru5jiEit+MIEBF5vd27d2PTpk0ICAiQvxo3bgzAPOpi0a5dO5vnbtq0Cf369UN0dDQCAwPx4IMPIiMjAzk5OQ6//tGjRxETEyMHPwDQtGlT1KhRA0ePHpWvxcXFycEPAERGRuLSpUtO3SsReQZHgIjI65lMJtx5552YM2eOzWORkZHyv/39/a0eO3PmDAYNGoTx48dj1qxZqFWrFrZu3YqxY8fCYDA4/PpCCEiSVOF1jUZj9bgkSTCZTA6/DhF5DgMgIvIqWq0WRqPR6lqbNm2wcuVKxMXFwcfH8V9bu3btQlFREd555x2oVOYB76+++qrC1yutadOmSElJwdmzZ+VRoCNHjiAzMxNNmjRxuD1E5D04BUZEXiUuLg5//PEHTp8+jfT0dJhMJjz++OO4cuUK7rvvPvz55584efIkNmzYgIcffrjc4KV+/fooKirCggULcPLkSXz22WdYvHixzetdv34dv/zyC9LT05Gbm2tTT9++fdGiRQvcf//92LNnD/788088+OCD6NGjh91pNyLyfgyAiMirPPvss1Cr1WjatCnCwsKQkpKCqKgobNu2DUajEQMGDEDz5s3x9NNPIzg4WB7ZsadVq1aYO3cu5syZg+bNm+Pzzz/H7Nmzrcp06dIF48ePR2JiIsLCwmySqAHzVNaaNWtQs2ZNdO/eHX379kW9evWwYsUKt98/EXmGJIQQSjeCiIiIyJM4AkRERETVDgMgIiIiqnYYABEREVG1wwCIiIiIqh0GQERERFTtMAAiIiKiaocBEBEREVU7DICIiIio2mEARERERNUOAyAiIiKqdhgAERERUbXDAIiIiIiqnf8H4aQf1u3V0uAAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def fixed_point_iteration(q, x0, tol=1e-6, max_iter=1000):\n", + " start_time = time.perf_counter() # Starting timing\n", + " x = x0\n", + " for i in range(max_iter):\n", + " x_new = q(x)\n", + " if np.abs(x_new - x) < tol:\n", + " end_time = time.perf_counter() # Ending timing\n", + " print(f\"Converged in {i+1} iterations.\")\n", + " print(f\"Convergence took {(end_time - start_time):.6f} seconds.\")\n", + " return x_new\n", + " x = x_new\n", + " end_time = time.perf_counter() \n", + " print(\"Did not converge.\")\n", + " print(f\"Total time taken: {(end_time - start_time):.6f} seconds.\")\n", + " return x\n", + "\n", + "def control_function(x):\n", + " return np.cos(x)\n", + "\n", + "x0 = 0.5 # Initial guess\n", + "result = fixed_point_iteration(control_function, x0)\n", + "print(f\"Fixed point: {result}\")\n", + "\n", + "# Collecting the values to plot\n", + "x = x0\n", + "vals = []\n", + "for i in range(50):\n", + " x = np.cos(x)\n", + " vals.append(x)\n", + "\n", + "# Plottingn the values\n", + "plt.figure()\n", + "plt.plot(vals, marker='o', markersize=3)\n", + "plt.title('Convergence of Undamped Fixed-Point Iteration')\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel('cos(x)')\n", + "plt.grid(True)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Anderson acceleration implementation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged in 34 iterations.\n", + "Convergence took 0.003145 seconds.\n", + "Fixed point: 0.7390848229131413\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHFCAYAAAAaD0bAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByLUlEQVR4nO3deVyU1f4H8M8As7AroCyKQGQq7uIG5pKGW9lyM61bmlul3Ey05eq1rqIWWf0M66blvRqlllSmaZGC5i655ZJpZrlgCiKogCIwzJzfHzgPDDMDMzjMM8rn/Xpxr/PMec6c58zyfDurQgghQERERNSAuMhdACIiIiJHYwBEREREDQ4DICIiImpwGAARERFRg8MAiIiIiBocBkBERETU4DAAIiIiogaHARARERE1OAyAiIiIqMFhAGSDI0eOYOzYsYiIiIBGo4GXlxe6dOmCt99+G5cvX5a7eGRnqampaNu2Ldzd3aFQKHDo0KFaz1m3bh0UCgX8/f1RWlpaL+WaPXs2FApFveTtjP72t79BoVDghRdekLsoEoVCgdmzZ8v2+m+++SbWrl1r93zPnDkDhUKBlJSUGtNt3boVCoXC7N/w4cMByFdH/fr1Q79+/axKW73svr6+6NevH77//nubXzclJQUKhQJnzpyx+dwLFy5g9uzZVv3GAJX1//XXX0vHdu/ejdmzZ+Pq1as2v7491VQOW94bR2AAZKX//ve/iI6Oxr59+/DKK69gw4YNWLNmDR5//HF89NFHGD9+vNxFJDu6dOkSRo0ahcjISGzYsAGZmZm45557aj1v6dKlAIDLly/Xyw2qocnNzcV3330HAFi5ciVKSkpkLpFzqK8AqC7lyMzMNPpLSkoCAGRmZmLChAkyl7B2w4cPR2ZmJnbt2oUPP/wQOTk5GDZsmM1B0AMPPIDMzEwEBwfbXIYLFy4gMTHR6gDInN27dyMxMdEpAiBL5Vi0aBEWLVrk+EJZ4CZ3AW4HmZmZmDRpEuLi4rB27Vqo1Wrpubi4OLz00kvYsGGDjCW8dTqdDuXl5UbX1pD9/vvv0Gq1ePrpp9G3b1+rzsnJyUFaWhr69++P3bt3Y+nSpRg5cmQ9l/TWFRcXw8PDQ+5imPXZZ59Bq9XigQcewPfff49vvvkGf//73+Uull3dzt+9li1bomfPnmafs3Tc2QQGBkpljY2NRUxMDO6++24kJyfjgQcesDqfJk2aoEmTJvVVTNnY8/chKirKLvnYC1uArPDmm29CoVBgyZIlZn+kVCoVHnroIemxXq/H22+/jdatW0OtVqNp06YYPXo0/vrrL6Pz+vXrh3bt2mHfvn3o3bs3PDw8cNddd+Gtt96CXq8HUNESoVKp8Prrr5u87m+//QaFQoH3339fOpaTk4Pnn38ezZs3h0qlQkREBBITE1FeXi6lMTRzv/3225g3bx4iIiKgVquxZcsWAMC3336LDh06QK1W46677sLChQvNdrsIIbBo0SJ06tQJ7u7uaNy4MYYPH45Tp07ZfJ0GV69exUsvvYS77rpLqruhQ4fit99+k9KUlZVh3rx5Uv02adIEY8eOxaVLl8y/gdWsW7cOMTEx8PDwgLe3N+Li4pCZmSk9P2bMGNx7770AgJEjR0KhUFjVbPvpp5+ivLwcU6dOxd/+9jds3rwZZ8+eNUln6M5Zvnw52rRpAw8PD3Ts2FFq6ajq+++/R6dOnaBWqxEREYF3333X7Gvb+l5s374dsbGx8PDwwLhx4wAAP/74I/r16wd/f3+4u7ujRYsWeOyxx1BcXCydf/nyZcTHx6NZs2ZQqVS46667MHPmTJPuPluusSbLli1DYGAgPv30U7i7u2PZsmVm0+3ZswfDhg2Dv78/NBoNIiMjkZCQYJTmt99+w5NPPonAwECo1Wq0aNECo0ePNiq7Nd8fS271u1dSUoKXXnoJnTp1gq+vL/z8/BATE4Nvv/3W6HUUCgWuX7+OTz/9VOq6qfr5tPYaLly4gBEjRsDb2xu+vr4YOXIkcnJyar1Oa1XtAhNCYOjQofD390dWVpaUpri4GG3btkWbNm1w/fp1Ka01n2UhBN5++22EhYVBo9GgS5cu+OGHH2653JGRkWjSpInRd7e23wzAfBeYNb99W7duRbdu3QAAY8eOld5TW7oPZ8+ejVdeeQUAEBERIeWxdetWKU1qaipiYmLg6ekJLy8vDBo0CAcPHjTKZ8yYMfDy8sIvv/yCgQMHwtvbGwMGDAAAZGRk4OGHH0bz5s2h0Whw99134/nnn0deXp7V5TDXBebo3xQjgmpUXl4uPDw8RI8ePaw+57nnnhMAxAsvvCA2bNggPvroI9GkSRMRGhoqLl26JKXr27ev8Pf3Fy1bthQfffSRyMjIEPHx8QKA+PTTT6V0jz76qAgNDRU6nc7odV599VWhUqlEXl6eEEKI7OxsERoaKsLCwsTHH38sNm3aJObOnSvUarUYM2aMdN7p06cFANGsWTNx3333ia+//lqkp6eL06dPix9++EG4uLiIfv36iTVr1oivvvpK9OjRQ4SHh4vqH5dnn31WKJVK8dJLL4kNGzaIzz//XLRu3VoEBgaKnJwcm6+zsLBQtG3bVnh6eoo5c+aIjRs3itWrV4spU6aIH3/8UQghhE6nE4MHDxaenp4iMTFRZGRkiP/973+iWbNmIioqShQXF9f43qxcuVIAEAMHDhRr164VqampIjo6WqhUKrFjxw4hhBB//PGH+PDDDwUA8eabb4rMzEzx66+/1vq+33PPPSI4OFiUl5eLTZs2CQBi9uzZJukAiPDwcNG9e3fx5ZdfirS0NNGvXz/h5uYm/vzzTyndpk2bhKurq7j33nvFN998I7766ivRrVs30aJFi1t6L/z8/ERoaKj44IMPxJYtW8S2bdvE6dOnhUajEXFxcWLt2rVi69atYuXKlWLUqFHiypUrQgghbty4ITp06CA8PT3Fu+++K9LT08Xrr78u3NzcxNChQ+t0jTXZtWuXACBeeeUVIYQQTz/9tFAoFOLUqVNG6TZs2CCUSqXo0KGDSElJET/++KNYtmyZeOKJJ6Q0hw4dEl5eXiI8PFx89NFHYvPmzWLFihVixIgRorCwUAhh/ffHcH2zZs2SHtvju3f16lUxZswYsXz5cvHjjz+KDRs2iJdfflm4uLgYfU8yMzOFu7u7GDp0qMjMzDT6fFpbjuLiYtGmTRvh6+srPvjgA7Fx40bx4osvSp+tTz75pMb3ZsuWLQKASE1NFVqt1ujPUh3l5eWJ5s2bix49eoiysjIhhBDPPPOMcHd3F0eOHJHSWftZnjVrlgAgxo8fL3744QexZMkS0axZMxEUFCT69u1bY/mrlvEf//iH0bHLly8LFxcXERsbK4Sw7jdDCCE++eQTAUCcPn1aOmbNb19BQYF07muvvSa9p+fOnau1/r/66ishhBDnzp0TkydPFgDEN998I+VRUFAghBDijTfeEAqFQowbN05899134ptvvhExMTHC09PT6LftmWeeEUqlUoSHh4ukpCSxefNmsXHjRiGEEIsXLxZJSUli3bp1Ytu2beLTTz8VHTt2FK1atZLez9rK0bdvX6P3xtG/KdUxAKpFTk6OAGD0Y1qT48ePCwAiPj7e6PiePXsEAPGvf/1LOta3b18BQOzZs8cobVRUlBg0aJD0eN26dQKASE9Pl46Vl5eLkJAQ8dhjj0nHnn/+eeHl5SXOnj1rlN+7774rAEgfdMOPcGRkpPTBNejWrZsIDQ0VpaWl0rGioiLh7+9vdNPNzMwUAMT//d//GZ1/7tw54e7uLl599VWbr3POnDkCgMjIyBCWfPHFFwKAWL16tdHxffv2CQBi0aJFFs/V6XQiJCREtG/f3iiYLCoqEk2bNpV+8IQw/YGpzfbt2wUAMX36dCGEEHq9XkRERIiwsDCh1+uN0gIQgYGB0o1XiIrPmYuLi0hKSpKO9ejRQ4SEhIgbN25IxwoLC4Wfn98tvxebN282Svv1118LAOLQoUMWr/Gjjz4SAMSXX35pdHz+/Pkmn09rr7Em48aNEwDE8ePHhRCV78nrr79ulC4yMlJERkYa1VN1/fv3F40aNRK5ubkW01j7/TFcX9Wbuz2+e9WVl5cLrVYrxo8fLzp37mz0nKenp3jmmWfqfA2LFy8WAMS3335rlO7ZZ5+1KQAy93fy5EkhhGkdCSHEzp07hZubm0hISBDLli0TAMT//vc/6XlrP8tXrlwRGo1GPProo0bpDEGzLQFQfHy80Gq1oqysTBw/flwMGTJEABAffvihTb8ZlgIga377DL9ftdW7gbnfp3feecfk9YUQIisrS7i5uYnJkycbHS8qKhJBQUFixIgR0rFnnnlGABDLli2r8fX1er3QarXi7NmzJp8jS+UQwjQAcvRvSnXsArMzQzfSmDFjjI53794dbdq0webNm42OBwUFoXv37kbHOnToYNT8OmTIEAQFBeGTTz6Rjm3cuBEXLlyQui8A4LvvvsN9992HkJAQlJeXS39DhgwBAGzbts3odR566CEolUrp8fXr17F//3488sgjUKlU0nEvLy8MGzbM6NzvvvsOCoUCTz/9tNFrBQUFoWPHjkZNr9Ze5w8//IB77rkH999/Pyz57rvv0KhRIwwbNszodTt16oSgoCCT163qxIkTuHDhAkaNGgUXl8qPvpeXFx577DH89NNPRt09tjAMfja8HwqFAmPGjMHZs2dN3nMAuO++++Dt7S09DgwMRNOmTaX6uH79Ovbt24e//e1v0Gg0Ujpvb+9bfi8aN26M/v37Gx3r1KkTVCoVnnvuOXz66acm3Q1ARReZp6enNMvHwPBZr36dtV0jAKPylpeXo+J3Drh27Rq+/PJLxMbGonXr1gCAvn37IjIyEikpKVL3we+//44///wT48ePN6qnqoqLi7Ft2zaMGDGixjEatn5/buXc6t89g6+++gq9evWCl5cX3NzcoFQqsXTpUhw/ftzia9elHFu2bIG3t7dR1z0Am8dXzZ8/H/v27TP6Cw0NtZi+V69eeOONN5CcnIxJkybh6aefNppAYu1nOTMzEyUlJXjqqaeM8o+NjUVYWJjRMcMYK8Nf9W73RYsWQalUQqVSoU2bNti9ezfmzJmD+Ph4u/xmWPPbV582btyI8vJyjB492qgeNBoN+vbta/Y387HHHjM5lpubi4kTJyI0NFT6bBrq2trPZ3X18ZtiCwZAtQgICICHhwdOnz5tVfr8/HwAMDsTICQkRHrewN/f3ySdWq3GjRs3pMdubm4YNWoU1qxZI42sT0lJQXBwMAYNGiSlu3jxItavXw+lUmn017ZtWwAw6qs1V8YrV65ACIHAwECTMlU/dvHiRSlt9df76aefTF7Lmuu8dOkSmjdvbpKu+utevXoVKpXK5HVzcnJMXreq2t4bvV6PK1eu1Pj65hQVFeGrr75C9+7d0aRJE1y9ehVXr17Fo48+CoVCIQVHVdVWH1euXIFer0dQUJBJuurHbH0vzF1/ZGQkNm3ahKZNm+If//gHIiMjERkZiYULF0pp8vPzERQUZDIWrGnTpnBzc7P5s33mzBmT8hpu0Kmpqbh27RpGjBgh1WdBQQFGjBiBc+fOISMjAwCkcV81fW6uXLkCnU5n1WfLlu/PrZxr7j345ptvMGLECDRr1gwrVqxAZmYm9u3bh3Hjxlk9+83acuTn55v9npv7vNXkrrvuQteuXY3+ahvM/dRTT0GlUqG0tFQaL1K1/NZ8lg2fNWu+H5GRkUb5zJkzx+j5ESNGYN++fdi/fz9OnDiB/Px8acylPX4zrPntq08XL14EAHTr1s2kTlNTU00+mx4eHvDx8TE6ptfrMXDgQHzzzTd49dVXsXnzZuzduxc//fQTANT5Wuz9m2IrzgKrhaurKwYMGIAffvgBf/31V60/ooY3KDs72yTthQsXEBAQUKdyjB07Fu+88w5WrVqFkSNHYt26dUhISICrq6uUJiAgAB06dMAbb7xhNo+QkBCjx9U/dI0bN4ZCoZC+MFVVHxwZEBAAhUKBHTt2mP3Bq8uMliZNmpgMFK8uICAA/v7+FmfdVf2vg+qqvjfVXbhwAS4uLmjcuLENJa7wxRdfoLi4GHv37jV7/po1a3DlyhWb8ja8F+YGpd7qe2FpDaHevXujd+/e0Ol02L9/Pz744AMkJCQgMDAQTzzxBPz9/bFnzx4IIYzyyM3NRXl5uc2f7ZCQEOzbt8/oWKtWrQBUtqglJCSYDGY2PD9o0CCpRaemz42fnx9cXV2t+mzZ8v25lXPNvQcrVqxAREQEUlNTjZ63ZT0pa8vh7++PvXv3mjxvz0HQ5uh0Ojz11FNo3Lgx1Go1xo8fj127dkktztZ+lg3fZUvfj/DwcOnx+vXrjeqw+nvRpEkTdO3a1Wx56+s3w5EM38uvv/7apHXMHHOfzaNHj+Lw4cNISUnBM888Ix3/448/bqls9v5NsRUDICvMmDEDaWlpePbZZ/Htt98adQ8BgFarxYYNGzBs2DCpa2HFihXSyH4A2LdvH44fP46ZM2fWqQxt2rRBjx498Mknn0Cn06G0tBRjx441SvPggw8iLS0NkZGRdfpSenp6omvXrli7di3effdd6TqvXbtmMtL+wQcfxFtvvYXz589jxIgRdbqm6oYMGYJ///vf+PHHH026aKq+7qpVq6DT6dCjRw+b8m/VqhWaNWuGzz//HC+//LL0hbt+/TpWr14tzfKw1dKlS+Ht7Y21a9caNZMDwP79+/HKK69g5cqVNi3k5+npie7du+Obb77BO++8I3XvFBUVYf369UZp7f1euLq6okePHmjdujVWrlyJn3/+GU888QQGDBiAL7/8EmvXrsWjjz4qpf/ss88AQJotYi2VSmX2xnP8+HFkZmbiscceM1tn8+bNw7fffov8/Hzcc889iIyMxLJlyzBt2jSzN013d3f07dsXX331Fd544w2LP6q38v251e8eUHHjUalURjeCnJwck1lggOX/6rW2HPfddx++/PJLrFu3zqgb7PPPP69T2a01a9Ys7NixA+np6fD09ESfPn3wyiuvSC2N1n6We/bsCY1Gg5UrVxp11+zevRtnz541CoDat29f5/LW129GdYbP7a20ClnKY9CgQXBzc8Off/5ptmvLGobrrv79+vjjj60uhzn2/k2xFQMgK8TExGDx4sWIj49HdHQ0Jk2ahLZt20Kr1eLgwYNYsmQJ2rVrh2HDhqFVq1Z47rnn8MEHH8DFxQVDhgzBmTNn8PrrryM0NBRTp06tcznGjRuH559/HhcuXEBsbKz0X8sGc+bMQUZGBmJjY/Hiiy+iVatWKCkpwZkzZ5CWloaPPvqo1hasOXPm4IEHHsCgQYMwZcoU6HQ6vPPOO/Dy8jJa7bpXr1547rnnMHbsWOzfvx99+vSBp6cnsrOzsXPnTrRv3x6TJk2y6foSEhKQmpqKhx9+GNOnT0f37t1x48YNbNu2DQ8++CDuu+8+PPHEE1i5ciWGDh2KKVOmoHv37lAqlfjrr7+wZcsWPPzww0ZfpKpcXFzw9ttv46mnnsKDDz6I559/HqWlpXjnnXdw9epVvPXWWzaVF6j4L6O9e/di0qRJZoO2Xr164f/+7/+wdOlSm1cynjt3LgYPHiytNaXT6TB//nx4enra/b346KOP8OOPP+KBBx5AixYtUFJSIk05N4zJGj16ND788EM888wzOHPmDNq3b4+dO3fizTffxNChQ2scu2ULQ+vPq6++ajJ2AqgIAjdv3owVK1ZgypQp+PDDDzFs2DD07NkTU6dORYsWLZCVlYWNGzdi5cqVAIAFCxbg3nvvRY8ePTB9+nTcfffduHjxItatW4ePP/4Y3t7et/T9scd378EHH8Q333yD+Ph4DB8+HOfOncPcuXMRHByMkydPGqVt3749tm7divXr1yM4OBje3t5o1aqV1eUYPXo03nvvPYwePRpvvPEGWrZsibS0NGzcuLEub5lVMjIykJSUhNdff126sSUlJeHll19Gv3798Oijj1r9WW7cuDFefvllzJs3DxMmTMDjjz+Oc+fOYfbs2TZ349WkPn4zzImMjIS7uztWrlyJNm3awMvLCyEhITW2OlZnCPQWLlyIZ555BkqlEq1atUJ4eDjmzJmDmTNn4tSpUxg8eDAaN26MixcvYu/evfD09ERiYmKNebdu3RqRkZGYPn06hBDw8/PD+vXrpa5oa8phrnXeUb8pFtVp6HQDdejQIfHMM8+IFi1aCJVKJTw9PUXnzp3Fv//9b6PZJTqdTsyfP1/cc889QqlUioCAAPH000+bTGvs27evaNu2rcnrPPPMMyIsLMzkeEFBgXB3dxcAxH//+1+zZbx06ZJ48cUXRUREhFAqlcLPz09ER0eLmTNnimvXrgkhKmeivPPOO2bzWLNmjWjfvr1QqVSiRYsW4q233hIvvviiaNy4sUnaZcuWiR49eghPT0/h7u4uIiMjxejRo8X+/fvrdJ1XrlwRU6ZMES1atBBKpVI0bdpUPPDAA+K3336T0mi1WvHuu++Kjh07Co1GI7y8vETr1q3F888/L81AqcnatWtFjx49hEajEZ6enmLAgAFi165dRmmsnQWWkJBQ6+yp6dOnCwDiwIEDQgjzU2+FECIsLMxkZs+6detEhw4djN4Lw/Tf6m7lvcjMzBSPPvqoCAsLE2q1Wvj7+4u+ffuKdevWGaXLz88XEydOFMHBwcLNzU2EhYWJGTNmiJKSEqN0tlxjVWVlZaJp06aiU6dOFtOUl5eL5s2bi/bt2xuVf8iQIcLX11eo1WoRGRkppk6danTesWPHxOOPPy78/f2l+hwzZoxR2a35/hiur/oMJ3t899566y0RHh4u1Gq1aNOmjfjvf/9r9v0+dOiQ6NWrl/Dw8DCZ9WTtNfz111/iscceE15eXsLb21s89thjYvfu3TbNAqvp+1G1ji5cuCCaNm0q+vfvbzSbSq/Xi2HDholGjRoZzRqy5rOs1+tFUlKSCA0NFSqVSnTo0EGsX7/eZKZRTSx9Tquz5jfD0iwwa3/7vvjiC9G6dWuhVCrNfr6qslT/M2bMECEhIcLFxUUAEFu2bDG6hvvuu0/4+PgItVotwsLCxPDhw8WmTZuMyuXp6Wn2NY8dOybi4uKEt7e3aNy4sXj88cdFVlaW2bJaKoe596a+f1NqoriZMZFFWq0WnTp1QrNmzZCeni53cYiIiG4Zu8DIxPjx4xEXF4fg4GDk5OTgo48+wvHjx41mBBEREd3OGACRiaKiIrz88su4dOkSlEolunTpgrS0tPrvjyUiInIQdoERERFRg8OFEImIiKjBYQBEREREDQ4DICIiImpwOAjaDL1ejwsXLsDb29vitgFERETkXIQQKCoqQkhIiMnK/NUxADLjwoULNe5oTERERM7r3Llzta6+zgDIDMOS3efOnTPZFfdWabVapKenY+DAgVAqlXbNmyqxnh2D9ewYrGfHYV07Rn3Vc2FhIUJDQ2vcGNuAAZAZhm4vHx+fegmAPDw84OPjwy9XPWI9Owbr2TFYz47DunaM+q5na4avcBA0ERERNTgMgIiIiKjBYQBEREREDQ4DICIiImpwGAARERFRg8MAiIiIiBocBkBERETU4DAAIiIiogaHARARERE1OFwJ2sno9AJ7T19GblEJmnpr0D3CD64u3JCViIjInhgAOZENR7ORuP4YsgtKpGPBvhrMGhaFwe2CZSwZERHRnYVdYE5iw9FsTFrxs1HwAwA5BSWYtOJnbDiaLVPJiIiI7jwMgJyATi+QuP4YhJnnDMcS1x+DTm8uBREREdmKAZAT2Hv6sknLT1UCQHZBCfaevuy4QhEREd3BGAA5gdwiy8FPXdIRERFRzRgAOYGm3hq7piMiIqKaMQByAt0j/BDsq4Glye4KVMwG6x7h58hiERER3bEYADkBVxcFZg2LMvucISiaNSyK6wERERHZCQMgJzG4XTAWP90FjTyURseDfDVY/HQXrgNERERkR1wI0YkMbheMG2U6TP3yMADg/Sc74YH2IWz5ISIisjO2ADkZbZW1ftoE+TD4ISIiqgcMgJxMWble+ndxmU7GkhAREd25GAA5GQZARERE9Y8BkJMp01UGQDe05TKWhIiI6M7FAMjJVG0BulGmryElERER1RUDICdj3AXGFiAiIqL6wADIyRh3gXEMEBERUX1gAORkOAiaiIio/jEAcjJVW4AYABEREdUPBkBOxngQNMcAERER1QcGQE6GXWBERET1jwGQkzFuAWIAREREVB8YADkZzgIjIiKqfwyAnAy7wIiIiOqf7AHQokWLEBERAY1Gg+joaOzYsaPG9CtXrkTHjh3h4eGB4OBgjB07Fvn5+dLzKSkpUCgUJn8lJSX1fSl2wS4wIiKi+idrAJSamoqEhATMnDkTBw8eRO/evTFkyBBkZWWZTb9z506MHj0a48ePx6+//oqvvvoK+/btw4QJE4zS+fj4IDs72+hPo9E44pJumdE0eO4FRkREVC9kDYAWLFiA8ePHY8KECWjTpg2Sk5MRGhqKxYsXm03/008/ITw8HC+++CIiIiJw77334vnnn8f+/fuN0ikUCgQFBRn93S7YBUZERFT/3OR64bKyMhw4cADTp083Oj5w4EDs3r3b7DmxsbGYOXMm0tLSMGTIEOTm5uLrr7/GAw88YJTu2rVrCAsLg06nQ6dOnTB37lx07tzZYllKS0tRWloqPS4sLAQAaLVaaLXaul6iWYb8LOVbWl4Z9BSXltv99RuK2uqZ7IP17BisZ8dhXTtGfdWzLfkphBDCrq9upQsXLqBZs2bYtWsXYmNjpeNvvvkmPv30U5w4ccLseV9//TXGjh2LkpISlJeX46GHHsLXX38NpVIJoKKV6I8//kD79u1RWFiIhQsXIi0tDYcPH0bLli3N5jl79mwkJiaaHP/888/h4eFhh6u13pyfXZFfqgAAeLgKJHVnKxAREZE1iouL8fe//x0FBQXw8fGpMa3sAdDu3bsRExMjHX/jjTewfPly/PbbbybnHDt2DPfffz+mTp2KQYMGITs7G6+88gq6deuGpUuXmn0dvV6PLl26oE+fPnj//ffNpjHXAhQaGoq8vLxaK9BWWq0WGRkZiIuLk4K2qu59ZxsuFlaURemqwLHZcXZ9/Yaitnom+2A9Owbr2XFY145RX/VcWFiIgIAAqwIg2brAAgIC4OrqipycHKPjubm5CAwMNHtOUlISevXqhVdeeQUA0KFDB3h6eqJ3796YN28egoODTc5xcXFBt27dcPLkSYtlUavVUKvVJseVSmW9fQEs5a3VCaN/K1xc4eYq+2S921Z9vodUifXsGKxnx2FdO4a969mWvGS7s6pUKkRHRyMjI8PoeEZGhlGXWFXFxcVwcTEusqurKwDAUkOWEAKHDh0yGxw5o6qDoAEuhkhERFQfZGsBAoBp06Zh1KhR6Nq1K2JiYrBkyRJkZWVh4sSJAIAZM2bg/Pnz+OyzzwAAw4YNw7PPPovFixdLXWAJCQno3r07QkJCAACJiYno2bMnWrZsicLCQrz//vs4dOgQPvzwQ9mu0xZVp8EDFWsBeWv4XyFERET2JGsANHLkSOTn52POnDnIzs5Gu3btkJaWhrCwMABAdna20ZpAY8aMQVFREf7zn//gpZdeQqNGjdC/f3/Mnz9fSnP16lU899xzyMnJga+vLzp37ozt27eje/fuDr8+WwkhTFqAOBWeiIjI/mQNgAAgPj4e8fHxZp9LSUkxOTZ58mRMnjzZYn7vvfce3nvvPXsVz6Gqjv9Ru7mgtFzPAIiIiKgecHStE6na/eXrXtHtdYOrQRMREdkdAyAnUrX7yxAAsQWIiIjI/hgAORFDAOTqooCXpqJ3kgEQERGR/TEAciKGAEjl6gIPVcX0/hJOgyciIrI7BkBOpExXEeyo3FzgrmQLEBERUX1hAOREysorZoGp3CpbgBgAERER2R8DICdimAVWtQvsRhlngREREdkbAyAnYhgDpHZzgTtbgIiIiOoNAyAnYgiAlK4ucFcyACIiIqovDICcSNVB0JVdYAyAiIiI7I0BkBORpsG7ucBddXMWGKfBExER2R0DICdSamYdILYAERER2R8DICdStQVICoC4FxgREZHdMQByIobd4CsWQuQgaCIiovrCAMiJlJVXHQRdMQaIXWBERET2xwDIiVRdCJHrABEREdUfBkBOpOpmqOwCIyIiqj8MgJyI2UHQ3AqDiIjI7hgAOZFSnWkAVKzVQQghZ7GIiIjuOAyAnIjxQogVAZAQlesDERERkX0wAHIiZUYLIbpJxzkTjIiIyL4YADkRbZUuMFcXBVRuFW8Pt8MgIiKyLwZATsTQAqS+GfhwIDQREVH9YADkRAzrACldbwZAnApPRERULxgAOZGqg6ABQMPFEImIiOoFAyAnUnU3eADcEZ6IiKieMAByItVbgDyUFTPB2AJERERkXwyAnEiZzjgAMqwFdIOzwIiIiOyKAZATMWkB4iwwIiKiesEAyIkY1gFSuxq3ALELjIiIyL4YADkRSy1ADICIiIjsy632JOQohgBIWgfo5nYYdRkDpNML7D19GblFJWjqrUH3CD+4uijsV1giIqLbGAMgJ1J9ELRGWgjRtjFAG45mI3H9MWQXlEjHgn01mDUsCoPbBduptERERLcvdoE5kVI7dIFtOJqNSSt+Ngp+ACCnoASTVvyMDUez7VRaIiKi25fsAdCiRYsQEREBjUaD6Oho7Nixo8b0K1euRMeOHeHh4YHg4GCMHTsW+fn5RmlWr16NqKgoqNVqREVFYc2aNfV5CXZTdosLIer0Aonrj0GYec5wLHH9Mej05lIQERE1HLIGQKmpqUhISMDMmTNx8OBB9O7dG0OGDEFWVpbZ9Dt37sTo0aMxfvx4/Prrr/jqq6+wb98+TJgwQUqTmZmJkSNHYtSoUTh8+DBGjRqFESNGYM+ePY66rDoRQkhdYIbNUN2Vtq0DtPf0ZZOWH6PXAJBdUIK9py/fWmGJiIhuc7IGQAsWLMD48eMxYcIEtGnTBsnJyQgNDcXixYvNpv/pp58QHh6OF198EREREbj33nvx/PPPY//+/VKa5ORkxMXFYcaMGWjdujVmzJiBAQMGIDk52UFXVTc6vYC42TBT2QVm20rQuUWWg5+6pCMiIrpTyTYIuqysDAcOHMD06dONjg8cOBC7d+82e05sbCxmzpyJtLQ0DBkyBLm5ufj666/xwAMPSGkyMzMxdepUo/MGDRpUYwBUWlqK0tJS6XFhYSEAQKvVQqvV2nppNTLkVz3fqgOdFUIHrVYLlWtFRFRcWm5VOfw9rHs7/T3c7H5dzsZSPZN9sZ4dg/XsOKxrx6iverYlP9kCoLy8POh0OgQGBhodDwwMRE5OjtlzYmNjsXLlSowcORIlJSUoLy/HQw89hA8++EBKk5OTY1OeAJCUlITExEST4+np6fDw8LDlsqyWkZFh9Pi6FjC8HZsz0uGqAP4oqDh26UoB0tLSas1TL4BGKldcLQMAc1PeBRqpgEvHfkLa8Vsr/+2iej1T/WA9Owbr2XFY145h73ouLi62Oq3s0+AVCuMbtRDC5JjBsWPH8OKLL+Lf//43Bg0ahOzsbLzyyiuYOHEili5dWqc8AWDGjBmYNm2a9LiwsBChoaEYOHAgfHx86nJZFmm1WmRkZCAuLg5KpVI6nltUCuzfBoUCeHDoECgUCvxyvgAfHNsDF5U7hg7tY1X+yvCLmLzqsMlAaMXN/533t44Y1DbQ9MQ7jKV6JvtiPTsG69lxWNeOUV/1bOjBsYZsAVBAQABcXV1NWmZyc3NNWnAMkpKS0KtXL7zyyisAgA4dOsDT0xO9e/fGvHnzEBwcjKCgIJvyBAC1Wg21Wm1yXKlU1tsXoHreelQ026lcXaBSqQAA3u4VZbqh1Vldjgc7NYebmyte/uoIrpVWdqsFNdB1gOrzPaRKrGfHYD07DuvaMexdz7bkJdsgaJVKhejoaJPmr4yMDMTGxpo9p7i4GC4uxkV2da2YKSVujiCOiYkxyTM9Pd1ins6i+iKIQN33AhvcLhhP9giVHo/sFoqd/+zf4IIfIiIiS2TtAps2bRpGjRqFrl27IiYmBkuWLEFWVhYmTpwIoKJr6vz58/jss88AAMOGDcOzzz6LxYsXS11gCQkJ6N69O0JCQgAAU6ZMQZ8+fTB//nw8/PDD+Pbbb7Fp0ybs3LlTtuu0hmENIHWVAMgwC6ysXA+dXti0lcX10sqgSePmwm0wiIiIqpA1ABo5ciTy8/MxZ84cZGdno127dkhLS0NYWBgAIDs722hNoDFjxqCoqAj/+c9/8NJLL6FRo0bo378/5s+fL6WJjY3FqlWr8Nprr+H1119HZGQkUlNT0aNHD4dfny2qL4IIVC6ECFR0g3mprX+7Cm9UjoS/dK20hpREREQNj+yDoOPj4xEfH2/2uZSUFJNjkydPxuTJk2vMc/jw4Rg+fLg9iucwWjNdYGo3FygUgBAV0+RtCoBKKsf/XCpiAERERFSV7FthUIWyctMASKFQwENp23YYBkYtQAyAiIiIjDAAchKlZlqAAMDdxtWgDQpLGAARERFZwgDISRhagJSuxm9JXXaEB4DCG5VdYNfLdLheZUo8ERFRQ8cAyEmYGwQNVNkQ9RZagAAgjwOhiYiIJAyAnIS5MUBA1bWArG/BKdHqpPz8PCsWVWQ3GBERUSUGQE7CsBCi2s18F9gNrfUtQEU3Z4ApFECYf8VeZgyAiIiIKjEAchKWWoCkAMiGLjBD95e32g2B3hoAXAuIiIioKgZATkJaB6j6GKA6zAIzTIH3cVeiiXfFfmJsASIiIqrEAMhJlFpqAVLa3gVmWATRR8MAiIiIyBwGQE7CnoOgK1uA3BgAERERmcEAyEkYBkFXXweoLjvCF1VtAfK6GQBxDBAREZGEAZCTsDgIug7rABkGQXMMEBERkXkMgJyEIQBS26EFyNAF5q2p7ALLu1YKvV7Yo6hERES3PQZATsLyNPg6zAIztABplPD3qlgIUasTKLihrek0IiKiBoMBkJMos7AZauVCiLYMgr45BshdCbWbKxp5KAFwHBAREZEBAyAnUWZxHaBbGAOkqWg9kgZCcxwQERERAMBN7gJQhcouMFej43XZDb7qQogA0MRbjZO51+oUAOn0AntPX0ZuUQmaemvQPcIPri4Km/MhIiJyJgyAnEStW2HUYS8wH01FANS0jjPBNhzNRuL6Y8guKJGOBftqMGtYFAa3C7YpLyIiImfCLjAnYQiAlK7GrSsaZR1agEoqF0IEIM0Eyy0qsXhOdRuOZmPSip+Ngh8AyCkowaQVP2PD0Wyr8yIiInI2DICchOXd4CuCGJvGAN0wbgGydS0gnV4gcf0xmJs0bziWuP4YdJxWT0REtykGQE6iti6w4rJyCFF7wFFWrpe6y0wCICtnge09fdmk5acqASC7oAR7T1+2Kj8iIiJnwwDISUgBkKvxIGjDLDC9qNwwtSZFJZVr/XhJs8A0AKxvAbK2q8yWLjUiIiJnwgDISVhcB0hZGRBZ0w1m2AneW+0mzdaytQusqbfGrumIiIicDQMgJ2GpC8zN1UVaG8iamWDVp8ADlQHQlWKt9Do16R7hh2BfDSxNdlegYjZY9wi/WvMiIiJyRgyAnISlhRAB2/YDM0yB99ZUrnDQyF0Jt5utQfnXa28FcnVRYNawKLPPGYKiWcOiuB4QERHdthgAOQlLLUBAlbWArOoCq9wHzMDFRYEAG1eDHtwuGIuf7mJSniBfDRY/3YXrABER0W2NAZCTqBwEbaYFSFk5E6w2lV1gxmtc2joOCKgIggJvnmfw/Yu9GfwQEdFtjwGQk7A0CBqo0gVmzRggMy1AQN0CIK1Ojws3p8N7qysCqmMXCq0+n4iIyFkxAHICOr2QFhW85S6wKjvBV1WXDVHPX7kBnV5Ao3RBn3uaAAAO/3XV6vOJiIicFQMgJ1B1Zpb5FqCK1hdrBkFX3wnewNbFEAHg7OViAEALPw90DPUFABxhAERERHcAbobqBAzdX4D5MUCGtYBu2DQG6Na7wLLyrwMAwvw90aF5IwDAkb8KrD6/Ku4qT0REzoQBkBOo2gJUfTNUwLYd4avvBG9QlwDoTH5FC1CYnwfaNfOFQlGxBYYhiLEWd5UnIiJnwy4wJ1B1ALRCYRoA2bIOkKELzNseXWCGAMjfA15qN9zdxAsAcOSc9a1A3FWeiIickewB0KJFixAREQGNRoPo6Gjs2LHDYtoxY8ZAoVCY/LVt21ZKk5KSYjZNSYnz7ltV0xR4QL5B0FmXK7vAAFTpBrtq1fncVZ6IiJyVrAFQamoqEhISMHPmTBw8eBC9e/fGkCFDkJWVZTb9woULkZ2dLf2dO3cOfn5+ePzxx43S+fj4GKXLzs6GRuO8+1bVtAgiUHUdoFufBl9cpsP10trHEun1wqgFCIA0EPqwleOAuKs8ERE5K1kDoAULFmD8+PGYMGEC2rRpg+TkZISGhmLx4sVm0/v6+iIoKEj6279/P65cuYKxY8capVMoFEbpgoKCHHE5dVZbC5BNs8AsLIToqXaTWpKsaQXKLSpFabkeri4KhDRyB2DcAiRE7a023FWeiIiclWyDoMvKynDgwAFMnz7d6PjAgQOxe/duq/JYunQp7r//foSFhRkdv3btGsLCwqDT6dCpUyfMnTsXnTt3tphPaWkpSksrg4LCworF/rRaLbRarbWXZBVDflXzLS4tA1AxANrc66lvbgh/vbTm8pTr9Lh+M0hyd4NJ2gAvFbIu30D21eto5quqsZx/5la08jRrpAH0Omj1Otwd4A6lqwJXirU4fakQoY09aszD38O6j5e/h5tD6pnsj/XsGKxnx2FdO0Z91bMt+ckWAOXl5UGn0yEwMNDoeGBgIHJycmo9Pzs7Gz/88AM+//xzo+OtW7dGSkoK2rdvj8LCQixcuBC9evXC4cOH0bJlS7N5JSUlITEx0eR4eno6PDxqvsnXVUZGhvTvkwUKAK4oKylGWlqaSdqTuRXPnz2fjbS08xbzvK4FDG/pzi2bUH1CmavWFYAC6dt/Qq5/zS04P918TQ/ddaMyBWlcce66Ap+t34bOATXnoRdAI5UrrpYBMLu3vEAjFXDp2E9IO15jVkZ5/lmoQKEW8FECkT4CNc2mr1rPVH9Yz47BenYc1rVj2Luei4uLrU4r+zT46rOehBBmZ0JVl5KSgkaNGuGRRx4xOt6zZ0/07NlTetyrVy906dIFH3zwAd5//32zec2YMQPTpk2THhcWFiI0NBQDBw6Ej4+PDVdTO61Wi4yMDMTFxUGprBins+NkHnDsZ/j5+mDo0BjTk37JwRd/HoFXI38MHdrNYt5Zl4uB/TvhqXLFsAcGmjz/Q+FhnP71IkJbtsXQni1qLOdvGSeBP0+jS6sWGDq0cmf4n8qP4Yt9f8G16V0YOrhVrderDL+IF1YdNvucAgrM+1tHDGobaPb56jb+ehFJab8hp7CytS7IR43XhrY2ycNcPZP9sZ4dg/XsOKxrx6ivejb04FhDtgAoICAArq6uJq09ubm5Jq1C1QkhsGzZMowaNQoqVc1dOS4uLujWrRtOnjxpMY1arYZarTY5rlQq6+0LUDVv3c2hWGqlq9nX8/aouMbScn2N5bk5AQzeGvPlDvSpGAieX6yt9brOXa0Yl3NXE2+jtJ1b+OGLfX/h1+wiq+rmwU7NsevUZXyx95zRcQWA90Z2woOdmtWaB1AxnX7yqsMmM8ouFpZi8qrDFneor8/3kCqxnh2D9ew4rGvHsHc925KXbIOgVSoVoqOjTZq/MjIyEBsbW+O527Ztwx9//IHx48fX+jpCCBw6dAjBwc674F7ts8CsGwRtaQC0gS2LIRpmgLXwM+4C7HBzJtjR84XQWzl93TAT7IluoUge2QlBPhoIAJevl1l1fn1Mp9fpBTL/zMe3h84j8898TsUnImpgZO0CmzZtGkaNGoWuXbsiJiYGS5YsQVZWFiZOnAigomvq/Pnz+Oyzz4zOW7p0KXr06IF27dqZ5JmYmIiePXuiZcuWKCwsxPvvv49Dhw7hww8/dMg11UWZriKwudV1gCxNgTewLQAyXgPI4O4mXtAoXXCttByn8q7h7qbeNeZzvbQcu//IBwCMvzcCLQO9cb2sHDPXHMXSnacxOiYMbhau28CW6fQxkf61Xps9V6bmFh9ERLcnWQOgkSNHIj8/H3PmzEF2djbatWuHtLQ0aVZXdna2yZpABQUFWL16NRYuXGg2z6tXr+K5555DTk4OfH190blzZ2zfvh3du3ev9+upq1pbgKSVoGtev8fSIogG1q4GfbW4DIU3t9So3gLk5uqCdiG+2H/2Cg6fK6g1ANpx8hLKdHqE+Xvg7qYVK0k/1qU5FqT/jvNXbyDtaA4e6hhSYx72nE5vWJm6enuPYWVqS11plvJytkDKmnzsmWbP6cs4kKeA/+nLiLm7ab2+VkNOw3qu+XvBur790tRWz44g+yDo+Ph4xMfHm30uJSXF5Jivr2+No7zfe+89vPfee/YqnkPUug6QlQshWtoJ3qCJV8UYoNpagAx7gAX6qKXgq6oOzRth/9krOPLXVTwW3bzGvNKPXQQA3N8mUBrcrlG6YnRMON7b9DuWbP8TwzoE1zjw3dp9x2pLV1tXmgIVXWlxUUG1fhltCaRq+0GwNpCyRz71k8YVn53c76DXashpWM/mvhes69s5jfl6dhSFsGZFuwamsLAQvr6+KCgoqJdZYGlpaRg6dKg0WOt/O05h3vfH8VDHELz/pOl6RfnXShE9bxMA4M83h1q8OS/I+B3vbz6J0TFhmPOwafdgTkEJeiZthquLAifnDYGLhXy+PXQeU1YdQvdwP3w50XRWmuH5TqGNsPYfvSxea7lOj25vbMKVYi2+eLanUffU5etliH1rM0q0enz+bA/ERgZYzEenF+j55uYaW678vVTY+6/7pboxV8+Zf+bjyf/+ZDEPA0NZLQUcOr3AvfN/tNgtpwAQ5KvBzn/2R8axnBp/ECwFUoZ3xhBI1fbDYk0+AJiGae6YNPzc35lpbjUIsuX+LXsLEBlvhmqOh6rybbqh1cFLbf5tMwyCrr4RqoG/V8VsMp1e4EpxGfy9TGe+AUBWtS0wqjOsCH0suxBl5XqL5f456yquFGvh665Et/DGRs/5earweHQolv90Fv/dfqrGAAioGNhdUwB0raQch85dQXSYn8U01nalnb9ajA1HyywGHL7uKqvGJP3nxz+QvOl3i61EH/69M+Z+f7zWFim9HvjH55Zbm6zJZ/a6XwEomIZp7og0ieuPoX/rwFpbdJ2t3ExT+/tqTQu8vTAAcgLa8oqPhKVAQqN0gUIBCFExDshiAFTLIGilqwv8PFW4fL0Ml66VWgyAztQSAIX7e8BH44bCknL8frEI7Zr5mk236XhF91f/1k3NDnSe0DsCK/acxZYTl/D7xSLcE2h+PNEnu07jz0vXoXZzgY+70qgLL8hHAx93N/x+8RpGL92Lz8Z3R6fQxib9y3ohsPNkntn8q3t97VHc0OpNjucUlGDiip8RGeBp5ixT5oIfoHLm2ourDqLc9GWM0mUXlOCFL0yDn6r5vPDFQdQ0iU0ARmsnMQ3T3O5psgtK0HVehjRW0RnKxDS3nsaWySz2wADICdQ2C0yhUMBd6YriMh1KyizfMWsbBA1U7Ap/+XoZLhWVorWFLdIMu8C38Dd/o1coFOjQvBF2/pGHw39dtRwAVRn/Y06YvycGtw3CD0dz8N/tp/DO4x1N0vx56Rre2XgCADBrWFuM7BZq0i1VVq7H+E/3Yfef+fj7f/fAS+2G/OtlMPQvB3ip4K5yxbnLNyzWi3RtgNngB6gMOP7Mu15rPlXTW1JT8FNVbTP0OYOfGqKagh+6fTlyb0hZN0OlCoZB0GoLLUBA5VT4Yq3lL31tLUCAdVPhDWsAhVtoAQKADs0rgp4j5wrMPv/npWs4lXcdSlcF+txjuXvr2T53AQDWHjqP3ELjD75OL/DSl4dRWq5H75YBeLJ7KFxdFIiJ9MfDnZohJtIfri4KuKtcsfSZbmgV6IXScv3N4KdS3rUynLt8Ax4qV4zrFQ4FTDfmMBybPOBui2WtykfjZnZzD0NevjUEoUR068b1Cpe7CFQPrJ30Yg8MgJyAIQBS1rAeTuVUeMszwWpbCBGoPQAqLitH7s3nwvwsd/UYxgEd/uuq2ecNrT897/KHdw0BWZcWjdEtvDG0OoFPdp8xem7J9lM4dO4qvNVumP9YhxpniqncXFBwo+ZN8Lw1bpj5QBQWP90FQb7GX7IgXw0WP90FkU28aszDYPjN2W/mAinA+h9nP09VjYGUn6d1gVRt+QT5qBHko2Eaprkj0gT7ajB9SBsE+zpPmZjGPu9r9wjL4zjtjQGQE6htEDRQORW+psUQi242Cd9KC1DW5YrWH193JXw9LOfT8eaK0Cdzr5ktk2H8T1xU7ft8Pdu7ohVoReYZ/PhbLr49dB5f7T+HBekVXV//HhaFkEbuNeax9/TlWvuYLxaWYu/pyxjcLhg7/9kfXzzbEwuf6IQvnu2Jnf/sj8Htgq3+r4+4qKAaA6kX+res9cc52FeDeTdn61kKpOY93M4u+cx+qC1mPxTFNExzR6SZNSwKKjcXzBrmPGViGvu8r45cD4gBkBMorWUhRABwV9W+HYbUBVbLGCDA8mKIZ2sZAG0Q5KNBE281dHqBXy8UGD2Xf60UB85eAWB5/E9V97cJRKC3GkWlOoxL2Ycpqw7hla+PQKsXaN/MR2ptqYmtiyWa60oDgO4RflYFHN0j/GoMpFxdFLX+OM8aFoWhHYJrDKSGdgixSz6D2wVjcDumYZo7Jw0ApysT09jnfXUUrgNkhqPXAXrh85/x3ZFs/PvBKIy7N8LseU8u+QmZp/Kx8IlOeNjMBqJ6vUDkzDQIAex/7X4EWJjhtfbgeSSkHkJspD8+f7anyfNLtv+JN9N+w7COIfjAzJpEVU34dB82Hc/F6w9GYXyVcn+1/xxe+foI2ob44PsXe9eYB1CxMNbEFT9bfP4jK74Ytq7xU1t5Jt0sT9UvhyHgqI/Voh21WKK902T+kYv0HXswsHcPrprLepYlDev69kxTWz3XFdcBus3UthUGUPt+YEWl5TCEspbWAQJq7wKTWoD8am4BAirGAW06nosj1cYBGbq/rGn9MazObIkC1q0NYWi5ySkoMTv7SoGK/8qwpn/Z8F8p1QOOoDqsVjq4XTDiooJq/UEwtEjVdz72TtMjwg/5xwV6WLgxObo8d2oa1nPNWNe3X5ra6tkRGAA5Aa01Y4BqGQRtGACtUbpA7Wa6fYVBbfuBGcYAtailCwyonAn2y1+VXWAlWh22/16x3o4143/stdGpoctp0oqfoYD5lhtb+petDTisYe2PuKPyISIijgFyCoZB0NZMg7+htRAAWTEFHqgcA3S1WIvSctO8ztzcBT7cwhpAVRlmgp3Kuy7NwMr8Mx83tDoE+WjQNqT27kN7bnRq7/5lS+OEiIjo9scWICdQ22aoQOV2GJa6wKxZBBGomN2ldFVAqxPIv1ZmNLtKq9PjwtWKQKO2QdBAxdTrUD93nLt8A0fPF6DX3QHIMHR/RTWtcdq6gb02OjUwtNzUV/8yERHdGdgC5ATssQ5QbTvBG7i4KKQB0tXHAZ2/cgM6vYBG6YKm3uYHUVfXoVkjABXrAen1Apul6e8WlpmuxpZZV9Yy9C9HB8jbv0xERM6LAZATsGoavGEdIAsrQUtrAFmxAnFTCwOhDd1fYX6eVrXeAMYrQv9yvgAXC0vhqXJFz7usC1isnS7OIIaIiOyJAZATsGYhRA8rB0HXtOqygaWB0LYMgDYwjAM68tdVafZX31ZNahyIXZ2zrQ1BRER3Po4BcgLWTIO3VxcYYHkqvC1T4A3a32wBulBQguWZZwEA/Vs1tfp8A3vOuiIiIqoNAyAnIE2Dr3EQdM3rAFk7CBqonAlWfWbVWUMXWEDtM8AMdp68BFcXBXR6gas3W6He2XgCXhq3Os+6IiIiqm/sAnMC1uwG7640bIVhfgyQtdPgAfu1ABlWTNbpjZcezC0qxaQVP2PD0Wyr8iEiInI0BkBOwKaVoLV6s89bsxO8gbkASK8X0hgga6bAG1ZwNrfqsuFY4vpjJsERERGRM2AA5ARsGQR9w0ILkDU7wRuYGwR9sagEpeV6uLko0KyWndcB21ZwJiIicjYMgGSm1wtodRWtJPZYB6imfcAMmnhVzLa6VFQKw164hu6vZo3d4VZDOQzsuYIzERGRozEAkpmh9Qewch2g2maBWTEIOsBbBQAo0epxrbSi5SjrZgDUwsrxP/ZewZmIiMiRGADJzCgAsmIrjGKtTmq1qUqaBWZFF5iHyg1e6or8DOOAbNkDDKifFZyJiIgchQGQzAwDoIGaAyBDF5hOL4yCJqCiG62oxPpB0IDpQOizNgyABriCMxER3d4YAMnMsAaQ0lUBlxqCBcMgaMC0G+x6WTkMk62saQECKtcCMgyEtrULDOAKzkREdPviQogys2YneKBigLRhF/fiMh0aVYlTCm/OAFO5uUCjtG4LiqotQEKIyi4wGxZBBLiCMxER3Z4YAMnMmjWADNyVrtDqynFDa9wCVGTDIogGVQOgq8VaaRq9LS1ABlzBmYiIbjfsApOZNTvBGxgGQlfvAqscAG19PFs1ADKM/wn0UVvdgkRERHQ7q3MLkFarRU5ODoqLi9GkSRP4+XG2T12USWOArAmAzK8FJO0Eb8UUeIOqY4CkPcCsnAFGRER0u7OpBejatWv4+OOP0a9fP/j6+iI8PBxRUVFo0qQJwsLC8Oyzz2Lfvn31VdY7ki1dYIbWmer7gdmyE7yBUQtQHXaBJyIiup1ZHQC99957CA8Px3//+1/0798f33zzDQ4dOoQTJ04gMzMTs2bNQnl5OeLi4jB48GCcPHmyPst9x7B2EDRgeUf4yn3A6jYGSAqArJwCT0REdLuzuslg9+7d2LJlC9q3b2/2+e7du2PcuHFYvHgxli1bhm3btqFly5Z2K+idypqd4A0sbYdRaMM+YAaGACj/ehlO510DwC4wIiJqOKwOgL766iur0mk0GsTHx9e5QA2N1oqNUA2kMUBaSy1A1neB+XmqoFBULKz464VCAGwBIiKihqNOs8AuXrxo8bkjR47YlNeiRYsQEREBjUaD6Oho7Nixw2LaMWPGQKFQmPy1bdvWKN3q1asRFRUFtVqNqKgorFmzxqYyOZI1O8EbVM4CMx4DZMtO8AZKVxf4eVTsCWaYiRbmxxYgIiJqGOoUALVv3x7r1q0zOf7uu++iR48eVueTmpqKhIQEzJw5EwcPHkTv3r0xZMgQZGVlmU2/cOFCZGdnS3/nzp2Dn58fHn/8cSlNZmYmRo4ciVGjRuHw4cMYNWoURowYgT179th+oQ5QasMYIHdpDJDxVhh1GQQNVHaDAYCvuxK+HtYHUERERLezOgVA//znPzFy5EhMnDgRN27cwPnz59G/f3+88847SE1NtTqfBQsWYPz48ZgwYQLatGmD5ORkhIaGYvHixWbT+/r6IigoSPrbv38/rly5grFjx0ppkpOTERcXhxkzZqB169aYMWMGBgwYgOTk5Lpcar2zZRaYh2EWmNbCLDAbBkEDxgFQOLu/iIioAanTOkAvvfQS7r//fjz99NPo0KEDLl++jJ49e+LIkSMIDAy0Ko+ysjIcOHAA06dPNzo+cOBA7N6926o8li5divvvvx9hYWHSsczMTEydOtUo3aBBg2oMgEpLS1FaWio9LiysGBOj1Wqh1WqtKou1DPkZ/v9GWcX/uykUtb6W2q1ie4nrJcblKiiu+LeHsvY8qgrwrAyYmjd2t/u1yql6PVP9YD07BuvZcVjXjlFf9WxLfnVeCPGuu+5C27ZtsXr1agDAiBEjrA5+ACAvLw86nc7knMDAQOTk5NR6fnZ2Nn744Qd8/vnnRsdzcnJszjMpKQmJiYkmx9PT0+HhUT8tIxkZGQCAX84rALgi9+IFpKX9VeM5Z2+m/f3UWaSlnZaO5+S7AlDg14P7UPyH9WUouOgCQyPgxewL+O77v3CnbeFlqGeqX6xnx2A9Ow7r2jHsXc/FxcVWp61TALRr1y48/fTT8Pf3x5EjR7Br1y5MnjwZ33//PT7++GM0btzY6rwUCuM7rhDC5Jg5KSkpaNSoER555JFbznPGjBmYNm2a9LiwsBChoaEYOHAgfHx8ai2LLbRaLTIyMhAXFwelUolTW/4Esv7EXWGhGDq0bY3n5v2Uhe+yfoN/02AMHdpROj778BYAWgy8rw9aNvWyqhwbf72IfQd/BVDRnbY/zwV/lbnjtaGtMait9YGss6pez1Q/WM+OwXp2HNa1Y9RXPRt6cKxRpwCof//+mDp1KubOnQulUok2bdrgvvvuw6hRo9C+fXv89VfNLRkAEBAQAFdXV5OWmdzc3FpbkoQQWLZsGUaNGgWVSmX0XFBQkM15qtVqqNVqk+NKpbLevgCGvHWiIjDTKN1qfS1vTcW1lpTrpbRCCGkdIH9vd6vKu+FoNiavOgxR7fjFwlJMXnUYi5/ugsHtgm28IudUn+8hVWI9Owbr2XFY145h73q2Ja86DYJOT0/HW2+9ZfRCkZGR2LlzJ55//nmr8lCpVIiOjjZp/srIyEBsbGyN527btg1//PEHxo8fb/JcTEyMSZ7p6em15ikXW9YBMrcQ4g2tDjp9RSjjbcUsMJ1eIHH9MZPgB4B0LHH9MSlPIiKiO1GdWoD69u1r9riLiwtef/11q/OZNm0aRo0aha5duyImJgZLlixBVlYWJk6cCKCia+r8+fP47LPPjM5bunQpevTogXbt2pnkOWXKFPTp0wfz58/Hww8/jG+//RabNm3Czp07bbhCx7FtN/ib0+CrLIRo2AnezUUBdyt2ct97+jKyC0osPi8AZBeUYO/py4iJ9K81PyIiotuR1S1Aq1atsjrTc+fOYdeuXbWmGzlyJJKTkzFnzhx06tQJ27dvR1pamjSrKzs722RNoIKCAqxevdps6w8AxMbGYtWqVfjkk0/QoUMHpKSkIDU11ab1iRxJWgjRtfbgxd3MXmBVp8BbM3Yqt8hy8FOXdERERLcjq1uAFi9ejNmzZ2Ps2LF46KGH0KZNG6PnCwoKsGvXLqxYsQKbNm3C0qVLrco3Pj7e4tYZKSkpJsd8fX1rHeU9fPhwDB8+3KrXl5tN6wDdXAm6aheYtA2GlYsgNvXW2DUdERHR7cjqAGjbtm347rvv8MEHH+Bf//oXPD09ERgYCI1GgytXriAnJwdNmjTB2LFjcfToUTRt2rQ+y33HMARAStfaW2/MdoHZuAhi9wg/BPtqkFNQYnYckAJAkK8G3SP8rMqPiIjodmTTGKAHH3wQDz74IPLz87Fz506cOXMGN27cQEBAADp37ozOnTvDxaVO46obLJt2gzesBF1lLzDDGCBr9wFzdVFg1rAoTFrxMxSAURBkCMFmDYuC6522IBAREVEVdRoE7e/vj4cfftjeZWmQbNkM1TAGqESrh14v4OKiQJHUAmT9Wzm4XTAWP90FieuPGQ2IDvLVYNawqDtmCjwREZEldQqAzp07B4VCgebNmwMA9u7di88//xxRUVF47rnn7FrAO51tY4AqB0rf0OrgqXaT1gDyVtu2jsLgdsGIiwrC3tOXkVtUgqbeFd1ebPkhIqKGoE79VX//+9+xZcsWABVbT9x///3Yu3cv/vWvf2HOnDl2LeCdzpZZYBq3yjSGgdDSIGgbWoAMXF0UiIn0x8OdmiEm0p/BDxERNRh1CoCOHj2K7t27AwC+/PJLtG/fHrt378bnn39uduYWWWZLC5BLlbV+DFPhpUHQVo4BIiIiojoGQFqtVto6YtOmTXjooYcAAK1bt0Z2drb9StcA2BIAAaYzwaRB0FbOAiMiIqI6BkBt27bFRx99hB07diAjIwODBw8GAFy4cAH+/lw92BaVXWDWvRWV22FUBD6FdRgETURE1NDVKQCaP38+Pv74Y/Tr1w9PPvkkOnas2Jl83bp1UtcYWaeyBci68Tce1VaDrlwIkS1ARERE1qpTs0G/fv2Ql5eHwsJCNG7cWDr+3HPPwcPDw26FawikAMiKQdBA1bWAKgKgohJ2gREREdmqzv0mrq6uKC8vx86dO6FQKHDPPfcgPDzcjkVrGGxZBwio0gWmNR4Ebc1O8ERERFShTl1g169fx7hx4xAcHIw+ffqgd+/eCAkJwfjx42vdp4uMaW0eBF0R6NwoK4cQwuaVoImIiKiOAdC0adOwbds2rF+/HlevXsXVq1fx7bffYtu2bXjppZfsXcY7WmldW4DKdCgt10stSOwCIyIisl6d+k1Wr16Nr7/+Gv369ZOODR06FO7u7hgxYgQWL15sr/Ld0YQQVcYAWdkCVGUMkGEAtIsC8FRZN4aIiIiI6tgCVFxcjMDAQJPjTZs2ZReYDbS6yq1IbV0HqESrM9oJXqHgKs5ERETWqlMAFBMTg1mzZqGkpHIjzRs3biAxMRExMTF2K9ydztB9BdiyDlBFo11xmQ4FHP9DRERUJ3XqAktOTsaQIUPQvHlzdOzYEQqFAocOHYJarUZ6erq9y3jHMnR/Aba3ABWX6aSd4DkDjIiIyDZ1unO2b98eJ0+exIoVK/Dbb79BCIEnnngCTz31FNzd3e1dxjuWIQBydVFYvRFp5V5g5dJO8GwBIiIisk2dAqCkpCQEBgbi2WefNTq+bNkyXLp0Cf/85z/tUrg7na0DoAHjWWC3shM8ERFRQ1anMUAff/wxWrdubXLcsEcYWcfWRRAB481QuRM8ERFR3dQpAMrJyUFwcLDJ8SZNmnA3eBvYuhM8YDwGiDvBExER1U2dAqDQ0FDs2rXL5PiuXbsQEhJyy4VqKGzdCR4wngXGFiAiIqK6qdPgkQkTJiAhIQFarRb9+/cHAGzevBmvvvoqV4K2gaEFSF2HFqASLccAERER1VWd7pyvvvoqLl++jPj4eJSVlQEANBoN/vnPf2LGjBl2LeCdzBAAKW1pAZJWgi6XdoL3ZgsQERGRTeoUACkUCsyfPx+vv/46jh8/Dnd3d7Rs2RJqtdre5bujlekqdnSv8xggqQuMLUBERES2uKU7p5eXF7p162avsjQ4dRkEbZgGf8NoGjxbgIiIiGxRp0HQZB+ldVgHyENZEbOW6wXyr1d0P3IQNBERkW0YAMnIsBlqXVqAAOBqMQdBExER1QUDIBnVpQtM5eYCt2rbZrALjIiIyDYMgGRUVm77IGjAuBVIoQC8VGwBIiIisgUDIBkZFkJU2zAGCKicCQYAXmo3uFi5kSoRERFVYAAko7qsAwQAHlVafDgAmoiIyHYMgGRUlzFAQOViiADH/xAREdWF7AHQokWLEBERAY1Gg+joaOzYsaPG9KWlpZg5cybCwsKgVqsRGRmJZcuWSc+npKRAoVCY/JWUlNT3pdistA67wQPGY4C4CCIREZHtZL17pqamIiEhAYsWLUKvXr3w8ccfY8iQITh27BhatGhh9pwRI0bg4sWLWLp0Ke6++27k5uaivLzcKI2Pjw9OnDhhdEyj0dTbddRVXVuAqo4BYgsQERGR7WQNgBYsWIDx48djwoQJAIDk5GRs3LgRixcvRlJSkkn6DRs2YNu2bTh16hT8/PwAAOHh4SbpFAoFgoKC6rXs9qCtw27wQLUuMI4BIiIisplsAVBZWRkOHDiA6dOnGx0fOHAgdu/ebfacdevWoWvXrnj77bexfPlyeHp64qGHHsLcuXPh7u4upbt27RrCwsKg0+nQqVMnzJ07F507d7ZYltLSUpSWlkqPCwsLAQBarRZarfZWLtOEIT+tVouSsopp8G4KYdPraKq0GHmqXOxexjtB1Xqm+sN6dgzWs+Owrh2jvurZlvxkC4Dy8vKg0+kQGBhodDwwMBA5OTlmzzl16hR27twJjUaDNWvWIC8vD/Hx8bh8+bI0Dqh169ZISUlB+/btUVhYiIULF6JXr144fPgwWrZsaTbfpKQkJCYmmhxPT0+Hh4fHLV6peRkZGTid5QLABX+ePIG0679ZfW5eTsV5AHDxrzNISztVL2W8E2RkZMhdhAaB9ewYrGfHYV07hr3rubi42Oq0so+gVSiM17ARQpgcM9Dr9VAoFFi5ciV8fX0BVHSjDR8+HB9++CHc3d3Rs2dP9OzZUzqnV69e6NKlCz744AO8//77ZvOdMWMGpk2bJj0uLCxEaGgoBg4cCB8fn1u9RCNarRYZGRmIi4vDD4XHgLyL6Ni+LYb2MD/myZxDP5zA7tyzAIDo9m0wNDbMrmW8E1StZ6WS3YT1hfXsGKxnx2FdO0Z91bOhB8casgVAAQEBcHV1NWntyc3NNWkVMggODkazZs2k4AcA2rRpAyEE/vrrL7MtPC4uLujWrRtOnjxpsSxqtRpqtdrkuFKprLcvgFKpRLm+Yi8wjcq21/GqMu6nkaeaX9Ia1Od7SJVYz47BenYc1rVj2LuebclLtmnwKpUK0dHRJs1fGRkZiI2NNXtOr169cOHCBVy7dk069vvvv8PFxQXNmzc3e44QAocOHUJwcLD9Cm8nddkNHqg+DZ5fUCIiIlvJug7QtGnT8L///Q/Lli3D8ePHMXXqVGRlZWHixIkAKrqmRo8eLaX/+9//Dn9/f4wdOxbHjh3D9u3b8corr2DcuHHSIOjExERs3LgRp06dwqFDhzB+/HgcOnRIytOZ2GchRNl7MYmIiG47st49R44cifz8fMyZMwfZ2dlo164d0tLSEBZWMaYlOzsbWVlZUnovLy9kZGRg8uTJ6Nq1K/z9/TFixAjMmzdPSnP16lU899xzyMnJga+vLzp37ozt27eje/fuDr++2pTVcSFEjbIy/Zm86+gR4Q9X7gdGRERkNdmbD+Lj4xEfH2/2uZSUFJNjrVu3rnHU+HvvvYf33nvPXsWrV9o6BEAbjmZj/g+Vizz+a81RfPDjH5g1LAqD2zlfNx8REZEzkn0rjIbM0AVm7W7wG45mY9KKn3H1hvE6BzkFJZi04mdsOJpt9zISERHdiRgAyciWMUA6vUDi+mMQZp4zHEtcfww6vbkUREREVBUDIBnZEgDtPX0Z2QWWN3QVALILSrD39GV7FY+IiOiOxQBIRoZB0EorusByi6zbzd7adERERA0ZAyAZldrQAtTU27rd7K1NR0RE1JAxAJJRmQ0LIXaP8EOwrwaWJrsrAAT7atA9ws9+BSQiIrpDMQCSkWEavNqKFiBXFwVmDYsCAJMgyPB41rAorgdERERkBQZAMinX6WGYsGXtOkCD2wVj8dNdEORr3M0V5KvB4qe7cB0gIiIiK8m+EGJDZRgADdi2EOLgdsGIiwrC3tOXkVtUgqbeFd1ebPkhIiKyHgMgmZSVV67XY+tmqK4uCsRE+tu7SERERA0Gu8BkYmgBclEAbjYGQERERHRreOeViWEGmDVrABEREZF98e4rE1tWgSYiIiL74t1XJmU2TIEnIiIi++LdVyaGNYBsHQBNREREt453X5mwC4yIiEg+vPvKxNAFxgCIiIjI8Xj3lQlbgIiIiOTDu69MDAshcgwQERGR4/HuKxNDFxjXASIiInI83n1lwi4wIiIi+fDuKxOuA0RERCQf3n1louUsMCIiItnw7isTqQuMY4CIiIgcjndfmXAMEBERkXx495UJF0IkIiKSD+++MqlcB8hV5pIQERE1PAyAZCKtA+SmkLkkREREDQ8DIJkYxgCpOQiaiIjI4Xj3lQnHABEREcmHd1+ZaDkLjIiISDa8+8pEagFiFxgREZHD8e4rk8p1gDgLjIiIyNEYAMmEY4CIiIjkI/vdd9GiRYiIiIBGo0F0dDR27NhRY/rS0lLMnDkTYWFhUKvViIyMxLJly4zSrF69GlFRUVCr1YiKisKaNWvq8xLqxNACpHTlNHgiIiJHkzUASk1NRUJCAmbOnImDBw+id+/eGDJkCLKysiyeM2LECGzevBlLly7FiRMn8MUXX6B169bS85mZmRg5ciRGjRqFw4cPY9SoURgxYgT27NnjiEuyWpmuYiFE7gZPRETkeG5yvviCBQswfvx4TJgwAQCQnJyMjRs3YvHixUhKSjJJv2HDBmzbtg2nTp2Cn58fACA8PNwoTXJyMuLi4jBjxgwAwIwZM7Bt2zYkJyfjiy++qN8LsgH3AiMiIpKPbAFQWVkZDhw4gOnTpxsdHzhwIHbv3m32nHXr1qFr1654++23sXz5cnh6euKhhx7C3Llz4e7uDqCiBWjq1KlG5w0aNAjJyckWy1JaWorS0lLpcWFhIQBAq9VCq9XW5fIsMuRXVq4DALhA2P01qLKeWbf1i/XsGKxnx2FdO0Z91bMt+ckWAOXl5UGn0yEwMNDoeGBgIHJycsyec+rUKezcuRMajQZr1qxBXl4e4uPjcfnyZWkcUE5Ojk15AkBSUhISExNNjqenp8PDw8PWS7PK1cJrABT4ed8eFJyol5cgABkZGXIXoUFgPTsG69lxWNeOYe96Li4utjqtrF1gAKBQGA8CFkKYHDPQ6/VQKBRYuXIlfH19AVR0ow0fPhwffvih1ApkS55ARTfZtGnTpMeFhYUIDQ3FwIED4ePjU6frskSr1SIjIwNuanegpAR97u2Fjs197foaVFnPcXFxUCqVchfnjsV6dgzWs+Owrh2jvurZ0INjDdkCoICAALi6upq0zOTm5pq04BgEBwejWbNmUvADAG3atIEQAn/99RdatmyJoKAgm/IEALVaDbVabXJcqVTW2xdAe3MavIdaxS9ZParP95AqsZ4dg/XsOKxrx7B3PduSl2wjcFUqFaKjo02avzIyMhAbG2v2nF69euHChQu4du2adOz333+Hi4sLmjdvDgCIiYkxyTM9Pd1innLhOkBERETykfXuO23aNPzvf//DsmXLcPz4cUydOhVZWVmYOHEigIquqdGjR0vp//73v8Pf3x9jx47FsWPHsH37drzyyisYN26c1P01ZcoUpKenY/78+fjtt98wf/58bNq0CQkJCXJcokXSLDBuhUFERORwso4BGjlyJPLz8zFnzhxkZ2ejXbt2SEtLQ1hYGAAgOzvbaE0gLy8vZGRkYPLkyejatSv8/f0xYsQIzJs3T0oTGxuLVatW4bXXXsPrr7+OyMhIpKamokePHg6/vpoY1gFiCxAREZHjyT4IOj4+HvHx8WafS0lJMTnWunXrWkeNDx8+HMOHD7dH8eqFXgA6PQMgIiIiufDuK4ObvV8AGAARERHJgXdfGZSLyn9zDBAREZHj8e4rg6otQNwMlYiIyPEYAMnA0AKkcnOpcYFGIiIiqh8MgGRgaAFSs/uLiIhIFrwDy8DQAqTkAGgiIiJZ8A4sg5uLQHMANBERkUx4B5ZB1TFARERE5Hi8A8ugXF8x8JkBEBERkTx4B5aB1ALELjAiIiJZ8A4sA8MsMLYAERERyYN3YBlwDBAREZG8eAeWgbQOEAMgIiIiWfAOLANDAKTkGCAiIiJZ8A4sAw6CJiIikhfvwDLgIGgiIiJ58Q4sAx0HQRMREcmKd2AZsAWIiIhIXrwDy6Bc3FwJmmOAiIiIZME7sAw4DZ6IiEhevAPLgAshEhERyYt3YBlwHSAiIiJ58Q4sA7YAERERyYt3YBlIs8DYAkRERCQL3oFlwHWAiIiI5MU7sAy4DhAREZG8eAeWgWEMEKfBExERyYN3YBmU67kQIhERkZx4B5YBu8CIiIjkxTuwDAxdYFwHiIiISB68A8uALUBERETy4h1YBpwGT0REJC/egWXAhRCJiIjkJfsdeNGiRYiIiIBGo0F0dDR27NhhMe3WrVuhUChM/n777TcpTUpKitk0JSUljrgcq3AaPBERkbzc5Hzx1NRUJCQkYNGiRejVqxc+/vhjDBkyBMeOHUOLFi0snnfixAn4+PhIj5s0aWL0vI+PD06cOGF0TKPR2Lfwt4BjgIiIiOQlawC0YMECjB8/HhMmTAAAJCcnY+PGjVi8eDGSkpIsnte0aVM0atTI4vMKhQJBQUH2Lq7dcDNUIiIieckWAJWVleHAgQOYPn260fGBAwdi9+7dNZ7buXNnlJSUICoqCq+99hruu+8+o+evXbuGsLAw6HQ6dOrUCXPnzkXnzp0t5ldaWorS0lLpcWFhIQBAq9VCq9Xaemk1Ki0rg05ULISoEHq7508VDPXK+q1frGfHYD07DuvaMeqrnm3JT7YAKC8vDzqdDoGBgUbHAwMDkZOTY/ac4OBgLFmyBNHR0SgtLcXy5csxYMAAbN26FX369AEAtG7dGikpKWjfvj0KCwuxcOFC9OrVC4cPH0bLli3N5puUlITExEST4+np6fDw8LjFKzWm1QOGat+6eRPcZW2Du/NlZGTIXYQGgfXsGKxnx2FdO4a967m4uNjqtAohhLDrq1vpwoULaNasGXbv3o2YmBjp+BtvvIHly5cbDWyuybBhw6BQKLBu3Tqzz+v1enTp0gV9+vTB+++/bzaNuRag0NBQ5OXlGY01sofLRTfQ4+2Kgd5H/z0AaqWrXfOnClqtFhkZGYiLi4NSqZS7OHcs1rNjsJ4dh3XtGPVVz4WFhQgICEBBQUGt92/Z2h8CAgLg6upq0tqTm5tr0ipUk549e2LFihUWn3dxcUG3bt1w8uRJi2nUajXUarXJcaVSafcvgHApk/7toVHDxUVh1/zJWH28h2SK9ewYrGfHYV07hr3r2Za8ZBuFq1KpEB0dbdL8lZGRgdjYWKvzOXjwIIKDgy0+L4TAoUOHakzjSGU3p4ApXRUMfoiIiGQi6wiUadOmYdSoUejatStiYmKwZMkSZGVlYeLEiQCAGTNm4Pz58/jss88AVMwSCw8PR9u2bVFWVoYVK1Zg9erVWL16tZRnYmIievbsiZYtW6KwsBDvv/8+Dh06hA8//FCWa6yuTFcRAHERRCIiIvnIGgCNHDkS+fn5mDNnDrKzs9GuXTukpaUhLCwMAJCdnY2srCwpfVlZGV5++WWcP38e7u7uaNu2Lb7//nsMHTpUSnP16lU899xzyMnJga+vLzp37ozt27eje/fuDr8+cwwtQJwCT0REJB/Z5yDFx8cjPj7e7HMpKSlGj1999VW8+uqrNeb33nvv4b333rNX8exOCoDYAkRERCQb3oUdzNAFpmQLEBERkWx4F3awyhYgDoAmIiKSCwMgB+MgaCIiIvnxLuxgWl3FupMcBE1ERCQf3oUdjLPAiIiI5Me7sINxFhgREZH8eBd2MM4CIyIikh/vwg7GFiAiIiL58S7sYJwFRkREJD/ehR2schA01wEiIiKSCwMgB+MsMCIiIvnxLuxg0jpA7AIjIiKSDe/CDiaNAWILEBERkWx4F3YwzgIjIiKSH+/CDsZ1gIiIiOTHu7CDsQWIiIhIfrwLOxhngREREcmPd2EHq1wIkesAERERyYUBkINpOQuMiIhIdrwLO1hZOdcBIiIikhvvwg7GdYCIiIjkx7uwg3EWGBERkfx4F3YwzgIjIiKSH+/CDsYuMCIiIvnxLuxghhYgJbvAiIiIZMO7sINVrgPEqiciIpIL78IOVrkOEBdCJCIikgsDIAeT1gHiGCAiIiLZ8C7sYOwCIyIikh/vwg4khOA0eCIiIifAu7ADaXVC+jdbgIiIiOTDu7ADGbq/ALYAERERyYl3YQcydH8BXAeIiIhITrwLO9CNMh0AQAGB/WevQKcXtZxBRERE9UH2AGjRokWIiIiARqNBdHQ0duzYYTHt1q1boVAoTP5+++03o3SrV69GVFQU1Go1oqKisGbNmvq+jFptOJqNRxbtAgAIKPD0sv24d/6P2HA0W+aSERERNTyyBkCpqalISEjAzJkzcfDgQfTu3RtDhgxBVlZWjeedOHEC2dnZ0l/Lli2l5zIzMzFy5EiMGjUKhw8fxqhRozBixAjs2bOnvi/Hog1HszFpxc+4VFRqdDynoASTVvzMIIiIiMjBZA2AFixYgPHjx2PChAlo06YNkpOTERoaisWLF9d4XtOmTREUFCT9ubq6Ss8lJycjLi4OM2bMQOvWrTFjxgwMGDAAycnJ9Xw15un0Aonrj8FcZ5fhWOL6Y+wOIyIiciA3uV64rKwMBw4cwPTp042ODxw4ELt3767x3M6dO6OkpARRUVF47bXXcN9990nPZWZmYurUqUbpBw0aVGMAVFpaitLSytaZwsJCAIBWq4VWq7X2kszac/oysgtKLD4vAGQXlCDzj1z0iPC7pdeiSob37VbfP6oZ69kxWM+Ow7p2jPqqZ1vyky0AysvLg06nQ2BgoNHxwMBA5OTkmD0nODgYS5YsQXR0NEpLS7F8+XIMGDAAW7duRZ8+fQAAOTk5NuUJAElJSUhMTDQ5np6eDg8PD1svzciBPAUA11rTpe/Yg/zjbAWyt4yMDLmL0CCwnh2D9ew4rGvHsHc9FxcXW51WtgDIQKEw3hRUCGFyzKBVq1Zo1aqV9DgmJgbnzp3Du+++KwVAtuYJADNmzMC0adOkx4WFhQgNDcXAgQPh4+Nj0/VU53/6Mj47ub/WdAN792ALkB1ptVpkZGQgLi4OSqVS7uLcsVjPjsF6dhzWtWPUVz0benCsIVsAFBAQAFdXV5OWmdzcXJMWnJr07NkTK1askB4HBQXZnKdarYZarTY5rlQqb/mNibm7KYJ9NcgpKDE7DkgBIMhXg5i7m8LVhTvE25s93kOqHevZMVjPjsO6dgx717Mteck2CFqlUiE6Otqk+SsjIwOxsbFW53Pw4EEEBwdLj2NiYkzyTE9PtylPe3J1UWDWsCgAFcFOVYbHs4ZFMfghIiJyIFm7wKZNm4ZRo0aha9euiImJwZIlS5CVlYWJEycCqOiaOn/+PD777DMAFTO8wsPD0bZtW5SVlWHFihVYvXo1Vq9eLeU5ZcoU9OnTB/Pnz8fDDz+Mb7/9Fps2bcLOnTtluUYAGNwuGIuf7oLE9ceMBkQH+Wowa1gUBrcLruFsIiIisjdZA6CRI0ciPz8fc+bMQXZ2Ntq1a4e0tDSEhYUBALKzs43WBCorK8PLL7+M8+fPw93dHW3btsX333+PoUOHSmliY2OxatUqvPbaa3j99dcRGRmJ1NRU9OjRw+HXV9XgdsGIiwpC5h+5SN+xBwN792C3FxERkUxkHwQdHx+P+Ph4s8+lpKQYPX711Vfx6quv1prn8OHDMXz4cHsUz65cXRToEeGH/OMCPSL8GPwQERHJRPatMIiIiIgcjQEQERERNTgMgIiIiKjBYQBEREREDQ4DICIiImpwGAARERFRg8MAiIiIiBocBkBERETU4DAAIiIiogZH9pWgnZEQFfu2FxYW2j1vrVaL4uJiFBYWcqfhesR6dgzWs2Ownh2Hde0Y9VXPhvu24T5eEwZAZhQVFQEAQkNDZS4JERER2aqoqAi+vr41plEIa8KkBkav1+PChQvw9vaGQmHf/boKCwsRGhqKc+fOwcfHx655UyXWs2Ownh2D9ew4rGvHqK96FkKgqKgIISEhcHGpeZQPW4DMcHFxQfPmzev1NXx8fPjlcgDWs2Ownh2D9ew4rGvHqI96rq3lx4CDoImIiKjBYQBEREREDQ4DIAdTq9WYNWsW1Gq13EW5o7GeHYP17BisZ8dhXTuGM9QzB0ETERFRg8MWICIiImpwGAARERFRg8MAiIiIiBocBkBERETU4DAAcqBFixYhIiICGo0G0dHR2LFjh9xFuu1t374dw4YNQ0hICBQKBdauXWv0vBACs2fPRkhICNzd3dGvXz/8+uuv8hT2NpWUlIRu3brB29sbTZs2xSOPPIITJ04YpWE928fixYvRoUMHaXG4mJgY/PDDD9LzrGf7S0pKgkKhQEJCgnSM9Wwfs2fPhkKhMPoLCgqSnpe7nhkAOUhqaioSEhIwc+ZMHDx4EL1798aQIUOQlZUld9Fua9evX0fHjh3xn//8x+zzb7/9NhYsWID//Oc/2LdvH4KCghAXFyft90a127ZtG/7xj3/gp59+QkZGBsrLyzFw4EBcv35dSsN6to/mzZvjrbfewv79+7F//370798fDz/8sHRTYD3b1759+7BkyRJ06NDB6Djr2X7atm2L7Oxs6e+XX36RnpO9ngU5RPfu3cXEiRONjrVu3VpMnz5dphLdeQCINWvWSI/1er0ICgoSb731lnSspKRE+Pr6io8++kiGEt4ZcnNzBQCxbds2IQTrub41btxY/O9//2M921lRUZFo2bKlyMjIEH379hVTpkwRQvDzbE+zZs0SHTt2NPucM9QzW4AcoKysDAcOHMDAgQONjg8cOBC7d++WqVR3vtOnTyMnJ8eo3tVqNfr27ct6vwUFBQUAAD8/PwCs5/qi0+mwatUqXL9+HTExMaxnO/vHP/6BBx54APfff7/RcdazfZ08eRIhISGIiIjAE088gVOnTgFwjnrmZqgOkJeXB51Oh8DAQKPjgYGByMnJkalUdz5D3Zqr97Nnz8pRpNueEALTpk3Dvffei3bt2gFgPdvbL7/8gpiYGJSUlMDLywtr1qxBVFSUdFNgPd+6VatW4eeff8a+fftMnuPn2X569OiBzz77DPfccw8uXryIefPmITY2Fr/++qtT1DMDIAdSKBRGj4UQJsfI/ljv9vPCCy/gyJEj2Llzp8lzrGf7aNWqFQ4dOoSrV69i9erVeOaZZ7Bt2zbpedbzrTl37hymTJmC9PR0aDQai+lYz7duyJAh0r/bt2+PmJgYREZG4tNPP0XPnj0ByFvP7AJzgICAALi6upq09uTm5ppEv2Q/htkGrHf7mDx5MtatW4ctW7agefPm0nHWs32pVCrcfffd6Nq1K5KSktCxY0csXLiQ9WwnBw4cQG5uLqKjo+Hm5gY3Nzds27YN77//Ptzc3KS6ZD3bn6enJ9q3b4+TJ086xeeZAZADqFQqREdHIyMjw+h4RkYGYmNjZSrVnS8iIgJBQUFG9V5WVoZt27ax3m0ghMALL7yAb775Bj/++CMiIiKMnmc91y8hBEpLS1nPdjJgwAD88ssvOHTokPTXtWtXPPXUUzh06BDuuusu1nM9KS0txfHjxxEcHOwcn2eHDLUmsWrVKqFUKsXSpUvFsWPHREJCgvD09BRnzpyRu2i3taKiInHw4EFx8OBBAUAsWLBAHDx4UJw9e1YIIcRbb70lfH19xTfffCN++eUX8eSTT4rg4GBRWFgoc8lvH5MmTRK+vr5i69atIjs7W/orLi6W0rCe7WPGjBli+/bt4vTp0+LIkSPiX//6l3BxcRHp6elCCNZzfak6C0wI1rO9vPTSS2Lr1q3i1KlT4qeffhIPPvig8Pb2lu57ctczAyAH+vDDD0VYWJhQqVSiS5cu0jRiqrstW7YIACZ/zzzzjBCiYqrlrFmzRFBQkFCr1aJPnz7il19+kbfQtxlz9QtAfPLJJ1Ia1rN9jBs3TvqNaNKkiRgwYIAU/AjBeq4v1QMg1rN9jBw5UgQHBwulUilCQkLE3/72N/Hrr79Kz8tdzwohhHBMWxMRERGRc+AYICIiImpwGAARERFRg8MAiIiIiBocBkBERETU4DAAIiIiogaHARARERE1OAyAiIiIqMFhAEREZEZ4eDiSk5PlLgYR1RMGQEQkuzFjxuCRRx4BAPTr1w8JCQkOe+2UlBQ0atTI5Pi+ffvw3HPPOawcRORYbnIXgIioPpSVlUGlUtX5/CZNmtixNETkbNgCREROY8yYMdi2bRsWLlwIhUIBhUKBM2fOAACOHTuGoUOHwsvLC4GBgRg1ahTy8vKkc/v164cXXngB06ZNQ0BAAOLi4gAACxYsQPv27eHp6YnQ0FDEx8fj2rVrAICtW7di7NixKCgokF5v9uzZAEy7wLKysvDwww/Dy8sLPj4+GDFiBC5evCg9P3v2bHTq1AnLly9HeHg4fH198cQTT6CoqKh+K42I6oQBEBE5jYULFyImJgbPPvsssrOzkZ2djdDQUGRnZ6Nv377o1KkT9u/fjw0bNuDixYsYMWKE0fmffvop3NzcsGvXLnz88ccAABcXF7z//vs4evQoPv30U/z444949dVXAQCxsbFITk6Gj4+P9Hovv/yySbmEEHjkkUdw+fJlbNu2DRkZGfjzzz8xcuRIo3R//vkn1q5di++++w7fffcdtm3bhrfeequeaouIbgW7wIjIafj6+kKlUsHDwwNBQUHS8cWLF6NLly548803pWPLli1DaGgofv/9d9xzzz0AgLvvvhtvv/22UZ5VxxNFRERg7ty5mDRpEhYtWgSVSgVfX18oFAqj16tu06ZNOHLkCE6fPo3Q0FAAwPLly9G2bVvs27cP3bp1AwDo9XqkpKTA29sbADBq1Chs3rwZb7zxxq1VDBHZHVuAiMjpHThwAFu2bIGXl5f017p1awAVrS4GXbt2NTl3y5YtiIuLQ7NmzeDt7Y3Ro0cjPz8f169ft/r1jx8/jtDQUCn4AYCoqCg0atQIx48fl46Fh4dLwQ8ABAcHIzc316ZrJSLHYAsQETk9vV6PYcOGYf78+SbPBQcHS//29PQ0eu7s2bMYOnQoJk6ciLlz58LPzw87d+7E+PHjodVqrX59IQQUCkWtx5VKpdHzCoUCer3e6tchIsdhAERETkWlUkGn0xkd69KlC1avXo3w8HC4uVn/s7V//36Ul5fj//7v/+DiUtHg/eWXX9b6etVFRUUhKysL586dk1qBjh07hoKCArRp08bq8hCR82AXGBE5lfDwcOzZswdnzpxBXl4e9Ho9/vGPf+Dy5ct48sknsXfvXpw6dQrp6ekYN25cjcFLZGQkysvL8cEHH+DUqVNYvnw5PvroI5PXu3btGjZv3oy8vDwUFxeb5HP//fejQ4cOeOqpp/Dzzz9j7969GD16NPr27Wu2242InB8DICJyKi+//DJcXV0RFRWFJk2aICsrCyEhIdi1axd0Oh0GDRqEdu3aYcqUKfD19ZVadszp1KkTFixYgPnz56Ndu3ZYuXIlkpKSjNLExsZi4sSJGDlyJJo0aWIyiBqo6Mpau3YtGjdujD59+uD+++/HXXfdhdTUVLtfPxE5hkIIIeQuBBEREZEjsQWIiIiIGhwGQERERNTgMAAiIiKiBocBEBERETU4DICIiIiowWEARERERA0OAyAiIiJqcBgAERERUYPDAIiIiIgaHAZARERE1OAwACIiIqIGhwEQERERNTj/D2Zw9TrBDjUfAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def anderson_acceleration(f, x0, m=5, tol=1e-6, max_iter=100):\n", + " start_time = time.perf_counter() # Starting timing\n", + " x = [x0] # List to store iterates\n", + " g = [f(x0) - x0] # List to store residuals\n", + "\n", + " for k in range(1, max_iter + 1):\n", + " x_new = f(x[-1])\n", + " g_new = x_new - x[-1]\n", + "\n", + " if np.abs(g_new) < tol:\n", + " end_time = time.perf_counter() # Ending timing\n", + " print(f\"Converged in {k} iterations.\")\n", + " print(f\"Convergence took {(end_time - start_time):.6f} seconds.\")\n", + " return x_new\n", + "\n", + " x.append(x_new)\n", + " g.append(g_new)\n", + "\n", + " if k > 1:\n", + " G_k = np.column_stack([g[i] - g[i - 1] for i in range(max(0, k - m), k)])\n", + " X_k = np.column_stack([x[i] - x[i - 1] for i in range(max(0, k - m), k)])\n", + "\n", + " # Ensure g_new and g are column vectors for matrix operations\n", + " g_new = np.array(g_new).reshape(-1, 1)\n", + "\n", + " # Solve the least squares problem\n", + " try:\n", + " Q, R = np.linalg.qr(G_k)\n", + " gamma_k = np.linalg.solve(R, Q.T @ g_new)\n", + " x_accel = x[-1] - (X_k @ gamma_k + G_k @ gamma_k).flatten()\n", + " x[-1] = x_accel\n", + " g[-1] = f(x_accel) - x_accel\n", + " except np.linalg.LinAlgError:\n", + " pass\n", + "\n", + " end_time = time.perf_counter() \n", + " print(\"Did not converge.\")\n", + " print(f\"Total time taken: {(end_time - start_time):.6f} seconds.\")\n", + " return x[-1]\n", + "\n", + "def control_function(x):\n", + " return np.cos(x)\n", + "\n", + "x0 = 0.5 # initial guess\n", + "result = anderson_acceleration(control_function, x0)\n", + "print(f\"Fixed point: {result}\")\n", + "\n", + "# Collecting the values to plot\n", + "x_vals = [x0]\n", + "x = x0\n", + "for _ in range(50):\n", + " x = np.cos(x)\n", + " x_vals.append(x)\n", + "\n", + "# Plotting the values\n", + "plt.figure()\n", + "plt.plot(x_vals, marker='o')\n", + "plt.title('Convergence of Anderson-Accelerated Fixed-Point Iteration')\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel('cos(x)')\n", + "plt.grid(True)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Anderson Acceleration converged in 11 iterations.\n", + "Anderson Acceleration took 0.002159 seconds.\n", + "Fixed point (Anderson): 2.0134438112396547\n", + "Fixed Point Iteration converged in 11 iterations.\n", + "Fixed Point Iteration took 0.000042 seconds.\n", + "Fixed point (Fixed Point Iteration): 2.0134438112396547\n", + "Anderson Acceleration Time: 0.002159 seconds for 11 iterations\n", + "Fixed Point Iteration Time: 0.000042 seconds for 11 iterations\n" + ] + } + ], + "source": [ + "def anderson_acceleration(f, x0, m=3, tol=1e-6, max_iter=100):\n", + " start_time = time.perf_counter() # Start timing\n", + " x = [x0] # List to store iterates\n", + " g = [f(x0) - x0] # List to store residuals\n", + "\n", + " G_k = np.zeros((1, 0))\n", + " X_k = np.zeros((1, 0))\n", + "\n", + " for k in range(1, max_iter + 1):\n", + " x_new = f(x[-1])\n", + " g_new = x_new - x[-1]\n", + "\n", + " if np.abs(g_new) < tol:\n", + " end_time = time.perf_counter() # End timing\n", + " print(f\"Anderson Acceleration converged in {k} iterations.\")\n", + " print(f\"Anderson Acceleration took {(end_time - start_time):.6f} seconds.\")\n", + " return x_new, (end_time - start_time), k\n", + "\n", + " x.append(x_new)\n", + " g.append(g_new)\n", + "\n", + " g_new = np.array([[g_new]]) # Reshape to (1,1)\n", + " x_new = np.array([[x_new]]) # Reshape to (1,1)\n", + "\n", + " if k == 1:\n", + " G_k = g_new - np.array([[g[-2]]])\n", + " X_k = x_new - np.array([[x[-2]]])\n", + " else:\n", + " G_k = np.hstack((G_k, g_new - np.array([[g[-2]]])))\n", + " X_k = np.hstack((X_k, x_new - np.array([[x[-2]]])))\n", + "\n", + " # Use only the last m values\n", + " if G_k.shape[1] > m:\n", + " G_k = G_k[:, -m:]\n", + " X_k = G_k[:, -m:]\n", + "\n", + " # Solve the least squares problem\n", + " try:\n", + " Q, R = np.linalg.qr(G_k)\n", + " gamma_k = np.linalg.solve(R, Q.T @ g_new.T).flatten()\n", + " x_accel = x[-1] + g[-1] - (X_k @ gamma_k + G_k @ gamma_k)\n", + " x[-1] = x_accel\n", + " g[-1] = f(x_accel) - x_accel\n", + " except np.linalg.LinAlgError:\n", + " pass\n", + "\n", + " end_time = time.perf_counter() # End timing\n", + " print(\"Anderson Acceleration did not converge.\")\n", + " print(f\"Total time taken: {(end_time - start_time):.6f} seconds.\")\n", + " return x[-1], (end_time - start_time), k\n", + "\n", + "def fixed_point_iteration(q, x0, tol=1e-6, max_iter=1000):\n", + " start_time = time.perf_counter() # Start timing\n", + " x = x0\n", + " for i in range(max_iter):\n", + " x_new = q(x)\n", + " if np.abs(x_new - x) < tol:\n", + " end_time = time.perf_counter() # End timing\n", + " print(f\"Fixed Point Iteration converged in {i+1} iterations.\")\n", + " print(f\"Fixed Point Iteration took {(end_time - start_time):.6f} seconds.\")\n", + " return x_new, (end_time - start_time), i + 1\n", + " x = x_new\n", + " end_time = time.perf_counter() # End timing\n", + " print(\"Fixed Point Iteration did not converge.\")\n", + " print(f\"Total time taken: {(end_time - start_time):.6f} seconds.\")\n", + " return x, (end_time - start_time), max_iter\n", + "\n", + "def control_function(x):\n", + " return np.sin(x) + np.arctan(x)\n", + "\n", + "# Initial guess\n", + "x0 = 0.5\n", + "\n", + "# Running Anderson Acceleration\n", + "result_anderson, time_anderson, iterations_anderson = anderson_acceleration(control_function, x0)\n", + "print(f\"Fixed point (Anderson): {result_anderson}\")\n", + "\n", + "# Running Fixed Point Iteration\n", + "result_fixed_point, time_fixed_point, iterations_fixed_point = fixed_point_iteration(control_function, x0)\n", + "print(f\"Fixed point (Fixed Point Iteration): {result_fixed_point}\")\n", + "\n", + "# Compare times\n", + "print(f\"Anderson Acceleration Time: {time_anderson:.6f} seconds for {iterations_anderson} iterations\")\n", + "print(f\"Fixed Point Iteration Time: {time_fixed_point:.6f} seconds for {iterations_fixed_point} iterations\")\n", + "\n", + "# Collecting the values to plot\n", + "x = x0\n", + "vals = []\n", + "for _ in range(50):\n", + " x = np.cos(x)\n", + " vals.append(x)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "tardis", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tardis/io/configuration/config_reader.py b/tardis/io/configuration/config_reader.py index bc67fe0324f..ec57a524ed7 100644 --- a/tardis/io/configuration/config_reader.py +++ b/tardis/io/configuration/config_reader.py @@ -411,13 +411,19 @@ def validate_montecarlo_section(montecarlo_section): ] = Configuration.parse_convergence_section( montecarlo_section["convergence_strategy"] ) + elif montecarlo_section["convergence_strategy"]["type"] == "anderson": + montecarlo_section[ + "convergence_strategy" + ] = Configuration.parse_convergence_section( + montecarlo_section["convergence_strategy"] + ) elif montecarlo_section["convergence_strategy"]["type"] == "custom": raise NotImplementedError( 'convergence_strategy is set to "custom"; ' "you need to implement your specific convergence treatment" ) else: - raise ValueError('convergence_strategy is not "damped" or "custom"') + raise ValueError('convergence_strategy is not "damped", "anderson", or "custom"') @staticmethod def parse_convergence_section(convergence_section_dict): diff --git a/tardis/io/configuration/schemas/montecarlo.yml b/tardis/io/configuration/schemas/montecarlo.yml index 1baa9f37dac..88c43f7e370 100644 --- a/tardis/io/configuration/schemas/montecarlo.yml +++ b/tardis/io/configuration/schemas/montecarlo.yml @@ -48,6 +48,7 @@ properties: convergence_strategy: oneOf: - $ref: 'montecarlo_definitions.yml#/definitions/convergence_strategy/damped' + - $ref: 'montecarlo_definitions.yml#/definitions/convergence_strategy/anderson' - $ref: 'montecarlo_definitions.yml#/definitions/convergence_strategy/custom' default: type: 'damped' diff --git a/tardis/io/configuration/schemas/montecarlo_definitions.yml b/tardis/io/configuration/schemas/montecarlo_definitions.yml index 1adf22599ae..deb58b5407b 100644 --- a/tardis/io/configuration/schemas/montecarlo_definitions.yml +++ b/tardis/io/configuration/schemas/montecarlo_definitions.yml @@ -110,6 +110,114 @@ definitions: type: number default: -0.5 description: L=4*pi*r**2*T^y + anderson: + title: 'Anderson Convergence Strategy' + type: object + additionalProperties: false + properties: + type: + enum: + - anderson + memory: + type: number + default: 5 + description: "Number of past iterations to use for computing Anderson acceleration." + minimum: 1 + threshold: + type: number + default: 0.05 + description: "Specifies the threshold that is taken as convergence (i.e., 0.05 means that the value does not change more than 5%)." + minimum: 0 + damping_constant: + type: number + default: 1.0 + description: "Damping constant used to stabilize the acceleration if necessary." + minimum: 0 + hold_iterations: + type: number + multipleOf: 1.0 + default: 3 + description: "The number of iterations that the convergence criteria need to be fulfilled before TARDIS accepts the simulation as converged." + stop_if_converged: + type: boolean + default: false + description: "Stop plasma iterations before the number of specified iterations are reached if the simulation is plasma and inner boundary state is converged." + fraction: + type: number + default: 0.8 + description: the fraction of shells that have to converge to the given convergence + threshold. For example, 0.8 means that 80% of shells have to converge + to the threshold that convergence is established + minimum: 0 + t_inner: + type: object + additionalProperties: false + properties: + damping_constant: + type: number + default: 0.5 + description: damping constant + minimum: 0 + threshold: + type: number + description: specifies the threshold that is taken as convergence (i.e. + 0.05 means that the value does not change more than 5%) + minimum: 0 + type: + type: string + default: 'damped' + description: THIS IS A DUMMY VARIABLE DO NOT USE + t_rad: + type: object + additionalProperties: false + properties: + damping_constant: + type: number + default: 0.5 + description: damping constant + minimum: 0 + threshold: + type: number + description: specifies the threshold that is taken as convergence (i.e. + 0.05 means that the value does not change more than 5%) + minimum: 0 + type: + type: string + default: 'damped' + description: THIS IS A DUMMY VARIABLE DO NOT USE + required: + - threshold + w: + type: object + additionalProperties: false + properties: + damping_constant: + type: number + default: 0.5 + description: damping constant + minimum: 0 + threshold: + type: number + description: specifies the threshold that is taken as convergence (i.e. + 0.05 means that the value does not change more than 5%) + minimum: 0 + type: + type: string + default: 'damped' + description: THIS IS A DUMMY VARIABLE DO NOT USE + required: + - threshold + lock_t_inner_cycles: + type: number + multipleOf: 1.0 + default: 1 + description: The number of cycles to lock the update of the inner boundary + temperature. This process helps with convergence. The default is to switch + it off (1 cycle) + t_inner_update_exponent: + type: number + default: -0.5 + description: L=4*pi*r**2*T^y custom: $$target: 'montecarlo_definitions.yml#/definitions/convergence_strategy/custom' title: 'Custom Convergence Strategy' diff --git a/tardis/io/logger/colored_logger.py b/tardis/io/logger/colored_logger.py deleted file mode 100644 index a03d88a1fac..00000000000 --- a/tardis/io/logger/colored_logger.py +++ /dev/null @@ -1,83 +0,0 @@ -import logging - -""" -Code for Custom Logger Classes (ColoredFormatter and ColorLogger) and its helper function -(formatter_message) is used from this thread -http://stackoverflow.com/questions/384076/how-can-i-color-python-logging-output -""" - -FORMAT = "[$BOLD{name:20s}$RESET][{levelname:18s}] \n\t{message:s} ($BOLD{filename:s}$RESET:{lineno:d})" -DEBUG_FORMAT = "[$BOLD{name:20s}$RESET][{levelname:18s}] {message:s} ($BOLD{filename:s}$RESET:{lineno:d})" - - -def formatter_message(message, use_color=True): - """ - Helper Function used for Coloring Log Output - """ - # These are the sequences need to get colored ouput - RESET_SEQ = "\033[0m" - BOLD_SEQ = "\033[1m" - if use_color: - message = message.replace("$RESET", RESET_SEQ).replace( - "$BOLD", BOLD_SEQ - ) - else: - message = message.replace("$RESET", "").replace("$BOLD", "") - return message - - -class ColoredFormatter(logging.Formatter): - """ - Custom logger class for changing levels color - """ - - def __init__(self, use_color=True): - self.non_debug = formatter_message(FORMAT, True) - self.debug = formatter_message(DEBUG_FORMAT, True) - logging.Formatter.__init__(self, self.non_debug, style="{") - self.use_color = use_color - - def format(self, record): - COLOR_SEQ = "\033[1;%dm" - RESET_SEQ = "\033[0m" - BLACK, RED, GREEN, YELLOW, BLUE, MAGENTA, CYAN, WHITE = range(8) - COLORS = { - "WARNING": YELLOW, - "INFO": WHITE, - "DEBUG": BLUE, - "CRITICAL": YELLOW, - "ERROR": RED, - } - - levelname = record.levelname - if self.use_color and levelname in COLORS: - levelname_color = ( - COLOR_SEQ % (30 + COLORS[levelname]) + levelname + RESET_SEQ - ) - record.levelname = levelname_color - - if record.levelno == logging.DEBUG: - self._style._fmt = self.debug - else: - self._style._fmt = self.non_debug - - return logging.Formatter.format(self, record) - - -class ColoredLogger(logging.Logger): - """ - Custom logger class with multiple destinations - """ - - COLOR_FORMAT = formatter_message(FORMAT, True) - - def __init__(self, name): - logging.Logger.__init__(self, name, logging.DEBUG) - - color_formatter = ColoredFormatter(self.COLOR_FORMAT) - - console = logging.StreamHandler() - console.setFormatter(color_formatter) - - self.addHandler(console) - return diff --git a/tardis/io/logger/logger.py b/tardis/io/logger/logger.py index 70fcabcc4e9..311f5976260 100644 --- a/tardis/io/logger/logger.py +++ b/tardis/io/logger/logger.py @@ -1,66 +1,7 @@ import logging -import sys - -from tardis.io.logger.colored_logger import ColoredFormatter, formatter_message - -logging.captureWarnings(True) -logger = logging.getLogger("tardis") - -console_handler = logging.StreamHandler(sys.stdout) -console_formatter = ColoredFormatter() -console_handler.setFormatter(console_formatter) - -logger.addHandler(console_handler) -logging.getLogger("py.warnings").addHandler(console_handler) - -LOGGING_LEVELS = { - "NOTSET": logging.NOTSET, - "DEBUG": logging.DEBUG, - "INFO": logging.INFO, - "WARNING": logging.WARNING, - "ERROR": logging.ERROR, - "CRITICAL": logging.CRITICAL, -} -DEFAULT_LOG_LEVEL = "INFO" -DEFAULT_SPECIFIC_STATE = False - - -class FilterLog(object): - """ - Filter Log Class for Filtering Logging Output - to a particular level - - Parameters - ---------- - log_level : logging object - allows to have a filter for the - particular log_level - """ - - def __init__(self, log_level): - self.log_level = log_level - - def filter(self, log_record): - """ - filter() allows to set the logging level for - all the record that are being parsed & hence remove those - which are not of the particular level - - Parameters - ---------- - log_record : logging.record - which the paricular record upon which the - filter will be applied - - Returns - ------- - boolean : True, if the current log_record has the - level that of the specified log_level - False, if the current log_record doesn't have the - same log_level as the specified one - """ - return log_record.levelno == self.log_level - +import re +from ipywidgets import Output, Tab, Layout +from IPython.display import display, HTML def logging_state(log_level, tardis_config, specific_log_level): """ @@ -71,13 +12,14 @@ def logging_state(log_level, tardis_config, specific_log_level): Parameters ---------- - log_level: str - Allows to input the log level for the simulation - Uses Python logging framework to determine the messages that will be output - specific_log_level: boolean - Allows to set specific logging levels. Logs of the `log_level` level would be output. + log_level : str + Allows input of the log level for the simulation. + Uses Python logging framework to determine the messages that will be output. + tardis_config : dict + Configuration dictionary for TARDIS. + specific_log_level : bool + Allows setting specific logging levels. Logs of the `log_level` level would be output. """ - if "debug" in tardis_config: specific_log_level = ( tardis_config["debug"]["specific_log_level"] @@ -89,7 +31,6 @@ def logging_state(log_level, tardis_config, specific_log_level): log_level if log_level else tardis_config["debug"]["log_level"] ) - # Displays a message when both log_level & tardis["debug"]["log_level"] are specified if log_level and tardis_config["debug"]["log_level"]: print( "log_level is defined both in Functional Argument & YAML Configuration {debug section}" @@ -99,36 +40,30 @@ def logging_state(log_level, tardis_config, specific_log_level): ) else: - # Adds empty `debug` section for the YAML tardis_config["debug"] = {} if log_level: logging_level = log_level else: - tardis_config["debug"]["log_level"] = DEFAULT_LOG_LEVEL + tardis_config["debug"]["log_level"] = "INFO" logging_level = tardis_config["debug"]["log_level"] if not specific_log_level: - tardis_config["debug"][ - "specific_log_level" - ] = DEFAULT_SPECIFIC_STATE + tardis_config["debug"]["specific_log_level"] = False specific_log_level = tardis_config["debug"]["specific_log_level"] logging_level = logging_level.upper() - if not logging_level in LOGGING_LEVELS: + if not logging_level in ["NOTSET", "DEBUG", "INFO", "WARNING", "ERROR", "ALL"]: raise ValueError( - f"Passed Value for log_level = {logging_level} is Invalid. Must be one of the following {list(LOGGING_LEVELS.keys())}" + f"Passed Value for log_level = {logging_level} is Invalid. Must be one of the following ['NOTSET', 'DEBUG', 'INFO', 'WARNING', 'ERROR', 'ALL']" ) - # Getting the TARDIS logger & all its children loggers logger = logging.getLogger("tardis") + tardis_loggers = [logging.getLogger(name) for name in logging.root.manager.loggerDict if name.startswith("tardis")] - # Creating a list for Storing all the Loggers which are derived from TARDIS - tardis_loggers = tardis_logger() - - if logging_level in LOGGING_LEVELS: + if logging_level in ["NOTSET", "DEBUG", "INFO", "WARNING", "ERROR"]: for logger in tardis_loggers: - logger.setLevel(LOGGING_LEVELS[logging_level]) + logger.setLevel(getattr(logging, logging_level)) if logger.filters: for filter in logger.filters: @@ -136,7 +71,7 @@ def logging_state(log_level, tardis_config, specific_log_level): logger.removeFilter(filter) if specific_log_level: - filter_log = FilterLog(LOGGING_LEVELS[logging_level]) + filter_log = FilterLog([getattr(logging, logging_level), logging.INFO, logging.DEBUG]) for logger in tardis_loggers: logger.addFilter(filter_log) else: @@ -144,23 +79,142 @@ def logging_state(log_level, tardis_config, specific_log_level): for logger in tardis_loggers: logger.removeFilter(filter) +log_outputs = { + "WARNING/ERROR": Output(layout=Layout(height='300px', overflow_y='auto')), + "INFO": Output(layout=Layout(height='300px', overflow_y='auto')), + "DEBUG": Output(layout=Layout(height='300px', overflow_y='auto')), + "ALL": Output(layout=Layout(height='300px', overflow_y='auto')) +} + +tab = Tab(children=[log_outputs["WARNING/ERROR"], log_outputs["INFO"], log_outputs["DEBUG"], log_outputs["ALL"]]) +tab.set_title(0, "WARNING/ERROR") +tab.set_title(1, "INFO") +tab.set_title(2, "DEBUG") +tab.set_title(3, "ALL") + +display(tab) -def tardis_logger(): +def remove_ansi_escape_sequences(text): """ - Generates the list of the loggers which are derived from TARDIS - All loggers which are of the form `tardis.module_name` are added to the list + Remove ANSI escape sequences from a string. Parameters ---------- - list_for_loggers : list - List for storing the loggers derived from TARDIS + text : str + The input string containing ANSI escape sequences. Returns ------- - list_for_loggers : list + str + The cleaned string without ANSI escape sequences. + """ + ansi_escape = re.compile(r'\x1B[@-_][0-?]*[ -/]*[@-~]') + return ansi_escape.sub('', text) + +class WidgetHandler(logging.Handler): """ - list_for_loggers = [] - for name in logging.root.manager.loggerDict: - if not name.find("tardis"): - list_for_loggers.append(logging.getLogger(name)) - return list_for_loggers + A custom logging handler that outputs log messages to IPython widgets. + + Parameters + ---------- + logging.Handler : class + Inherits from the logging.Handler class. + """ + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + def emit(self, record): + """ + Emit a log record. + + Parameters + ---------- + record : logging.LogRecord + The log record to be emitted. + """ + log_entry = self.format(record) + clean_log_entry = remove_ansi_escape_sequences(log_entry) + + if record.levelno == logging.INFO: + color = '#D3D3D3' + elif record.levelno == logging.WARNING: + color = 'orange' + elif record.levelno == logging.ERROR: + color = 'red' + elif record.levelno == logging.CRITICAL: + color = 'orange' + elif record.levelno == logging.DEBUG: + color = 'blue' + else: + color = 'black' + + parts = clean_log_entry.split(' ', 2) + if len(parts) > 2: + prefix = parts[0] + levelname = parts[1] + message = parts[2] + html_output = f'{prefix} {levelname} {message}' + else: + html_output = clean_log_entry + + if record.levelno in (logging.WARNING, logging.ERROR): + with log_outputs["WARNING/ERROR"]: + display(HTML(f"
{html_output}
")) + elif record.levelno == logging.INFO: + with log_outputs["INFO"]: + display(HTML(f"
{html_output}
")) + elif record.levelno == logging.DEBUG: + with log_outputs["DEBUG"]: + display(HTML(f"
{html_output}
")) + with log_outputs["ALL"]: + display(HTML(f"
{html_output}
")) + +widget_handler = WidgetHandler() +widget_handler.setFormatter(logging.Formatter('%(name)s [%(levelname)s] %(message)s (%(filename)s:%(lineno)d)')) + +logging.captureWarnings(True) +logger = logging.getLogger("tardis") +logger.setLevel(logging.DEBUG) + +# To fix the issue of duplicate logs +for handler in logger.handlers[:]: + logger.removeHandler(handler) + +root_logger = logging.getLogger() +for handler in root_logger.handlers[:]: + root_logger.removeHandler(handler) + +logger.addHandler(widget_handler) +logging.getLogger("py.warnings").addHandler(widget_handler) + +class FilterLog(object): + """ + Filter Log Class for Filtering Logging Output + to a particular level. + + Parameters + ---------- + log_level : logging object + allows to have a filter for the + particular log_level + """ + def __init__(self, log_levels): + self.log_levels = log_levels + + def filter(self, log_record): + """ + Determine if the specified record is to be logged. + + Parameters + ---------- + log_record : logging.LogRecord + The log record to be filtered. + + Returns + ------- + boolean : True, if the current log_record has the + level that of the specified log_level, + False, if the current log_record doesn't have the + same log_level as the specified one. + """ + return log_record.levelno in self.log_levels \ No newline at end of file diff --git a/tardis/simulation/base.py b/tardis/simulation/base.py index 7efa377a314..e717aff9a64 100644 --- a/tardis/simulation/base.py +++ b/tardis/simulation/base.py @@ -517,7 +517,7 @@ def log_plasma_state( plasma_state_log["next_w"] = next_dilution_factor plasma_state_log.columns.name = "Shell No." - if is_notebook(): + if False and is_notebook(): #TODO: remove it with something better logger.info("\n\tPlasma stratification:") # Displaying the DataFrame only when the logging level is NOTSET, DEBUG or INFO diff --git a/tardis/simulation/convergence.py b/tardis/simulation/convergence.py index 3891b50374e..3b79a4a86de 100644 --- a/tardis/simulation/convergence.py +++ b/tardis/simulation/convergence.py @@ -1,5 +1,6 @@ import numpy as np - +from astropy import units as u +from astropy.units import Quantity class ConvergenceSolver: def __init__(self, strategy): @@ -21,9 +22,16 @@ def __init__(self, strategy): self.convergence_strategy = strategy self.damping_factor = self.convergence_strategy.damping_constant self.threshold = self.convergence_strategy.threshold + self.memory = self.convergence_strategy.get('memory', 5) + + #Initializing the history for Anderson method + self.x_history = [] + self.g_history = [] if self.convergence_strategy.type in ("damped"): self.converge = self.damped_converge + elif self.convergence_strategy.type in ("anderson"): + self.converge = self.anderson_converge elif self.convergence_strategy.type in ("custom"): raise NotImplementedError( "Convergence strategy type is custom; " @@ -32,10 +40,68 @@ def __init__(self, strategy): else: raise ValueError( f"Convergence strategy type is " - f"not damped or custom " + f"not damped, anderson or custom " f"- input is {self.convergence_strategy.type}" ) + def anderson_converge(self, value, estimated_value): + """ + Anderson acceleration method for convergence. + + Parameters + ---------- + value : Quantity or float + The current value of the physical property. + estimated_value : Quantity or float + The estimated value of the physical property. + + Returns + ------- + Quantity + The accelerated converged value with the same units as the input value. + + Raises + ------ + LinAlgError + If the QR decomposition or the linear system solution fails, the method + should fall back to returning the current estimated value without acceleration. + """ + if not isinstance(value, Quantity): + value = value * u.K + if not isinstance(estimated_value, Quantity): + estimated_value = estimated_value * u.K + + original_unit = value.unit + if value.unit != estimated_value.unit: + value = value.to(estimated_value.unit) + + value = value.value + estimated_value = estimated_value.value + + x_new = estimated_value + g_new = x_new - value + + self.x_history.append(value) + self.g_history.append(g_new) + + if len(self.x_history) > 1: + m = min(self.memory, len(self.x_history) - 1) + G_k = np.column_stack([self.g_history[i] - self.g_history[i - 1] for i in range(-m, 0)]) + X_k = np.column_stack([self.x_history[i] - self.x_history[i - 1] for i in range(-m, 0)]) + + g_new = np.array(g_new).reshape(-1, 1) + try: + Q, R = np.linalg.qr(G_k) + gamma_k = np.linalg.solve(R, Q.T @ g_new) + x_accel = self.x_history[-1] - (X_k @ gamma_k + G_k @ gamma_k).flatten() + self.x_history[-1] = x_accel + self.g_history[-1] = estimated_value - x_accel + return (x_accel * original_unit) + except np.linalg.LinAlgError: + pass + + return x_new * original_unit + def damped_converge(self, value, estimated_value): """Damped convergence solver diff --git a/tardis_anderson.yml b/tardis_anderson.yml new file mode 100644 index 00000000000..39eb6713307 --- /dev/null +++ b/tardis_anderson.yml @@ -0,0 +1,58 @@ +# Example YAML configuration for TARDIS +tardis_config_version: v1.0 + +supernova: + luminosity_requested: 9.44 log_lsun + time_explosion: 13 day + +atom_data: kurucz_cd23_chianti_H_He.h5 + +model: + structure: + type: specific + velocity: + start: 1.1e4 km/s + stop: 20000 km/s + num: 20 + density: + type: branch85_w7 + + abundances: + type: uniform + O: 0.19 + Mg: 0.03 + Si: 0.52 + S: 0.19 + Ar: 0.04 + Ca: 0.03 + +plasma: + disable_electron_scattering: no + ionization: lte + excitation: lte + radiative_rates_type: dilute-blackbody + line_interaction_type: macroatom + +montecarlo: + seed: 23111963 + no_of_packets: 4.0e+4 + iterations: 20 + nthreads: 1 + + last_no_of_packets: 1.e+5 + no_of_virtual_packets: 10 + + convergence_strategy: + type: damped + damping_constant: 0.2 + #memory: 10 + threshold: 0.05 + fraction: 0.8 + hold_iterations: 3 + t_inner: + damping_constant: 0.2 + +spectrum: + start: 500 angstrom + stop: 20000 angstrom + num: 10000