{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Cosmic shear\n", "\n", "This example simulates a galaxy catalogue with shears affected by weak lensing, combining the [galaxies](../1-basic/density.ipynb) and [lensing](../1-basic/lensing.ipynb) examples with models for the intrinsic galaxy ellipticity and the resulting shear." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "The setup of galaxies and weak lensing fields is the same as in the basic examples.\n", "The setup for angular matter power spectra matches the definition from the [Matter shell definition](../1-basic/shells.ipynb) example." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import healpy as hp\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# use the CAMB cosmology that generated the matter power spectra\n", "import camb\n", "from cosmology import Cosmology\n", "\n", "# almost all GLASS functionality is available from the `glass` namespace\n", "import glass\n", "import glass.ext.camb\n", "\n", "# how many arcmin2 over the entire sphere\n", "ARCMIN2_SPHERE = 60**6 // 100 / np.pi\n", "\n", "# creating a numpy random number generator for sampling\n", "rng = np.random.default_rng(seed=42)\n", "\n", "# cosmology for the simulation\n", "h = 0.7\n", "Oc = 0.25\n", "Ob = 0.05\n", "\n", "# basic parameters of the simulation\n", "nside = lmax = 256\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", "pars = camb.set_params(\n", " H0=100 * h,\n", " omch2=Oc * h**2,\n", " ombh2=Ob * h**2,\n", " NonLinear=camb.model.NonLinear_both,\n", ")\n", "\n", "# get the cosmology from CAMB\n", "cosmo = Cosmology.from_camb(pars)\n", "\n", "# shells of 200 Mpc in comoving distance spacing\n", "zb = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n", "\n", "# linear window function for shells\n", "ws = glass.linear_windows(zb)\n", "\n", "# compute the angular matter power spectra of the shells with CAMB\n", "cls = glass.ext.camb.matter_cls(pars, lmax, ws)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Matter" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# apply discretisation to the full set of spectra:\n", "# - HEALPix pixel window function (`nside=nside`)\n", "# - maximum angular mode number (`lmax=lmax`)\n", "# - number of correlated shells (`ncorr=3`)\n", "cls = glass.discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)\n", "\n", "# compute Gaussian spectra for lognormal fields from discretised spectra\n", "gls = glass.lognormal_gls(cls)\n", "\n", "# generator for lognormal matter fields\n", "matter = glass.generate_lognormal(gls, nside, ncorr=3, rng=rng)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lensing" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this will compute the convergence field iteratively\n", "convergence = glass.MultiPlaneConvergence(cosmo)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Galaxies" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# standard deviation in each component of galaxy ellipticity\n", "# this is very small so that the galaxy density can be small, too\n", "sigma_e = 0.01\n", "\n", "# galaxy number density per arcmin2, over all shells\n", "n_arcmin2 = 0.01\n", "\n", "# localised redshift distribution with the given density\n", "z = np.arange(0.0, 2.0, 0.01)\n", "dndz = np.exp(-((z - 0.5) ** 2) / (0.1) ** 2)\n", "dndz *= n_arcmin2 / np.trapezoid(dndz, z)\n", "\n", "# distribute dN/dz over the radial window functions\n", "ngal = glass.partition(z, dndz, ws)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation\n", "\n", "Simulate the galaxies with shears.\n", "In each iteration, get the shears and map them to a HEALPix map for later analysis.\n", "\n", "In addition, generate the galaxy ellipticities, drawn from the intrinsic normal distribution.\n", "The standard deviation is much too small to be realistic, but enables the example to get away with fewer total galaxies.\n", "\n", "Finally, apply the reduced shear from the lensing maps to the galaxy ellipticities, producing the galaxy shears." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# number of HEALPix pixels in the maps\n", "npix = 12 * nside**2\n", "\n", "# map for galaxy numbers\n", "num = np.zeros(npix)\n", "\n", "# map for sum of shears\n", "she = np.zeros(npix, dtype=complex)\n", "\n", "# simulate the matter fields in the main loop\n", "for i, delta_i in enumerate(matter):\n", " # compute the lensing maps for this shell\n", " convergence.add_window(delta_i, ws[i])\n", " kappa_i = convergence.kappa\n", " gamm1_i, gamm2_i = glass.shear_from_convergence(kappa_i)\n", "\n", " # generate galaxy positions uniformly over the sphere\n", " for gal_lon, gal_lat, gal_count in glass.uniform_positions(ngal[i], rng=rng):\n", " # generate galaxy ellipticities from the chosen distribution\n", " gal_eps = glass.ellipticity_intnorm(gal_count, sigma_e, rng=rng)\n", "\n", " # apply the shear fields to the ellipticities\n", " gal_she = glass.galaxy_shear(\n", " gal_lon,\n", " gal_lat,\n", " gal_eps,\n", " kappa_i,\n", " gamm1_i,\n", " gamm2_i,\n", " )\n", "\n", " # map the galaxy shears to a HEALPix map; this is opaque but works\n", " gal_pix = hp.ang2pix(nside, gal_lon, gal_lat, lonlat=True)\n", " s = np.argsort(gal_pix)\n", " pix, start, count = np.unique(gal_pix[s], return_index=True, return_counts=True)\n", " num[pix] += count\n", " she[pix] += list(map(np.sum, np.split(gal_she[s], start[1:])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "Compute the angular power spectrum of the observed galaxy shears.\n", "To compare with the expectation, take into account the expected noise level due to shape noise, and the expected mixing matrix for a uniform distribution of points." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACaGElEQVR4nOzdd3hT5fvH8XfSXboolF32XmXvDYKAgCAOXLhQ/KJs/QEiKCqIIiAKIiggioiI4AIEkYoMQfbeFJBC2d275/dHmtOTNm3TJmnScL+ui8uTk5NzniLQT59xPzpFURSEEEIIIUSxp3d0A4QQQgghhG1IsBNCCCGEcBES7IQQQgghXIQEOyGEEEIIFyHBTgghhBDCRUiwE0IIIYRwERLshBBCCCFchAQ7IYQQQggX4e7oBhR3GRkZREZG4u/vj06nc3RzhBBCCOFiFEUhNjaWChUqoNfn3Scnwc5KkZGRhIaGOroZQgghhHBxly9fplKlSnleI8HOSv7+/oDhNzsgIMDBrRFCCCGEq4mJiSE0NFTNHHmRYGcl4/BrQECABDshhBBC2I0lU75k8YQQQgghhIuQYCeEEEII4SIk2AkhhBBCuAgJdkIIIYQQLkKCnRBCCCGEi5BgJ4QQQgjhIiTYCSGEEEK4CAl2QgghhBAuQoKdEEIIIYSLkGAnhBBCCOEiJNgJIYQQwiaeeeYZHnzwQbs/56233qJJkyZOcx9nIsFOCCGEEDbx8ccfs2zZMkc3wyydTse6detMzo0fP54tW7Y4pkF24u7oBgghhBDCNQQGBjq6CQXi5+eHn5+fo5thU9JjJ2wiKSmJgwcPkp6e7uimCCGEsLMffviBRo0a4ePjQ6lSpejRowfx8fE5hmK7dOnCq6++yujRoylZsiRly5Zl8eLFxMfH8+yzz+Lv70/NmjXZsGGD+plly5YRFBRk8rx169ah0+lybc+///7LfffdR+nSpQkMDKRz587s379ffb9q1aoADBw4EJ1Op77OPhSbkZHBtGnTqFSpEl5eXjRp0oSNGzeq70dERKDT6fjxxx/p2rUrvr6+hIWFsWvXroL/JtqJBDthtZs3bxIWFkbTpk155513HN0cIYQQdnT16lWGDBnCc889x4kTJwgPD2fQoEEoimL2+q+++orSpUuzZ88eXn31VV5++WUefvhh2rVrx/79++nZsydPPfUUCQkJhW5TbGwsQ4cOZfv27fzzzz/UqlWLPn36EBsbCxiCH8DSpUu5evWq+jq7jz/+mI8++ohZs2Zx+PBhevXqRf/+/Tlz5ozJdW+88Qbjx4/n4MGD1K5dmyFDhpCWllbo9tuSDMUKq6SmpvLII49w+vRpAL7//nveeustxzZKCCGKqRYtWnDt2rUif265cuXYu3evRddevXqVtLQ0Bg0aRJUqVQBo1KhRrteHhYUxefJkACZOnMj7779P6dKlGTZsGABTpkzhs88+4/Dhw7Rp06ZQ7e/WrZvJ60WLFhEUFMRff/3FAw88QEhICABBQUGUK1cu1/vMmjWL//u//+Oxxx4DYObMmWzdupW5c+cyf/589brx48fTt29fAN5++20aNGjA2bNnqVu3bqHab0sS7IRVxo8fz9atW9XXp0+fJjk5GS8vLwe2KneXLl3C09Mzz7/YtpSSkoKnp2eRPEsIUfxdu3aNK1euOLoZeQoLC6N79+40atSIXr160bNnTwYPHkzJkiXNXt+4cWP12M3NjVKlSpkEwbJlywJw/fr1QrcpKiqKyZMnEx4ezvXr10lPTychIYFLly5ZfI+YmBgiIyNp3769yfn27dtz6NAhk3Par6l8+fJq+yXYiWJt6dKlzJs3z+Rceno6J0+eJCwszEGtyt327dvp0qULiqLw8MMP83//9380bdrU5s9RFIXw8HBmzZrFhg0beOyxx/j2229t/hwhhOspqh86rXmum5sbmzdvZufOnWzatIlPPvmEN954g927d5u93sPDw+S1TqczOWecO5eRkQGAXq/PMaybmpqaZ5uGDh3KrVu3+Pjjj6lSpQpeXl60bduWlJQUi7+ugsir/Y4mwU4Uyu7duxk+fLj6unnz5uzbtw+AI0eOOGWw+/XXX9XFHatWrWLVqlX07NmTCRMm0KVLlzwn5loiLS2NH374gVmzZqm/FwArV65k0qRJNGzY0Kr7CyFcn6XDoY6m0+lo37497du3Z8qUKVSpUoW1a9fa5N4hISHExsYSHx9PiRIlADh48GCen9mxYwcLFiygT58+AFy+fJmbN2+aXOPh4ZHnAr+AgAAqVKjAjh076Ny5s8m9W7VqVcivpujJ4glRYFevXmXQoEHqT0Ivv/wy06ZNU98/evSoo5qWp7Nnz+Y4t2nTJrp160br1q358ccfC/UTV1xcHB9//DE1a9ZkyJAhJqHOaOXKlYVqsxBCOJvdu3czffp09u7dy6VLl/jxxx+5ceMG9erVs8n9W7duja+vL5MmTeLcuXN8++23+dbGq1WrFl9//TUnTpxg9+7dPPHEE/j4+JhcU7VqVbZs2cK1a9e4c+eO2fu89tprzJw5k1WrVnHq1CkmTJjAwYMHGTVqlE2+tqIgwU4USHJyMoMGDSIyMhKAjh07MnfuXJP5EkeOHHFU8/JkDHbu7u58+umnVKtWTX3v33//5aGHHqJ+/fp8+eWXJCcn53u/q1evMmnSJEJDQxk9ejQXL15U32vatCkLFy5Erzf8FVu5cmWuK8aEEKI4CQgIYNu2bfTp04fatWszefJkPvroI3r37m2T+wcHB/PNN9+wfv16GjVqxMqVK/NdlPfll19y584dmjVrxlNPPcXIkSMpU6aMyTUfffQRmzdvJjQ0NNdpOCNHjmTs2LGMGzeORo0asXHjRn7++Wdq1aplk6+tSCjCKtHR0QqgREdHO7opdpeRkaE8//zzCqAASqVKlZRr166p7wUGBiqAUrlyZQe3NKeMjAylRIkSCqDUqlVLURRFSU1NVVauXKmEhYWpX5PxV4UKFZQPP/xQiYmJyXGvY8eOKc8995zi6emZ43O9e/dWtmzZomRkZCiKoig9evRQ39u1a1eRfs1CCCFcQ0GyhvTYCYt99tlnfPnllwB4e3uzbt06dTWTTqdT55BdunSJ6Ohoh7XTnKioKOLj4wGoWbMmYOi5e+yxxzhw4AAbNmygS5cu6vWRkZG89tprVK5cmTfeeIOoqCjCw8N54IEHaNCgAUuWLFGHoj08PBg6dChHjhxh/fr1dOvWTZ2v9/jjj6v3lOFYIYQQ9ibBTlhk27ZtJnMMvvjiC5o3b25yjXZxwLFjx4qsbZbQzq8zBjsjnU7H/fffz9atW/nnn38YOHCg+t7du3eZPn06FSpUoGvXrvz222/qewEBAbz++utcuHCBZcuWmV0cMWjQILX0y6pVq5ymgKUQQgjXJMFO5OvSpUsMHjxYDSXjxo3jiSeeyHGddp6dsy2gOHfunHqcPdhpGRdRHD9+nGeffVZd0q5dVBEaGspHH33E5cuXmTlzJhUrVsz1foGBgeoqraioKJOaf0IIIYStSbATeUpISODBBx/kxo0bANx33328//77Zq/V9ljZZQFFRgZcuwb79sHPP8OqVbB+PWzbBgcOwNmzEBUFCQmQbaFCXj125tSrV48lS5Zw/vx5xo4dS3BwMM2bN+ebb77h3LlzjB07loCAAIuaLcOxQgghiopOUWSpnjViYmIIDAwkOjra4m/0xYWiKDzxxBNqGKlevTr//vsvwcHBZq+/desWpUuXBgwbPxeodyopCa5cMfz677+sY+25q1fB0qFMvR78/MDfH/z9OXv9Opdu3yYW6PzwwwQ1awY1akD16ob/Zttw2pYSExMpW7YssbGxBAQEEBUVhbe3t92eJ4QQwrUUJGtIsLOSKwe7Dz/8kNdffx2AEiVK8M8//+RbZLdChQpcvXqVUqVKcePGjdyL/ioK/PEHfPopbN8Ot2/buvkFU7KkadDT/rdSJXBzs+r2Q4cOZfny5QD8+OOPJvP4itrFixeZM2cO3bp1o3///g5rhxBCCMtIsCsC8+fPZ/78+aSnp3P69GmXC3a///47ffr0UeeWWRpGevXqxaZNmwBDnbcc29TEx8PXX8O8eXDihOUNCgmBihWzflWqBL6+EBcHsbFZ/zX+0rxWMn8Vet6BhwdUrQr16kGTJoZfYWFQrRpYuFvF77//zv333w/A4MGDWb16dWFbY7UBAwbw888/AzB58mSmTZtm9a4bQggh7EeCXRFyxR67s2fP0rJlS+7evQvAlClTePvtty367Lhx45g9ezZg2NXhvvvuM7xx8SLMnw+LF0PmfVVlykCtWoawpg1vxgBXvjxkriwtDOMQsS/wQMeOrJo9G86dg/PnDb+Mx5cvG+bxWSogwBDwwsKyAl+DBmBmmDUtLY0KFSpw48YNvL29iYqKcsifl6SkJEqWLElSUpJ67oknnuDLL79UV+8KIYRwLgXJGrJXrDARGxvLgAED1FA3YMAApk6davHntUO1R48c4T4vL/j4Y1i3Lmdo6tgRRo2CAQPA3X5/FI0rYhOA4AYNoEULw6/sUlIMAdQY9LKHv8w6eKqYGPj7b8MvIzc3qFs3q1evSRNo3hz34GAeeeQR5s+fT1JSEuvWrePpp5+215ecq+3bt5uEOoAVK1bw33//sXbtWkqWLFnkbRJCCGf3zDPPcPfuXdatW+fopuRLgp0wMW7cOI4fPw4YVoYuX75c3RbLEo0aNcILeAwY9N57OefOeXrCkCGGQJfLli62ZvGKWE9PQ8+hua1jMjLgwgU4eBAOHTL89+BBQy+fVno6HDtm+LViRdb5OnV4s3p10oB/gO+++cYhwc44TA6Gf6hWrVpFYmIif/31F+3atWP9+vUmW60JIURxFB4eTteuXblz5w5BBVgcFxERQbVq1Thw4ABNmjRRz3/88cfFZltICXbChLEAr4+PDz/99FPBhgsjIwn74QcuAWXANNSVKwcvvwwvvQSZu1UUFW2wq1GjRuFuotcbFlPUqAEPPZR1/vbtrKBn/O/x45Caavr5U6coe+oUCzNfxm/eTEr79nh27AitW0ObNoYhZzszBjudTseHH37Iyy+/TL9+/bh+/TonT56kTZs2/PLLL7Rq1crubRFCiOIiMDDQ0U2wnN02NrtHuNpescb9XuvUqVOwD06dqiju7opiWO+q/spo2VJRvvlGUZKT7dJeSzz11FPqfq1Hjhyx/wOTkxXl4EFFWbZMUUaOVJTWrRXF0zPH702OX5UrK8ojjyjK7NmKsnOnoiQl2bRZV69eVX8fWrRooZ4/f/68UrduXfU9Hx8f5ccff7Tps4UQriU9PV2ZPn26UrVqVcXb21tp3Lixsnr1aiUjI0Pp3r270rNnT3XP7Fu3bikVK1ZU3nzzTUVRFGXr1q0KoPz6669Ko0aNFC8vL6V169Y5/n3++++/lQ4dOije3t5KpUqVlFdffVWJi4tT309KSlJef/11pVKlSoqnp6dSo0YN5YsvvlAuXLiQYx/voUOHKoqiKBs2bFDat2+vBAYGKsHBwUrfvn2Vs2fPqvfM/rnOnTsriqIoQ4cOVQYMGGDy7FdffVUJCQlRvLy8lPbt2yt79uxR3zd+jX/88YfSvHlzxcfHR2nbtq1y8uTJQv1+FyRrSLCzkqsFOw8PDwVQmjZtavmHVq0yCShpOp2yEpQ2oJw9c8Z+jbVQ27Zt1b+k8fHxjmlEUpKi/POP8t/rryvfgnI+v5AHiuLtrSidOyvK5MmK8vvvihITY1UTvv76a/X3YdKkSSbv3b59W+ncubP6vk6nU+bMmWPV84QQruvdd99V6tatq2zcuFE5d+6csnTpUsXLy0sJDw9X/vvvP6VkyZLK3LlzFUVRlIcfflhp1aqVkpqaqihKVuipV6+esmnTJuXw4cPKAw88oFStWlVJSUlRFEVRzp49q5QoUUKZM2eOcvr0aWXHjh1K06ZNlWeeeUZtwyOPPKKEhoYqP/74o3Lu3Dnljz/+UL777jslLS1NWbNmjQIop06dUq5evarcvXtXURRF+eGHH5Q1a9YoZ86cUQ4cOKD069dPadSokZKenq4oiqLs2bNHDWRXr15Vbt26pShKzmA3cuRIpUKFCsr69euVY8eOKUOHDlVKliypXm/8Glu3bq2Eh4crx44dUzp27Ki0a9euUL/fEuyKkCsFu5SUFPUbe4cOHSz70H//KUrJkllhZNQo5cNRo9T7rF271q5ttkRISIgCKBUrVnR0U5SMjAylQYMGCqCUBeXa558ryoQJitKli6L4+uYd9NzcFKVFC0UZM0ZR1q5VlBs3CvRsbc9leHh4jveTkpKUJ5980uSn1VdffVVJS0uz0VcvhHAFSUlJiq+vr7Jz506T888//7wyZMgQRVEU5fvvv1e8vb2VCRMmKCVKlFBOnz6tXmcMPd9995167tatW4qPj4+yatUq9V4vvviiyf3//vtvRa/XK4mJicqpU6cUQNm8ebPZNhqfcefOnTy/lhs3bpiM5hh7+w4cOGBynTbYxcXFKR4eHsqKFSvU91NSUpQKFSooH3zwgcnz//jjD/Wa3377TQGUxMTEPNtkTkGyhsyxE6qEhAT12NfXN/8PZGTAM8/AnTuG148+CnPmUHn1asNKWAx7xj744IO2b6yFYmJi1O3QLNlKzN50Oh1Dhgxh8uTJRAFLb99mwowZhjfT0gyLLv75B3buNKy2vXAh68Pp6bB3r+HXnDmGc/XqQadOhhXGHTtC5cpmn6soijq/rkSJErRt2zbHNV5eXixfvpxq1arxzjvvAPDJJ5/g6enJrFmzbPZ7IITIQ4sWhq0Ti1q5coZ/Wyxw9uxZEhISsspZZUpJSaFp5qK4hx9+mLVr1/L+++/z2WefUcvMojTtv0PBwcHUqVOHE5n1TQ8dOsThw4dZoVmEpigKGRkZXLhwgSNHjuDm5kbnzp0L9GWeOXOGKVOmsHv3bm7evKnWar106VK+BfiNzp07R2pqKu3bt1fPeXh40KpVK7X9Ro0bN1aPy2fOo75+/TqVc/m32hYk2AlVgYPdvHmG3SPAUHPus89Ap6NRo0bqJXbZM7YAjKVOwDmCHaAGOzDsHTthwgTDG+7uWXXxXnrJcO6//7JKqmzbZgh+WidOGH59/rnhddWq8OCD8MgjhgUZmYWHjxw5QlRUFABdu3bF09PTbNt0Oh3Tpk2jWrVqDBs2jPT0dD766CMaNWrE0KFDbfnbIIQw59o1wzaKTiwuLg4wLLarWLGiyXvGepgJCQns27cPNzc3zpw5U6hnvPTSS4wcOTLHe5UrVzZZFFcQ/fr1o0qVKixevJgKFSqQkZFBw4YNSUlJKdT98uPh4aEeGwvBZxSkXmohSLATqnhNnbZ8g93Ro2AMJABffWXYlgtDgPL09CQlJYWjR4/ao6kWs8mKWBurXr06rVu3Zvfu3Rw+fJhjx47RoEED8xdXqmQoDzNkiOH1rVuwY4ch5P39N+zbZ+jJM4qIgLlzDb8qVzYEvEcf5fc//1Qv6dmzZ75tfPbZZ0lKSuJ///sfAC+++CK1atWiXbt2hfuihRCWyb5bjxM+t379+nh5eXHp0qVce8zGjRuHXq9nw4YN9OnTh759+9KtWzeTa/755x+15+rOnTucPn2aevXqAdCsWTOOHz+e6w/kjRo1IiMjg7/++osePXrkeN/4w2u65t/HW7ducerUKRYvXkzHjh0BQ23P/D6XXY0aNfD09GTHjh1UqVIFgNTUVP79919Gjx6d6+eKTIEHeoUJV5pjd+jQIXVu1bBhw3K/MClJUcLCsuZ+jR6d45KwsDAFUNzc3JQkG6/uLIjp06erX9P333/vsHZk9/HHH6vteuONNwp/o9hYRdm8WVGmTFGUrl0VxcPD7Py8Kz4+ynRQwkA5eeKExbd/+eWX1XaWLVtWuXTpUuHbKoRwGW+88YZSqlQpZdmyZcrZs2eVffv2KfPmzVOWLVum/Prrr4qnp6eyb98+RVEUZeLEiUqlSpWU27dvK4qSNf+sQYMGyh9//KEcOXJE6d+/v1K5cmUlObOCwqFDhxQfHx9lxIgRyoEDB5TTp08r69atU0aMGKG24ZlnnlFCQ0OVtWvXKufPn1e2bt2qztH777//FJ1Opyxbtky5fv26Ehsbq6SnpyulSpVSnnzySeXMmTPKli1blJYtW5rMB09NTVV8fHyUd999V7l27Zq66CL74olRo0YpFSpUUDZs2GCyeCL716id43fgwAEFUC5cuFDg329ZPFGEXCnY7dq1S/0mPmrUqNwvfP31rNDQoIGimJkI+sQTT6j3OnTokP0anY/nnntObcf+/fsd1o7srl69quj1egVQqlevrpYFsNqdO4qydKmi3H+/2fIzCigZtWsbVtpaEPBSUlKUrl27qr+HTZs2NSk3IIS4N2VkZChz585V6tSpo3h4eCghISFKr169lPDwcKVs2bLK9OnT1WtTUlKU5s2bK4888oiiKFmh55dfflEaNGigeHp6Kq1atcrxvWLPnj3Kfffdp/j5+SklSpRQGjdurLz33nvq+4mJicqYMWOU8uXLK56enkrNmjWVJUuWqO9PmzZNKVeunKLT6dRyJ5s3b1bq1auneHl5KY0bN1bCw8NzLPRbvHixEhoaquj1+lzLnSQmJiqvvvqqUrp06TzLnUiwK4ZcKdht2bJF/QY+ceJE8xeFhyuKTmcICZ6ehnptZrz//vvqvb755hs7tjpvnTp1UtvhbP+PevToobbtn3/+sf0Dbt5UlMWLlRtNmyppua20bddOUb780tDzl+ttbirVq1dX2/rwww/bLogKIe45lq5YFVkKkjUs3ytKuDztHLsSJUrkvCA6Gp5+2hAJAN591zDR3wyTPWMdOM/OOMeuTJkyBdtFowg8/vjj6vG3335r+weUKgUvvMCMrl0pDwwHrjdooC6oAAyrb59/3jC/5vnnDa+zbZtTqlQpfv75Z/z9/QFYvXq1umpWCCGEc5FgJ1T5rop95RW4dMlw3LkzjB2b6720wc5RK2MTEhKIjIwEnGdFrNbAgQPVibqrVq0iLS3NLs/ZtGkTN4DFej3u27ZBZKShXIp2aX98PCxZAu3bQ/368OGHkLmKFqBBgwZ8++236qquqVOnsmbNGru0VwghROFJsBOqPIPd99/DN98YjgMCDKtg3dxyvVflypXVHh5H9didP39ePXbGYBcUFETfvn0BiIqKIjw83ObPiIyMVH//W7ZsSXBwsKF3bvRoOHwY9uwxlFbR9maePAmvv25YkTtoEGzZAorCAw88wAxjzT3g6aef5uDBgzZvsxDCtXXp0gVFUQgKCnJ0U1ySBDuhyrXcyZUrMHx41uv58yFziXdudDqd2mt38eJFYmJibNpWSzhjqZPshhjLmGCf4djNmzerxznKnOh00LIlLFwIV68awnqnTlnvp6XB2rXQowc0aAALFvD6//7Hk08+CRh+EOjfv79aH08IIYTjSbATKm2PnTrHztzuEk88YdH9tIWKj2UvrFsEtMHOGXvsAB544AH8/PwA+PHHH0lKSrLp/Y27TUA+9et8fQ3zJ//6C06fhokTIbNKOmAogjxiBLpKlVgaGMhDmdXUL1++zMMPP4ySbV6eEEIIx5BgJ1Rmh2Jz2V3CEo5eQFEcgp2Pjw+DBg0CIDo6mg0bNtjs3hkZGWqPnb+/P61bt7bsg7VqwfTpcPGiYQg+s5AnADExuM+fzw+HD7PFy4vewPa//2b//v02a7cQQojCk2AnVDmCXR67S1jC0VuLFYdgB6bDsStXrrTZfQ8dOqTuk9utWzeTrW0s4uEBDz9s2OXi4EF44QXw8VHf7paczHrgBHBr5kxITrZZ24UQQhSOBDuh0s6x8/PwgCefzPpmPXo0dO9eoPs5usfOuE9syZIlDYsGnFT37t0JCQkB4JdffrHZfESLh2EtERYGixcb9q794APDnrSZ6gA9V6+GatUM70VHW/csIYQQhSbBTqi0PXaVv/wSDh0yvGjQADSrIS1VunRpymXuP3jkyJEinYeVnJzMpczSLM7cWweGTaIffvhhAJKSkli3bp1N7mvTYGcUHAyvvQZnz5Kxdi1/u2u2m756Ff7v/wx71E6YYHgthBCiSEmwEypjsOsElFqyxHDS0xNWrABv70Ld09hrd/PmTa5fv26LZlokIiKCjIwMwHlXxGppixXbYjg2Pj5e3dy6WrVqtv89cHND/+CDzH7gAVoBawDFOPcyJgZmzjT06r34Ipw5Y9tnCyGEyJUEO6GKj48nAFgO6CzYXcISjppnV1zm1xm1bduWypUrA4YSJca5cYW1bds2UlJSAENvnc7CBS8F1blzZ/4FBgPfTJpkmIeXWXSZlBTD8G3duvDUU3DqlF3aIIQQIosEO6FKSEjgfUCtUJfP7hKWcNQ8u+IW7PR6vbqIIj09ndWrV1t1v99//109ttkwrBmdNHXvfjp50hDkIiIMBY6NRY8zMgzFrevXN8zbPHnSbu0RQoh7nQQ7oUpISKBf5rHi45Pv7hKWkB47y9ly71jj/Do3Nze6detm1b3yEhYWRmBgIGDoJVQUxVD/buZMw/Zz771nmJcHhoC3YoUh4D3xhKE2nhBCCJuSYCdUCQkJ+Gce6ypXznd3CUvUr19fHQaUHru8NWrUiPr16wOwY8cOLl68WKj7XL58mROZoal169Z23bbHzc2NDh06AHDjxg1OanvjAgNh0iRDD96MGVCqlOG8osC33xoW5QwZIkO0QghhQxLshCo+Lo4SxhclSuR1qcVKlChB9erVAcPuE8YFDfZmLHXi5+dHmTJliuSZ1tLpdDZZRJHnNmJ2oB2O/euvv3Je4O9vWCV74QK8/75pwPvuO0PAe/FFQykVIYQQVpFgJ1RpCQmoxStsFOwga55dfHw8ERERNrtvbtLS0rhw4QJgWBFrr4UD9vDYY4+px2+99RbvvPMOyQUs/GuXMid56Ny5s3q8bdu23C/09zeUQ4mIMAzVli5tOJ+ebpibV6uWYW7e7dv2bbAQQrgwCXZCpdPUsbNHsIOimWd36dIl0tLSgOIzDGtUo0YNBgwYABhq8U2ZMoXGjRvz559/WvT59PR0tccuMDCQli1b2q2tRs2aNVP3Fv7rr7/yr1fo52cIcBcuwDvvZC2ySEqCDz+E6tUNW5ppCmYLIYSwjAQ7AWD4Zqz9RmrDYKddQFEU8+yK4/w6rRUrVvDaa6/hlrlw5fTp03Tv3p2nnnqKqKioPD974MABbmf2eHXv3h13bQFhO/Hw8KB9+/YAREZGqsPg+fLzg8mT4dw5GDcOvLwM56Oj4Y03oEYNWLQIMkO6EEKI/EmwEwCkpqbirZ3/Vox77Ip7sCtRogQffPAB+/fvp127dur5b775hrp167Jw4cJc5yoW9TCskXaeXZ7DseaULg2zZhkKGb/wAugz/1mKioKXXoImTUBTvkUIIUTuJNgJwLAi1iTK2TDY1a5dW92AXnrsLNe4cWP+/vtvFi9eTMmSJQG4e/cuL7/8Mu3atePgwYM5PuOoYKedZ2d2AYUlQkMNc+2OHYOHHso6f+wY3H+/4ZcD9hwWQojiRIKdAOwb7Dw8PKhbty4Ap06dUndEsBftUGBxDnZgKFz8wgsvcOrUKYYOHaqe3717N82bN2fs2LHExsYCEBsby86dOwHD112tWrUia2fLli3xztx2rtDBzqhuXfjhB9i+HVq1yjr/+++GXVBefBGuXbPuGUII4aIk2AnAsGLVXsEOsubZpaWlccrOdcuMPXbe3t5UqFDBrs8qKiEhISxbtoytW7eqITkjI4M5c+ZQr1491qxZQ3h4OKmpqUDR9tYBeHl50bZtWwAuXrxY6Bp8Jtq3h127DDXvMrdbIyPD0KtXu7Zh+LaQPyScOnWKefPmFen+xUIIURQk2AnAvj12UHRbi2VkZKg9dtWrV0evd60/4l26dOHQoUO89957ag/ZlStXGDx4MM8++6x6XVEHO7Bynl1u9HpDEeOTJw1Fjv0zS2jHxsJrr0GjRrBxY4FumZ6eTq9evRg1ahQjRoywTTuFEMJJuNZ3PVFo9g52RbW12JUrV9S6b8V9GDY3np6eTJo0iWPHjtG7d2/1/K1btwBwd3ena9euRd4um8yzy42Pj6HI8dmzMGwYGGsTnj4NvXtD//6G9yzwzz//qD2KxqFrIYRwFRLsBGD/odii6rFzlYUTlqhevTq//fYbP/zwg8mQc7t27Qgw1oYrQm3atMHT0xOwLNgpipJ/zbvsypQxlEDZuxc0K4b55RfDDhaTJuVb/+63335TjyMjI4mXenlCCBciwU4A9u+xq1KlCn5+foB9e+zupWAHhm3IHnroIU6ePMn//d//0aVLFz788EOHtMXHx4dWmYsdzp49S2RkZK7XJicnM2jQICpUqMCOHTsK/rBmzQyLK1asAGOoTUkxDNc2bAjr1+f6UW2wAzh//nzBny+EEE5Kgp0A7B/sdDqd2msXERGhruS0NVdaEVsQ/v7+vP/++2zdulUNV45g6Ty7mTNnsm7dOq5du8bs2bML9zCdDh5/HE6dgokTIbO3kIgI6NsXHn4YsoXLy5cvc/jwYZNzFhdUzrR8+XKqVq3KvHnzCtduIYSwIwl2ArD/UCyYzrM7duyYze8P916PnbOxZJ7dqVOneO+999TXJ06csO6hfn6GLcgOHwbt3MIffjCUTvnkE8N+tOTsrQPTPzOWeOedd7h48SJTpkwp+FCyEELYmQQ7Adi/xw6KZp6d8Zu0u7s7oaGhdnmGyF27du3UrdDM9dhlZGTw4osvmtQyPHPmjFqmxSp16sCWLbB8uWE3CzCsnh05Etq0gf37rQ526enpREREABAdHc3du3etb7cQQtiQBDsBFE2ws/fKWEVR1G/S1apVK5J9UoUpPz8/WrRoAcDx48e5ceOGyftLly7NEfjS0tIKPByaK50OnnrKMDz7wgtZ5/fuRWnZkp4bNuAHBAUFqW8VJNhFRkaSptm79tKlSzZotBBC2I4EOwG4Ro9dVFSUusJRhmEdJ7d5dlFRUbz22mvq644dO6rHJ0+etG0jgoMNhYz//tuwWhbQZWTwano6x4F3mjdXt2krSLDLXnjZJoWYhRDChiTYCaBo5tiFhIRQpkwZwD49djK/zjlo59lpg92YMWO4c+cOAE888YRJcWCr59nlpkMH2L8fpk8nJXOIOBR4ZcsW1ioKlTD0uhlrH+Yne5CTHjshhLORYCeAoumxg6zh2Bs3bth8OycJds6hQ4cO6DILCBsXUGzcuJGVK1cCEBwczOzZs9Wt0cCOwQ7A0xNlwgS6ly3LBs3pznfvcgx4TlG4YGHJEwl2QghnJ8FOAKbBTnFzyyodYWPa4Vhb99rdq6VOnE1gYCBNmjQB4PDhw/z333+8/PLL6vuzZs2iTJky1K5dWw2ANh+Kzeb48eNsj4ykD/Buo0ZQrhwAAcAXQKknngALhlUl2AkhnJ0EOwGYDsUqvr5ZWzbZmHYBha3n2Wl77GrUqGHTe4uCMQ7HKorCgAED1JWkXbp04ZlnngEMBY2rVasGGIKdtaVDTp8+net8uV9//VU99n/+eThxgjOaOX4hBw4YCht//jnk0Q6ZYyeEcHYS7ASQbSjW19duz7Fnj53xm7per6dq1ao2vbcoGO08u/379wPg5eXF559/rvbSAepwbGxsLFeuXCnUs65fv85TTz1FnTp1qF27Np9//nmOa7RlTvr27QtBQVyfMYP7gcvGN+LiYPhw6NEDLlww+yxb9tjduHGDxMTEQn9eCCHMkWAHzJkzhwYNGlC/fn1Gjhx5TxYdNQl2mVt/2UODzBWKYNseO0VROHPmDACVK1fGy8vLZvcWBadd8Wo0efJkateubXKuXr166nFBh2MzMjJYtGgRderU4ZtvvgEMfw6GDx9usivE7du32blzJwB16tRRh+lr1qzJ70BDYEOlSlk3/vNPaNQI5s+HjAz1tKIoOYLd1atXTWryWWrt2rVUqVKF0NBQbt68WeDPCyFEbu75YHfjxg0+/fRT9u3bx5EjR9i3bx///POPo5tV5LTBTmfHYOfn56cOvx09epQMzTdOa9y+fZvo6GhA5tc5g1KlSpn0ztavX5/XX389x3XaYFeQBRSHDx+mQ4cOvPTSS2qRYB8fH/X9UaNGqXvm/v7776Rn7jzRt29f9ZoyZcpQokQJYoCR3t6waRNUrmx4Mz4eXnkFunWDzLmbN2/ezNHDpihKgXsaz5w5w9ChQ0lMTOTWrVts3LixQJ8XQoi83PPBDgwFUpOSkkhNTSU1NVUtyXEvSYqLw9jHpbPTilgj4zy7+Ph4m81RkhWxzqd///6AYZ/gRYsW4WlmQU5BV8amp6fz2muv0axZM3bt2qWef+qpp4iIiGDKlCnquddff51333035zBsJp1Op/5ZiYiIIK1rVzh61DAca/TXX9C4MXzyCRdzGZ4tyHBsYmIigwcPNtkr2djTLIQQtuD0wW7btm3069ePChUqoNPpWLduXY5r5s+fT9WqVfH29qZ169bs2bPH4vuHhIQwfvx4KleuTIUKFejRo8e9OfE+s7AvYLdSJ0b2mGcnK2Kdz+TJk5k9ezZbtmyhffv2Zq8p6FDsvHnzmDVrltoDV6dOHf7880+WL19OmTJlePvtt032oX3zzTf57rvvAAgICKBDhw4m9zP+WUlLSzMENH9/+Owzw9ZkxnmaCQkwciShw4ZRMfNz2h/+CvLDyciRIzl8+LDJudOnT1v8eSGEyI/TB7v4+HjCwsKYP3++2fdXrVrF2LFjmTp1Kvv37ycsLIxevXqZ1Ehr0qQJDRs2zPErMjKSO3fu8OuvvxIREcGVK1fYuXOn2T0uXV4RBjt7rIyVHjvn4+Pjw5gxY+jatWuu1wQHB6shyZIeu59//lk9fvvttzl06FCO+0+aNIlZs2apr40hsGfPnjl6DbV/VkxW1HbrBkeOgKaIctnDhzkCPILpHEJLe+yWL1/OF198ARh+b/R6wz+/0mMnhLAlpw92vXv35t1332XgwIFm3589ezbDhg3j2WefpX79+ixcuBBfX1+WLFmiXnPw4EGOHj2a41eFChX4448/qFmzJsHBwfj4+NC3b997co6dLiEh60UR9tjZI9jdkz2uxZhxOPbatWvqfDlz0tPT+ffffwEIDQ1lypQpuS6SGTduHJ9++qnJOe0wrFGuwQ4Mi4g+/RQ2b4aKhr66ksAq4N2ICIIyL7Mk2J06dcqklt/ChQvVldunT5++JxdsCSHsw+mDXV5SUlLYt28fPXr0UM/p9Xp69OhhMv8mL6GhoezcuZOkpCTS09MJDw+nTp06uV6fnJxMTEyMya/iTlEU9NpJ4XYOdrVr18bDwwOw3VCs9pty9erVbXJPUTQsXUBx7NgxdS/g1q1b53vfESNGsHjxYnx8fKhRowaDBg3KcU2ewc6oRw84coRtFSuqp+ru28dhoBuWBbvFixeTkPnD0wsvvMDTTz+trhCOjY0lKioq33sIIYQlinWwu3nzJunp6ZQtW9bkfNmyZbl27ZpF92jTpg19+vShadOmNG7cmBo1aqiTvs2ZMWMGgYGB6q/Q0FCrvgZnkJKSgo+2x8DOwc7T01MNzydPnixUuYjsjN+UK1asiK8d6/AJ27N0nt3u3bvV4zZt2lh07xdeeIFr165x6tQpAgICcrxvUbADKFmSUSEhDAHuZJ4KBbYAQ/bsgXzq0RkLNAO88cYbANSqVUs9J8OxQghbKdbBzlbee+89Tpw4wbFjx5g3b55JAdXsJk6cSHR0tPrr8uXLuV5bXBTVPrFaxuHYtLQ0qyePR0dHc+PGDUDm1xVHlq6M1QY7S3rsjAICAnBzczP7XoUKFdTh3DyDHYZFEt8BfStVgu7d1fND79xBadECDhzI9bNXr15Vj8uXLw9gUtNPFlAIIWylWAe70qVL4+bmlmMYIyoqinKZe0HampeXFwEBASa/ijvtdmJAkQQ7Wy6gkBWxxZulQ7HGYOfm5kazZs1s8my9Xq/OyTx//nyudRVjY2O5c8fQV+dZowZs2sSCOnVIynxfd/w4tG4NH35oUtTYKDIyEjAsFjEGSW2wkx47IYStFOtg5+npSfPmzdmyZYt6LiMjgy1bttC2bVsHtqx4cWSPHVg/z06CXfFWqVIlSmT+mcttKDYmJoZjx44B0LhxY5sOtxv/zCQnJ+dabFhb0qRKlSqg13OgY0eaAfuNb6SmwuuvQ69eoOmhUxRF7bEz9taB6VCs9NgJIWzF6YNdXFwcBw8e5ODBgwBcuHCBgwcPqhOWx44dy+LFi/nqq684ceIEL7/8MvHx8Tz77LMObHXx4ohgZ8seO1kRW7zp9Xp1zuX58+dJSkrKcc3evXvVlaOWzq+zlCXz7HIEOwxb150A2gCnBw0C4xSOP/4wFDX+9VcA7t69S3JyMmAY+jWqXLmyWn5FeuyEELbi9MFu7969NG3alKZNmwKGINe0aVO1wvyjjz7KrFmzmDJlCk2aNOHgwYNs3Lgxx4IKkbv4+HhM+j+KINhVqVJF7aWxtsdOatgVf8bh2IyMDLMhp7Dz6yxhTbADSAU2du5sCHTGHrmbN6FfPxg5kmuahRPaHjs3Nzf1B5EzZ87kub3etGnTaNasmbrnrRBC5Mbpg12XLl1QFCXHr2XLlqnXvPLKK1y8eJHk5GR2795t83/4XZ0jeuz0er06HHvhwgXi4uIKfS/psSv+8lsZ68zBDjJLnnTrBocPg3ZV/SefUGnwYIxfnTbYQdY8u+Tk5FwXYt28eZO33nqLAwcO8MEHH1j8dQkh7k1OH+yE/Tki2IHpPDvj/KnCMH4zLlOmjEssZrkX5bUyVlEUtWh4YGCgyaIDW9D+MFCQYGf8L2hq2ZUuDevWwfz5kLlIwv/8efYCLwLlsy3qsqTkSUREhDoM/d9//1n0NQkh7l0S7App/vz51K9fn5YtWzq6KVZzVLCzxTy7hIQEdcWhDMMWX3mtjL106ZK68r1Vq1bqVly2UrlyZdzd3QHThTha2mBnrF1ZsWJFtTSSSZFinQ7+9z/4919o0AAAX+Bz4KGVK+HOHfVSS0qeaHvytFslCiGEORLsCmnEiBEcP35c3eKoOHNEuRMw7bHbtGlToe5x/vx59ViCXfFVs2ZNtdZc9qHYwhQmLgh3d3eqVasGGHrszG3vZQx25cqVw9vbGzCUPjKWVdIGP1WjRvDvv2xv3Fg9VWnPHmjWDPbsASwreaINjTdu3JDtx4QQeZJgJxzWY9e6dWuCgoIA+P777/nzzz8LfA+ZX+caPD091f9/p06dMllIYM/5dUbGHwri4+Nz1MVMTk5Wy5Voh18ha57dtWvX1JWvJnx8+LRePR4EbhvPRURAhw4wdy61ND+MWNJjl5SUpG6rJoQQ5kiwEw4Ldn5+fsycOVN9PXz4cLOlLvIiK2Jdh3E4NjEx0aQHzDi/DgxDsfaQ1wIKbbDKHuy0r3Ob/xYZGclPQBMg3RhMU1NhzBjKjxhBxcyafLkFu+x70cpwrBAiLxLshMOGYsGwl2e7du0Aw1DU+++/X6DPS7BzHeZWxqamprJ/v6EEcPXq1QkJCbHLs/MKduYWThjlWBlrhrG3LyYwELe//zYUMc6k++kn9qSm0grD6vDU1NQcn8++Wta4fZ4QQphToGAnQwCuSdtjp+h04ONTZM/W6/V8/vnn6uT1GTNmcOrUKYs/L8HOdZhbGXv48GG1F9ce8+uMbBHszM6zA9NdJzw8YOZMQ/Hi4GAAKqSmsh0YmZ7OBc2cUSPpsRNCFESBgl2nTp3U440bN9q8McIxtMEuw9s7q4J+EWnYsCHjx48HICUlheHDh1s8Qdz4TbhkyZIEZ36jFMWTuZWxRTG/DkznZ2ZfGWtNj11sbKz6A7FJDbu+feHgQcjsrfYAZgMBQ4fCbXU2HqmpqWowNJIeOyFEXgoU7LQTmidNmmTyXvfu3W3TIlHkTHrsbLgHZ0G8+eab6srE8PBwvv7663w/oy3qKr11xZ+2x844FKudX2fPYFe1alW1jEpBeuzM1rLT0Iay7MWJCQ2F8HCO9u2rniq3ezc0bQq7dgFw5cqVHD/kmAt2slJWCGFUoGCn0/TkZP+H5Lbmp0xRvJjMsSvC+XVavr6+zJ8/X309btw4bt26lednIiIi1B82JNgVfwEBAepeqtl77Dw9PWnSpIndnu3l5aX2vp05c8bk3zdreuyMNRbBTLAD8PAgbvJk+gA3s24EnTrBrFlcNjO8m30o9sUXX6Rs2bJ8++23uX15BWbc7aKwZYiEEI5ToGB348YN1q1bx4ULF0xCHpDjtSg+TFbFOijYAfTu3ZtHHnkEMHxjeV0zydwcKXXieozDsbdu3eL06dPqStGmTZvilbmTg70Yd4GIjo7mu+++U88bg11QUFCOnU1Kliyp7nlsbo6dtsfOGFrNPXcDhlWzRwMDDSfT0uC116g2ZgxB2a7X9thFRkayePFibty4wdChQ9m+fXv+X6gFpk+fzttvv82DDz7IHU1BZSGE8ytQsBs7diy//PILjz32GOfPn6dt27YMHTqU9957L9/eFVfjSjtPJGp67PR+fg5ty9y5c9VvnkuWLGHbtm25XisLJ1yPdp7d8uXL1eOi2P/5mWeeUY+HDRvG8ePHSU9PV4f7s/fWgeEHWmOv3aVLl3KMZOQ5FJupVKlSBAcHcwXo5+8PEyeq71U6cIC9QJjmem2PnbaXMC0tjYceeijXPWcL4tChQ4Ch9ExuhZOFEM6pQMHuxRdf5Msvv2T37t3cvXuXb7/9lsGDB6PX6+nQoYO92uiUXGnnibTYWPVY5+/vwJYYvvlNnz5dfT18+HBSUlLMXivBzvVo59kVdbAbMmQIQ4cOBQzTEwYPHsyZM2dIS0sDzAc7yBqOTUpK4ubNmybvWRLsIKu3MOK//0h8803YsAFKlQKgBrALGJp5rbbH7sqVKyb3uX79OgMHDiQxMTHvLzYf2nbbIigKIYqOVatiq1WrRr9+/Zg4cSIrVqyweeNEEdGUsdE5cCjWaPjw4Woh2hMnTvDhhx+avU6CnevR9thpA0VRBDudTseCBQvUPYxPnDjBo48+qr6fW7DLawGFpcEux9Zi998P+/ZxOnNo1gdYhmG/2WjNzhjaosjGxR/79u3jlVdeyfVZlpBgJ0TxJatihUmwc+QcOyM3Nzc+//xzde/Qd955J8dKRcgKdn5+fpQpU6ZI2yjsQxvsjEqXLk316tWL5Pm+vr788MMP+Gf2XB8+fFh9L78eO8gZ7LSLJ3KbYwdQv3599dhYkJkqVXg8NJTPNNe9CKyKjES5cAEwDXZz5sxR5/stXbqU6OjoXJ+Xl8TERO7evau+lmAnRPEiq2IFeu2wjRMEO4AmTZowevRowFDW5H//+5/Jn7m0tDQiIiIAQ2+dLN5xDeXKlcuxQKF169ZF+v+3du3aLF26NMd5S4Jd9gUUxp6vEiVKqGHRHG3xZW2Jl/NXrvA/YFzp0iRn9sg1VxRo1gw2bDAJdr179+bxxx8HDP8+G1cWF1T2unkS7IQoXmRVrECXkJD1wkmCHcBbb71FaGgoAJs3bzZZqXjp0iV17pOsiHUdOp0uR69dUQzDZvfQQw8xZswYk3OWDMVm3zXFZNeJPLRo0UIdSt2VWcMuLi5OXZG6r0ED3uzRA2O/te7uXejbl/t27MD4L2/FihVp0KCBes9jx46px6mpqYwaNYqRI0eSnJycZ1sk2AlRvNl8Vezo0aNZunQp+/bty/cfEOF4iqLglrllE+BUwc7Pz49PP/1UfT169Gj1G53Mr3NdzhDsAGbOnKnuYxwQEECdOnXMXte0aVM1lO3cuVM9n5CQoA6H5hfs/Pz8aNy4MQBHjx4lNjbWJFCFhoaSVKcOLYCfjCcVhecuX+Y3oEZQEL6+vrkGu9WrVzNv3jw++eQT5s6dm2dbtMPHIMFOiOKmQMFuzJgx+a6K7datG1euXGHmzJk0a9aMhg0b8uijj/Lee+/xyy+/2OWLEIWXnJyMyV4TThTsAPr378/AgQMBw4q/iZmlICTYuS7tylhAXUhT1Dw8PNiwYQOzZs3it99+yzFEbOTv76+GsiNHjqhhzpIadlrG4diMjAz27t1rMl+vcuXKhISEEA0MBI4/9RRKZpjsDWyLj4d9+0yC3fHjx9Vj7fDup59+Smpqaq7tyN5jd/XqVbV3XAjh/AoU7LIztyq2f//+TJ48me+//55jx46xd+9eXnvtNcqXL8+ff/5pk0YL2zHZdQLAQVuK5WXevHn4ZdbX+/zzz9m5c6cEOxem7bGrW7cuQUFBDmtLQEAA48aNy7ecU/v27QFDD7gxRFm6ItYo+zw7bbALDQ1VFwgpwK7Onbnz/fcYK9pVSE2F9u0pt2mT+vul7bFTF2RgWHCxdu3aXNuRPdhlZGTk6MUTQjivQge748eP8+abb5qsnjLH29ubFi1a8NxzzzFnzpzCPk7YicmuE+B0PXYAlSpV4p133lFfv/TSSyYTwyXYuZYmTZqoc3a1JZacmTHYAezYsQOwPthph0CNPXZGN27cIKJaNZoC6uBvcjK6Z55hkY8P7hgCXHR0NBkZGRw8eNDkWXkNx5oLcTIcK0TxUehgN2PGDI4ePWr2p+mkpCR1E2/h3IpDsAN45ZVXaNasGWCYg7Rx40bA8IODJcNcovioXLkyn3/+Oc8++yxTpkxxdHMsYotgV6tWLUqWLAmY77HTBrvr16/z33//EQl0Bfa1aKG+9/DVq/wBhGCoxXfmzBnitSWNMCzQ2LNnj9l2ZO+xAwl2QhQnhQ52//zzDyNHjjT7nre3N8OGDWPGjBmFbpizc5UtxYpLsHN3d+fzzz9XJ6kbVa9ePcc5UfwNGzaMJUuWULFiRUc3xSKVK1emUqVKAOzevZu0tLQCz7HT6/XqQpHr16+bbKdXuXJlk1qNN27cUEudpACH//c/WLQIPD0B6AzsA6J++40DBw6on9MWQv7444/NtkOCnRDFW6G/I/733395DoENHz6cn3/+ubC3d3qusqVYjjl2ThrswFASIntFfRmGFc7C2GsXHx/PoUOHTIY0LemxA2jbtq16fCGzCLG/vz+BgYFme+yMKlasCMOGwV9/kRwcDEAo0GfGDPSardnee+89SmVuVfb999/n2JIMZChWiOKu0MEuODjY7E92Rq1atTK7W4BwLsWlx87onXfeMen9kGAnnEX24diCDsWC6Tw7I2MB5KCgINzd3QHTHjtA7S2kTRvubtnC9szzHunpPLJhA/MAd6Bdu3a89NJLgKHI9xtvvGHyrOTkZLXYvLGGJEiwE6I4KXSw69SpE8uWLcv9xno9Sdr6aMIpFbdgFxAQYFLbrlu3bg5sjRBZcgt2Xl5eFq/sNVfaxRiw9Ho9pUuXBvIIdkCZsDAGBQayQHOPV4FtHh6U1+t59dVXCczcg/arr75S56sCXLt2TT1u1qyZOs1Bgp0QxUehg9348eNZvHgxixYtMvv+rl27imx/R1F4xWko1mjgwIH89ddf/PLLL/Tp08fRzRECgMaNG6t7tWqDXYUKFSzemScoKChHgWbtlmXGeXbXr19Xw5a/v79JjT2dTkfthg0ZATwPGMvEt01NRdeiBeUuXWLWrFnq9S+99BKxsbGA6TBs5cqV1Z5GCXZCFB+FDnbNmzdnwYIF/O9//+O+++5j3bp1XLp0idu3b/PTTz/xf//3f+q+hcJ5FbceO6NOnTrxwAMPyFZ2wmm4u7urix+uXLmi7sZj6TCsUfbhWO2QqHGeXUpKijoHT9tbZ1S/fn0AlgCdAHUm3ZUr0KkTz+v1am/3pUuXmDRpEpBzJa/x2devX5edhIQoJqxaTjhs2DDCw8OJjo5m0KBBVKtWjZCQEAYOHEijRo1y7LUonE9xDXZCOCPtcKxRQYOddgEFmPbYaRdQpKenA+aDnXYHij1Ac+CmcUeP5GR0zz/PjxUq4O/tDRh2o9i7d2+uwQ4wGfoVQjgvd2tv0KFDB/bs2cPJkyfZv38/CQkJNGzY0OwkYOF8JNgJYTu2CHZ59dhpS54Y5RfsAKKAmLVrKf3JJ7DAMPsu8JtvOFa9Ok3Pn+cWsGLFCnUoGQxDyNkXUNSoUaNAX4sQouhZHeyM6tatm2OPR+H8tHPs0j08cHNzc2h7hCjO2rRpg06nQ1EU9VxBC2jXr18fPz8/4uLigNx77IzyGoo1CgwMpFqdOjB/PrRoAS+/DMnJhJ4/z79Af2D9+vUmwdTZeuzu3LnD4MGD8fX1ZfXq1Xhn9jYKIUxJZdd7nLbHLsPHx6FtEaK4CwwMpFGjRibnCtpj5+bmpm6l5u/vb7ri1cIeu/Lly5usxNVu08azz0J4OJQrB0A1DNuS1T19mu3bt5vcQ3tvWy2g+Omnn5g6dWq+21Fmt3r1av78809+/fVXfvnlF5u0RQhXJMHuHqcNdoqvr0PbIoQryD4cW9BgB4a9XF944QVWrFiBl5eXet7SHjudTmcyHGvcjk/Vpg38+y80bw6AP/AT8PCZMwB4eHhQqlSpXGvZ3bp1izfeeKPARegvX77M4MGDmTZtGh9++GGBPnv9+nX12FxhZSGEgQS7e5x2KFaCnRDWs0Wwq1WrFosXL6Zfv34m580Fu9y2XdMOxzZt2jTnBZUqwd9/c7d3b/XUe8BKoFrZsuh0ulyD3YwZM5g+fToPP/wwN27csPCrMpTBSktLA2Dv3r0Wfw4gJiZGPb5582aBPivEvUSCXSG5zF6xmmCn8/NzaFuEcAW2CHa5sXQoFmDw4MGAYXi4V69e5m/o40Pgr7/yfkAAGZmnHgPW3bkD//1H2bJl1d0uLl26pH5sz549gKHsypnMXj5LaPetLcjnALXWHlCgMCnEvcYmiyeSkpI4fPgw169fJyMjw+S9/v372+IRTmfEiBGMGDGCmJgYtYp7cZQaF4dxuYQEOyGsV6VKFSpWrMiVK1fw8vJS92a1hew9dt7e3gRn7g2bXc+ePTl37hxBQUG5XgOg0+u58NhjPLhoESswDMvWi4+HFi1wW7uWKlWqcO7cOc6dO0dGRgZ6vZ7Tp0+rn9fuVpGfgwcPqscXL14kJSUFT09Piz4rPXZCWMbqYLdx40aefvpps3/RdDqdWmtJOCclc+UdgN7f34EtEcI16HQ6Jk+ezGuvvcbIkSPVbblswbhfrHE4s1KlSnkW6bZ0958+ffrw4KJFtAV+BqoDREVBly68Ur8+YzBM27h8+TIlS5YkKipK/Wxee4Znp+2xy8jIICIigtq1a1v0WW2wkx47IXJn9b84r776Kg8//DBXr14lIyPD5JeEOueXIcFOCJsbPnw40dHRvPfeeza9r06nM+m1y20YtqC6d++Oh4cHx4CWQES1aoY3UlIYffAgH2L4ZnHixIkcQ6iWBrtr166ZBEKAs2fPWtxG7VCs9NgJkTurg11UVBRjx46lbNmytmiPKGK6+PisYylOLITN2LKnTks7z85Wwc7Pz4/OnTsDcBvYOmEC/O9/6vvjgV+As/v35wh2lg7FanvrjAoyz0567ISwjNX/8gwePJjw8HAbNEU4gi4hIeuFBDshnJ49euwAHnroIfW4YdOmhmLGn31GRmbR8j7Ag7NmEZW5cMLI0h477fw6o4L02GmD3e3bt2VESIhcWD3H7tNPP+Xhhx/m77//plGjRnh4eJi8P3LkSGsfIexIn5iY9UKCnRBOz17BbtiwYSQmJlKyZMms1f7Dh5NUuTJJffsSDFS6c4dnPvuMlcDuzM9ZGuzM9dgVNthlZGRw9+5dmy5MEcJVWB3sVq5cyaZNm/D29iY8PNxkIq9Op5Ng5+TckpKyXkiwE8Lp2WMoFgw7XowZMybHed8+fehVtixfRkVRGwhMSiIceAZYRcGHYr29vXFzcyM+Pr5AQ7HaOXZgGI6VYCdETlYPxb7xxhu8/fbbREdHExERwYULF9Rf58+ft0UbhZ0oioJbcnLWCQl2Qji9nj17AobtxrLXzLMX37Aw2gBbM197A98BbwJR167lOywaGxur9s41btyYWrVqARAREUFqamq+z09JSSFJ+0MosoBCiNxYHexSUlJ49NFH7TZRWNhPUlISJlFOgp0QTq9Pnz4cP36cs2fPUrp06SJ5Zr169bgD9AK+1JyfBixXFG7ks4/soUOH1OOmTZtSs2ZNANLT07l48WK+z8/eWweygEKI3FidxoYOHcqqVats0RZRxLTbiQES7IQoJurVq2d2Fwp7MW5Plgq8ALwG6k4VTwB+Dz4Imr1cs9POr2vSpInaYweWzbPTzq8zkh47Icyzeo5deno6H3zwAb///juNGzfOsXhi9uzZ1j5C2ElCQoIEOyFEvurVq2fyehYQ6evLosx/Q/wOHYLWreHXX6FBgxyf36NZSdu8eXOT3SbOnDnD/fffj6IoHD16lNDQUIKCgkw+b67HToKdEOZZHeyOHDmibjB99OhRqxskio4EOyGEJYw9dlox3brR8ddf+QWoCBARAe3awfffQ7a9aXfvNqyh9fLyIiwsjETNanxjj91XX33Fs88+S7ly5Th37hy+vr5ZzzLTYydDsUKYZ3Ww27p1a/4XCackQ7FCCEuUKlWKkJAQkzDVqVMnXv/1V1oB+8qXp9zVqxATA337woIF8OKLgKHmnHH1a9OmTfH09DQ7FPv9998DhlW2Bw8epF27duo1MhQrhOWsDnbTpk3L9T2dTsebb75p7SOEnUiPnRDCUvXr1+evv/4CwN3dnbZt2wIQCXzYvz8fRUXBunWQng4vvQQXL8K775oMw7Zu3RqAcuXKUaJECeLj4zl27BiKorBv3z71uitXrpg8W3rshLCc1cFu7dq1Jq9TU1O5cOEC7u7u1KhRw2WD3fz585k/f36xrn4uwU4IYal69eqpwa569eqEhoaq7128eRPWrIHXXgPjvOrp0+HiRfZVr65e16pVK8DwQ3/z5s3Ztm0bFy9eZN++fVzXLL7IHuxkjp0QlrM62JmrJh4TE8MzzzzDwIEDrb290xoxYgQjRowgJiaGwMBARzenUCTYCSEspZ1nV6tWLZP9wa9evQp6PXz0EVSpAqNHg6LAihX0Cw5mFnCXrB47gPbt27Nt2zbAsIORlvTYCVF4dik+FxAQwNtvv+2yvXWuQubYCSEs1ahRI/W4Xr16eHt7U7JkSSDb7hMjR8KPP4KPDwCNb99mO9CkZEmqa3rvtMWVV65cafKsyMhIk9cyx04Iy9mtqnB0dDTR0dH2ur2wAemxE0JYqlOnTjz99NO0bduWESNGAFC+fHnA0GOnKErWxQ8+CFu3kh4cDEADYEtCAjrNCI9xjh4YCt1r5dVjZyyVEh8fb7K6VghhYPVQ7Lx580xeK4rC1atX+frrr+ndu7e1txd2pA12GW5u6DW1pYQQQkuv1/PVV1+ZnCtfvjzHjx8nMTEx57SU1q1Z/+ab1BkzhtpAcHIydOpkKIfSpw/BwcHUq1ePEydO5HhWXnPsqlatyunTpwFDr512rp8QwgbBbs6cOSav9Xo9ISEhDB06lIkTJ1p7e2FH2qHYdG9v+3XfCiFcUrly5dTja9eu5ZhvvCUigmeAn4AOAPHx0L+/Wg6lffv2uQY7RVHQ6XSAaY9d9erVJdgJkQerg92FCxds0Q7hANoeu3RvbzzyvFoIIUwZh2LBMBxbp04dk/f37t3LbaAHED1gAF4//ZRVDiUigvZt2/LFF1/kuG9iYiJ3795V5/Bpg12NGjXU4+t5bGMmxL3Kqk6a1NRUunfvrhafFMWLNtgpmirvQghhCW2wu3DhAsuXLzepW2f83hBcvjxeP/4I48dnfXjGDAatW4d2AkgJzTxf7XCsNthptzeTjgUhcrIq2Hl4eHD48GFbtUUUMZPFExLshBAFpB2KfeWVVxg6dChdunQhKiqK2NhYtUetZs2ahnIoH34In3xiOAYCfvmFPzw8CMi8x3333afeTxvsjHPsfH19aaDZi/bkyZP2+tKEKLasnlb15JNP8uWXX9qiLaKIJcXGZv207OfnyKYIIYohbY9dQkICYBhG3bNnD+fOnVPf0w6f8sorsHatWg6lY2oq4UDL0FC6deumXmauxy4gIIC6deuq5yXYCZGT1XPs0tLSWLJkCX/88QfNmzc36UoHmG2sQi6cTrpmeEMnwU4IUUDaYKd1/PhxkxImNWvWNL2gf38ID4c+feDWLZoCO4C/3LO+JZkLdv7+/pQtW5bAwECio6M5deqUrb4UIVyG1cHu6NGjNGvWDEBdqWRkXNEknJMSF6ce6/39HdgSIURxVLlyZYKCgrh79676X4ATJ06Y/Ptv0mNn1KoV7NgBvXrBxYt4XL5M5zfeoDmwj6xgpyiKSY+dTqejbt267N69m4sXL5KQkICvTCURQmV1sNu6dast2iEcQBvs3AIC8rhSCCFy8vX15ZdffmHbtm0MGTKEGjVqoCgKx48fVwsJg5keO6M6dWDnTrj/fjhyBI87dwgHBpIV7BITE8nIyAAMwQ5Qg52iKJw5c4awsDA7fpVCFC9Wz7G7dOmSacXxbO8JJxYfrx66SY+dEKIQOnTowKRJk6hWrRpVq1YFDD12Z8+eVa8x22NnVKECbNtmKF4M+AG/AY2OHgVMV8Rqg52RzLMTwpTVwa5atWpmN2O+desW1apVs/b2wo50mZOdAdlOTAhhtfr16wMQFxfHP//8A0DJkiXVenS5CgqC33+HgQMB8ASmR0TA3Lkmwc4/8wdQCXZC5M7qYKetDq4VFxeHt7e3tbcXdqTX7rMowU4IYSVjsAPUfVxzHYbNztsbVq9mbUhI1rkxY/B/7z31pfTYCZG/Qs+xGzt2LGBYIPHmm2+aTF5NT09n9+7dNGnSxOoGCvuRYCeEsCVt8WCjPIdhs3Nz4+v27Tm0bh1vZZ4qv3w5S4FhZAW7GjVq4O7uTlpamgQ7IbIpdLA7cOAAYOixO3LkiMlEWU9PT8LCwhivrTIunI57cnLWCwl2QggraXvsjCzusctUsVIl3gauAZ/pdOgUhWeAEOCYlxdgKI5v3DP21KlTZGRkoNfLbtdCgBXBzrga9tlnn2XevHnq3Id7xfz585k/fz7p6emObkqhZGRk4JGamnVCgp0QwkpW99gBTZs2BeBzwL1sWebeuIF7ejp9gaZffw1jxkBgIHXr1uX06dMkJiZy+fJlqlSpYoOvQIjiz+ofcWrVqsXq1atznF+yZAkzZ8609vZOa8SIERw/fpx///3X0U0plMTEREyinAQ7IYSVAgICqFixosm5gvbYPfXUU7Rq1QqA+deu0T09nejM9yqcPw9du8L169SpU0f9jOxXLkQWq4PdokWLTCayGjVo0ICFCxdae3thJyb7xIIEOyGETWQfji1oj52HhwcrV65UR4G2AV2B225uhgsOHICOHWkYGKh+5sKFC9Y0WQiXYnWwu3btmtltZUJCQrh69aq1txd2IsFOCGEP2uFYX19fypUrV+B7VK9encWLF+Pm5oZOp6P7+PGU2LcPKlUyXHD6NA/Pm0etzOsjIiI4ffo048ePZ8+ePQV+3q1bt+jWrRsDBgww2QpNiOLI6p0nQkND2bFjR46adTt27KBChQrW3l7YiQQ7IYQ9aHvsatSoUeitJR999FGaNm2Kt7c3lStXNpzcvh169ICzZ/G5fp2/gV4YeuxefPFF/vrrL7766isuXbqEj4+Pxc8aN26cOm98wYIFjB49ulBtFsIZWN1jN2zYMEaPHs3SpUu5ePEiFy9eZMmSJYwZM4Zhw4bZoo3CDuLj4yXYCSFsLnuws0bt2rWzQh1AlSrw99/QuDEAZYFwwO/QIf766y8Abt68yapVqwr0nF9++UU93rt3r1VtFsLRrO6xe+2117h16xb/+9//1C5sb29v/u///o+JEyda3UBhH9JjJ4Swh5YtW1K1alUiIiIYNGiQ7R9QrhyEh0PfvrBrF0HA3OPHuQhsyrxkwYIFPPPMMxbfUjv8qi3dJURxZHWPnU6nY+bMmdy4cYN//vmHQ4cOcfv2baZMmWKL9gk7kWAnhLAHb29vjh49ypkzZ3jqqafs85CSJWHzZvZkLqDwBX4BHsp8+99//y1QxYJUTeknCXaiuLNZRUc/Pz9atmxJw4YN8cosIimclwzFCiHspUSJEgUuc1KIh/Bpr16syXzpCawCns18/fnnn1t8K+mxE67EJsHu77//5sknn6Rdu3ZcuXIFgK+//prt27fb4vbCDrQ9dopOZ9inUQghipHQmjV5FFia+doNWAK8DBw8eNDi+yiKYvO2CeEoVge7NWvW0KtXL3x8fNi/fz/JmdtURUdHM336dKsbKOxDG+zSPD2hkCvXhBDCUapVq0Y68DwwV3N+AfDQxYuFumd8fLz1DRPCgawOdu+++y4LFy5k8eLFeHh4qOfbt2/P/v37rb29sBOTYCe9dUKIYqhq1aoAKMAYYIGmaPHEmzehEJ0LcXFxtmmcEA5idbA7deoUnTp1ynE+MDCQu3fvWnt7YSfaOXYZBaj3JIQQziJ7/dT17dvzaZkyWSfeeAOmTIFchlojIyNJSkoyOSfBThR3Vge7cuXKcfbs2Rznt2/fTvXq1a29vbATkzl2EuyEEMVQaGioSQHk+g0a8F2tWrymveidd+D//i9HuFu8eDEVK1bMsSWmBDtR3NmkQPGoUaPYvXs3Op2OyMhIVqxYwfjx43n55Zdt0UZhB4lxcRjjnCIrYoUQxZCnpyeVjNuMYdijPCAggFnAq9oLP/wQRo2CjAz11IsvvgjAxWxz8STYieLO6gLFEyZMICMjg+7du5OQkECnTp3w8vJi/PjxvPrqq/nfQDhEanR01gsJdkKIYqpq1apcvnwZMOx6ERAQAMCnwNT33qP05MmG3rpPPoGkJFi4kKQ89oOVYCeKO5sUKH7jjTe4ffs2R48e5Z9//uHGjRu88847tmifsJOM2Fj1WO/n58CWCCFE4XXt2hUwTAtq2LChGuwArvTtC8uWgT7zW93ixfDss+zbvTvX+0mwE8Wd1T12Rp6entSrVw+g0Js+i6KjaP7x0vv7O7AlQghReJMmTaJZs2Y0btwYHx8fAjUrY2NiYuDpp8HLC554AtLTYflySh06hDuQZuZ+EuxEcWeTAsVffvklDRs2xNvbG29vbxo2bMgXX3xhi1sLO9EGOzfNT7hCCFGceHl5MWDAAHWFrLbHLiYmhgsXLnClQwdYvRoyS3LVPXSI1Rh2q8guLi4u34LFd+/epV27dnTo0IFYzeiHEM7A6mA3ZcoURo0aRb9+/Vi9ejWrV6+mX79+jBkzRvaLdWImwU7zE64QQhRn2mD3+++/U7NmTUJDQ3n255+59eWXKJlbXj4I/EDOcJeRkZGjBEp2a9asYdeuXezYsYOffvrJpu0XwlpWD8V+9tlnLF68mCFDhqjn+vfvT+PGjXn11VeZNm2atY8QdqBLSFCP3WQoVgjhIrTBbs2aNWRkroRdtmwZ4eHh7Fi8mKCnn8YX6AesAR4CtMsp4uLi8MmjDNTNmzfV41u3btm0/UJYy+oeu9TUVFq0aJHjfPPmzUlLMzeDQTgDfWJi1gtZFSuEcBHaYBcZGWnyXkREBJ+fPUsfwLhx2APAj4CX5rr85tnFxMSoxzIUK5yN1cHuqaee4rPPPstxftGiRTzxxBPW3l7YiZt2qEGCnRDCRQTmM7Xkq6++4i8wCXd9MQ13+YU1bbDTHgvhDGyyKvbLL79k06ZNtGnTBoDdu3dz6dIlnn76acaOHateN3v2bFs8TtiABDshhCsKyGcxmLEg8Tbgo27dGPfnn5TAEPTWAgORHjtRvFkd7I4ePUqzZs0AOHfuHAClS5emdOnSHD16VL3O1UqgzJ8/n/nz55Oenu7ophSKh7ZApwQ7IYSLyC/YGYWEhFDywQfp/eefrAf8gN4Ywl3C7dt5flaCnXBmVge7rVu32qIdxc6IESMYMWIEMTEx+Xb9O5uMjAw8tfMfJdgJIVyEuWBXsWJFEhMTua0JbE2bNqVq1ar8jSHQbSAr3EVNngzdu0MuCygk2AlnZvUcu8TERBI0KywvXrzI3Llz2bRpk7W3FnaSkJCAr/aEBDshhIswF+zKlClDWFiYyblmzZpRuXJlALYD9wPGiFb20CEYMAC0i8w0JNgJZ2Z1sBswYADLly8HDEUbW7VqxUcffcSAAQPMLqoQjpeQkIBJlJNgJ4RwESVKlECvN/3WFhISkiPYNW3alCpVqqivd2Aa7ti8OddwJ4snhDOzOtjt37+fjh07AvDDDz9Qrlw5Ll68yPLly5k3b57VDRS2J8FOCOGqdDpdjl673IJdUFAQn332Gffddx+TJk1iJ4Zwl+KZWbZ482bo3x80o1IgPXbCuVkd7BISEvDPLHC7adMmBg0ahF6vp02bNurqI+Fc4uPjJdgJIVyWuWDXpEkT9bW/vz81atQAYPjw4WzatInOnTsDsBP46vHHwVi4/Y8/oF8/k3AnwU44M6uDXc2aNVm3bh2XL1/m999/p2fPngBcv37d4tVJomhJj50QwpWZC3b16tVTOyHatm2bY7jWz89PPT5RsiSrnnsONb79+aca7tLS0kzmlRcm2F26dEndEUMIW7PJXrHjx4+natWqtG7dmrZt2wKG3rumTZta3UBhexLshBCuLHuwK1OmDF5eXqxcuZLnnnuOuXPn5viMNtjFxcUx8eef6QlEG0/++Sc88ACxUVEmn4uLiytQSJs5cyZVqlShd+/eFn9GiIKwutzJ4MGD6dChA1evXjWZw9C9e3cGDhxo7e2FHSQkJBCkPeHrm8uVQghR/GQvQRUSEgJA37596du3r9nPaIPd6dOnuXDhAheAnsB2X188EhJg61a8H3sMb0BT4p24uDizI1RpaWkcOnSIsLAw3N0N324nTJgAGDo/oqOji125LOH8rO6xAyhXrhxNmzY16dpu1aoVdevWtcXthY1p59ileniA3iZ/DIQQwimYG4rNjzbY/fXXX+rxHmDpY49B5j19tm/PsbdsbsOxL774Ii1atOCxxx4z+372vWyFsAX5jn4P0g7Fpnl55XmtEEIUN9YGu+zWRkbCxo2QeU1v4Acgc+1srsFu6dKlAKxZs4YjR46wefNmk/evXLmSb7uEKCgJdvcgbbBLl2AnhHAxhQl2Pj4+ORZUGP37778obdrAhg2keXsD8ADwPeCB+Vp2t7NtS9a0aVN1caGRBDthDxLs7kHaodj0XLbMEUKI4kob7Dw8PCyax6bT6XLttbt16xYXLlyADh3Y9vrrGNfEDgC+A+Lu3MnxGePe6Ubm9hWXYCfswapgl5qaSvfu3Tlz5oyt2iOKQIIm2CkS7IQQLkYb5EJCQtDpdBZ9LnuwMy54AEOvHcC5SpV4ADDuRzEIqDNtGmTuv52amsoXX3zBp59+mu/zJNgJe7Aq2Hl4eHD48GFbtUUUkZSYGPV/vCKlToQQLkbbY2fJMKyRdsGfXq9n/Pjx6uu9e/cChmHXrUB/slbGVty5E556CtLS+Prrrxk2bJi61WZeJNgJe7C63MmTTz7Jl19+yfvvv2+L9ogikK6dDyLBTgjhYgob7L799ltWr15NSEgIbdq0wcPDQ/3eZuzEiI42VLb7A3gQ+InMFbLffQfu7ry9davFz5NVscIerA52aWlpLFmyhD/++IPmzZtTIltQmD17trWPEDaWFq2W3ESfx0owIYQojgob7MqWLcsrr7yivlYUhdKlS3Pz5k0OHToEmC6U+B14CPjJzQ239HT45hs+KlGCRwDFgudJj52wB6uD3dGjR2nWrBlgKOqoZem8BlG0lLg49Vhv3A9RCCFcRMWKFdXjqlWrFvo+Op2Oxo0b8+effxIVFUVUVFSOFbC/ASsHDuTJdesgLY3B8fEsAl4k/3B37do10tLSTObyCWEtq/80bS1At7NwDibBTvbzFUK4mPr16zNlyhSOHj3KiBEjrLpXWFgYf/75JwBHjhwxW9rko7Nn8Ro8mEGrVuGmKLwApAEv53PvjIwMoqKiTIKoENaySbmTv//+myeffJJ27dqpXctff/0127dvt8XthY1pg5279NgJIVzQ22+/zZo1a6wOTY0bN1aPDx06ZDbYHTx4kEe++47HFIW0zHPDgU8suL8MxwpbszrYrVmzhl69euHj48P+/ftJTk4GDBNMp0+fbnUDhe3pEhLUY/egIMc1RAghnJw22B0+fNhssDP6AXgKMFasewWYk8/9b926ZWULhTBldbB79913WbhwIYsXL8bDw0M93759e/bv32/t7YUd6BMT1WM36bETQohc1a9fHzc3N8DQY2dcyWo8l913wFBAyZxjPhr4MI/7z5s3jyeffJI1a9aYLWIsREFZHexOnTpFp06dcpwPDAzk7t271t5e2IFbUlLWCyl3IoQQufL29qZOnTqAIdhdvnwZgI4dO5pcN3PmTPV4BXD3o4/UcDcemJH5XqNGjXjuuefUazdu3MiKFSsYPHgwL774ot2+DnHvsDrYlStXjrNnz+Y4v337dqpXr27t7YUdSLATQgjLde3aNce5Pn36cP/99wMwatQoRo8ebfJ+0OjR7NMEtQnAO8C+vXt54IEHzD5n1apVKIolhVKEyJ3VwW7YsGGMGjWK3bt3o9PpiIyMZMWKFYwfP56XX85vTZBwBPeUlKwXEuyEECJPTz75ZI5zvXr14ueff+b48ePMnj0bT09PVq5cSf369fn888/R6XRE9evHS5rPTAY8pk8nKJe5zfHx8Vy6dMkuX4O4d1hd7mTChAlkZGTQvXt3EhIS6NSpE15eXowfP55XX33VFm0UNuYpwU4IISzWunVrQkND1WFYMAyp6nQ66tWrp5577LHHeOyxx9TXvr6+LMLwjXa+8eTbb1Nn+PBcn3X8+HGqVKli2y9A3FOs7rHT6XS88cYb3L59m6NHj/LPP/9w48YN3nnnHVu0T9hYeno6XtoJuhLshBAiTzqdjgEDBqivu3btalEBfuNOTAuAkZrzFRYuZGIunzl27FjhGyoENqpjB+Dp6Un9+vVp1aoVfrJNldNKSEjAJMpJsBNCiHxNnDiRUqVK4ebmxrvvvmvRZ3x9fdXjT4D3SpdWX08HXjfzGQl2wlpWB7unn36aJUuWcO7cOVu0R9iZBDshhCi4ChUqEBERwdWrV2nXrp1Fn9EGO4DvK1SAD7OKn8wExmX7zPHjx61sqbjXWR3sPD09ef/996lVqxahoaE8+eSTfPHFF5w5c8YW7RM2Fh8fL8FOCCEKwc/Pj5CQEIuvzx7sSpQoAePHw/vvq+dmYah1Z3T8+HH27NnDnj17rGusuGdZHey++OILTp8+zeXLl/nggw/w8/Pjo48+om7dulSqVMkWbRQ2JD12QghRNMwGO4D/+z9mavbpngPMKFcOgLi4OFq3bk3r1q1Zu3atek1iYiJr167l6tWrdm+3KN5sNseuZMmSlCpVipIlSxIUFIS7u3uBfrIRRUOCnRBCFA0fHx+T19qg93VoKG9q3ptw7Rojsn3+lVdeITY2FoDx48czaNAg7rvvPql1J/JkdbCbNGkS7dq1o1SpUkyYMIGkpCQmTJjAtWvXOHDggC3aKGxIG+zS9HrQbAMnhBDCdjyy/ftaQvODdFBQEO8Cb2ne/xTQVn+NjIzkyy+/BGDBggWAYXGFdn/ZO3fu0KpVK1q3bk10dLRtvwBRLFldx+79998nJCSEqVOnMmjQIGrXrm2Ldjm9+fPnM3/+/GK3t198fDwVMo9TPT2t/wMghBDCItoeu8DAQADeBtxA7b1bAFSpVo0JFy4AEBERkeM+V65coXTmCtuJEyfy77//AjB16lTmzp1rn8aLYsPqHrsDBw7wxhtvsGfPHtq3b0/FihV5/PHHWbRoEadPn7ZFG53SiBEjOH78uPoXqrgw6bHz8nJoW4QQ4l6SvcfOaArAxKzKdv934QIvZB4nJyfnuM9///2nHu/evVs93r9/v62aKooxqztswsLCCAsLY+RIQ/nFQ4cOMWfOHEaMGEFGRkax69Fyddpgly7BTgghioy2xy5As3gCgPfeg/R0+OADABYDOvIPdhkZGeqxm5ubTdsriierg52iKBw4cIDw8HDCw8PZvn07MTExNG7cmM6dO9uijcKGtOVO0rNN7BVCCGE/2h47vT5rwEyn04FOZyiDkp4OH30EwCJg6eHDOe5z5coV9Vgb7LT3FPcuq4NdcHAwcXFxhIWF0blzZ4YNG0bHjh1z3eRYOFZSTAzG6byKBDshhCgy2lWyZgOZTgcffkhsXBz+n38OwLP79sGcOSb30fbYaUfFpMdOgA2C3TfffEPHjh1zdisLp5SmWTWlZKuxJIQQwn60PWrakiUmPW06HSnvvsu7n3/OZOO5sWN5Hfgg8+XSpUtp1aoVw4cPlx47kYPVwa5v377cvXuXjz76iBMnTgBQv359nn/+eXXVj3Ae2mAnNeyEEKLo6HQ6s+ezl0Xx8vbmTSAFmJZ5bibgCRh3qX355Zfp27evBDuRg9V/Cvbu3UuNGjWYM2cOt2/f5vbt28yZM4caNWrICh0nlB4Tox7r/Pwc2BIhhLh3jRs3Tg16K1asMHnP29sbgHeAhVWqqOffISvogaEUigQ7kZ3VfwrGjBlD//79iYiI4Mcff+THH3/kwoULPPDAA4wePdoGTRS2lJFZxRxA7+/vwJYIIcS9RdtjV6tWLXbu3Mmvv/5K//79Ta5zd3dXQ9qXISHEvv22+t6bGHrvAKKjo/OcY5eamsqkSZOYMmWKVKi4h1g9FLt3714WL16Mu3vWrdzd3Xn99ddp0aKFtbcXNqbExanHbjIvUgghioz2+yRAmzZtcr3Wy8uLxMREkpOTuT10KBOnTuXTzPdexzAse+f2bZMeu+xbjS1cuJAZM2YAULVqVZ577jlbfBnCyVndYxcQEMClS5dynL98+TL+0iPkfOLj1UMJdkIIYV/GYVY/Pz+eeOIJiz9nHI5NSkoiKSmJ+cCLmvdHA9VmzeJaZKR6LjEx0eQeixcvVo+XL19e0KaLYsrqYPfoo4/y/PPPs2rVKi5fvszly5f57rvveOGFFxgyZIgt2ihsSJeQoB57SEkaIYSwqyFDhrBjxw6OHz9eoDJgXpkF5JOTk9XAthgYGxSEsY+uw5EjzE9LwzjAmz3YaeW2cEO4HquHYmfNmoVOp+Ppp58mLS0NMKzwefnll3n//fetbqCwLZ3mL74EOyGEsC+dTke7du0K/Dljj11ycjJJSUnq+b9r1uSpvXtZjmGP2RcAD+A5UK9TFIXvvvuOI0eOmLRD3BusDnaenp58/PHHzJgxg3PnzgFQo0YNk61ThPNw0wQ7WTwhhBDOydhjl5SUZNITV758eb4FUoFvMXwTH4oh3L2fOSLz119/8fjjj+d671u3bhETE8OqVasYOHAgderUsdeXIRzA6mBn5OvrS8OGDQH5ycCZuWl+8pM6dkII4ZyMwS46Opo///xTPV++fHkAVmMId6swLKR4HAi6eBFSU5mTbacKrS+++IJhw4apr7/99lsOm9m2TBRfNil68+WXX9KwYUO8vb3x9vamYcOGfPHFF7a4tbAx95SUrBcS7IQQwikZh2IB3n33XfXYGOwA1gGDgOTM130SEuDRRwk0s12kTqcjISHBJNQBHDlyhNTUVBu2XDia1cFuypQpjBo1in79+rF69WpWr15Nv379GDNmDFOmTLFFG4UNeUqwE0IIp2fsscsuODgYT09P9fVvwABAHYtZu5Zxu3bhme1zSUlJLF261Ow9o6KirG2ucCJWD8V+9tlnLF682GQFbP/+/WncuDGvvvoq06ZNy+PToqh5aH8yk2AnhBBOSdtjp+Xj44O/vz+3bt1Sz/0OPAD8AvgAYZcu8SPwEFm9edHR0fz+++9m73nt2jUqVapku8YLh7K6xy41NdVsIeLmzZurq2SFc0hNTcVHU8xSgp0QQjin3HrsvL29zX5v3QL0BZTMhYt9gbWAMR7GxMTkus3n1atXrW6vcB5WB7unnnqKzz77LMf5RYsWFagYo7C/xMRETKKcBDshhHBKuQU7Hx8fk/InWluB5B9/JNnDA4DewE8YevEuX77MlStXzH7O2mAXHR1Ninaaj3Aomy6eeOGFF3jhhRdo1KgRixcvRq/XM3bsWPWXcKyEhAQJdkIIUQzkNhTr7e1NcnKy2fcAElq25NO+fTHuCt4TwxCttgBZ//79qVixovrammC3Zs0aKlWqROXKlbl7926h7yNsx+o5dkePHqVZs2YAah270qVLU7p0aY4ePapeJyVQHC8+Pl6CnRBCFAN59djlJTExkaNBQfTEMPcuAOgOrMcwPBuPYceoN998k5YtWwKGOXaFcfv2bQYPHgxAXFwcW7duZeDAgYW6l7Adq4Pd1q1bbdEOUQS0PXYZOh36XP7hEEII4Vh59djlJTExkcTERP4BegCbgCCgM7ARw/Bs27ZtTVbWFrbH7tChQzmeLRzPJkOxonjQBrsUd3eQXlQhhHBKefXYzZgxI9fPJSYmqnPw/sXQW3c7870OwPYSJagWHEzZsmXVkbTCBrvsc/0k2DkHCXb3EO1QbIpn9ipHQgghnEVeq2JHjhzJBx98QL169XK8n30Lsv0Ywp2xOEpYfDz07Il7XBzBwcEAXL9+vVBtzD7XT4Kdc5Bgdw/R9tilyTCsEEI4rbzq2Pn6+vLaa6+xb98+hg8fbvK+cShW6yDQDbhhPLFnD/ToQZWAAMCwqjU3t27donfv3jz88MMkZO5FayQ9ds5Jgt09RBvs0iXYCSGE08qrx87Ix8eHzz77jDfffFM9px2K1ToMLHv6aQgJMZzYt4+vo6IIBO7evUuGtsapxurVq9m4cSM//PAD48aNM3kv+3NyK8MiipZVwS41NZXu3btz5swZW7VH2FFiXJxarDI9n5VVQgghHCevHru8zmXvsVu/fj1jx45l6tSpjP7iCwgPhzJlAKifkMAmIEBROH36NIMGDeLBBx8kPj5e/fzu3bvV44ULF5oUR7ZkKDYiIoJvv/3W5J7CvqxaFevh4cHhw4dt1RZhZ8m3b6vHigQ7IYRwWgUpd6I9p51jV6pUKXr37k3v3r2zLq5fH7ZuhS5d4MYNWmEoi9K6Xj1iMi+ZNm0aM2fOVO+hdezYMcLCwtRnaWUPdunp6bRs2ZKbN28yfPhws5sZCNuzeij2ySef5Msvv7RFW4SdpWnmURi3nRFCCOF8cgt2nmYWvuXWY5drzbvMcBeT2SvYmqyadwAbN25UL42NjTX56OXLl9Xj7D122YNeVFQUN2/eBAy9faJoWF3HLi0tjSVLlvDHH3/QvHlzSmQrejt79mxrHyFsRBvspDixEEI4r9yGYs0V+/f391ePb926pQasPGveNWjAF0OG8OTSpZQB2mCoc9cLTLYHi4mJMfmYNtjl12OX27w9YV823Xni9OnTJu/JbhPOJV37F9TPz3ENEUIIkafceuzMqV27tnp87Nix/HvsMiXXqkV34E8gBGiLoefuZU1PXF7BLr85dhLsHEN2nriHZGi61N00P+EJIYRwLuZ627p372722vr166PT6VAUhWXLlqnn8wt2JUuW5Cio4a40hnD3xZUrEBsL/v45hmL/++8/9Ti/VbHahRai6Ei5k3uIEhenHusl2AkhhNPK3mO3YsUKvvvuO7PX+vr6UrNmzRzn89t+LCgoCIAjGOrc3cw83yIlBe6/H2Jjreqx0w7piqJjdY+d0fHjx7l06VKO/5H9+/e31SOEtTTLzd0DAx3YECGEEHnJHuwef/zxPK9v2LBhjtJj+fXYGYMdGMKdseeuFMDOndC7N+l375p8piBz7CTYOYbVwe78+fMMHDiQI0eOqF3BkDW/Lj093dpHCFvRBDsPzV9oIYQQziW/3rbsypYtm+OcJUOxWofJCnfBADt2sADoDRi/e1y5ckW9Pr9VsRLsHMPqodhRo0ZRrVo1rl+/jq+vL8eOHWPbtm20aNGC8PBwGzRR2Ipe89OUBDshhHBeBVk8AfDggw/mOFeQHjujQ8D97u4ome91BNaDumtRUlKS2oEjPXbOyepgt2vXLqZNm0bp0qXR6/Xo9Xo6dOjAjBkzGDlypC3aKGzETfOXUIKdEEI4L72+YN+ee/bsySeffGJyztI5dtntTU/nyJw5GEvadwJ+IyvcpaamsmvXLn744QeTz0mwcw5WB7v09HS1hk7p0qWJjIwEoEqVKpw6dcra2wsb0muCnSyeEEII55Wamlqg63U6Ha+88orJufx67IKDg82eVxSF035+9ADuZJ7rDPwK+AIDBw6kXbt2OT4nwc45WB3sGjZsyKFDhwBo3bo1H3zwATt27GDatGlUr17d6gYK2/HQzoeQAsVCCOG06tSpox4/99xzFn/OV7OrkLm9W7U8PDz4888/zb536tQpDoBJuOuCIdyFr19v9jMyx845WB3sJk+erBYhnDZtGhcuXKBjx46sX7+eefPmWd1AYTse2r9kEuyEEMJp+fr68u+//zJv3rwC7eCk3dv11q1b+V7ftWtXunXrluO8scNmP3AfcNd4PfALYK4vUHrsnIPVq2J79eqlHtesWZOTJ09y+/ZtSpYsKTtPOBlPbde+BDshhHBqLVq0oEWLFgX6TKlSpdSSJJYEOzDt5TMyBjuAfRjC3WYgCEPNu1+AfoA2yhkXVhi/90uwcwy7FCgODg6WUOeEPLVVwCXYCSGEy6latap67OnpadFnzM3Fy75F6F6gJ2Dccbw78DM5e+60JVAk2DlGoXrsxo4da/G1BelCFvaTmpqKb+YSdUCCnRBCuKCPPvqI3377jbS0ND788EOLPpPfIov169ezYsUKVqxYQU9gExCIYf7dT0B/wDi77pNPPuHVV1/F29s712B37do1+vXrR1BQEL/++muBS7uIvBUq2B04cMDW7RB2lpCQgEmUk2AnhBAup3r16pw/f57k5GRq1Khh0WfyCnZDhgyhd+/erF69GoA9QC8M4S4AwxDtT8AADOHu9ddfB+C1117LNdiNHTuWvXv3AjBnzhwmTJhgUTuFZQoV7LZu3Wrrdgg7i4+PNw12+fyEJoQQoniqVKlSga43N8fOqH379oDpsO5uDOHudwzhriewDniQrHBnLthlZGRw9epVVq5cqZ4zBjxhO1Yvnpg2bVqu7+l0Ot58801rHyFsQNtjl+TmhncBi18KIYRwTR4eHrm+Z6xXl/2af4D7MYQ7fwxBby2GcGecZZc92KWlpfHII4+YnJN5eLZndbBbu3atyevU1FQuXLiAu7s7NWrUkGDnJBISEjCWJE7x8KBguxAKIYRwVbktdvTy8qJRo0aA+fC3i6yeO38MQW8tMBBDiDMX7Hbu3GlyToKd7Vkd7MzNt4uJieGZZ55h4MCB1t5e2EhCQgLlMo9TLVwpJYQQ4t5VsWJF3N0NMSG3Fba7yOq58wN6Az8Cl8+cyRHaWrdunePzEuxszy7jcQEBAbz99tvSW+dEtHPs0iTYCSGEyJRbj12FChXU4+w9du+99x6fffYZHTt2ZMzq1dwPxGW+1wc41bAhM995x+QzR48ezfGMgm6dJvJnt4lW0dHRREdH53+hE5g1axYNGjSgYcOGfPPNN45ujl0kxMVlBbt8NoYWQggh8gp2AQEBDB8+nG3btjFo0CB2YOitM4a7+zMyWAPk142grXsnbMPqodjs24YpisLVq1f5+uuv6d27t7W3t7sjR47w7bffsm/fPhRFoWvXrjzwwAMEBQU5umk2lXz3rnqcLsFOCCFEPipWrKgeZx+K1b7W6/U899xzLFmyhD7ABqAE8ADwAzAYyG3AVYKd7Vkd7ObMmWPyWq/XExISwtChQ5k4caK1t7e7EydO0LZtW7wzw05YWBgbN27ksccec3DLbCvlzh31WJFSJ0IIIfKRV49d9qLCX3zxBWPGjKFRo0b0AdZjCHf9gNXAw5gPdzExMRa3Z/Pmzaxbt47Ro0dTq1Ytiz93r7F6KPbChQsmv86dO8c///zD9OnT8ff3z/8G+di2bRv9+vWjQoUK6HQ61q1bl+Oa+fPnU7VqVby9vWndujV79uyx+P4NGzYkPDycu3fvcufOHcLDw7ly5YrV7XY2qZoeOyWPmkVCCCHuLfXq1TN7vnz58upxfsFOp9PRoEEDALYBfYH4zPf6A98D5oqq3NF0OuQlLS2Nnj17smDBgmIxGuhITl/MLD4+nrCwMObPn2/2/VWrVjF27FimTp3K/v37CQsLo1evXly/fl29pkmTJjRs2DDHr8jISOrXr8/IkSPp1q0bgwYNok2bNri5uRXVl1dk0rU/FcmuE0IIITI98cQTDBgwgCZNmpicDwgIUI/zGoo10i7C+AvDUGxC5usBmA930dHRFg3HJiUlqcfnzp3L9/p7mdVDsbntG6vT6fD29qZmzZoMGDCA4ODgQt2/d+/eeabz2bNnM2zYMJ599lkAFi5cyG+//caSJUvUbUoOHjyY5zNeeuklXnrpJQBeeOEFl+zi1QY7nZ+fA1sihBDCmbi5uamjYUFBQerCx9DQUPWa/HrsjMaOHavuER+OIdz9CvhiKF78NTAE0OxcTvXq1alfvz4hISF89dVXZmvmZS+LMnXqVMaPH2+TkUFXY5M6dvv37yc9PZ06deoAcPr0adzc3Khbty4LFixg3LhxbN++nfr161vdYK2UlBT27dtnMpdPr9fTo0cPdu3aZfF9rl+/TpkyZTh16hR79uxh4cKFuV6bnJxs8tNFQeYHOFJ6bKx6rJe/CEIIIcz45ZdfeOKJJ+jevbtJD172sJVbXbtp06YRFhZG6dKl6du3L1sxzLP7FfABHgVOAG9rPhMZGUlkZCRgyBQzZszgwQcfNLlv9l69adOmER0dzdy5cwv+Rbo4q4diBwwYQI8ePYiMjGTfvn3s27eP//77j/vuu48hQ4Zw5coVOnXqxJgxY2zRXhM3b94kPT2dsmXLmpwvW7Ys165ds/g+AwYMoH79+jz55JMsXbpULchozowZMwgMDFR/aX+icWZKXJx6LMFOCCGEOR07duTSpUssXbrU5Hz2IJdbj12JEiV4+umn6d69u3ruTwwrYzMyX7+V+dqckydPMnDgwBwjbeYKGX/88ce5fh33Mqt77D788EM2b95sMhYfGBjIW2+9Rc+ePRk1ahRTpkyhZ8+e1j7KbgrSuzdx4kST4eeYmJjiEe40wc49MNCBDRFCCFHcWNpjZ5Q9+K0HXgdmZb5eBpwGDufy+eXLl9OkSRPS09P55ptvTObNi7xZHeyio6O5fv16jmHWGzduqMOUQUFBdtk2pHTp0ri5uREVFWVyPioqinLlyuXyKet4eXnl+pOKU0tIUA8l2AkhhCgIS+fY5eUjoDHwNIZSKD8BLYGbZq719PTkzp07rF+/nmeeeabAz7qX2WQo9rnnnmPt2rX8999//Pfff6xdu5bnn39eHSPfs2cPtWvXtvZROXh6etK8eXO2bNminsvIyGDLli20bdvW5s8rzvSaYOdZsqQDWyKEEKK4sXQoVuuLL77Ice5FwFiQrCqGGnfmephmzpxJcHAwTz75ZAFbKqzusfv8888ZM2YMjz32GGlpaYaburszdOhQtXhx3bp1zf4PtkRcXBxnz55VX1+4cIGDBw8SHBxM5cqVGTt2LEOHDqVFixa0atWKuXPnEh8fr66SFQZ6zVJxCXZCCCEKoqBDsQDPP/88Z86cYebMmeq5ZGAgcNLfH//YWLoAHwMjbNnYe5zVwc7Pz4/FixczZ84czp8/DxiWLvtpSmpkr41TEHv37qVr167qa+P8tqFDh7Js2TIeffRRbty4wZQpU7h27RpNmjRh48aNORZU3OvcNcHOw8W2SxNCCGFfhR2KLWGmbmokcPDNN+kweTK6lBT+BxwCFlnfTIENCxT7+fnRuHFjGjdubBLqrNWlSxcURcnxa9myZeo1r7zyChcvXiQ5OZndu3fTunVrmz3fVbhrlorLqlghhBAFYUmB4oJQWrcmPrPeHcCnQI/C3EdR8r8o09q1a2nSpEmOFb+uxuoeO4AtW7awZcsWrl+/TkZGhsl7S5YsscUjhJXctYtXZOcJIYQQBVDYHrvU1FSz59u0aYOubVs+euUVxmHYkeIHoANwtADtmj59Om+88YZF1w4aNAiA5557zqWna1ndY/f222/Ts2dPtmzZws2bN7lz547JL+EcPLV/uSTYCSGEKIDCzLGDnMGubdu2HD16FE9PTzw8PJji5cW6zPcCgd+A8lhu8uTJJq83bNjA2rVrC9ST52qs7rFbuHAhy5Yt46mnnrJFe4SdeGYubAEk2AkhhCiQwvbYZS91tnPnTpPXJQICeOLGDcIxlD6pDPwCdAbiC9jGnTt30qdPHwDWrVvHgAEDCngH12B1j11KSgrt2rWzRVuEnSiKgrcEOyGEEIWUvYfOzc3Nos9pg51enzNy+Pv7k4Bh27GIzHPNgZUUPKC899576vFrr70GQHp6OqtXr+aPP/4wuTb7tDFXYnWwe+GFF/j2229t0RZhJ6mpqfhqT0iwE0IIUQDZe+wspQ125oZvjbtWRQF9gbuZ5/sBnxXgOXfv3mX9+vXqa+NQ7OrVq3nkkUe47777cm2Xq7F6KDYpKYlFixbxxx9/0Lhx4xz/82drVr0Ix0hISMAY5VL0ejwt/ElLCCGEAPsFO20VjePAIOB3DIspXgRuAJNzfCqncePGmbw29sgNHz7c7PU9evQgLi6On376iSpVqljwhOLD6mB3+PBhtU7d0aOma1l0Op21txc2EB8frwa7ZHd3rFukLoQQ4l5T2PIm2mBnbl6ej4+PyeutGLYcW4FhSPENDOHu43ye89VXX5m8NvbY5ZZDduzYAcDLL79s0tPnCqwOdlu3brVFO4QdaXvskt3dkSp2QgghCqKwPXbly2etca1bt26O97MHO4DvgFIYatsBzAVuAd/k8Zz09HST18Zgd/fu3Tzbt3379jzfL45sUscO4Pjx41y6dMkknet0Ovr162erR4hCSkhIoEzmcaqVRSWFEELcewob7CZPnszq1atJSEgwWxjYXLADmI8h3L2d+XopcAdDOZTsRo8eneOcoih89913+bbP19c332uKG6uD3fnz5xk4cCBHjhxBp9Pl6P7MnqJF0YuPi1N77NIk2AkhhCigwg7FBgUFce7cOTIyMszew9vbO9fPTgNKA69iCCurgfuBbdmu+/jjnAO1iqIwZMiQfNtnbsuz4s7qVbGjRo2iWrVqXL9+HV9fX44dO8a2bdto0aIF4eHhNmiisFZidLSa4NPy+EskhBBCmGNpeRNz3N3dcw2GufXYGY0CjHU3fDD02FlSYC02Ntaitrlij53VwW7Xrl1MmzaN0qVLo9fr0ev1dOjQgRkzZjBy5EhbtFFYKUWzA0i6BDshhBCF0Lt3bwAmTZpks3vmFezuu+8+zpw9y+7hw/k185wfsAFolc99Ld35yhWDndVDsenp6fhnbipfunRpIiMjqVOnDlWqVOHUqVNWN1BYL1UzeVSRYCeEEKIQfv75Z86cOWN2EURh5RXsAgICqFGjBsHlyjEY+AnoBQRgKInSA9hn5fNdMdhZ3WPXsGFDDh06BEDr1q354IMP2LFjB9OmTaN69epWN1BYT9tjp7jgfAIhhBD25+7uTr169WxayiyvOXbG4sWBgYEkAw8CWzLfCwI2AU2sfH54eDjXrl3LcT46Oprk5GQr7+4YVge7yZMnq4UAp02bxoULF+jYsSPr169n3rx5VjdQWC8tOjrrhQv+dCKEEKJ4yqvHrly5coAh2AEkAf2BvzLfDwY2A02tbMPTTz9t8vrQoUNUqFCBypUr51suxRlZHex69erFoEGDAKhZsyYnT57k5s2bXL9+nW7dulndQGG99JgY9VinqfIthBBCOFL2YGfc+kuv1/P8888DpoWNEzBsPWasPlca+BNobUUbNm/ebPL6ySefJCEhgevXrzNjxgwr7uwYVgc7c4KDg2XXCSeSoVkdJMFOCCGEs8ge7BYsWMDUqVM5evQoNWrUAAxbl2rFA33ICndBGHruOtmoTZcvX1aPr1+/bqO7Fh27BDvhXJS4OPXYLXPOghBCCOFo2efY1axZk7feeot69eqp5wYMGKAGwBUrVgAQi2EhxR+Z1/gDG4GeNmiTtrRLcazFK8HuHiDBTgghhDPKr44dQKlSpTh+/Djbtm0zKTqcADwAaikUH+BnDPPwCkob4CTYCaenS0hQj92DghzXECGEEEIjr1WxWlWrVqVjx445pnklA4Mw7EoB4AWsAZ4qYDu026Hq9VnRyLg4dPPmzfTu3Zs1a9ao55yVzfaKFc5LG+w8JdgJIYRwEoXdg1YrFRgCJAJPYwg2y4FywIcW3iMlJUXtPTTXY9ezp2GQd+PGjdSvX599+/ZZHEqLmvTY3QP0iYnqsWfJkg5siRBCCJHFmq3KtNKBZ4BPNOc+AOYAlizlNNasy8jIIDIyMuu+6els2bLF5Nrjx4+zePFi6xpsRxLs7gFumiKLXsHBDmyJEEIIkcXd3XYDhwowEpioOTcaWAmY36k2i3Eo9tlnnzU5n56eTo8ePXJcf/Xq1cI31M4k2N0D3DVLxaXHTgghhLOwZbAzeh9D711a5utHgXBvb/JaOmjssVu+fLnJ+dzm0znzPDsJdvcAD82kUKljJ4QQwlmULl26wJ955JFH8r3mKwyrY+MzX7dNSmIHUCWX60+cOMGkSZNynM9tVawEO+FQnqmpWS9kr1ghhBBOolatWrz88stUqVKF8PBwiz6zYMECgi2YVrQB6AZQqhQADYE9QBsz1/br18/sLhMnT540e29FUSxqqyNIsLsHeKWlZb2QYCeEEMKJLFiwgIiICDp37mzR9aVKleLq1at8/PHH9O+fd9W6PQC7dnEhc/VtGWArhlW0ljh//rzZ89JjJxxKgp0QQghX4unpyciRI3nqKQsq1tWqxZBq1TCubfUGvgXexrIVs+ZIsBMOoygK3pl/ANMAPPNbGySEEEIUf3q9njVr1gAQ6+7O/cDnmvenAKuAwnR3OHOwkwLFLi4lJUX9Q5vk5oafrrA/nwghhBDOJbe5bqmpqdy5c4eQkBDAEPLSgOHASeAjDD1bDwP1gIHA2QI815m3GpMeOxcXHx+vBrtkOywrF0IIIZyNu7u7GurAdJuwuUA/IDrzdUPgX6BvAe6fpCkj5mwk2Lm4hIQENdilSLATQghxD8q+x+x6oCVwLPN1EPArMBXL5t0laLbqdDYS7FycSbCT+XVCCCFciKVlR7Q9dkZnMJQ+Wa059xbwM5BfKf+9e/c6bckTCXYuLiEmBq/M4zQJdkIIIVyIdrg1L+aCHUAc8AjwOob9ZgEeAA4AbfO435kzZ9i3b5/F7SxKEuxcXNKtW+pxmre3A1sihBBC2FaXLl3o168fISEhjB49mvbt2/Pzzz/nuC77UGx2HwL3AzczX1cBtmEIfLl9csOGDYVutz3JpCsXl3LnjnqcLsFOCCGEC9HpdPz8889kZGTk2isHuffYaf0BNMFQ464ThoA0E+gKPA3cyHb98ePHAUhLSyM1NRUfH59CfAW2Jz12Lk4b7DKc5A+dEEIIYUv5Bbf8euyMrmDYhuwdwFip7n7gEIaAp3Xx4kViY2OpVasW5cuXV4Oeo0mwc3HaYKdIsBNCCHEPsjTYgWGu3RTgPuBa5rnywJ8Y6t8Z561fvHiR2bNnExERQXR0NI899pgNW1x4EuxcXHpMTNYL2U5MCCHEPSi/Hr2WLVvmOPcnEAZs0pwbC+zDMGQbGRlpMs/u5MmT1jfUBiTYuThtsNP5+TmwJUIIIYRj5NVjt3v3bnr27Gn2vesYhmLHAMaSxA2APcBEYO/u3eq1zrIbhQQ7F5cRG6seS7ATQghxL8qrx65FixYkJyebfc/X1xcFw24VzTGUQQHwAKZjWDlbI/Ocs+wfK8HOxWmDnVtAgANbIoQQQjhGbsHuoYceQq/X57pF2KJFi9Tj40Br4D2yat61Aw4Do3GeQOUs7RB2osTHq8cS7IQQQtyLchuKXb3asO9Ebj12TzzxBH37Zu0imwpMBjoC5zLP+QJzMPTeERFhk/ZaQ4Kdi9Npgp17YKADWyKEEEI4hrkeu549e6qBL7dgBxBgplNkF4aFFfM052oD+Ppa1U5bkGDn4nSajYo9goIc1xAhhBDCQZ5//nmT13Xr1jUZZjU3FPvhhx8Cuff2xQOjMBQzPgu8ClCmjE3aaw0Jdi5On5ioHnsFBzuwJUIIIYRjPPbYY8yePZsZM2aQmprK8ePHqVKlivr+E088oR63bNmSuLg4xo8fb/ZeU6dONXn9N9AQWGWPhheCbCnm4tw03cueJUs6sCVCCCGEY+j1esaMGZPr+/369WPu3LnExcXx+uuv4+HhYfa6kJAQvM1sz5n7QG7Rk2Dn4tw1wc67VCkHtkQIIYRwTjqdjlGjRll07e3bt+3cGuvIUKyL85AeOyGEEMJm0tLSHN2EPEmwc3EeqanqsRQoFkIIIQom++KJkSNH4ubm5qDW5E+CnYvz0v5kIXvFCiGEEIWmKApVq1bl+PHjLFiwwOS9/v37O6hVpmSOnYszBrsMQO/j49jGCCGEEMWMuXIntWvXJiUlxeTcO++8U1RNypP02Lk4n8xNiZP0eshjE2QhhBBC5KQddnV3z+oPC8pWG9bT07OompQnCXYuTFEUfBQFgMQ8NkAWQgghhHnTpk1Td65YuXKler5UtkoTuZVIKWoyFOvCkpOTMc6qS3aX/9VCCCFEQVWuXJlTp04RGxtL06ZN1fM+2aY3SbATdhcfH68GuxQJdkIIIUSh1KxZM99rZChW2F1CXBzG7YhTnOQPnBBCCOGK3J2kA0WCnQtLvHVLPU6VYCeEEELYVGBgoHqcfWjWUSTYubAkTbBL9/JyYEuEEEII17Nt2zYGDBjA0qVLKeEktWKdo99Q2EXynTvqcbqZTYuFEEIIUXiNGzdm3bp1jm6GCemxc2GpmmCX4SRdxEIIIYSwHwl2Liz17l31WPH1zf1CIYQQQrgECXYuTBvsZJ9YIYQQwvVJsHNh6TExWS8k2AkhhBAuT4KdC8uIjVWP9f7+DmyJEEIIIYqCBDsXJsFOCCGEuLdIsHNl8fHqoZumiKIQQgghXJMEO1eWkKAeukuwE0IIIVyeBDsXptcEOw8JdkIIIYTLk2DnwnSJieqxV3CwA1sihBBCiKIgwc6FuSUlqceeJUs6sCVCCCGEKAoS7FyYe3KyeuxdqpQDWyKEEEKIoiDBzoV5SLATQggh7ikS7FyYR2qqeixDsUIIIYTrk2DnwrzS0tRjnWwpJoQQQrg8CXYuzBjsEgHc3BzaFiGEEELYnwQ7F+adng5Aol7+NwshhBD3AvmO78J8MzIASJLeOiGEEOKeIMHORSmKgm/mcbIEOyGEEOKeIMHORSUmJGBcLpHs4eHQtgghhBCiaEiwc1EJt2+r/3NTJdgJIYQQ9wQJdi4q6dYt9TjV09OBLRFCCCFEUZFg56ISb95Uj9O8vBzYEiGEEEIUFQl2Lirlzh31ON3Hx4EtEUIIIURRkWDnorTBLsPb24EtEUIIIURRkWDnolLv3lWPFV/f3C8UQgghhMuQYOei0qKjs17IPrFCCCHEPUGCnYtKi4nJeiHBTgghhLgnSLBzURmaYKfz93dgS4QQQghRVCTYuaiM2Fj12E2CnRBCCHFPkGDnquLj1UO3gAAHNkQIIYQQRUWCnavSBDv3wEAHNkQIIYQQRUWCnYvSJSaqxx5BQY5riBBCCCGKjAQ7YP78+VStWhVvb29at27Nnj17HN0kq+kTEtRjz5IlHdgSIYQQQhSVez7YrVq1irFjxzJ16lT2799PWFgYvXr14vr1645umlXckpLUY6/gYAe2RAghhBBFxd3RDXC02bNnM2zYMJ599lkAFi5cyG+//caSJUuYMGGCg1sHR7/4grt79xb4cyU1wVSCnRBCCHFvuKeDXUpKCvv27WPixInqOb1eT48ePdi1a5fZzyQnJ5OcnKy+jtEWAraD27Nn0+nECavu4VO6tI1aI4QQQghndk8Pxd68eZP09HTKli1rcr5s2bJcu3bN7GdmzJhBYGCg+is0NLQomlpoFzw88K9Y0dHNEEIIIUQRuKd77Apj4sSJjB07Vn0dExNj13AX8NJL/JVL72F+9F5e1Bk3Dr27/G8WQggh7gX39Hf80qVL4+bmRlRUlMn5qKgoypUrZ/YzXl5eeHl5FUXzAGgyahSMGlVkzxNCCCFE8XVPD8V6enrSvHlztmzZop7LyMhgy5YttG3b1oEtE0IIIYQouHu6xw5g7NixDB06lBYtWtCqVSvmzp1LfHy8ukpWCCGEEKK4uOeD3aOPPsqNGzeYMmUK165do0mTJmzcuDHHggohhBBCCGenUxRFcXQjirOYmBgCAwOJjo4mICDA0c0RQgghhIspSNa4p+fYCSGEEEK4Egl2QgghhBAuQoKdEEIIIYSLkGAnhBBCCOEiJNgJIYQQQrgICXZCCCGEEC5Cgp0QQgghhIuQYCeEEEII4SIk2AkhhBBCuAgJdkIIIYQQLuKe3yvWWsYd2WJiYhzcEiGEEEK4ImPGsGQXWAl2VoqNjQUgNDTUwS0RQgghhCuLjY0lMDAwz2t0iiXxT+QqIyODyMhI/P390el0dnlGTEwMoaGhXL58Od/Nf8W9S/6cFEzLli35999/Hd0MizmyvUXxbFs+wxb3suYeBf1sQa6Xv+cF4yp/zxVFITY2lgoVKqDX5z2LTnrsrKTX66lUqVKRPCsgIED+Iot8yZ8Ty7i5uRWr3ydHtrconm3LZ9jiXtbco6CfLcyz5O+5ZVzp73l+PXVGsnhCCHFPGjFihKObUCCObG9RPNuWz7DFvay5R0E/W9z+LBYnxe331hbtlaHYYiAmJobAwECio6OL1U8eomjJnxMhXJ/8PRf5kR67YsDLy4upU6fi5eXl6KYIJyZ/ToRwffL3XORHeuyEEEIIIVyE9NgJIYQQQrgICXZCCCGEEC5Cgp0QQgghhIuQYCeEEEII4SIk2Dm5+fPnU7VqVby9vWndujV79uxxdJNEMTVw4EBKlizJ4MGDHd0UIYSNXb58mS5dulC/fn0aN27M6tWrHd0k4SAS7JzYqlWrGDt2LFOnTmX//v2EhYXRq1cvrl+/7uimiWJo1KhRLF++3NHNEELYgbu7O3PnzuX48eNs2rSJ0aNHEx8f7+hmCQeQYOfEZs+ezbBhw3j22WepX78+CxcuxNfXlyVLlji6aaIY6tKlC/7+/o5uhhDCDsqXL0+TJk0AKFeuHKVLl+b27duObZRwCAl2TiolJYV9+/bRo0cP9Zxer6dHjx7s2rXLgS0TjrBt2zb69etHhQoV0Ol0rFu3Lsc1MmwvRPFly7/j+/btIz09ndDQUDu3WjgjCXZO6ubNm6Snp1O2bFmT82XLluXatWsOapVwlPj4eMLCwpg/f77Z92XYXojizVZ/x2/fvs3TTz/NokWLiqLZwglJsBOiGOjduzfvvvsuAwcONPu+DNsLUbzZ4u94cnIyDz74IBMmTKBdu3ZF1XThZCTYOanSpUvj5uZGVFSUyfmoqCjKlSvnoFYJZyTD9kK4Nkv+jiuK8v/t3GlIVN8bB/CvM40t465pZpmaUbY5ZY1ZQ9kqhFK9iIgyC6EFzNTU9gV9YWgu2e+FZC8yCFqhglYTjLBSwyYiWlWyoGg3ZyydmvN/1f07P0tnTH+j1+8HfDH3nHvu85zLhWfOuSPWrl2LefPmITY21l6hUi/Awq6XcnR0RFhYGMrKyqRjZrMZZWVliIiIsGNk1NtYu22/YMECLF++HJcvX8aIESNY9BH1EdY84xUVFTh16hTOnz8PjUYDjUaDhw8f2iNcsrMB9g6A/iwlJQVxcXGYNm0atFotCgoKYDQasW7dOnuHRn3QjRs37B0CEfUQnU4Hs9ls7zCoF2Bh14utWLEC79+/x969e/H27VtoNBpcvXq13bc26t+4bU8kb3zGyRbciu3lEhIS8PLlS7S0tKCyshLh4eH2Dol6GW7bE8kbn3GyBVfsiPoAg8GAFy9eSJ/r6+uh1+vh4eEBf39/btsT9XF8xqm7OAghhL2DIKKOlZeXY+7cue2Ox8XF4dixYwCAf/75Bzk5OdK2fWFhIVd4ifoIPuPUXVjYEREREckE37EjIiIikgkWdkREREQywcKOiIiISCZY2BERERHJBAs7IiIiIplgYUdEREQkEyzsiIiIiGSChR0RERGRTLCwIyIiIpIJFnZEREREMsHCjoiIiEgmWNgRUZ8VGRmJpKQke4fx1/p6Hv91/KmpqVi6dOl/dj2ivoSFHRER9Sl6vR4ajcbeYRD1SizsiKhfa21ttXcI9Bsd3ZcHDx6wsCP6AxZ2RGSzq1evQqfTwc3NDZ6enoiOjkZtba3UHhkZicTERKSnp8PDwwPDhg3D/v37LcZoamrCqlWroFar4evri/z8/HZbegEBASgoKLA4T6PRtBvL2rh+xZaQkICkpCR4eXkhKirqt2NFRkZi8+bNSEpKgru7O3x8fFBcXAyj0Yh169bB2dkZwcHBuHLlinROS0sLEhMT4e3tjUGDBkGn06G6utpiXKPRiDVr1sDJyQm+vr7Izc1td22z2YysrCwEBgZi8ODBCA0NxdmzZ38bZ9t4O5pza+ayKzn/8uPHDyQkJMDV1RVeXl7Ys2cPhBBW52TtfXn9+jU+fPiA0NDQDueDqL9iYUdENjMajUhJScG9e/dQVlYGhUKBZcuWwWw2S31KSkqgVqtRWVmJ7OxsZGRkoLS0VGpPSUlBRUUFLl68iNLSUty6dQs1NTU9Htev2BwdHVFRUYGioqI/jldSUgIvLy9UVVVh8+bN2LRpE5YvX46ZM2eipqYGixYtQmxsLJqbmwEA6enpOHfuHEpKSlBTU4Pg4GBERUXh06dP0phpaWm4efMmLly4gOvXr6O8vLxd3llZWTh+/DiKiorw6NEjJCcnY/Xq1bh582aH+Xc259awNee25w0YMABVVVU4dOgQ8vLycPToUZtysua+6PV6uLq6IjAw0Ka8iPoNQUT0l96/fy8AiIcPHwohhJgzZ47Q6XQWfaZPny62bdsmhBDi69evQqVSiTNnzkjtX758EUOGDBFbtmyRjo0aNUrk5+dbjBMaGir27dsnXadt/87i+nXOlClTOs3p3zn8+PFDqNVqERsbKx178+aNACDu3LkjDAaDUKlU4sSJE1J7a2urGD58uMjOzhZCCNHU1CQcHR3F6dOnpT4fP34UgwcPlvL4/v27GDJkiLh9+7ZFPPHx8WLlypVWxyuE5Zx3NpddybnteSEhIcJsNkvHtm3bJkJCQqzOydr7kpmZKWbPnt1pP6L+iit2RGSz58+fY+XKlQgKCoKLiwsCAgIAAA0NDVKfyZMnW5zj6+uLd+/eAQDq6upgMpmg1WqldldXV4wdO7bH4wKAsLAwq8Zrm4NSqYSnpycmTZokHfPx8QEAvHv3DrW1tTCZTJg1a5bUrlKpoNVq8fjxYwBAbW0tWltbER4eLvXx8PCwyPvFixdobm7GwoUL4eTkJP0dP3683bZyR/EClnNuLVtybmvGjBlwcHCQPkdEROD58+f4+fOn1TlZc1/0ej23YYk6MMDeARBR3xMTE4NRo0ahuLgYw4cPh9lsxsSJEy1eeFepVBbnODg4tNsS7YxCobB4TwsATCbTX8UFAGq12qrr/y6Htsd+FTK25tURg8EAALh06RL8/Pws2gYOHNjhuR3NubVz2RM5W5uTNfdFr9dj8eLFVl+bqL/hih0R2eTjx494+vQpdu/ejfnz5yMkJASfP3+2aYygoCCoVCqLHxY0Njbi2bNnFv2GDh2KN2/eSJ+/fv2K+vr6Hovrb4wePVp6P+wXk8mE6upqjB8/XuqjUqlQWVkp9fn8+bNF3uPHj8fAgQPR0NCA4OBgi7+RI0d2OT5b5rIr2uYEAHfv3sWYMWOgVCq7LaempibU1dXxF7FEHeCKHRHZxN3dHZ6enjhy5Ah8fX3R0NCA7du32zSGs7Mz4uLikJaWBg8PD3h7e2Pfvn1QKBQW23nz5s3DsWPHEBMTAzc3N+zduxdKpbLH4vobarUamzZtknLy9/dHdnY2mpubER8fDwBwcnJCfHw80tLS4OnpCW9vb+zatQsKxf+/Yzs7OyM1NRXJyckwm83Q6XRobGxERUUFXFxcEBcX16X4bJnLrmhoaEBKSgo2bNiAmpoaHD58WPrFb3fl9ODBAyiVSkyYMKHb4iaSGxZ2RGQThUKBkydPIjExERMnTsTYsWNRWFiIyMhIm8bJy8vDxo0bER0dDRcXF6Snp+PVq1cYNGiQ1GfHjh2or69HdHQ0XF1dkZmZ+cdVpu6K628cOHAAZrMZsbGxaGpqwrRp03Dt2jW4u7tLfXJycmAwGBATEwNnZ2ds3boVjY2NFuNkZmZi6NChyMrKQl1dHdzc3DB16lTs3Lmzy7HZMpddsWbNGnz79g1arRZKpRJbtmzB+vXrpfbuyEmv12PcuHGdbkkT9WcO4t8vXRAR2YHRaISfnx9yc3OlFS4iIrINV+yIyC7u37+PJ0+eQKvVorGxERkZGQCAJUuW2DkyIqK+i4UdEdnNwYMH8fTpUzg6OiIsLAy3bt2Cl5eXvcMiIuqzuBVLREREJBP8dydEREREMsHCjoiIiEgmWNgRERERyQQLOyIiIiKZYGFHREREJBMs7IiIiIhkgoUdERERkUywsCMiIiKSCRZ2RERERDLBwo6IiIhIJljYEREREcnE/wCnHSylxtmpkQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# compute the expected number of galaxies in each pixel\n", "nbar = ARCMIN2_SPHERE / npix * n_arcmin2\n", "\n", "# normalise the maps by the expected number of galaxies in each pixel\n", "she /= nbar\n", "num /= nbar\n", "\n", "# get the angular power spectra from the galaxy shears\n", "cls = hp.anafast([num, she.real, she.imag], pol=True, lmax=lmax, use_pixel_weights=True)\n", "\n", "# get the theory cls from CAMB\n", "pars.NonLinear = \"NonLinear_both\"\n", "pars.Want_CMB = False\n", "pars.min_l = 1\n", "pars.set_for_lmax(lmax)\n", "pars.SourceWindows = [\n", " camb.sources.SplinedSourceWindow(z=z, W=dndz, source_type=\"lensing\"),\n", "]\n", "theory_cls = camb.get_results(pars).get_source_cls_dict(lmax=lmax, raw_cl=True)\n", "\n", "# factor transforming convergence to shear\n", "ell = np.arange(lmax + 1)\n", "fl = (ell + 2) * (ell + 1) * ell * (ell - 1) / np.clip(ell**2 * (ell + 1) ** 2, 1, None)\n", "\n", "# the noise level from discrete observations with shape noise\n", "ntot = nbar * npix\n", "nl = 4 * np.pi / ntot * sigma_e**2 * (ell >= 2)\n", "\n", "# mixing matrix for uniform distribution of points\n", "mm = (1 - 1 / ntot) * np.eye(lmax + 1, lmax + 1) + (ell + 1 / 2) / ntot\n", "mm[:2, :] = mm[:, :2] = 0\n", "\n", "# the shear pixel window function for HEALPix\n", "_, pw = hp.pixwin(nside, lmax=lmax, pol=True)\n", "\n", "# plot the realised and expected cls\n", "plt.plot(ell, cls[1] - nl, \"-k\", lw=2, label=\"simulation\")\n", "plt.plot(ell, pw**2 * mm @ (fl * theory_cls[\"W1xW1\"]), \"-r\", lw=2, label=\"expectation\")\n", "plt.xscale(\"symlog\", linthresh=10, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", "plt.yscale(\"symlog\", linthresh=1e-9, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", "plt.xlabel(\"angular mode number $l$\")\n", "plt.ylabel(\"angular power spectrum $C_l^{EE}$\")\n", "plt.legend(frameon=False)\n", "plt.tight_layout()\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "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.12.7" } }, "nbformat": 4, "nbformat_minor": 0 }