{ "cells": [ { "cell_type": "markdown", "id": "71990519", "metadata": { "papermill": { "duration": 0.002836, "end_time": "2024-11-19T22:31:12.550591", "exception": false, "start_time": "2024-11-19T22:31:12.547755", "status": "completed" }, "tags": [] }, "source": [ "# GE Discovery MI (Listmode Reconstruction; With Time of Flight)" ] }, { "cell_type": "markdown", "id": "7c513fdf", "metadata": { "papermill": { "duration": 0.00225, "end_time": "2024-11-19T22:31:12.555794", "exception": false, "start_time": "2024-11-19T22:31:12.553544", "status": "completed" }, "tags": [] }, "source": [ "This tutorial demonstrates how to use PyTomography's algorithm for VOI based uncertainty estimation in PET reconstruction. The first part of the tutorial is the same as the introductory PET reconstruction tutorial for the discovery MI data. The data here can be downloaded from https://zenodo.org/records/8404015" ] }, { "cell_type": "code", "execution_count": 1, "id": "ea48e0e5", "metadata": { "execution": { "iopub.execute_input": "2024-11-19T22:31:12.563036Z", "iopub.status.busy": "2024-11-19T22:31:12.562415Z", "iopub.status.idle": "2024-11-19T22:31:16.653254Z", "shell.execute_reply": "2024-11-19T22:31:16.652547Z" }, "papermill": { "duration": 4.096959, "end_time": "2024-11-19T22:31:16.655415", "exception": false, "start_time": "2024-11-19T22:31:12.558456", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "import pytomography\n", "from pytomography.metadata import ObjectMeta\n", "from pytomography.metadata.PET import PETLMProjMeta\n", "from pytomography.projectors.PET import PETLMSystemMatrix\n", "from pytomography.algorithms import OSEM, BSREM\n", "from pytomography.priors import RelativeDifferencePrior\n", "from pytomography.likelihoods import PoissonLogLikelihood\n", "from pytomography.transforms.shared import GaussianFilter\n", "import matplotlib.pyplot as plt\n", "from pytomography.io.PET import clinical\n", "import torch\n", "import itk\n", "import numpy as np\n", "from pytomography.callbacks import DataStorageCallback" ] }, { "cell_type": "markdown", "id": "d558cd1d", "metadata": {}, "source": [ "This code is from the introductory PET tutorial for discovery MI data" ] }, { "cell_type": "code", "execution_count": 2, "id": "2d2bd0c3", "metadata": { "execution": { "iopub.execute_input": "2024-11-19T22:31:16.699534Z", "iopub.status.busy": "2024-11-19T22:31:16.699042Z", "iopub.status.idle": "2024-11-19T22:31:55.254439Z", "shell.execute_reply": "2024-11-19T22:31:55.253351Z" }, "papermill": { "duration": 38.561221, "end_time": "2024-11-19T22:31:55.256256", "exception": false, "start_time": "2024-11-19T22:31:16.695035", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "scanner_name = 'discovery_MI'\n", "info = clinical.get_detector_info(scanner_name)\n", "\n", "# Get listmode events\n", "detector_ids = clinical.get_detector_ids_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/LIST0000.BLF', scanner_name)\n", "# Get scatter/random correction term at each event\n", "additive_term = clinical.get_additive_term_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/corrections.h5')\n", "# Get multiplicative weights (attenuation/normalization) at each event\n", "weights = clinical.get_weights_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/corrections.h5')\n", "# Get ALL valid detector-pairs and corresponding multiplicative weights\n", "detector_ids_sensitivity, weights_sensitivity = clinical.get_sensitivity_ids_and_weights_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/corrections.h5', scanner_name)\n", "\n", "# Define object space reconstruction matrix\n", "object_meta = ObjectMeta(\n", " dr=(2.78,2.78,2.78), #mm\n", " shape=(192,192,71) #voxels\n", ")\n", "# Get time-of-flight metadata and define projection space metadata\n", "tof_meta = clinical.get_tof_meta(scanner_name)\n", "proj_meta = PETLMProjMeta(\n", " detector_ids,\n", " info,\n", " tof_meta=tof_meta,\n", " detector_ids_sensitivity=detector_ids_sensitivity,\n", " weights_sensitivity=weights_sensitivity\n", " )\n", "\n", "psf_transform = GaussianFilter(4.5) # 4.5mm gaussian psf\n", "system_matrix = PETLMSystemMatrix(\n", " object_meta,\n", " proj_meta,\n", " obj2obj_transforms = [psf_transform],\n", " N_splits=8,\n", ")\n", "likelihood = PoissonLogLikelihood(\n", " system_matrix,\n", " additive_term=additive_term/weights,\n", ")" ] }, { "cell_type": "markdown", "id": "bb8c61cc", "metadata": { "papermill": { "duration": 0.002528, "end_time": "2024-11-19T22:32:00.698885", "exception": false, "start_time": "2024-11-19T22:32:00.696357", "status": "completed" }, "tags": [] }, "source": [ "For uncertainty estimation, we to create a `DataStorageCallback` instance to store the reconstructed image for every iteration of reconstruction" ] }, { "cell_type": "code", "execution_count": 3, "id": "b83a54d6", "metadata": { "execution": { "iopub.execute_input": "2024-11-19T22:32:00.705301Z", "iopub.status.busy": "2024-11-19T22:32:00.705056Z", "iopub.status.idle": "2024-11-19T22:32:50.069910Z", "shell.execute_reply": "2024-11-19T22:32:50.068854Z" }, "papermill": { "duration": 49.370579, "end_time": "2024-11-19T22:32:50.071796", "exception": false, "start_time": "2024-11-19T22:32:00.701217", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "recon_algorithm = OSEM(likelihood)\n", "data_storage_callback = DataStorageCallback(likelihood, torch.clone(recon_algorithm.object_prediction))\n", "recon_OSEM = recon_algorithm(n_iters=2, n_subsets=34, callback=data_storage_callback)" ] }, { "cell_type": "markdown", "id": "bb2d9091", "metadata": {}, "source": [ "We can plot the reconstructed image:" ] }, { "cell_type": "code", "execution_count": 4, "id": "10f76174", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAADfCAYAAADSk77/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB+hklEQVR4nO39abBlyXXfh/4y9z7n3KmmrqoeCj2gGyCaBCCToggSgCiCFKlwyJRE0QrziRptOcxBfO85nvkUYSviCdQnB0MOx3M4RBOkZdJBgSaDCpkSGY6QDfPxiaMIAiBAECABohuNnqprvvM9Z++d6Q+ZK3Nlnn2rC2B1DV1ndVTfe8/ZQ+7cuXJN/7WW8d57VrSiFT0QZO/2AFa0ohXdOVox/IpW9ADRiuFXtKIHiFYMv6IVPUC0YvgVregBohXDr2hFDxCtGH5FK3qAaMXwK1rRA0TtrR5ozC0fuqIVregukPf96x6zkvArWtEDRCuGX9GKHiBaMfyKVvQA0YrhV7SiB4hWDL+iFT1AtGL4Fa3oAaI3hOHX1tb4gR/4fn75l/81L730AoeHe9y4cZXPfvbT/ORPfohv+ZZved1rfOd3fie/+Iv/kpdeeoGjo322t6/x/PN/zG/+5q/z4z/+Y/zAD3w/a2trxTkf/OA/wvv+lv6dOnUqnfeBD3xg6furVy8xm81Gx/bf//f/3dLxP/VT/+xPNmkrWtGdIH+LBM0t/Xv/+/+cf/HFF1/3ev/r//qL/sSJ06PX+PEf/9Atjempp54pzvvgB//xrT6OP3XqoXTeBz7w50eP+Zt/828vjW19fctfv3596dif+qmfvuU5Wv1b/Xsj/t0K3VY0zXvf+15+5Vf+j0IyXrx4kY997OOcOHGC9773m5hOpwD81b/6XXzkI/87f+7PfYDFYpGO/57v+Y/4/u//vvT33t4eH/3o77K9vc3Zs2d597vfxZkzZ25pPB/96Ed54YUvjX6n73kcfd/3/Wd8+MM/W3z21//6/43Tp0/f0v1XtKJ7jm5VIr7e7jKdrvsXXnihOOe//W//v75tZ+mYp556xv/+7/9+ccx//V//aHGdX/qlX07ffeELX/CnT58tvjem9e973zf7D33oJ/xjjz1efFdL+L/7d/+TW9oZj5Pw3nv/7LPvLI797d/+d6PHrST86t/d/ncrdNts+L/1t/4mTz75ZPr73/7bX+O/+C9+mL7PcL8XXniBv/bXvoeu69JnP/RDP8jJkyfT329/+9vS75/61O9z48aNeoPit37rt/j+7/9BXn311ds1/IJeeeWV9Pv3fd9/ln7/2q/9Wr7pm75x6ZgVreh+odvG8H/lr/zl4u9/+k9/bPS4z33uc3zkI/9n+vvEiRN867d+IP2tN4Pv/M7/gB/5kQ/yrne963YN85bo53/+F9jf3wfg7/ydv5XMkB/4gWxq/E//00/f0TGtaEW3g24bw3/91//p4u/f+q3fPvbY+ruv//qvT7//5m/+Vvp9MpnwwQ/+f/j0pz/JjRtX+chH/nf+4T/8r3j66advaUw/9EM/yC/8ws8v/fsH/+D/fdPztre3+bmf+3kAzp07x3/4H343m5ub/I2/8b0A7Ozs8L/8Lz93S2NY0YruJbptTrtz584Vf1+8ePHYY+vvzp/P5/7oj/4Tvud7/qMlx9ypU6f49m//83z7t/95/vE//iD/9J/+GD/8w/+AYRiOvc973vMe3vOe9yx93rav/9gf+tBP8p/+p38PCGr9iRMnkunx4Q//bNIAVrSi+4nuCvDGGFP87VVp/Oeff573ve+bC7W/prZt+c//8/8nP/IjH3zDxvjRj36UT3ziEwB827d9K//wH/6X6buf+In/8Q2774pW9EbSbWP4K1euFH8/+uijxx77yCOPVOdeLf7+oz/6I/7CX/j3edvb3sHf//v/d37u536e1157bek6f//v/8DS5qHpP/6P/x7GtEv/vvu7/9qtPBIf+tBPpt/f+ta3AvA7v/NRfu/3fu+Wzl/Riu41um0M//GPf6L4+33ve++xx9bfffzjHx897rnnnuN/+B9+nO/93r/Jo4++hb/8l/8qBwcH6fuHHnqI8+fP/wlGfXP68Id/lt3d3eKzD33oJ96w+61oRW803TaG/6Vf+uXi7x/8wR8YPe6rvuqr+I7v+Pb0997eHr/6q///9Hct/TX98i//Mr/2a79efKbDfreb9vb2CuecduataEX3I902hv+Zn/nnvPjii+nvb/3WD/Df/Df/hKZp0mdPPvkk/+Jf/DyTySR99mM/9uPs7Oykv3/2Z/85v/iL/5K/+Bf/YnGunK+jAZcuXeLatWu36xFG6UMf+kmuXLnClStX+Gf/7KcKDWNFK7rf6LZ56ReLBd/7vX+LX/mV/yPFrX/4h/9f/I2/8df52Mc+ztbWFu9///vSdwAf+9jH+Ef/qHS8WWv5ru/6K3zXd/0V9vf3+eQnP8Xly5c5ffo0733vNxWw3Z/5mQ/fdEw/9EM/yF/6S985+t2P/ug/4Xd/93df97k+/vGPc/788f6IFa3ofqLbiqX/jd/4Db7jO/59fu7nPsyFCxcAeOyxx0aZ7pd+6Zf523/77zKfz4vPtcd+c3OT97//fcfc6zf54Ad/5KbjOS4sB/DP//mHuQV+X9GK3lR020vR/tqv/Rpve9s7+Ht/7z/hL/2l7+Rrv/bf4+zZsywWC1599VV+/dd/g5/5mQ/zq7/6q6Pnf9d3fTff+q0f4Du+49v5hm/4Bh599BEefvhhptMp165d45Of/BS/8Av/gp/+6f/5pjH4Fa1oRctkvL+1/vCrMtUrWtG9Tasy1Sta0YoKWjH8ilb0ANFKT39T0nHoQ72/u1u4zi1Zeyu6j2jF8PcdjTGzVd+q741N3xtTKnPea4Y/hvnjMT4xfn3cakO432jF8Pc01cw9wthLTB3+NsZi5HhjsaZNf3tcYniPS4ytKR+jv1fHeTeyEaw2gHudVgx/z5EweWRcTGTqzMjhZ2TgyNhG/dTHCVkzSZ8nZk9MP4x8FhjeKyY3agyuqNrrwPdqA4BSG1htBPcKrRj+niHDGJOHDD87ytzWthia8FN/ZwIk2fsBF0M1hvGNAIjM69I/CJtAuIZL51kziT9bnO9xrmdwcwa/KJk/XDX8WGkC9xStGP6u0ZgkD+m7NqXyRlXcWAxNwayBseXczOhGO+YMGO8CQ+PC8V6f36QROHIM13uHoYnHNen4sMFYrJ1gmUADzq3hcQzDgsHNC3MBRFPocW6BmAfjPoHVJnAnaMXwd5yOl+TWTmnMlMbOEkMnhjfN0pW8X1bFZRMIjDYk6e18lzYNa1vF+JGhsUqNH+JnTRqH3my8HzCmoTVTaEJuxGCDtHeFlCfduzdHSQsofAJpA3CsmP6NpxXD3zEKjG4wGDsFAjOLQ62xgdHbZh1r2sQoWtUeY3wtUbV0F/u7uA4uHOOixI6vXzN9rQVYGzQOCKq/83363pkOaydps5Bz0hPH3xum2GFSjCP5B0QD8Qu8H1gx/htLK4Z/wykzOlGKWzOlsdPCLm7slIldp7EzGjOhd3MWziVbOdntBAYV6e69w7lopxtbhNqd6wsNIEl/+cw4Zds3mHiM1hYMNjn8tLYAwQywvqexOQNSXx/IGkpjgbXCCShj790hg7NB7S+cfyub/3bTiuHfEMr2uUh0Y1oauxYYu9lMklGoMW1Q5WlwkakGN0+2cZC0E4yNtrpI+hFekOta24bNwIO1KLvcLh0vm4V3yyG6MUefpqx9BO1ANg3vXdgUaIPdb1paMys2Fec7OndI1x/Qu0OcWyS7v/T+ryT/7aAVw99W0tJcPOxtkuhBZV9j2mwpmznay6YpwmBDlOzOd3iCJDS2yaozNjjlIqNptd7aIJG9HwITe8fYqw6M3tCYNt1fmD957E15fPhZmhXaYejocrxfzBETxtaYCa2dBYdfvK6jyybFYBmiOROeeRFNiAV4s2L820Arhr8tFBndNAWT2yi1xQ4WKQfB9h58F9XknOYri9+5Lv5tMb70yKfPsXgTmFqr9QZLY1qccVg/GdlYgjlg7YTWTJNZ0dowrtrxlscWHX7Fk4fnludqzKy4p77O4DsaP8OY/Jnzemw2OSwhbBj9cMjgFzh3BL5fMf6fkFYM/ycizejTFE5rbGAiYfSlReznUYJ12YGl7WYFnrFMMDaGz1S83YoUNZMlr3iQ/A0NDZ4Bh8MrUI6P2kFjgpkQ/gWpbmNYUPsCtMOwVu51KNDGuRDp7f3AEEOCQmL/W2OT6eLQTkmbtJQQDbBB8mNxMd6/YvyvnFYM/xVRyeiNXWPanjhG+kUPtid5qQe3SOEv7UwDML5J1xmLtYs0zjZ8X2wmHpekq9j61oMzywypkXlyP+ujSeBIjJo3pWEp1i/3Q+7hyZuItQzRFyFM7YgmireFFiD2v+AFtF/BGBu9/NPoxBRVv69s/DCiFR1PK4b/sigzOtgQN7drTJpNJu1mssuT5I6MAAG5ptFpRZgMUcWXY+3F9wkoYws1uAblaJLjhNlNJaNFIxiU+l2H5+Qa8r1sSBabJHjB9AQJLra6Nw58h1P+ANlMXPxpzSTZ+vJcQeOIkODG4nyw7wdB+sW4voB6clgv3GFFy7Ri+FuiSnVXAJmmmao4dZeYWktu+elcX+DWhbE00Ka2nWVT0Cp/YhJK6ZyuqT6vj01jib6D4ZgqKcKAmtmFqb13ydHXAH2MoXsGDA09RL9DjgoYgqahNx2PC+OqUILJcUkX57GJ95rFsa+leQ7RjEWU+gLxXcXzj6MVw78uZWa3do3WbiRHnCDhIKPexpxdNT5dmE7H4DVCzvkcd68pIeawhcTUXvt83/Gaf6KeD75PZkbenMa1AG2Hy3GNmeHokC0ozMGAAwY3D+o4TfynQnfqui4yfcjmawrb3mIZyJGDxoRQpvMtDoezPc51NHYWQpiR8Yfo4AvzvWJ6TSuGvymZoGbHGHrbbDCbnArST1R3XzPC8qKWn6KyW5uTUMRzbyPbON9nJFxUZ8c83hhSzHsJYWcyI2kKTD5kezpqJPEm5bUpk2/q55OYujMBueewDIQNxLmOzoT7tGZWOAZlDnp/lJ7JGWgIYbuGCc4MdFFbYOT+IRLggsPSTOO8rSfcQje0gend0YrpK1ox/CiVMFhrprTNBpN2k2mzBUDvQnltryTQEuzVDzjKENcSRj563MW55tW1cjw8MoY46EQ1rhxuGbY64KMGIPceNMKtwtmn51DXDuNfxukHT31U0WmWknK06u98jytMiQrwYzJDh41qwMVYvvM9Q8QgxAEGcyjNTwxLxm3Nm0mOMMT79UTgzkrFT7Ri+IJM/H+TM9fslNauMWk3mTSbNGYWpUZk+KgOWzNRAJagSg8xA01nomlmh6jC47Aav07uzCMqfhFnr0JzLkp7Gc/gF0qdLr3h2lEoxy+FBJ1KwtE4fjLU1hqLUym0jmW1HyKACEfDpHDC6fmAYKr0DmCe0HdL/gwbtoOWWbHZhXdmac2MppmEzbC3MSJyBCwItZlXTL9i+ETHg2esyWr3sWcrx5OE4gxBjTcpNl4WqRAStbamwOhdYvYwyhIdZ7Dgu6C4Rmy6MaFAhY2P5f2QQoEyVi0JZQwaYy8hMj3GpOrjcN6B6Rh8R+/naax5jNnUEOyBMZam2nDK5+2SWVPn8YfrBD0lGCuTJQ3HRNxBa6b4ZsiYf2/J+foPtrRfMbzGvSsAjbVTpMiEjq9rL7smr5g6XFUt6pH1nb3fQ2GP5/sEH4Go3UmyG1ug48K9sgoexjdgfM5gq4tb1LF9b1yRaCP3H3NCpnCa73Aeej9PabF5RjNeQDv8UliPbP4kLUXuF59bj08/AxAQhOp68n3YABsaOwtahV3E3IQQvx/cUcDq+44HlekfcIYvM9l08QlNstgG32N8hKYuea1F5c2eZmr7vQo75esvM73cY9n+L/0Ezvf04qH22Ykozrsyji/RgbZI3hE1v3bQyXkBhpt9DhKC03NTbyYag5BAQ5Tj11EEq3wOGWWXjxVTSXwI5RwEPL6EBeV61rSYxtIwjdrPOt2wz+AO8G6RnuFBYvwHmOGzB14YHaXqgpJyynuepdUy0wPJE20NDBVwBcOx6iwoySvebIWi09LfRXPBexeZfR5i0srOrzeKxkzj7wH33ihb2vsB17QJEKSmKJ0vpo3GAYinf0z1T5sYJI9+a0McfdDRgTTGJqL0wjidGdGi0uYVNRpyeBHf01pU+K9JG4fWhpphyryz9BzAA5iD/4AyfJNsdWumFHXjKJlLq5SOfqlSs5acybYlq9ZFgYpKYoXrlrHyQgrakESjEWhix+r0UhdjzsWGodTpxACR0RszK+Ld4dldyqxLyTae4tzWzJCkH62qF6q3gGbcIjxCxOUbu56kryFoPxoUpGv0eQZMtLcd5calSWMJAHoJZ3oFTorOVEk9lnt1dko3HOCcVOIZgHHcwpuJHkCGN+E/02LtGo2ZFkwbjtCMWQNFlkExSa1lYEje+5g4ohhapF5ZISZXm9EUEGhtZMLsYXe+D05AuaZa2ClExjIeAEKce2LXaVmLYat5YpaQYTdBJ9vIRiM+A2sm6fjC1tb3kXCcoArjNbyNkp/su8g2eXBzNGZCyxoDHT1HJdxXmQQJwkxZPUd8Hc6U7y8w/SQwfLOVpL81Exa9VZ78N7+kf8AYPqvx1gR4rJYqkNFvhZdZSTLIKrNmXFFtjc/naEaNH8YfehOIVWai/enicUW9ukpl1uWkkkc8FsbQabDlk2uTYMB6Ua27QkXX96vDizp+L/kCcm0dux9zajrf0fuj4j7yeXhui/OTUbOnRgFqX0J9pzFtSrQSg1PvuaVt1tLxg7yiN7mK/wAxfCweaSwSX5cqMkB2MCkEnGboMUrAFR0Hp2Jy1AI35WJPzO5Lp5nY4kvx8bGnMjaVxaqLaiwxsh8YfBclrWPw84SDL8Yp49O2usmFKnRarz6u/l0n9Qy+NIcK6ewdvQXjG3plQmnSuQC6oq5ssKNzk8KCYgbNi3mxZkLTTFWREXizx+wfEIYvY+zShWUMBKKRcOHM5YSWZDdXjrEcVioZdmk01aIeheIqLSDlmyvATjgvVI5t7IzWrAFRatMl1TtcS9B3DucczvQYQuqqvq9+ptrxV48TyDF/+iLmrsOIsnl6XLG5pGOF6Z1oRHmzkgSfoIFZGgIbWgX9HZtb8TnYQvPIkOJaYwsqv4s+E3gzM/0DwPA1s0+TdB9zziWgDUpFVp7pumlDbc9LPFsz3PKIlp1P2m9QMrvKODOlDSuqtKGhiei8ASkwqW3bnNjjTbCr9fOKE1DGbWhKTUSSfnwuW13DbgttRTG7ZAImEJEvMwlR9w7Zea7YPMN9Q6hu8H1i+rFoR5qTiELUuQhynyJZKJoP1rbhGR0jTP/mcuS9yRl+Oa1Vl4SWaqt6cYtXt7UbZVUWhSarJaEwovac19J/jPKGkZk6qPfjzR+E6utjSarw4Of0bl7G5JVDCx+mRUwAuaczHZ2xDG5Z29Dmjv4nx9TAoDDz5d8yjgI1eMy8pOq4uhovxM2jzOYLczRJZoaYRSH1NicQDb4vAEKWNjnvcIR8e9MH/rbgnCPUzFf22puA3sQMX6a11qWhdby9cAoVXl2LU5h5zewJjEK5oOV6x6nzx45WedeFgVItPCWp6g0nLPIOAaUIs4/VpUsxcjXOxkwYkCYVTTF2PRZgafNJYUdTqshybg2ykUiGHLsMKqqkfvSHWCalbyBtKiF2j2kx1qbONymSYXPsPlcCyv6UXN03An1sHJMDb/qgCXnRcN4cTP8mZfgyrbWxa7nMs6ikCgkGJCedrWLKksgRSijnmG+Cs8brlIi6hsY2jMFwR2P85Ji3LPSiXrx8Hj3kArTxOKwLz+WbUhVOm4ZqIlFITSyNGXAmF8IQ3H4eQ9l5RjYfnSSkc+rTfIRf4o8RrL44Oqt5MAQHpPMBeyCMPhYGBfFpzNIcLiB3v8EtpRkn9b3eaLApDNmq8QymDcUzEyrv/mf6NyHDB2aXtFYb1XcNDYW844czLKlgYkJxDTjvElZ8GBaJgcdUVWDpO+9Dqmeh/tMUjJ2vkR1Wci0hgZ3KxhI2nj47vKKEXfJLoGrg+SaX1hJbNoYRa80h+R+UNjumfgf12inmrrz1nqWyVbKx1ddMkRETEmOcKZF4WUIHLIOAiCwTLA0DHYOZ40yX5ynWxBeIsaj+mPHnkTElTdC19IQCmm+W3Po3GcMbjJkUKrxUjxUJpb22SRqSMeeSSy6VZYcKp647tgqNeaaNacLiciSmL3q7QdIYtF0LpTqcykB7cCZLXM1kQYW14FRl22jbTux6cprJWPOGNKRNTJeMSplq3mKipHWuT0g8ZwRw0xQhunqz0b+nwphmJD8gfi+edalyUzjbZI4stEwxZsbEbNAQTBLNh4V/wbi0mduo+qfr6tBiBfKZ2JAOPTRzumGfRb8TUXkd9zPTv4kYPkr2qvFDqgArjJMWVZ8YMTF+rBSDz55ubbMLadu9zirTnvQxsjZ4roXckJFrMs50HQUyMSZ0gvNWIdRG7HStxSwxnYIQF+2qjsEaZCkdYcJDxagjvoql0KEai9j+dWZdrbLr59COUFy8XpNbWzkCtmBQAKLiGdS7rYtzyjPWm5bBhoYZpsXb9agZ9HReHHn3r6R/EzF8gKJKnN2aSVLNNEKriQ474xs650B5dZ3rCkePSF9di67wA1QbAUBtMtSbgTVtkLpY5jGxQwpkpKwydQ/xKUicvbEzBjtnMewzuAX9cFiBXUofRO/nWFVQUjsrj6MacJRwBSyWjstazaT4vTjO5GIdg7XBjKgwCLK5DJWppdtlpfckvhV/gKGh90dl5ARb+AgyDiE06CDes0zJ7ZMGJw69idnAmobWzhKMuMPh3Zz7NUb/JmF4E73x0qgxVJPVDSBSGSmVKWZNG8sq9UuMUteQK0JSCtAxrkoOSQ0tqtSaPIZw/wnW9thYflnb+PnYCTNzgg1zmjW/xWA6unbOXnOVo+E6R8bSD4dLfgUZx+CHYAp4OyrJj5Pux2koWiOwpi1scnmmRhyFWiUfKfCR520EBizHe5KjTUt9R5eSZaQclnYeyjXGNLj0/n2O6d/Mpje+iWHcGb1rMabnfgXmvAkYvlTl22aNtlmjsbPCftUSXqqotmYGFqyPXtzKyy0eYh0WKmLNVV05Ie38kqyxWs1OjB0dazrxRdPEbnDCnOdx9xjnZjOsgcXgebV7jIvty1w3L3DINSWFXeXxD1pMnfsOJWimCJkpp5oGBCUIauUPEJ9GY2dJFQawMgfRF1KCecrKNzpn3okjUPLjTZu67+hxGJU+nEwvsuqekmySgzNI+ZBZ2FT3bZXaL4U257GaT5fXgplGR+T9Ccy5zxm+9MiLk05LY6GsMoYkCiHJEBNyFd79OFtVw12TSq4camPQWxnHUuisChXJ9VszZcOc4YJ7lHef3uAdJx0bjedwsHx25xSf257xuSYs5PmwO2rPL4XL1FwUcemKdKss7QSTsKCcJw0yZ81J1pszzNiiZULHnAUHzN0euIOle2VzoCmYPd9HRzeyViQbk3NdsSGNFyUpMxHFiTf4HqPew9KY4kYh2lnv50mlD3b99L4F5tzHDF9lvtk1FX4LL3qIlWWFcry77n6q/q4WpfojfSZqayhEkWvW6SjAWLwXwsLs/Ryr7GRr2iyJ0iYwYdps8bC7wDtObPDnHznivW97hY1He/o9w2///lvYbDc4uvJWujZkoUnhR+3llyeoGf84ZtcSu23WE1ZfMgqd71kMe0nKt3adjfYs53maC/4cJ9qWxhhudB2X/TbX7EUOjcUNO0ulvupKOBLX14g4mbcmgY9yjYExH0p6jd4tvWeZf88QsPtjzG60JtGnYh3pXdk2tNQ2U5zlvovR36cMH1F0KtZel2AK1tVyrDz+Ea+SVVeRLpKHHjzTy8wvseLUrFHCUT5LogQ2sa7IEBMJiSPlbBdFMdR4WjNl05zj8fYk/94Zzze/60VO/N2vwT/+KLPL1/im//GTbP/eU7x6sMH24VvomkNcVJ1r0sxeM0nhm/A2JcSIg9GaCa1ZS2OTkJxUkllvznDBvY0/c/Ih3nt24LH1YMJ8anudT18/zx8dzHjVQm/mdErl1kU4k1oOOHe4tBl5P9C7uXKUDsU7XTKX1Oaq57Y22YREi0kAHSNryKX7Zd9Ak4p64GC4z9B49ynDWyTe3tq1Qo0Hjt35j/Oah2yuvAgbmohAKz3GcrxORRWJmNBfWFBON2dKmKvEs42SxMUmEa85bbY45c/y1ImW95y9zonveZrhm98fHuStz7B15QbfcPEyf7j7Fi4enmHbnKEzBzjTLzXI0M+tk3oSLkAh58bMl/BceWOa2HDuxKxz3j3Bs+tn+O7HD/nAX7uC/donwHm+7qef51/9wVMsXjvB3vwsR3aH3s+Tf8HGXvGpiw4lbFnH9Qe3SBl+x42vls5hHUSglPYVaDNHvPLkOghCuuvtcuhxkoRLGN8RwaC/95149yHDB3BNaOK4Qdusq28qL3Vl09WL3vgG02SUXThm2Qnjo8o+5rkOqrlNvdYlxm59sOy9GZYWXB1+EySdNF9s7Yx1c4Yz7iTnZvDI6V38295d3Ne952t55Ff+FY++6DgxmTDt15f8BIUPY0QFHlN50/PKOSaWtDIBeSh2bWtnNEw4ySaPb1r+zLOv4P8ff4ehCUvqjPsV3vNPrvMH2+eYHE2SxuRxNLRM7AatmSnHWpldOIYa1HMYHzTPqS4Rrjdq7WCVTdpnaK/gAzRc2KlN/jhobzB7AGJasoMAyrm3nXj3GcMHu72xa0zardCiWbV90hI+HK1ivJGKRZ+KNMbCDtWOXkuD7HAbinOFbYq8dDqs72JILG8Yteoo9qv0crdmwtRuseVPc6KZstF4rPWwKBnCnztP+8iUzcYxtZbWt4VHXcasPe43s3eX5oYcnhyYM/iy+KTHsdacYmZazkw9m++cQJOXk3vfN3DhkX/J5pfOhXBYwYRNKGrJDEtDx1H4XhBxlDXzj6veU8B2tZdfYMhVTYHGtMFhF00wTNbyElgqtrx2KkLgBQBUO3RjzoTQ4GT+7l0pPx58vGcp2O1ts8Gk2SzUsEIaiA1qJ2H3jqGiUe+9tglVzLb4W33ufLdUq066m4I4lUYKM4gksbNQVFE1tkgtnc2Eid1gZrZomHDkBl48MHzipUfofuY3sH/8+XS95pf+Da/8+oTn9luuLzrm5jDfx2QoscSPGztN2ASpT1drBDJ3yTtfSLVc2so5QbYNdH5grzdc/6jDfuYz4eBugf+Jf8XvffFRXt537JltBj8f1Z6cOOE0c1ahz1BAc0prpqEmX5xDaS2lx1nWAchNLeS91STH5IKjSsLHDUSKYJYdg/J7DiHJNaxdi5vuTUAHd5nuIwkfwTUCrImgGs2wxdEqlAPZgRdSHpvinQiQ47ia6SIpZZeHsLCs8vhLjFcKUIRr5rEJA7Ypu6tcXA1NdJDNmLGBw7HtD/j8dsPg11j8n0/x3i/8Fme/8bexmxOe+yXLv3npAn9wfeBVf41FZPgmevwtGQRjTBNquvms1taU4tcsb6DZgTUUm6HzPQfMefVgg4/88RO8/7/6JGvrv8v+/ozffe0Cv3N1wucPttk1l1MvPh0iG0yIb/d+nuZLQ35lDFLMIvlLYo+5BIE2yiEaqU4jrkVb2uRMOHawFsNyA84aNShJPVpLtGZC087wvQu96r3nXlXt7yOGtxjlkYfc47wGXOQzlBMG6TM+jqe+mRamvda6WONAhuYCWV0kx5Iz6KdRoaXAPFLySa7ZxM3AMbAwhyzMIYfuiO3rJ9lZbPK5vbfy1HMdM+v5xI0pn7ne8/xwmSvmZTofJbxpAohESSYgOsYsgy0Z/rgkIqElM0d9v/AHXLdXeX5vjcGv84e7b8F52O3g4kHPy90NXrUvsBj2oi3cFOdadIGRLm3e2sSINy2YXcJl+pgaHVk7aJOzlFq7U9766LyVzb0GIJka1KO0osZOae06zi7uadX+PmF4kyS2VKlJC6RSR2tnFeRFG640DnRJjq7KE5zMAPWi8/HjxR/TfVUsWaSnMTaqxLl6jjE2YvjDgu3MPJRq9nN2gGt2xtWDR/nj/VOcsDOMgdeGq1yyL3HIdfphvrTRjUFZxd+g58NTpvAeNycl0jDiHPycHX+JFy1cPdji0/tTDsweB2aXud9jwQF9f1Co0hLyEr9AH4tWjN07/hK1qbwpjZ1Tt7bW66BwWJplppd3FbIDl0O1ENdNNMs8LmVB6jlpmimNW0uC6F6U8vcBw2forDFtCoHVHVnHfpeFpj3AulRUId2qBSLXEulcM3qd8ZWuORbrtqXmIckfusEjQGcCvLO3R3H8gREX7LGwe1wza7TMwMOR2WbR76VnDCGuMn48UOIA0hxITF050Uz1/Lr+ffKCx3kQrPzgOxZujxu2Y4cJ3jg6d8DC7acMtBpRJ+PyUbsZaxpZk/M9fRy3zF0pgUtEnTEWY7NpIPBmOSZt3DbXIEzfK0cqRMShCc4852M5LGsZmJdJQNF30jbrDH4Rn+/ek/L3AcMHIIx0ctWL8riwTfg+53qXsXRl25s2VY+Ra2pv+lioJ9VO0/XWRkJ2etMIT2GXGCp5/AkLr+cQb6UbTFltddHvF9fXDCUlmjAjxStrcIoZZ6w6jq3nTUBFwuwSGcEH+7sfchmw3h3SD4f57cWsxTiIpBXVacdjmlRW7aVrbDk2jXvXjtPU6ddA7xcjMOd4TqW25/dUAoLC0Bsa0+DtLPiAAOMtg1uouW0ivLsNRTOw3GtS/h5neBMncY222Siy3/IRluPsTyg90OJw0xJuOTZdXT+BP6SIpUqXpUyeqYEa2eMdveI0KREjHK9y2+lSEYvQBz0AfCQnX29Mheki7aDJ+PA89tyYIiWM+HqsuXFFwqvHRS/48XI+msjwQ9FySswT6XFXAKFknqOk1XOuK+Lm4h1RGzDabCoZxxzDmHV3GpHC6X7G3fR6+fe4RmKXXJH0Tax3aKOwkDkRv1AGcwW8/b0m5e9hhs+qfIi7b1ZqdXbOLJ2ppGP4eXz8ufYByM9Cjfcl02XJCj7afWPMriUW5CyvOvW2GI8ps7q8d8u15nzZSknO0+sqbYSFf2NIG8Rxm5xuN4UlORWPSyFdYnav7FrdxnqEyceeXzZlyVXw6TnKsJ08o5bM4TulwRWO2ZBHIfMpY6+fvxgLkeFxWN8V6MDWNKHunlkO9YWIQgv3oJS/hxk+qPKN3WA6OcFaexoIdmM/HBY2tJAwaircCEUfNojqoSk3i8JRowo1hmu4VAhDasmFg/MCFYBG8TmZ6SXsJwUfB7dYKtestQU5XjOTeP3DtVV/O0csEtGNqOVjqMEhMWd6ZtriPCmpZc00eshLvL94ypNXnWXVXLLtShxD1ozGtKSSkRva6BcYfIeGHo+R5LlbkwFIS449Bry3RZmtrGGFHIfeQl3eOqT4gjVhs29M7PQzgmMQzcuaFmda8D33UkbdPcrwueFjY6dMmk2mdiuqkENaVLdCYxK8vFPpmQ3hM8UAarGOZ1fZwrcgzqBagtfMXveyk7JZdfy7uKfyJaT1U+EJimev7HCZD+d6pY7aiO23BY4ASvU4zWUyg5b71xfPE+fEUdnPqYNsWd+vGLd6ZgHWONVgo7hfdORpj7r0E6ivlebEj8xpdATLPEiUo44GybMViVljOBDx/Rh7T8Xl71GGt2As1k5jiuYstiqu0Fgju75uMJhV8NzZRcfDLZbeEtIdyZ5WW8TL5TtJcgl/S8GH1NRCATTGGkEAS446kYKC4tJSRc5roh1uTZu0Cy2dtZQJz7/8M2ka5FyC0Q4wRgGMCBENS24T3Y9k4ontWlelkbLf8jqd65c2tBC/zvX9tKkx+C4lKWlUoFWbWNp0K21PMP/9CLpPpLyUPxNG1/OlzaPjojg1ok+/j/wvRJbupbz5e5DhVTtnJTmlUOGgwji64aOoqt6H5AiRZInxDOjHFdy08XHH9i4mUUzUrq+83CZWvjG5XntW/0tVEht6lTtKp1deoDlnvrEzpnZTbTJDjPG6FPsV73sdZ9Y0xugybn3/IpdcFrLLjCl2bhhLsONdtbGWiSxNdvAV+QJ5A0vXXnrTr79ph+NG/DTJmx8AMRm5SPpbYLV6jooNzUzCO9LMrarc1hBr0aYEYl1vKEkbsioSY1q86e8ZKX8PMjxIh1dZDL2fB/sq1ogXySDllAY3p2eBV7ZxCsnVDjmUI81n+1ns9ox26/PmYmxuS1Wp2Pm6TZJIAFjonILWygalwB1Shmtqt4KzzDt6jnCmp6H0GIsVLdevG1bq5y6mUpkeepOQsJo8qzxWiiR4l2rt0aiwnjC04mCHw1irHHZVf3Ym1FVyk/9jJLqS3putzJLq/pJdmJJ8fE9dlLRwFqa1UIYf9c9ifqg2IrL2oTUlg005EtZMcF4y7+LGavp7RsrfYwwfWzorpvQ+1kMzLoXGdAilWMRiax/zEjWlXT46BwWvLQwzKI+6hKJALaSkHpfqZGuyeu/NQM+ijC0fEzoUqKjzy4s7bEjBLtUq45j3e2zhJiZLWo8sQLHv+/CbA8gNLnBxAzBd0qZMDE1pkkw/rfKGy5VSO5tHSu0dE/36nBEshUBgw99lBxzI4UitnekUaD03erya9AYnkl1jA2qzKDkqk7nocitq3zPcI1L+HmP4qNBHCQ9ZNdJhH0HbaRx6XfQBKBxCqeMLGdiSstTiBpK84uQMKqGMLgse8mFYRHV4gbUtg2mZNiGs1RDNApvV8xT7N3mRSiWXML6s2g6+S6o8MY1z8KXGEo4tu8yGkl+lY8rR4V2ozCuL1LlF3Khko7DBmywlnGUDMVZ50zPopmbs8L5yyyfBuuvIQi2ltURM/gFTS99KI4jz412E+9oyFFm//4JZY8JUGQ1ZroOXzhWfgQI1CeV5zBGb+tyGCd6u4xtlWvijuy7l7zGGD9LW1IUTlfqaQ2ER5VQ5gYRB6jCXfsnOZHVWp6iWdcpVZZhKRU3xXB/sP+e7WC++I+S0r9MSHF3BAbTcJEFLoTFnmIt52TrMF/wUZcWa46iQXnF9hU1MHFR681AS2eeIRWHrFwyrwmpKA5E+7jLWnn6JYYNnfpLMscZMQiaftIh6HQbUz6dDkRpxqDdzfb72Lyytm3r+lAMySfDCM5+r5ejQa31PQd5Z0zLcAx77e4jhtbNuGUabjpKFVdnn+h9VAYXCA0swDcZKGo02RUSlTxZMkJ2E4rDLrYgbWmYMpmMwXWLSMdVR1Nb6e2snwdTwGRGoJUs9J/mCpRofPhOmy8yej9caQzSnjMXHvnCy4VrTKqfVsucbpWIb9AaVIbRajde94TCEZBrcEqNC2Tyj9gGkyErSoCoHW7VG6iiMrn5TPlPeKCXrMEV6/DwdkzYYHL1fFKWzw9gnsUFmj/OLu27L30MML0CbtRTqGo+1q6oqYz3ZyBIq7fRGLZT4IqUGeUr/rJhEL1JxDKFi2bX0k2aQvZvT23lUw4PE1x72OoWzxupDdiwGR96wdM/6mdN4VZ84ccpl6Zp/MrLxhIu4NFcpnGRko2np/QKLS51sxKY1xkYfhPgh+lgkpALd1P3waMB0OdWV5efLz6krC1Mwqc6XD6XFchltnekn1ykaZdicl5Aur3wixjRItl7w5VisWU+RHTHBchGNXBlZ1qK1LY2fMjjx2D/wEj53jgmNJNZp7Gx0txfSXtwSANMUTE/hC1jGhst3GkhyM6ffzRal8z29n3PktqO3tszzvlmfumSOQCEJgzAYH195b10EUoej9LNkZvcjEsZgItPn+DTeYpxlYFGouYUNrha+1qJkrLobbdAYOgZhOuleWzjY3JKUTN9hi1Ca3ggk98BGZF7dWKN41ogmlPz345Kfgpd9CK2krQvMTqjk68yAo8M4GyvyqtoM3uUW5EbyKaah77wfuFsFL+8Rhg9OupAks56AKJgM5lhKeDEuAWbqMsP5qmUixuBt8TkEiStRXy0N0iL1ywumMB88xWJxPjQdNGZe+AMKrSAujlTtNjqrxmxPodojn8KBxiYVX9vovmL44lqvt9CSmhs2CAnbBWenTeaLb9YRGKrUk5dIgNY4lp7DhLFKBmEqU2UyCEeSVaSqbRhNh/OOwc8LqTrEVOnBy0bl0jW8cYUEF/NszPtRbzRhrkptDkLVXkuDXMWa8L519eKs7ke4dywb5kyLYXHX3Hb3AMOLdF+jteu0dr2EtiowB8T2RRGXnpwwsdUwZNW9zqZK6pXJveVCXTbA9/l+vtzdJRMqdZ2V+yrknTABkBxJdaOD47zrUKp/qOQWR5c2uPoa9cLMzB7t9DrXQHnfvxwKjEms7qPi4DFG39hp6rsu+f21z6I0O/LzSMHO1syYEGrfW4KHe+bXmfgJLQ0OT0/Pwhwxt4ccssPc7bFwe6F3XtSkJBqihYRW253vg6nl5sGCNGWCUflO1Jh9CeLxfoAmVBm20WEZ7nU8Badii/FBkHGXylrfAwwfpHtrQ0+4iV0vUWco7zwk9U8YDaA10yQdRVrLDqv9I9pRFDYFx+APC2edNS2mCVITk9FwrQk+BUHJeQJCTWeH5QjCSH04oyrg+GaE6SMDx2qrOo02zVTspKoXaEITFscuq+3G15/p4236ziRvV2AzEyG9WYvJDj2hMt493iNPP6M1LWvNaU7aR3nUvYUzZp3NtqGxwayYNYaN1rDWwKwJqMX5ADcWjhuLgcv9AVfsJW6Ylzn0N+hjx51s3oU0YGMnWB9DmkpSW9vifBfWQY2mNFUiUeXXkc3CujltfEw3tmGIVlJFVYKUt+DNXZHyd5nhs2e+bdaZxiQZS8NAsMMyICVPXCnpmuSwEWktarHkiNdJIOH6pZMrjEYDKHKDgtbOmNiN4OQzXbrOwLzomZ6fyqaFVyTn+ABBFaz2zcAnMi4dFkJh1rV9PLjo/R2171WsvfpsjOpNQpjeK4dbGEuP8zYgB1UoL4SqskNL5li0oKaZMm02OW0v8Fb/FO8+s84TG57zs4HGeAZvmFrHRuPYaHum1jF4w0Hfcnk+4eK85fndk3xxb8YXjWVougTIGvwia2MxRFfX5JPfg+Doyo2CMttSf57mO+YEONMFE9GEtakrGeu06jyvTm0obVwDd17K330Jbyyh6+s6rd3IKC4fdk7d3wuWva5ph/b595ROqrqVBOjsEFW5RjG7blpokzTXJNVkw70nqZOpx2Grem0SDxYTQIf/BjcPyTrRzq0lodSis5QY73DtJjsCKbvTZifdzZTK4767CdINryScxjKgHHDB2547uERmSSE9Ei7C2pZJs8mp9gnead7Gex6e8F1PXOapr77O7K1reOcZLi8YDj1egGnOYBqPsTDMDfvXp3zm1fN84sY6H738FH/oWy5OPsdhdy2bFPTBv2NJztLiHVUmBpDem9YWBSQ1VhnHEbrXpjVURAIUhqS6t41NVMTsutNS/i4zfFDnbfpXhpgGP0+ZZ6W0K2G1vZunJgPJ/lUMVSDufEPPUek9lu+MDZl5KoPLeZfGZWlwQMOEfgTAkXqxRYx8y1pk/DDeudnDuL2iaQbk8Fy4R8X0MjbZuBD/w0J55HPNtnzNr3QpLav64TNd2JOk1osUhzYxjcGmjUlrJm2zzkZzlsfcU7zz3IRvOb/Ls9/dwTe/F//Wt2L295m+8CJcuoa/sY/fOcIf9pi1BnNihtlcY/Oo4+TvvMSZT59nPpxhfuUxDv0+QzNnHjfSsc00jFvBZY+hXC0nS+vknK0dejEM1/tFAmqNrc/kU4jmYkYj9ndcyt9VhjeYALKxuSii83vJCSMTWQMt9AsTqKzzOQYsKbHeZ0+4FL3sHRjTJa+/qHApbo9Dt1QK4aYA6MhosmEpE8uaNiTCNFusmzNs+dNs+C3W/ZQGS8fAVXOVbfsavZ3jhnJhhtDPcuup2o7PeIEhx9vR6qOOs7+ek65e+LeiIYidDhLnFylvGoWBUJJREGez5gTneIpn1k7yDWcWvPOZ1zBvfwZ/+jTM55i9Pdg/hK4H56Ab8Ac9LOJGs7mGOXuCtfd43rl+hf3fbPGc4OjS0wxtxzYvMmc3wZ7H/Ak5whIjIiqmLziKGrAla6T2VTgVZnS+L8LCxtjUZkyHMtN79C7WsL+zUv4uMrxBsuLESTe4efJKi6qdPaPLalWa+Cj9NTMc5zwS9T8wVNx1Y632YKuFEsqSkCEvuWdeYL+lGIf3AWnX2BkbzTlO8QgPu/Ocm6xxatqw3hpmDRwN8KW9NV5ya8yb3eBoin3OvXcBF24Dkxs7SQkptZdex/Plsy+fjpdwo+RDbD6QYnyCWhrkf5uclRa7FN+e2HU27Tke9Wd5csvyjtPbbD5jYH0G+/uYa9cxV67hX7kGe3N8N+D2etyhw7Q+lBM7XASmf/gUk3fD17x4mZ1uwqsHa+ztPsW82YtCYrxAJojEzm2l0rzKuqCjhixr27yOvCTgVQw36iSdYM7lrMeJDX0QPY7BLOKmYLmTQJy7yPBKnY+OqCFNXpkhBRqgIn3P7RIiOeW4kxMlJP6ps72AlMpqcKHKS9xmg1lQaRaRhuhchFzMQuz+reZhnvRfzdvXT/L2k5anNgbOz444Oek4MV2w3034nWsn+cTVR9ibbzO3u/TDIRJ+w5XqoIRxdBhycPNQMUf1LNee8RzJEDXxy6Hx45Pn3i9rGtBHYJPAb2P0w0jThgZLkJBrzRnOu8d569Yaz2wOnDl5gJlauHIDs38Ih3P81T2GV/fxncPYwI123WI2Wuxm9O3sxoq4/cCJpwf+1P4VXjq8wEF/ku2jCyzsHoOdF4U0i7CuhlRX/efgy9tAaxyG1iTF9NK+nNbOoonYxoKdd77u3V234YWRNTCiVlO1MywVhSDaW7Uk9xQvOfVXU5qEXFPuISp0rVnUFWqxUMRiCS90vT3NOZ7k2Y2TfMNZeNfJPR4/vcOJrTlrp3pmb7H4wbP1mx1wlhdeegvX7Ass2AWptjKWcRWRY5DbF6cyWTHFVeZRM6zBRK/8rSzer0RDCOcFz30Ib+o6cUByvjqCdD9pHuZhc5KnTxie2pwz2+zwvcVf3YPZHLoed/2IYXvAdx7TgF0z2M0Ws97CxgScw+8ewmLAD47m7IyzTx3w7NUDXjva4oXDc2zb15ibHeq+BaKWF3n8yimnqWZ62ThSuNXkxCwNwipANyZrO/p4idnrNervoE5/lxjepEkINnEp1WuAiTjDjMkTCIRqNZXk1tcQUI7E9lszSw600NklNHzQWXJQ2slpLDFDL1Vxjf9aM+W0eQvvaB/hWx/x/IW3v8ipJxbYLYtfeOxmQ/u+p/EPnearz36ayf92iU9eO88Xjx7mqLmBhn+ONUcMDsMJgtyTfHbv+1uQRmNMf2vhufy9LePzQknNd9V86XTaCRbL1G5xzj3MkyemvPvkEW976DpN63H7A3Z3jnEeWotpDcaC62GYQzN4sAPNxhCcL73DH3a4neCvsaenTC9MeMfLV7m2mPJH2xtc7M5yYK8UZo8eTxp+VZhEfwagY/Fp5nSprWLjiPF1yuIbjZ3hcDSCTIzjkWpG4VpynTuTUHOXGD48qKTB6sICx5FIPM0M2nuq8801Ak4ccmHxqTryPtvEtfmQrh83I001BqCxMzbdSR5eb3jnqeuc+5YWc/ok7M8ZXjvAtAZ/7iH8M09j+4EnX/wETz7/CCcPz3PdvkA/HKVnSd7daE8OdLF8VpYK4QCVCHPcfGGiaj/G9Deb6yoz5abkyKCdZehvMLFaJmadU2aNxzYMX3XmBg+/bR8zCcztjwaYDJjWYmYtdrMLIbVDj3fgjjxmL6D9/MLhDweGvQiS2vKYtYYTDx9x4cohp6cn2VhsBcbEZqeqyeMR4eIqwRLG3BQbgZbassa0lpjgs7JpOIr1In6hoYoGyXrLDr07B8K5ayq9rlk3ppKXOdRNsaMmO5xJqrKqARa6LRLEpBaOik20cwf0fl4UutDhE2wuvKgdZBp1JiOaMmGzhUdO7WGffUs49ouXcPsO9h3N4SFubQ337NuZ/ZmXePIjA4/4h7jYnKR3h9nejFdMCyTya2vWEhgoz9XxSTwCoy2ZPpzz+sz+5ZFsPEbwD8IEZpaiGmtscXLScmFt4Il37jB7z8P47UPc7iIw/GGPWWthraV5bBO7tWDYXuAOPL7zDNf7oOoPHr8A14XYfLPXYdYa2g04s37EyekpNvxGEVaVMYqvpK4NCMvOOMmos2ZS7H1SvxCkxuI8lUHTzt/iuiPh37tJd4HhBV1XOlOAxKigADWRQhYUS5pPwktT7pzhepmBvB9y1xFC0YlhpClh7QMY44EkyaJNP7dzDnrYP5rC7gFsrsH6FLsZz3c+2J8bW3D+FOemC05PpqwNJzmyN+Ks1CZJrsaDFaRW1FCKcWUp+/pUe9pr+pP1Na+RfoI1n/gZG63l9LRj8pY1uHAO017H9DuB4fscdmM2wZyYw2wfczXY9G4ObvAMc4uLWnc7Jdj6U49dM8wmPRMLLTmvQZtj4anFN6OLm2hvvZbuDUvJWFpLiOaLrLdQoSiDcvR7HOufp1F83EFP/d2T8GPMjk2JM6J26bDcECVJKkmlJLNmEonXJwjlSMukokRRGkOEhqo01Tze7PHX6mvv5mw3V3lp/xyfuvoQT/zb55n96TOYc6do3m4Cs8+mmN1daFtY9Gy0AyenUzYOT3FgN+n0QsAVbaWc72k1Witlj2lVvd603J8AeLM0U4yXnBWzzKZNqAQhhTQYa1omfsasgZl1MLUwm8L6DLPWYuwCrIFpC6c24PRJmC+waxNMuwP+AH99YDiwdAcNgzO0rYPpEGz+tQazcEwmQ9qujs1TYAR5h112/MbPrbG0rAFEqPeQQrI55z9nxCWG9l0IqFY96sVBLWMyMSxN8T7fWOX+7jC8WiA16imXK5KQlC3i4YGJLTobrN5JLRS79lie+2isXjLlVHNICfEJszd2VpgbHsfc73G1P+Rzuyd512fO8PYzO7SntzBntgLD39jBbn8GdvZY/LtX+NLBM+x2DhcLP0pbJZE+OkHGuT50Kh0BkEgiSyYNurm9lBx3S9iGqA2RtaLwLA4XgScOFxJgnA1AmsOjAKyxSqMwBqyFtoGhgaaBjSn2TE/LAglb2bmlmTkmJz3Nw+twah1r95hM4gapGEaDbDKgqSk96tqeroVQ5fMB6P1R8v2Ee4yXGSvi8xozUa2ttGHeoWSau+a0kxRTKCdNdswm2k+DDxlq0hNcVzJdCoUkL2ivssNcwUBj5auPqw0nGw+0qQqPNGXQlWsW/oCr5iqf39ni311+iNOfOuDhs1exF06H8fzRRY7++IiXnzvF7199it+6bHlhcY09c7XQYGqJIIshAYtEUpg2eOmxGFXzLiQKMIK0+zKBNrdEtpBSshnqLLGeMP8Lc8h+77ixaOheWdBc2U7M7p3HOMB7WPQBaXd4FMJvgDmzQbM+wZ5YYC/PGfZ72hOG5rENzLNvgRObwIu002vxiT26uq81La2ZJvW8BbwJRT1rZgybQpucfDXp3A4onca60KXTUSfvksYhfQiTA9GGUuDcQbX+LjB8qdLoDiVBXcrFDKQSq8GlycXkTaFm9vrl3YzyxOfEj7Hzi6y5iLMfYmPBkG0VMP/79gYvHz7Mp25scuJLj/F1/9s1Tp28DMDLl0/xxd1H+ezulOd2PZ/Zv85r9ksshr2yA6sC84hGIePRXWsSkxFU/pLp7ywtq6iZ4SFssHNzwM5i4MpiysGlltnVPcyZTXAe3zkYDBx1cDiH1oafXQ/WYmYTmE1gbULjPHatp3l4HfPW87hn3wZNg/3CS8wPWxYOBoZga0vc3U6KEDBAY10o6FExfXygMG7v8GZQ5bdCyrIIntoed2qjrTft5YKauWGF1AscRvIh3gi6KxI+A2naBIYoCzUGZJkGyzRmkhBxKeYu7YSLwo1xcj2F5M7AibJccbLLfQkeKcImqd9cTPE0obChSPrezTkw13nZvoK78RjX5xt89NqjKZf7ypHn0lHHq8NlrtnX2OUSi35vKUIg907aRCzz1fsFbigXkTVtgGW6mzH9GyDZU/xdptlhZExq44weGKy3zM0eV/0Brx1OuXx1i5Ov7dJsxrLWC493Hru/gLUJprWwGIIphAuq/9oEM22wBwt8azBPnsN/1Vvxj12AboHfPuTq9ib7PSwiNLY101R4NM2GCZDfJjap1Fpf/CWp4mEN5l57em0mh5vSTLWfSL8nyKi7tO7jmjfG4poeM8SUWT/wRiPu7jDDmwjgWLbbpb6YpBo63+NMizUqwaSKiULIeQ/45+xVH2tHXIIwSmbX7ZIK4I+yt+Qa4rzRnzk6enfAtrlIbzt2Ds/yxcMAwe3o2THb7Jsb7HOFRbefogNasxGGLxBc2PJeVW+9oPHcGo7h1ug4L/7NIwF1GSjJ8pOy1b2fs2f2uDo/yZd2TnDhhW223tIlsI3vk5c02O7WBYY/6vFHPXQDNBbfR7t/bQJNg7lxHfPiSxx+cpfnd55iZ+EYot+g9rCH9+QKjTE0tHBFTYNkWjVtqECrNRikDFkZMh7DcmifQI1ClPMaGoYUan0jzK5luvv58JEEfw0y6T3O9Fi7nBbbkCW/MQ29y2r4cQu/bk+cQ29NkvDWTpIkTbu4FUCPA9WNJqXiBrQFnTsMzil7xJ65hKGh90chBz5i82tVL9l8qgaceG5E4xB7cHm+SmfS7aUxsI5y3IVRU0Br4zsLCJnIIDKnvmPf3uDS0Tk+v7fGky+d4qu2jzDnN7GbTcgjWJ9gNqZwYgPsEWwf4A97/GLAHHaYWRvK3rTx/nt72OdfZPj4F/nsZx7ms7szrs8XDHQxOjBR76pLfe/zE8YEn6hp6oiNCBw5LgB5GqyAusxQ4DGKtlK+rFev+w8ITiHX82vS+ns9E/R20Z1n+MrWK9RwsuRKcc4h27QSPnEmLDK9i2oHSlpwWjNQoTYdHpHvtbpVl29KAAplpxUYAd8nD3/ofRfhu9K9VmzvY5yDZdSgSWqlrgUg5kg5b3FsX/ZL0IvrOMfe6y/AItwVY9PO5e8C0wet7dBvc9Xv8vzeGp+7cYqnvnSD2Ykp9rGtYKef2YKHTuPPnMJcugJXtvONWgsbQa0H4GCO+eIr9J98hYsfnfKJ6yd5btez7Y7obGjt1RAYfEgm04BU/ZH6d0VpK82kBLvbAZgW3cxSPP0SJdLCQvwCGq5ba2kyFmccLRlMFdbyG++pv8MMH2O3KvzhxdZTu6WoT6LeWybBuSGJNjG8oWksoWHM6aXVeC05U2zfZscLlCq9XKv+PSO5+mKjyXn5y6m9xfkj5ZSlSIbGC+jNyvke726HvVfHn/PE1rH8Ek9fPo9IeBM3yxDfDu+k9ws6f8B1e5nnd0/z8PqMr/3CJhfO7NN81UNw/hScPol/+HzIjV90MLkIUxv8A5tTzIn1gKWfd/hL2wyv7vPSb6/xsdfO8+lty5f2j9gxe/R0o+8sbNbB1JAwbw2jTc+Znil3sdXNMOScMAvBwy4tyHT+e/IbmaZcK2I2Grn2RPlF3lhM/R2X8OJhzmiouPMKQ1BO2nHXGP1cvbjC+2+CV9abnLsOoOu8SSNAY22unoNbyrAaKsmQQoGFM0eSgsp2TDXTy/ni1KmfLUj1zORhE4uSI5ZktiYm0rzOvH8ltJQwAxXgRsa1DJGWBS02cWcOOTDXedXt8Ic3zvBrLz/K1x1e562vXGHt6Rs0bz8Xlvp68H2YaQuzFhoXvPTW4K/v07+0z+4XLC+9dpqPXz3FH+42vLC74LLf5sju0xPh0qqqTGoJplGckfHqohUy9rHcCoghZGmJ5XOvwTxny3iJm/lX7Ous9dtNd5Dhs8NuSaU3uZY4aHWwKf8Wm5s8QbomWWak8gVqVV/w0AapVZ9jqk38zjPQuUO8H1QFXcEAjNi2ymEVGLhNkQC59vJzWJmW11lguaimzhTsHcG/4SSGq+zuaD+HcY9kusm4j9sm0gJW10zeecXskeFHcyLIzy/owbnZ4UrzKl84bJldPslr87O8c2+Lt3/pBm956VXWru5hD49CbvzuIf6wg95jOMRvHzL//B6XPr/Jx187z2d3pjy363j18IhLXGfPbtMxxzGETEg3X4JP1xVoxcaWOL0wp6yL13OEije/jgbJd2OJXvJ7rVEEn84bH4u/KxK+Jm1flUkzJaMIuegoknOFxjyzck2IUl8znM9Y/GUUVHCqaSddUc5a4rsqbl6DebRdKNcUBpYNBqBoVaQWpRRVrDu21ptc3kBCMQqtLhrFqAXz3iKZaKNSvIvA7CmdVz1TfhelI8p7R+cO2TWXwMKwd4GLh1t8YXfGo9cf5a0Xz/PMJw945vwXmK31eGfoe0vXNTjnOVxMeGHnCT63t8Znt+FLe0dc8bvs2R0WHOLIjrHOHdC5w6J8dR57KXFF45uYDSwNCw4Y/HxpnrNKPl4wQydsLdVziGup7pEnyTwmCbA3XsrfWYZPzonS2YOqc64LYNQhKogTnmzzMs2xZoolwANVmikq0SHu6jpnOpfI7gt8e2On2YaLzpnwGOXCqlU1kSiCL0i9xE0bvPmxnZN+ZqmUIuMRLLc49vKxwbkUPgv/akYNUNye5Wy6MbLpnFoaWTNNPpJUp0CZHSI1U68A3zP4Oc51LMweOwzMzR4X/TrP753m1O5pzrcbPLaxxVuvbXF26jjRDngM88FwMFh2esOL+4aX9ju+1F/jqn2FOXt4P6RioeE1DLFIyLxAWNabkXzWmNA2SldLznDs3I48FWg5Jmoi17VkKLjeMOoxiO8qzalt42ZkeSNj8XeQ4ZedIkI6waUugAEUCymdo2qOG6Py4SvVuAZJOBzE9syO3ANNhihOts4dJnUwhZzIzNzaWdrNJZZaxPrNMhLQUjY5kPRRmICFfliUsVxx/kXQh6TzaryANDcQE0We2fueUJAhxNDzeKKkfx2m1zZ64V9RjC6LVC9oa0LBkYndCNV9OQohy6gdyXx35pDGTNg3V7lht7jizvDKzkM8v7vJicmEzVYKm8Ji8BwMPdeGQ66ba1w3rzAfdgLiEYuzjpawKfYxZXVs8x0rpNKYvJk6utDNRq2JZNLF9aa/0+tSm1wNWZDI+y/nVq0BLNgZjV3E+X5j21DdlbAcqPAH45OiPdvh+FzrTUjUce1Uk2tprHt9Xd2dJoW8vE0VbYFUP64MmYm9KkymxmSX1W3d6EC0CAnJhPxxR2OaIq+/zs0O4xzCQpaadmrzAfHykj5zvqMfjrBqXkWtLBJuCtx9dNIl9V9MhVx3MBzTFAxUaFXeYSK6rTWzwIQeBtOhY9VDbLc8mDnWzOndAUd2mx1ziUtmi1m/wbSTgo+R0UzPgd1m4fc4GnZCBxjAmwnWd6GbC6TITs6DKE2p1s5oVK6+UOcP4vyGOdbPJ+tL+2o0Zc2xDLsWm77SFmpqpNZBWsdvnKf+LoTlyLtt7A6yxJhGLSCTq5SMlbmycV50DnKNmLPymHWVJgVrzWG4mJpbl0gypac9lCXODCBx38L5SAnISOMTXLwZ8HaWmbWwyaNGIh5m3ydm1zF9LWVF2oi2lMphJadViEyYhBcPc5JquiUVXph8GpJPmvUl5pZ71aXJ+uEwqLWmizX5NFJwHF3orcM5x8Cchdljr2JGG/0tvZunja8WCFJaXOZR1pZ24DZR+5iaLZqYq7/ggHmU6kGrW8R7Rts+VppNjTJVvoP2+QBLjr5wzxDxOcZvWoxPskjfSMfdXax4Y5PUqXe+1JBRMX1N2o46jtmlkoy0djruWuF6MY4+opbX3na5v09SUzlm1MZctneqxm6kqEdkfEmKUQvdxRLQMifFvFGm7OpNcIg4/4FFCt8Fm7HH0qrNUDM6HJcBJwwv85RyGOr5i0i7zh2mLkKFE7aKnNRzOeAKdVzOkRoJgwIhad9P0YYMl9KNUWp2kvJMaJgwQTcbCT0QRKML8yrdho4HS9W+IF1cIzD7bEnF188m2mO5Gbyxjrs776Wn3r1znHu8zlilrvtczlqDG2pm10AHXbI422ZRNYesjivv9ti9NbOlSIEullE1NQjOmFi9NTUqyAuwqT3bWAo0mJYmZM1IX0PbywA9RyFU10BjZ8l55XzoYy6tmqGn8Oqn51US3kxCd1id4KT9LSMmk8x17w7yOyvs6dy+W2t3yyZdLmYSfsR8drUOPC506TUD1vTJrxIYW0VLjE0agKOj84epqmznD7K/Rm1Mkiw11vKaagxSy17mOtw/mxPWtoHpY2qt9MIDQvto0yTN4o2muybha+YKeOWQnCLMM+bcC+cO1SLN39/s5RRklFSmTJWsw4Cj55Ml9XGgiYDCCvXxjC9tPNk8mojTbgw409NAShHW+O66tr5gvNfsKSZmHWnAOcQe7i0zHBZjY+NDZ9M1UwgotoxaCqfJPWxW47UEq73PY3M9+B6rcg/ytNticyzDZhkYI3+P3aNONvLGxU074wGE0bNwcalb0GC6sDHSFA5Fke76ObR3vqbaf1F09k2NTkIR0nzNgV49s2iVzumy428c3RWGFxQULoa1Ym532sUVTJb4uWa4VDNeMeeYl3zpvr6UNE1UUwf64uWm6yqVvdYAtKooL3VITrpgyzcmdJAJz1TGZmUclgj28Q5ngh0rmoN2joXa+jlkaE3LxGywZc4y9et0Zs4RQ1J7DaFBhvUTem+jhjFJTnpHdA0ZVy40Y4s5Tc425aEWxqht+tprL5h1CWuGZ26W/FFJPS+0mcrMUx1afYxaBHMrzGGYEx3qzMwuBSc7d5i0E2n4qaWtXjOD75c0zrKMdfabiA9IVHTZnAuQV9QIM3Q8z50TNGB+u7xRobk7yPAxNqydNSaoh8bapGIFMIvNOepKfRfKDrCQYKOhqV8ppaiBUuu0U7E2HQYbKtyL2udU2etUxIPcw06TXuA6UUc/n4wpqdaqhTWESrZTs8HET7FYOubM/R4Lt5fMh0ZDdsU+9w3GuGTLy2Ks7eZ0zsi85zlbrlik7wUkrUC0FWsEhTiuvS05p02Gr9aAFiGxuWWexKmZn21IqDsRFDbO0VL4Tm1WehNa0s7SZlx2+k1RDDUHKbEoFlMVwJY833Fa6RtBd1bC+yBNtApkjA1hYZPrzufiEuoFSyNJk8NCpeNn3CE3pp4L1WqalujSrUYWnHdlPN7jcKalsY6GYI7oBphyb6cWio4KSNpo5w4YjAB8YnVdsp+hDiW5OCdNdD45HL05ZN9f5aC/yrzbxvkuOdqkBp+Mx9o25hWoCAX9EqOFOHuJ8Ku1rjRGlpF2WQqXEQ/nexqzXO8/nSfvXVcmVn6L4n2Z2N7brjOx60t9CwSkNMSsw2FYZJs8QhKEOevKSyXWIWzbeg7yRhwBP6ZP5kBaOwhYSuZgQJB1IkwK0zVpqqNTc1vojqv0tTTJXt8xj+Xy7ifYbe1o006t4865mYofzm3UsRkNl1SxeI1UTTZ6pK2fJEmuN5CUL63ypst7D/QssDisIPrUJieLSVR5kVgpG81YHANHZp+OIw6H68y7bbphP92jqeqzp3mw4H3OOjTeFvNubbtUjUjOl41Xw37bCE0VjSJ09ZnT+wi6US275NnrbMVSda4878al9Oj8LKquQXy3ofRYrTnkUGnIe6d0zKpnKaolVxgKqYugJf7SBqi6GxvTxDWQUZ8FKlTWkyrecpxP5HbSnWf4JbVQqVCAFFPQ1WbTsVaqj+aElnqDiBdZukfRxndEhdJADS2xCudgKliotBTfYZmE9Et5xsqpVC8wPQ9FD/gR7/vErMcYv0Pw4kD67MDvsHB7idmdD/ngztsUdathr8GG72WwYRPx+djGztRGYwmFP3LEIRSGnKXxbXCKqV9n4kPpqI4Fh2afI3ZYmP04Z7dmk8qchk0i+gtMT+Nnhd9AtJUEqU2wY2W+0BQbqEj0Mt04MHsowNmEFs++ZzDzQqUf62JcrB+TGR1UM5G4nuTc+hzLJGoR4+G/2013NQ6vkVDJnk1hrpyYklX7nGQj6DLtuQ8/3dKmIvcQNQtTNq9I46ntStl4tEd4SRKV3UJHn1UzWzp3uZeeqI6tnTKxG2yYM2z606z7dQRt1jMkadmZBT3z6KhTTTjlmlWIMdvcpTlUhv1K1JdgvmVBt2bGmj3FhjnNlj/FhttgnQlrtmWttTgPnXfsDx3XzQ60sAt0w35poo1QMbfkuU1ZaV4cmNPCr5PwCj6/Z0suXml9aBbamFlRrUtLYyCaTAB9YToeKyDUeqnNjdpBDCxBbu8G3b0iliYXqACJl3dFpRDtSCuYS9lOEENGyIZQSVdK55hG6DnFILU61fscl9Xq5VgoJi2gSk0tPezleD02QUxrLaaNzH7OPcbD9gSnZmIPehaDp/OOhXMcuAXODizsAUdcz9dQJkwteer7heNjqSWbw1pi/wpZLK3dYNOe5WH3OA+bEzy8MWG9NUwtbLQws6Hk3ODhxmLKpcN11o7WeWUyY8e8zGLYL5ywaUxeb9wVQyjHbZ7TMmW5DhWaJrCiMLDUTAwPMsN41UNOSWunNt68GRwvnbMZmLXSMK6M8JPrZedel6HdWntEbdTH+DduB91dpJ1KJ9RSPBxQHq+Zy5iGGkMfrjOuNqYUV0KcW9R1I95TylhrLWnyNY6XTMHxtVxZ1trsdBPnjfWhkMZYKHFi1zlpH+Up9xRv29rg7SfhifWexng6Zzlyhv1+wo3OcOlojc2dKV+0E+aTgC/XMzCmeuraerX/o4TPhmUpErO1G5y0j/K4e4JnT2zxjpPw7NacjTY0emysZ2IdbawYtLeY8OrRjI9fP80f3ngXf9Suc9U8z7zfKZhT0oqNz9V9DDbgFrD4SksRn4HBxkSZ0ilobRslvsV5KXNVdh46zqcT3vNy6K04V94reZ70z7wmSoyGdK0pSmuptR+ciZLs9MbR3W01JRNf2cYwrvYlzz7B8SQLMl1TmQH6HMjAGtn9k/ZAZHzV2krj9g0R3VVh6+tx1ZJJpKU1LY0J9nCwLxXmW4dnCLbxuj3DBfc47zy1yfvPdfyZh6/y2DM72HUTGiseGo62Wy5eOcEXdk6y1a7RXn+CvclVenfIotstCiwCJRP73Pd+uaBnxuKH0GEfYvlmwqY9yxPuSb7m5CZ/9tzANz56mSe+8QCz2eKPBvzcgfOYqYUmiPn5xYHH/+AxHpqdwL32NF0z57rv6Ib90jFqHIP4HCyxjmEIHTqfUZJhLtscsRi6BGLR60bCpsbPyVGCSuqKt13b17UJRBMiCiN15QMselxFz5t+Tq313qHrHoQDKL4P0Yk3IcMn25XKZjdj/eKWJ8DFMFn98mpQg1wnxYBtAFpYa6ONl9VDUbcEHy9hwnRPuqVNSa5f03HSAyRU1C+Ffqxtmdh1HuYp3nXiBN/28IJv+1MvsPWNW5ivejb0XgNoGzaPFjz0/EWe+cQXeejjF9hot7j+2rPMJ3uhNZWAXGiWNr9lZ2bJ7LWz1GBZb87wlHuGbzy7xXseWvDNT7/MQ3+2xf7Zr8PPZpidXbh4FffqNgwOM2ng1Dqbzxq+4aGLnP/4Pp17lOHSs/xBO2fHvxLn1KEBVz4CkHTo0PjsszluXuvPxeYfYiSi6C2YG9AVGyEjCWphgynDdb1KohGB43Fp7op4vTyXEiZ1TkTK+/BSqswvD+Q20h1leI8PkhKXOm6kMI8kZ1QqvtDrLdxc1CDg42VDKUI8Md4/KEdbSHUVSR9LZXuWgC7O5lh8GOdI6eiRUJG+f1iEuZmBzAOExbXWnOEt5ixfe8bzZ595mRPf8zT+ne/APfRQfvb1jXCvr99h/dnP8g3t7+N+5wJf2DnBpfkF5u0O8045jCrVc3nMtlj4RSjKBHX+pD/Pk+ubvP/snPe97WXOfPsm5uvejvvqZ8Fa/LVr2EWHubyL2w9lqeyFNbhwnsmFczxz4QU+8HO7XJ+f4NL2Exw123TuMM1JrhaT2z2ZxiZPvTh4s8QM1WjlPR4XcpWNr143EjkxKuRofWl7C6NP7Aatybn2qTGIQuyluaxKq6UYvWL2pLGoiFFez2+8M+8OMnzsJIILO1kKl+SMppArrkIfkr45wuxLapRxEY+f7cIM56x7vNfdZCcjalmICCTtQTYlK1pDyfQpdXekjFHvLQ25TpowO5AiFbP2JGe4wBNbE959apdT7/b4p5/Ab27C/j7mxg3YP8A8dAZ/6hT0Pf6h08z+1Gne+fJlnr78NM8fXuCafYGF2R3ZMJs0rzoMpLWs9DdZ+jd2xkl/ikfWLW8/vc3JZz3m8bNBff/c5zGHR3B9F//SVfov7dHvekwL09k17NoEf+Fh7Lue4J1P/jHP7W3wuZ2HuGLOpHr90hhUh75kvrU3PptXDtw8MlGX1lBSj5X2mIFay2p5sZ6iE9UYSxM37dbMsExYN6dokOw76MwE/CKP1SnzYCR9O9wk+4HGSEdY3nw2vM9S1/vgeEu16GJBBo14q3H05YtVdqjyCQi80pkuMP9Q2mfitR/UAhPHXrb5c4OFdOuoPaQYdrxWqlFnpwmsUyDMXNiQ0j3USxXM+7o9w1l/lsfWDU8+tE3zli2YzTD7+5hLl+GlS/hru5gT65iHTsKpLby1mAsPcert2zzxuYGz7TrTYYsDcyVDQylNn8JRJhKGerPLxT4mZoMtv8bJCZw+cYg9txbaPd/YgS9dZrh+hNvu6K875ruWYd5iG4+dHDBZu4p56DSc2GTr7Y6venGfh9e22Dg6xb65jHc5xCh1CLQjUeZH297eO7qq5lxIQ7Vpgy3emUhSysIiddnptPHFGgcta0zMOlPWsd4ymJ6BWRIqMhapQy9rVjSD47z7+pn0fAsK9Y2mu4C0C2q9ZC4V3+ksOHKcXts+kGP4RTfOSq0zpglIMwc0eeLr4opFJR0D0mDAe5eTYdRmoL3aWgORfnDBOdcF+86LrZeftYjPRmafNptsmNOcMDO2Ws9srQc3xVy5GlorX7yGe20Hv72AyweYi9uYU+uYE2swmWDPTHl41nF2NmN2uBWSQpiTOuVGGzNsTGWpZu2xD+9HFdeQBppYWgvOGfxhD9v7+N1D+ud2WFyBxV7D0Df0fZibCY5ux2Nf3qM99WroBW8Np9fmnJ6e5MTBabbtVoK9ynjSvAJjDrZkBok3v44sVL4HCd+KY1Duof0BuhusRXIX8vU65gV6r4nX7GWDUanXx2ExJO1ZayLyPLq0252gO8zwojsqpFp03OiQhYSC8uIsVaK0axul7gk+W/WWE0dQcW48NmPAQ6mq0RbUURMp8tJ9FRmgSfaltvtrx56+ZlIBoye/NWuhKIOxNKLBHHaYGzuwe5CY3R3GJou7Pfagxx51mPMnMK1hs+3ZnKwxO9xYjlKYPGeNqv7T+wXeuZQvrx4qHW9paOP1hsHiD+b4/SP89iH9tmex19J1FjeEgdvGY6zHLQzD9oB9bQdzcg0MbMwWbE5gnVlUmctcg3x7FSqMv6fQaUzEEede+azZ4ZgFgbSdGiLTlRlsPUdjCxUXy4o5M2BVIs6SdC7agy3H9tMzxQ2hgJUrP9OdorvgpXfgDV4qwfo+KciF9FN14wZTlqzKYZpSfZI0zIF5sj/TRMcdXq4rrZ9tLHXUc5RMjLrmnU7Fbew0jUHIpI1DLT7GVeWl/nExXDOYjiM3sD9M2N+fcvqgI8K+Ak1jAofzmNZgT0xhfQqLAbfTcdC3CVYqmWO1emzVd0AqK6WLPBQaTHpjYRxN4zBTi7EGpg3NCcN627NO0EZdB25hEvN7H/9nDPb0lK0Tc6Y2VMurKbXlJkv4PG9DsflChgDrrEqpWAOk/IywTiyGgMyUWnue3BhUh8ockgcR1pDxTUxzztEcHUnSpAWU+CQ0hn55PSyjLd9oumtxeFHpnetGq/qY6DW3apdPk2JKj6h8P0jTPmlL5VWjPlWNBqK6qsobS9MLKFXIWyWr/Ad1G+jyucYx0wM9R77jxmKNS7tbPHL5Ek0/wPoMe/4E/nAe2ii3FrM+Db3YrMW/dIX5qwPXFi37XVg8TWW2FPdPDNCkedPfSfxZyDHQeUcvUz9tYHMtdOg5c0RzwsM0Xmfh6K/3DPsOY8HODGbSYFqLt2HyBw99dGDqpBxdX6DWUAJjukLTCpt28Jn0hBLYhReectMN92hSLTuHaHAZgxFOUu8rRms8Q2wg6VKOv8xXDiuPw7J1ya6xsmDaQ398yfDbR3eB4T2hG0qOP0o1mAxNjCWQYp75cTSm8iVTQQF0spqtrh8qqCWpLOEyTdIp1Niszum0Rj2OpNIzidKhwUeJU/sZsvNMUHhBeuxzxMv7m3zi+gke+aM9Hn/3DcyzT+CffgLmc8zBEf7EFu7USTixhblyBfexF/ji5x/ic3stFw+P6Mw84eHrppgBLBJ/6sKZhBx50YhkjM73LPweu/6QnW6do/kE2h7/+COY+Rw77+Cog/VpaA3lHObSLvbqHGMNZrOFjSm+dwyvHXH56kmuzmHb3EjloOMg0nxLZmB6J9H/Mbg5wxALTNpQWHPabAW72YX01Fo7ca6jR0JkjtbMGOiSui7lrXQfwZpxB0nxNVlaS70Bidhkcy3UZ5Dn0SXBxuoNSOELHbV6I2PwcJdLXHmvIZ6qR3tsClF767V9XScvHCeNJbymbbukqptcuCB1FK1UW3H4COkxabv89bLBivCezXn/UqgipLoecGm+xWe217nw2jnOfewF1k9twCMP40+eDEPd3MSvrWEODzCXr7H9Kc/Hr57h+d2BS36bzhxFp+akAHnIOCU7T5KQxAGKy0UxtaOyd3P27C7bizNc3tvg8auXaQG/Lk7DBtoG1mcwbbHWYmb7MG1gNsG0DX77gP0vGp7fPsmlw4E9c6PIVRgroCE0+K4ooBHeSdgYWtaUVmKXGD6o+bHVVLykMJ50p5Hvjy3M4WNZc2kuocpv1REbiIKnWhe1OarzBqSyMHcAVgt3jeHzBDrf0ZA7uQSmDla0i9JHZ2stS8smd2IRdc8KiCdXdBGSiR9iUUcrHW+ifagBFXUCz5iqLzagAEl0wQlrAoCkXIxNkUcuiL+Bjl1zDecd3fZ5YIv5b72Vb7v6EmeOOvj6r8E/8QQ4h33pJfiN32P7I9f41595ml+7ZPjDo8tcta/Q+3kIK5kpg83PnjWZPkp6G+e1iRGG3NJqcDk11PmOPXODVw8e49PbJ3jsE3s8/vbnMBceCuG53uEXfWj5PJvCow9hHn0INtfx1sKXLjK8us8fvnCej15f48XFNfa5ugSKcTiM74N0N6EoZHonFVbemhIUJeum3oxTiDWaWBM7MERNMlSpzd2EQj39Zuk9j0loXfdPnKJJs1Sw5johSLcpz2jABXcCYSd01yQ8kCU8KtyV4rEh2cWgNIBKImgnSdp5scfs1CNllaklYI6daxVeM6yEWGTTCdlagviSXPN8bhpTJFH/Cynkc2imZ05n59jtJ/BsAo/zrf/6Jc5sfxz7NZdhvsA/9xoX/03Hr3/pKX7jsuEzh1d5zTzP3O2mhpjYjdCgQdmc2kGUJPyIZJSFD2GDmPs9rvoDPrd7mkdeO8ep/9+LbH71LvZsqNmO87A2CZJ+NsVPJ7C+Dkdz/I099j4Pn9nZ4nPbA5ftRRb9XjEX4RJZUroYbk01EtQm39hp8FEotb+w3cUhmp5zqPwzue+ALjo5JhRkTWhEpDEZel1gGWJM3ju3BL5Kz6AiSrnc1p2T7nDXGN7HeHzwrDZuBk2u8uqI4RfGO3Omq9ykx9eYl7ccQQ4JSnqo3kD64TA4tiIKrLGzJJnFyz0wD8pKVM16DvE2J95ICKwAjqCljlLz1CbT2yOes3P2th/n+vw0n997hj/9/CFPnXyFed/yyv4Zfvf6Gn9wvecz/QtcJmShAcyak8zsBmvmJEMsySyVZ7RKqkNv5ZwumyedP+CKvchnb6zTuRk3uqf56ud2eeYt11g/39M+1NBMG8xDkUHnC9g/gKs7HPy76/zec2/hUzcsfzy/zA4XkYKaGno6uJj5ZkoVOb1PWxarSEVDpVactpFF2kdyvgsgH5+r0tamm6wTvcno45bnqJTeog0ZGrytuyZJGHaSTKa8Bu6cdIe7KuGjWp8cFw5jJ2gYYgiTLUt0ULaVOlZojMHC52XIR0v5OoUySUM3pPCfa1oak3HS3sZQEOKRtxErrUEjObwoIbmMuIve4bqKq3M443jFzjmcX+Daa+d4+WCNx66vs3BwfQF/vHPEC+ZlrvjnOOpvhFpxsaRVw4Q1tvA4DkwwOzoO03NpaaNLWJVzmDUR7x2H7PCyv0q3fYb9bo1Xjk5z8XCdJy7u8/DpXU5duc70tf2QEN87/MLRXRr49B88widubPD87pwr9hUWfWgAqbHxqc5hXBOgcAACkjEkEA1I7nrtNW+KuRdHm36nqY5ivSaoUIeqarD8XppyFc5CgZzEH1KEROuwbfJh3TnpDneV4YOUd17aJ2XJImqyODjGJeRy/24NI61Jq49F9pNfQEyGKMA52NyqSdWRqwslSqy2rtnmCVVWrAciesvF0KBz3VKHFaEgAQZ6f8DC7XFkt7luT/HS7iOc2DmFx3FojrhhXmNvuBSaKqomBtaEtlczH9RtZ4boWNopHE6SUx4G3BUw45xAkpsvdv6AK3yJbXOZVw9O8/z+Of54Z4PHN2c8ee0UD7/ccfITHYM3HA0NB0PDjUXD5/cantvp+aJ/hV13kX7IG49seNIOq3ifBfhnGV4rITdJOZX3qguNWNtHE0vVx1va1GyxHsawIN67oiyYrNVaAOmwcX4MhRxMeBHVjOIOxd+F7qoNL4i7ITpQZP3VHVdqqiuhFmpZdLIcG4euvP26IIRcS5jYmKyuaZCMJWTcSex3eXyhYaIxsca+zzaxTo8dy7MOTqBw/cHN6YdDjswNDpqrTO0WELSRxbCXGh+mcROq2g6mCw09RmL+EsMWTzdA57ITVR+X042DY7T3c47YZt9cZc/eYG//cS4fneDVgwknp2usNWt0Do4GOOw9+73j6vyIS1znOq8ks6KeBz2PSw4z34NpwbsiXVnHuIX5pNAIkPwX3jher4x5vfEWERUsGIu3rsiWTJWT7KRM+jE5d7+4h4oWpR6Bd1idh7vO8IDvcW5B1+/jm2Voq7Gll1vH2UVKLJVIjmpy7SjREmMgJ7hoaGlC58VQlaaM5GPJ6VPjs0v7PI9Rh7yWpsKH9lQpfqvG1vUHHNprS/NQLzDnO+Z+L2kiEnfWMGPpiS5pn86ETUKHPPVPiUT0sdwzQG8POGp2uOLP8KX9h1jfW6OhoaNnbubMzSELDunMEZ0/ZD5sFziHZNLEZ5S5Kr3YPW7ILbpqXIH2zaQ21SakD4cYe0RY6vNG7O9axa+jKkCuwKOiHtpxl7WFplDn9b0ksacfDhn8AskevZN0lxneR7u2p3cB06xLKycG0cUtVfhOFtCYpJfPNbPn68bc65ECEVKFNjFSlA7ZK18uFv1dEe+WxS1xftMvORPr5KGQ/tkV95Tr9f4Q48qFpZk9e7tDrftDa1NHGzGL5M4S+z/O/MlvRzSbIWkc+bMgoRd2j317hdaEmLiUpx6GOVKGLBVxjMlMNdMdf994jMnONp38UwNmUsqvdzi9JpTnfuyeS+snmhNiMsj7rEOcVsw2SdCpNo86AuBcHwuULHBuETf9Oyfd4a4zfCDvO5w7ShAGLS21ip7yje0kVrwJsVpitZKa4cf+ttikqnn53bZ4pxszZKbXsdNauurrNnaqpLdbWtTWTGIkonIaulId1M4sa9rCay0xXkub1G2rHFmgWiTF8FCTpPosdKmNqnGYiwbdzrmep2Kxe1doLhKCGlxo8Zw2nMpEmrSbodS2ndFQSjrtM0lrQc+fJDRhU/krObLQAk3WRAI6L6MmC0CXui7VdfTGIM+H70Mkpjq2HrMUANXfL2eCumSCDe4InyT8naV7gOFDYQzvFjgsAzaBVfIRORbf0EKEvAa7qraDlde1mlBJHhHJLUyTOnlWDRPGbbFSqhf3jIum3kCClhLVzyQ9qh5pWpuIzidncrjHub7AF2hml/bIOjSUnoGY/GEaUAkkUzYisw/o5BlNukpLrdGkuTAUm5aeQ9mIDCEjEGAw80JSary73ihlLpe+Uz4AnZoKoU21YU4RWpOoiNpIUsSkfq9RcGg1vX5mfT/3OgxbOwODhI9Am7sg3eGeYHgQjz1+gfdtSH5pyjZHmiwW7Cyp+8XOq2PLx8ynxN19xNI7LDRhMRZNKpV3uI4K5JJcZcqkNZNUh023INbOHlE/i0WskntEm2lMblRYXsNW/xpSOWdhTKvHGbK95JkmrDFhRh8s7sKvYM0kheqc61IGms6TT/NYqa9aeorJkTYl3fCRDPzRbaM1WjJd3wXJHTYlmz3llZblvaPzEgEYEpLRMyQUZDF/lCi5dE3UZjPG7CyvSZ0Wq/0gtUY4+AXOL7iTQJua7hGGB0mbDRNiafwsqsJteimpGqnJQAaNghLnjlSd0ZQYA7Ltr+K3ltjWOUodsYuTFNLMbifqc13mSsp1ZQ93iAXnDC9YRutlG7PSVgiVXV1MydXSs5Cqpi+YJnwfmnQ636cW3E1ls/fMo809LyIbTknRWgWuzZxCAo8wuzyXmA5hHJP8HJH5JaqhnWaycUlbrJrJ5Xp5rEpTUe3Aw4aYtYHC1+IomF4+TxuzGpOmsZCqjEc7ioGgvfouS3fujnSHe4rhswNPcuWNydjocER2qoWGk8tJE0vlpVHhFcokGbFVQdT9Ka5plxaVJrm2bjUtJA6q5B6zElePnmTVsFAYR2za9IxiK0sXVCngYKrFz8DgllVx7XiSyq16/iRjDEJRRslac77H0qaiIrUPol78WosqM8Ea5YtQEk6wDFIo1ERJHCMmmqnSRmty/Luw6ZXTUtPYXGiJKxtV6gpMmFNh+gI+K++oamRZdyKq71tHhUQgDYUqf/foHmL4QAF0s2AwLda3NEzRbX2FJHQFmQnr8Fzh0SfUER8Usg8y3FGYtVXFMhxjEqUpFoC0D9LjFwrq7LQYs1zD+hgCxOKZFuqjZKk1JjB8E1MxdbMGYmWfeuHpvx19qqc3mI7GTOjNUUpBzRlji9zkwy/Xf6sZqyatqZSbQjYBBMnoYnxcYuiayeWd5c6rfaEByb1CZ9vcpiwjMVUUhVzEVKC52fwKwqC01Zd9QfJPO0UTqlKEhoy/YnQZSzCL5gzuSKnyd0e6wz3H8LGyre9x7oghtmkKXulyAYXFPRQvsbiS2GNKFRNm12CVkL/cFi8qaQdJwt18VxZP+FKHUJPRWgCkevHhPi0kjSKrxlHVFVwAErVo8L6lT80ie5yrHGjEBVlVoRV1eWCunJYleEhs5FrjENJhwOJzxehLPoYUDclhyt7PE54BKEwnGVdmfldsIvK7biIpc6ZzIcTZKOafgGckTBjmZTnMWoy92rT0sfp4NCKwqlrb+bCZ3umMuJvRPcbwkFR7t2AwR3S9pW3Waew0LYx0pKid4l1NyLiyZJN4v3V9uvRddFIdF5POdmy2jb0JSR7yoptY84xYdTehvwpvebl4INbFj+qtkO5OEmq2hyUafBmhis9g2mzyqHnI0mdYeh4J6WkEoIxF5sgY1fooepXFWRdKSuWWynUOe72OtSakY9r9cFiUopJQZWac8R54pcRt0zhESi9lQyrzT3f2rT3rtV9AM274PhS0CNGhcn6K+VMFM8N8h2fthn1cId3vLt2DDA/C9Do2b4yNlXGUk0ar3BGjnr4TlZ1c+aUO++idOd35JvY7ZMaRllYJymmUI0qFsbwZgl0cx53j4hPlHBspfTTiwNPVU5dQZ6oqq2w0wsD5msqeRZiyWa6Oo5hdO7ASWrBS+Ytxyrxii6KevZurOHSObRsXguuNUunHYccjzrAlH0PcjLRQUFrCceZZeY9lMJbWHoxf1iRRiT3GBO2r94u4uR3ddUedpnuU4QF8tOePGLDQE0ZrKXZ30OpsRjQlstx0YxWvulxHAzYkXVfsv6Cy53x5S67vHhh5ElI2zUBnMmx0YJFr61nJnGpoVZsrqZem69bXUtRFh560O4IY65ZwnrfFs9cqaPZbtAj6LvxdYs2lzqA4s25GtTYVVO5JzFnX8xoRZsd24g2mVXHtER+KDmUWYTNPsu2XzlHHau1ENqbU2iplPWYhokn8C96WIKT6Pothj0W/Sz8c4NwRdyvmPkb3MMNDSK4hODwA69qk4kLlpCliqbq5QYmllvOyY8slBJjHJRu8ITvXxP5zgFGSsrgeAblmjMXbIfoLMk5c4KHWZ+ayOq5tOlxqTjiCs4/SVreqyuq4TQ43yBteVneXGSd8n1sr6YiFa0ITx57DwPRqfrSTSoNgZB4kMSdHS5plxlFa1ahfAFuYJTWIqgDRKI+6ttfl3dYmnp4X2Zxw5A2uoiWnHE3SGGoSc6zrD+iGgyjdB+4FVV7oHmf47MQb3BEMYZLbZj3lF4u3FirJXlEBj5WiiTbgn0XFlFRNCA4+YXrIaa7eDBm0wfKLtxHZ1tpZYEzrRsdVhOekQ4wpAR/eh+Sh0Ogyjy9lE2KxfoKxE1rTMHiFGoyIO237jknLOJgU5dA2qNQXDL8TagBaKkZczvjTtQYlLdga8QGIVJ0UjCwOvsTAMXynS4aFMQ8Fw+q51JTSbkdKQXvvAhQ6hiAbO0vmog7Nia9gqXlJ8XzhmJzdeMRi2GUYDvC+42575Wu6xxkexJ43Hpxb0BN32SarZELppYyUDtbe/Nr5p3PB9bWWHFPHjdA7nMnJFU6YUcXPoVwo4dicNis2rzN9QPwpxtULXmPC5TpaTR2LRRfSuHJQjZEw6XHfl89Q3itUeHWgIL55cwztlzXGXyIb4tnW49PSGkhalh15juShV05WXeOg0ACVttIkhrVJ0kN2IOqx6ucuNYaABOyHo+SkuxeZHe4LhofsuT9KVW6c70JCRpQakBmqLmNUk1ZhtVQPTCodacp01ywxlhnG+Y7FUIayJH5uGpsgtoXNagYGuSdNUu+tCd1UpR/94BdLTiT9u3MdnSmdheE6k9SsIRwrhTDyuHXIUMJng58z+FwrTnvwNXBIM1T9bJJCXLxBuVasuZer4w45VKaqAOV7K/UfaRTSK0ds2fdd/DBjNrZ+b9JrXgO1rLVFYRIN2JL57V2ugS9OUWnR3Q9HDO4A7yQx5t5idrhvGB4C0/sQrkMq4QTPfRPBLQlKqdTl47zJkKVBvkPw6DLCYLWDqAyllQUtAnY8Mn2MOxXNKQjpm1IbfiBX100OLD3OKqNOfpbOxn4kjDZLDD1E77VINgEcSU58dkqWYchkl9vscddJNUtz6vvCLNFjtghMOhSD0O+h3nzjoHAmMHrS5kwGO6WQrHp3kkdf338sjKa99LKxjZUckzn2fkioRO1DKLLgUsmze4/Z4b5ieAiZXQbcEQIPyRKtLRZpHaKpGXvZns0x7OCBD6Wcayo8vFEFrdVHGUNLifwbYmxd4K7Ol0zrcQn5VjukClV3DO7pcuZdY5vkPGtN8CUMalzCxBO7nrPYfEfvD5S9OyTGaM00aQOp2KTPDR1rDUTGoW1tEze+5FhUTC5/14Ac2bxqaR+up5p9qndXvqtlcy9hJ6IPSG9uxXqInw0+VE7Oz5aLkuR3H5Ji7nVmh/uO4SFJ+lgpZ4gvUGLdmhE0c5Zx5rFqM2W2losJJxo4ItfSjCjn6usgdcpj4Q6Nv873GtJxkD28Ra30sQWr1HENExXm9CYzaqM85QV01YZmjrqpYr5mjlujrltfRz9vHS2pQ2fGW6wFr5yKlkkyz8JbdYWqrN9V0lpMNt9uxlO1pld8F9+pI7QD16ZLuQ6kyu0k3SvXBYhZhDikEOX9wOxwXzI8wID3AAsGF15QY6ZFkwJtg2s7bPB9Uis1aekCUCSBKCbH6JZIMW12zHGIxo9nFTjZfQRsP5TSTaRHYKIyClGTllral5Ht0zZ5u7VkCwy/ls0IhKls6qUm49JUS9D0+UhItDwvIvfokmOsLky6rCGVWAuZewEgBdW7bGZRJN2ksY33+MvRDu31X66mo1tHyUYWzg1Va0im3L1ps9d0nzI8SIwegvce43C+pZGa4DZIEfHCyuJuyEwvNAbiIcJ1k9qqIgKeUo0FZdvLNeMihrL+najTUtU2nY9GktnCw61DREVl3+L+TYFWS4sXDeQRjaWEu0okQI9Bq6w9pJJYGbeeQ31jkQEZx3GOs+J4mc/oG9Hj0F1/g1khTkSXGnfq62iqmT1rHCrqkTSOsvyYnKe9/INb0LuAnnMu5LV7CR3fB8wO9zXDx0o53oF3DMZiTEvqZtM42vj+2gK2OWSvLApQUzHs0t2Uqh6KZ2Qbr7YBtUqri216Hyqo6tr1CexiVM03M8ZEoX2UM47aaabj2sUGp9TTdB0s1pTMLinDdQnu1OqJDmf7nFVYhUPD26hhvXHDMcsprWOqtoQwCyScHGtUZlzEOQD05mgUCxDewbDEsDfz2tcbeL1JDENg9n44wHupR3f/MLrQfczwIJMdMsRiJxvjkqfW+Q5rJgx2XniZIav5UoxC+r5rWpLilLZlrfIWtmzF7Pp64mjyKs5so8Svyx7rewvyz48wjFy/VpVlnDqE2LvAYJrZZZy1Ay4xjSOAgI6B2moMetZSpEd66SGX+RNp3ijtpWeRxq59CTnLLbcUyxuI7jy8HH/PzzEsjbOeY60NhbBnHxNhjoJj7j5mdrjvGV7Ikzz43uOGHmdCmCQgvNZCzL6ZJvvembyIUqy5KbPbxpM4AkMJWEc7vPSiKxY5y+q/aAw6ZVM2H2taBLMvY8kAmqyi6+/033pccl0I98rNJLv02VJoUj4n141bStapIhLpWSWMpxJJxlp9QWB6KdopZo5R9Qp03D/Y7CHC0ZLDcUFLI2lddc8CyM629P5GJHo9Z2F8uaT0MBxEFf7ed8zdjN4kDC8U1XyIEr9nMAEXPvgFjZsWjO8Vgkqr+dp5M1ACMXTByJrxtNddS7mlEJoCbeBJFV/0WCRmLX+L02uoqqbUUl0+kzVZLGzf0dNHWzlvMoMwt3Ik1irwcRtLjoyUUFQ9j3KVOpqhk3ZCYcms+hdvVSIZUdOYW0KasJ/HsFnGF4xpKZryOMewFlJYtC9i66E01b2JnPty6U3G8JBfiEj8Ae87vF/gTMvgpjQuS3ydFKEXbLiSeIPDFRMyCxf7m2egjTCZZjYNbhFpp1sjhXuMQV6bhC7Tm0oG+GjJr8KGsSaAjEMzjhTGqOHFPoankgda6gdWWW0SBxfsgUQwZBzyrOIHkE1KtJa6Ll46D1swqhS7rM2WELloYtppVPtV66yMnszzVuMt9OakN+Bw/T6p8FKhJsfWgwb5ZqA3IcNrEubXzj2XJL51ObNLCjwEQa1t4CpLLDJ0wagm5rkb8L5NlXUgp9+mWm4ysrEiEmhgyHihzFT2qX7SaBOH0lhtcT+tpqcxKTNhiGm5ySmpTZ1YXUaPtbUzvM+ajoxXwoABi+9SpCM9m8n49GTG6Kw6YGInWN/mIiDKfND1DmTuA5OTNyQ7YfDZYQngB4XIi952HWEInwmCrk8/3wwqfE1vcoYX0s69UuIHMEdL46fByda4ot1vQn5VOeOQQzo1KCeAOko1W45LDjttU1bqsjCPeKOli4o42JInvJaCsfKsxxXFQmBZtQ3nj1eXIVZ6NUmSVw7PVEDU5TBa2qQkJyAgC9Ml43NKE8vSgZk3TekX0HOEg+Sk1LXkdIuvMGEUG1JjJgmFVz/fUMXuk2SXuPp9GGr7cugBYXihMYnfY0ybHHy9OyxsvCTB/Axvw6KTVMjkhY9MLNIqVYGVxBQfpH8Og+V89cZUnVwrssTUV6kzF02JVIWGvMEUaD7xSOvIgpb0Xmf3ZVNBmwmyydU2uaj8hYRPDRtz+y9PbhgpG9jErsc57Irr4Illw+ZJq9EOTRl/Br90xbM61zPYeaGdiYc9dyeOUn8pXNe/6Rld6AFjeKFa4g8YTOjD4kI8H6XyYmLraBs1AIL3VhaSleYZKkvP+6ApCOTXSsXZKvFC+wFkTGUoq8HE5hapB3la6DmsV0cIChMkfR7PN0OoaRfXtTyLqPCCo4dsZ+t7DGQYsNBQoRszDiFWELK2AjsFP4hu1jGYeYp2yIbp1CaTVO847ymFeOiBo7KwJXqeBALrSJt9YvwsBN7MjC70gDK8kJL46aO8AXiIDGXxNtj9fexiIr29fZRDGIsZSiltXYu1bYbU+lwrDogtqW2o92ayZNRltgYzp9NFKURD8DCwoJf4eqXmZkhwZiDdJbeYhWIDChuV7mk3VCq4lKoq+9K32Vch11EblMCeBe+gJbAeRxoPktWXN1Gn7Ov8rMqpWWhJy+p8yeTaxHnzM7rQA87wQvULjx5+AD8AFucCCNyJB99XKavFJSJzmDZUoXH9kuSDjPzTWWf1ptCrFFede730BEqayXWT1qC8z6K+j5HxNlW40SAiXeQzJY5E5tNzINpQjhqo0KUTW385JJaO8dqz7sKxTiP/ejXv+ZzXo1KKh08eVDLe+1t6emMe5L1BdFGL4Rio2dIpQTMwxc9Ay6Giam4r5tX5+QUTpdGVdq5mvCxt3YhkHMkHOO566m+tHgc6BliU7hOfaWkuLMuS2GUmrqMRhQou9PoMry7wZRx7/5EWQsfRiuG/LLoZs4+kYibkjFXfu1Iq1d99RTSGHFtWb0smHTsuDqlg2JrcLUlVGJOs+X5jG2fJzHpDOu5+xy1dc5Pv3ry0Yvg7TmMbwpj3/fWZ7su/861pHrWEfL3zliWqplrCH8ek9TX0PW9ls3vwmPcroVth+BUX31bSC1MW9a1Ip+NQXMdJqjHp+JXRcYk4mW5Fmrtjfg93GLtrIHPM9VcM/kbRiuHfMLodi/a4a3wl15ZNot6UvlIz4nY+34Opgt8NWjH8A0M3k7S3SmObxu2gFbPfKVox/Iq+DFox5v1Of3Jv0YpWtKL7hlYS/oGjW1HLtVNwJdXfTLRi+AeOXo+BbxFYtKL7klYq/YpW9ADRSsKvqKKVCv9mppWEX9GKHiBaMfyKVvQA0YrhV7SiB4hWDL+iFT1AtGL4Fa3oAaIVw69oRQ8QrRh+RSt6gGjF8Cta0QNEK4Zf0YoeIFox/IpW9ADRiuFXtKIHiFYMv6IVPUB0y1VrV7SiFd3/tJLwK1rRA0Qrhl/Rih4gWjH8ilb0ANGK4Ve0ogeIVgy/ohU9QLRi+BWt6AGiFcOvaEUPEK0YfkUreoBoxfArWtEDRP8XftM+w7eAU3wAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1,1,figsize=(3,3)) \n", "plt.sca(ax)\n", "plt.imshow(recon_OSEM[:,:,20].T.cpu(), cmap='magma', interpolation='gaussian', vmax=0.5)\n", "plt.axis('off')\n", "plt.text(0.03, 0.97, 'OSEM', color='white', fontsize=15, fontweight='bold', transform=plt.gca().transAxes, ha='left', va='top')\n", "plt.xlim(35,157)\n", "plt.ylim(45,152)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "54d3ac97", "metadata": {}, "source": [ "Lets load masks of each of the NEMA spheres:" ] }, { "cell_type": "code", "execution_count": 5, "id": "14974d6c", "metadata": {}, "outputs": [], "source": [ "mask = np.transpose(itk.GetArrayFromImage(itk.imread('/disk1/uncertainty_spect/fdg_scan/fdg_spheres.seg.nrrd')), (2,1,0))\n", "masks = [torch.tensor(mask==i).to(torch.float32).to(pytomography.device) for i in range(1,9)]" ] }, { "cell_type": "markdown", "id": "b564f1d9", "metadata": {}, "source": [ "We can estimate the uncertainty on the largest sphere, corresponding to the first mask:" ] }, { "cell_type": "code", "execution_count": 6, "id": "d6a98233", "metadata": {}, "outputs": [], "source": [ "uncertainty_abs, uncertainty_pct = recon_algorithm.compute_uncertainty(\n", " mask = masks[0],\n", " data_storage_callback = data_storage_callback,\n", " return_pct = True,\n", ")" ] }, { "cell_type": "code", "execution_count": 7, "id": "5062ddae", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated uncertainty in large sphere: 1.63%\n" ] } ], "source": [ "print(f'Estimated uncertainty in large sphere: {uncertainty_pct:.2f}%')" ] } ], "metadata": { "kernelspec": { "display_name": "pytomo_install_test", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.10" }, "papermill": { "default_parameters": {}, "duration": 682.085087, "end_time": "2024-11-19T22:42:33.818138", "environment_variables": {}, "exception": null, "input_path": "t_GE_HDF5.ipynb", "output_path": "t_GE_HDF5.ipynb", "parameters": {}, "start_time": "2024-11-19T22:31:11.733051", "version": "2.6.0" } }, "nbformat": 4, "nbformat_minor": 5 }