{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "0", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "import logging\n", "import warnings\n", "\n", "import numpy as np\n", "\n", "np.set_printoptions(precision=4)\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "\n", "logger = logging.getLogger()\n", "logger.setLevel(logging.ERROR)" ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "# Temperature dependent average observable\n", "\n", "This is a more difficult case because the derivatives are more complicated. However, this case has also been coded in the library. Extra data will need to be provided, but then everything is handled by flags in the code for creating the data and models. The catch is that we cannot foresee every way that an observable may depend on the extrapolation variable. But if the derivatives of the observable being averaged with respect to the extrapolation variable are known up to the desired order, we can incorporate these into our computation of extrapolation and interpolation coefficients. That means these must also be provided with the input data." ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "For our averaged quantity that depends explicitly on temperature, we select the average dimensionless potential energy of a single ideal gas particle, $\\beta u = \\langle \\beta a x \\rangle = \\langle \\beta x \\rangle$ where $a=1$ for simplicity. The average is over all particles and configurations. \n", "\n", "Everything can be set up in pretty much the same way as Case 1. EXCEPT that we now have to provide not only data for the observable, but also its derivatives in $\\beta$. The derivatives must be supplied up to the desired maximum order. This is a bit cumbersome in general, but hopefully your explicit dependence is polynomial (as it is here) or exponential in $\\beta$. In these cases, derivatives of arbitrarily high order are trivial to compute. \n", "\n", "So now the observable data can have 3 dimensions - for each configurational snapshot, for each derivative order starting at zero (i.e., the observable value) and going to the maximum desired, and for each element in the observable array. Note that the second dimension represents the derivatives of the observable with respect to the extrapolation variable. So for the first derivative, this is just $x$ in this example and for all subsequent derivatives up to the desired order, we pad with zeros." ] }, { "cell_type": "code", "execution_count": 2, "id": "3", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "# random number generator\n", "import cmomy\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import xarray as xr\n", "\n", "# Import idealgas module\n", "from thermoextrap import idealgas\n", "\n", "rng = cmomy.random.default_rng(seed=0)\n", "\n", "\n", "# Define test betas and reference beta\n", "betas = np.arange(0.1, 10.0, 0.5)\n", "beta_ref = betas[11]\n", "vol = 1.0\n", "\n", "# Define orders to extrapolate to\n", "orders = [1, 2, 4, 6]\n", "order = orders[-1]\n", "\n", "npart = 1000 # Number of particles (in single configuration)\n", "nconfig = 100_000 # Number of configurations\n", "\n", "# Generate all the data we could want\n", "xdata, udata = idealgas.generate_data((nconfig, npart), beta_ref, vol)" ] }, { "cell_type": "code", "execution_count": 3, "id": "4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Size: 6MB\n", "array([[1.0138, 0.9448, 0.9595, ..., 0.9873, 0.986 , 1.0242],\n", " [0.181 , 0.1687, 0.1713, ..., 0.1763, 0.1761, 0.1829],\n", " [0. , 0. , 0. , ..., 0. , 0. , 0. ],\n", " ...,\n", " [0. , 0. , 0. , ..., 0. , 0. , 0. ],\n", " [0. , 0. , 0. , ..., 0. , 0. , 0. ],\n", " [0. , 0. , 0. , ..., 0. , 0. , 0. ]])\n", "Coordinates:\n", " * deriv (deriv) int64 56B 0 1 2 3 4 5 6\n", "Dimensions without coordinates: rec\n" ] } ], "source": [ "# You can use numpy and just wrap with xarray after...\n", "xdataDepend = np.array([xdata * beta_ref, xdata])\n", "xdataDepend = np.vstack((xdataDepend, np.zeros((orders[-1] - 1, xdata.shape[0]))))\n", "xdata = xr.DataArray(xdata, dims=[\"rec\"])\n", "udata = xr.DataArray(udata, dims=[\"rec\"])\n", "# Note the naming of dimensions here and see \"Data_Organization\" for a full explanation\n", "xdataDepend = xr.DataArray(\n", " xdataDepend,\n", " dims=[\"deriv\", \"rec\"],\n", " coords={\"deriv\": np.arange(xdataDepend.shape[0])},\n", ")\n", "print(xdataDepend)\n", "\n", "# # Or you can accomplish this with xarray if you really want to...\n", "# xdata = xr.DataArray(xdata, dims=['rec'])\n", "# udata = xr.DataArray(udata, dims=['rec'])\n", "# xdataDepend = (\n", "# xr.concat([xdata*beta_ref, xdata], dim='deriv')\n", "# .assign_coords(deriv=lambda x: np.arange(x.sizes['deriv']))\n", "# .reindex(deriv=np.arange(orders[-1] + 1))\n", "# .fillna(0.0)\n", "# )\n", "# print(xdataDepend)" ] }, { "cell_type": "code", "execution_count": 4, "id": "5", "metadata": {}, "outputs": [], "source": [ "import thermoextrap as xtrap" ] }, { "cell_type": "code", "execution_count": 5, "id": "6", "metadata": {}, "outputs": [], "source": [ "# Create the model and data, here with the data object just created inside the call for the model\n", "# The data is always accessible directly as model.data\n", "order = orders[-1]\n", "xem_dep = xtrap.beta.factory_extrapmodel(\n", " beta=beta_ref,\n", " data=xtrap.DataCentralMomentsVals.from_vals(\n", " uv=udata,\n", " xv=xdataDepend,\n", " # by specifying deriv dimension, trigger doing coef calculation with explicit derivatives\n", " deriv_dim=\"deriv\",\n", " order=orders[-1],\n", " central=True,\n", " ),\n", ")" ] }, { "cell_type": "code", "execution_count": 6, "id": "7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Model parameters (derivatives):\n", " Size: 56B\n", "array([ 0.9792, 0.0167, -0.0127, -0.0042, -0.6049, 1.8554, 2.1772])\n", "Dimensions without coordinates: order\n", "\n", "\n", "Model predictions:\n", " Size: 32B\n", "array([0.6949, 0.7366, 0.7752, 0.8106])\n", "Coordinates:\n", " * beta (beta) float64 32B 0.1 0.6 1.1 1.6\n", " dalpha (beta) float64 32B -5.5 -5.0 -4.5 -4.0\n", " beta0 float64 8B 5.6\n", "\n", "\n" ] } ], "source": [ "# Check the derivatives\n", "print(\"Model parameters (derivatives):\")\n", "print(xem_dep.derivs(norm=False))\n", "print(\"\\n\")\n", "\n", "# Finally, look at predictions\n", "print(\"Model predictions:\")\n", "print(xem_dep.predict(betas[:4], order=2))\n", "print(\"\\n\")" ] }, { "cell_type": "code", "execution_count": 7, "id": "8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bootstrapped uncertainties in predictions:\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " Size: 32B\n", "array([0.0997, 0.0824, 0.0667, 0.0527])\n", "Coordinates:\n", " * beta (beta) float64 32B 0.1 0.6 1.1 1.6\n", " dalpha (beta) float64 32B -5.5 -5.0 -4.5 -4.0\n", " beta0 float64 8B 5.6\n" ] } ], "source": [ "# And bootstrapped uncertainties\n", "print(\"Bootstrapped uncertainties in predictions:\")\n", "print(xem_dep.resample(nrep=100, parallel=True).predict(betas[:4], order=2).std(\"rep\"))" ] }, { "cell_type": "code", "execution_count": 8, "id": "9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True extrapolation coefficients: [ 9.7922e-01 1.7150e-02 -1.3566e-02 1.0069e-02 -6.7219e-03 3.6244e-03\n", " -9.3213e-04]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\t With N_configs = 10: [ 9.7620e-01 6.3382e-02 -1.5816e-01 8.3563e-01 1.7880e+01 1.0777e+02\n", " -3.1199e+03]\n", "\t With N_configs = 100: [ 9.8067e-01 4.3129e-02 1.6502e-01 -1.1619e+00 1.0416e+01 6.0027e+01\n", " -3.1348e+03]\n", "\t With N_configs = 1000: [ 9.8033e-01 1.9356e-02 3.7339e-02 1.8383e-01 -8.0818e-01 -4.5850e+01\n", " -6.0897e+01]\n", "\t With N_configs = 10000: [ 9.7936e-01 1.7196e-02 2.1037e-02 2.1810e-01 -4.0756e+00 2.2653e+01\n", " -2.0604e+01]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\t With N_configs = 100000: [ 0.9792 0.0167 -0.0127 -0.0042 -0.6049 1.8554 2.1772]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHBCAYAAAAPRE2qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAACwl0lEQVR4nOzdeVxUVf8H8M/s7AMIsiiLioKgLAopUqLhBuVW5pqplUtabqnPQ9avtHrsMS1N01wqH5fM1Nyy3HDfVzQXcINAFgHBGdZZ7++PkSvDLGzDMMD3/XrdFzP3nnvPGcSZz5x77rkchmEYEEIIIYSQRo3b0A0ghBBCCCF1R6GOEEIIIaQJoFBHCCGEENIEUKgjhBBCCGkCKNQRQgghhDQBFOoIIYQQQpoACnWEEEIIIU0Av6Eb0NDUajUyMzNhb28PDofT0M0hhBBCSDPGMAwKCwvh6ekJLrdmfW/NPtRlZmbCy8uroZtBCCGEEMJKT09H69ata7RPsw919vb2ADS/PAcHhwZuDSGEEEKaM6lUCi8vLzaf1ESzD3Xlp1wdHBwo1BFCCCHEItRmSBhdKEEIIYQQ0gRQqCOEEEIIaQIo1BFCCCGENAEU6gghhBBCmgAKdYQQQgghTQCFOkIIIYSQJoBCHSGEEEJIE0ChjhBCCCGkCaBQRwghhBDSBFCoI4QQQghpAijUEUIIIYQ0ARTqCCGEEEKaAAp1hBBCCCFNAIU6QgghhJAmgEIdIYQQQkgTQKGOEEIIIaQJoFBHCLF4ycnJcHd3R2FhYbX3WblyJQYNGlSPrbJsn332GUJDQy3mOISQ+mdRoe7kyZMYOHAgPD09weFwsHv3bqPljx8/Dg6Ho7MkJSWZp8GEkFpLT0/HO++8A09PTwiFQvj4+GDGjBl48uSJTtn58+dj2rRpsLe3BwCUlZVh/Pjx6Ny5M/h8PoYMGaKzz8SJE3Hp0iWcPn3aaDvGjx+v931kwIAB1X4tvXr1wsyZM6td3lLpe9+dM2cOEhISGqZBhJAasahQV1xcjJCQEKxcubJG+yUnJyMrK4td2rdvX08tJISYwsOHDxEeHo67d+9i69atuH//Pn744QckJCQgMjIS+fn5bNlHjx5h7969mDBhArtOpVLB2toa06dPR58+ffTWIRKJMHr0aKxYsaLK9gwYMEDrPSQrKwtbt26t+wutgGEYKJVKkx7THOzs7NCiRYuGbgYhpBosKtTFxsbiiy++wGuvvVaj/Vq2bAl3d3d24fF49dRCQogpTJs2DUKhEIcOHUJ0dDS8vb0RGxuLI0eOICMjA/Pnz2fL/vbbbwgJCUHr1q3Zdba2tli9ejUmTpwId3d3g/UMGjQIu3fvRmlpqdH2iEQirfcQd3d3ODk5AdCcERAKhTh16hRbfunSpXBxcUFWVhbGjx+PEydOYPny5WwvX2pqKnsm4eDBgwgPD4dIJMKpU6fw4MEDDB48GG5ubrCzs0NERASOHDmi1R5fX198/vnnGD16NOzs7ODp6akTTtPS0jB48GDY2dnBwcEBw4cPx+PHjw2+xkuXLqFv375wcXGBWCxGdHQ0rl69qlUnAAwdOhQcDod9Xvn0q1qtxsKFC9G6dWuIRCKEhobiwIED7PbU1FRwOBz8/vvv6N27N2xsbBASEoJz584Z/TcghNSdRYW62goLC4OHhwdiYmJw7Ngxo2VlMhmkUqnWQggxn/z8fBw8eBBTp06FtbW11jZ3d3eMGTMG27ZtA8MwADTDMsLDw2tVV3h4OBQKBS5evFjr9pafWh07diwkEgmuX7+O+fPnY926dfDw8MDy5csRGRmJiRMnsr18Xl5e7P7z5s3DokWLcOfOHQQHB6OoqAhxcXE4cuQIrl27hv79+2PgwIFIS0vTqvfrr79GcHAwrl69ivj4eMyaNQuHDx8GoOn1GzJkCPLz83HixAkcPnwYDx48wIgRIwy+jsLCQowbNw6nTp3C+fPn0b59e8TFxbHjFC9dugQA+Pnnn5GVlcU+r2z58uVYunQplixZghs3bqB///4YNGgQ7t27p1Vu/vz5mDNnDhITE9GhQweMGjWqUfZUEtKoMBYKALNr1y6jZZKSkpi1a9cyV65cYc6ePcu89957DIfDYU6cOGFwn08//ZQBoLNIJBITvwJCiD7nz583+v/7m2++YQAwjx8/ZhiGYUJCQpiFCxcaPN64ceOYwYMHG9zu5OTEbNiwwej+PB6PsbW11Voq1imTyZiwsDBm+PDhTFBQEPPuu+9qHSM6OpqZMWOG1rpjx44xAJjdu3cbrLtcYGAgs2LFCva5j48PM2DAAK0yI0aMYGJjYxmGYZhDhw4xPB6PSUtLY7ffunWLAcBcvHiRYRjNe11ISIjBOpVKJWNvb8/s27ePXafv36XycTw9PZkvv/xSq0xERAQzdepUhmEYJiUlhQHArF+/Xqdtd+7cMfJbIIQwDMNIJJJa5xJ+QwRJU/H394e/vz/7PDIyEunp6ViyZAl69uypd5/4+HjMnj2bfS6VSrW+VRNCGhbzrIeOw+EAAEpLS2FlZVXr41lbW6OkpMRomd69e2P16tVa65ydndnHQqEQmzdvRnBwMHx8fLBs2bJq11+5l7G4uBgLFizAH3/8gczMTCiVSpSWlur01EVGRuo8L6/3zp078PLy0nrvCgwMhKOjI+7cuYOIiAidduTk5OD//u//cPToUTx+/BgqlQolJSU69RojlUqRmZmJqKgorfVRUVG4fv261rrg4GD2sYeHB9uGgICAatdHCKmZRh3q9OnevTs2b95scLtIJIJIJDJjiwghFfn5+YHD4eD27dt6r1pNSkqCk5MTXFxcAAAuLi4oKCiodX35+flwdXU1WsbW1hZ+fn5Gy5w9e5Y9Xn5+PmxtbatVf+Vyc+fOxcGDB7FkyRL4+fnB2toaw4YNg1wur/JY5UGXYRj2cUWG1gOaq3xzc3OxbNky+Pj4QCQSITIyslr1GmqHsXoFAoFOebVaXeO6CCHV1yTG1FV07do19lshIcTytGjRAn379sWqVat0LmDIzs7Gli1bMGLECDYIhIWF4fbt27Wq68GDBygrK0NYWFid2vzgwQPMmjUL69atQ/fu3fHWW29pBRShUAiVSlWtY506dQrjx4/H0KFD0blzZ7i7uyM1NVWn3Pnz53Wel/dyBQYGIi0tDenp6ez227dvQyKRoGPHjgbrnT59OuLi4hAUFASRSIS8vDytMgKBwOjrcHBwgKenp840MWfPnjVYLyHEfCyqp66oqAj3799nn6ekpCAxMRHOzs7w9vZGfHw8MjIysHHjRgDAsmXL4Ovri6CgIMjlcmzevBk7d+7Ezp07G+olEEKqYeXKlejRowf69++PL774Am3atMGtW7cwd+5ctGrVCl9++SVbtn///nj33XehUqm0rmy/ffs25HI58vPzUVhYiMTERADQulLz1KlTaNu2Ldq1a2e0PTKZDNnZ2Vrr+Hw+XFxcoFKpMHbsWPTr1w8TJkxAbGwsOnfujKVLl2Lu3LkANFeOXrhwAampqbCzs9M6dVuZn58ffv/9dwwcOBAcDgeffPKJ3h6sM2fOYPHixRgyZAgOHz6M7du3Y//+/QCAPn36IDg4GGPGjMGyZcugVCoxdepUREdHG7yoxM/PD5s2bUJ4eDikUinmzp2rc6GKr68vEhISEBUVBZFIxF4BXNHcuXPx6aefol27dggNDcXPP/+MxMREbNmyxejvmBBiBiYe31cn5QOLKy/jxo1jGEYzoDk6Opot/9///pdp164dY2VlxTg5OTEvvvgis3///hrVWZcBiYSQ2ktNTWXGjx/PuLu7MwKBgPHy8mI++OADJi8vT6ucUqlkWrVqxRw4cEBrvY+Pj973i4r69evHLFq0yGg7xo0bp/c4/v7+DMMwzIIFCxgPDw+tdu3evZsRCoXMtWvXGIZhmOTkZKZ79+6MtbU1A4BJSUlh388KCgq06ktJSWF69+7NWFtbM15eXszKlSt1LrTw8fFhFixYwAwfPpyxsbFh3NzcmGXLlmkd559//mEGDRrE2NraMvb29swbb7zBZGdns9srX+Bw9epVJjw8nBGJREz79u2Z7du3Mz4+Psy3337Lltm7dy/j5+fH8Pl8xsfHR+9xVCoVs2DBAqZVq1aMQCBgQkJCmL/++kvr9QFgfzcMwzAFBQUMAObYsWNG/y0IIXXLJRyGeTYquZmSSqUQi8WQSCRwcHBo6OYQQvRYtWoV9uzZg4MHD1Z7n5s3byImJgZ3796FWCyux9aZnq+vL2bOnNkk7lJBCKmZuuQSizr9Sggh+kyaNAkFBQUoLCxkbxVWlczMTGzcuLHRBTpCCKktCnWEEIvH5/O17jJRHf369aun1hBCiGWiUEcIIRZG39WwhJCGk5OWC0leoc56sYs9WnobnzLJnCjUEUIIIYQYkJOWi/EBM6AoU+hsE1gJsCFpucUEOwp1hBBCCGl0zNV7Jskr1BvoAEBRpoAkr5BCHSGEEEKaHnOErcbUe2ZOFOoIIYQQYhLmClvm7j1jRAJAqCcyyZUmq8MUKNQRQgghzYA5etDMGbbMFbTyC4qhjuwE8PTcWVWlRn5BsUnrqwsKdYQQQkgDamqnK80RtswZtIpK5PrrAQAeV7PdQlCoI4QQQhqIOU9XyhkA9tY62+Rypcl60MwVtswZtOwcbeq03Zwo1BFCCCGVmPPKSgpbls3J3QkqIaDmc3S2cZUMnNydGqBV+lGoI4QQ0miY61TluOA5UOi5NbqAw8H/biwxWV1NMWwZC0CmYudoA/iUgeus0tmmzueZtPcst6gYRVEqCG11Tx/Li/nILSpGB5PVVjcU6gghhDQK5gpbKfeyIevirzcEyVRqpNzLNlmoKyqRQ+bCgVJPBuGXcEzbs2UkBJmKTMQxGoBkIt2wVxv/iHLxwvt/g8fV/VtQqTn4R/4yuqGjSeq6V5yKvnFXDdZ1r6g7otDOJHXVlUWFupMnT+Lrr7/GlStXkJWVhV27dmHIkCFG9zlx4gRmz56NW7duwdPTE/PmzcOUKVPM02BCCCEAgNvXUpD1KF9nvUdrZwSGtTFJHSn3siEb6guuq1pnmyyXa7KwVVQiB9rKDQYgUwath2VPIR1eBpFANwSVKPh4WPYUMSao564y3WgIupsSihh0rXM9OTYFRgNQTlH3OtcBAI/lT8Cz0t/zx+MyeCx/Ur0DqdWAUml0UZRkgOdguC6FqKS2L8PkLCrUFRcXIyQkBBMmTMDrr79eZfmUlBTExcVh4sSJ2Lx5M86cOYOpU6fC1dW1WvsTQkhT9uvJI7j/5JHOer8WrTGyZx+T1XP7WgreX7EBXBfdsKXO42LlB+NNEuyuSh/ghRk3DQaGqylh6IbOda7HXAEIAHI4uXjFz/Bryrlvmt6mIqtivXUAmmBSJCgEpNIqA05VS2HBXfDCDNdTmHgSOHulzvVYdXIB3jT8eq02bQDe+rDqY+np9a3MbmgE8Flrw9vtrKo8hrlYVKiLjY1FbGxstcv/8MMP8Pb2xrJlywAAHTt2xOXLl7FkyRIKdYQ0IWce/Y20wmyd9d727ohqXfcPcXPXs//WBTwsyNRZ39bJE68EdTNJHb+ePIJEh9XgOep+aCWqOcBJmCzYnUi9jhemGw4mJ1KumyTUPcFTo8HkCZ7W/KB6emqelGUbr+fJQ+DsWa191HI5VHI5GIUCjFJp8CdHpYIdj8fuJxQUgdfBcF2qy0fxeMs2MCoVrNRqOHK5WvVmFhdD9ewxo1KxP5lK6/ivvADMEhr8NZRt+R+SB77FPq/YIj9oh4VcADmVypX/LHzjJSDMxXA9J48hcfspnToqP3cD0KrCczWAqxWeP7WKAGAkaOVmASkpAIACAMkGSz7HAcD+7+PxAD4fD7lcpAkN/94sjUWFupo6d+4c+vXrp7Wuf//++PHHH6FQKCAQCHT2kclkkMlk7HOpVFrv7STEnJpaADrz6G/sy/tM74fs9TwOgM9MUp+56tl/6wJOKRbrPZ2ToeAAt+YZDXYMw0CpVEKhUOhdlEoluFwu7j95pDfQAZqwcPzGBfCy86FSKqGUyaBSKKCSyxHQrh16dOmiFRyWrl0LuUwGlVIJlUIBpUKheaxUQqlQINXbCQHtDNf155GdOPnlF1CpVFCrVFCr1ZqlwuO9Q4aAXyFg/fD339h87x7UDAO1Wg2VWg2bgT3x4ocGBvsDwKHDwKff4NWUFNwsKwPDMFA/+52V/2QAqJ/9nA9gZoXdcwG0AdB+ZBzigkQGq/ktIQFxEz9ElwrrtgB4y9AOFbg8q6ec1cgY4BUHg+V/+echPt2dAAAYA2Bzpe0hAPKqUe+MYilsYThsbVLJ8X8Gtj12cEBLoRDg8wE+H+uLivDR06d6y7ZVlmGkkXYctOZjXjXa+6+XXsJXr7zC1qlgGETMmvW8HsBoPZg8CfjXCoDPx6nTpzF45swq6xQIBJCXlgJcLsDRjP1bOm0aDmRdMl6XBWnUoS47Oxtubm5a69zc3KBUKpGXlwcPDw+dfRYtWoQFCxaYq4mklpKfPsTjslyd9W5WrvB3bGuyeswVTMxVV1MLQACQVmi81yStMBtRtT3lxjDAs56M22n3wLMzXM/tvy8jKk/Fho5TV65AIpVCVlamWWQylD37yS5yOQZ26oTuPj7sfhdLMsF71XA9Xy7+Ah/sOgGFWo3TAwbAx8qK3Xfl3bv44O+/q3xZfgIBxr01BJhuuMzZM2ex5tePddZPBdCj0rp4APrvEaDReWQcAmA4BP2Tn4+/r1wxcgRA9c03Wh9I/wA4U7keuQIwUg+ePgX+/huZz/avir6RUMUAFNAdS1eRGiowrVsDdnZs6OBIJGzPkDGMSASMGsXuZ+UjApBmuC5VhfF7L7wATJrE7gs+X/O8Gp0Twpd7ArhtcLtKWWZ453v3gJYtnz9ftAj46KMq69TLpWXVZQCgRw/gX/96/lwmAyqEuip5+wAdn305ysio/n483YtGVAqZnoLP2Yh0O5AaSqMOdQDA4WhfScM8Oz9eeX25+Ph4zJ49m30ulUrh5eVVfw00k6YUgpKfPsTPqXMNnsqZ4Pu1SV6TOYOJueqqcQCqEGpqsqTl3wDP3Ug9Zw4jSnG9zuNmoFSiyNcaeM3wa87/4UdcPTobxXI5ShQKFCuVKFEoUKJSaR4/+6lUq7FYJNI69lcqFbZB80EuHtkf/eINT4PwdOM24Ne32efjATw03CyWxx9/oOLQcPnIOHBfNRxMiuRqpBRqpuyQ7dypta261ygqFAqguBhGA5AB7HB9gYANDrzCQqOhzhTUM2YA1tZsndxjx4ATJ9jtXA4HDEd3zF5FophewLuz4fLvf8Pt4UNwuVxwOBxwOBzNYy6XXcflcuE4ZQowcSJbp6CoCEExMRA720L35OBzLVqKYXPoENDx+Xi3locO4cXPP2frM7SIxWLg55/Z/ezuHAYUqw3W5de+DQIHOwMAukZHA++8o7V90KFDkEql7GdexZ8VHzv5+KDISKgL7RKGF+07an12lj+2stIeMxYSEoJ3KrSjYj1FvjwAjw3W07ljANpXuJix8md1+fPu3bUvqODxeHj//ffZ54XOHKjUGQY/J7zt3dnnbdu2xYwZM3TKVa6bpyfQvfzyy8BtAVTqdIN1BTh766xvKI061Lm7uyM7Wztg5OTkgM/no0WLFnr3EYlEEIlq/kZXG+bqBWpqISgpP81oMEnKTzPJ66kyAD1KRpTCrmYhxEBASuM/Bq+Lkbp2/A9Rj1D3ENStFTDZyEf/558De4dqt7cWit54CfjY8Kmcop27gWfjZiqSA5BA88ZTcbpOJYBVz7ZVXuQjY9ADhk9PPczJw/s3b1bZZi6A/5aVoeLbeBaAxGePO1cRmWQiG8DDgw0AoowMQF71VZBlgYGAv//zsNLKDsZOmAm5XLja2kLA5QKTJ2vV6XnzJrodPAgBn/98EQjYn/xnP1u2aAHEhEO3r+u54I4BGPfZC+AJBOALheAJBOAJBAjs1Ano2VOr7K979gDQfOjx+XzweDx24fP52HHnAoDTBut6JS4Oexd9Dy6XCx6PB+6zcFVxsXJ0ZE95AcDCjz/GQoANYQCw7tI+pOJn/ZUAcAwJACJicOjSJYNljHG0tcXNmzex9c5h3DQStKZMeRcdO2pfwNCvXz+d4UDV4WZrj5tPDW+fM/k9vOxt+IrRH3/8sVr1aN6/Dxj8nJg9/v1qv3/HxcUhLi5O77aqPo+mDHsL/u/W/P2bz+djxYoVWuuq+xnbqVMndux9Tb3++ut4/fXXzdZxUleNOtRFRkZi3759WusOHTqE8PBwvePpzMmcvUB1CkEMU/3emcfXwXM2EkwuHUfUpft1DiYlngxg5HqZknU/AlcX1T0ADQwD/s/TcEWLFwO7avfhoGNoBNDF8KBenDxpkrqK2toARsbNFCnxrAenClyu9imeCkspl4s8R/1fmsrtthPjN2dnSFQqSJRKdil7FiJnBwVhaY8ez4MOj4eZy5bp7RsJquJtignpAvyaUOVLUgOQ3bgBK1tbtl7b//4Xwh9+gI21NaytrGCsdwaDBwEbfmefvr9qFfLz8yESiWBlZcV+Yay8dOjQAahwNsBq588owz59NQAAhgx7HR9v0b998LOlOr7YtcHo9oDOnfHh0PHVOtbgwcZrvVyWwQ6c18fd0xW+vr7Vqqucvl6TQA9fPMjjGAwMgR41q8MQb3t3XDdST8VeoLrq4t4Ox57yoX7eR8rigo8u7qaZ/0zzefNZvXc0+Du2xQTfr80SgKJad6790Isa8ndsC39YTngzxKJCXVFREe7fv88+T0lJQWJiIpydneHt7Y34+HhkZGRg48aNAIApU6Zg5cqVmD17NiZOnIhz587hxx9/xNatWxvqJbCq7AX6Yxuicneb5PRUSZeWwATDbSlZ+AXw1wj9+6uNn87QUsVl3di0yTQhaGgEEGuknuRk4JAJ6qnqtfN5gI2NwXBTk6WooysA3Tm8yhVFvgS0713nesoK7gE4brCesrcmAP/9H1ueeXaFF6fCqba/b9/GXwcPIjc3V+9SXFyMoMIgDIbhU5X/lCpwK9/w65V07w6sXcs+5wKw/+knvRcucavoQRO28sbUqVNhY2MDW1tb2NjYaD2uuE4QGKg1ZubLZcvwn+XLAWhCUCn2Gq2roqlTp1a7bEUcrmkmX62KX4vWSFQbDiZ+LYz8H6uhcP+O2GckBIX7m2ZaDnMFE3PVAwCOQld8GPA9ipW6f/u2fAc4Ck13SzJzhaDGEoCaIosKdZcvX0bv3r3Z5+Vj38aNG4cNGzYgKysLaWnPB5S2adMGf/75J2bNmoXvv/8enp6e+O677xrHdCaHDpmuF6iF8Uu7IZUCOca+R+vx7MNea7G3M76PlxcQJaxzMIG7EsBdw/W8/jrw2vTqH0/fa+HzUfTgDIDfDFZTNPcj4LeBNfu9GVB2eheATYa3R74IvDi07hWdNjLYGcDhS5fwx4pftULatWvXtE4lXbx8Gf+qOEBZj6qCVvl2DocDBwcHiMViraVzZ90Plg0bNkAgEGiVc3BwwC9HjyFdvcVgYPCxa433vn9fZ1t1VBxT0wKO+MdICGoBx1rVUVkrkQuSjdTTSmS4p7UmRvbsA5yEWeapM2cIMlcwMWcvkKPQ1aThjTRfFhXqevXqxV7ooM+GDRt01kVHR+Pq1au6hS1deDjQsotJeoHAyQBwwnBd094HFtagLh5Pa2wLq4oBvZgyBejYt86/mqJL+2As1BV1CAQi6h62yv4xPrayTKZ7OqShMQyDvLw8pKenay29e/fGgAEDYFdqBZWd4cDw99HrSD19Vmt9bm6uVqhzddX9cOFwOHB2doarqytcXV1hZ+0OlVr/nGEqNQdvDhiJKWu2ws7ODlyukSkoKhg6VH+oDRD6YOfKzgZn9+872DQXOnVxaIftyzvpvVuBOpeLoW+a5jTYK5264dexu6Fupft74WaoMXuTaeapA0w3D111mDMEEUL0s6hQ16wMHWqSAATgWdgyEuo8PYCOnepcTVFRmdGL6YqKjPcSVZdAZgOVwHAwEchMc6NmN2ELZBjpMXETGh83VhNVhS27Uv0zkv/www84d+4cG94ePXqEsjLd37NKpcKAAQPgXdYSm4wEoNLrtwBoxiu5uLjoDXAvvPACfvvtNzbAubq6okWLFlpjnBL2XcGXK/cZrGf+4EA4OBi+uKEm7GyEwEMh1P/ov+G5nY1pJgZt094dol2peu8rKuJw0GaBacZRtfR2xaZNi+r9pvSEkOaHQl0TwC0VQcU1HBi4paa52lfMiKEyEoLEjNgk9bS39cWKP7sYvCH0wFhfk9QT6tQBP/zXcAB6+18dTFIPgCrD1qHTX2OV/f/hUqWr9g4cOIA9z648NCY9Pf35k3+soDYwSdeSr5cg7o0ecHR0NNiD5u7ujjfeeMNofeYKWgDg7GQL7rmbgFDP25VcCWcn00wL2tLbFf+7scQsYaultyuFN0KIyVGoqycOagejAchBbZpeDABozWmFxUZDUCs9e9WcuUKQq50t7M7woObr/nkKlQxc37A1ST0F2QVGA1BBdgHQvu69M3/99RcWzF8OG7deBuvKzMxEoSodDMNojfOqOIeivb09vLy89C7t27cH8CxsqdQAT3/Yau/nA2dn5zq/JnMFLUATqIQcQFFYqrNNYCWA2MXeZHVR2CKENGYU6upJe6EPln5rLAD5mKwuc4UgAEZDkKkUZBeAJwd4cv3jK00Vtoqe6ptPvnrbGYbBo0ePkJycjKSkJCQnJ7PLtm3bdCbOvH/3NoJdehoMWxwo0apVKxQXF8PO7vkFKbNnz8bkyZPh5eWlmbS0CuYKW+YOWhuSltPpSkIIqQKFuvpkhgBUzlgIamzqErZqoqperYqnEBmGwRdffIHbt2+z4a2kRH87kpKStEJdQEAAFPJio2HrxOmD6NBVdyB+mzY1vwk6R6YAZPU7/7+5gxb1oBFCSNUo1JFqK8guqHq7CXrQahK26qImvVocDgc//fQTUlNTjR7T0dERxZUm9/X19cXZs2cws8cnBsOWodva1ZTYxR4CKwEUZbr10KlKQghp2ijUNQHmCluQK42GLchNMwVIfZ9CLCgowJUrV3B0z8ka9WoFBAQgNTUVPB4Pbdq0QUBAAPz9/eHv788+dnV11XsvQ1fPFmYJW3SqkhBCmi8KdU1BEwlbFdXXKcRvv/2WndTaHo7oxjE8j1flORMXL16Mb775Bu3atYNQWLPeQnOGLepBI4SQ5olCXT1RFpVqbkOlb9oItVqz3USaQtgyheLiYiQmJuLy5cvssmPHDgQFBbFlyq8SBQA5ZFAxKvA4undHEIj4ELtoX6Gs7y4INUFhixBCSH2iUFdPBAwD7lnDQUtg5M4ZtWHJYaumqjMurLS0FDdu3NAKcLdv34a60v1cL1++rBXqwsPDERkZifDwcISHh6NtKz+4OLQEr1IvJ52qJIQQ0thQqKtHTSlomVPlU5UV524rD1utWrVCZmam0eMIhUJkZ2vfi9Ld3R1nz541sAchhBDSeFGoI9VmrisrZTIZbj+8hYSEBCQkJMDW1haHDx/WKhMSEqIV6vh8Pjp37sz2wIWHh6NTp041HvtGCCGENFYU6poAc4Wt+hrsr1KpcO3aNTbEnT59GqWlz8cc2tnZQaVSad17dNCgQXB3d2cDXHBwMKys9N8/lRBCCGkOLC7UrVq1Cl9//TWysrIQFBSEZcuW4aWXXtJb9vjx4+jdu7fO+jt37iAgIKC+m2oxGuuVlbdv38bHH3+M48ePo6DA8LQsnp6eyMnJgYeHB7tuypQpmDJliknaQQghhDQFFhXqtm3bhpkzZ2LVqlWIiorCmjVrEBsbi9u3b8Pb29vgfsnJyXBweH6loqtrww9wN+cksIDlX1mZkZEBAGjV6vl9aK2srLBr1y6dsq1bt0ZMTAxiYmLw8ssva+1DCCGEEP04TOXJuBpQt27d0KVLF6xevZpd17FjRwwZMgSLFi3SKV/eU1dQUABHR8da1SmVSiEWiyGRSLSCoSnkpOU220lgCwoKcPz4cRw5cgQJCQlITk7G3LlzsXjxYq1ybdq0gVQqRe/evdkg1759e5PdYYEQQghpTOqSSyymp04ul+PKlSv497//rbW+X79+VV6tGBYWhrKyMgQGBuLjjz/We0q2nEwmg0wmY59LpdK6NdwIS+89M6WSkhKcOXOGHRd39epVnelFEhISdPY7efIkWrVqBa6++fwIIYQQUm0WE+ry8vKgUqng5uamtd7NzU1nWopyHh4eWLt2Lbp27QqZTIZNmzYhJiYGx48fR8+ePfXus2jRIixYsMDk7W/OVqxYgTlz5kAul+vdzuPx0K1bN/Tt21drehIA8PLyMlczCSGEkCbNYkJducqn3SqHgIrK77tZLjIyEunp6ViyZInBUBcfH8/eJgrQ9NRRsKgemUyGw4cPIzo6Gvb2z8cE+vj46AS6zp07IyYmBn369EHPnj21yhNCCCHE9Cwm1Lm4uIDH4+n0yuXk5Oj03hnTvXt3bN682eB2kUgEkUhU63Y2N6WlpTh06BC2b9+Offv2QSqVYuvWrRg58vmtx6Kjo+Hn54devXqxFze0bNmyAVtNCCGEND8WE+qEQiG6du2Kw4cPY+jQoez6w4cPY/DgwdU+zrVr17SmviA1V1JSgr/++gs7duzAH3/8gaKiIq3tO3bs0Ap1YrEY9+7dM3czCSGEEFKBxYQ6AJg9ezbGjh3L3p9z7dq1SEtLY+cji4+PR0ZGBjZu3AgAWLZsGXx9fREUFAS5XI7Nmzdj586d2LlzZ0O+jEbr8OHDWLduHfbv34+SkhKd7WKxGIMHD8aoUaMaoHWEEEIIMcaiQt2IESPw5MkTLFy4EFlZWejUqRP+/PNP+Pj4AACysrKQlpbGlpfL5ZgzZw4yMjJgbW2NoKAg7N+/H3FxcQ31Ehq18+fPY/v27VrrnJycMGTIELzxxhuIiYmh224RQgghFsqi5qlrCPU5T50lkkgk2Lt3L3bs2IGvvvoKHTt2ZLfduXMHgYGBcHFxwdChQzFs2DD07t0bAoGgAVtMCCGENB9NYp46Un8KCgqwZ88e7NixA4cOHYJCobnLRXh4OD755BO2XMeOHXHq1Cl0794dfD79aRBCCCGNCX1yN1FPnjzB7t27sWPHDhw5cgRKpVKnzIkTJ7RCHQC8+OKL5moiIYQQQkyIQl0T9Pnnn2PBggVQqVQ621q1aoVhw4Zh2LBh6NGjRwO0jhBCCCH1gUJdI1c+JLLiBM3t27fXCnTe3t5skOvWrRvdkosQQghpgijUNVISiQRbtmzBmjVr8M033yAmJobd9uqrryIwMBCvvPIKhg0bhoiICIN35SCEEEJI00ChrhFhGAaXL1/GmjVrsHXrVnYuuTVr1miFOjs7O9y6dauhmkkIIYSQBkChrhEoLCzEL7/8gjVr1uDatWs627Ozs6FSqcDj8RqgdYQQQgixBBTqLNj9+/exePFi/PLLLyguLtbaZm9vjzfffBOTJ09GSEhIA7WQEEIIIZaCQp0Fy8/Px7p167TWRUREYPLkyRg5ciRsbW0bqGWEEEIIsTQU6izEtWvXUFhYiJ49e7LrIiIiEBoaivv372PMmDGYPHkywsLCGrCVhBBCCLFUFOoaUHFxMX799VesWbMGly5dQmhoKK5evcpeqcrhcPDLL7+gdevWsLe3b+DWEkIIIcSSUahrADdu3MCaNWuwefNmSKVSdn1iYiIuXbqEF154gV1X8d6shBBCCCGGUKgzk5KSEvz2229Ys2YNzp8/r7M9NDQUkydPphBHCCGEkFqhUGcGpaWl8PX1RW5urtZ6GxsbjBo1CpMmTaIJggkhhBBSJ3S/KDOwtrZGdHQ0+zw4OBjff/89MjMzsX79erzwwgsU6AghhBBSJ9RTZybTpk2DnZ0dJk+ejG7dulGII4QQQohJcZjyO8I3UwUFBXB2dsatW7fg4ODQ0M0hhBBCSDMmlUoRFBSE/Px8ODk51WjfZt9Td//+fQBAUFBQA7eEEEIIIUTj/v37iIiIqNE+1FP3rKcuPT2deuoIIYQQ0qCkUim8vLyop642eDweAMDBwYFCHSGEEEIsQnk+qYlmH+oIIYQQYjqPcySQSEp11ovF1nBrKW6AFjUfFOpIs0ZvPoQQYjqPcyQYN34t5AqVzjahgIf/bZhE7631iEIdabbozYcQQkxLIinV+54KAHKFChJJKb2v1iOafJg0W9V58yGEEEIaCwp1hBBCCCFNAIU6QgghhJAmgEIdIYQQQkgTQKGOEEIIISYhFltDKNA/v5pQwINYbG3mFjUvdPUrabbK33wMXf1Kbz6EEFIzbi3F+N+GSTRVVANp9rcJk0qlEIvFkEgkdEeJZojmqSOEEGJJ6pJLqKeONGtuLcUU3gghhDQJNKaOEEIIIaQJoJ46QsyETvUSQhoSvQc1fRTqCDEDuiUZIaQh0XtQ80CnXwkxA7olGSGkIdF7UPNAoY4QYvGSk5Ph7u6OwsLCau+zcuVKDBo0qB5bZdk+++wzhIaGWsxxCCH1j0IdIaRBpKen45133oGnpyeEQiF8fHwwY8YMPHnyRKfs/PnzMW3aNNjb2wMAjh8/jsGDB8PDwwO2trYIDQ3Fli1btPaZOHEiLl26hNOnTxttx/jx48HhcHSWAQMGVPu19OrVCzNnzqx2eUvF4XCwe/durXVz5sxBQkJCwzSIEFIjFhXqTp48iYEDB8LT01Pvm0tlx48f1/tmnJSUZJ4GE0Jq5eHDhwgPD8fdu3exdetW3L9/Hz/88AMSEhIQGRmJ/Px8tuyjR4+wd+9eTJgwgV139uxZBAcHY+fOnbhx4wbefvttvPXWW9i3bx9bRiQSYfTo0VixYkWV7RkwYACysrK0lq1bt5r0NTMMA6VSadJjmoOdnR1atGjR0M0ghFSDRYW64uJihISEYOXKlTXaLzk5WevNuH379vXUQkKIKUybNg1CoRCHDh1CdHQ0vL29ERsbiyNHjiAjIwPz589ny/72228ICQlB69at2XUfffQRPv/8c/To0QPt2rXD9OnTMWDAAOzatUurnkGDBmH37t0oLTU+XkgkEsHd3V1rcXJyAqD58igUCnHq1Cm2/NKlS+Hi4oKsrCyMHz8eJ06cwPLly9kvlqmpqeyXzoMHDyI8PBwikQinTp3CgwcPMHjwYLi5ucHOzg4RERE4cuSIVnt8fX3x+eefY/To0bCzs4Onp6dOOE1LS8PgwYNhZ2cHBwcHDB8+HI8fPzb4Gi9duoS+ffvCxcUFYrEY0dHRuHr1qladADB06FBwOBz2eeXTr2q1GgsXLkTr1q0hEokQGhqKAwcOsNtTU1PB4XDw+++/o3fv3rCxsUFISAjOnTtn9N+AEFJ3FhXqYmNj8cUXX+C1116r0X4tW7bUejPm8fTfd46QhkL3Q3wuPz8fBw8exNSpU2Ftrf263d3dMWbMGGzbtg3lN7s5efIkwsPDqzyuRCKBs7Oz1rrw8HAoFApcvHix1u0tP7U6duxYSCQSXL9+HfPnz8e6devg4eGB5cuXIzIyEhMnTmS/WHp5ebH7z5s3D4sWLcKdO3cQHByMoqIixMXF4ciRI7h27Rr69++PgQMHIi0tTaver7/+GsHBwbh69Sri4+Mxa9YsHD58GICm12/IkCHIz8/HiRMncPjwYTx48AAjRoww+DoKCwsxbtw4nDp1CufPn0f79u0RFxfHjlO8dOkSAODnn39GVlYW+7yy5cuXY+nSpViyZAlu3LiB/v37Y9CgQbh3755Wufnz52POnDlITExEhw4dMGrUqEbZU9lU0HtQM8FYKADMrl27jJY5duwYA4Dx9fVl3N3dmZdffpk5evSo0X3KysoYiUTCLunp6QwARiKRmLD1GhkFEubWo2ydJaPA9HURy5f9+CmTfDdLZ8l+/LShm2ZW58+fN/r/+5tvvmEAMI8fP2YYhmFCQkKYhQsXGj3m9u3bGaFQyNy8eVNnm5OTE7NhwwaD+44bN47h8XiMra2t1lKxTplMxoSFhTHDhw9ngoKCmHfffVfrGNHR0cyMGTO01pW/P+3evdto2xmGYQIDA5kVK1awz318fJgBAwZolRkxYgQTGxvLMAzDHDp0iOHxeExaWhq7/datWwwA5uLFiwzDMMynn37KhISEGKxTqVQy9vb2zL59+9h1+v5dKh/H09OT+fLLL7XKREREMFOnTmUYhmFSUlIYAMz69et12nbnzh0jvwVS3+g9qHGQSCS1ziWNep46Dw8PrF27Fl27doVMJsOmTZsQExOD48ePo2fPnnr3WbRoERYsWFDvbct8KsUrSzdArtQzJxCfh/0fjoenI91rtjmhW5JVD/Osh47D4QAASktLYWVlZbD88ePHMX78eKxbtw5BQUE6262trVFSUmK0zt69e2P16tVa6yr2+gmFQmzevBnBwcHw8fHBsmXLqvtydHoZi4uLsWDBAvzxxx/IzMyEUqlEaWmpTk9dZGSkzvPyeu/cuQMvLy+tHsHAwEA4Ojrizp07iIiI0GlHTk4O/u///g9Hjx7F48ePoVKpUFJSolOvMVKpFJmZmYiKitJaHxUVhevXr2utCw4OZh97eHiwbQgICKh2fcS06D2o6WvUoc7f3x/+/v7s88jISKSnp2PJkiUGQ118fDxmz57NPpdKpVpvjKbytLhUb6ADALlShafFpRTqSL2w9Fnj/fz8wOFwcPv2bQwZMkRne1JSEpycnODi4gIAcHFxQUFBgd5jnThxAgMHDsQ333yDt956S2+Z/Px8uLq6Gm2Tra0t/Pz8jJY5e/Yse7z8/HzY2toaLV/x2BXNnTsXBw8exJIlS+Dn5wdra2sMGzYMcrm8ymOVB12GYdjHFRlaD2iu8s3NzcWyZcvg4+MDkUiEyMjIatVrqB3G6hUIBDrl1Wp1jetqDiz9/yxpPBp1qNOne/fu2Lx5s8HtIpEIIpHIjC0ixHwaw6zxLVq0QN++fbFq1SrMmjVLa1xddnY2tmzZgrfeeosNAmFhYbh9+7bOcY4fP45XX30V//3vfzFp0iS9dT148ABlZWUICwurU5sfPHiAWbNmYd26dfjtt9/w1ltvISEhAVyuZliyUCiESqX/S1xlp06dwvjx4zF06FAAQFFREVJTU3XKnT9/Xud5eS9XYGAg0tLSkJ6ezn4pvX37NiQSCTp27Giw3lWrViEuLg6AZkqZvLw8rTICgcDo63BwcICnpydOnz6t9cX57NmzeOGFF6p45USfxvB/ljQeTS7UXbt2je3qb04yn0rxtFj3m56jrTX1CDYj1Zk13hI+IFauXIkePXqgf//++OKLL9CmTRvcunULc+fORatWrfDll1+yZfv37493330XKpWKvQjq+PHjeOWVVzBjxgy8/vrryM7OBqAJVxVPm546dQpt27ZFu3btjLZHJpOxxyjH5/Ph4uIClUqFsWPHol+/fpgwYQJiY2PRuXNnLF26FHPnzgWguXL0woULSE1NhZ2dnc4FGxX5+fnh999/x8CBA8HhcPDJJ5/o7cE6c+YMFi9ejCFDhuDw4cPYvn079u/fDwDo06cPgoODMWbMGCxbtgxKpRJTp05FdHS0wYtK/Pz8sGnTJoSHh0MqlWLu3Lk6F6r4+voiISEBUVFREIlE7BXAFc2dOxeffvop2rVrh9DQUPz8889ITEzUmSeQVE9j+T/b3DWWz1iLCnVFRUW4f/8++zwlJQWJiYlwdnaGt7c34uPjkZGRgY0bNwIAli1bBl9fXwQFBUEul2Pz5s3YuXMndu7c2VAvoUHQ+D3S2LRv3x6XL1/GZ599hhEjRuDJkydwd3fHkCFD8Omnn2qFori4OAgEAhw5cgT9+/cHAGzYsAElJSVYtGgRFi1axJaNjo7G8ePH2edbt27FxIkTq2zPgQMHdL4M+vv7IykpCV9++SVSU1PZOfDc3d2xfv16DB8+HH379kVoaCjmzJmDcePGITAwEKWlpUhJSTFY17fffou3334bPXr0gIuLC/71r39BKpXqlPvwww9x5coVLFiwAPb29li6dCn7+svn8fzggw/Qs2dPcLlcDBgwwOicfD/99BMmTZqEsLAweHt74z//+Q/mzJmjVWbp0qWYPXs21q1bh1atWuntQZw+fTqkUik+/PBD5OTkIDAwEHv37qWppEiT1Zg+YzlM+ahkC3D8+HH07t1bZ/24ceOwYcMGjB8/np3/CQAWL16MtWvXIiMjA9bW1ggKCkJ8fDx7eqE6pFIpxGIxJBIJHBxM949yO+Mx3lj5i8Ht298fjcBWbo2uLmLZ7t7LxpSpGwxu/2HVeHRo726+BpnIqlWrsGfPHhw8eLDa+9y8eRMxMTG4e/cuxOLG1dPh6+uLmTNnNom7VBDjmur/2abE3J+xdcklFtVT16tXLxjLmBs2bNB6Pm/ePMybN6+eW1U7jrbWEPJ5BpO9oy3NCURIdU2aNAkFBQUoLCxkbxVWlczMTGzcuLHRBTpCSPU0llOi5mRRoa4p8XR0wP4Pxze5Pzj6T0QaAp/P17rLRHX069evnlpDCGlojemUqDlRqKtHno4OTeqPiv4TWb7yWeMNXUlHs8Y3DvrGspGmif7P1g5NG6YfhTpSbfSfyPK5tRTjfxsmmWXOK5pbizRl5vr7Nuf/WVJ7jEANfbdr4FjYne8o1NWjzGIJ8mW6M9k7i2zgaWu6/6hNcfweneatPXPMGk9za5GmzNx/303tTg9N7f1bxlVC1rEM4OrZqNZstxQU6upJZrEEff5cDZla901BxOXhSNx7Jgt2TW38Hp3mtXw0txZpyujvu/aa4vu3QMTVH+gAgPtsu4WwnJY0MfmyEr2BDgBkapXeHrw6ETBQW6t1FggsZsaaaqvOaV5CCCGWx5zv34xA9zNPba0GI2i+t6OjnromwJy9gk1NUztNQAghhjSl97vGdErUnCjUNQHV6RU0RahztLUG35oDOfR0q6Pxjd1riqcJCCFEH3O+35njooLGdErUnCjUkeoTMJB3LINcX4Dk8hrdqV5zX83blL4lmxtdaUvKNcW/BXO8N5jr/S6nrMhoD1pOWRECQXc4qi8U6ki15ctK9Ac6AHIT9ggCjefy8epqar2C5pxbi660JeXM9bdgzr/vptaDJlUaCHQAwH22vZFxFtlAxOUZHOLkLLJpgFbpR6GunjSmPwJLY86xEuYKj+bsFTTHt35zzq1FVyKScub6W3BrKcbi70YjI1eis62Vq2mnH3laXAoZRwFGT06UKdWNrgfN3lpUp+2WyNNWjCNx75llirK6olBXTxrTH4GlMddYiaZ4msCc3/qb2txahJTLfCrF2F+2Gxw//Nf0t032/8hc70Pm6kFztbOt03ZL5WkrbhSf2xTq6pG5/gioV7B2zH2awBy9gk2tR9DcmuJ4LXMx1+9O5sRAacPRWc8vMd2Y3rtPclHUvljv+4Ncrdluqr/xpni60hzoc08/iwp1J0+exNdff40rV64gKysLu3btwpAhQ4zuc+LECcyePRu3bt2Cp6cn5s2bhylTppinwRaCegVrx5ynCczZK2iO8Jj5VIrY734yS08GAKiEgJqv+0HOVZrug5zG7tXe4xwJxn7yHdRiuc42rkSITZ9PN8nv7k5BDqTDyyAS6P4xlyj4uFOQgw5wr3M9UmUZbIRyiPi69ciUfJMHLRu+4bpMpamdFqXPPf0sKtQVFxcjJCQEEyZMwOuvv15l+ZSUFMTFxWHixInYvHkzzpw5g6lTp8LV1bVa+zcl5ugVbGrfjMx5msBc38bNFR7N2ZORW1SM3DA+GKFuqOPIGeQWFaODCeqRSEpRaKfU2wskK1GadOze/lsX8LAgU2d9WydPvBLUzSR1AMCvJ4/g/pNHOuv9WrTGyJ59TFbPn9fPost7ieBxdUO2Ss3Bn9fPYkLf2DrXk1GWiVf8bhqsJ0MaCSC4zvXImKd4pZ3hemRMvzrXYe66XO1sjYZHU73fOYts4ChUgsPVDfiMWmjSz4nGckrUnCwq1MXGxiI2tvr/8X/44Qd4e3tj2bJlAICOHTvi8uXLWLJkSbMLdeZA34xqz1zfks0VHs15yiiXKYYsSGYwqOYyxSapx1y9QPtvXcApxWLwHHQ/xDMUHODWPJMEu19PHkGiw2rwHHXrSVRzgJMwWbDLKHoMnpv+XlMel0Hq0wxIpVKoVCqo1Wq9S+vWrcHhPA/UOTk5yM3NZberVCo8kabrfT3l9QisFQCA8+fPQyKRgGEYqNVqrZ8VHwcFBcHf3589RllZGXbu3Ikb3EzwgwzXc+7iCfR28YeLiwu7Pjk5GYcPH2brMLRYW1tj2rRp7H5CoQI8ueG6jh39A3d+SwDDMOjcuTNee+01rTJfffUVioqKAECrnorPASBkSJTR8Lj8h8Wwz3v++y/fDwAWLlwIe3t79vnhw4exb98+rXLlPzm2SvQfkQOOnnTB4wiwZcP3ePB3mtb6inWV69+/v9ZrVSgUmDp1qtY+XDsVuNaMzjEmjXsPXTt2Z59fu3YN3333nW6DKrePx8P69eu11m3cuBEJCQng2avBtWbg5+eHsWPHsttt+Q5wFLpWeWxzsahQV1Pnzp1Dv37a32L69++PH3/8EQqFAgKBQGcfmUwGmUzGPpdKpfXezqaEegRrp6kOHjYHri3P6Kkwri3PJPWYqxfowZNHRoPJ6b8vQZgphUKhQK9evWBjYwOo1YBSiVvXr+P06dNQyGRQlJVBIZdrHpcvcjmUCgWcbW2B8A7gRRqu58ctP2L9lOlQqlRQlS9qNYZ17oy5PXsCKhWgVAIqFQJWrkSZUgmVWg2VWg3ls58qhoFSrUb7N/oh7t+Gv5js+n0XFg6faPT3IuvVC0IAYBhArcY3Dx/ivxkZWmXaD43CG5+1NHyQ9T8C78zHe/fuIbGs6i8W/3FzQ3zL58eTKpV4884dtB0agZFBrQ3ut2f3brz51Wa42NsDHA7A5eJiXh4+uHevyjpdhEJMO3CA3c/J3xoYY7j8vp3b8ejoAwDAGH9/vHbqFMDnAzwewONh6bJlyCup+raTS60mgPe64b+HP/dtR+b5DL3b40eNgr2Hh6ZegQCXT5/GihUr9JZ1CxDjnTEv692mYhQ4cf4o9m9OqLK9Tk5OWqFOrVZrBS4Hd2tM2dsXfJHu///d8qVoJ1/Fhq309HRs2LChyjoFAoFOqLtw4QJ2H9rO1lWIG1h1fy67nc8RYJb/SosJdo061GVnZ8PNTftUkpubG5RKJfLy8uDh4aGzz6JFi7BgwQJzNZHUgrl6BJtieDTnuBlzjAMC6vf01KNHj5CZno7iggI8zP4bonDDH3qqtDvArlJAJgPkckz9+Wek5+VBJpdDplBAplCgTC6HTKl8vqhUWBwWhre9vNj9Hvu5AEbyzb49f+CrXz8AANwTCuGnVGpCHYBjAD6oxutqC2BCURwQafjf+7FUir/v3NFZ3zUtDdi/X2tdKgCZTsnnVCYY2qg+flzrub6OWR6MB/iy3ALg77+rfVNz5vFj4PFjnTpVCmOvVrOdqfS70z1pb6BOuRz44w/2ubh7K2DMCwbLq0sqhNPkZM1Si3rtf/sdT17XH7YAQPW0yPDO3bR7jqtbp17JSdUrd+oU8MUXgFCoWTjatVo7CvUGOgBQc1QofpwCR3sBYGWl+aJQB8bqUjIKFCulFOpMhVPpH5rtAubo/7OLj4/H7Nmz2edSqRReXl7110BSK+boEWyKp5PNNW7G5EGLYQCFAigpAUpLoSoshDQ3F5LcXDwqSgIv2HDYevT7FkzeNhXFZWUokclQLJejRC5HiVKJYoUCJSoVipVKqBkGBW3bAmVl7PJlcTF+ePae4T8yBq+HGxkHuH8/8OvzHoYjAKrulwEKz58Hzp9nn6tHxoELw2GLWyG4KOXa45KE1agPABR8PuDgAONRTD91+/ZAjx5sTxB4PLht2QKZSgUehwM+jwcelwsehwMelws+jwexlfEvC24ODvAIDASXywWPywWXywWXw9H8fLZwpk4FRCK2Byv84kW8c/WqVpk8XxcARsLHoEHA1E8wcd8+ZOXna44LzecBl8PR/Hy2jsvl4qVOnYDAQHZ3O7kc3/35JySeKpTiscFqpkW9CJ+R8zS/Y7UaYBhEZWVhU1ISOAyjqbP8Z/ny7LmIxwNCQjT7qdVwsZEAOGewrm+jomHvxwVUKnja2AAuLppe1GfL9rQ0KJRKcNRqQKXS/FSrwXm2lD928rNFqpF/o59btYK4RAyOUsn20JY/dmIYzTqFAmAYjAMQU2FfToWfBdB8+TDk48IifKFn38rPXc+eBc6eZdcLAPxdoVw+gING6kG/fkCSZp7B3gDuCASacCgSPV/KA2PF58OHa4KgSARYWeH/GAZvvvEG/sATY7VZjEYd6tzd3ZGdna21LicnB3w+Hy1atNC7j0gkgkjUuK7yIfXHnNPOmGPwcLEy32jYKlbmA9DtwTZKrdaEreJioKgIKC6GMPceeK6Gg5Zyz694+OUPkEilkBQWQlJUBElJiWYpK0M4j4cBPB4b4uTFxQhjGEgASKD9sd12aARGBhs+FVZ4+TLW3rhRvZdy/75WL07FiKtW6f7bVGTl6gxERbEfBKKzZ4HCQq0yIh5Ps/D57GLXpw/Qsye7Hz/vIdRaH0/aAtr4YvCUKRCIRGgxbhzg7s6e9opOScFPly9DIBJpFoFAa+Hz+RAIBLCxscEfKYkA9hqsZ+iQITi7bht4PB74fD54PB64XP19XP+sXm30d/P1H5vxFL8b3N5v4CuYu2ar0WNU9tqwYXit0rqtp/biJjYY3MelXXvgpT6Y0qd2YwWtAHwQF4ebD69ga9GXBsvFDn8bLm27aq1r82ypsZIHwH3DoS5q4UK0smlncHt0NavJKHkAVDhtWFnw3r1G62Gp1fBQKOChUDwPehUeZ8hScUz5vcHdvRYvRqtCO7bnml0qPtezjSuXo1OFbRluauPtrNCxYw8goLydxTUbg+sGQBkgBoYZ7uW0JI061EVGRrKDNcsdOnQI4eHhesfTEVLZU3kuipW64ypNPfjVRiDHK363oWIUOtt4HAFsBMYDRbXI5bAqLdAb6ABN2LL6YweQuUvzxlZ5eRbYKi6yoiLklpUhF9Ba0rq3gvUaw6eM5uzeZXB8DgC8D2BAhecCaHq9dH87VZ8Ks+8YAOy6pLNewOPBRiSCrZUVbJ4tsvXrYe3goPkmbmWFF48fh/zCBdjY2+NJGwGAuwbrsXtvItCxL/s8IScHHA4HIpEIVlZWEAgEBs8QVCTatQGlRkJdcEQEPh46Xu82/7Aw+IeFVVkHgGehzjCBUAA7O7tqHasqQa3a4KSaY/DLRFCrWsUdHQG2rXBdDvD0dFmq5JrtpqB4CqhVAFfP93+1DFCYcHogW74D+BwBlHreG/gcAWz5FjbnI5f7vGdLnxIA943sHxMDVCc8VqWKkIrLlwGBt1bPPLvIZDVbb1cIIKfubTYDiwp1RUVFuH//+V9DSkoKEhMT4ezsDG9vb8THxyMjIwMbN24EAEyZMgUrV67E7NmzMXHiRJw7dw4//vgjtm6t2TdC0jw9lefi2+T3Db6ZmnLwa7FSqjfQAZrBw8VJV+BYaK/p+Skq0vys6WO5HC4BYmCb4W+Udl9/i9QkiVZAGwKgYl/lNgAfPdtWqHMEDbenRXjHyOs1Oj4HgKRfP+DzzwEbG8DaGhwbG7Tq3h1ypRJisZhdHBwcYNfRHpqTOvpFTpuG6yPmwMbGBra2trCxsYGNjU21vtgNadcOQ97RvJKjaeeR8HSxwbJutvZaz1u2NDJg34hWIhckGwlArUQuevaqOb8WrZFopB6/FoZ7P2sqLiwKuAY8yNWdPqWdc2vNdhOwgxP+HsQF31F3m/Ip8M4fTiar58Zgw/W8baJ6AMBR6IpZ/ivr/culucKjRYVUgUCz2NtXXdaYqgKkBbGoUHf58mX07t2bfV4+9m3cuHHYsGEDsrKykJb2/FLoNm3a4M8//8SsWbPw/fffw9PTE9999x1NZ9IEmKMHrVgp1fvGAxgZ/MowQGkpIJEAUunzparnDsVAvLPhxowYyY7/qA01DM8wUlEUoDNS6OqcOQgLCABsbQFbW6jOn8fD//yn1m0BgKioKAhDHbQCWsWlTZs2QKdOWvuk/POP3mNllDzQutqsMi9HF7TyrPs3f2+1E9RygKunF0gt12w3hRfdA7FlmgM4Prqhk/lHgXdWBOrZq+ZG9uwDnIRZ5qkDYLLgZozYxR7MUyFKsnX/3wqsBBC71PHDuwJ5Ngfy7KrLmYKj0LXeB9qbKzw2tZBq7rrqyqJCXa9evfTOV1NO3yXJ0dHRuHr1aj22ilRkjrBVbz1oDKPp1Xr6VLMU3gWMnX2aMUMTtCoHNZX+23AZFSAGYGRMhrMz0M5F843Szk7z89njUmtrPGIYpMvlSJfJkF5SgnSpFOlPnyI9Px/pOTmYO3ky5n/0EcDLAVLja9S03L59NYOKn2lpbQ2n1avh6upqcLHyAK5ip8Fjfvfdd9Ubn1MN5npDZfIFuGGkF4j5QwD41r0esYs9hNdkUJzT7c0UmjiYmDq4NbSW3q7YkLQckjzdPmSxiz1aelvGFYg1lZOWa5bXZI7waK56zBUezV1XXVlUqCOWzVynK6vsQXuUDMfCzOfhrHyRSHTXVV6vrjC4topTlTh92nDvGYejufJNLNb8LF8qPy9f56oE8JfBqq4tXoxbxQ468y4OGjRIZ9yoPumFhYCTE1CSb7Rc165dYdPJWSugVZx8FQD69OmD/Hzjx3kqz8WN5L1m+eZqzjdUc/TONNVgYi4tvV2b1O8oJy0X4wNmQFGmv/dxQ9LyJvV6TclcIdXcddUFhbomwmJPV1amVmt6uwoKgPx87Z/lj/m5wEgjx3jl1TqdqoRAADg6At5V3B3g0081A231hTZbW515kwwpLi7GnSRjF/kD7777LgRF9khPT9da7+BgPBwJhUK0bt0ajo6OAKru1dry068m+Xsw9zfXxvKGWl1NLZg0NWIXewisBAaDlil7UyV5hXrrAQBFmQKSvEL6WyHVRqGuCTDngH+j9u4FMtT6A1v5z8q9ZfoEiIGRRnrQhELAzU0TzPQtYrHxbdbWmkBW1eDXQYOqfZWWSqVCWloakpOT0a1bNzg5PR9/tX37dsyIn2pw9nOlTIXSp3Lk5mRCqVSCz3/+3zIkJAQpKSnw8vLSu7i6umpNQ2HuUxJNKWgRUo56U0ljRaGuHplruow69aCpVJqg9eRJ1YtdIfCFkYmaFyyofg+atbXmdKGTk2Y8WcWfviIAyYb3PXfOJJfE12asllQqRXJyMpKSkpCcnMw+vnfvHnv7ub/++gsDBjyfsMPf3x/S7FL8MOgwrB2fj8J3cXGBm5sbXOzdMGHEJHh7e+uEurlz52Lu3JpdddUUw5a5xhwRUo56U0ljRKGunlhM7xkA/Pgj8LBUN6SV955V9xYqAWIARkJdVBTQtVJQ0xfanJw084QZYqbLxw31aqnVajBlXK1/H4Zh0KFDB60pdwxJSkrSCnUdO3bEO++8A39/fwQEBMDf3x++vr4QCqt7f4DmzVxjjsx5yo0QQuoDhbp6UufxZwyj6UHLywNyc5//rPi4/Kd9EfBtgOFjrVpVdQ+agwPQooVmcXZ+/rji4qYGjMwaj+++a7AetNpSSjm4eeU+Ll++jGvXriEpKQn3799Hv379sHfv89n4ORyOwYlaBQIB/Pz84O/vD39/f7zwgvakvI6Ojjo3iSbVZ64xR3TKjRDS2FGoayiHDwHpKsOhLS9Pc9uV6ggQAzAS6oYMAdRu+oNaeYirxkSttvJc8JP31XvYqs9xYbdu3cIff/yBy5cv4/Lly0hNTdVbLilJ96bTPXr0gL29PRveynvd2rRpo3XKlDRedMqNmJs5e4hpGEPTR59EDeXf8dUbf2Zvr7mBs6ur7s/yx25qABsNH+OTT0zSg9aYBuEXFhbi6tWrCA8Ph63t8zt8Hj9+HP/+978N7lfe69apUycwDKN126fvvzd8P0NCCKkNc/UQ09QpzQOFuobStSsQ5KQdzir/dHExPvbsGU0P2lazzRlmaYPwi4uLkZiYyPa+Xb58GcnJyWAYBseOHUOvXr3YsuHh4exjGxsbdOnSBeHh4QgPD0fXrl3h5+dHvW6EELMyRw8xTZ3SPNCnV0NZu9Y0NzVG45rt2lRWr16Nixcv4vLly7h9+zbUBqZJuXz5slaoCwkJwYYNGxAeHo6AgADweLpTjBBCCCGNEYW6emLue8VZYg9aXTEMg5s3byI3Nxcvv6w9b93KlStx+/ZtvfsJhUKEhISwvW8VWVlZYdy4cfXWZmJ6dFUqIYRUD4W6etIce89MISUlBQkJCUhISMDRo0eRk5ODjh076gS48PBw3L59G3w+H507d2ZPoYaHh6NTp040XYgZmGvQNV2VSggh1UOhrh41xd4zU8vJycHRo0fZIJeSkqJTJikpCYWFhbC3f94j8+GHH2Lq1KkICQmBVTXGHRLTMvega7oqlRBCqsatuoh5rVq1Cm3atIGVlRW6du2KU6dOGSx7/PhxcDgcnUXfdBTE8pw4cQJubm4YNWoU1q9frxPo7O3t8eqrr2LJkiVgKk2QHBwcjG7dulGgayDVGXRNCLEc5cMY9KFhDE2HRfXUbdu2DTNnzsSqVasQFRWFNWvWIDY2Frdv34a3t7fB/ZKTk7VufO7qSt/oLYVMJsP58+eRkJCA7t27Iy4ujt3WtWtX8Pl8KJ/NxycUCtGjRw/ExMQgJiYGERERdCUqIYSYAA1jaB4s6hPzm2++wTvvvIN3330XALBs2TIcPHgQq1evxqJFiwzu17JlSzg6OpqplcQYtVqNxMREJCQk4MiRIzh16hRKS0sBAKNHj9YKdXZ2dnj77bfh5OSEmJgYREVFwcbGpqGaTgghTZq5hjHQJMcNx2JCnVwux5UrV3Qmhu3Xrx/Onj1rdN+wsDCUlZUhMDAQH3/8MXr37m2wrEwmY2+8Dmhu0E7q5p9//sGff/6JhIQEHDt2DPn5+XrLHT16VGdC3zVr1pirmYQQQuoZTXLcsCwm1OXl5UGlUsHNzU1rvZubG7Kzs/Xu4+HhgbVr16Jr166QyWTYtGkTYmJicPz4cfTs2VPvPosWLcKCBQtM3v7mbNu2bfjXv/6ld1urVq3Y06kvv/yyVqAjhBDStNAkxw3LYkJducof+pV7dioqvwdnucjISKSnp2PJkiUGQ118fDxmz57NPpdKpfDy8jJBy5u2goIC7N27Fzt27MDSpUvRoUMHdltMTAz72NHREb1792aDnL+/PwW5JojmjiOEEMtjMaHOxcUFPB5Pp1cuJydHp/fOmO7du2Pz5s0Gt4tEIohEolq3szl58uQJ9uzZg+3bt+PIkSPsBQ2RkZH46KOP2HKhoaH4+uuv0atXL4SFhdFdGpoBGnRNCCGWx2JCnVAoRNeuXXH48GEMHTqUXX/48GEMHjy42se5du0aPDw86qOJzUJubi52796N7du34+jRo1CpVDplKk8zw+PxMGfOHHM1kVgImjuOEEIsi8WEOgCYPXs2xo4di/DwcERGRmLt2rVIS0vDlClTAGhOnWZkZGDjxo0ANFfH+vr6IigoCHK5HJs3b8bOnTuxc+fOhnwZjdann36KL7/8Um+Q8/b2xrBhwzBs2DB069atAVpHqouuPCOEkObJokLdiBEj8OTJEyxcuBBZWVno1KkT/vzzT/j4+AAAsrKykJaWxpaXy+WYM2cOMjIyYG1tjaCgIOzfv19r2gyiX2ZmJlxcXLRup9WhQwetQOfr64s33ngDw4YNQ0REBI2NawToyjNCSEOi8bYNi8NUnqq/mZFKpRCLxZBIJFoTGDdF6enp+P3337F9+3acPXsW+/btwyuvvMJul0gk6NGjBwYNGoRhw4ahS5cuFOQamXtXH2JquP4rkQFg1eX/on2XtmZsESGkuaGzBXVTl1xiUT11xPT++ecf7Ny5E9u3b8f58+e1tu3YsUMr1InFYty6dcvcTSSEENKE0HjbhkOhrglKTU3Fb7/9hh07duDSpUt6ywQGBiIoKMjMLSOEEEJIfaFQ1wStWrUKX3/9tc76zp07sxc7BAYGNkDLCCGEENOg07y6KNQ1YgzD4MKFC+jYsSPEYjG7ftiwYWyoCwsLw7Bhw/D6669rTdRMCCGENFZ0UZh+3IZuAKk5iUSCVatWITQ0FJGRkdi0aZPW9oiICCxfvhz37t3D1atX8dFHH1GgaybKrzzTh648I4Q0FdW5HVlzRD11jQTDMLh06RLWrl2LrVu3oqSkhN32ww8/YNq0aeyVqhwOB9OnT2+oppIGRHd6IISQ5otCnYUrLCzEli1bsGbNGiQmJups7969OyZPngy1Wk235yIA6MozQghprijUWbBDhw7htddeQ3FxsdZ6BwcHvPnmm5g8eTKCg4MbqHWEEEIIsSQU6ixYly5doFA8HzMQERGByZMnY+TIkbC1tW3AlhFCCCHE0lCoswDXrl3DmjVr4OLigi+++IJd7+LiggkTJgAAJk+ejLCwsIZqIjEBuvyeEEJMg25Hph/dJqyBbhNWXFyMX3/9FWvWrGEnCHZ0dERmZiasra3N1g5iHnT5PSGEmFZT/aJMtwlrRG7cuIE1a9Zg8+bNkEqlWtsUCgWuXbuGHj16NFDrSH2pzuX3jflNiBBCzI0uCtNFoc4MFAoFewVr5fuvAkBoaCgmT56M0aNHm7W3kBBCCCFNB4U6M+BwOPj444+RkZHBrrOxscGoUaMwefJkhIeHs3PMEUIIIYTUBoU6M+Dz+XjnnXewcOFCBAcHY/LkyRgzZozWrb0IIYQQYpkay/g9CnVm8t577yE2NhbdunWjXjlCCCGkkWhMF7o1+1CnUqkAAI8ePar38WytW7fWOgVLmo9iRSGUQgWUMt03Bb5IgGJFIR49etQALSOEEGJMyq1/UFgq0butrBRIvnUXcq7MZPWVX0RZnk9qotlPaXLp0iW88MILDd0MQgghhBDWxYsXERERUaN9mn2oKygogLOzM9LT0+nKU0IIIYQ0KKlUCi8vL+Tn58PJyalG+zb70688Hg+A5n6qFOoIIYQQYgnK80lNNPtQRwghpPFoLFchEtIQKNQRQghpFBrTVYiENARuQzeAEEIIqY7q3G6PkOaMQh0hhBBCSBNAoY4QQgghpAmgUEcIIYQQ0gRQqCOEEEIIaQIo1BFCCGkUxC72EFgJ9G4TWAkgdrE3c4sIsSw0pQkhhJBGoaW3KzYkLad56ggxgEIdIYQ88zhHAomkVGe9WGwNt5biBmgRqayltyuFN0IMoFBHCCHQBLpx49dCrlDpbBMKePjfhkkU7AghFo3G1BFCCACJpFRvoAMAuUKltwePEEIsCYU6QgghhJAmgEIdIYQQQkgTQKGOEEIIIaQJoFBHCLF4ycnJcHd3R2Fh9W/YvnLlSgwaNKgeW2XZPvvsM4SGhlrMcQgh9Y9CHSGkQaSnp+Odd96Bp6cnhEIhfHx8MGPGDDx58kSn7Pz58zFt2jTY2+tOLnv//n3Y29vD0dFRa/3EiRNx6dIlnD592mg7xo8fDw6HA/8OHjh2JJ5drl/7iS0jFPAgFlsbPEavXr0wc+ZM4y+4EeBwONi9e7fWujlz5iAhIaFhGkQIqRGLCnUnT57EwIED4enpqffNpbLjx4+Dw+HoLElJSeZpMCGkVh4+fIjw8HDcvXsXW7duxf379/HDDz8gISEBkZGRyM/PZ8s+evQIe/fuxYQJE3SOo1AoMGrUKLz00ks620QiEUaPHo0VK1ZU2Z4BAwYgKysLN/5Oxukz13H6zHUcOLAXP6wajx9WjTfJdCYMw0CpVNbpGA3Bzs4OLVq0aOhmEEKqwaJCXXFxMUJCQrBy5coa7ZecnIysrCx2ad++fT21kBBiCtOmTYNQKMShQ4cQHR0Nb29vxMbG4siRI8jIyMD8+fPZsr/99htCQkLQunVrneN8/PHHCAgIwPDhw/XWM2jQIOzevRulpcanIxGJRHB3d0fnTh0Q1SMYUT2CEREegA7t3ZGZkQSv1q44deoUW37p0qVwcXFBVlYWxo8fjxMnTmD58uXsF8vU1FT2S+fBgwcRHh4OkUiEU6dO4cGDBxg8eDDc3NxgZ2eHiIgIHDlyRKs9vr6++PzzzzF69GjY2dnB09NTJ5ympaVh8ODBsLOzg4ODA4YPH47Hjx8bfI2XLl1C37594eLiArFYjOjoaFy9elWrTgAYOnQoOBwO+7zy6Ve1Wo2FCxeidevWEIlECA0NxYEDB9jtqamp4HA4+P3339G7d2/Y2NggJCQE586dM/pvQAipO4sKdbGxsfjiiy/w2muv1Wi/li1bwt3dnV14PJ7BsjKZDFKpVGshhJhPfn4+Dh48iKlTp8LaWvuUpru7O8aMGYNt27aBYRgAmh788PBwneMcPXoU27dvx/fff2+wrvDwcCgUCly8eLFabctJy8W9qw+1llYO3pj09iSMHTsWEokE169fx/z587Fu3Tp4eHhg+fLliIyMxMSJE9kvll5eXuwx582bh0WLFuHOnTsIDg5GUVER4uLicOTIEVy7dg39+/fHwIEDkZaWptWWr7/+GsHBwbh69Sri4+Mxa9YsHD58GICm12/IkCHIz8/HiRMncPjwYTx48AAjRoww+NoKCwsxbtw4nDp1CufPn0f79u0RFxfHjlO8dOkSAODnn39GVlYW+7yy5cuXY+nSpViyZAlu3LiB/v37Y9CgQbh3755Wufnz52POnDlITExEhw4dMGrUqEbZU0lIo8JYKADMrl27jJY5duwYA4Dx9fVl3N3dmZdffpk5evSo0X0+/fRTBoDOIpFITNh6Qogh58+fN/r/+5tvvmEAMI8fP2YYhmFCQkKYhQsXapXJy8tjvLy8mBMnTjAMwzA///wzIxaL9R7PycmJ2bBhg8H2jBs3juHxeIyNjQ3DA5/hgcfwwGPaIpDpwxnG9OEMY/pbDWc6B3Vmhg8fzgQFBTHvvvuu1jGio6OZGTNmaK0rf3/avXu3kd+GRmBgILNixQr2uY+PDzNgwACtMiNGjGBiY2MZhmGYQ4cOMTwej0lLS2O337p1iwHAXLx4kWEYzXtdSEiIwTqVSiVjb2/P7Nu3j12n79+l8nE8PT2ZL7/8UqtMREQEM3XqVIZhGCYlJYUBwKxfv16nbXfu3DHyWyCEMAzDSCSSWucSi+qpqykPDw+sXbsWO3fuxO+//w5/f3/ExMTg5MmTBveJj4+HRCJhl/T0dDO2mBBSFeZZDx2HwwEAlJaWwsrKSqvMxIkTMXr0aPTs2bPK41lbW6OkpMRomd69e2PPL/vQDX3QDX3RDX3hBT92u0qmxn8/+xo7d+5EaWkpli1bVu3XU7mXsbi4GPPmzUNgYCAcHR1hZ2eHpKQknZ66yMhIned37twBANy5cwdeXl5aPYLlxysvU1lOTg6mTJmCDh06QCwWQywWo6ioSKdeY6RSKTIzMxEVFaW1PioqSqfe4OBg9rGHhwfbBkJI/WnU93719/eHv78/+zwyMhLp6elYsmSJwTd7kUgEkUhkriYSQirx8/MDh8PB7du3MWTIEJ3tSUlJcHJygouLCwDAxcUFBQUFWmWOHj2KvXv3YsmSJQA0QVCtVoPP52Pt2rV4++232bL5+flwdTV+A3hbW1v4ePnChmNnsMzV61fZ4+Xn58PW1rZar7dyublz5+LgwYNYsmQJ/Pz8YG1tjWHDhkEul1d5rPKgyzAM+7giQ+sBzVW+ubm5WLZsGXx8fCASiRAZGVmteg21w1i9AoFAp7xara5xXYSQ6mvUPXX6dO/eXWdsByHEcrRo0QJ9+/bFqlWrdC5gyM7OxpYtWzBixAg2CISFheH27dta5c6dO4fExER2WbhwIezt7ZGYmIihQ4ey5R48eICysjKEhYXVqc0lTBEWffMl1q1bh+7du+Ott97SCihCoRAqlf77xlZ26tQpjB8/HkOHDkXnzp3h7u6O1NRUnXLnz5/XeR4QEABA0yuXlpamdabh9u3bkEgk6Nixo8F6p0+fjri4OAQFBUEkEiEvL0+rjEAgMPo6HBwc4OnpqTNNzNmzZw3WSwgxn0bdU6fPtWvX2K5+QohlWrlyJXr06IH+/fvjiy++QJs2bXDr1i3MnTsXrVq1wpdffsmW7d+/P959912oVCr2IqjKAeLy5cvgcrno1KmT1vpTp06hbdu2aNeundH2yGQy5OblQsaUses44EDIEYFhGNzCRUR1fxETJkxAbGwsOnfujKVLl2Lu3LkANFeOXrhwAampqbCzs4Ozs7PBuvz8/PD7779j4MCB4HA4+OSTT/T2YJ05cwaLFy/GkCFDcPjwYWzfvh379+8HAPTp0wfBwcEYM2YMli1bBqVSialTpyI6OlrvRSXl9W7atAnh4eGQSqWYO3euzoUqvr6+SEhIQFRUFEQiEZycnHSOM3fuXHz66ado164dQkND8fPPPyMxMRFbtmwx+js2lcc5Ekgkulczi8XWdZ52hpDGzqJCXVFREe7fv88+T0lJQWJiIpydneHt7Y34+HhkZGRg48aNAIBly5bB19cXQUFBkMvl2Lx5M3bu3ImdO3c21EsghFRD+/btcfnyZXz22WcYMWIEnjx5And3dwwZMgSffvqpViiKi4uDQCDAkSNH0L9//xrVs3XrVkycOLHKcgcOHNCalgMAbGCPHuiPFNxBKUrw+XxN0HR3d8f69esxfPhw9O3bF6GhoZgzZw7GjRuHwMBAlJaWIiUlxWBd3377Ld5++2306NEDLi4u+Ne//qX3KvwPP/wQV65cwYIFC2Bvb4+lS5eyr798Hs8PPvgAPXv2BJfLxYABA4zOyffTTz9h0qRJCAsLg7e3N/7zn/9gzpw5WmWWLl2K2bNnY926dWjVqpXeHsTp06dDKpXiww8/RE5ODgIDA7F3716zTCX1OEeCcePXQq7Q7U0UCngmmU+QkMaMw5SPSrYAx48fR+/evXXWjxs3Dhs2bMD48ePZ+Z8AYPHixVi7di0yMjJgbW2NoKAgxMfHIy4urtp1SqVSiMViSCQSODg4mOqlEEJMaNWqVdizZw8OHjxY7X1u3ryJmJgY3L17F2Jx1R/0OWm5GB8wA4oyhc42gZUAG5KWo6W38bF5puLr64uZM2c2ibtUmNLde9mYMnWDwe0/rBqPDu3dzdcgQupBXXKJRfXU9erVC8Yy5oYNG7Sez5s3D/PmzavnVhFCGtqkSZNQUFCAwsJCvbcK0yczMxMbN26sVqADgJbertiQtBySPN37y4pd7M0W6AghpLYsKtQRQog+fD5f6y4T1dGvX78a19PS25XCGyGk0aJQRwghFkbfWDZCCKlKk5vShBBCCCGkOaJQRwghpFEQi60hFOi/t7dQwINYbK13GyHNhUlOvxYUFMDBwYGdQ4oQQggxNbeWYvxvwySap44QA2od6u7fv489e/Zgz549OHfuHOzt7REXF4fBgwcjNjYWdnaGb7dDCCGE1IZbSzGFN0IMqNHp1/T0dHz00UcICgpCWFgYTp8+jbfffhvZ2dk4duwY/P398dVXX8HV1RWxsbFYvXp1fbWbEEIIIYRUUKPJh8t75oYMGYJ+/frByspKb7lHjx5h165d2LdvHw4dOmSyxtYHmnyYEEIIIZaiLrmk1neUqMkkoJaMQh0hhBBCLEVdckmtr3596aWXkJ2dXdvdCSGEEEKICdU61IWHh6Nbt25ISkrSWn/t2rUa3XuVEEJI45eTlot7Vx/qLDlpuQ3dNEKajVpf/bp+/XosWLAAL774Inbv3o2WLVvi448/xs6dOzFo0CBTtpEQQpqUxzmSJjUtR05aLsYFz4FCz2geAYeD/91YQrdfI8QM6jRP3aeffgqhUIi+fftCpVKhf//+uHTpErp06VKr4508eRJff/01rly5gqysLOzatQtDhgwxus+JEycwe/Zs3Lp1C56enpg3bx6mTJlSq/oJIaS+Pc6RYNz4tZArVDrbhAIe/rdhUqMLdin3slHSPQBqke7JH4VMjZR72SYLdZlPpXharBuIHW2t4elI46JJ81brUJeVlYVFixZh/fr1CAwMRFJSEkaOHFnrQAcAxcXFCAkJwYQJE/D6669XWT4lJQVxcXGYOHEiNm/ejDNnzmDq1KlwdXWt1v6EEGJuEkmp3kAHAHKFChJJqclC3a8nj+D+k0c66/1atMbInn1MUgcApBVIUdRTDaGtXGebvJiPtAIpupmgnsynUryydAPkSj2BmM/D/g/HmyzYUXgkjVGtQ13btm0REBCA7du345VXXsHBgwcxfPhwPHr0CP/6179qdczY2FjExsZWu/wPP/wAb29vLFu2DADQsWNHXL58GUuWLKFQRwixXD5l4DrrBhN1vunuyvPrySNIdFgNnqPuKdFENQc4CZMFu3RkoW/cVfC4unWp1Byk3/c3ST1Pi0vBExZDbK/U2SaX8fG0uNQkgcuc4ZFUD8MwUKvVWs8N4XK54HKf9xqr1Wooldp/M4b2F4lEWs8VCgVUqud/B1wuF0KhsEZtN6dah7qff/4ZI0eOZJ/3798fx44dw6uvvop//vkHq1atMkkDjTl37hz69eunta5///748ccfoVAoIBAIdPaRyWSQyWTsc6lUWu/tJISQcokFd/HC+38bDECJBXfRAe51ruf+k0d6Ax0A8LgMrqclI+JBGyiVSqhUKnZp0aIFvLy8tMonJCSwH27lC7ufQoHbuffh1sFwXYl3zuPrjxOhUiigVqmgVqs1PyssC957DzwOB1CrAZUK+06cwOELFzRl1WqolEoUu9qgzwdK8Hh6fncqDorPuAECJ8z5+Wc8fPyYDQIVfzIMA/WznxNefBEjIyIAhgEYBgWFhRi8ahUU4haQ9x2h9/XIlSq8Et0bW2P7IqBlS81KhsG+W7fw0V9/sXUwgPbPZ4+dra1xYcoUTZ3P9p3555/YdecOyl9VednyxwDAAHitY0esfPVVrfZ0/O475JWUsGUq71v++IdBgzAyOJjd73p2Nl5cu5bdj62Lw9GqFwBSFyxASwcHgMsFOBwsPnQIn+zbp9M+VDwOgFAfH1z+z38AHo9dXv7kExz7+2+9v9uK/jV+PL6aPZvdT65UwqpTpyr3A4CjR4+id+/e7PM//vgDgwcPrnI/gUAAuVy7t3nmzJlaeWbAgAH466+/qtWOhlDrUFcx0JXr0qULzp49a7arX7Ozs+Hm5qa1zs3NDUqlEnl5efDw8NDZZ9GiRViwYIFZ2kcIaVyS81KQW5yns97V1gX+Lm201uXn56O0tJT9oqhvKSsrQ2hwMPy8vACZDCgrQ/LdG+AFGw5Am7asx28z/gWFQoE148bB3dYWUCgAhQLbLl7EfxMSoFCpoFAqNT/LF7UaCrUaSrUavra2GDVpFNDO8Gvdv/9PfDX2fZ3108RirHR2BlQqQKkEVCrEPn4MhZHfW+eRcRjYR2Rwe+LtO9j0659GjgB8sm4dKvZTngWwolIZ76j2eJOn/4Odx2Mg++4r4Mw9HAVwzWhtGtE3bwI//MA+VwI4BcDKrRRtjeyX9PgxSv77X611TwHcrEadLgDw8cda654ASKvGvk8vXgQuXtRal/dsqYrit9+A335jn6sBFFVjPwDARx9pPVUB0D3RrkuVkgKMGqW1rtoT427YoFlq4+WXAaEQEAg0S4UePqOUSmDwYM0+5ftfuqRdJjkZWLsWmDSpdm2rZ3W6UEIfX19fnDlzxtSHNYjz7FtFufJvCJXXl4uPj8fs2bPZ51KpVOdbKSGkcVGr1SgpKUFJSQmKi4t1HhcXFUFeUoIxgwcDxcVASQlQUoJfdu1CwoULmnJiLgI/d9HfC5TPwYQvnOCfW8qGs/CjR5FSqjvmqrKVAPwqPC8eGQdesOEA9E/OE/x9/jwAYMmVK1p9dk9QvbBSIpEA6WkADNdjiFIiASQSrXU8wGioU0P/GMHqbgcAdcuWz3t0uFxw9bSjSn5+AOMKzvXrmn/nKjDt2wPt22t6p7hccORy4ODBKvfjcTjAgAGAs7NmXw4HgtRU2Jw7Bw6HAw6g+/PZY0eRCBg+nO0RA4cDxzNn4J6SwpZDhf3KHwNAC19fICrqWeM1f6e+u3bBoayMLVOxrorPHV54AfDxYfe1KihA0NGj7GviVPoJhmEf86KjNQHnWY+me0oKwlJS2Dbo3ZfDQQdrayAwUPMF4dnS8cEDlJaUsMdiy1d4DIaBt7U1YGXF7sdRKtGj0n6VP+XLn4sBQC7XLABaAOiJqvEZBti7V2tdewC9KjwPTkkBNm1qGqFu//79eO+991BUVITQ0FDMmTMHcXFx+PTTT3H+/Hm89NJLmDhxok7vWX1xd3fXmQA5JycHfD4fLVq00LuPSCTSOWdOCGlYKpUK586dg0QigeTJE0geP4YkJweSvDxICwogefoUEokE0qIifBMbi65OTmww+yMpCQMrfDgZwgcwptK6cwB+evbYO6o9OvP0v2/weAxyL5+E/5l77LrqvovIKj0XQQ3dEWHPMarnfSCKHj0AR0e2x0GUlgbB5csQcLnPFx6PXfhcLgR8Plo7OgIhoQDuGKzHy6UFOkdHg8/ng8fjgcfng8/n46XQUCA29nnA4vPx0bp1UDEMeHy+ZhEIwBcIwBMIwOPzcUdcCmP9VLGv9MXCNyaBy+WCx+OxY54qLsJevTT1PTMlPR1DsrO1ytwueoib2GT4l7dgAeATgb9ycqBUKsHlcjXhiMNhH1f8KRKJgAqfBy0YBjKFAsnZeRi5epvBai5evozAVtqfcyOfLbWxArq9ktV1ac2aWu3XEdXrWdRnwrOlNmo7MEsIQKe7iGHY0/VQKtke7cpLlEKBE3K5/u1VrJ/5bNFa36aNnhZahhqFujlz5mDYsGGIi4vDX3/9hddeew2vvvoqDh48iLFjx+Kvv/7C6tWrcezYMXTo0KG+2syKjIzEvmfn9csdOnQI4eHhesfTEUJMh3n2bZylVuPQnj24e/MmpHl5mnCWn68JalIpJIWFkJSUQFJairfbtcOCtm01PTESCRiJBC/dv1+teh+v0P74s65me5XQnDISCgSAjQ1gawvbkhLg6dPqHWDkSOB1J00IEInQZ9s2dJRIYGVlpfmyaGUFkbU1RFZWsLKxgejZ0jM6GujWTbOfQADHQ1uRhx0Gqxk1ajje+uZnCAQCzZdT/vO36XeeLdXxxa4NMBbqIl9+GR+vGF+tY32ywnjk2PPwOC4WGY4IL/boisFte1WrrnJeXl46Z1GepJXi5lPD+8g5mh7BluVj3WqIw+FAKBSCx+PB2qYMQpH+CzKIheBwnn/5sOCLF8ypRn+daWlpmD59Onx9fdGnTx8EBARgypQpWLZsGT744AMAmkGF8+fPx/bt22vcmKKiItyv8MaekpKCxMREODs7w9vbG/Hx8cjIyMDGjRsBAFOmTMHKlSsxe/ZsTJw4EefOncOPP/6IrVu31rhuQohGdnY2bl6/jtyUFOSmpSE3IwO52dnIzc1FbkEBcqVS5BYXw1MoxI2AACAvD3jyBCgsxAqGwR/VqeP6deD6dfY5H4AtgKpPmAHFvXppTunY2gI2NnB5+hTd9+yBjbU1bG1sYGNjAxtbW9ja2cHGzg62YjFsHBxgIxYDU6Zo9ntmemYmxhYUwMbGBrcLH+Isvjdc8cCBgE8E+3TFO9WNV9r4fONXuNo7OKB169a1OnZFfi1aI1HNMXhBhl+LutdRLsDZG+ekhusKcPY2ST18of5hNdXdXl08UQn6vHoFPJ7uWCyVigueqLZ9coTUrxqFOl9fX1y4cAG+vr4AgDfffBOTJ09GZGQkW2bq1KmIjo6uVWMuX76sdcVK+di3cePGYcOGDcjKykJa2vPhpG3atMGff/6JWbNm4fvvv4enpye+++47ms6EEGiu9H78+DFyHz9GbmqqZklPR25mJnJzcpD75Alynz7FLz16wFcm04SzvDzsycjAlGdX0xnDlcl0BhFXNbuaDZ8PsZUVbENDgTFjALEYcHAAxGLM/eUXcKysIG7ZUrM4OcHBwQFisRg8ezUEdlzY2dlBIBAg49nxbPkOCBG64tyzaY1qytPTE56engCAjH/ygBoO4aoNd4EzMoyELXeBs0nqGdmzD3ASZpmnzt+xLSb4fo3HZbq3BHOzcoW/o7HLDqpPxLWDysjvTsS1M0k9djZqvYEOAHg8NexsqjnwnhAzq1GomzdvHt59910kJSUhLi4OwcHBOHv2LDp27MiWKR+cXBu9evUyOvfMBj1XwkRHR+Pq1au1qo+QxkqpVCIrMxPpf/+N9Js34W9jg1CRCMjIADIykH7/PrxPnKjWsTJ37IBvhedVzftvKxDA1c4O7s7OYL75BhxXV6BFC0AsxsQbN9AvMxNiR0eIxWKtxcHBweiwiE9ffFHv+qfyXHyb/D6UTxSaKwUq4HMEmOW/Eo7CxnMLqmjfUKwf5QC1u26vEjebweytoSary5TBrSr+jm3hb/Sa0bqz5Ttj/4NOEPF1T4vKlHy86W2aQGxuNNExMZUahbpx48bBwcEB33zzDT7//HNwuVwEBAQgLCwMYWFhCAgIwBdffKHVc0cIqYWyMiAzE0f37sXf164h/Z9/kJ6VhfQnT5BeWIhMuRwV+wriAYRWeO4OzZVg1Zk+IH/sWODFFwEXF8DFBZ2KivDvAwfg6uUFV3d3uLq6ai3W1oZHsUX37VubV2tUsVIKJaP/2kslo0CxUmqSUOdq6wJVPsfgHGiuti51rgMAWnq7YtPWxZDkFepsE7vY0z1Sq1CiFKJE2XTGT9FEx8SUajzic+jQoRg6dCiKiopw/fp1JCYmIjExEVu2bMGtW7dQVlYGT09PvP766wgODkZwcDCGDh1aH20npFF6mpqKfy5dQvqtW0i/fx/p//yDtOxsCMvK8FOLFpretjzNzFNLAFRnmst0JyegRw+gVSugVSsIWrXCq+vWQWBnB1cPD7i2bq35WSmgubi46MyO3gHAIjPNNWlJ7K2dcSg1GFyu7gxcarUQ0zuYrheopbcrhTcCQHOXDH2BDtBMdGyqu2SQ5qHWl/HY2dkhKioKUeVz5kAzLUFSUhIb9E6fPo1Vq1ZRqCPNC8MAjx8DDx4A9+/j4okTWHfyJJIeP0ZycTFyDQwxcATw06MK459EIngLhUChdo+Oq4MDvNzd4eXjAy8/P3j5+iI8PFwz4WYFe2s5kL+5ypeVQCLnw9DbYr6sBJ62prknK6kdZ5ENRFweZGrdECTi8uAssmmAVhFiOWoU6nJycmBtbQ17e3u923k8HoKCghAUFIQxY8bg7t27ZpnahBCzU6uBR4+Qf+0aks6dQ/KNG0h+8ADJ2dn4UqFAYIVJaTMBrK/GIZ8CKNq5E3Z+fpoeN2dnjD17Fj0ePGCnd2jdujWsrKzq6UURYtk8bcU4Evce8mW6F/I4i2xMFrpt+Q7gcwR6T/vzOQLY8qnnjFimGoW6EydOYOzYsejVqxcGDx6MwYMHs1eOlbtw4QL27NmDXbt2IT09HUVF1b4RCSGWRaEA/vkHuH8fzP37+OPwYSTdu4fkzEwkSaVIZhi9t+cZASCQywW8vQE/P/g7O7O35/FwdYV/QADatm/PBrWKi52d9tV7lXvDm6O8IuMXXuUVFaMVddA0G5624nrvMXUUumKW/0oUK3XvDW7Ld2hUF+aQ5qVGoe6NN95A9+7dsWfPHuzcuRMzZ85ESEgIBg4ciIyMDOzduxcqlQqvvPIKvvrqK/Tr16++2k2I6Tx+DNy4gScXLiDp6lUk370LcW4uXn/yRDNTOTQXHUwGkFWNwyV/8AGwZAk7GaafQoGLc+bA398fDg70Db+m5GUCqFSGL2CQl9FE48T0HIWuFN5Io1PjMXVeXl54//338f7770MikWDfvn3Yv38/PD09sW3bNkRFRYHL5dZHWwmpG5kMuHMHuH4dTy5cwJXz53H57l1cLi7GZQDpFYq+BOB1ALC21txPsl07+F+5gqz056U8PT3h7++PgIAArZ/e3t5Ahf8DAoEAERERILVjy3PGkT/CDc7u/9bbjXMaC0IAzbQlQj7P4NWvjrbVvWcKIXW4UAIAxGIx3nzzTbz55pumag8hdccwQGam5o4FN248X5KSsFmlwv8BSKniEMlOTsDNm4CHB3vj7dn79uFdqRT+/v7o0KED9bqZUWmJFUqrng+5TmgQPmkIno4O2P/heJqnjpgE3cSONG4lJcCtW2xwK7x6FdeuX8flwkJcBrAIgE+F4kJbW6TomRzbzs4OXbp0QVBQENvrxnh4aN3bdODAgfX+ckjDMdcgfEIq83R0oPBGTIJCHWkcGEZz0UKFnreSxEQk3r+PywyDywAuA0iC9oS7Q6Oi4DNwIBAcDAQHI7ysDDbBwQgLC0N4eDi7dOjQgYYNELMMwieksqfyXLoog5gEhTpimTIzgTNngLNnNfcX/ftvQPr8Te9FAOcAVHUHxisvvog3/vUv9nkbhoFEIgGfT3/6jQWNOSJNGXsbPAPTpzS22+CRhkWfbKThqVSa8WtnzrBBriA1FccBJADgAvgO0FxNGhgIBAeDf/o01A8fah1GIBAgJCREqwcuMDBQqwyHw6FA18jQmCPSlJnrNnikeaBPN2J+hYXA+fOaXrgzZ4Dz51FSWIgz0IS4BABX8bwXzsXeHsvPnAEnIAB4dkP47v/+NyQHD2oFuE6dOkEkEjXMa2qGzHnKiMYcEUJI1Swu1K1atQpff/01srKyEBQUhGXLluGll17SW/b48ePo3bu3zvo7d+4gICCgvptKqoNhgLQ0rV443LgBqNV4BGADNCHuLADdO25q5BcXI9PZGa0Ez+cjW7RoEb766qt6bz7Rj04ZEUKI5bGoULdt2zbMnDkTq1atQlRUFNasWYPY2Fjcvn1bM/eXAcnJyVrTS7i60odJg1EogMTE571wZ84AmZlgAMgAsDe48vXFk8BAfPLnn3oP07lzZ8TExCAmJgY9e/bUmT6k4lWpxPzolBEhhFgeiwp133zzDd555x28++67AIBly5bh4MGDWL16NRYtWmRwv5YtW8LR0dFMrSRaCgqAc+ee98JdvKiZZgRAKp6dTuVwcJTHw7SICHwyezbQowfg6YnOajVc3NyQl5eHNm3asCHu5ZdfRsuWLRvyVRFCCCGNjsWEOrlcjitXruDf//631vp+/frh7NmzRvcNCwtDWVkZAgMD8fHHH+s9JVtOJpNBJpOxz6VS3TFBxIiiIuDoUeDAAeDECeD2bXZTLoBjABKEQiTw+XjwLNyBYQClEglCIT4ZNowtz+Vy8euvv6Jt27Zo06aNeV8HIYRYAFu+A/gcgcGhDLZ8GktKqs9iQl1eXh5UKhXc3Ny01ru5uSE7O1vvPh4eHli7di26du0KmUyGTZs2ISYmBsePH0fPnj317rNo0SIsWLDA5O1vshgGSE4G/voL+PNP4ORJQK49+u03d3f8RybD9YICzQq5XKeMnZ0dWrRoAYZhtE6dxsTE1PtLIIQQS+UodMUs/5U0Tx0xCYsJdeUqj5WqHAIq8vf3h7+/P/s8MjIS6enpWLJkicFQFx8fj9mzZ7PPpVIpvLy8TNDyJqSkBDh2TBPi/voLSHl+U61HAFr6+kIYFwf07Qv06AHF4cO4XulWcQKBAJGRkejTpw9iYmIQEREBgYBuvE4IIZU5Cl0pvBGTsJhQ5+LiAh6Pp9Mrl5OTo9N7Z0z37t2xefNmg9tFIhFNe6HPvXvPQ9zx40CFU9T/CATY2aYNdsjlOJeair9WrcKA2Fh2+8svvwwul4vQ0FB2XNyLL74IW1vbBnghxBzolBEhhFgeiwl1QqEQXbt2xeHDhzF06FB2/eHDhzF48OBqH+fatWvw8PCojyY2LaWlmvBWflr1wQOtzQ89PbHT2xvbnzzBpXv3gLt32W07du7UCnUeHh548uQJXazSjNApI0IIsTwWE+oAYPbs2Rg7dizCw8MRGRmJtWvXIi0tDVOmTAGgOXWakZGBjRs3AtBcHevr64ugoCDI5XJs3rwZO3fuxM6dOxvyZViuBw+eh7hjx4CysufbBAI86NoVvzk5Ycc//+Dq7duaW3VV0rlzZ3Tq1ElnPQW65odOGRFCiGWxqFA3YsQIPHnyBAsXLkRWVhY6deqEP//8Ez4+PgCArKwspKWlseXlcjnmzJmDjIwMWFtbIygoCPv370dcXFxDvQTLUlamubCh/LRqhd42AEDr1kBcHBAbC8TEYNVnn+Gbb77ROUxoaCiGDRuGYcOGaY1hJIQQ0nhkPpXS7faaOA7DMExDN6IhSaVSiMViSCQSnQluG6XU1Och7uhRds44AACfD7z4Im6FhWGHTIZ34+PRqnVrdvO5c+fQo0cPAEDXrl3ZIOfn52fmF0EIIQ0rs1iCfFmJznpnkQ08bcUN0KK6yXwqxStLN0CuVOlsE/J52P/heAp2FqIuucSieupILT1+DGzbBmzZopn8tyJPTzADBuDvoCDsePwYO/btw51vvwUAOAcE4IMPPmCLduvWDcuXL8fAgQNp3jhCSLOVWSxBnz9XQ6bWDUAiLg9H4t5rdMHuaXGp3kAHAHKlCk+LSynUNQEU6hqrwkJg1y5NkDtyBFCrNeu5XCAqCkxsLBLbtsWO69exY+dO3P3pJ51D7NixQyvUcblcTJ8+3VyvgBBCLFK+rERvoAMAmVqFfFlJowt1pHmgUNeYyOXAwYOaILd3r+YK1nIvvACMGQOMGIGf9u/Hf/7zHzyodEUroJkH8MUXX8SwYcPw2muvmbHxpKmh8TmEEGJZKNRZOrVac1/VLVuA7duB/Pzn2zp00AS50aOBCuPe5HK5VqDjcrno2bMnhg0bhqFDh8LT09Ocr4A0QTQ+hxBCLA+FOkv199+aILd1K1Dhil+4uwOjRoEZNQoXlEqsWbsWs0pKEFxh16FDh2L69Ono2bMn3njjDQwZMqRGEzgTUhUan0MIIZaHQp0l+ecfTYj75RdNqCvn4AC8/jowejQkXbpgy6+/Ys277+LGjRsAAGtra6xatYot7ubmhry8vKZxNS8hhJA6c7S1hpDPM9i77mhr3QCtIqZGoa6hPXmiOa36yy/AqVPP1wuFmjnkxowBExeHSzdvYu3atdg6eDBKSrQvs9+/fz9UKhV4PB67jgIdIYTUjrPIBiIuz+DVr84imwZoVd14Ojpg/4fjaRxsE0ehriGUlAD79mlOrx44ACie3T+TwwGiozXj5F5/HYV8PrZs2YI1UVFITEzUOUz37t0xefJkDB8+XCvQEUIIqT1PWzGOxL3XpOapAzTBjsJb00ahzlyUSiAhQRPkdu0CioqebwsN1QS5kSM1d3l45tSff+K9997TOoyDgwPefPNNTJ48GcHBwSCEEGJ6nrbiRhveSPNFoc4c4uOBn34CcnKer/P1fX7lamAgCgsLkZeXh4pT/vbv3x/e3t5IS0tDREQEJk+ejJEjR8LW1tbcr4AQLTQ+hxBCLA+FOnNIS9MEOhcXYPhwTZiLjAQ4HFy7dg1rpkzBli1bEBkZiUOHDrG78Xg8rF69Gh4eHggLC2vAF0CINhqfQwghlofu/WqOe79evQpkZwN9+wICAYqLi/Hrr79izZo1uHTpklbRe/fu0b1WCSGEkGaK7v1q6bp0AQBcv34da9aswebNm1FYWKhVxNbWFqNHjwafT/8khBBCCKk5ShBmUFpaipdffhnnz5/X2RYaGorJkydj9OjRNA0JIQZkFkua3JWIhBBiahTqzMDa2hrW1s8HjtvY2GDUqFGYPHkywsPDweFwGrB1hFi2zGIJ+vy52uCcYUfi3qNgR0g1PJXnolgp1Vlvy3eAo9C1AVpETI1CnZlMmjQJT548weTJkzFmzBiIxfQhREh15MtK9AY6AJCpVciXlVCoI6QKT+W5+Db5fSgZhc42PkeAWf4rKdg1Ac0+1KlUmg+LR48e1evpzx49euCPP/4Ah8NBYWGhzpg6Qoh+jwtyoHwiMbw9IwuOJfpDHyFEI7s0FfnZhv8f3bdKhru1zIwtIoZIpZre1PJ8UhPN/urXS5cu4YUXXmjoZhBCCCGEsC5evIiIiIga7dPsQ11BQQGcnZ2Rnp5OFyoQQgghpEFJpVJ4eXkhPz8fTk5ONdq32Z9+Lb9nqoODA4U6QgghhFiE2tzTnVsP7SCEEEIIIWZGoY4QQgghpAmgUEcIIYQQ0gRQqCOEEEIIaQIo1BFCCCGENAEU6gghhBBCmgAKdYQQQgghTQCFOkIIIYSQJoBCHSGEEEJIE0ChjhBCCCGkCaBQRwghhBDSBFCoI4QQQghpAijUEUIIIYQ0ARTqCCGEEEKaAAp1hBBCCCFNAIU6QojFS05Ohru7OwoLC6u9z8qVKzFo0KB6bJVl++yzzxAaGmoxxyGE1D8KdYSQBpGeno533nkHnp6eEAqF8PHxwYwZM/DkyROdsvPnz8e0adNgb2/PrmMYBkuWLEGHDh0gEong5eWF//znP+z2iRMn4tKlSzh9+rTRdowfPx4cDkdnGTBgQLVfS69evTBz5sxql7dUHA4Hu3fv1lo3Z84cJCQkNEyDCCE1YlGh7uTJkxg4cCA8PT31vrlUdvz4cb1vxklJSeZpMCGkVh4+fIjw8HDcvXsXW7duxf379/HDDz8gISEBkZGRyM/PZ8s+evQIe/fuxYQJE7SOMWPGDKxfvx5LlixBUlIS9u3bhxdeeIHdLhKJMHr0aKxYsaLK9gwYMABZWVlay9atW033gqEJoUql0qTHNAc7Ozu0aNGioZtBCKkGiwp1xcXFCAkJwcqVK2u0X3Jystabcfv27euphYQQU5g2bRqEQiEOHTqE6OhoeHt7IzY2FkeOHEFGRgbmz5/Plv3tt98QEhKC1q1bs+vu3LmD1atXY8+ePRg0aBDatGmD0NBQ9OnTR6ueQYMGYffu3SgtLTXaHpFIBHd3d63FyckJgObLo1AoxKlTp9jyS5cuhYuLC7KysjB+/HicOHECy5cvZ79Ypqamsl86Dx48iPDwcIhEIpw6dQoPHjzA4MGD4ebmBjs7O0RERODIkSNa7fH19cXnn3+O0aNHw87ODp6enjrhNC0tDYMHD4adnR0cHBwwfPhwPH782OBrvHTpEvr27QsXFxeIxWJER0fj6tWrWnUCwNChQ8HhcNjnlU+/qtVqLFy4EK1bt4ZIJEJoaCgOHDjAbk9NTQWHw8Hvv/+O3r17w8bGBiEhITh37pzRfwNCSN1ZVKiLjY3FF198gddee61G+7Vs2VLrzZjH49VTCwkhdZWfn4+DBw9i6tSpsLa21trm7u6OMWPGYNu2bWAYBoCmBz88PFyr3L59+9C2bVv88ccfaNOmDXx9ffHuu+9q9fABQHh4OBQKBS5evFjr9pafWh07diwkEgmuX7+O+fPnY926dfDw8MDy5csRGRmJiRMnsl8svby82P3nzZuHRYsW4c6dOwgODkZRURHi4uJw5MgRXLt2Df3798fAgQORlpamVe/XX3+N4OBgXL16FfHx8Zg1axYOHz4MQNPrN2TIEOTn5+PEiRM4fPgwHjx4gBEjRhh8HYWFhRg3bhxOnTqF8+fPo3379oiLi2PHKV66dAkA8PPPPyMrK4t9Xtny5cuxdOlSLFmyBDdu3ED//v0xaNAg3Lt3T6vc/PnzMWfOHCQmJqJDhw4YNWpUo+ypJKRRYSwUAGbXrl1Gyxw7dowBwPj6+jLu7u7Myy+/zBw9etToPmVlZYxEImGX9PR0BgAjkUhM2HpCiCHnz583+v/7m2++YQAwjx8/ZhiGYUJCQpiFCxdqlZk8eTIjEomYbt26MSdPnmSOHTvGhIaGMr1799Y5npOTE7NhwwaD7Rk3bhzD4/EYW1tbraVinTKZjAkLC2OGDx/OBAUFMe+++67WMaKjo5kZM2ZorSt/f9q9e7exXwfDMAwTGBjIrFixgn3u4+PDDBgwQKvMiBEjmNjYWIZhGObQoUMMj8dj0tLS2O23bt1iADAXL15kGIZhPv30UyYkJMRgnUqlkrG3t2f27dvHrtP371L5OJ6ensyXX36pVSYiIoKZOnUqwzAMk5KSwgBg1q9fr9O2O3fuGPktEEIYhmEkEkmtc4lF9dTVlIeHB9auXYudO3fi999/h7+/P2JiYnDy5EmD+yxatAhisZhdKn6jJoQ0POZZDx2HwwEAlJaWwsrKSquMWq2GTCbDxo0b8dJLL6FXr1748ccfcezYMSQnJ2uVtba2RklJidE6e/fujcTERK1l2rRp7HahUIjNmzdj586dKC0txbJly6r9eir3MhYXF2PevHkIDAyEo6Mj7OzskJSUpNNTFxkZqfP8zp07ADSnn728vLTev8qPV16mspycHEyZMgUdOnRg3/+Kiop06jVGKpUiMzMTUVFRWuujoqJ06g0ODmYfe3h4sG0ghNQffkM3oC78/f3h7+/PPo+MjER6ejqWLFmCnj176t0nPj4es2fPZp9LpVIKdoSYkZ+fHzgcDm7fvo0hQ4bobE9KSoKTkxNcXFwAAC4uLigoKNAq4+HhAT6fjw4dOrDrOnbsCEAz1qzi+0J+fj5cXV2NtsnW1hZ+fn5Gy5w9e5Y9Xn5+PmxtbY2Wr3jsiubOnYuDBw9iyZIl8PPzg7W1NYYNGwa5XF7lscqDLsMw7OOKDK0HNFf55ubmYtmyZfDx8YFIJEJkZGS16jXUDmP1CgQCnfJqtbrGdTWUzGIJ8mW6XwacRTbwtBU3QIsIqVqjDnX6dO/eHZs3bza4XSQSQSQSmbFFhJCKWrRogb59+2LVqlWYNWuW1ri67OxsbNmyBW+99RYbBMLCwnD79m2tY0RFRUGpVOLBgwdo164dAODu3bsAAB8fH7bcgwcPUFZWhrCwsDq1+cGDB5g1axbWrVuH3377DW+99RYSEhLA5WpOdgiFQqhUqmod69SpUxg/fjyGDh0KACgqKkJqaqpOufPnz+s8DwgIAKDplUtLS0N6ejr7pfT27duQSCRsuNVX76pVqxAXFwdAM6VMXl6eVhmBQGD0dTg4OMDT0xOnT5/W+uJ89uxZrSuPG7vMYgn6/LkaMrXu70LE5eFI3HsU7IhFatSnX/W5du0a29VPCLFMK1euhEwmQ//+/XHy5Emkp6fjwIED6Nu3L1q1aoUvv/ySLdu/f3+cO3dOK2z06dMHXbp0wdtvv41r167hypUrmDx5Mvr27avVe3fq1Cm0bduWDX6GyGQyZGdnay3lgUelUmHs2LHo168fJkyYgJ9//hk3b97E0qVL2f19fX1x4cIFpKamIi8vz2iPlJ+fH37//XckJibi+vXrGD16tN7yZ86cweLFi3H37l18//332L59O2bMmMG+/uDgYIwZMwZXr17FxYsX8dZbbyE6OlrndG/Fejdt2oQ7d+7gwoULGDNmjM6FKr6+vkhISEB2drZO72i5uXPn4r///S+2bduG5ORk/Pvf/0ZiYiLbtqYgX1aiN9ABgEyt0tuDR4glsKhQV1RUxI5nAYCUlBQkJiayYz7i4+Px1ltvseWXLVuG3bt34969e7h16xbi4+Oxc+dOvP/++w3RfEJINbVv3x6XL19Gu3btMGLECLRr1w6TJk1C7969ce7cOTg7O7Nl4+LiIBAItKb94HK52LdvH1xcXNCzZ0+88sor6NixI3799VeterZu3YqJEydW2Z4DBw7Aw8NDa3nxxRcBAF9++SVSU1Oxdu1aAJordNevX4+PP/6Yfa+aM2cOeDweAgMD4erqanSc2rfffgsnJyf06NEDAwcORP/+/dGlSxedch9++CGuXLmCsLAwfP7551i6dCn69+8P4PkkwU5OTujZsyf69OmDtm3bYtu2bQbr/emnn1BQUICwsDCMHTsW06dPR8uWLbXKLF26FIcPH4aXl5fB3s3p06fjww8/xIcffojOnTvjwIED2Lt3L00lRYgF4DDlo5ItwPHjx9G7d2+d9ePGjcOGDRswfvx4dv4nAFi8eDHWrl2LjIwMWFtbIygoCPHx8ezpheqQSqUQi8WQSCRwcHAw1UshhJjQqlWrsGfPHhw8eLDa+9y8eRMxMTG4e/cuxOLGdarM19cXM2fObBJ3qWiMbuZnYfDhnwxu39P3bXRypjNCpH7UJZdY1Ji6Xr16wVjG3LBhg9bzefPmYd7/t3fvcVGW+f/HX8MMDDBySEEEAQ+tgsdSsVLzlGaptdp2sjK1o3Y2d9tq6/vtsJnbybUyTXf7bQctK8391uZaVqaVpkJipqJ5BE+BB0AGmWGY+/fHwMgIKBoxA7yfj8c8mPu672vuzzjivL3u+77uP//5N65KRPztzjvv5OjRoxw7dsznVmGnsn//ft5+++0GF+hERM5WQIU6EZHqWCwWn7tM1MawYcN+o2pE6la+Mw+7q7BKu80SSXTIqa/cFqlMoU5EJMBUdzWs1J/m1nCsQeYar35tbg2vs33lO/P4+9Z7cRmlVdZZTME8mDJTwU5qTaFORESkkgRbFF+MuKte5qmzuwqrDXQALqMUu6tQoU5qTaFORKTc/vxC8u3Hq7RH28JIiNaFVE1Jgi1Kc9FJg6NQJyKCJ9CNfOlNnK6qh9xCLGY+/eMEBTsRCWgBNU+diIi/5NuPVxvoAJyusmpH8EREAolCnYiIiEgjoFAnIiLiJzZLJBZTcLXrLKZgbBYd8pfa0zl1IiIifhIdEsuDKTM1T53UCYU6ERERP4oOiVV4kzqhw68iInimLQmxmKtdF2IxE20Lq+eKRETOjEbqRESAhOhIPv3jBM1TF+D22wvqZVJgkYZIoU5EpFxCdKTC21nKzc6j4NCxKu1RMRG0TK6bQ4v77QUMXTK7xtt3fTHirjoLdrofqzREARXqVq5cyQsvvEBGRgYHDhxg8eLFjB49+pR9VqxYwZQpU9i0aRMJCQn8+c9/ZtKkSfVTsIiIkJudx019/4QjxKiyzuo0MX/Vi3US7I44iqsNdAAOdxlHHMV1Eup0P1ZpqAIq1Nntds477zxuueUWrr766tNuv2vXLkaMGMEdd9zBvHnz+O6777j77ruJjY2tVX8REX/YemgXefZDVdpjbTGkxLSrs/0sWPkF2w/vrdL+uxaJjBkwtM728/OBgxx8KRlrqLvKuqMlQfx84GCdjdaFW5xYLa4q7Q5X3X2d6X6s0lAFVKgbPnw4w4cPr/X2r7/+OsnJycyYMQOATp06kZ6ezosvvqhQJyIBaeuhXfwr50+YzVVHtcqOmLiFF+sk2C1Y+QWZkbMxR1fdT6bbBCups2CXeWgHI1M2Yw6q5j25TWTu20E/uv2qfRiGQVHpYUae+1ON+7G7jgDx5OfnU1paimEYuN1un5+Vn0dHRxMdHe19DZfLxbZt2zjs3nfKWrZs2cI5neIJDw/3tuXl5fHzzz9791HTIzg4mP79+/u83o8//sj+/ftxu93e92oYhvd5xc/ExER69uzp03fJkiWUlJT49Kvcp+L5hRdeSNu2bcEwwDA4cuQIS5cu9ayv2K7S9gBGeT3XXn01YaGhFR8EP27cyA/r1/vWeHLNhkGL5s256oorvP0A/rN0KfsPHPB5b5z0XgHO79KFi9LSvOvKysqY89ZbJ2o7qW+lFYwaNozE+Hjv8s49e/j0yy9P2syo0jfIZOLeceNOtBsGy7//nh+zsrzLbVq3ZvSVV0JSEoEooELdmVq9ejXDhg3zabvssst44403KC0tJTi46oSODocDh8PhXS4srHrOhIjIbyU9eytmS9VQAmA2G6z6eSMxpkhKS0uJiYnBYjnxz3Rubi579uyhtLS0+ofTicvpJDwkhO1BRdUGOgBzkMFHSxfzwxvv4CotpaysjDKXizKXiwHdunHDwIFQVgYuF5SVMeaFF3A4nZ7tyspwlf+seG4Z1JO+99e8r/kfvMX0sbfgNgzchkGZ2+19XvHYP348ISYTuN1QVsbj69bxwk8/edYDbsMg8ZJzGff37jXux/G3/4VNRQxevZrMY1XP7zvZs4mJPNqqledL3O3miNNJl02biEuN4rb3L6mx39ixY1lyrBk9Q0K8IWlpURHjDh8+7T5jTCbyWrTw9sMweL6oiPml1Y8MVnaT2cy8iu+18r7jnU6qjvlW9TbQttLyHuCmWvQDuPz226l87fenwF9q0e984KqT2l4Cvq5F34eBiyotu4B7atEPoPPjj5NYafkn4P5a9AsG7n3qKZ+2hcCsSsuXA6MXLICVK2tZTf1q0KHu4MGDxMXF+bTFxcXhcrk4dOgQ8RVJvZJp06bx1Ekfmog0XS6XC7vdTnFxMblF+zh6/BAlx0twOByUlpbidDqJi4pnSO+h4HB4HiUlvPbPf/LLL7/gOH78xMPhwFFS4nk4nZQ4HNw3YAC/T0nx9tttPgq/r7me5155hVsXjAJge8+enGuxQGkplJbyfm4u9+fmnvY9nQtMGDMCOllr3CZr124+XLCkSrvp44+5YepUn7bFgPMU++vWOwWoeV9FTic78/NPWbN7zhzf5dPss1obNsBXOzDVcnNj717Ye+Lw9JnM8WXs8x3Nq/U+DQMO+caw2valrMzzkFMLCvI8TOV/sobh+Q9Kbdhsnn4VfUtKPL9/FSwWaNasbuutQw061AGYTL6/DhVDsie3V3j00UeZMmWKd7mwsJCkAB1GFRHP73RxcTEFBQUUFhaSmJhIs0r/qO7ZsYP/fvwx9vx8igsKsBcWUnzsGMV2O/aiIoqLi7EXF1NWWsry228Hux2Ki8Fu5+6vvmL2tm0ARLYKY9LHl2KxmiEMKg9NbHSU0Sv5aqIPnpjuZAawvRb1X7l5s8+yacwQ+H3trrB1/fCDz3JIrXrB6cd8alZ2zjmeQ0tms+cLzGzGvHatZwTtLIUEBRETGkqQyYQ5KIggk+nEo3zZuOsusFq9X8gJa9bQIz3duz4oKIiIiOhT7idq0j0wIYYB8+fTOjeXoKAgTCYTpvL+J//sdPHF0LevNwCEOhxMmDWLkKRTfzVed8klxE4dC/Hx3gDQKSuLBz/7zLO/Svv1eW4yYQsPh4kTTwQHk4mrli3j3G3bvNsA3n4VzwG6duwIl112InCYTPzvm29yvKSkSr+Tl9MGD4aUFG+/xMOHeWXxYu97MpXXUvm7s2K/zW66CSoONZtMXL5hA80r/d2sqW/zc86B0aN96v3z0qWMq+Y/JhX7qtCtWzdIS/P2s5SV8fa77/r2qeF7vtPQodCqlXe55969zK9mZO3k/kFBQXD99T5tt6Sn03/7id/0Vq1awaBB1e43EJgM4+QD0oHBZDKd9urXAQMG0KNHD15++WVv2+LFi7nuuusoLi6u9vDryQoLC4mKiqKgoIDISE1lIFKX3G43RUVFFBQU+DySWremW7t2UFAABQWUHj7MrU895Vl/7BgFRUUUFBdTUFJCodNJWaV/ppa1bs1Qw/AGsyWlpYysZT1l+I7GTAH+Xv78dIfc7r7+K1pnFXgWgoLoahhsqsU/n8937cpD3bp5AovVyvOtQygYVfXihQpr51gI+W4bwcHBvHL77SQnJEBwMAQHsyori/eWLyc4JIRgq9X3Z2golvLl6ObN2X2OmeMd/1Pjfo6t7sM1nS7EYrFgNpsxm81YLBZiYmJITEz02Xbfvn0EBQV5t6nYvmJ5+rIPKExcVOO+ovddzUPDa3uwr2b7incwa/tDNa6/+3cv0Dr83F+9H139Kv70a3JJgx6p69OnD5988olP2+eff05aWlqtAp2I1I7b7ebIkSPk5eV5Hrm5mI4f5w99+kBeHhw+DAUFPPjGG3yxeTOFx497A1l1secBPCNdFSzAe3hC1+kUn3TYy3YG76Nk7FjCIyM9h1jCw/ndhg30+/FHwsPCaJF6mqkwli+HiI6ecGax8NqKFZSUlBAaGorVavV5VG4LCwvzjHiVC/t2MQW8U+Nufn/DDdz32slnInn0HTKEvvfU7syiZxa/ecr1LVvGMmDAgFq9VuvWrU+5vkVUFKc6O/mcyIY1KbDuxyoNVUCFuqKiIrZXGubctWsXmZmZNG/enOTkZB599FH27dvH22+/DcCkSZOYOXMmU6ZM4Y477mD16tW88cYbvPfee/56CyINQsV5p3l5eXTo0IFQqxWOHoW8PL787395/YMPyDt8mLz8fPKOHeNwSQnuk0al2gN/OOl1s/GclHw6BRVPgoMhKgpTVBSRu3dztPx8oZCgIKKsViJDQ4kKDyeqWTOiIiOJiooiYcwY6NnTG8xS7Xb+tXIl4ZGR2Jo1Izw8nPDwcGw2m8/PsLAwn4sOAO4uf8DpR4GIjITwExFy4MCBtXinVQXZzZQ1M9V4BWeQvfpblZ2p37VIJNNd835+1yKxml5np0tsG7bbTdVf0VtmoktsmzrZj80SicUUXOMIms1Sd0dbdD9WaYgCKtSlp6czePBg73LFuW/jx4/nzTff5MCBA2RnZ3vXt2vXjiVLlvDggw/y2muvkZCQwCuvvKLpTKTJKywo4O05cziwYwd5+/eT98sv5B054glpRUUcrXQF+I8tWtCtoMB7IvE+PFd8nU4eeIJVbCy0aAHR0URt307o3r1EWa1EhYV5AllEhCeQRUcTdc45RMXE0DMtDa69FkJDvefbpO/cic1mIyoqitCK6RNqIQ6YUHGuUAPQ1hTH69u6YrVWM9eaw8LfouOq6XXmxgwYCiupl3nq+nTsBtueZO/RX6qsSzwnzrO+DmgETeTUAirUDRo0iFOd4vfmm29WaRs4cCA/nHQysUhjYxgG+fn55OTkkJOTw96dO8nZvJmcHTvIyclhbEoKt8bGQk4OZGfjzM7mPru9Vq+dV3kahogIYm02OHgQgHCLhdjwcGIjI4lt3pzY2FhiW7UiNjGR2DZtMCZN8jnZ+J9uN/8v6EyuITyhffv2Z9WvoekYE0uLW3ZT1rzqZQ/Njjjp+EndBZO6DG6n4wludRPeTkUjaCI1C6hQJ9JU2e12jhw54rkS2+WCAwcgO5sHp05laUYGOUePYj/FPFZdKybHLNcCCAVKTtouIjiYWJuN2Kgob0hrceut0K+fZ8TNamWA3c6uvDxiY2Ox2c7kjLXyq8caqPo8tBd8oJjgA1VvSi8i8mso1InUA6fTyfbt29m6dSvbMjPZs3kzObt3k3PgANmHD3O0pIR2Vis7W7aE/fu9c1HtA7JO/dIA5LVvD+PGQXIyJCVhSk5mXno6ETExntG18ofVWvNcYhVsNtsZh7nGQIf2RKShU6gTqSOGYZCbm+s5Jyw4GLKzYcsW5s+fz/j33vOZlqM6ex0O3Dk5nik3LBZITCSptJSwgwdJiooiKTaWpMREks49l6ROnUhKSSEpOZmkpKRqL3u/umPH3+aN+kG+M69ewpYO7YlIQ6ZQJ3KGHA6Hd9Qta+NGtqanszUri6ycHAocDr5q357B+/d7ZiIHWlPzVB3BJhOtIyJIiokhOSGB43/9K7aOHSEuDsxmpjmdvBgcXOMkm01BY5szLComguDQYEpLqr6f4NBgomIi/FCViDQGCnUip2MYGHl5/OGaa9iYlcWuQ4eqTO9R2dadOxkMEBICHTuS2rYtPdLTSWnXjtTu3Um56CLad+5MUlIScXFxpzwPLSSktvcQaLzsrsJqAx2AyyjF7ipsUKGuZXIsb2a9TMGhqvcmjYqJoGVyw3kvIhJYFOqkycvNzSUjI4P1GRlkZWSwdfNmukVH889u3SArC7ZswXTkCFnAjlO8TptmzUhJSKDl1VfDbbdB27ZgNtMK0PXZUlnL5FiFNxGpcwp10qTk5+ezbt060tetI33lStLT08muPKVHueMAa9eeaDCZSA0LI8fhIKVFC1LbtiWla1dS+/QhJS2NDh07El5xb0QRERE/UKiTRqugoACbzea5i8ChQ7BuHf+YPp0/f/HFKfuZAFdkJMbkyZg6dYJOnaBDB95xu7HZbE36/DYREQlcCnXSKBQVFbF+/XrS09NJX7OG9NWr2ZadzZphw7hgxw7Y4TlwmnZSv2ZBQfSKjyftvPPoNXQoXQYNokNqqud+nSdvWw/vQ0RE5Gwp1EmDtHr1ak+AW7eO9FWr2LJzZ7V3I0n//HMuqFhISaFnjx7cf/AgaYMH03v0aDp27dqgJ8xtCupzUmARkYZMoU4CmtPp5ODBgyQnJ3saDh6EtWu5cdw4dhcU1NjPCpwfHU3kkCEwcSL07u25Nynwcr1ULnWlPicF3m8v4Iij6p0emlvDSbBF1dl+5Oztzy8k3368Snu0LYyEaAV8adoU6iSguN1uMjMz+fLLL/nyiy/4ZuVKOsXEkN6nD6xZ45nQF89h1N3lfYKB7s2akda+PWkXXUTaqFF0GTqUYE0H0mjUx6TA++0FDF0yG4e76qyC1iAzX4y4S8HOz/bnFzLypTdxuqp+RiEWM5/+cYKCnTRpCnXiV4Zh8PPPP3tC3Jdfsvyrrzhy9KjPNj/u3Yvjww+xAphM0LkztyYkcElMDL1HjqTb6NFYm+BtraRuHXEUVxvoABzuMo44ihXq/CzffrzaQAfgdJWRbz+uUCdNWsCFulmzZvHCCy9w4MABunTpwowZM+jfv3+123799dcMHjy4SvuWLVtITU39rUuVX2nNmjVcc8017N27t8ZtEkwmBrduTeGECcRecgmkpUFEBMPrsU4REZGGIKBC3fvvv8/kyZOZNWsW/fr1Y86cOQwfPpzNmzefOKeqGlu3bvW592VsrCb1DCT5+fl8/fXXJCcn07NnT8jNhc8+o93ixVUCXTQwOCqKIRdcwJCxY0kZMwaTDqOKiIicVkCFuunTp3Pbbbdx++23AzBjxgw+++wzZs+ezbRp02rs17JlS6Kjo+upSjmd48eP891333kPqWZkZOB2u5nYsyevm0yQkQFAS+AiIMJiYUhqKkOuuooed9yBOSnJr/WLiIg0RAET6pxOJxkZGTzyyCM+7cOGDWPVqlWn7NujRw9KSkro3Lkzjz/+eLWHZCs4HA4cDod3ubCw6hV1cmZcLhcZGRneEPfdd9/5/BlX+PKHSjfL6tEDhg9n1eWXY+rTBywB81dRRESkQQqYb9JDhw5RVlZGXFycT3tcXBwHDx6stk98fDxz586lV69eOBwO3nnnHYYMGcLXX3/NgAEDqu0zbdo0nnrqqTqvvymbOnUqTz75ZI3ruwJDrVaG9O6NceutmC6/HOLjAc/dG0QCQXNrONYgc41Xvza36jZw/hZtCyPEYq7x6tdoW9VJw0WakoAJdRVOvgWTYRg13pYpJSWFlJQU73KfPn3IycnhxRdfrDHUPfroo0yZMsW7XFhYSJIO952WYRj89NNPLFy4kLvuuotWrVrB/v2wdCkDv/rKZ9s2wFBgSLt2XDJqFHHXXgsXXKDROAloCbYovhhxl+apC2AJ0ZF8+scJmqdOpAYB8y0bExOD2WyuMiqXm5tbZfTuVC666CLmzZtX43qr1YrVaj3rOpsSwzDIzMxk4cKFLFy4kG3btgEQt3kzd+/a5T03rg9wAzAoPJwhl1xC+2uv9YzGtWzpv+JFzkKCLUrhLcAlREcqvInUIGBCXUhICL169WLZsmVcddVV3vZly5YxatSoWr/O+vXriS8/tCdnzjAMMjIyvEFuR/k9UytbuHAhd4NnzrjevbEOH867l1/uuWuD2VzvNYuIiEgAhTqAKVOmcPPNN5OWlkafPn2YO3cu2dnZTJo0CfAcOt23bx9vv/024Lk6tm3btnTp0gWn08m8efNYtGgRixYt8ufbaLDmzp3LtGnT2L17d5V1QcAA4Brgqm7d4M474brrNBonIiISIAIq1F1//fUcPnyYp59+mgMHDtC1a1eWLFlCmzZtADhw4ADZ5beJAs8Vs3/605/Yt28fYWFhdOnShU8//ZQRI0b46y00GG63G5PJ5HO+YmlpqU+gCwIGUx7kYmKIGz8exo+Hbt3qu1wRERE5DZNhGIa/i/CnwsJCoqKiKCgo8JnAuDEqKytj1apVfPjhhyxatIgPPviAfv36QWkpLFnC/tmzaffZZwzCE+RGBwcTO2qUJ8hddhkEB/v5HYiIiDRuvyaXBNRIndS9srIyVq5cycKFC/noo498LkRZOHs2/RYuhPnzIS+PBCAPiExLgwkTYMwYaNHCX6WLiIjIGVCoa4RcLhdff/01CxcuZPHixeTm5lbZJsRkonj+/BMNcXFw881Ejh8PXbvWY7UiIiJSFxTqGqEpU6bw6quvVmm3BgUx3DC4xjC40jCIDAmB3//eMyp32WWaR07OSL4zD7ur6h1ZbJZIokN0/2URkfqmb/EG7ujRo9hsNkIq3fT+iiuu8Ia6MIuFEUFBXON0MtLtJgI8U49UHF5t3twvdUvDlu/M4+9b78VllFZZZzEF82DKTAU7EZF6FuTvAuTMGYbB6tWrmTBhAgkJCXz00UcnVubmMnjjRsafcw4fAHkuFwudTsa0akXEQw/BTz/B2rVw990KdHLW7K7CagMdgMsorXYET0REflsaqWtACgoKmDdvHnPmzGHjxo3e9tdff50xvXvD3/4Gb71FcGkpbwKEhMCoUZ5RuWHDdHhVRESkEdO3fIAzDIN169YxZ84cFixYQHGx730poyIiOO+XXyjr2BGz2+1pTEuDW27R4VUREZEmRKEugK1Zs4ZJkyaRmZlZZd1F3bszMTSU69auJTwry9N4+eXwP/8DffvWb6EiIiLidwp1ASwmJsYn0EVGRjL2ssuYeOgQ3ZcvP7HhqFHw2GOeCyBERESkSVKoCwDHjh3jvffeIyQkhAkTJnjbzz33XC699FLy8/OZOHQoYzIysH34oWelyQTXXOMJc+ed55/CpcmyWSKxmIJrvPrVZmncd2cREQlEuk2YH28Ttn79eubMmcP8+fMpKioiKSmJXbt2YTabvdvYP/sM24svwhdfeBqCguCGG+Avf4HOneu1XpHKNE+diEjd023CGhC73c6CBQuYM2cO69at81mXk5PDN998w6CBAz0h7q9/xfbNN56VFguMGwePPAIdOvihchFf0SGxCm8iIgFEoa6ebNiwgTlz5jBv3jyOHTvms85ms3HjjTcy8c476fXLL9CnD6xZ41kZEgK33goPPwxt29Z/4SIiItIgKNTVg+PHjzNgwAAKC30PVfXo0YOJEydyw/XXE7l8OUycCD/84FkZGgp33gkPPQSJiX6oWkRERBoS3VGiHoSFhTF27FgAwsPDue2221i7di0Za9cyMSqKyP794Q9/8AQ6m80T5HbtgpdfVqATERGRWtFIXT2599576dKlCzfddBNRNhu8+y7cfDNs3erZIDIS7rsPJk+GmBi/1ioiIiINT5O/+rWgoIDo6GhycnJ++6tfnU547z2YPh127/a0RUfDXXfBpEme5yIiItJkFRYWkpSURH5+PlFRUWfUt8mHur1795KUlOTvMkRERES8cnJySDzDU7CafKhzu93s37+fiIgITCbTr369ioRdLyN/clb0GQU+fUaBT59Rw6DPKfCd/BkZhsGxY8dISEggKOjMLn1o8ufUBQUFnXESro3IyEj9AgU4fUaBT59R4NNn1DDocwp8lT+jMz3sWkFXv4qIiIg0Agp1IiIiIo2AQl0ds1qtPPHEE1itVn+XIjXQZxT49BkFPn1GDYM+p8BXl59Rk79QQkRERKQx0EidiIiISCOgUCciIiLSCCjUiYiIiDQCCnUiIiIijYBCXR2bNWsW7dq1IzQ0lF69evHNN9/4uyQpN23aNHr37k1ERAQtW7Zk9OjRbN261d9lySlMmzYNk8nE5MmT/V2KVLJv3z7Gjh1LixYtCA8P5/zzzycjI8PfZUk5l8vF448/Trt27QgLC6N9+/Y8/fTTuN1uf5fWZK1cuZIrr7yShIQETCYT//73v33WG4bBk08+SUJCAmFhYQwaNIhNmzad8X4U6urQ+++/z+TJk3nsscdYv349/fv3Z/jw4WRnZ/u7NAFWrFjBPffcw/fff8+yZctwuVwMGzYMu93u79KkGuvWrWPu3Ll0797d36VIJUePHqVfv34EBwfz3//+l82bN/PSSy8RHR3t79Kk3HPPPcfrr7/OzJkz2bJlC88//zwvvPACr776qr9La7LsdjvnnXceM2fOrHb9888/z/Tp05k5cybr1q2jVatWXHrppRw7duyM9qMpTerQhRdeSM+ePZk9e7a3rVOnTowePZpp06b5sTKpTl5eHi1btmTFihUMGDDA3+VIJUVFRfTs2ZNZs2bxzDPPcP755zNjxgx/lyXAI488wnfffaejEAHsiiuuIC4ujjfeeMPbdvXVVxMeHs4777zjx8oEwGQysXjxYkaPHg14RukSEhKYPHkyDz/8MAAOh4O4uDiee+45Jk6cWOvX1khdHXE6nWRkZDBs2DCf9mHDhrFq1So/VSWnUlBQAEDz5s39XImc7J577mHkyJEMHTrU36XIST7++GPS0tK49tpradmyJT169OAf//iHv8uSSi6++GK+/PJLtm3bBsCGDRv49ttvGTFihJ8rk+rs2rWLgwcP+uQHq9XKwIEDzzg/WOq6uKbq0KFDlJWVERcX59MeFxfHwYMH/VSV1MQwDKZMmcLFF19M165d/V2OVLJgwQJ++OEH1q1b5+9SpBo7d+5k9uzZTJkyhb/85S+sXbuW+++/H6vVyrhx4/xdngAPP/wwBQUFpKamYjabKSsrY+rUqdxwww3+Lk2qUZERqssPe/bsOaPXUqirYyaTyWfZMIwqbeJ/9957Lz/++CPffvutv0uRSnJycnjggQf4/PPPCQ0N9Xc5Ug23201aWhrPPvssAD169GDTpk3Mnj1boS5AvP/++8ybN493332XLl26kJmZyeTJk0lISGD8+PH+Lk9qUBf5QaGujsTExGA2m6uMyuXm5lZJ3+Jf9913Hx9//DErV64kMTHR3+VIJRkZGeTm5tKrVy9vW1lZGStXrmTmzJk4HA7MZrMfK5T4+Hg6d+7s09apUycWLVrkp4rkZA899BCPPPIIY8aMAaBbt27s2bOHadOmKdQFoFatWgGeEbv4+Hhv+9nkB51TV0dCQkLo1asXy5Yt82lftmwZffv29VNVUplhGNx777189NFHfPXVV7Rr187fJclJhgwZwsaNG8nMzPQ+0tLSuOmmm8jMzFSgCwD9+vWrMhXQtm3baNOmjZ8qkpMVFxcTFOT79W42mzWlSYBq164drVq18skPTqeTFStWnHF+0EhdHZoyZQo333wzaWlp9OnTh7lz55Kdnc2kSZP8XZrgOfn+3Xff5f/+7/+IiIjwjqpGRUURFhbm5+oEICIioso5jjabjRYtWujcxwDx4IMP0rdvX5599lmuu+461q5dy9y5c5k7d66/S5NyV155JVOnTiU5OZkuXbqwfv16pk+fzq233urv0pqsoqIitm/f7l3etWsXmZmZNG/enOTkZCZPnsyzzz5Lhw4d6NChA88++yzh4eHceOONZ7YjQ+rUa6+9ZrRp08YICQkxevbsaaxYscLfJUk5oNrHv/71L3+XJqcwcOBA44EHHvB3GVLJJ598YnTt2tWwWq1GamqqMXfuXH+XJJUUFhYaDzzwgJGcnGyEhoYa7du3Nx577DHD4XD4u7Qma/ny5dV+/4wfP94wDMNwu93GE088YbRq1cqwWq3GgAEDjI0bN57xfjRPnYiIiEgjoHPqRERERBoBhToRERGRRkChTkRERKQRUKgTERERaQQU6kREREQaAYU6ERERkUZAoU5ERESkEVCoExEREWkEFOpEREREGgGFOhEREZFGQKFORKQOvPXWW3Tu3Jnw8HBSU1P5z3/+4++SRKSJUagTEfmVFi9ezD333MPjjz/OTz/9xPDhw5k0aZK/yxKRJsZkGIbh7yJERBqyiy++mEsuuYSnn34agGXLlnHttdeSn5/v38JEpEnRSJ2IyK9w7NgxVq9ezciRI71tS5cu5fzzz/dfUSLSJFn8XYCISEO2YcMGTCYT3bt3p7i4mPnz5/Pqq6+yaNEif5cmIk2MQp2IyK+QmZlJamoqmZmZ9O3bF4CrrrrKZ+RORKQ+6PCriMivkJmZSY8ePejatStr1qxhxowZfP755zzxxBP+Lk1EmhiN1ImI/AqZmZnceOONREREcMEFF3DBBReQlZXF999/7+/SRKSJ0UidiMhZcrlcbNq0idTUVJ/2DRs20L9/fz9VJSJNlUbqRETOUlZWFiUlJTzzzDPEx8cTHh7O7Nmz2bVrF3fccYe/yxORJkahTkTkLGVmZhIfH4/NZqN///7YbDYuvvhili9fTnx8vL/LE5EmRqFOROQsZWZmcuGFF7J48WJ/lyIionPqRETOVmZmJt27d/d3GSIigEKdiMhZ27Bhg0KdiAQM3ftVREREpBHQSJ2IiIhII6BQJyIiItIIKNSJiIiINAIKdSIiIiKNgEKdiIiISCOgUCciIiLSCCjUiYiIiDQCCnUiIiIijYBCnYiIiEgjoFAnIiIi0gj8f+AqRNOjbt8gAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(len(orders), sharex=True, sharey=True)\n", "\n", "nsampvals = np.array((10.0 * np.ones(5)) ** np.arange(1, 6), dtype=int)\n", "nsampcolors = plt.cm.viridis(np.arange(0.0, 1.0, float(1.0 / len(nsampvals))))\n", "\n", "# First plot the analytical result\n", "for a in ax:\n", " a.plot(betas, betas * idealgas.x_ave(betas, vol), \"k--\", linewidth=2.0)\n", "\n", "# Next look at extrapolation with an infinite number of samples\n", "# This is possible in the ideal gas model in both temperature and volume\n", "for j, o in enumerate(orders):\n", " trueExtrap, trueDerivs = idealgas.x_beta_extrap_depend(o, beta_ref, betas, vol)\n", " ax[j].plot(betas, trueExtrap, \"r-\", zorder=0)\n", " if j == len(orders) - 1:\n", " print(f\"True extrapolation coefficients: {trueDerivs}\")\n", "\n", "for i, n in enumerate(nsampvals):\n", " thisinds = rng.choice(len(xdata), size=n, replace=False)\n", "\n", " # Get parameters for extrapolation model with this data by training it - the parameters are the derivatives\n", " xem_dep = xtrap.beta.factory_extrapmodel(\n", " beta=beta_ref,\n", " xalpha=True,\n", " data=xtrap.DataCentralMomentsVals.from_vals(\n", " uv=udata[thisinds],\n", " xv=xdataDepend.sel(rec=thisinds),\n", " deriv_dim=\"deriv\",\n", " central=True,\n", " order=order,\n", " ),\n", " )\n", "\n", " out = xem_dep.predict(betas, cumsum=True)\n", " print(f\"\\t With N_configs = {n:6}: {xem_dep.derivs(norm=False).values.flatten()}\")\n", " for j, o in enumerate(orders):\n", " out.sel(order=o).plot(\n", " marker=\"s\",\n", " ms=4,\n", " color=nsampcolors[i],\n", " ls=\"None\",\n", " label=f\"N={n}\",\n", " ax=ax[j],\n", " )\n", "\n", "ax[2].set_ylabel(r\"$\\langle \\beta x \\rangle$\")\n", "ax[-1].set_xlabel(r\"$\\beta$\")\n", "\n", "for j, o in enumerate(orders):\n", " ax[j].annotate(\"O(%i) Extrapolation\" % (o), xy=(0.4, 0.7), xycoords=\"axes fraction\")\n", "\n", "ax[-1].set_ylim((0.0, 2.0))\n", "ax[-1].yaxis.set_major_locator(mpl.ticker.MaxNLocator(nbins=4, prune=\"both\"))\n", "\n", "fig.tight_layout()\n", "fig.subplots_adjust(hspace=0.0)\n", "\n", "for a in ax:\n", " a.set_title(None)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "In each of the figures, the true behavior is shown as a black dashed line, the analytical result (infinite sampling) for each order of extrapolation is shown as a red solid line and the results with 10, 100, 1000, ... 100000 randomly sampled configurations are shown with purple (fewer samples) to green (more samples) circles. For higher order extrapolation, the analytical, infinite sampling result matches very closely with the true temperature dependence of $\\langle x \\rangle$. However, the finite-sampling results are in practice very poor due to difficulties in accurately estimating the higher-order moments of the potential energy distribution. The higher orders are actually quite accurate (if you try zooming in) close to the point we're extrapolating from, but the error grows very quickly as we move further away." ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Tags", "hide_input": false, "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.8" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }