diff --git a/2 - Driftcorrection.py b/2 - Driftcorrection.py index e49ce78..df059c6 100644 --- a/2 - Driftcorrection.py +++ b/2 - Driftcorrection.py @@ -97,8 +97,14 @@ def plot_stack(images, n, grid=False): n=widgets.IntSlider(original.shape[0]//2, 0, original.shape[0]-1, 1, continuous_update=False) ) +# Get extent, by any means necessary. Example with center and square of fftsize +center = [dim//2 for dim in original.shape[1:]] +extent = (center[0]-fftsize, center[0]+fftsize, center[1]-fftsize, center[1]+fftsize) +fftsize = max(extent[1]-extent[0], extent[3]-extent[2]) //2 # used in calc half matrices + # Step 1 to 3 of the algorithm as described in section 4 of the paper. -sobel = crop_and_filter(original[Eslice,...].rechunk({0:dE}), sigma=3, finalsize=2*fftsize) +sobel = Reg.crop_and_filter_extent(original[Eslice, ...].rechunk({0: dE}), extent, sigma=sigma) +#sobel = crop_and_filter(original[Eslice,...].rechunk({0:dE}), sigma=3, finalsize=2*fftsize) sobel = (sobel - sobel.mean(axis=(1,2), keepdims=True)) #.persist() sobel diff --git a/StatRegistration_multiNLP.ipynb b/StatRegistration_multiNLP.ipynb new file mode 100644 index 0000000..e7bce53 --- /dev/null +++ b/StatRegistration_multiNLP.ipynb @@ -0,0 +1,859 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1bb538aa", + "metadata": {}, + "source": [ + "This is statregistration.py but then for multiple .nlp's. Edited from Tobias' script:\n", + "\n", + " Calculates shifts from CPcorrected and create a shift corrected stack. The calculations implements the drift\n", + " correction algorithm as described in section 4 of the paper.\n", + " https://doi.org/10.1016/j.ultramic.2019.112913\n", + " https://github.com/TAdeJong/LEEM-analysis/blob/master/2%20-%20Driftcorrection.ipynb\n", + "\n", + " It uses cross correlation of all pairs of\n", + " images after applying digital smoothing and edge detection filters to align\n", + " Low Energy Electron Microscopy images with each other. When applied correctly,\n", + " this allows for sub-pixel accurate image registration.\n", + "\n", + " Config parameters:\n", + " SAVEFIG, boolean whether to save the figures\n", + " stride: A stride larger than 1 takes 1 every stride images of the total dataset, This decreases computation time\n", + " by a factor of stride**2, but decreases accuracy\n", + " blocksize: dE is the blocksize used by dask, the number of images computed for at once.\n", + " fftsize: The size of the image for which the drift correction is calculated\n", + " startI: Starting frame for which the drift correction is calculated\n", + " endI: Ending frame for which the drift correction is calculated\n", + " sigma: the gaussian width over which the images are smoothened\n", + "\n", + "Added napari to choose a rectangular patch on which to perform drift correction." + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "52f39001", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\wgstam\\Anaconda3\\envs\\pyL5\\lib\\site-packages\\distributed\\node.py:160: UserWarning: Port 8787 is already in use.\n", + "Perhaps you already have a cluster running?\n", + "Hosting the HTTP server on port 62459 instead\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "
\n", + "
\n", + "

Client

\n", + "

Client-ae1d0953-8ffb-11ec-88dc-2cf05d1358fc

\n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "
Connection method: Cluster objectCluster type: distributed.LocalCluster
\n", + " Dashboard: http://127.0.0.1:62459/status\n", + "
\n", + "\n", + " \n", + "
\n", + "

Cluster Info

\n", + "
\n", + "
\n", + "
\n", + "
\n", + "

LocalCluster

\n", + "

43df7dec

\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + "
\n", + " Dashboard: http://127.0.0.1:62459/status\n", + " \n", + " Workers: 1\n", + "
\n", + " Total threads: 6\n", + " \n", + " Total memory: 19.87 GiB\n", + "
Status: runningUsing processes: True
\n", + "\n", + "
\n", + " \n", + "

Scheduler Info

\n", + "
\n", + "\n", + "
\n", + "
\n", + "
\n", + "
\n", + "

Scheduler

\n", + "

Scheduler-920f837d-3a34-4324-85ae-dedc20bba490

\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Comm: tcp://127.0.0.1:62460\n", + " \n", + " Workers: 1\n", + "
\n", + " Dashboard: http://127.0.0.1:62459/status\n", + " \n", + " Total threads: 6\n", + "
\n", + " Started: Just now\n", + " \n", + " Total memory: 19.87 GiB\n", + "
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "

Workers

\n", + "
\n", + "\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 0

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:62474\n", + " \n", + " Total threads: 6\n", + "
\n", + " Dashboard: http://127.0.0.1:62475/status\n", + " \n", + " Memory: 19.87 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:62463\n", + "
\n", + " Local directory: C:\\Users\\wgstam\\PycharmProjects\\pyL5\\solidsnakephysics\\dask-worker-space\\worker-loam3fyf\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "\n", + "
\n", + "
\n", + "\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "\n", + "
\n", + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import dask.array as da\n", + "from dask.distributed import Client, LocalCluster\n", + "import scipy.ndimage as ndi\n", + "\n", + "import os\n", + "import time\n", + "\n", + "import napari \n", + "%gui qt\n", + "\n", + "from pyL5.lib.analysis.container import Container\n", + "from pyL5.analysis.CorrectChannelPlate.CorrectChannelPlate import CorrectChannelPlate\n", + "import pyL5.lib.analysis.Registration as Reg\n", + "\n", + "cluster = LocalCluster(n_workers=1, threads_per_worker=6)\n", + "client = Client(cluster)\n", + "client" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "cc591321", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_masking(DX_DY, W_n, coords, dx, dy, shifts, min_normed_weight, sigma):\n", + " \"\"\"Plot W, DX and DY to pick a value for W_{min} (Step 7 of algorithm)\"\"\"\n", + " extent = [startI, endI, endI, startI]\n", + " fig, axs = plt.subplots(1, 4, figsize=(12, 3), constrained_layout=True)\n", + " im = {}\n", + " im[0] = axs[0].imshow(DX_DY[0], cmap='seismic', extent=extent, interpolation='none')\n", + " im[1] = axs[1].imshow(DX_DY[1], cmap='seismic', extent=extent, interpolation='none')\n", + " im[2] = axs[2].imshow(W_n - np.diag(np.diag(W_n)), cmap='inferno',\n", + " extent=extent, clim=(0.0, None), interpolation='none')\n", + " axs[3].plot(coords, dx, 'x', label='dx')\n", + " axs[3].plot(coords, dy, 'x', label='dy')\n", + " axs[3].plot(shifts[:, 0], color='C0')\n", + " axs[3].plot(shifts[:, 1], color='C1')\n", + " axs[3].set_xlabel('frames')\n", + " axs[3].set_ylabel('shift (pixels)')\n", + " axs[3].set_box_aspect(1)\n", + " axs[3].legend()\n", + " \n", + " axs[0].set_ylabel('$j$')\n", + " fig.colorbar(im[0], ax=axs[:2], shrink=0.82, fraction=0.1)\n", + " axs[0].contourf(W_n, [0, min_normed_weight], colors='black', alpha=0.6, extent=extent, origin='upper')\n", + " axs[1].contourf(W_n, [0, min_normed_weight], colors='black', alpha=0.6, extent=extent, origin='upper')\n", + " CF = axs[2].contourf(W_n, [0, min_normed_weight], colors='white', alpha=0.2, extent=extent, origin='upper')\n", + " cbar = fig.colorbar(im[2], ax=axs[2], shrink=0.82, fraction=0.1)\n", + " cbar.ax.fill_between([0, 1], 0, min_normed_weight, color='white', alpha=0.2)\n", + " axs[0].set_title('$DX_{ij}$')\n", + " axs[1].set_title('$DY_{ij}$')\n", + " axs[2].set_title('$W_{ij}$')\n", + " plt.show()\n", + " return min_normed_weight" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "05476d21", + "metadata": {}, + "outputs": [], + "source": [ + "folder = 'D:\\\\20220210-36-CuGrKalbac-old\\\\growth'\n", + "#folder = 'D:\\\\20211130-27-Si111SbPLD\\\\PLD2_100mJ_400C\\\\growthIVs'\n", + "names = [f.name for f in os.scandir(folder) if f.is_file() and f.name[-4:] == \".nlp\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fed2f705", + "metadata": {}, + "outputs": [], + "source": [ + "for name in names:\n", + " script = CorrectChannelPlate(os.path.join(folder, name))\n", + " script.start()" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "id": "d551cc5a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Array Chunk
Bytes 248.48 MiB 3.65 MiB
Shape (68, 1096, 1748) (1, 1096, 1748)
Count 68 Tasks 68 Chunks
Type uint16 numpy.ndarray
\n", + "
\n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + "\n", + " \n", + " 1748\n", + " 1096\n", + " 68\n", + "\n", + "
" + ], + "text/plain": [ + "dask.array" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "conts = [Container(os.path.join(folder,f)) for f in names]\n", + "#original = da.stack([cont.getStack().getDaskArray() for cont in conts])\n", + "#original = da.stack([cont.getStack('CPcorrected').getDaskArray() for cont in conts])\n", + "original = da.image.imread(os.path.join(folder+'\\driftcorrected01\\*'))\n", + "m = 1\n", + "subfolder = 'driftcorrected%02d' %m + 'it2'\n", + "#original = original[:,m]\n", + "original" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "534a8b02", + "metadata": {}, + "outputs": [], + "source": [ + "# config\n", + "SAVEFIG = True\n", + "stride = 1\n", + "dE = 20\n", + "fftsize = 256\n", + "startI, endI = 0, -1\n", + "Eslice = slice(startI,endI,stride)\n", + "sigma = 10\n", + "min_norm = 0.4 #minimum" + ] + }, + { + "cell_type": "markdown", + "id": "cf9cc677", + "metadata": {}, + "source": [ + "## Step 0: choose area\n", + "Choose the (rectangular) area on which to perform drift correction." + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "id": "29c1506f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(292, 804, 618, 1130)" + ] + }, + "execution_count": 62, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "center = [dim//2 for dim in original.shape[1:]]\n", + "extent = (center[0]-fftsize, center[0]+fftsize, center[1]-fftsize, center[1]+fftsize)\n", + "extent" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "id": "2b532e70", + "metadata": {}, + "outputs": [], + "source": [ + "viewer = napari.view_image(np.swapaxes(original, -1, -2), name='original')\n", + "\n", + "# create the square in napari\n", + "center = np.array(original.shape[1:]) // 2\n", + "square = np.array([[center[1]+fftsize, center[0]+fftsize],\n", + " [center[1]-fftsize, center[0]+fftsize],\n", + " [center[1]-fftsize, center[0]-fftsize],\n", + " [center[1]+fftsize, center[0]-fftsize] \n", + " ])\n", + "shapes_layer = viewer.add_shapes(square, shape_type='polygon', edge_width=2,\n", + " edge_color='white')\n", + "shapes_layer._fixed_aspect = True # Keep it square" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "id": "697eb9c4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The extent in x,y is: (214, 726, 646, 1158) pixels, which makes the largest side/2 256 pixels.\n" + ] + } + ], + "source": [ + "# load the outer coordinates of napari\n", + "coords = np.flip(np.array(shapes_layer.data).astype(int)[0])\n", + "extent = np.min(coords[:,0]), np.max(coords[:,0]), np.min(coords[:,1]), np.max(coords[:,1]) #xmin, xmax, ymin, ymax\n", + "fftsize = max(extent[1]-extent[0], extent[3]-extent[2]) //2 #This is basically for \n", + "print('The extent in x,y is:', extent, 'pixels, which makes the largest side/2', fftsize, 'pixels.')\n", + "viewer.close()" + ] + }, + { + "cell_type": "markdown", + "id": "322a68b7", + "metadata": {}, + "source": [ + "## Now starting the steps of the algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb6da5b0", + "metadata": {}, + "outputs": [], + "source": [ + "def crop_and_filter_extent(images, extent, sigma=11, mode='nearest'):\n", + " \"\"\"Crop images to extent chosen and apply the filters. Cropping is initially with a margin of sigma,\n", + " to prevent edge effects of the filters. extent = minx,maxx,miny,maxy of ROI\"\"\"\n", + " result = images[:, extent[0]-sigma:extent[1]+sigma,\n", + " extent[2]-sigma:extent[3]+sigma]\n", + " result = result.map_blocks(filter_block, dtype=np.float64,\n", + " sigma=sigma, mode=mode)\n", + " if sigma > 0:\n", + " result = result[:, sigma:-sigma, sigma:-sigma]\n", + " return result" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "id": "d94f7854", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Step 1 to 3 of the algorithm as described in section 4 of the paper.\n", + "sobel = crop_and_filter_extent(original[Eslice, ...].rechunk({0: dE}), extent, sigma=sigma)\n", + "sobel = (sobel - sobel.mean(axis=(1, 2), keepdims=True)) # .persist()\n", + "\n", + "# Step 4 of the algorithm as described in paper.\n", + "Corr = Reg.dask_cross_corr(sobel)\n", + "\n", + "# Step 5 of the algorithm\n", + "weights, argmax = Reg.max_and_argmax(Corr)" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "id": "686cf0bf", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "memory backpressure assumed\n", + "39.07800000003772 seconds\n" + ] + } + ], + "source": [ + "# Do actual computations\n", + "t = time.monotonic()\n", + "W, DX_DY = Reg.calculate_halfmatrices(weights, argmax, fftsize=fftsize)\n", + "print(time.monotonic() - t, ' seconds')" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "4877c08f", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAADGCAYAAACn6KM+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACMK0lEQVR4nO2debxkRXm/n7e77zL7MAvbLMwwbALiMDOKIiAiATQqQUUBRVQSwIDgz5goIYlEgyGJGlBUFjdcEBExoiDIIoyIbDPAsA0wG8zCMvt6t+5+f39U1enuun26+95e77313E9/6p46dc6p013Vp6u+9b6vqCqBQCAQCAQCgUAgEGg+iWZXIBAIBAKBQCAQCAQChjBACwQCgUAgEAgEAoEWIQzQAoFAIBAIBAKBQKBFCAO0QCAQCAQCgUAgEGgRwgAtEAgEAoFAIBAIBFqEMEALBAKBQCAQCAQCgRYhDNACgUAgEAgEAoFAoEUIA7RAIBAIBAKBQCAQaBFSza5AoHaIyG7AJmCnzdoEPAJcqqpP2zJzgCXAfqr6is37KPA/wBGqurrhFQ8EBkm5Nm/3vwy8V1XvzzvuJ8AY4IOqqg2u9ohmPxHdVabMK3Cnqp7UkAoFWh4ReRL4d1W9xW6/CXgCOEtVf2zzDgEeAmao6pYmVXVYc+JJh+nGDTti9y9atDL020CgRoQB2vBiLrBJVScDiMhk4O+Bh0VkvqouVdXlIvI74LPAF0TkbcBVwIlhcBYYgsylTJsXkWuB/wfcb8v8K3AwcHQYnDWeLuD8MmX+BaY0oi6BIcMWYELe9kXAVi/vM8CPw+CsfmzYsI2/PPxvsfs7UmeHfhsI1IiwxHF4MRczqwiAqm5U1a8Ai4Gz88r9F3CuiBwK3AKcp6qPNLCegUCtmEv5Nv8N4EQRmSMipwLnAO9TLSvkBOqAYB48pV6BgMdWYDyAiEwF3gf8EDtAE5GJwEeBbzWpfiMDBdVM7CsQCNSOoKANLw4n78dqHkuBaW5DVReLyCPAw8B/qeovGlO9QKDmlG3zqrpWRH4OfBuYj1GL1zWshoF+hAdPYIBEAzTgXOAXmKXLe9m8TwF/VtWlTajbiEHJksl2N7sagcCIIExWDi/mUvzH6gRgvdsQkQSQAbIYNY28fX8vIgeIyN9ahS0QaGXmUkGbx6powPmquthlhvbeeIKCFhgEW4HxItIGnIdRyrYBE+zz7HzgSsj1aft/6Nc1RVFNx74CgUDtCM/CYYKIdABvAJ708pPAkVj7G8vXgYnAi5hlIRGq+h1VfUFVv+cciwQCrcgA23w70INZ0hsR2nvjEYyCVupV9hwiPxCR10Xk6by8S0VkrYg8YV/vydt3sYgsE5HnReTEmt5QoBE4Be1DwJOq+jx2gAb8NZAG7oBcn7b/h35dS1TRTHfsKxAI1I4wQBs+HIpRxJ7z8s8DeoHfAojIucApwN9g1LN/FBFxhUXkTza9r+41DgSqo6I2b3kT8LR607yhvTeHGihoPwKKeYv7X1Wda1+3A4jIwcBpwCH2mO/YQXxg6OAGaBdhlTLMAG08xjnIt5zDH9en7f/3Nbaaw52goAUCjSKYAgwfDgeeUdU+ABGZgVmrfx7wHlXtE5Hjga8Cx6rqayJys90+Gfg/EZkCvC4i44DtTbmLQKByyrb5vLJz8ZZChvbeHJyCVg2qulBEZlVY/GTgRlXtAVaKyDLgLcBfqqxGoHFsBY4GtqnqH2zeNky/HgV8AHJ92v4f+nWNUc2ima5mVyMQGBEEBW34MBc4TES2i8hm4G5gN2CBqj4iIgcBNwJnqupTAGrcLn0D+II9x2GYGGmHAmFZSKDVmUuJNu+VdXGT8gntvQlUuMRxiog8lvc6p8LTXyAiS+wSyN1s3jQgP4TIGvKcJgWGBFsxn9mVeXnbME5CfqSqLjiX69MQ+nXt0Syku+NfgUCgZgQFbZigqhcAF5TYv5QisYVU9dsY73aQe7i9EXiqDtUMBGpGuTbvlT22SHZo702igpnBDaq6YICn/S7wFUBt+nWMdz8pUjbEvxtCqOqNmAnG/Lxn6P/Z5g/QQr+uB2EpYyDQEIKCFsjnjYQfrIGRQ2jvTaAWTkKKoaqvqWpGVbPAdZhljGAUsxl5RacDIczC8MT1afd/6Nc1RDSLpLtjX4FAoHa01ABNRE6yXraWicgXi+wXEfmm3b9EROY1o57DFVU9W1WXYx5sIZ5MAwhtvnmE9t486uFmX0T2yts8hdzytluB00SkQ0RmA/sD/hLYuhD6d2PJ69MQ+nUdUMim41+BQKBmtMwSR+tV69vAX2FmPB8VkVtV9dm8Yu/GPFz3B47ALGk5otF1Hc6IyP8B93kOFgJ1ILT55hPae+OphZMQG3j8WIyt2hrgS8CxIjIXs3xxFcZhDKr6jIjcBDyLccd+vrW/rSuhfzeP0K/rhCqS6Wl2LQKBEUHLDNAwy1GWqeoKABG5EeN9K/9hdjLwY+tO9yERmSgie6nqK42v7vBEVf+m2XUYQYQ232RCe288LlB1Najq6UWyv1+i/GXAZVVedqCE/t0kQr+uF4oEpSwQaAittMSxEk9bwRtXYDgR2nxgRFIPG7QWJPTvwPBCs5DuiX8FAoGa0UrPwko8bVXsjcu6ZT4HYMyYMfMPOuig6moHkM2adNcuAF5anfdcHT0agA0bClfOjBuXtPWBqVP3Yfz4wlMmWmCI7G5LFVavfgnJZthn772hszPKB3MPw5FVq1axYcOGZtxdzdp8sfa+bVvpi7vP1aWZTGHaY5+3mzebtK9vKzbEEFOmtAOFbQcg6YX/3WfatMIdthG58u54l7pr+3XzU/8e4rYrKbt8uTFVmj9rL4riN/x6d4QGdLRFy5dvUNWpdb9QEWqhoA0R6vZMA+YDzJthD8/aNKGF267f4JVTk2bSpl/u6hpldrePiq63zS4O7EjEtEeBvl19TJ+2vnTt1UvFzy9+ftXCe8gd55fPXdB1nYJj20azbvM4erv6mO3qGleXcsR9ev754spH+UXuWbT4vphralw5gZe3ZNmwM1uXLxLJ1n118LBiypQpOmvWrGZXI9CiLFq0KPZZ3EoDtEo8bVXsjUtVrwWuBViwYIE+9sjAbcKz9mdEAvvrsdt6KXra2p4/8ADn/vKXAFz70Bx7VJtNzVu7fbtbDpBk27ar2Lmzs+i10t6qgVSZT8btjzvO31/qfO620mm46KJzae/byTX//d/0Ttm7dCXK1KncPcQRvd+WbAU/5/xjKiVLgre8ZaDevGtGzdp8QXs/4AB97KtfheOPB/q3Y7ftPieXbtli0g0bTLpqlUlvvtmkP/mJW3X1T3xggS00YQIA28Sk45M7C+p1zfvfb/4ZM8ak48aZ1DUOOwkQpRMnFm77DdpPPUq1Fb+vuO0xY94MwGNHv6GwgH+NuAY92IbeDGxd5VvfeqlZVaiFDdoQoW7PtPkzE/rQ59tpv8g04nT2egCy1gW6M7Fz2+m+LWa72w5Q+kzYsPZ1xuFh1/dfBeCpxYcx9sB9Afj9y3sC8JYppqzYEUAyYb5HEqJsfW4V77n4d+aco8wkJe1m8kYTdlLGtjlt7zB1GD3ebtsJwFR7sdvNHS+FraWUDZT0FgZRTvR2I6MO5EvnHcrKJ1fz489dZ3a4x4Udw2jWfm9kCr8/1O3vs/eQ9mag/AFxonA0lbXHxQ6g88tYMr1tBdd2+92ALJE0lU/3thUcl7V1F1FO/N0y6oIqku6tz7mHKbNmzeKxxx5rdjUCLYqIxD6LW+k5+Siwv/WytRY4DTjDK3MrJhDpjRhD6q2tt1bfxdF0M5HugVH8IRQY0QyLNn/Nl7/MitQBAOzb6f22fOCBJtQo0OoMUzHeZ1j0783PrQIg4Q3QklKp7BQYPigEBS0QaAgtM0BT1bSIXADcCSSBH1jvW+fZ/VcDtwPvAZYBu4BP1rNOg1VlAETsLJh2lSlpqFYxi1MHKpnczy9TasllnELmX3uwDOT9ruazaRUa1eb998ptp1KFH3aycMK7SNtxs8dSuEYxnSblhGG/MfR5TtTiGosv51XaqLxKlmoX7n77941EVdesWQeIO38tqXVdB4GQW2cwnKlr/84K2Z4U6azxi5JKnAXklLRIdbLdIZGwalXbOFu3tNthkqT50S2i0TJBt7IxKXZA5qlD4g/Q7Lki5SvKTxamHv3KRxdIFd2vCaPUFV1qlyrMywLS1tG/XPTVZ+5B7Bul3qpASeSXop9C1q/KCXcep2bZ7QHYMkTva8Kd02xrRgr2++9/lF+mjtUgQUEL1IMHroBp82DtYpPOPgZWLsxtr10MR3222bVsOC0zQANQ1dsxD6z8vKvz/lfg/EbXqxTXnHpqtMwxEBgoQ7HNBwLVMIKWOIb+HRhehAFaoB5Mmwe//AQc9TmTTlsAL/8FZr4N/vR1OO1nI3LANlKek7Wl6Cy0m7k3M3huEtwXEAarOpVTxCqddM9XGCqx7aqkDoO1PYuzNYtTfOrBcFDiakVmICtXZs406bZthR94pR9+pfZdw5laS9BDiBHiJKR+JBRpz+RszKxy5itpJhQbSMK2NSlMc3ZeRnVJiEa2ZtGl3D4v3ylrzuYsWn7h2rU9tzobtJRdVWJtz5xDEk16ClfCOz4R891QxN27u59IXUskIU9Bk6T3fe+cFjm1isIvwUgJa3PKfoza595Gq16562TtTyzJOlWsvw1aP9XObou33y0pFasSJlWKHp9IZvqrmzVDkWx4ZgZqzOxj4NQfwS8+BpleePFOQGDFHyHZAQ9eVThgO/QDJi+bLlTehtnAbQT+KgoEAoFAsxhJClogMKxQhaCgBerB7GNg5pHwwu9h9GTYtdEM1sAM2CQBK++HZDtIEv5wCRzwbrj/v+BNp+fSn34I9j0WVtyXSyUB4/YyDx/FpNteMWEj/LL7HtsyA7/wnKw7bkbOTLENdrK8nEpVqXo1UNWskjqU2660LnFqVr29OAYMvov8gdDvsy7X0OOMJkcSI/GeLSPESUj9yAramwTrrdEpYvFKWukv42J2S/2jS1hbKKfc+J+iU9AGaINGjELWTzlzx1t1rJiyJvTYfbasJpFi9l8xHugjxSzhPU+yZZ4Tkct7+x7FebjPFle9Sp/TbmpxZa3fYSrUs4eJBgUtUAdWLjQDMDCDszeeCs/8xgzSOidA91ZId0OqEx77PkyYbgZzUw402/ufBI/9AA44Cf7wL4XpC783xzmj0mzGnLdY2T/8i1li+cevwuEfMwO/d14Cv70oflAXNwDc58iqBnphgBYIBAKBhjFSnIQMVXa+sJwxB8wpXzAw8lCFdF/5coHAQFi50NqezYdVf4K/+g9Y+N+QtEOU7m1w8N/A0t9B9xZAYOsas2/D8yZ98Q6jrr3we2gfa9KOcbnt3h0UBC4cPdnsG7unSSfMMOkeh8LaR41K99j3zXX/+FUzoMv2wayjjXq373Em3f9Eo/AlO8xATSQ3AJwwHR78Zs627tQfDehtCQO0geBLQ9Fa+wR5gVUKDmlr2w2Avj4XH6r0T5PB2prF5bsYZ8X2F1Ph8peXu/3uHHF1jTvvQG3S4tSsSuzm4uzXyjEsbdBUjfFjmQ8gd++F72mFIccM1ouj4+W0iZ03k5cLy8UZuJWThP1KVbq/CO5+223R3rTXliqVhuOolR3dCFDWgm5dJQral8x5Y7Rd2dmc+UpaT/pasz/hbM/sc8jZbJXw/Jf1lJx+9k1+tHnfu6LbjnPP7tuSJZzdlo3p5pSySty7+6paNkPBz5xscVXKR8spZjHEqlrOrqzY+6zetfw6eXXxrxEpmp7CVhdUC38kBAK1YO1iM3i5699g1CTY+00mf9bR1vbsnbD8XjMAIwHZXphxBKx+GMZPg21rc8siXdo50QzmXDpqN+jabM7bPsaUSXXCjlfNebeuNvtee9qkbnXCs/9XWNdVfzLpintN+uKdJvXjMyZSZoDXOcHc177vMMs4B0B4TgYCgUCgYTgbtFKvQCDQiiiS7ot9DTVEZIaI/FFEnhORZ0TkIps/SUTuEpEXbbpb3jEXi8gyEXleRE5sXu2HEUd91gxexu0F4/c2A7bTfmaWCJ72MzjyAlNu1tHQ1mmWIK5+xNigbVsHe801A6695sKuTWbpYffWwrRrs1G5Eu3Qu8vkpXtMmunLDZ6mzTfpvu80ZR1TTKxXJu9n0t1mm3TCDJNOnJUrO3ZPM/E0eoq5/qR9Ye/DB/y2hGdhJfixmXbsMKmTltJpcmPd8UVPMXXqGABGdzrvTgMbGw9W6Wmv4BN28aE6O60dUq9COh1dc3RniYOrvPZAGZaKV61RNYpVhUpMu1V+2tpy7QByglBH5ATNqb+JSBHrGb0bm/Y8GNdEnNfS1zBeHvfosbNKbfZYZ+jWFqMkuz7V2Vk8P07WG0Tgv/aUpyD6/bxSJa0R3hiHkZdLIcwMVosiaDpJpm8LkItz5hQyZ3PmlLOO1DkA7Oz+BuCv84CsVZOzcfZl5Lw5ujK+stb/gCqMWqvFeYBsH422DfIBBv1t1crds1/OV7ucDZqvmhU7RYxyJtaOJlLObLno3JKA+jlxHG4KWhr4B1VdLCLjgEUichfwCeAeVb1cRL4IfBH4gogcjAk4fwiwN3C3iBygqhXIu4Gy7NoEoyfl7LTcoOmBK8xAbe1iM1hbuxj2ehM8+XNY8CmTHvBus0TxgHfDC3dY27M7itugicCK+70yd+SOXXA2PPXLnCFuIgUbXrQ2Z/fHpPeZAaAq7Hgtt++wj8Cyu832AAnPyUAgEAg0lGSZVyAQaEGcDVrca4ihqq+o6mL7/3bgOWAacDJgo75zPfA39v+TgRtVtUdVV2ICzL+loZUezrglij5OYctPJ86EM36RS2e+FU64zChXJ/xHYbr/iXDYaXDYh+GNH4a5p8P+JxQve8J/mKXR2bSZTVxwtpnwSbbnBnX90vvMADCRNHZzruyCT5nBmbNBW7lwQG/H8JmWrSdxs/YFNmiV4ZSzgXo+LOfxsFIbrGLKXb5woIqZNRhGM/YjjkwGtmwxLyjfyOz2pIkT7aZpI3aT3aLFHWNt2gY9XQB0jE4yaWK2X7uOcHVwO3yFzE9dOV/Gc6mf791DPyoxvnTKoKtrXLlyfaKV+kwr1cUjOAmpAVkh091Btns9ANo2zuR73hqdouaUszGdnyvYzoyaCECyw/y4bkvaOJ6JTKShjUqZfe2pGJGg26rk9jnYT2NqtxK8VV7cE8gFPHbx0Hw01V54Pt8TpfZXrKMgynnqnWQnk2UiWUD7fFXKeXWz6lPGez7a/dkeG8Mt1nOivXcpfP5qJlmQFiPbV7gv01d4nxkbe81X1LLu3F7dMypFFdBaICgyTG1kRWQWcDjwMLCHqr4CZhAnIrvbYtOAh/IOW2Pz/HOdA5wDMNPFDA2UZ9dGY4NWCb7K5lK3HNJPS+GXfeAKM5g79INGrfvozfD0r4xnRjeIW3FfLt3/xOJeHCfONLZ1zsZu7eIB2aG17lM8EAi0KN8BPtnsSgSGKGGJYyAwRBl+SxwBEJGxwK+Az6rqNvFjTOQVLZLXb0Gpql4LXAuwYMGCei04HV5kM8ZOrJiC1mjyXeH7A0DHQAaAcecoQxigFcF5d4sEMqdKuYyxYwu3EwlyC3PMLGOfp/Zv325m4ArM1ihv1hJP3E+cwrpXMtmVr6BlMpBKpaCzM3ofypncDHayvty9u/c9xCobIE5B27DBbJdTlzxVarxt3xMnmvd9wgRXLE/3sA35mu98B3bsIBHZdXnX8lUp1zHilDRHnMfUOLeklShoscfYdPPmwv2+nZwfKC7Ojq6RtLBSVopqay0iPwDeC7yuqofavEnAL4BZwCrgw6q62e67GDgb8wV9oareWWUVmosK2peCPmMPHXlz9FQm563R/aT2lbStr5qJlpRV0BJWARJREs78wtqeJSXmh3mX97DzfsBLr1G1pNOUU7tf7ZK4yLmE89Lo7MdSxfuXi3EmRbw6unNpvoKmO0lgnozZLqvmJVwdCvu09hVvmRmnoJVQwqC/l0Zn21daQSu8ZjZTaGPmFDXfW6NT3px6546ThNbRBk0hXZ25lYicBFyJ+dH0PVW93Ns/AfgpMBPzVfE1Vf1hVRctXZ82zODsZ6p6i81+TUT2surZXsDrNn8NMCPv8OnAunrVbUTRtYXI/X0ACBOZgUAgEGggTkEr9aqAHwEneXlfxBj27w/cY7fxDPtPAr4jzh99oCSrl6xpdhUCrUY2G/8qg+133wbeDRwMnG77Zz7nA8+q6puAY4Gvi0g7dUCMVPZ94DlV/UberluBs+z/ZwG/ycs/TUQ6RGQ2sD/wSD3qNuLYtdGkYYAWMTSnXwdINgu7uss/9uOcuDmbnNFult8paC5tawPcLGLxeGfd3WbWafv24h6l/DBR/mR9HG4S3wkT/qS+n+8re/6102noaDM2aM5ZpX9MuWvGXbtScqJA5Z9ZNUJCXIiuIUtPDyxbBrNmme24xuQ+GPfmOaMzm+5rj581y3wO+++fd+xO287XrCk8h/9BLFtWeO442zNXR1du1KjC8/neH/3rlbvHYkR1HW2S5Y8V7vdVvXKdshUUtSFCtaMjVV1obUbyORnzgw6MYf99wBfIM+wHVoqIM+z/S5XVaBqZdJKd6ycycd0Sk2Htv9T3nOhtO5szp5xN2NMIE73fMuXGjTH9unNUV3SMW+3l4p/5scMyG027T3SYL+OcUmaPc92swyo8nda2rN3+5m6zqfuB72y6Y77U/XssUNLcA8GdI5GAcROBN5j73Dyu8FyeTVk2nSq6P9PrbNCKP5Pce5KwNnyunK+GFcNd072/GU/VS6cLFbSEVf+yJa6R9W3pakX1CtpbgGWqugJARG7E9M9n868CjLODp7HAJoy3xXrwduBM4CkRecLm/TNwOXCTiJwNvAycCqCqz4jITba+aeD84MGxRkQDtApt0EYAI2KAFggEAoHWQKhogDZFRPJHzNdau45SVGXYHwgEKkBLKmXl+u00YHXe9hrgCO8cV2GUqnXAOOAjqqUvOlhU9QHiY3u/K+aYy4DL6lGfEU3XJpMGBS1iRAzQdu2CJ57oPwEep1o5ccAxxoQwY9w4Myt18EEHmQxnD3PbbRiTBwDnRtOfVDEnWbTo3QW5lSpAvmBQaSimUvZivmKYThuzpTEz22DsWJ64r/ixcXZt1YSmahTF6rhrV3PqUjc2bYJf/rK//ZePb4PmKWhOgfvIqR8C4Pnnc7Oya9eYZfhnf/1KvvWta6Ku4FLXl/a4+ebCc/qKma+kOZeRcV4c49K4eyuVF22bmbtl27cX7E7428WvVPH+Shnu685doOoybFDVBTW8pM+QNtzf1TWKxU8fwtu+/yiQU26wKoyzhXKps4Vy3hqdzZlTzto/Y47vvdKOY6f2wNL9mXHYdHa8Zn44pZ1i49lCrXv4UADa2sw5kynvAWHrkLReIFPt1t6tLV2Q+sN2SRaex91Loq2rID+bZ9/l25BJso/EtDTb0mk29aV58fFDTL6nBkb2XjEKWVdf8RV2/vE+mci7Y3H1EaDPU7vSWXM/SauUuTpl7LHt9rP281OJ3Bimu++5ovWpGlVIlxwrleu3lfTFE4EngOOAOcBdIvInVd02kKoGhhhhiWM/hvtvgcAA+M53zm12FQKBhnLu/fc3uwojkhrYoBXjNWvQTzDsDwAkpr2x2VUYfmQ1/lWeSvriJ4Fb1LAMWAkcVJO6B1qXsMSxHy2kbdSPHTvgT38qbyLiJtS77ASdb+7iTM66rT3bvLlzTcbkybwH4/zn9ii2oa+gmRm4hx56N5UQN9kfp0pVqlblq4hOQXTHrjchddhnJpx73nlMm164oqhS740D3e/qFGcLVsr0Z6D2Y8XO5SumQ51d6TSPvfoq+/7kJ0D5H7zuY3IWVymnXjmjs6VLAfi3f/mX6JhzzrOztn1djN7xOqMjO7ZCu62F1hCx3TYuNw/t+U+M6tgZk18uJWa7GPnHmtqdCcDZPGT/G9i5GznLNRxm1CpU0AaDM+y/nP6G/TeIyDeAvRkGhv09mSQrtk1k7OLDgJxC4zwu+jjVy8U5c94anc2ZU87aLzJfzvPmZXjfRHPsxh7TKzv6ch4e83l21Wwg5+Ux6VQ7ZzPlVD2btltlLJUXc60Ybr9TndzxKU+hy2ZzX+oZzx5rzPZ9AdhGN3vMncoTr/YWvZYjTgnr9dTD6Nruuu5e/Xxrw+beE/9eALoziYJjM3aXs/1z+X322NFJc/Y+e253rbaoDUB3GW+Tg0UVNF2V+PwosL91sLEW47znDK/My5jlhX8SkT2AA4EV1Vw0MATYtdEEe24b3eyatAwjYoAWaDy33npubPxu3xbcPYg05ns/PiRJ/DFx+Of6wAeuGdgJAoEY6uYHOo9ahJ/9VA3OUQ21iIMmIj/HOASZIiJrgC8RDPsDdWTqITNjB2jZmAGaE5US3jMu/znUawda4h0TOWex+Wl7jU67lNFtu0dgyg7Q1j+Tb+JVY5TczQ3mcNW0iFwA3IlZ0/oD2z/Ps/uvBr4C/EhEnsLc/hdUdUO1VQ+0OLs2meWNpX7wjTBaYoBWKn5NXpkZwI+BPTFfEdeq6pWVnH/z5pe45RYzYDj11PCDvBi33tp/eaOf5wZU/gDLZxjGsaw59W7zXwXG2/+vqEmNhx8XFsn7Sd7/lT4mwuNk4FT74FHV02N2tYRhf737d9uodqYeMpOxu41hx/NBXCjG6AP2wx9NTD54n5LHxE34pXU4aNc1osqpDVW9Hbjdy7s67/91wAnVXSUw5Ni1MSxv9GiJARq5+DWXi8gX7fYXvDJp4B9UdbGIjAMWichdqvqsfzKf0aNh7puU7uuu48mHcsv2/HGE+24WL3VfzW5hxc9tejfvsP99HM3Y7xfnm96t48uLTH3uxRfz2qvnlKtuUbZed13RfN9pv7ungS5CeEP+xuFnAzB/3vB/KH3hH7P88pdNuXTd2vzo+fOZOm8eL682Lfg02xj2mVm8VaQzplxPj9l2TdYt/ZyxZi0Act550THXXm3a+4pVCVbs6B+A3fGzT1m9ptK4EZbNW+zyHdug3Q+nKZNr59vBrXj+oO1bkz/60dIH+DMScTMUtcSXl/3ZRf8XZQWzjx+78krOdNHHm0Adlzi2EnV9pm3rg3tfh2Xb94RRe0YKjSOuGUTPtEiWseltfw3Al44zmy89uYZ/X3wxABvPM8sE+3rMAmUXIBkgNesN3LhoPgBJb+lj2ipDPbYPu6V7vWV+4Perup/hfcfkX9V5zU8KsCxPXcLYLSzdVvjEzHg/AvzfAO59ysR87cQN6DLe/lLd0ilmcSZcLjs3MWq+d3xfHa7civUpdvTVadpIQdNhSmrI8sAVMG0ezD4ml7dyIaxdDEd9tlm1MjgFLRDRKr/AT4bIeOt64G/8Aqr6iqoutv9vB54juEoODF1Cmw+MSGoUqLrVCf07MDzJSvwr0NpMmwe//AQ8+xsze7ByodmeNq/ZNbMKWhig5dMqE5lx8WuKYgOUHg48XMnJsxs2sP266xjl5fsTVrnZtuL57uvHet1nPsYD3Pu5n207jFIw3nkSceQrae3t8OqrZttFcbazYbtuv51ixNUpOr1Xtzjc/vx79u9/jK3LT39uZkg/9s619uL26v7axri1jOXWQPr45fzzljqPPzXpvIY4xcZNYcbV/b77wHOn3iDq1+Z7e2HnTiZMMG3R3erWbZ79hM3fbbSRzlKjzQx5W5sLiGqvnTGtrCvfhfWqVUAumHWvdePdLwSD/Rycsbzb39FnlWbnicd9XjYdb9dn+h9fT2/xexjIslp3zrHtxmFApAe8bp3++caRjnLtuVZr5+utzD3xRH3PXwHJcu/VQI1LW4+6PtM6EsLsMW3Mn2T6kXP64TsJ8W2mOqyDjpzjDpO/q8/0w6295jtg76m7+PqpRjmbfLVZQrmr95sApHu3ROdLJPZk748ZJ3wJ+0RJ9us2pm4pcWmh05A4Ihf23tPNKU3FjvbLOlavMfc3sa3wqLjfAP4qGvXSOHJ2ZIUOPBy+IxDIU/GcUxBx5zCp88nhtlNuf0wdekYL7fUKHaOCpofJFMpIZPYx8N4r4KYzYc674JUn4NQfFSpqzSIM0PrRsAGaiNyNWWvvc8kAzzMW+BXw2VJxMUTkHOAcAOcT5lteGf8Lzm37PmRcvl3FxdM2vbvCOjuuufLK3I8jNyiwP1BfjRmgxX0Ju/w4b3Y+vvEy5AZ3bt+/lzlHYGA0ss3nt/eZ03KT8F/+8jVlPX8evKcNEGljlW3bYVrEmjV2/yzztH9te17P2FmZ3cs1dimkG8C5FcCTdrxs/nETGl48NH/A56cON88xEG+e7hyXfemTlR8UqB0i5d3Nug+2hWnmM21i50QANj+3Csh5UCw3QBuV6iso5wZoO+zAzHlsrJRs9lU2PPOyrYMWpA4Xz6vN1rEtUdwbpI9Ggxx/gOZijxWpT17ZSWXszQKDJMbTZWCIsOrPJl1+D7z9s2Zw1uxljpk0dG+BUcEGLZ+GDdBU9fi4fSLymojsZWca8+PX+OXaMA+yn6nqLWWudy1wLcBk+yTwXXyX9btvSdsfC+64rJdCfPDmfj9E3A9R7wB/wIW37Q+mYutq0+kuMLBP/j07wyGHNUKK1Aj//fGVBd+uyD0x4/J94mbRB2Kv5Cti/rFxdXV0dtbNa1Aj23x+e1/wxjcq2SwTZBsHpFZAZ1yAZkdn0d1RMdtO2vLd33pGZ+1RYa+FevtTqURBftzIqz0V1yOKs4cUfftKDgSuufJKzr3oov6z4nEKbpxM5/bXSvWph5edfGOY2C+rBiGS+x6MYwgM0Jr5TJs+frpCnmrlgjhTXJ1yA7J2GyzaDehcOReE2rnSh5zNmVPORrcbtzpdfd8BIJvpLjiHP/BSz7Nh5JWwTEiACDeI9IJBR8cV+dp2zjxSkmXn0hVMfMMsU1f7hHWDx8iuyzunc5efs0ErDBlQbCVKsfJ+vk8yb4/zP+KXdUpawpP1ik225pMSrZ8jPAVN18eFf6BBZPJCTTxyDYyZCg98wyhpzaLL+k8KCloBraJVu/g1UBi/JkJEBPg+8JyqfqOBdQsE6kFo84GRiVPQSr2GPqF/B4YfWsL+LNigleeBK4xalc/KhSa/EaxcCE/dZP6XBGTT8IdL4KjPNXeZYwhSXZRWeRIWjV8jInsD31PV9wBvx0SUfUpEnrDH/bN12VqSxLhxjHvrW+m96y4AovkDO0vrZqLGVqioFcMJCm5plvuNkXBndzPG/aNe23Ixdbepr97F4cqv27KlYHvPqVP7lV1jl1lGaptVAvpN4Pv2XG56zl9X5sr5J4gibsbcZTnPdKXKlqtjnA2aI5VqVtyNurV5be8gM2Y83d2wrnPf2CWOe6ftMkNPxUql2gvL2zZqV1UVPcbR7l+ku/Di/X57xylptmClStomCk18JqVf71e/TandC/cBqOZmwytVznwleaCKVzkbs3raX6k2X0FLJMoraEOfuj7TwAgqbvmgjx8Q2V/q6JfLFgmg7Lw1Opszp5yNavt7AHrSOY/I+UQKT7TksXD5ZSJmKWQcua/ywnspFlQ6XpUr/R0f9/74lKtxpecpde5qn0ZZZMBenAeC1ikI9ojAOelwdl/OSUej1Ku1i2H+WfCXb8PuB8NrT8Okfc1ADZq31LHLmlkEBa2AlhigqepGisSvsfEw3mP/f4AQcigwTAhtPjCiGR4qWSyhfweGI6qgQSkbPLOPMYOxG8+APQ6FDS801knHUZ+FB/7X/L/5JUikYNMK2LqmcLDY6IFapKCFAVo+w/sp6RCBRCL2Zt2MX7dV1Do9tSlllaasUxJcfl4ZP+yZo7PTnL2fshCd3CoFMXVyxDk08cv7+e6qW9av73dpd80NVm1zql6my53Us63xVSlfCSinlPnlyu2vJIiMu5av5vnHNiJmVYtgm3ssrq2+zEwAZhY3j8lRwx/TZU/lFcjGKGflRCCnlhWYMtn/X7Nq2x5sKjyonBJWaaR2H7/9VmNjVqk9XKlrNNu+a2QoaPVFjI1S0vOI6NSq/p4PS3/nFlOjfJzNmVPOOlImpud1B32yoA6RjVmkypnjc/ZxFJQrVyd3XFa87SLHJLxTuvchQWFdcl4ai3t1jPZLYblKvTj6Nm7+eTXv3hPeZ+U/y31bNNfFUzFdvC0R58uyFgiaGTnP0row+xgzOHv5L3D05xu/tPDVp0x62s8g3QM3fBgevQ4e/wl81AaFbaSqB2GAFsPIGKAFAoFAoDWoxItjIBBoPRTIhgFaVaxcCK8sMf8/9n3Y9x2NHaRtXgWjdjPXBTjsI7DkRkh3wyPXwUt/brzr/WCDVpQR8ZTcZ84crrmliIMs9yPB2d64fH9218pi7dZH+IJHHzXpwyZkzbm33cZl/3EO11x9NQsfsN6wOsyhObM2kz9v7tyidZw0ED/h9eK88wA4e4vxqHPNjb9sZm1qSynJZbQfWGFok07D5Zdfw113wZ/+VN4FfXe3UZSckOrSZctMevPNTo36RXTMpz51LwDXXXdNdE0w3/H5jLZ9ybnu77LqbKZzX3tts71jg02XUpDvm6a5fF8E8u+xEjOrdHoSz6/q4H6rJJ7T0cE1F1yQK+DbpJZRwQMVUokXx0BJ9thtHf/vtH+FsR2FO9q9tRi+ktptPPXSZTpQZqNp4+sePhSAZ1fNBmD8QbO5eePHWPfUambNXwX0V7yccvZ3S38I5MVJ6zGdWdV2wj4bUsbGPkz07jTn67XBumw5ydpObFNJ53mbK9ifLdjWRLxNVKLd1Pm+8w8D4PQr/lxYQAu/KKJrOOx2orfLyy98X/265vK9LyrvHgCkt6fwGLftvsRcWfWCPvbac6czhWkqyT3/6n0R1wwhGxS0weOWER51Efzxq/Cufyu0SWsE46eB7YOsXAjL7oI3nwOPXgvP3QqHnZ6rS6OWOu7aBG2joc2PVjyyCT2thpxrBzhDFRezCuCaSy9tXkUCLc510X/nnXduE+tRPR//+DXR/wWDs0D9cEscS70CLcHeb5zR7CpUhe5aGv2/6snVTazJMMG62Y97BcqwdrEZjM2xETrG7W221y5uXB26txiPX/k2Zwe/D1J2cLTkRnj21tz+afPqX6cQpLooI3vqt9w0e1x03BK4QLzukDFjBlGvZpPNNt/TW62Ji9I8DEkmjTnhuHGly7m3wDnA9N+SCRNM6px8rl/vZrcuJJX8nfk3myWR7o1sLNv9bxTbIcZb+0YXB82/tv/xuHynjLnyfhhBv87ltvPJqXD2Yjt3Fi/oVy6usoHKCEsca4MqjLLqvzNM8tM4D6Rus8N0kLY20xmcDVtSsv0CT7d5HiOdolYuTlomYTqttpkvpKwYe2hJGfVP0kYx0khdsopaylOWPIWtEhK9CWhrRxWyCtK3q/g5Hd622y+93SWv3U8pi/KdCuYpZyUUNPqsGucraGnv2r12f8Z5Krb5mbS52TqgVGavGIjBKVFb15h0x6sw/xONXU7YtRnGT88NFsEMxD56E7xwJ/zlKrjp49A+Fk6/oTGBrHdtDMsbixCekoFAIBBoHGGJYyAwNFEJSlktGGPDwmx/tfHX7tpinJS4wdYDV+SWWM4+BtY8AqsfgWQbzDq6MaEAdm0KCloRRsYAzS2p8V0slsPN8vqzvSWMevyiPT3F8wNNJv8DaU4ctLohYuKHjR1bqFbF4d6KHs/0YZQVzHK/pYv8qE6njUrm9xWXuj5nt51NmosX6OPCA/qHO7XLMxstqEYxSt17bp/9weHkb0elSlkV8RNHJEFBqx4F0kQ2Z5EdVhSA08Wk9BQh94+z4+q0ylnKdIakdRWYyHOH6HuKdD3XbTubMz9OWm/m+6ZcwtQpkzb9K9vWZets8jVpvhMk46lU2UL7Ol9hK0Y/u7REEmnP2Rgnuour5LHqXKSg9dhyMa4TnYLmK5fufXb7iyhokWLmcHaCTgVzypl/6T7rWdK9HfkhI+saSnF4PS+bQqodRk+p3QDtgSvMUsR8JS5O9eraYpyEOPL3r1wIG5ebGGmvPwvXHQdbXqq/jdyujTBpdv3OP0QJNmiBQCAQaBzBBi0QGJpYBa2VbNBE5G0i8m0RWSIi60XkZRG5XUTOF5EJZY79gYi8LiJP5+VdKiJrReQJ+3pP3r6LRWSZiDwvIidWVfFxe8KO16o6RYQLgL1yodmOsx/L9EHvdmOD5pOvlH36QRi7B6xbDPu+s/5LMIOCVpSRNY3ppswrnb31DV18w5kinHjiSvufm4lzZe06/MwelV27mSQSYYZ7KKMK6TSpVKFHtzj1yeU7z6NOtXK/k3NNwY/WBzv6Onh5x6RYpWrfzuJ9yNmspa1NmhOhfAEu7re6b6MWR34zjrvvSFcYrF2iq0SwRauMGihoIrIK2I75ok2r6gIRmQT8ApgFrAI+rKqbq7pQi9PPg2GcR8OYfPVUmShulwppLyCxU06yXpwz563R2Zw55aw9eTaQi5umCVsuYTu750FR7TNS6LbbhVSi2/jHaCKJ5Nnd9bM5cwzCvm1QVBKnMCoTU5cmTqu3mg2aiPweWAf8BrgMeB3zY+sA4J3Ab0TkG6p6a8wpfgRcBfzYy/9fVf2ad62DgdOAQ4C9gbtF5ABVHVyjGbcnbH9lUIf2wwXAvunjMG0BrF0EH76+/8Cqe6tJ8xU0h7NJczZnXVtM/nPWYUi9bNEyfdCzFUYFGzSfoKAFAoEB8vVmVyAwlHE2aNUraO9U1bmqusBufxG4R1X3B+6x24FAoFaokM0kY19N4ExVPVtVb1XVdaqaVtUdqrpYVb+uqscCD8YdrKoLgU0VXutk4EZV7VHVlcAy4C2DrvnYPWF7jRQ0MAOomUcat/kH/XVx1avLzleNmth/31GfzQ3CfvkJOO5fzcROx3i46Sx48Kr6eHXcZd/+4CSkH0EmySd+in0QDOHZ9OHoxXGg6ulQJpuF7m46O/srXvn4scSczZmzRXO/kzsiUxBna5Wb12lrM3Zj7lz9FC3PBi3Cs0lz3h39OpXrkm6/f1wx4j/6ZGFdo+xB2pyNhDZWDfWzQTsZONb+fz1wH/CFelyo6ShoL/1szrSMDRrt1nNir7F7EvcxONszUcYdNBtQeqzIM9F6b3S93nl6jEx3bZwz563R2Zw55awjdU7BdiZlDE2dDVpkU5bpKp5vlTbF1t3lF1WYvLhw6RSayout1C+2mqdkxdiSRZ4WY5QvKWOD1u+LLJO33et9ebltZ4PmvDS6Q+wlIgEynVM9zT1pXW3QyLaOgqaqGwBEZAzQpapZETkAOAj4var2uTID5AIR+TjwGPAPVomfBjyUV2aNzeuHiJwDnAMwc+bM4lcYt4dZ4pjN9vOuOihWLoSV95n/n7kFDvtw/0GaU8WKKWiOfCWtZxss/G+YehD88TI44xe1V9KiINVhiaNPUNCq5Nzbbmt2FQKBQGDoUJkN2hQReSzvdY53FgX+ICKL8vbtoaqvANh098bdVCAwAlDIZhKxryayEOgUkWkY9fyTmOWLg+G7wBxgLvAKuSUjxUamRYfCqnqtqi5Q1QVTp04tfpVxe4FmYNdgxo8eTvU6/ONme993FtqkOSIFrcQAzSlpAO/4gll6uH4pTD2wUGFLpIxzkmoJA7RYRtZUb7lZ27j9cd4cU6kBeAAcQoracLRBG273UwpVq6CNB/qLh3Fioq+Yuf2joglo14ZzbT6ZNMfF2bexYUfxi3l9ytmkjRqVKFUsOr8vzFXy8fp1y6lu9prbtxcWGDWqcNtdJM7wzSlulch5I5nK3OxvyFu6WIy3q+o6EdkduEtElpYoO/xQ0J5kpIjlFLS2gu3+XgkLvTfSYdp+MmW9FYpGvdsJN22eYpaIvDnaAn2mj7s4Z85bo7M585W0nd3fsFWxXhytItZPOcvaOvq2atH+Iqs8XFnJfSFosh0RSEiRmGOWSEmLe7/KKGhRfpyClvHqmv9l5Mc3S9v3Va2XRif65XtpBLQvUZhvlS2toxdHRZq1lLEcoqq7RORs4Fuq+t8i8vhgTqSq0bpDEbkOsEE/WQPkR26fjrF/GxxjrT+C7a/C2Crnkpzq9dozZnvXplwA7HwVrXuLSYs5CSnGyw8CCm2jYN3jcM07YOtqOOpz8MA3auN6v8stcQwDNJ+goAUCgUCgcbgljqVeZVDVdTZ9Hfg1xhbkNRHZy1xC9sI4DAgEArVCzVLKuFcTERF5G/BRwC1rGtSsrPsOsZwCOA+PtwKniUiHiMwG9gceGWR9jYIGtXG171SvjJ3UeO1pE8Osn4v9ChQ0h1PKPvxj+PuHIdUJrzxhJkb+9PVChyLVKGlBQYtlBMkKjcL33hjGwC1J/gym1nPBfhNIp2HLFiZO3D3aLIUTM3xhaOLEwhTG2DRng5VKKqM7c7PKLr5Z9Bt7qRdbzA905ipnKzHe7u/sNOdxSplva+YLMNXFQbN12enFSOrxZtvLDRyCF8fKcEscB4m1N0mo6nb7/wnAlzE/oM4CLrfpb2pQ29ZEQNqyZEaPL8jW9tLvq3saaaSkGTvVVLvpWO3JNG1J0zF67fdBW8IpaDk7NbNtz9lr+k1W1pvUxjlz3hqdzZlTzsZ0fg7Is0nLWq+NTknzvELmbq5wf1E8VU1TnWhyAknJkhQhsWtHzHFF4pPlb7tYZf7+cl4Zoxhm/dxl5v7v9b54+7znURTnzL7hKbvfV86ypi7SL2BabWnyUsY4PgtcDPxaVZ8RkX2BP5Y7SER+jrFbnSIia4AvAceKyFyMDrkKOBfAnvcm4FnMp3L+oD04grFBA9hRw2DVTnXu2Wbil+02q3C/s0HrLBl9wOB7dWwbDZI0niBTnbBzY6Fr/sHapUUDtOAkxCcM0AKBQCDQOKp3ErIH8GsxI4QUcIOq3iEijwI32WVOLwOnVl3XQCAQoQiqrTdAU9X7gfvztlcAF1Zw3OlFsr9fovxlGHf+1RMtcayhJ8dM3vL6V58uMkDbbLwyJiv4/nUDLTcIO/ofzLLGiTOMTdrNnzBLiY+/1JQb7EBt1yZoHwupjvJlRxgjc4BWzhinHM6LWzIZTSNee91Gr5BbOO486dU5vkqgNMPNK2UpnA3aFLMZZxLlFDMnZrjUE7XyvDi6tpyALVvMv+PHmwMiWzJv9jbOWCyuz9mLOps0rJIWZ0fn7s0Xr9y9lerauSaRKKyrI85uLm57JLWxaqjMBi0W++PrTUXyNwLvqqJmQwcFzUg/xSzbPirmAINYL4ZqvQ5Ku+nTiTbTdlPJDMlEYR/2lbNkpKDZtHeXSe0PrMiWzCledtvZnPk2aS5umiZMp3WKmlg7Ml8x66es5RNJhM7zI5DoALrswTEKmcP1YaeMOdsxl5/1V1uUea5XoqBlCs8ZxaazqTrFym67hYT9AkPbcgpQr+WGWr2CJiInAVdilhl9T1UvL1LmWOAKzHKNDar6jphz/ZYSFneq+v6qKlsvHrjCuKtvGw3rFpntRApW3Af7HmuU4GnzzEDHpTEDnqvvX85h0ydw5JwpsOrPgKDAYw8v5M1veG/hgKl7S+X2Z461iwttzmYfA987AdY8bPrZ3f9mlLXjL4WnfwVP3wKn/Sx33TL1Z9fGoJ7FMDIHaIFAIBBoDvVzsx8IBOpMNbZmIpIEvg38FcbpxqMicquqPptXZiLwHeAkVX3ZOgKK42sl9rUu0+YZxaljHCy7G8btDY/9AA44Cf5wCRzwbnjwm2Zg9MtP2AHSFUUHOYdNn8AFNzzOVWcczpFjp6IoOxjNgTsfK7IEcVHxGGilOOqz5tr5yx03LYM3nmoGZJo1r7v+FRJtZoL4mf+DJ38Obzod7v8vk/70QzGDz0XB/iyGkfmUrNStXRx53tyued/7OPfWYkHq3ezW0FHOsiRQBBlpXhwr9sQ5RFCFdHrAH2GcyJUL+5XnxdF5O7PX6ncSR9qbgXYqlS/XxeyPlLQYW85Y75EDQgZ3kpEUW6+WVGmDFgBU0HQSTRkFTJ0Xx6Tz6uh5RLQ4hS2K29XmKWiJ/gpa+bqk7Tl7bB28z9bZltnUKWROOWtPnl2wnbTx1LJSaHPmjk+6OGpFbNFyebk+KfkKWlwf970tumV8URyzMs9xp4glvGdJdHwJO2enlLm33a9iutDGjIRT5QqvpXZb6mj3rirVKmhvAZZZFRwRuRETv/DZvDJnALeo6svmmhrr7McubcSeaxQwU1Wfr6aCDWH2MWbA85MPmN8fj30fdn8DvPB7mDDDpPu83ahWZTwmHjlnCledcTjn/WQRP9yrjYPpYIz0kli/GH5+Bpx+gyn4y08Yb5EDHaBB/+WOri4v3GliC2b7TMiATA8k28397DbLpAe8O2/w+S8wbYFxMjLzbSYdu7t5D5yqWKv4asOAllhMLCKTROQuEXnRprEuZkQkKSKPi8jv4soEAq1OaPOBEU2VXhxbndC/A8MV1UTsi/LxC6cBq/O2iwV7PgDYTUTus3EOP16uTiLyPuAJ4A67PVdEis2ctw6zj4FDToGMNYd5/TmTbrVvz0t/hnQPLPyfsh4TR7en2Nad5unVGyHVSeKAEwGF3u3wl6vgho+YgV42ayZXH7hicN4XneMQMAO1034G7/pXY9OW7DRLHTO9ZsC1eZUp98LvIdlm0tGTYO2j0LcLVt5vJnh3bTT2comUOee0eQOr0zCmVZ6EXwTuUdXLReSLdvsLMWUvAp4Dxsfsj6dcnLNKcbNwzvBFhNwy6LiZtqGjpJHNBnua+lO/Nl+hgub2d3UVbrtm3d9cLM8zqfN8Wa6tOCOxOKXMr4xfzuIraeWUs2L3Hl/NChW0OMXMP24YDDDqyshQ0Or/TMtWofxnfduowXsgFXsutalYGzHFuYc1XzDOFi0X98yUj1PS3Oyxq6EkPZu0IkKfb7emkgJJIu58cd4WnXxVzkatnw2aX86b844rn+c5WMsJlgnPa2O/fKucuW2po1diFTKlFbRy8QsrCfacAuZj7ElHAX8RkYdU9YUS570Uo87dB6CqT4jIrFIVbTorF8Lye+Cwj8CSX8D0N8OaR2HGEbD6YePMo2ebGbhkevurV3n85om1AMyZ1EHvtgypVX+h7eBT4NlfG5VrwnRY+N/GicjmVXDge2LPVRKnbLklj2DUvY/eDK8sgXu/DNkkZHth4j7Gk+SoSSbW2ajdch4bs2nzaktA13aYckChjVsAaBEFDSNxX2//vx74m2KFRGQ68NfA9xpTrUCgboQ2HxiZ1CAO2hAg9O/AsEOpOg5aJcGe1wB3qOpOVd0ALKSIUyCPtKpurfA2mo8bbB31OWODtuBsWPOYWQ64+hGT9myHqW8wA5mffgh+cWbRAcyDyzdww8MvAzCzcyfj2cX5fRfy9LQPGu+IJGDrGuMev28XTN6v+sGQi7uWr6g98A047t+grdMsY9zyMuw11wzO9pprXPzve6wpO+ddxl6tzzgYYsML5j0Ig7MCWuVJuIeqvgKgqq+UMAq9AvgnYFy5E1pp/RyAmTNnmsxqbUZKHudiRA2DWEjD0Qat9eyFatrmC9r7HntAd3fFt+pszJyo4cQtX0krmM+ZbI16e3tLK08ultgo613OL+srZjEKmsMpadmazi0N0gYtjnLnaZ022BxGhoJW12fajIlS3FNfIlV6O8q3nTsR348q1uciNS5dsC1YJc2zh3PqlrNFczZnvpKWzprxrai1r7Mhp3xFLR9fXZNECpEk6rLi5KpycdCyWnx/VM6JQP7xUcWK5xc5pJ8y6raj1MtXZ6PmlLTiVawV2WxVF3gU2N8Gel4LnIaxOcvnN8BVYuTQduAI4H/LnPdpETkDSIrI/hgX+w9WU9G64gY2+ekJ/2G8OC74lHGwccJ/mD7z/O1GUZs4s+gAZsmarcycNJoXX98B3VuRcXtw4Vv3Y859F8AZPzeFfvcPsPEFmDATXn8Gjvmn2gyGfEVt7WJ4xxfMYM3dxwHvNssbD3g3vHCHGYg99Uuz7BExatuhHzL2arOPDoO0PBqmoInI3SLydJHXyRUe/17gdVVdVEl5Vb1WVReo6oKpU6dWVfdAYDA0ss0XtPdcZOlAoPVwbvZLvYYATX2mjamfY6NdLyyr27kDQxtVIZNJxr7KH69p4ALgTsyy3ptsAOjzROQ8W+Y5jC3ZEuARjCv+p8uc+jPAIUAPcAOwDbN0uDVxClR+euQF8LGbzUDsjF+Y7Wnz4NWnTGDpV5+CFdYnSp792DlH78vaLdZOYcwU6BjHobqcUWf8JDfY6dpo1KutL5sllY9935yjHveTTZvBmruPmW+FEy4z+Sf8R25CZ9bRRm3b/0RY8cecx8pa1muI07CpXFU9Pm6fiLwmInvZmca9gGJee94OvF9E3gN0AuNF5Keq+rEBV2awXhz9ci7YUsHsWNyX1DBQ1oYyTVAtmt3mB9qc47w3FrVBK0c5ZcyvXJx3R4dXPmGnkVOp0nUp5WAyhxTWweEPFAZj+FauUiORYeJmv5n9W1XI9qYi741Y2yunjDnVSjwvjpHXx1Rh55ZkLg4aQDKRs5mOlq55tk0u38VWk5RVy7M2HporGClrxh7VxTFz9mLOW6PryU45SyXOKtiOsfIqIBKjEq7KRkHLqhW5yvW9SCnzbMajeGalD4/u2pXLxNiPFZzb88aY8bZdvDOnlNlt9W3Bsi4OWrZEZLDqqcbNvjlebwdu9/Ku9rb/B/ifAZx2lqpeAlziMmwstfsGW8+m4XtMfOclcN9XAYVffCynUNmlhas372JXr22fmbTxouifw3mCPOGynGdIZ4NWa8XKXdud16VHXmDSB64wzkXWLjZ5zvlJvpoYVDSgdWzQbgXOsv+fhZG4C1DVi1V1uqrOwsji9w5qcBYItAZDuM33iysaCFSOW+I4xBW0Mgzh/h0IxKCQySZiX03kJhH5JzGMEpFvAf/ZzApVjRuwHHkBfOgHJm/MFPjjZQUDq+de2Z47JttXuKzZncOpWkdeULi9dnFDbqUAXz2E/mpiAGgdG7TLMR3sbOBl4FQAEdkbI2+/pyZXqdaLoz+DXuDF0ffSaF2nDsF4aIGGUL82b704VqoyxRHvxTGPSj1+Om+OcfHP4jwjxqlWNt+/x8E5VIx5fwZqS1ZPZWywnmZbValq1XrVjvo+0xSjqDivhZ5yFhXztstpH5E6kqeSZO1R4hQzuytRqUzj7MK8uGWRp0WbRrZl1uasnJJWjH6eH60Xx/4FvThn/fY7e68qpai4r978y3reGMshLk6dvQetxpvnAFEqdgbSaI4A/gtjdzYO+BlGoR665A9UDjgJ9jgUXnsajryoQGFa+uq26H/J9EFnW/9z5CtSs4/pr24FWpKWeEqq6kaMS1U/fx3Q70GmqvcxFKXrQMAS2nxgxDICnISE/h0YrjRZKYujDxONfBRmufBK1bIBDIYOKxfCFhsf7dHrYP/jufrlaRw2fQJLX9nOjEmjWL2pi97eHtZlhb2bW9tAjWiJAVrDKGdDUu64AbDbbgcCsHlzMLgONJdEzJyzrzrF2aKVxCpor6lxUrcHmwr3O+XMUakXx0r7aoySVoyy3b7S74XBfo/4xw+EwV6rFe3dhokNWlMRSse7crZpvi2VFH/fXQwtcecscm6nnKi1uc5aKc3FQesXD80dR4d3Ik9JszZquThnhd4a45S0Yvh2appIIaSM4Jiv+MUpZ5Xif80Mn6FAWVpUQXsUs4z4zcBk4BoR+ZCqfqi51aoBzo7s0FPgpb+Y4NM3ncVRR32Tb/5sGc9k9+FN++/D6k1dbNmxi8l7TW52jQM1IjwlA4FAINA4nBfHQCAwpFCVVlXQzlbVx+z/rwIni8iZzaxQzciPNfbUzdC7A46/lEM33s1/p+5g7ravM2OHMamZMkqYOHZ08+oaqCkja4BWzi7DeXGLsy0ZlF2HW//ePoBjAoEqsDZo5Ri4jVoVnkgzdbLB9PpknFpoisTdZ4UzwoMzcCt/nlIMR6UpKGjVo4JmkkjGeE7UhPejzFfO/MMjmzWTJtqMm+5UqtCbI0DWKiYJq6pF25HHQhcrzIuHZnGeJNWLg+a8OSZTYwvz7XG+PdlAlDQfdV4cs4X2W7Fxz/D2+/HM/Lc3Oo0UlnPmdyUUtv7eGL3vIy/OWeStMWZbSFDxd9ogaCUFTUTGq+o2YIWITPJ239aMOtWcfFu0D/8YfvpBePyn0LWZF4/9HtzaxUMrN9GWFDokY7w4BoYF4SkZCAQCgcYxAmzQAoHhiNJyNmg3AO8FFmGqlz96VGDfZlSqbuz3LtjjEOMs5Oh/5Lfrdwde4uy3z+IHf15FJu15cQwMaUbWJ1lu1rbS/SVnwAu/vKZOnQ3A+vXFwuAEAq3PQMSO1/rMJOYebdYWbbDKWa3Uqmoop5g3wlPiQD1J+se1olJVgyWOInIScCVG1v2eqo7o2A/iqVe+98YIq1JJGYWt4BBnuhWJQ+afyFNiuXPFKGu587s4afZ8ydJtdiBKmkgKifM+Ww3iqUiJ0t4exf4sKKakObu/VlKmYmmxJY6q+l6bzm52XRrCyoWwaQUAfQ9fy90792XiqCn86/sO4VeL19Lb2836rixTm1zNQG1onZ4WCAQCgeGPW+JY6lXycEkC3wbeDRwMnC4iBzeg5oHAiMY4W5HYVzMRkQ+IyDdE5Osi8jdNrUw9cM5Cjv93AB6fdgZjdSeHTzWj/jEdKcamlE3dTaxjoKa04PRqC+PPSrfZeBOJ/HFutqBIIDBUcGH9BuTF0eK6QuS0sZy3xmqVn0EoRM4+LVtp3LNKz91MtaqcwtaKXhyrX+L4FmCZqq4AEJEbgZOBZ2tQuyGCIolspZHIYnHqVzZjOn82a1OtOMpZZMfWaHwlLT8vEqo8j5GFz+oi+EqYKx8ZxEnhNnH5OR+W+RUqNoRR50Ez67a9AlZhcypcvzho9uRiLy6pTGkPn1XSSgqaQ0S+A+wH/NxmnScif6Wq5zexWrXFOQuZtgDu+AKH7z2GZUv35qQOszqrLSm0pTMcuNduza1noGaU/UUhIl9V1X8WkWOAJaq6pf7VCgQCgcCwpHonIdOA1XnbazCBagOBQF1p2UDV7wAOVRsDQkSuB55qbpVqTL6zkKkHseulRWR1Hoe89QQA2pIJkpqGZFvx4wNDjkqmQm6y6d8B94vIiyJyq4h8RUQ+JCIz6li/QCDQcny+2RUIDGFUoTedKPkCpojIY3mvc/JOUVSIaEztA4GRi6pR0OJeTeR5YGbe9gxgSbmDROQHIvK6iDydlzdJRO6yv3XvEpHd8vZdLCLLROR5ETmxpndQIVffv5zXxx5I++tLAOXQaRN4cPkGtnb1kdB08OI4jCg7jamqT9j0TAARacOs+38T8DbgYhG5U1X/uY71bAwDXdLk1nMVGBybeBTOY//27c/b/OC1rKVoZScK1aJq2uYg79EvPpDDI+cglQaF9tNy5eMqNYB7TfRz5F2GcsGzmxG4eghTYRSIDaq6IGbfGswPMMd0YF0NqjZkEAFJZqHXuMcnZYNE2+WGgnG/73t0k3SvTe2zy34Q2mfKZazL92wm10uydjyctmvvnLv9OL8YUcDqKMcGqo6cgdjUOc5wjkskVbAd527fJ99JiO84JGPfh4RA0nfsUYqEc2Gvhdv425TOj27C7S92A4Vl+tXSfbZuAJQqDG8g7r7cNZL1jZqdbbKtWQyTgedE5BG7/WbgLyJyK4Cqvj/muB8BVwE/zsv7InCPql4uIl+021+wdq6nAYcAewN3i8gB6iKrN4iTtt7IvSu7OE03csCo7azasIPrb/gpkxPTSGk6eHEcRgz4k1TVPuBJ+wJARBYBQ3+AFggEAoG6ks3mJrAGyaPA/iIyG1iL+dF0Rg2qFggESqC0lhfHPP5tMAep6kIRmeVlnwwca/+/HrgP+ILNv1FVe4CVIrIMYw/7l8Fce6Bcff9yDps+gSPfeDTTlnwX0nBoYhXXXn893xv9bc4Y9V2SW8ISx+FErYbaR9foPPXB/0UwUHf6fuqUM+dCXPOnEgsNpnfb7UAANm9eM8BKN54EWaQVVwoNRpmI+wyLqR3agvdcS+LeP/te+D5v4nFtO28G1ZNDIjf70vphJXJfA1XOCFerxtbCkUc5RbHFFONqbllV0yJyAXAnplH+QFWfqVHVhgSqQrYv1a/lRuqV77jDbXsBqn2nGc6+SIs4p09JZcpMrNMQ8dpe5MCjUDnzpbNyend+rXzHIb2Z7wNtZBUyA/mez2rx7ay33U/e88qpt12s0t4+38YrUs6cI5K0FG57gaylzs65W8kGTUREDfeXKjPA0+6hqq8AqOorIrK7zZ8GPJRXbo3NK3bNc4BzAGbOzK28dIOsJWu2RmkyAX9etpG37zeZPy/bSEJgj/GFK65e29bNHuM7+fYfl/Hew/YiuccVfHnN3/KO3vt5V/uztJ32U9ruaiO1JSxxrAVX37+clzbu5H1v2jv6rH775Dpe29ZNVok+q3Kf2VtmT+73eZ/3jjkV16MmT2tV3VWL8wQCgUBgeFMDBQ1VvR24vSYVCgQCFaEKmRYaoAF/FJFfAb9R1Zddpoi0A0cBZwF/xCxlrJaKbV9V9VrgWoAFCxZEZQ6bPoELbnicTx+7Lxfc8DjvPHAKv1q8jiP3ncRlty3lHftP5t4XN9KeFLOMWYRMJks6C2ccMQNV5eZFa+jLjOIzHRN4f+IvXJf9EIdmD6Yz+YK5SCJ+pnWgA49yA5CBDFbizp/JwksbdxaU32fymH4D10bW6aWNO/m/x9dyy+K1fGj+NP73rhdQVfoyynEHTeWrty0tSO9Zup6OVIKE/czSmSx9GWWvCZ1cu3BF9HlfdcbhsZ9NMVpjOrXeVGj0EJVxLqDj7GOccpZOc+5ttxXsamvbD+jvqnz27Om2RH3Xh9eEbLa27rkHeq6BuDuPO3e5/Py02qClrYaqaaMVurRPRPsTBcVdmowmxN3MbO451ds+lm0TZ/b7ItmEmXScFKc+xyk+ldqkVYO9dsLviwO9Ztw9NINK+0ELuN2v9Os4UAKFbF+Ktl4z0o1aslOv1KZOUWsfXfw8VkGTZJ/dzDJq//2A8i78nS2Ss2tL9DqX9J5Kly7sJ5rqdLdQmG8VNrE2NKI2dbZpvj2dU+TyXOk7mzOjnEF78mzmzZvHCWNPNWJWd6+to2dj5khniu/vi1HCXN2dG31PvIqqVuIRoz2FiqOmvW1rH+gUsoQLFN7r7t8paPazTGSjsvUg21oDtJOATwE/t0uet2AM/pPAH4D/dX4UBsBrIrKXVc/2AtxykKptX4+cM4Wrzjicc368iM62BL9abA5/cIWx3b7/xY0A9GZcuzTpzEmj+dnDqxnTkaQvo+zGdibLdmS3ffnkrns5/2c/oWvKYeaQEoHeD5s+gW//cRm/W/IKF75rP879ySJ6+jKxA4+4AUhfOkNfFo47cEpB2WMPnMK9SzfQnpTCwUoWjtl/MpfdtjRK32WPOeOIGfxuyStksoqqkkwI73/TXvz8kTW844ApBcccd6DZPvaAKdz7fO462IFsX5ZosOvSt+27G39ZsZk2VyeEdNYMet+8z0Quu20pC7z0hIPN75dMVvnZw6tJiPkqOGjPcdyzdD2H7DWee5au52Cbzpk6huXrd5JMCKpZsgrTJo7iZw+vZu6MifzvXS/y6WPncOScKQNpLiNkgFYHzv3Od2D9+n75mcy5AOy0EwJuBUluxvi79a9cIFAX/l/R3IsuOjf6319F1LHqea75+MfrWanAEEO1egUtUHsS097IqC379cufdPA+ALR5SxxzgYlHUPi5EY4iZPoFamseqtoNfAf4jnVgNwXoqjIc1K0Y5e1ym/4mL/8GEfkGxknI/sAjRc9QgiPnTGHujAk8sGwje0/oZN3WbmZPHs3KjbuYsdsoVm/uYtbk0azaaBamjW5P8vKmXQiwsyfDGxIvc1Pnf7I0O5tZiVGMPe1Kvn3jWZy1/TJzgRJLHI+cM4Vvnj6Xs3/0GJfdtjTKTyWEe5auj9K2ZP/tnnT/WYZ7n98AwD1LzW/hP9ptf4AJsNAOPl16z9L1TBzdxs8eXs1+u49l5YadZLOKAjc8soakCPe9sKHgGHc9l1/sOm6w69K/rNgMQF+Rso++tAWAx7z0D88Wmmi4eZqlr24H4JlXtgHwrE2Xrzc/+DN5Ez5rtxgHTk+sNudc+MJ6PnPc/gyEETFAW/HSGk475x8rNr9wPzLdymV3nNtub1OYOdO8LGcfbnYmb/hJ4UmcOrPLzMx9+LTCbIcbyLn8uHialZYbDNksXHrptRx8UJZzzzuPrdukLtfyV4SXMw0YyArySs0M3Dkvu+waetLNCbRaL15dvpwjTj6ZcRWWdx9r0kvdgoBtNs2Fg4Vrr3uvV7r/lx/APO5n/v3397Nec6mfH1cuLvWppKn0879mFXFJHujt8R9I5RpX3P4Yj29NY15Tr16LJY4jnXU79uArf76IxJ8HdlycPdc2K2luI/fBpDA/MFavabPHuqOtkmPT+/7xw0Dl3+PJmE7q+UkkGznHc14he8pexxe+vnjMf7LyydVcvvVi5s2bx79ee2nROg0U/7haakpxDiQz3n2798E3c3t5yWpWbf9FDWtUSIstcYywDuxeGcgxIvJzjEOQKSKyBvgSZmB2k4icDbwMnGrP/4yI3ISZkUgD5w/Gg+ODyzfwzLptnHL43vzf4+s4ar8p/HnZhoL0gWUbaE8lUFV29WZ4276T+cuKjbSnEvxV6ileete32X3ZzWxd+SBLsgdz5GnXM+s3S6GLkkscAfaZNKZAOJ46tp31O3rZfVwHr2/vid3eY1wHr203fdANLKdN7GTtlm72mtDJK3nbbqAJMHO3Uby8uYuZk0bx8qauaN+Use1s2GEU7WWv72DquA6mjm3n2Ve2R+ebPrGTNVu6o2PddV3+PpNG8dKmLntfo3lp0y5mTxnNyg27okHvvlNGs2LDLvabOoZldiDl/t9/9zG8+PpODth9LC+8viOX7jGWF17bQTJhvL9mVMlkYe6MiTyxeguHTZvAkrVbedP0CTy5ZitvmTWJR1Ztoi0pqCrpLLx9zmT+vHwjJx26Jw8t38jZR88eaFOpszVpIBAItDDnNbsCIxC3xLHUKxCoBSufXF2+UKBiFLPEMe411FDV01V1L1VtU9Xpqvp9Vd2oqu9S1f1tuimv/GWqOkdVD1TV3w/0eg8u38AFNzzO379zDve/sIEzjpjBn5dt4LiDpkbpA8s20JFKkBQTfLo9KTy0YiNzZ0ygI5Vg3Lv+gY/f20m6cxJ7pnawZM1WmH0Mnbtb5xMlljgCPLTSqFFtSaEtKWzY0ctR+01h/fYejtpvSuz2a9t7aE8laEsKr2zt5qj9prBui0lf9bZXb+6Kyq7e3GXyNpl0jd3eaM8LcPLcvenuy7B6cxenHL53dL617nz22Fe8/Jc25a7z8qZdHLXfFFZtsOlGk66028vW74zKLl+/0+S9btIXX99RkL7w2g5SCSGbhVMXTGd0e4r2pPDk6i2866CpPLV2K+86aCpL1pj0kVWb6EglSCWE9lSS9qTw4PKNfPSIGTyychPnHzeHf77laR5cvmFA7WVEKGj77rsPN97w3fJP/jK2OuzYAcC5bsnWq69GRZLzzKz02uPOBPp7xHPq069ubN0ljlkSnHfeuci2rVzz5S/Dfv2Xuww/snR0NLsOtaVt1nzO//fHuPtus12pOd6WLYXpsmUm3b59MXAZsCw69qMf3RuA1z1njU5tdbO67e3n8JnPXMN2szIgWvrrrumUFNu1otTlxzlO9RWYgZoiun0/Bf7LTuk7Jc1hAyYP6JyV0AoDkAkTmvdjKixxrJ5pk9bxHx/5F7SvsI2KHwPLs0Vy5bNd5kuvd7PR2dfuOBmA5zfuHpVdus1o6BPbTGdOSqGE41Sk068wMp70mSVZiW7TyZ1HSVzstV4bk2yX7eTuy8J1iOjLw8v3iVt+UozuXpg0l9mX55Q0KIydBpDVwms5YcTlq78/W5gftx9/v8vP9uUKZwo7g/RsNWm6pzDfbkvGpjYGnotp52wBSSQ58mOFThdqRus5CRlSLFmzlavOOLwg/ee/Pog/L9sYpccdNLWsR8Crzjic1Q/8kb3TuzjvSPMs7nR9v8QSxweXb+A/fvccAB+cN43bnnqVHjLR4PDepeuLpgU2aIkEQjb2mIGUvXfpej56xAxufTJng3b3c69zxhEzuOHh1U2r09wZE3j2le185M3TuWXxOj53wv6sWL+T17Z1k1Eq/sxmTBrT7/MeiB1aSwzQRGQS8AtgFrAK+LCqbi5SbiLwPeBQzGTOp1S1ITEo8rnmc58z/zzwAOfeemvFx11yyTUMCSchgboz1No8XALcYP+/q+KjPvOZa+pSm8DQZSQ4CRlq/XvXi2byZeOrXUy2NmeV8OrTQaEaSSit5WbfISL/papfKJfXbJyLdfcj3aV/d/ScgrQU0Q/8LQfASmDXJpgwjU5xDm7ilzguWbOVI2ZP4o/PryeRSHDNmfMjL47FBh7lBiADGazEnT+ThfcetldB+RmTxsTWpVF1+qeTDmLJmq18/xMLWLJmK//5gcMKjh/MZzZUnYQUjdxepNyVwB2q+iHrRjXGNZVHpb8I3LRuZ+EHORAPc74I5zzgDalQW7X24tgKtF5cqPq2+QooFhJuoPg2kYHyxDXB4dbl4hghNmh1fqZh5vrijMrKefFLFHZYkfIPKL9Ew36mxxleV/KlY1zJRZt+nDS3nbAeIX0lzeVntI6d03mnzHrXSBR64ux/nLPY7SvYjo1DVxOkVRW0v6J//3p3kbzhw+jJJt21ASZMoyNh20mJJY7nvWMO+00dyx+fX89pb57BYdMn9hs0+AOPSgYgjSrb6DoNdmBVK1rFBu1kTMR2bPo3fgERGQ8cA3wfQFV7q/TUEwg0k9DmAyOSEWKDFvp3YNjRajZoIvJpEXkKOEhEluS9VgJLGl6hRjLaDhp2GZuydqeglQlU3W1DSXS2DS/naMORVlHQ4iK357MvsB74oYi8CVgEXKSqRRdb94vknkrVTkXJDxDluXPyJ/nyPUI2T6wZIInEEKpshfj30/z7q2mbz2/vkycb76LlbtFv/n7qhOTt233/jrk273sz8ye3y73t/rWcslKu7n6oQp9BdW3vZO3ewdkK57MqHWC0XpNsDCPEBq2+z7RJmOnVSCnzXA+7XO9Hs4uV5WJtuXhZxZatZZxJmN32f4ALarzBWXVJYpQeiVO6KpXd/XIDkeuzCqo26G8uu6yS5p0mmei0p7Od234VCC5GmWeDJjG2Zy4/P6ZbtM8oYZqynSPrKWhOIbOfXb9vI2uDpokkOlj3lBVQrYImIidhlOMk8D1VvTym3JuBh4CPqOrNMadbArwP43kxXy3bnu/cY1jiFLSdZoDWYQdoWUmVfFJ19Zpyo8IAreVpmIImIneLyNNFXidXeIoUxj/0d1X1cGAnZtlIUVT1WlVdoKoLpk6dWoM7CAQGRiPbfH57HzcutPdA6zJcFLRmPtOmjG3JZWaBYY5qdQqaiCSBb2OWHx4MnC4iB8eU+y/gzjKn/KaqrgIOUNWX8l7De3AGeUsc7QDNLltOS+mBV3dfUNCGCg2bs1XV4+P2iUhc5PZ81gBrVPVhu30zJR5mJRmo8YcvLbjymUw/4zKnILgJLJeOBMP4QCHNavOVqrW1tH0qZxbivJr29BTfX2uK3X+j+l+lyuVIZbjYoDX1mWZt0JwiJlbzKRtDOMY2rRLHDwlrp1buh3ikpHmKWqSkxSliGrO/hjgFzCllg7VJczX0vTfGIQmntNntvH2KU+Xsl6S9hia8sk5pc+dytmbF0tZV0N4CLFPVFQAiciNmObAf7fwzwK+AN5c5X5+I/BCYJiLf9Heq6oXVVLalGTURJJG3xNG0j4yUjoPW3WcaYWdbq1g4BeJolUU1cZHbI1T1VRFZLSIHqurzwLvo36kDgaFCaPOBEckIWeIY+ndg2JEa3cG0ow6J3d/xUsc0EXksL+taVb02b3sakO/6cw1wRP45RGQacApwHOUHaO8FjrdlF5W9geFEIgmjdjNOQsjZoPVpklElDusKCtqQoVUGaEUjt4vI3pg1yu+x5T4D/Mx6u1oBfHJQV4sz/qg0YJSjhA2a897oQiuVCtUSaALN9+JY9zY/UBs012Z9wTi3Ejpng+bac5xy5vIrtbOKs4Pz6+pIxjxbSqlTcdeu1MFqws6bp1LVxUdrpq1ZK6h3I2Q1QWOeaRlrU2aVMSFjt4u3UXXl+6wSlLZeCp1tmkq0KMT1C5f6yonaVz/FzE9tnK5cvhfnLOPFQYvSmJ7Z78umRA9OZyDdX5EbsJLmHR+9u26/Z3Pm26apZ7umeSeUZGdhmZQzsDWzGOp55XPx0DQZE8AzkaRePjYlkWDXph2x+0eNGrW9u7t7QalTFMnzP8ArgC+oakbKKIGqugG4UUSeU9UnSxYejoyeEilobWKXOFJ64NXVl7FBqsOP0lanJQZoqroRM3vo568D3pO3/QRQqvMHAkOC0OYDI5XhssSxFKF/B4YjSs55zCBZA8zI254OrPPKLMAMugCmAO8RkbSq/p9/MhH5J1X9b+BvpUisiGG9xBGMHZp1EtJul8r2lflZ392XoTMV1LOhQEsM0BpOtTZogaHPMP4sRQoVpkbaRMUpauWUs0q7oKOvr7Lz51P2PiusRMKbT3feHQfbpBqpJrVCsx8hClpTEeswIE5Jq+gcMWn+/op0mrjlI5Gr45j9kfdC/yol7smVdepbQgqOVy20i4vCyJVR0nzEOmLwzyfOfsyzTXP5xW7Bt0tTV9YpaZEKZ/NtvCtJeEql3dY626BVGdL1UWB/EZkNrAVOA84oOL/qbPe/iPwI+F2xwZnlOZs+FrN/eDN6EmxcDkCbW+JYyQCtPQzQhgIt8LgOBAKBwEhhhNigBQLDkkwVyydVNS0iF2C8MyaBH6jqMyJynt1/9QDP91ubFh9JD3fGTIHVjwB5AzQts8SxNxNc7A8RRtYALc6wZbBx0Erg1vE7JWNIxUELjGgqaadxk+Ju4rZam8tKlTXf1nMgnilr1R/jFLVKKWdvN9yot4ImIpcCf4eJMQbwz6p6u913MXA2kAEuVNVybrxbl6zkbM8i26Y4NcopS4X7NcarI+SEJ4m2C7UTceKUVXDivDf2ty2zaaztme/l0dNs4gxei5bVgjzf9iwuP05J8706OiXNf9ejGkafi7VNi1HYzM6YLyQXMy3rFDRnz2ZVPOe90f0wjztPjajE42fp4/V24HYvr+jATFU/Uck5ReQA4PPALPJ+16rqcYOt55Bg9GRjg5bN0maXOPaWGaB192WDB8chQhgyBAKBQKBhNMgG7X9V9Wv5GTbe0mnAIcDewN0icoD669QCgUBRlHjfLU3ml8DVwPeAkdOfR08BzUDPVtpwSxzLOwkJCtrQYGQN0PxpW3+6faBT6kXioPleHN12ItG6M+MDnfEPDA18dSmOypt9EriNYpYn5WzPfE+RA/V4ONC+U+y4WsZ9K4WvqPmU62/NjOHWKJp0PycDN6pqD7BSRJZh4jL9pSm1qRJVibwxuqeQtDlVqngby/aYGEkZl/aatKuvHYAJb5hF2p4s43lz9OOfCUpWIdHbZbZ7u21aGOzQbUf5fb0mjbwj+4qbFs+P9nvb+bfq26D1KdENkVOuMjZNJjpLnrqcTZpT1HybNN/7oyQLlbMEueuq5wFSU9ZJerZ4J/HHR+LnJ1I5SbUOZOvkIbJK0qr63WZXouG4YNU7N5KyA7ReLf1A7+rLBBf7Q4Twy7wKzr311mZXIRAIBIYUTkEr9QKmiMhjea9zBniZC0RkiYj8QER2s3nFYjBNq/6OAoGRgarpv3GvRiMik0RkEvBbEfl7EdnL5dn84csDV8D2V8z/uzbSpsZzVu+zt5U8rCcM0IYMI0tBq9SVXFw5Px5aEU9JQzHemZvxD0rayMJv9m1txfdTYsmE6wI6wGUvcfHRKlW5XF19b46lrhPnvLXiqtcodt5gbNYGGtPOz28lKrRB26Cqse7nReRuYM8iuy4Bvgt8BfPRfgX4OvApKovBNDRQIAuatn3TxQJLx/RVq345mzPNWMXHU9p8lcxdCvorNYr35kW2aL4NWYyNWTF7sYL9xW+l35dNgR2dOwe5tMQnHGeT5u8vGyetnJIW7Xd2ZPkXKbxmZKfWz7DNHevUuBivjnWMgwbSagraImw4Prv9j3n7FNi34TVqFNPmwS8+Zv7ftZHkdjP31LvbASUP6+rLsFcYoA0JRtYALRAIBAJNpRZeHFX1+ErKich1wO/sZiUxmAKBQAx2XqBlyHfJP+KYfQz89TfgV2fDoutpW/EEAD27v7HkYV19GUYFN/tDgjBAy6fcNPUA4qG5Sb46hiOpOem0mbwMXTfQKAYaB62hilCNFLNyVOsFEqo3p20k9XYSIiJ7qapd+8MpwNP2/1uBG0TkGxgnIfsDj9SvJsOLWCGqUWvbXLdwpmmRy9j8Mp4NWl7RhOS8KfrnzBZuRkqYz0CVNL/qjmyxnZHBmjmX/9Mhsk1zHiA9r46a6rDb9Y2D1opOQkTkVOAOVd0uIv8CzAO+oqqPN7lq9eWAk0z64h2k9jwGXn2V3mzpL//gxXHo0MKP8UAgEAgMNxoQqPq/RWQuZkyxCjjXXFefEZGbgGeBNHB+8OAYCAyMVhygAf+qqr8UkaOAE4GvYbw6HtHcatWZtYtMOvNtJNc+Q1oT9JX5gLp7gw3aUCEM0OrMULRJCwQCQ5M42z5HK9ik1VtBU9UzS+y7DLisfldvEIJRXhK1+bXsYltlC/Jqcur+JHwprB7XKJ7txx8TvLhkUZwya9dllbBK46T5StpAEGtb5lL1bdPcPYhni2Y9QgqmU2kihdbJTqyF3ey7xvTXwHdV9Tc2HuKgEZFVwHZ77rSqLrCOR36Bibe2Cviwqm6u5jqDZuVCuPmT0DYG9ppLsmMP+l64nd5Xl1LK91Fwsz90CAO0QCAQCDSUVhgoBgKBgVO3gXt1rBWRa4Djgf8SkQ5q46X8naq6IW/7i8A9qnq5iHzRbn+hBtcZOGsXw6k/gl9/Gnq2kxw7lT6S9G16KfaQvkyWdFaHrILW19fHmjVr6G5AIM1a09nZyfTp02nzvbGVIAzQqmEA31QirWsX0s/mpZUrGxgSDCXbyyFU1bLEeapspe7coEDVwxvnrSHrtV73u6uIN0YA9VwDOuUsUtBijoP4fiLpXptal6rpPm+/3fa9ODr6eW+Mea5mvdSphyW8OGrW/C+AIP0UNN8mLSuF275N2kCVNL/qvnfHwn0lPD2Sp6gl/B94VjlzXh0lVb8vYG3ZKNAfBk4CvqaqW0RkLwo9OtaKk4Fj7f/XA/fRrAHaUZ81acdY6NlGoqONNEn6Zr0z9pDuPvPpDVUFbc2aNYwbN45Zs2YhQ+hHhqqyceNG1qxZw+zZlfu1CQvwAoFAINAwnA1aqVcgEGg9tMyrafVS3aWqt6jqi3b7FVX9Q7WnBf4gIovy4jDu4RwQ2XT3YgeKyDkuhuP69eurrEYZOsYZBU3T9JGiNxPvtKfLDtA6h6gXx+7ubiZPnjykBmcAIsLkyZMHrPy10LxqCxM3/VyikQyx9lMSp7D53uYCgcDQoJUGPbVwsz/iEZCkRiqSJLI2zY9Slot7Fh0mheUSyUI9pNiMrTtDQswxTmUTYlQ13/DabccZZCeiK9g00pSKl/evWsyLI+59YUDT0P1UrCh/cDZp/WruxUkr3JeKubbvedLuz/Z524W2afWixO//4cbbVXWdiOwO3CUiSys9UFWvBa4FWLBgQX3Hrh3joHsbibFmgNaXib9cT5/58DpTQ1ebGWqDM8dg6j10P6VAIBAIDDmCghYIDF1aUUGrB6q6zqavA78G3gK8ZpdPYtPXm1dDS8c46N1BIpsmrUl60+UVtBAHrTZceumlfO1rX6vb+YOCVglxvxiKzAj6Wf6geUioUQ3wgx0Ynrj2X6/QSHHx0JL2eZMZhIFE9MMizmCrFQ25YoirYlPiyMUQbNBqg+bZi2nWdDxJZu229+CJs0nLFj6w8rut60pxdmmxP8r72ZjF2J6VI654P9vvIvWLOda3OVPP5sz36hiZfVXp3bEUcZ4ffSUtmk63tmbqbNGiSuYraPXz4lhCoBk2iMgYIGHjqo0BTgC+jImleBZwuU1/07xaWtrNEkfJ9pKWFH2lljj2Dm0btIFw9f3LOWz6BI6cMyXKe3D5Bpas2cp575jTxJpVTlDQAoFAINAwgoIWCAxdVONfw4g9gAdE5ElMMPvbVPUOzMDsr0TkReCv7HZzsTZoZNLGSUiJAZpzEjJUvTgOhMOmT+CCGx7nweXGCeeDyzdwwQ2Pc9j0CVWd97LLLuPAAw/k+OOP5/nnnwfgzW9+M/fddx8AF198MZdccklV13C0xJRwpbElROT/AX+Lmch5CvikqjZvLrZeMkGL0dJq3xClGW2+GpVpRDGEFLNytOItjAQbtLr3b+vFMdvn7JbMd3Q27pFuFTXNWOUnbeZmsxmTZqx3x3wHiu4Ht/v29zWZrApZlchLo7h+U05BS2cKU997ozvcSTX+NHK/AXxepXMuEU2SptDtoO/FsYxXR3dtX0krR5ySVio+WqTSRXHQYuzhvDho/WzR6sxIeHyo6grgTUXyNwLvanyNShAN0HrJSGkbtK4RNEA7cs4UrjrjcC644XE+dsRMfvrwy1x1xuEFitpAWbRoETfeeCOPP/446XSaefPmMX/+fH70ox/xoQ99iG9+85vccccdPPzwwzW5h1ZR0Fxsif2Be+x2ASIyDbgQWKCqh2IcCp/W0FoGArUjtPnAiMQtcSz1GgaE/h0YdrhA1XGvQBPoGAcodG8hXcaL41B3sz9QjpwzhY8dMZNv3ruMjx0xs6rBGcCf/vQnTjnlFEaPHs348eN5//vfD8AhhxzCmWeeyfve9z5+8IMf0N7eXovqt4aCRuWxJVLAKBHpA0YD6xpRuWgauoLpaDcZ6NQKNwsZZ4vmiFOpGmmzlk4XX6YQV4chYU/XujS8zY905axiO6xy/bzOClux/tQvVuEgaQVFbYSYuNa/f2clp4xZ40/xfiX7tmhOQYtS55HRemhM5BUv53RM8Z4XWWuf5cVBiz7sjKewRYqa99DxpTv/ujZffKePkBcbLW9f3v5+HhGdSuVUK8+LY6RiJQvtxPzeGPcEjFXS8u/H8+jox0HzVT1X11hbNHOSmBpVSRiItR4d40y6axNZ6aCvhJOQbufFsa1VtJn68uDyDfz04Ze58Lj9+OnDL/PWOZOrHqTFeWN86qmnmDhxIq+99lpV58+nVT6lsrElVHUt8DXgZeAVYGsN4lwEAs0itPnAiMQtcRzmClro34FhSVDQWgw3QOvaZJc4Bi+OkLM5u+qMw/ncCQdGyx2dTdpgOOaYY/j1r39NV1cX27dv57e//S0At9xyCxs3bmThwoVceOGFbNmypSb30LD5VBG5G9izyK6KrOlEZDfMrORsYAvwSxH5mKr+NKb8OcA5ADNnzhxMlWvKUJ41DgrZ4Ghkm89v71OmNL+9D2mGkQ1aq6I69L9TmvpM2w2jFiW8X8WRgmSVsUTxeGg+TkkrtoIibhZXsEJN1rcls2lc3LMaESlp+ZnlLunbmFkvjlFz7OeFubiy1q8unnfH6HL2OF9J6818P+8a8bHRIE8xy7prlYuL1kY9vTiGgViLEQ3QtpBJ7lvaBm0EeXFcsmZrgc2Zs0lbsmbroFW0efPm8ZGPfIS5c+eyzz77cPTRRwPwxS9+kXvuuYcZM2ZwwQUXcNFFF3H99eU9uJajYb9AVPX4uH0i8pqI7KWqr5SILXE8sFJV19tjbgGOBIo+zBoaKDAQKEIj23x+e58zJ7T3QCujQG+zK1E1zXymzZ/hj8wCgcYw9KdWhhlugIaSTZSxQUuPHCchxVzpHzlnStVLHC+55JJ+Xho///nPR/9feOGFVZ0/n1aZIq4ktsTLwFtFZDTQhfGk81jDaliMCuKgOUQKJ+QrtTlrlO0ZmPqJmxbNq2ycrVlQ1qqibm1eBNractvJmO/iym3S4r/M/bhnw8zVcmkaqLRV2tdqZatWX5QR4Auuvs80FTSdjJSyCGeT5tueaWXtoph5hZuQT3qqjIqa/u4rZ37qbM8iL49amOYqafO9CgzmMRN3TLbQPs55THR35itpGtmcdRbk+6f346TFVccpZ+3Js6N9Ls8pacVM6/J3iLM9E6/uqVG5svWyQUPRYReSeogTDdBMmyjpZt8qaB2pofCcCLTKp1Q0toSI7C0itwOo6sPAzcBijDviBHY2MRAYgoQ2HxihKNBX5jXkCf07MOxwgarjXoEm0D42+jeTaCtrgzaqLRnr6CLQWrSEghYXW0JV1wHvydv+EvClBlYtEKgLoc0HRi7DX0EL/TswLFHQEbVMYgjQMT76VyVFn+8VNY/uvuyI8eA4HGiJAdpwptqJiuDKPtCKtOoEXAjGPRSorw2aiJwKXAq8AXiLqj6Wt+9i4GzMCPFCVb3T5s8HfgSMAm4HLtIh8EvUBap29HMG4i1tzPYlC47Lpk3aZwNW9+Yd71Ygugl5f5VkAqFPBentAYhS+rzPttcqoi4wtV1mFUku7vz5rvEB0rYunrmd2rq60ACav9/5z3Dn7EkivXlLtDPWRWjkWMM6CcG50+8srIu7prjlmTbDey/iljxGx3tORPKdhLjljr4L/jhX/s5RScLV1TkNcfcmKaSuTkKq6xYichJwJWbt/PdU9XJv/0fJhaTYAXxaVZ+s6qLDmY6cglbOBs0paIGhQRhKBwKBQKCB1H2J49PAB4CF+ZkicjAmEPQhwEnAdyTndu+7GA+J+9vXSdVWIhAYjmRUY1/lsP3t28C7gYOB022/zGcl8A5VPQz4CmHZb2lSHZDsAExMvHJLHDtHgIv94UJQ0JrEQJWxZgSLjnM4EJcfVL6hyWB8XLSKtuCcofQ102yphdzyx30/tBb1XeKoqs9B0YCiJwM3qmoPsFJElgFvEZFVwHhV/Ys97sfA3wC/r1slq0UUSVT/fetUqGg77//Ih4cU7qt11y8bcSEmREAUZLtYoOpyOHf7ku+aPs+dvivnudnPHW93xwS4jgJal1HSIN4Ff5z7fd/tfnRu5zQkkYI6KWhQ9ef/FmCZqq4AEJEbMf3y2ej8qg/mlX8ImF7dJUcAHWNhV0/ZAVpPX4bOVBigDRVa8ekdCAQCgWGLW+JY6sUUEXks73VODS48DVidt73G5k2z//v5gUAgDzO1orGvnp6eCWX6bVwfjONsWnmipFWwnhw1kSodB60vMyKCVDeKSy+9lK997Wt1O3/zp3yHEn6k6Wz/mQo3aeu7208kWmKCvSTpdKEy4mbj/Vn4uPzAcCeBm5ktZ4NW5xi1gQpp3b5aVv3ZoKoL4naWChKtqsVc2kNxWUFL5Lcu1s1+ptcoP04Jk342aIXbGWd7Zu24MmnzYy2dtbZpxQJVe6d0mwn76md71m233ZdAZINmT97ngmfbE3mP1ShUQCbmS8YPIZDXtCOVzwXeTichkyAh5j6kZ2vhsamc3RYAKWvX5Stf1oV9FKi6QK3K5fuBrDVRPD+/N8a54C/rft/LcDZp9iDqRbZE1+jo6Ni6a9eu2H7LAPqaiLwTM0A7akAVHIlEA7R2etMlljj2jqAB2gNXwLR5MPuYXN7KhbB2MRz12WbVakC02lM7EAgEAsOaihS00mdQPV5VDy3yihucgZmtn5G3PR1YZ/OnF8kPBAJ5KEq2xF8FxPXBAkTkMOB7wMnWI2qgFNaToyTLxEHry44cJyHT5sEvP2EGZWDSX37C5FfBZZddxoEHHsjxxx/P888/TyaTYd683DlffPFF5s+fX9U1HC2u6bQYA5DAnLiW71Wuu7tIuRrZog2GcrPqA7VBiysXbNNGHkXE5UDA0jQ3+7cCN4jIN4C9Mc5AHlHVjIhsF5G3Ag8DHwe+1YwKDpTIa6MLYhzt8FSmbOF3tnr7k9aerZjw4lQ1X0nL2le0sqRfIOpM4bZvuNrA74ismpeknbpnH8xWOYxUuLR9SKc6vROkC8rFBbZ2OMVMna1bCa+Pbpfv4XGwSppR9eo3917lx/YosL+IzAbWYpz2nJFfQERmArcAZ6rqC9VdboTgYqEl20p6cezuy9A5UgZos4+BU39kBmULzobHvm+28xW1AbJo0SJuvPFGHn/8cdLpNPPmzWP+/PlMmDCBJ554grlz5/LDH/6QT3ziEzW5haCgBQKBQKCB1NeLo4icIiJrgLcBt4nInQCq+gxwE8YhwR3A+ZrzwPBpzIz9MmA5we4lEChKNQqaGm8qFwB3As8BN6nqMyJynoicZ4v9GzAZ42X1CRF5LOZ0AUfeEse+Ikscr75/OQ8u32C8ONoB2oPLN3D1/csbWs2GM/sYMzhb+N8mrWJwBvCnP/2JU045hdGjRzN+/Hje//73A/C3f/u3/PCHPySTyfCLX/yCM844o8yZKiMoaDXGxV9yylml3u4GqzbVRZ1ShXQ62KCNUCoRilvFi2McIR5aIa3VZ+vuxfHXwK9j9l0GXFYk/zHg0LpVqtaoiQfm4plJwtl1+TJX4eftbM6cLVraxkHL2HL5R7su7kzHfAeJCYG+LDm53KVp77O1TgvUrVy1opKLZ9Yv3lk2z36sGOobxeVVLGWu7RRD7UtBOhnZoEVE6p71quhOnSyMjxbFS6sQ39ujr6RJifPFxUrzlbSydcimqZcJpUKlSxnjz6F6OybWYH7e1Xn//y3wt1VdZKRhB2iSbCtwEnL1/cs5bPoEDps+gQtueJzedJatXb1cfMsS7nzmNa464/Bm1bgxrFxolLNj/smks4+uepBWxDswH/zgB/n3f/93jjvuOObPn8/kyZOruoajFZ7WgUAgEBgxVG+DFggEmkHVNmiBemAHaKu39tKbyfLd+5bx4PINHDZ9Auf+ZBG/fXIdH5o/jR09ae5d+jq/W/IKV51xOEfOmdLkitcRZ3N26o/guEtyyx2dTdogOOaYY/j1r39NV1cX27dv57e//S0AnZ2dnHjiiXz605/mk5/8ZC1qDwQFreYkvUk/N9huRcXBzaqnUmac7jupHOh5Ks0PBJqJUwgH295jaaF4aHG0Tp8M0mYtiLwVWo+HSasm+TZmTlmLlCV3nE0zNk3nHRcJY148NEdC7XNNPeXMb2Kun7l8p/JlXd0Sxff7Spl3Lw5J5O9z55JcWRUyagIsRzZo0U04W7TCLwNVp6R5fdluR14do+t6lYyb+i5isxZ3aJyS1pM2cZtdjDV3vJKv0jXHi2OgSdgB2tjRowHzu/OCGx7n08fui6ryq0Vr6LXKWiYLnzxy1vAenIHx1phvc+Zs0tYuHrSKNm/ePD7ykY8wd+5c9tlnH44++uho30c/+lFuueUWTjjhhOrrbmndXxKBQCAQGIY4G7RAIDCUUCAjYXKl5bADtMkTxgBw9f0r2HfKGC67bSmdbYlocNaREs49Zg4/ffhl3jpn8vAepBVzpT/7mKqXOF5yySVccskl/fIfeOABPvWpT5H0VZoqCAO0GuOUMl85E6lsYn2wNiKV2rAVO3+lSkLcNcrZqMVtB4YuIQ5aEXxpLk6qa2GFrTFkCcsYq0SM3Vki6bwv2jhoKc+mKuviFhb34piw3hvbk+a4zkTuuzlhO3Gq0FFk3n5IJoh32ZqISVNaUEcS3g9+p6xZGzRJeOf37OoK9ttzO7u2hEr0ngBIplBBi+zcXLyyyMujPSZhvTdGdU25A03i4psl2rz9qYL9zvbM2aiJ5r4DxCsT4xgyUs46Uib2s1PWonJ5cdDqp58pmRGufovIScCVQBL4nqpe3uQqRQO0RLIdgHEdKRa/vIWEGNf6yQSkEkJ7Kslb50zmrXMmc8ENjw//ZY4N4pRTTmH58uXce++9NT3vSPwZFQgEAoGmki3zCgQCrYYCKtnY13BHzLrTbwPvBg4GTheRg5taqQeu4NEnlwAgKTNAW725i3eOWUlWoT2VIJlI8PkTD+SaM+dzwQ2PA3DVGYezZM3W2NMGKufXv/41S5YsYcqU2g52R/pUbs3x45+57UTCTKSXM1OJ2+/nOzUqLn8gxClog72Gv7/cdlDUWoNaCjuNioPWF7NSrpT3xprbnvn4b+SIV8x8nJOQwGBRFbK9KdK9RrlxClpS82yvvPIA2Yx5MGX7bGrVKOfFMd8GzXnrjkzDvDqkssaehV7bCXttx+ortFFyjg21r9DWLFKvolhunhdH5+XRk+P7earM3x8Zztn77U0h6STq4qD1dtljCmOKqd3WZAcAEilohcuV+tmm+V4Zs332fMUVteg8kttOOOXLk8762abZY+K8O0bh0LSeX3BKhnp/gbY0bwGWqeoKABG5ETgZE7qjOUybx+H3fw2Apa/vMlmygZ27upg7YwLL1+/kwnftx3fvW8FVZxweDczOe8ecoJ61OEFBCwQCgUADcW72S70CgUAroiX+RgDTgNV522tsXvOYfQypk0zkkJ41TwAwIdHNX73nA/zf+UdxzZnzyWRzitmRc6Zw3jvmNLHC1aGt6HGvAgZT7zC1OxjcjHgRQxs/a6C2OAOdbK+0fKk4SOW82g1WAIizVYtTzIKNWmCgtNmJal9JKxUHre5eHCu90IhV1oKTkHoRp5xF3hudOqVmO2tVKufFcVA/fZz3RhfvzG/uLs6Z76XRJ8Z7Y7nYbpr3vIhiFGleqoJaj5OSdu3Oa3+eUubiozllrZ+9XGSb5uy+um1+tMzEpsUVtYJ4c/4bFmOE5meX9O5Yzkh4kGhQ0Iq9sf26jYicA5wDMHPmzHrXCeafBS8/xDlP/ILVycnscczf8XdHm0HYkXOmRErZUFfMOjs72bhxI5MnTy4aj6xVUVU2btxIZ2dn+cJ5jNRfCIFAIBBoCvUNVB0IBOqHjuwJ1DXAjLzt6cA6v5CqXgtcC7BgwYL6Sz4rF9K39A6uT5zCpaNu5rMPT+fB/aYM+QGZz/Tp01mzZg3r169vdlUGTGdnJ9OnTx/QMS0xQBORU4FLgTcAb1HVx2LKtZ73HI9Stjf5s/yVhkyq1Sx/Ke+N+Wk2Sz+Xk5XaxcVds5ySVksbtaGiwg2nNg+tGecvjrrboPnEeXccsQx/G7S692/jrSFSwJwNWr9ink2as+tyx7n9KesJMZV3noq7tB//zPvqjWKyebHYyBSqemJt0CLFLFt8+Um/OGj5zwvvHJpNQFYi1zOStu3OKWb+dnQee09qbdPctWy+pqytmlXOYm3T4hQ1G0cNgIw9h7Ux62eT5u7bKVd2v7tr37uj264PI15BexTYX0RmA2uB04AzmlqjlQvpu/Eszu+7kE987ExGJ57l2zeexfk/Az565rAapLW1tTF79uxmV6NhtIoN2tPAB4DYEN8t6T0nEBg8oc0HRjDD3gYt9O/AsMMMVEeu91U161EvAO4EngNuUtVnmlqptYv5/Rv+k0+4wdjsY2g77Xr+ee6u4KVxiNMSU7mq+hxQbk1p63jPKTEF7xS0VpskL6ZWxSliPnH3MtB7jFO1aql2tbpy5mhEm4+zz6oHfty/VmbQNmgDVcDiLjDi46INfxu0uvdvAURz8c+scpSw8cycahXFQ5NCWzNXvthp43CtNLKFcnHQUs7g07Zrb9pXslbVSrhrmDNEqpQ7ILoXW+eY73LxLpAf5wwXF86WkUQWSShJgaQQKWXqpXhplC++IuaXS8WUi+nTzhat8AZsXQuPiWKnZdNF9ztvjc67Y6GSVq9JDiVTVy+RrY+q3g7c3ux6RBz1Wd7v580+hlmzj+G8ZtQnUDOG0i+DYt5zjmhSXQC45oMfzG18/vMA3HufeTA4W8AaBhVvCFdeeQ3jN6xodjUChpZr8/BeAK688t8BeOIJk1tsuexQ4eRmV2DEEQJVW1quf7/+TK46K9bbH/6j7TJIcfqJ2W6LGeS1GvrKk6x88tBmV2OYoHV24x8IBBwNG6CJyN3AnkV2XaKqv6nkFEXyYp8Q+V50gB5JJp+u4BqD5wtfqKTUFGBDXetRGaEehRxYj5M2ss377f2DH5S6tvcJE75cadFW+YzL1uNfJ0xoiXo0iLq0+coZ+ssYm/1MG/2lzU/D5gouUwl/jt+1reSBU/7nYy3Rnsv0q28BcOP8ZtejYdSlfyuMFHf6NWPRokUbROSlIrtapa3kE+pUGbWs0z5xOxo2QFPV46s8RUXec/KuF3nREZHHVHVBldevmlCP1q1HPc7byDbfiu0dWqcuoR7969G8qw+PJY7hmRbq0cr1qM+ZlezIdhIyYFR1arH8Vmkr+YQ6VUaj6tQqTkIqIfKeIyLtGO85tza5ToFAPQltPjAMCYGqLaF/B4YcqtnYVyAQqB0tMUATkVNEZA3wNuA2EbnT5u8tIrdDi3rPCQQGSWjzgZGLc7Nf6jV4RORUEXlGRLIisiAvf5aIdInIE/Z1dd6++SLylIgsE5FvSpVRUEP/DgxH1Nqgxb0CgUDtaAknIar6a+DXRfLXAe/J2x6s95x6BgYZCKEehYzYetS5zbfK+wqtU5dQj0KaXI+6zrY7F/fXFNm3XFXnFsn/Lsa+6yFMfzsJ+P1gKxCeaQ0n1KOQutUj2KDVjFZpK/mEOlVGQ+okOhT8YgcCgUBgWCCyh8JHy5T630XVrvEXkfuAz7sg0SIyC/idqh7qldsL+KOqHmS3TweOVdVzq7l+IDDcmDhhsh580Ftj9z/3woMvbt68+YAGVikQGLa0xBLHQCAQCIwkytqgTRGRx/Je58SeamDMFpHHReR+ETna5k3DOOxwrLF5gUDAQ8nGvgKBQO0Y0gM0ETlJRJ63dgNfLLJfrD3BMhFZIiLzKj12kPWZJCJ3iciLNt2tSJkZIvJHEXnO2klclLfvUhFZm2cj8R7/+FrVw5ZbZe0unsj3+lTp8bWoh4gcmHe/T4jINhH5rN1Xq/ejqE1KkXJF20St3o9aENr84Othy4U2X1iuCW3eeXEs9WKDqi7IexUsKRGRu0Xk6SKvUmHtXgFmqurhwOeAG0RkPAN0d19PQv8efD1sudC/C8vVtH8rSpZM7CtQGfXoq4OoQ9F+W9/v/orqlRQzifa7FqnPRBG5WUSW2vfqbY2q05AdoIlIEvg28G7gYOB0ETnYK/ZuYH/7OgdjZ1DpsYPhi8A9qro/cI/d9kkD/6CqbwDeCpzvXft/VXWufQ02Wn0l9XC8014r/0t+IMdXVQ9Vfd7dLzAf2EWh7UYt3g9nk7IwrkCZNlGr96MqQpuvuh6O0OZpZpuv3oujqh6vqocWecXGH1PVHlXdaP9fBCwHDsAoZtPzipZ0d18vQv+uuh6O0L+pY//WTPwrUJY69tWBEtdvm/175yKM0yRHs+tzJXCHXQL/Jlu3htRpyA7QgLcAy1R1har2AjcC/uzpycCP1fAQMFGMvUElxw6Gk4Hr7f/XA3/jF1DVV1R1sf1/O+bDrvVymrL1qPPxgz3PuzBG/MWCOg4aVX1OVZ8vU6xUm6jV+1Etoc1XUY86Hz/Y84zANl9fL45xiMhU++MIEdkXM8hZoaqvANtF5K0iIsDHgUoCTdea0L+rqEedjx/seYZZ/1ZUM7GvQEXUq68OiBL9tmm/d0RkOvDXwPfysptZn/HAMcD3AVS1V1W3NKpOQ3mANg1YnbddzG4grkwlxw6GPezDHpvuXqqwGKP1w4GH87IvELN05QdVyKaV1kOBP4jIIim08RjQfdSgHo7TgJ97ebV4PyqhVJuo1ftRLaHNV1+P0OZzNKnN193N/ilSxMU95kG7RESeBG4GzlPVTXbfpzE/CpZhlLVBe3CsgtC/q69H6N856tC/FSUT+wpURL366qDx+m0zf+9cAfwThW5+m1mffYH1wA/tssvviciYRtWpJdzsD5JK7Abiygza5kBE7gb2LLLrkkqOzzvPWOBXwGdVdZvN/i7wFVuXrwBfBz5Vx3q8XVXXicjuwF0islRVY5dM1LEeiAnU+n7g4rzsmrwfpZY95Z+iSF6ruTgNbT60+Yrq0dptXqmnm32Nd3H/K0z7K3bMY8ChxfY1kNC/Q/+uqB5N7d9VBqQWkZMwy8aSwPdU9XJvv9j978EsD/2EU3qGCS31W8Pvt1JdCMhq6vFe4HVVXSQixzalEv1JAfOAz6jqwyJyJQ1cYjmUB2hrgBl528XsBuLKtFdwbFFU9fi4fSLymojspaqv2GUnr8eUa8N0iJ+p6i15534tr8x1wO/qWQ81MXlQ1ddF5NcY6X0hUNHxtaqH5d3A4vz3oFbvR4WUak8DuY96Etp8lfUIbb6AJrV5p6AFPEL/rrIeoX8XUJf+XY1SJjn7q7+y9XtURG5V1WfziuXbWR6BGdQeMeiLth6V9POGENNvm/V75+3A+8U4zukExovIT5tYHzCf1RpVdSsCbsYM0BpSp6G8xPFRYH8RmW1nqk4DbvXK3Ap8XAxvBbZaObKSYwfDrcBZ9v+zKGLHYGeHvg88p6rf8Pbtlbd5CsYQuF71GCMi49z/wAl51yt7fK3qkcfpeEtBavh+VEKpNlGr96NaQpuvrh6hzRfSpDZfvZOQYUro39XVI/TvQurQvxW0xKs81dhZDhfq1VcHRIl+25TfO6p6sapOV9VZmPfkXlX9WLPqY+v0KrBaRA60We8Cnm1YnVR1yL4wEvgLGJuBS2zeeRjbAjBS8rft/qeABaWOrUF9JmM8urxo00k2f2/gdvv/UZhfKEuAJ+zrPXbfT2w9l2AawF51rMe+wJP29Uz+exB3fD3qYbdHAxuBCd7xtXo/TsHMhPQArwF3xtSjaJuo1fsR2nxo86HNK8AdwGNlXnc0q4818xX6d1X1CP27zv17ypQpD02cOPGFuNeYMWNWeP34HK/eH8Isa3TbZwJXeWV+BxyVt31PfjsfDq969NVB1KFov61VX6mybscCv6umrdawLnNtW14C/B+wW6PqJLYCgUAgEAgEAoFAXRCRU4ETVfVv7faZwFtU9TN5ZW4D/lNVH7Db9wD/pCY0RiAwYhjKSxwDgUAgEAgEAkODauwsA4ERRRigBQKBQCAQCATqTTV2loHAiGIoe3EMBAKBQCAQCAwBVDUtIhcAd2Lc7P9AVZ8RkfPs/quB2zG2UMswbvY/2az6BgLNJNigBQKBQCAQCAQCgUCLEJY4BgKBQCAQCAQCgUCLEAZogUAgEAgEAoFAoCQicqGIPCciP2t2XYY7YYljIBAIBAKBQCAQKImILAXeraor8/JSqppuYrWGJUFBCwQCgUAgEAgEArGIyNWYgPC3ishWEblWRP4A/FhEZonIn0RksX0daY85VkTuF5GbROQFEblcRD4qIo+IyFMiMseWmyoivxKRR+3r7Tb/HSLyhH09LiLjmvYGNJigoAUCgUAgEAgEAoGSiMgqYAFwAfA+4ChV7RKR0UBWVbtFZH/g56q6QESOBf4PeAOwCVgBfE9VvyQiFwGzVfWzInID8B1VfUBEZgJ3quobROS3wOWq+mcRGQt0jxS1LrjZDwQCgUAgEAgEAgPhVlXtsv+3AVeJyFwgAxyQV+5RF8tORJYDf7D5TwHvtP8fDxwsIu6Y8VYt+zPwDWvzdouqrqnXzbQaYYAWCAQCgUAgEAgEBsLOvP//H/Aa8CaM+VR33r6evP+zedtZcuOQBPC2vAGf43IRuQ0TG+8hETleVZfWqP4tTbBBCwQCgUAgEAgEAoNlAvCKqmaBMzGByAfCHzDLJgGwShwiMkdVn1LV/wIeAw6qTXVbnzBACwQCgUAgEAgEAoPlO8BZIvIQZnnjzjLlfS4EFojIEhF5FjjP5n9WRJ4WkSeBLuD3NatxixOchAQCgUAgEAgEAoFAixAUtEAgEAgEAoFAIBBoEcIALRAIBAKBQCAQCARahDBACwQCgUAgEAgEAoEWIQzQAoFAIBAIBAKBQKBFCAO0QCAQCAQCgUAgEGgRwgAtEAgEAoFAIBAIBFqEMEALBAKBQCAQCAQCgRYhDNACgUAgEAgEAoFAoEX4/2JVtzQtkS/fAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "shiftshape (68, 2)\n" + ] + } + ], + "source": [ + "# Step 6 of the algorithm\n", + "w_diag = np.atleast_2d(np.diag(W))\n", + "W_n = W / np.sqrt(w_diag.T*w_diag)\n", + "\n", + "# Step 7 of the algorithm\n", + "nr = np.arange(W.shape[0])*stride + startI\n", + "coords2, weightmatrix, DX, DY, row_mask = Reg.threshold_and_mask(min_norm, W, DX_DY, nr)\n", + "\n", + "# Step 8 of the algorithm: reduce the shift matrix to two vectors of absolute shifts\n", + "dx, dy = Reg.calc_shift_vectors(DX, DY, weightmatrix)\n", + "\n", + "# Interpolate the shifts for all values not in coords\n", + "shifts = np.stack(Reg.interp_shifts(coords2, [dx, dy], n=original.shape[0]), axis=1)\n", + "neededMargins = np.ceil(shifts.max(axis=0)).astype(int)\n", + "\n", + "plot_masking(DX_DY, W_n, coords2, dx, dy, shifts, min_norm, sigma)\n", + "\n", + "print(\"shiftshape\", shifts.shape)\n", + "shifts = da.from_array(shifts, chunks=(dE, -1))" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "c6222f17", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Array Chunk
Bytes 313.09 MiB 92.09 MiB
Shape (68, 1231, 1961) (20, 1231, 1961)
Count 129 Tasks 4 Chunks
Type uint16 numpy.ndarray
\n", + "
\n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + "\n", + " \n", + " 1961\n", + " 1231\n", + " 68\n", + "\n", + "
" + ], + "text/plain": [ + "dask.array" + ] + }, + "execution_count": 56, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Step 9, the actual shifting of the original images\n", + "\n", + "# Inferring output dtype is not supported in dask yet, so we need original.dtype here.\n", + "@da.as_gufunc(signature=\"(i,j),(2)->(i,j)\", output_dtypes=original.dtype, vectorize=True)\n", + "def shift_images(image, shift):\n", + " \"\"\"Shift `image` by `shift` pixels.\"\"\"\n", + " return ndi.shift(image, shift=shift, order=1)\n", + "\n", + "padded = da.pad(original.rechunk({0: dE}),\n", + " ((0, 0),\n", + " (0, neededMargins[0]),\n", + " (0, neededMargins[1])\n", + " ),\n", + " mode='constant'\n", + " )\n", + "corrected = shift_images(padded.rechunk({1: -1, 2: -1}), shifts)\n", + "corrected" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c236f861", + "metadata": {}, + "outputs": [], + "source": [ + "# Optional crop of images TODO" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "efaeabad", + "metadata": {}, + "outputs": [], + "source": [ + "# Save as png with dask\n", + "from pyL5.lib.analysis.stack import da_imsave\n", + "os.makedirs(os.path.join(folder, subfolder), exist_ok=True)\n", + "da_imsave(os.path.join(folder, subfolder, 'image{:04d}.png'),corrected, compute=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "7054cea6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Saving stack as .png\n", + "Now saving movie stack.mp4\n" + ] + } + ], + "source": [ + "from pyL5.solidsnakephysics.helperFunctions import save_movie\n", + "save_movie(corrected,os.path.join(folder, subfolder))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.8.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/registration/registration.py b/registration/registration.py index 328a815..f072460 100644 --- a/registration/registration.py +++ b/registration/registration.py @@ -39,6 +39,16 @@ def crop_and_filter(images, sigma=11, mode='nearest', finalsize=256): result = result[:, sigma:-sigma, sigma:-sigma] return result +def crop_and_filter_extent(images, extent, sigma=11, mode='nearest'): + """Crop images to extent chosen and apply the filters. Cropping is initially with a margin of sigma, + to prevent edge effects of the filters. extent = minx,maxx,miny,maxy of ROI""" + result = images[:, extent[0]-sigma:extent[1]+sigma, + extent[2]-sigma:extent[3]+sigma] + result = result.map_blocks(filter_block, dtype=np.float64, + sigma=sigma, mode=mode) + if sigma > 0: + result = result[:, sigma:-sigma, sigma:-sigma] + return result def dask_cross_corr(data): """Return the dask array with the crosscorrelations of data