{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a80dfd37",
   "metadata": {
    "papermill": {
     "duration": 0.003296,
     "end_time": "2024-11-19T22:25:19.153105",
     "exception": false,
     "start_time": "2024-11-19T22:25:19.149809",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Multi Photopeak Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "afbded82",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:19.160631Z",
     "iopub.status.busy": "2024-11-19T22:25:19.160019Z",
     "iopub.status.idle": "2024-11-19T22:25:23.199202Z",
     "shell.execute_reply": "2024-11-19T22:25:23.198483Z"
    },
    "papermill": {
     "duration": 4.044966,
     "end_time": "2024-11-19T22:25:23.200988",
     "exception": false,
     "start_time": "2024-11-19T22:25:19.156022",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from pytomography.io.SPECT import dicom\n",
    "from pytomography.transforms.SPECT import SPECTAttenuationTransform, SPECTPSFTransform\n",
    "from pytomography.algorithms import OSEM\n",
    "from pytomography.projectors.SPECT import SPECTSystemMatrix\n",
    "from pytomography.likelihoods import PoissonLogLikelihood\n",
    "from pytomography.utils import print_collimator_parameters\n",
    "from pytomography.projectors import ExtendedSystemMatrix\n",
    "import pytomography\n",
    "import matplotlib.pyplot as plt\n",
    "import pydicom\n",
    "import torch\n",
    "import shutil"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92d0f7a7",
   "metadata": {
    "papermill": {
     "duration": 0.002842,
     "end_time": "2024-11-19T22:25:23.207339",
     "exception": false,
     "start_time": "2024-11-19T22:25:23.204497",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Change this to the path where you downloaded the PyTomography tutorial data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f9a1f317",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:23.214508Z",
     "iopub.status.busy": "2024-11-19T22:25:23.214013Z",
     "iopub.status.idle": "2024-11-19T22:25:23.217499Z",
     "shell.execute_reply": "2024-11-19T22:25:23.216912Z"
    },
    "papermill": {
     "duration": 0.008696,
     "end_time": "2024-11-19T22:25:23.218974",
     "exception": false,
     "start_time": "2024-11-19T22:25:23.210278",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "PATH = '/mnt/mydisk2/pytomo_tutorial_data/SPECT'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "737c685a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:23.226402Z",
     "iopub.status.busy": "2024-11-19T22:25:23.225801Z",
     "iopub.status.idle": "2024-11-19T22:25:23.229696Z",
     "shell.execute_reply": "2024-11-19T22:25:23.229195Z"
    },
    "papermill": {
     "duration": 0.008913,
     "end_time": "2024-11-19T22:25:23.230969",
     "exception": false,
     "start_time": "2024-11-19T22:25:23.222056",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "path = os.path.join(PATH, 'Lu177-NEMA-SymT2')\n",
    "path_CT = os.path.join(path, 'CT')\n",
    "files_CT = [os.path.join(path_CT, file) for file in os.listdir(path_CT)]\n",
    "file_NM = os.path.join(path, 'projection_data.dcm')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc3bc703",
   "metadata": {
    "papermill": {
     "duration": 0.002651,
     "end_time": "2024-11-19T22:25:23.236308",
     "exception": false,
     "start_time": "2024-11-19T22:25:23.233657",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The idea behind dual peak reconstruction is that we will construct two seperate system matrices (one for each peak) and then combine them together. Most of this is similar to the standard DICOM tutorial, but now we're constructing two different system matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "19525b31",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:23.243648Z",
     "iopub.status.busy": "2024-11-19T22:25:23.243170Z",
     "iopub.status.idle": "2024-11-19T22:25:23.259978Z",
     "shell.execute_reply": "2024-11-19T22:25:23.259312Z"
    },
    "papermill": {
     "duration": 0.022416,
     "end_time": "2024-11-19T22:25:23.261474",
     "exception": false,
     "start_time": "2024-11-19T22:25:23.239058",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index 0:   Name: \"Window Group A\", Energies: [187.19999694824keV, 228.80000305176keV]\n",
      "Index 1:   Name: \"Window Group B\", Energies: [166.39999389648keV, 187.19999694824keV]\n",
      "Index 2:   Name: \"Window Group C\", Energies: [228.80000305176keV, 249.60000610352keV]\n",
      "Index 3:   Name: \"Window Group D\", Energies: [101.69999694824keV, 124.30000305176keV]\n",
      "Index 4:   Name: \"Window Group E\", Energies: [85.879997253418keV, 101.69999694824keV]\n",
      "Index 5:   Name: \"Window Group F\", Energies: [124.30000305176keV, 146.89999389648keV]\n"
     ]
    }
   ],
   "source": [
    "dicom.print_energy_window_info(file_NM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c86e2dcc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:23.269339Z",
     "iopub.status.busy": "2024-11-19T22:25:23.268939Z",
     "iopub.status.idle": "2024-11-19T22:25:28.327155Z",
     "shell.execute_reply": "2024-11-19T22:25:28.326316Z"
    },
    "papermill": {
     "duration": 5.063979,
     "end_time": "2024-11-19T22:25:28.328882",
     "exception": false,
     "start_time": "2024-11-19T22:25:23.264903",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Given photopeak energy 208.0 keV and CT energy 130 keV from the CT DICOM header, the HU->mu conversion from the following configuration is used: 208.0 keV SPECT energy, 130 keV CT energy, and scanner model symbiat2\n",
      "Given photopeak energy 113.0 keV and CT energy 130 keV from the CT DICOM header, the HU->mu conversion from the following configuration is used: 113.0 keV SPECT energy, 130 keV CT energy, and scanner model symbiat2\n"
     ]
    }
   ],
   "source": [
    "object_meta, proj_meta = dicom.get_metadata(file_NM)\n",
    "photopeak208 = dicom.get_projections(file_NM, index_peak=0)\n",
    "photopeak113 = dicom.get_projections(file_NM, index_peak=3)\n",
    "scatter208 = dicom.get_energy_window_scatter_estimate(file_NM, index_peak=0, index_lower=1, index_upper=2)\n",
    "scatter113 = dicom.get_energy_window_scatter_estimate(file_NM, index_peak=3, index_lower=4, index_upper=5)\n",
    "amap208 = dicom.get_attenuation_map_from_CT_slices(files_CT, file_NM, index_peak=0)\n",
    "amap113 = dicom.get_attenuation_map_from_CT_slices(files_CT, file_NM, index_peak=3)\n",
    "collimator_name = 'SY-ME'\n",
    "intrinsic_resolution=0.38 #mm at 140 keV\n",
    "psf_meta208 = dicom.get_psfmeta_from_scanner_params(\n",
    "    collimator_name,\n",
    "    energy_keV =  208,\n",
    "    intrinsic_resolution_140keV=intrinsic_resolution\n",
    ")\n",
    "psf_meta113 = dicom.get_psfmeta_from_scanner_params(\n",
    "    collimator_name,\n",
    "    energy_keV = 113,\n",
    "    intrinsic_resolution_140keV=intrinsic_resolution\n",
    ")\n",
    "att_transform208 = SPECTAttenuationTransform(attenuation_map=amap208)\n",
    "att_transform113 = SPECTAttenuationTransform(attenuation_map=amap113)\n",
    "psf_transform208 = SPECTPSFTransform(psf_meta=psf_meta208)\n",
    "psf_transform113 = SPECTPSFTransform(psf_meta=psf_meta113)\n",
    "system_matrix208 = SPECTSystemMatrix(\n",
    "    obj2obj_transforms=[att_transform208, psf_transform208],\n",
    "    proj2proj_transforms=[],\n",
    "    object_meta=object_meta,\n",
    "    proj_meta=proj_meta,\n",
    ")\n",
    "system_matrix113 = SPECTSystemMatrix(\n",
    "    obj2obj_transforms=[att_transform113, psf_transform113],\n",
    "    proj2proj_transforms=[],\n",
    "    object_meta=object_meta,\n",
    "    proj_meta=proj_meta,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db0730c9",
   "metadata": {
    "papermill": {
     "duration": 0.002893,
     "end_time": "2024-11-19T22:25:28.335255",
     "exception": false,
     "start_time": "2024-11-19T22:25:28.332362",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "From these two system matrices, we build a single system matrix that simulatenously projects to both peaks at once\n",
    "\n",
    "* The system matrices need to be scaled by the relative calibration factor in CPS/MBq in order to compensate for differnces in count rate between the two peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fd270bfe",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:28.342514Z",
     "iopub.status.busy": "2024-11-19T22:25:28.342274Z",
     "iopub.status.idle": "2024-11-19T22:25:28.346055Z",
     "shell.execute_reply": "2024-11-19T22:25:28.345476Z"
    },
    "papermill": {
     "duration": 0.009064,
     "end_time": "2024-11-19T22:25:28.347252",
     "exception": false,
     "start_time": "2024-11-19T22:25:28.338188",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# These calibration factors are approximate and come from SIMIND\n",
    "calib208 = 10.2419 # CPS/MBq\n",
    "calib113 = 11.3 # CPS/MBq\n",
    "\n",
    "system_matrix = ExtendedSystemMatrix(\n",
    "    system_matrices= [\n",
    "        calib113*system_matrix113,\n",
    "        calib208*system_matrix208\n",
    "    ] \n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b31911a4",
   "metadata": {
    "papermill": {
     "duration": 0.002881,
     "end_time": "2024-11-19T22:25:28.354020",
     "exception": false,
     "start_time": "2024-11-19T22:25:28.351139",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now if we create an object and forward project, it will project to both peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e5b3a753",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:28.361391Z",
     "iopub.status.busy": "2024-11-19T22:25:28.360845Z",
     "iopub.status.idle": "2024-11-19T22:25:28.384306Z",
     "shell.execute_reply": "2024-11-19T22:25:28.383655Z"
    },
    "papermill": {
     "duration": 0.02874,
     "end_time": "2024-11-19T22:25:28.385589",
     "exception": false,
     "start_time": "2024-11-19T22:25:28.356849",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "test_object = torch.zeros((128,128,128)).to(pytomography.device)\n",
    "test_object[50:70,50:70,50:70] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a8f500ef",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:28.393096Z",
     "iopub.status.busy": "2024-11-19T22:25:28.392506Z",
     "iopub.status.idle": "2024-11-19T22:25:29.712722Z",
     "shell.execute_reply": "2024-11-19T22:25:29.712092Z"
    },
    "papermill": {
     "duration": 1.325524,
     "end_time": "2024-11-19T22:25:29.714151",
     "exception": false,
     "start_time": "2024-11-19T22:25:28.388627",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([2, 96, 128, 128])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_projections = system_matrix.forward(test_object)\n",
    "test_projections.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4c6576b9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:29.725319Z",
     "iopub.status.busy": "2024-11-19T22:25:29.724692Z",
     "iopub.status.idle": "2024-11-19T22:25:29.799367Z",
     "shell.execute_reply": "2024-11-19T22:25:29.798772Z"
    },
    "papermill": {
     "duration": 0.083189,
     "end_time": "2024-11-19T22:25:29.800870",
     "exception": false,
     "start_time": "2024-11-19T22:25:29.717681",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAADiCAYAAACfgLRQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIEJJREFUeJzt3Xl8lPW1+PEzS5LJRhAEQQTUFuqr1QqCBXFpFdyq1uutiguKWmtr5daK1HrduGr784e4L607uLxYxOsCxY3gcgGxFFGoCooWCcgalpA9mZlz/6CZS/ieh8wQIN8kn/fr9fxz5vssM89z5swkZ75PSFVVAABAiwq39AEAAAAKMgAAXqAgAwDgAQoyAAAeoCADAOABCjIAAB6gIAMA4AEKMgAAHoimOzAUSnsoAADYiWp8l4/zDRkAAA9QkAEA8AAFGQAAD1CQAQDwAAUZAAAPUJABAPDAHi3Iffv2lVGjrpEJE56WJUs+lvr6GlGNy80337TL9Q466CC56qpfyuOP/0UWLvyb1NRUimpcnnzy8V2ud8wxx8ijjz4sH3wwV1avXinV1RVSUVEmn322RB566AHp3bt3xs9h7NjbRDUuY8felvG66bj88stENS6rV6+UcLjpl79fv36iGpf6+hrp3r37XjkmtH3RaFROOukkufvucbJgwXzZsqVU6uqqZe3a1fLaa6/IT3/60ya3MXToUJk5c4Zs3LhOqqrKZenST+WPf7xT8vPzA9eJRCJy9dW/ljlz3pfNmzdKXV21bNy4ToqL35ZLL71EQqFQRs9j5MhLRTUuEyY8ndF66TrxxBNFNS5VVeVSVFTU5PguXbpIbW2VqMbl6KOP3ivHhHZE0yQSaXK5//4HzHVvvvnWXa537bXXmes9+eRTu1zvzjv/pKqq33zzjc6e/Y5OmjRZZ858XdeuXauqquXl5frjH5+U1rE3LGPH3q6qqmPH3p7Reuku+fkddNu2baqqevrpZzY5/qGHHlZV1enTZ+yV42FpH8vQoaek8mrNmjU6Y8ZfdcqUqbpkyZJU/LHHHg9c/3e/G62qqolEQt97732dOvVFXbNmjaqqLl26VDt37uqsk52dq++9976qqtbU1Ghx8WydPHmKzp07TxOJhKqqvvzyKxk9j5EjL1dV1QkTJu611+rrr79WVdWrr76mybGjR49RVdUlS5a0+Dlm8X9pss42OaJhYBo7+8Uvfql3332PXnjhxfq9731fn332OVVtuiD/7Gf/pg8++JCOHHm5HnFEv1ShbaogH3bYD7R370OdeFZWLPXhoKSkRMPhrLRfsL1dkEUi+uSTT6mq6rRpL+1yXHZ2rpaWlqqq6tlnn9PiFxNL611OPHGYTpv2kh533I+dx84//wKtr69XVdVLLhnpPN6v3wBNJBJaX1+vp512Riqem1ugs2YVB17LDUV8xYoV2rPnwY0eGzDgR1pWVqaqqsOHX5j289gXBfnmm29VVdW///3vTY79xz/+oaqq1157XYufYxb/l31akHdeJkyYmFZB3nlpKIpNFeRdLdFojlZVVamq6uGHH5nxvvdmQR48+FhV3f6toVOnLoHjzjtvuKqqrlu3TiOR7Ba/mFja7tLwIXHWrGLnsalTX1RV1SeeeNJ5rFevQzQej6uq6ve+9/1Gj02fPkNVVceMucHc5+OPP6Gqqg888GDax7kvCnKPHr1SH1B29d5x9NGDU3ls/YWAhWXnpSlttqlLVSWZTIqISG1t7R7Z5oABA2TNmlUSj9fK6NHXNXqsT58+8thjf5avvvpCqqsrZOvWTfL+++/KxRdf5Gznww8/lM8++0xycnJkxIiLA/d3xRWXi4jIc8+9IIlEYo88B8Dy8cefiIhIz54HNYpnZWXJGWds///ypElTnPVKSkpk3rwPRETknHP+rdFjNTU1ae27tHRThkdrO+SQQ2Tp0k9FNS733Xdvo/9Pd+/eXe699x75/PN/SGXlNtm2bYssWDBfrrnmNxKJRBpt59tvv5W33npbRP4vBy1XXHGZiIhMnz5DNm3aM88B7VubLMjhcFjGjr1N8vPz5bPPPpOvvvqq2ds866wz5f3335GioiI577zhct9996ceO/fcn8vixYvkV7+6Surq6uT119+QhQs/kqOO6i8vvPCcPP30k872nn56goiIXH75SHN/PXr0kJNPHvavsc80+/iBXenT57siIrJ27bpG8b59+6aathYuXGiuu3DhRyIi0r9/v0bxN954U0RErrnmaunZs2ejx4466ii54ILhUlVVJc8//0Kzj3/QoEHy4YfzpE+fPjJq1G9l9OjrRVVFROT444+XTz9dLKNH/05isZjMmlUs8+Z9IN/5znfkkUcekpkzZ0g02niu/oacGzHiIucxEZFYLCYXXDD8X2MnNPv4ARHZs01dOy/76k/WPXserBMmTNQJEybqa69N15KSElVV/fLLL/X73z9it/a945+sR436rcbjcV2/fr0OGjSk0fjDDz9Sq6urtaqqSs855+eNHuvV6xBdvHixqrr/m9t//wO0trZWVVX79x/oHMdNN92iqqpz585r8T+zsLTt5YADDtQtW7aoquqoUb9t9NiZZ56tqqqbN28OXL/hf8ULFixoFA+Fojpx4rOpP+sWF8/WSZMm65w5czWRSOgnn3yigwcfm9GxWn+y/vd/P1crKyu1oqJCzzrrbOe5bdy4UROJhP7617/RUCiaeqxTpy5aXDxbVVVvvXVso/Wi0Rxdv369qqqT1yIRveiiEaqqunLlykbbZGHZ1dJknW0LBfnII49yjnfhwoX6wx/2z/iYdyzIoVBU77tve3PYsmXL9NBD+zjjJ0+eoqqqo0ePMbc3cOAgVbUbRKZNe0lVVR9++BHnsS+//FJVVa+44soWv4hY2u4SiWSnGrMWL16sWVmxRo9feOHFqqq6atWqwG1ceeVVqRyxHh89ekzqw2eDiooKveee+7RLl24ZHe/OBfn663+viURC165dqwMG/MgZf9dd41RV9aGHHja3d+CBPbW2tlbXr1/vPDZ+/L2qqjpjxl+dxxoK+R13/LHFzyFL61ma0ib+ZL148WIJhaISCkWlR49ect55wyUvL08++miB/Md/jNqtbebl5cl///c0ue66a2XOnLlyzDHHyT//+c9GY0KhkJx++mkiIjJ16ovmdhYuXCjl5eXSv39/ycnJafTYU09t/7PYRRddKNnZ2an4CSecIH369JHy8vLA7QJ7wmOP/VmGDRsqpaWlcu65w6W+vn6PbbuwsFD++tfpMn78OHnkkT9Lnz6HSV5eoRx++JHy6quvyfXXXycLFsyXgw46qOmN7SQSicijjz4s99xztyxbtkwGDz5WPvroI2fcGWecLiIiU6dOM7ezZs0aWb58uXTt2lX69OnT6LGnntr+W+fTTjtVunXrlor37t1bTjzxJ5JMJmXChIkZHzsQqMmS/S+782mgJbusi4o66erVqzUej2f0Tblh33V1daq6/feF2dm55tj99z8g3ZdPVVUPPLBno/VDoaiuXLlSVRv/9KPhz3xPPfV0i3+iY2m7ywMPPKiqqps2bdJ+/QaYY5rzJ+snnnhSVVUfeeRRc7033nhTVVUnTnw27WNu+IbckJ/r1q3Tjh07B46vqKhIOz+HDDneWX/u3HmqqvqHP/xnKvZf/3WHqqoWF89u8XPI0rqWprjdCm1EWVmZvPLKqzJq1DXys5+dJUuWLMlo/ZkzX5fjjjtWjjjiCLnxxj/IHXfc6YzZcaatiROfbXKbO3d7q6pMnPic3HbbLXL55ZfJ1KkvSkFBgZx77s9FhGYR7D333DNerr32t7JlyxY55ZTT5ZNPPjHHffPNNyIist9++0lBQYFUVFQ4Yxoatr75ZmUqFg6H5ZJLRoiIyOTJU81tT5o0WU477VQZNmxoxsc/Z85cOfjg3nLooYfK+PHj5Kqrfp1q4tpRQ45Om/aSVFZW7nKbVqf0008/I8ceO0Quu+xSGTfubhERGTnyktRjwJ7UZguyiEhlZZWIiHTt2iXjdT/5ZLHcfPOtMmvWm3L77WOlsLBQfv/7GxqNKS0tlaqqKsnLy5MxY27YrZ8+PPPMBLnllptk2LChctBBB8mpp54i+fn5snTpUpk/f37G2wOaMm7c/5frr79Otm7dKqeccrr5p94GX3zxhVRWVkp+fr4MHDhQ3nvvPWfMwIEDRERk0aKPU7GuXbtKLBYTEZFt27aZ2y4rKxMRkU6dOmX8HEpKSmTEiEuluPgtufLKX0hBQYGMGHGp8/PAVatWSd++fWXcuPG7fJ5BXnxxmjz44P1y2GGHyZAhQyQ3N1cOPvhg2bJli7z88isZbw/YlTbxP+QgJ530ExER+fLL5bu1/ueffy7HH/8TWbFihYwZM1r+8pdHG/22MZlMyqxZxSIicv755+3WPlauXCmzZ78jkUhELrtsZOp3j888M3G3tgfsyl13/T+54YYxsnXrVjn55NMCf8rUoL6+XmbOfF1ERC666ALn8V69esmQIceIiMgrr7yaim/atCn1O+RBg35kbnvw4MEiIrJixYqMn4eIyNq1a+WEE06URYsWyQUXDJeXX36pUS+GyP/99Gp387OyslKmTNn+Df+KKy5L/fZ40qTJe2x+AyAl3f+v7M7fy/f2/5BvvPEm3X//A5x4x46dU3NAb9myRbt27Z7xvnf82VOPHr106dKlqqr6/PMvNJo1q3//gVpTU6MVFRV66aWXmT+B+MEPfmj+dKJhGT78QlVV3bBhg6pu//9YJsfMwpLO0jAl7ebNm3XgwEFpr9e//8DU1JmnnvrTVLypqTMbZvhau3atHnFEv0aPnXTSyamZ9DJ5f7B+9tShw346Z85cVd0+01heXmHqsR49eunmzZu1vr5eR48e43SRi0T04IO/oxdffEngPgcNGqKqqtu2bUsds/VTRRaWppYm6+yeLMj9+w/U+fM/TC0NBaakpKRRvFu3Ho3W69atR6PHV61apaqq69evbxTfOQlUVevr63XRokU6bdpLOnnyFP2f/5mj5eXlqWI8bNipGb1gQVNndunSTT/++GNV3T4h/o6NXueee36qeaSkpETffPMtff75F3TmzNdTv4mePHlK4D53nLO6YfstfeGwtK3lrLPOTl1fCxYsSP1uf+dl/Ph7zfV3vLnEO++8q1OmTNVvv/1WVYNvLnHggT31q6++SuXpnDlzdcqUqbpgwYLUsRQXz9acnLy0n0fQ1Jl5eYX69tuzVFV13rwPtKioU+qx44//Seq9aN26dVpcPFuff/4FnT59hi5fvlxVVefP/3CX+/30009Tx7xo0aIWP58srXNpss42OaJhYBo7+/GPT0prWzvfEKJ370PTWm/nOzf95jejdPLkKbps2TLdvHmz1tXV6aZNm/SDD+br2LG379a3zF3NZV1U1EnnzftAVVXfeuttzc0taPQc7r33fl2yZImWl5drVVWVrlixQt9551294YYbzd8w77g8+OBDqed55plnZ3zcLCy7WhoKWVNWrFgRuI2hQ0/R119/Q0tLS7W6ulq/+OIL/dOf7tKCgqLAdQoKivSmm27Rv/1tgW7dulXr6+u1tLRU3333Pf3lL3+V0Y1fdnwe1lzW2dm5+sorr6rq9qK541/PunTpprfffqcuXLhQy8rKtKamRktKSnTu3Hk6duztTc53f91116deo50nT2FhSXdpSkjVaE00hEJtuv8LAIC9SjW+y8fbdFMXAACtBQUZAAAPUJABAPAABRkAAA9QkAEA8AAFGQAAD1CQAQDwAAUZAAAPUJABAPAABRkAAA8wH+ZeEWp6SKvn22e5pBFLa1ZYtHmtNR99yzGRkPFaatp5ZuVokPaZu/6dcQAA2iEKMgAAHqAgAwDgAQoyAAAeoKkrkNUI4n5+sZocJGR9ztnzn31C5n6aK71tprvvUND2mnPs6jaHJI37jKrW2atrwohaMfgt4kTSzcf07+++b76zNDeXA/PMGWi9h6W/b+s41chHNRq4kkk3H8nRxviGDACAByjIAAB4gIIMAIAHKMgAAHiApq6AWXxCoSwjlu3EIuGYE4tG3FjIakAJaORIu2GqGY1Vwft2jzNsNMCEjW2GjdfMiomIRIx4ROyxO4tLjROrjm91Y3Wl5vp1xthkstYY2fabSPyTXjOliJ2jkUieE8syYtFIrrsX4zq38iFIus1RVu5Z+RS073Tz2cqxsJFj1rig/YSN9zFLnVY5sarEJicWlKO19VudWDJZbYxsWznKN2QAADxAQQYAwAMUZAAAPEBBBgDAA+2sqcttGLEaQ0REImG3ESQ7WuTE8nO6OrG8SGcnliMFTiwa0MRkNlNomo0gRtNF1IhF1D711tgs4zLJMWJRo7EkFrGbQGIRd2xW2Dg/xrrVCffWbOu00omtyF9m7ntj1VJ3m3XrnZg1A1F7vS3cvmM1LNl5Eo24OZWX08WJFWX1dGNygLuuWjmafmNV2LharVgkzXUjAY2XWUbcjBn5lB0xYsa4oLFGSNRIifJ6N7i2xm3KWp7n5qKIyIaqz5xYe8hRviEDAOABCjIAAB6gIAMA4AEKMgAAHmhnTV3p34YtaszuUxg70Il1C/V1Y0m3saRD1J3ly2q6EBEx+p3M28pZDRZRY5tRY3tBjRw5Rg9LtrF+rjEuL+o2UxQaMRGR/IjbjFEQdW+hGA6562+pc8/ZlxUdnFis9Ahz3zWxMidWF3dj8US9sXbbmhmoZVlNfEZjVNjNHRGRWPZ+Tqxz9ned2CHJPk6sR8zN74457oVuXfsidu5ZMSvNsoxtRtMcJyKSE3ZzIhYxYsa4vKibd/kR+5rOM/IxZoxNqnvwG2rc2Qo/3ZbvxLI2HG7uuzq2xYm1hxzlGzIAAB6gIAMA4AEKMgAAHqAgAwDgAQoyAAAeaMNd1ul1cAZ1WWdF3Y7AgrA7TWZ3dWO98t0Ow/1y3H3HAm4tmm4Hp9U9HTW6ko3mUYkZXc4iIrGwG0+3M7Mou86Nxaz7DIt0yHfvaZxb6K4fMp73ts3uvWwL1rmd7Ruq3fMgIrK8yj1nW8MrnFgiUeHEWu+kfK2Eda9g4z7kIiI5kUIn1jnp/hKiV67bUf3dDu5+uuQY135Antj56F4dETMf3W1aeZsb0P0ci7rx3KjbbZyXY8Ry3RzL7WB1Kotkd3SPM5Lvvm4ad4+9ssR9c4v9s7sTW1sVkKOV7tj2kKN8QwYAwAMUZAAAPEBBBgDAAxRkAAA80IabutJj3ddURCQSznFiMeOexoUR916tHY0Grs7u5szp7kTshpEso+nDagTJNqfVc5szcgOaVaxGkoIsdwq9AqOBq0Oe26iV39Fu6op1MZ7PfsblaMwpmL3avffxgdvcc1OUHTDlYpU7NhRw71m0vHDQ9LbGPcvzxW34s6bEtBq4useMJiijgUpEJBJy188yci/LauAyYjFjmsocI+9ERGJGs1Ys1x2bXeAee1Ynd3uRznaehDu7r2Wo0I1J3N1POHeTE+ux3s3bomxjeyISq3CbattDjrb9ZwgAQCtAQQYAwAMUZAAAPEBBBgDAA+2+qcuaGUjEbvaKqvtyWfc0tu6hat3D1IqJ2A1cZqOXsX620TBiNXAFztRlNHVZswDFso1YnhvLLgiY6ajQncknVOQ2l4SM6ciiZe5+rGO07u0sYp/H9tAw0jqkfx7CxjnLUvek2/fztpof3WvfionYzVpWo1e2sX7UyD2rgSsnO6CpK80Grqh7i3CJFLnXfng/e7as0H5uY5V0cBsipd7NvXDHcieWazyfoByNiNss2x5ytO0/QwAAWgEKMgAAHqAgAwDgAQoyAAAeoKkrAyExGpGMZiujz8scZ4Qy26axbrrjrAYUEfsWctbxRIxxYaNRJmCSJQkZnTahbKPDw3hC1rrW8QS+vsZ5hL+Cmnms8xg2zrqdE+ndKjEoT6y4eQ2auZPeuuGAps+QcTvUkJl76d6zNeB7WZbbWCXRNHMnYtxG03jeQTnaXvENGQAAD1CQAQDwAAUZAAAPUJABAPAATV0B9sWsMFazSRCrCSVsNKFYR22Na+5+QuluM5OX0XpBMnmRAM+YOZrmJR0KaOpKl/UWFjJzLCBJm5GP1n4yeR9qr/iGDACAByjIAAB4gIIMAIAHKMgAAHigfTV1tYPbd7UaNGuhBaV7+WXUeNmMpqXmrCvSwm9t5PIeQ4UCAMADFGQAADxAQQYAwAMUZAAAPEBBBgDAA+2ry3ofMe9TnMH6++JTUib7SH+azPTuB5uRNG/w3NwuVbRdba0HuNk51RzmNJsJY5x1kLs/jWh7wTdkAAA8QEEGAMADFGQAADxAQQYAwAM0dWG3cG9TtDUZ3Te8NTQjmfczNsYFPRk6rvY5viEDAOABCjIAAB6gIAMA4AEKMgAAHqCpqxUzJ7FqLc1WZsMJTSRAg1bTOGnO3oXdwSsJAIAHKMgAAHiAggwAgAcoyAAAeICmrmZqTh9S2Lgd2d7YTyaretdIYjWMpPli0CPW/oQ8m0LLarLcK42XreCrlWenxkut4DQCAND2UZABAPAABRkAAA9QkAEA8ABNXRkI8/lll8ymDTqrgLYj3XxuZt631/fa9vmsAQDwDAUZAAAPUJABAPAABRkAAA/Q1NVOZTJbkDV7VyjcvNmGQnv49ovezTAGr6V7pYUymU0vg7HOuhlcv6F0v0ZZ4zLJO26ruM/xigMA4AEKMgAAHqAgAwDgAQoyAAAeoCADAOABuqz3AqtpMZN7gTLb5A6MF87s0AaaybdvJ63m/sHk4x7j2zUIAEC7REEGAMADFGQAADxAQQYAwAM0dbUxzf2E1ZxGEmtKv7Sn+QP2kr3Rc+TbVK37LM/MHSXckPGiW9P10g7WGG+XAAB4gIIMAIAHKMgAAHiAggwAgAdo6gKADCXVbUeKtGCHkiY921HSbeBS4zVDY3xDBgDAAxRkAAA8QEEGAMADFGQAADxAU1cGkpJeQ4PRzyBqxJIB89QkjcHWSGs/CaNxwppTKKjBIpF041YDi1rjjAl7NG7PaKRx97UMGbF017WOEe2PGrlj5qOxbtzMHfu6MvdjjLXyzIpZ16+ViyIiESPPzLcmK/esWCIg7+LGjtKNmU1d7rCg+c7Sfa9ta/iGDACAByjIAAB4gIIMAIAHKMgAAHigfTV1ZTCdjRpjE1LvxOqNTdYZsVojFgq412HCCGcZMbvhxHqOEXff5p7tW6RF6rPcWNjdT6TSjWXF6sz9hPPi7r6jte7AqNHsUu7upzbuXsq1VvOL2OcR/rJyUUREjdv+xY3r38rHGiPJahJungTdZjEr7H6XiSTdHcWNPLGataxYJjNbhUJuPomxb4m6+RjJDvheFnHjIauBK+HGkmVuLlfX5Tkx69yI2Oe2PeAbMgAAHqAgAwDgAQoyAAAeoCADAOCB9tXUZVA1miFEJJF0mxJqQpVOrLzebRDaWmc0QxhtVDG3h0RERCJGs1fU+OgUNRpOrGaTnLA7Lha2dx6LuJdEntF5UVjnNnp1qM1xYkXVNeZ+Cre68VgH9zUPGc+7vDTbiX1b6TaMlNXZDTm1oWonlkza1wH2NWsWNvvc1Cfd81hhnNuyunwntrHOaHQ0eqhi9fZbZNTIqYiVj2aOus8x24jFzCm5RGJRN56b5b4P5W02Yhvcpq7YhipzP1kd3dcynG+8byTc51i10o2tqShwYptr7RytNt5r20OO8g0ZAAAPUJABAPAABRkAAA9QkAEA8AAFGQAAD7ThLmvjfpxWLKiDM2F0VCc3OLG1oU7uNivc2Oao+1JnR+yp8ayObGtoNJzeuGxjCjxrKs7tY90uSqsbPDfidlTnRdzXtzBqz41XmOV2iuYZXaXWYW41Ol+Xlbtd3ysq7A7vbbrOiSXU7T61rhfsZcY0mUnj3IiI1MS3OrHSnDVOrKTS7bKOJ2NObF2Oe6EHzippXJhGOprjsoxtWjHrVxQiIjnG2JiRezGjc7vAyMf8qP0emGd0c8ci7ljrXs4batzXd0mZG1tZYUyXKyLbkkaOJt18bms5yjdkAAA8QEEGAMADFGQAADxAQQYAwANtuKnL4jY0BDZ1xSucWHntt05sdbY7Pd3mSGcnFku408ZFEm4jkohIWN3PSWHjs1PIuM+xNS6i7rhowKnPMuJZxn6yQm4sx4jFjIYyEbt5LDvs7jtszGdYFXfP47pa93x9HV5m7ntb9SonlkhY0wemf/9s7I70Gi+TSbupq7a+zIltCn/txOJZbuPQupoDnFhetZGjxrUvIhI22g3TjUWsHDXmiI0E3LU8aoy1GjytaXStHMuJ2O9DRo+bOa2v1VZVbky3u6bazbHl4S/MfVvvtVZTV1vLUb4hAwDgAQoyAAAeoCADAOABCjIAAB4IqWpaU52EQm21/8tu2ggZDUrhsDvTTMSMuffrDYfcxomQdbPfAOmODVnNX8a61vPLZP2I8XzCxjVijRMRCYuxfsC52Flc3OaOmsQ2J1ZdV2quX1u/1YkljXvritj3o8XeZDUy2dd+yLi2IhH3vthZRiwayXX3Yly/meSJOc7IHavx0l43oKHMOk4jd6zcs447KEcjRo6m+7zr1c2nqoSbj5V17uyHIm03R4OaiBvwDRkAAA9QkAEA8AAFGQAAD1CQAQDwAE1dgdJrLrFulShmA9ae/+yTSVOYq3nHYzaKZbLNNNdXayYe8xZ9brOEBty2T9VqBGldzSHtS8C9QtPOR6tZK91r1be8zSDPmtEIKtK847Ty0Zpxrb3lKE1dAAC0AhRkAAA8QEEGAMADFGQAADxAU9deEdSE0pb49lkuk9uwpXXJo81oDfnoWz4FNMgZrFtm2shRmroAAGgFKMgAAHiAggwAgAcoyAAAeICCDACAB2id3ivaZodgY61/Gju0F60hH/3Lp9bwqrU1fEMGAMADFGQAADxAQQYAwAMUZAAAPEBBBgDAAxRkAAA8QEEGAMADFGQAADxAQQYAwAMUZAAAPEBBBgDAAxRkAAA8QEEGAMADFGQAADxAQQYAwAMUZAAAPEBBBgDAAxRkAAA8QEEGAMADFGQAADxAQQYAwAMUZAAAPEBBBgDAAxRkAAA8QEEGAMADFGQAADxAQQYAwAMhVdWWPggAANo7viEDAOABCjIAAB6gIAMA4AEKMgAAHqAgAwDgAQoyAAAeoCADAOABCjIAAB6gIAMA4IH/BXlvMFjChB+LAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "IDX_113 = 0\n",
    "IDX_208 = 1\n",
    "test_projection_113 = test_projections[IDX_113][25].cpu().T\n",
    "test_projection_208 = test_projections[IDX_208][25].cpu().T\n",
    "fig, ax = plt.subplots(1,2, figsize=(6,3), gridspec_kw={'wspace':0.0})\n",
    "ax[0].imshow(test_projection_113, cmap='magma', vmax=80)\n",
    "ax[0].set_xlim(45,90)\n",
    "ax[0].set_ylim(40,80)\n",
    "ax[0].axis('off')\n",
    "ax[0].text(0.03, 0.97, '113 keV', color='white', transform=ax[0].transAxes, ha='left', va='top', fontsize=16)\n",
    "ax[1].imshow(test_projection_208, cmap='magma', vmax=80)\n",
    "ax[1].set_xlim(45,90)\n",
    "ax[1].set_ylim(40,80)\n",
    "ax[1].axis('off')\n",
    "ax[1].text(0.03, 0.97, '208 keV', color='white', transform=ax[1].transAxes, ha='left', va='top', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e79464ab",
   "metadata": {
    "papermill": {
     "duration": 0.00328,
     "end_time": "2024-11-19T22:25:29.807503",
     "exception": false,
     "start_time": "2024-11-19T22:25:29.804223",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We can create the likelihood by stacking the photopeak/scatters as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "eceda3a1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:29.817544Z",
     "iopub.status.busy": "2024-11-19T22:25:29.816924Z",
     "iopub.status.idle": "2024-11-19T22:25:29.821226Z",
     "shell.execute_reply": "2024-11-19T22:25:29.820627Z"
    },
    "papermill": {
     "duration": 0.010203,
     "end_time": "2024-11-19T22:25:29.822461",
     "exception": false,
     "start_time": "2024-11-19T22:25:29.812258",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "photopeak = torch.stack([photopeak113, photopeak208], dim=0)\n",
    "scatter = torch.stack([scatter113, scatter208], dim=0)\n",
    "likelihood = PoissonLogLikelihood(system_matrix, photopeak, scatter)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2afa2e2",
   "metadata": {
    "papermill": {
     "duration": 0.003309,
     "end_time": "2024-11-19T22:25:29.829093",
     "exception": false,
     "start_time": "2024-11-19T22:25:29.825784",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Then we reconstruct as normal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b79df128",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:29.838562Z",
     "iopub.status.busy": "2024-11-19T22:25:29.838067Z",
     "iopub.status.idle": "2024-11-19T22:25:56.775065Z",
     "shell.execute_reply": "2024-11-19T22:25:56.773705Z"
    },
    "papermill": {
     "duration": 26.942928,
     "end_time": "2024-11-19T22:25:56.777040",
     "exception": false,
     "start_time": "2024-11-19T22:25:29.834112",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "recon_algo = OSEM(likelihood)\n",
    "recon_OSEM = recon_algo(n_iters = 8, n_subsets=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "041e4dfa",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:56.787444Z",
     "iopub.status.busy": "2024-11-19T22:25:56.786939Z",
     "iopub.status.idle": "2024-11-19T22:25:56.942580Z",
     "shell.execute_reply": "2024-11-19T22:25:56.941955Z"
    },
    "papermill": {
     "duration": 0.161889,
     "end_time": "2024-11-19T22:25:56.943983",
     "exception": false,
     "start_time": "2024-11-19T22:25:56.782094",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAGOCAYAAAC3yitIAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAr0RJREFUeJzt/W2sLclZHgxfVd1r731mPGeIH+IZY4ZgiQhCADvhw4wjBSM5GVkIxYqE/PDHlvl4FQkiHEeKQhRBlD/+kTg4EiiOhXgsIiETQuJIYJE4tgARHCGDLcVIWCLyG1uImYCeZAbPnL1Xd1c9P6qqu6r6rurqj7X2Wnvf19E6e63u6qrqr7rq/iyhtdZgMBgMBoNxspC33QEGg8FgMBh5MFkzGAwGg3HiYLJmMBgMBuPEwWTNYDAYDMaJg8mawWAwGIwTB5M1g8FgMBgnDiZrBoPBYDBOHEzWDAaDwWCcOOrb7gCDwWAwGEtwfX2N/X6/up6LiwtcXV1t0KPDgcmawWAwGGeH6+trvP71r8Pzz/+/q+t6+umn8YUvfOGkCZvJmsFgMBhnh/1+j+ef/3/x///CR/Dw4WOL63nppVfwta//v7Hf75msGQwGg8E4BB4+fAwPHz5+2904OJisGQwGg3G+UMp81hx/BmCyZjAYDMb54p6QNYduMRgMBoNx4mDJmsFgMBjnC63NZ83xZwAmawaDwWCcL5ReqQY/D7JmNTiDwWAwGCcOlqwZDAaDcb64Jw5mTNYMBoPBOF8wWTMYDAaDceK4J2TNNmsGg8FgME4cLFkzGAwG43yhV0rW+jwkayZrBoPBYJwthFYQKwh3zbHHRDFZC8G8zmAwGIwyaN3edhfuFJiBGQwGg3G+uCcOZkzWDAaDwThfKL0uCxlnMGMwGAwGg7EFWLJmMBgMxvmC1eAMBoPBYJw47glZsxqcwWAwGIwTB0vWDAaDwThfaL0usQmvZ81gMBgMxoFxT9TgTNYMBoPBOF9w6BaDwWAwGIxTAEvWDAaDwThfsBqcwWAwGIwTxz1ZdYvV4AwGg8FgnDhYsmYwGAzG2UIoBbFCsl5z7DHBZM1gMBiM84XW62KlzyTOmtXgDAaDwWCcOFiyZjAYDMb5gr3BGQwGg8E4cdwTsmY1OIPBYDAYJw6WrBkMBoNxvrgn6UaZrBkMBoNxvrgnanAmawaDwWCcL5ReSdbnIVmzzZrBYDAYjBMHS9aMOwRh/5bOlMV0kR5+nXOOi49lMBib4p4kRWGyZpwB5pLj3PKHrnPuJILBYBSDbdYMxm3jEKSbQ2wV2volZtJmMBjLwGTNOEEsIemU+0VMuHPcNLZw6aAIn0mbwdgMemXoFqvBGYy5SJH0GtIsO1ZkJgg6QarUMeOyrv0caa/FeQw2DMZBwGpwBuNYoEirnKBzRFteSaI9rWbVnyZwv/5DqdeHFhkMxt0CkzXjFpEn6VkknCLbImSOFTnJWCa2e/DIfpC6U7bxpecQ94HV7Ix7BJasGYxDYVrdHRB1MRHPJzuxum4JrVMvuxr6TpL2VN3pCUtYR0pqZ9Jm3ANwulEGYy1KJeOlJE3vX0/Ac6DI9gyBe1K5kIAldf8cfdJdo24fS+1M2gzGXQKTNWNDzLUdhyTXE1BPftNkmibm/LHlhJ5GSMg+BgIPyviXhyDuoXMz+xbZ1QcbeUo9bkoxGHcCWvXv0+LjzwBM1oyVmEPQCUmYlKblagk5dbwokaipY+OXOnHqfiSIEG4bIWlnMYOwo0mAgIic2pZ4ojOZM84ErAZnMHJYH2Y1kipHRG1+h1JqdMgcQvbKFhG214bWA8FqOKlY0rNy77QccQti21Sb5fCujzAStbD9Wu6Jzh7mjDMBO5gxGCnEA/mYWIrtrwEp+eRsCNsnVZ/ANNQkIYf9CYnf3zbZRSEBMUwWNNTwndJik2nEB0c0ER2jNWX3jkwEGfIe6jVErfXQbkja43ppcCIXBuPUwGTNKESZJJ2WlqeQIOqZqmxqO0XQOfKL9wlU0OgCCV/DfGLu0pGKmyLuULL2iZqemKTOyy/rpH2tzWTAtDE4tqUd0ijkJHEmbcaJgdXgDAaQt22uCbVKkWpI1HPU1WE90TaPpN1+gaq4HhGomu1GT4J1++L+9pK3vy+hEs9K/RPXtJ8UiJiww/ayjm2uCBkWxqTNOFHck/WsmawZCawh6XkES0q5HlEv9dyOCc8RtSPpUul6ZDP3SVq7P2NVtk/evrpcZI4ZCL1csvbrF7ZvjrBNP1ypCce2bGjZFGmb0gwG4zBgsmZEKCNpUzIOtTJlRPR7wNTsN22nXiJh+/Alauo3VT5uVwgJLYy9Wil4au1Qou4d0Pz6PR4LiHXUcPp8s5MWNwEQIWGb41w/04eTRJ4MCbPlxz0cd4jBODRYDc64PygnaFOaJmkgVGWn6kpn/PLbkCPCLHGyGtVD2J8DVfikyjxUlWvdQaE1x0VqcNcP2vHNL2T/UM5pE+eatNU7u7UeCNttG0jbOKDRIDKxESFhpvslHuZM3IxjYWWc9ea5+g8DJut7jQVSdL8hVHnHoVakc5jntay9F8T3Zo6PExhLxDFob+qou1b93avCExLsQOIVSZwKxsNayhpKtQFh5/ohAm/wgVipflL9SpULYPsiIHsHOJ+0DdLpUcPwskRc+Kxc533NUScZDMZcrNMtMs4UAmmiHkjXlBRju3SGqKWoIdwHMvyIuj/e7XfHl9ilfcew4FN4/Ki+qH0nSTuiFkKiEjUkJCqxg4SEFDWk3PX7Yxv4lGrdJ+OR1B/ty33IOoV/rcNr7t8Pc4+oOodP2I73THj3X3j/QoTPUHzVGYxN4dTgaz4z8L73vQ/f/u3fjieeeAKvec1r8Pa3vx2f//znJ4/7pV/6JXzDN3wDrq6u8M3f/M342Mc+NqtdJut7hS1JOh7kabKRog7IQXqkAY9cJnseeXCHntN5spxSdQ/lqmCfREzeldkWkVkvtUf9mCLbmJT9bXE91DlRpE3W7V3vfrIUT6RGfSkk7Yi4Q6RIO/ccMhgzcWSy/o3f+A38yI/8CP7bf/tv+PjHP46mafA3/+bfxMsvv5w85rd/+7fx/d///fjBH/xBfOYzn8Hb3/52vP3tb8fnPve54naF1iX5lGBm6IwzxRbq7rCsP4AHqluPkP1yvup1iAl2scqxrTSWDtNe3GG9Hbk9kMi9SUOwz7NNuzZ7svbOW+kGnW6hdAOlW3MOejifuC9xeJjbHhznZUTL2dKpc/Lb8OsO/oLqX3TNE/dnALVt2B7XYeqhhpb4eFaL31Vo3R60/pdeeglPPvkk/s9H/j4ePna5vJ5XbvAV//f78eKLL+Lhw4ezj/+TP/kTvOY1r8Fv/MZv4K//9b9OlnnHO96Bl19+Gb/yK7/Sb/vO7/xOvPGNb8QHP/jBonZYsr6zEMhLMGOpJ3Ac66WmsZQVS17m2EhiiyRJKWtIWYfStUjYthPwybjvc+LYFKEPdUXEosd10/U6wq9ASdh9OU/SDo73tgVSMDkRmofYLBBL1oH6O5KqSambkKpzqvOUtB1OAOPnjqVsxkq4dKNrPivw4osvAgBe/epXJ8t86lOfwlvf+tZg23PPPYdPfepTxe2wuHznMDXwUSrh9GpXIXGMt/tkOyJpiqxENZBuH2LkLSGZIKoiJzJUJKHnUFKv6ffYMUsISXqET7aHajjvqL5Y3Q/QkxTX7+yEwN1Wv++eZ3gwYcmFlgXnJwkJf8ie5nXCqzDllBYv48lSNmMBNgrdeumll4LNl5eXuLzMS+xKKbznPe/BX/trfw3f9E3flCz3/PPP46mnngq2PfXUU3j++eeLu8mS9Z1BiRQ9lqQpoqYl6Gm7c6C2jiS8ULLeQYpdaE8N+hBiKtSrhGyNyr1LquPN3zKiH+zWg407tjOX9pWyS09lVkvVO9JmiF0vbbvr7e7BSLqO7N3mC23fjh3UcpqXWOIOJW1KymYwZmIjm/UzzzyDJ598sv+8733vm2z6R37kR/C5z30OH/nIRw59lixZ3w1Qg1yGMJJ26VjtGat2x6Qax0PHKuGcDdaFM1FSq0NywY7kuU1L1+HCF3SbRmLOq9JdHX1sdSaEK7ViWEnM+RoE/SAl5Sg2nMjMRk2i8jnPE+t5R7nKaSmbJWzG7eBLX/pSYLOekqp/9Ed/FL/yK7+C3/zN38RXf/VXZ8s+/fTTeOGFF4JtL7zwAp5++uni/rFkffaY9sAV0b9hh0wSNSV9hfZISRI15QgWEH6kJh9JcwXYkuBy0rWCSkrbvmNafA7kxCSWqEkTQfoaUE5sKVBSdtBnz7YdtzsKZyO8xk2B0MYd1jM8g2O7NwKb9gCWsBkLsZHN+uHDh8EnRdZaa/zoj/4o/sN/+A/45Cc/ide//vWTXXz22WfxiU98Itj28Y9/HM8++2zxabJkfbZIL1OZXZ5yRAjjgVaIaFDu6y0Z3KsROfvlHFxebAEJoWWRQOUW0RCoAul1KXn3Ur2V3gcpv7PtyL6vFFzyEXeeh5aS436bvnWgnNZi9BoA75r5tvNw+3DelIZglFpVSNK+TUva5ne8GtiQFU1FFbGUzZiA1lO5dKePn4Ef+ZEfwS/8wi/gP/7H/4gnnniitzs/+eSTePDgAQDgne98J173utf1qvQf+7Efw3d913fh/e9/P77ne74HH/nIR/DpT38aH/rQh4rbZcn6LEET9cgGTX368qH047b3RB3ZK52dUhLJNXo7aAFp+PtH0nWBZAnQ4U9zER9Hhjdp1UvX/kdF5B2nJC1rL62qD0LQEuFoVF/jbZRkTUrZhObD7CMk61hqB/28hHbw+FmLbNnwJ5jxM8BSNuO08K/+1b/Ciy++iLe85S147Wtf239+8Rd/sS/zxS9+EX/8x3/c/37zm9+MX/iFX8CHPvQhvOENb8C/+3f/Dh/96EezTmkxWLI+O6SJ2nyJVI6pWgJi9AZzz2nMLxdn3PIHe7rOiTY1AqlUaElKcjEGaTiUKEuh0QF66EucAtRJ2Gaj6WeX6VcsXfvn5feZ+p46t76fmWOmvNgpRzb/uNhXgLoedCfjhsLt1OpisU1b9J7wVpIe2bF9CZtohMHwceSFPEpSk/z6r//6aNv3fd/34fu+7/tmteWDyfpsUKD27qUW+3OC+OKyFFGTTmMZgp5uMyzr1M6O6ACC7EZq6oGw4/643yUhXCP1b4qwTSeS5+mcyyj7dilRl/Q1uW/GxCUm6mGC0fX75/TJTVb64wjHNIq0RyuBCRCEDTBpMyZxT1bdYjX4WWAeUY+delKf0JFoiqiTSTcideu49+Pto3oJ9WsOqcxlVBvkMVTWMSLjV6/+tgtjUHWVoKS/VL2Uuj7VD9LjPHYMQ3jth3LjCU/8cceMJ2/RcyOiNr1yw3GE6WOkEgfoIYoTqTDuH1iyPmnkQ7JSRD3yBp9STxPlSaKm1OAzJDGq3ViKDlSzxEBNScNrMSVhD9J/NSuUjM6clldvU8fHkwiy3zPNArEKvHc4y5XzNBakKp1wuOud1bxEK04Sd85m/ve045kDJWmbkox7Cr0yC9mRHEPXgsn6ZFHo7T0iUEpqmiCXyC6dI+otSHrUvFOBj5ZopEGT6zIbdq7+kUq877D744jSSeDKfqel41KiLpHa/UlDcRa2zOTIJ123Pf7uCH2KsPv+BcfJwHN8krCBgLRNnSXEzaR973BP1OBM1ieHuSSdVi1OZQYb6h2rvn37dIqoU5OAVKgT2bbfN2s7jZ2dyDYWSNWUWjwlpfuEDSApZef6F/ydIVHPAdVvh7XXx22jJPBYlT8ieyr8y3NCGxM2rOOZqyFKxpIl7pi0z2PwZTDmgMn6pDDh6Q0kSDr8nbcV5lrPE3VKWh95EhMS3BT8Af2QhJ1DirBNt1Tfv9SxTqpOEfUqBzPiesbSNcTY9hyUyUjVfv9SToOkTTwi875dypThYrohe8L2t5vQQXd83JCET8qiv7YUaTNh3ysorJSsN+vJQcFkfTKgnWrGIVnDfp+oKWmaCr/K9mBCoh76lNiesF86pBKZxKQReyevIbm5x6ckVZ+0fSk7Ppf4uC2Imqp3ziSldOJUpEonpGvqGpNajN5rfiBseNtNOXc8cR6x1E3atpmw7x1YDc44HjJEPUHS/u8+OQVCQi2WrEU6VWhM+Km6ShzEgvJe6E9M2IPEdTiUSufZeGxgJFWbbfT60j62dpCbi9RkJmUySN2PqUlR4HDm39eedMMQMCqUlQz3ApCO0WY79n2AVhp6BeGuOfaYYLK+dYyJOpXgJEvSXvk5RB2qtsuIuhS+vbcEOeI4NGn7yJFfbMv2t/vfS4j60Ch2PFtK9JHtOlfPyBFNyMCGDSBYHGRkfgCAXqPh6gTKk6qcx4DMYKTAZH2rCFMsjm3TPlmObcQpkvbLBOUpByIiwQlVfm64UhALPBHulLNxp+yqc5AjkRwpp/anVOVOsj4kUecmLs5uDT2UPUS7ofakLAENZSIJpG0AICZB/b5Is+FIO59Uhb3F7wWOnBv8tsBkfWsQ0S96XelSKdov428fS8dR8osJUs8lNImRdQabEZ9MOZstsdVugZR3dSxJu785oi69BqRdP4phntv/kbc2OXFL3NcJm3cpYfv9GGU+K+yDhjISuCXqIHUpZ0G7n2CbNeNwCFXfueQmc7y7qeQmwd9IzZ36PoeggTFJpwb3ErIakYolgpRkfUzynvJMnyJq/2/qOiw9nymVNxX7Xao5IZOdEPeJbDdSk/vHlj5nsVe+WW+8DQjblENA2KY8h3gx7gaYrI+ORCpFj5zdvp6oCdtz6UBXQtRUnXOcroJtRCrMQCKdKWHnPI8pe/uxQXl9UyRdgtHkK0GMc/pWZLNe6VOQr5uWukcOhUg/c0G7TqqGDCTsILEK6XgGsC37joIla8b2SDiTRTZpiqhLvbGpfak0ksHfQo/xHKhlJt3vmLCnMJbCQsK+TYw1CXR41txY8ynMiV8nHboSZYI2Eqr+lHQdl5s8B6JPKRNNeKDfcfShXzFhG1hS9hzPzGEck30nwWTN2BZTixMQqm+PqFNLVOZbzEubOaJOrdFMrixVsMDFFuFF7nwC1WriHOfYUedgiqgpB7sSFEmV8b6EJoRKvTplr6a2l5JxrlxKqh5ti3LPExW5gsOz5LKg9eFgkUrcljcbclI2Ezbj9MFkfRRMxVHLJFFLUZttG6h8iwg+QdL+foqwAdrpym+7xEmsJI73tlXeQEhCJUSdO/ey+5K28+b6WypRF7ctkJSq56rJ40iEWNL2r2vKG92P0x4k7DRpjzOfMWGfPViyZmwDeik/KuFJjqhzkke8fYvwppxXrk/Yseq7RLIu7V+JB3OqLWA76TolTfv7ppzPKKQkzLitKRIsnfzEKu4lhJ37Pmo7tXRpaTw/9epYtbfQsufXkkxoZTHZTNjnCK1XJkXh0C3GeLTx7NTetmDQ9FXf9iPFzuwqJGsfi4g7Uq+mbJe5tlKS9VxVeMqGOnXMlrHNPqaIOpcsZQqkbwFhqy/SUCRU4Wvt/XGMtOtHCXmb48dETZlj4vCurKbGpS4Fgnzj5pxd+UjKZsJmnBmYrA+GCaIOnMqiEC14A5c3uJUOtFtI2n49UxIx5VBGlSlVhU/1J4e55zo35Cnevpaoi5yrUGZGGN13R/Je3cUe4iIm+q6IREuQIurYJ8Nvi6zHV8u7cq6+xCIhcUw24w6A1eCM5SghavuVUD9T6kFaZZofdGPv3dH+SDU95cBG2S7z7Y/Jwm87h9XOaDOIJGfbDX4fwGEthdhGn1OHT12rWMqmvPOn49+rIsKeo/3ov2ccLkk1vvaebU/13acvjUmbJGwYLZZWLF2fO5isGfNRsBZ1P0iFpOwvwkFhRKYTcdaUQ1DatjgerEskMOeRO0XcjizmqGFTUmSJFEwRHdUf6tjcOcRlc1L1GlAajTg5zFzkznl8jY/jwDc1OSw5Np44mI1DuTC8yydsYLw2IhP2WeKekPU6AxbDQmApUZfcgjnq1Ph3LKnHUrrWLp+1WzWqyxJPMdkmpFDXXilm2XzFWJWa61+ppEwR9aGQ0mgA06ryKSzRDIT+C+PEOmvbykn07rn0f4+O959zUI6Zg8Mm2W/hZRAEEL6PtHMog3EbYMl6NeaRNICAqOMBJYdAbe1JwwMZ0+pmSm1d4mW9xmN4VBehio2PTYHqB7n8JsJzKvEKz6nq/f1Un6h2D4VB/Uzc4wUkTF3TlPbG/aUW2VgSTx4QaaZ/1Pe4ruQzGjiTxRI2RtJ1GIMdtOJVxjhJ8EIejDyoWXcpUcek7dUaSb3+WDHXo9qR95Sn7mDDQzDIuTWbU/ZyKu/zsL9M/Tzqg1df8rxSNs4Vjk8xaVMESF67CS3EWuc+v92UZ3yqz+RxmXtQ1q90nP3ksakENgscIovNAoSZxtwXj7CLVuw6jwH9PkKrdb6C5+JnyGrwRaCk6WVEPWWr9hHbSlPpPekep6Vn93fsUJWvX4hqUkJaijmOSiXOSlT99DnT6vElTlRUX8ntRN72+DjKzDHWNmxja15Kxtu0TWssSBV45nxJTUzuWe3fx3DZWr8GBuM2wZL1bKxTe7vfuaUuKfhOYsFfTx1uvHXpPOCUdO3/9dXrgyRf9fWXh/yMnaL6fSvJJOdsNtf7O/6dMwWkfs/JcR5si535Imk3R9pxP3LnPceEMWeBlWQd/SRyHtnnNBaTUQMZs4D/LIbHOLNQpA4XckLCZpwk7omDGZP1LJSGZKVJ2pQPpelJ251A702d7JZ2AzdN2FMI6g5U7yFh9+UJ+6WvbszZWKdQakOnws18Gz11bkmb+oSqOXfsFOZ4wJebOBJJSSaytm2ZqnXWimIz2s3Fq89RkQ+/ZWC7jvWeTi3OGc7OFEzWjBBT+b2H7SUknVp72kco8Q42VXLAyowjw0A+dgiK66Uk7L5+ot+u7zHZTYWLUXVsUWYtsjZ9f9uRJK7Uc+H2zSHsEuk95cQ4SLq0Xbz/nWl36v6tyf42CyKzvGaSsO2+fhw4jwGecXfAZD0baaIuVXfPJZ1gQA7CiAb1qdYKUtrbuWAc8WOhexASezKblJAIklUknOFSsb05jOoIVP/LUmiWOrSl7Pi3idSEKEfYSyXq3KIk/u+kKj5qN3R2G08G4vq38BUI2nbP6IiwfZU4MHY6G9UGJuzTwH1xMGOyLsLYuYTO7z0tSeccslIDY8r21ndNh+SYsj/GA6K/vGM/EYg8z2PCNvsyjj0TiVjmqHmPCUozMCpzZKLeQiWeqzcfLjc2p/hS9VZLgC5Vec/1IA8SpxAZ0MzxREpSIJHljHEy0CvV4By6dRdB2ajHRO1L06NY6IRk3dtaZyIO73JtBYPR1PGgJwUxYQ/bEqspEaSRG1TXSNdxnVuR+7FJesq5jfJ0p9LIpr772MJ2niPFWAU+V5rPEXXOKXBKFe/7MvjvRbxaF0nYozhstmEzbgdM1pMQ0a9xfu9+n79NhFmTqJzfFKYkUWo/Rdj9PkJymRpwfXW4X3+pV3ic/jFH0EtIltQaRM5v7jwpApijTj2WNO17hE9d44B0MoSdAzlhjDQyyXXLZ4QMlvRlVH+hY98Stfho8Y8SwuZVuk4bCusUHmeiLGGyziKRhnBkp46W+fOIepTuMyVtFIRxppyfcoQ9B746fIqw3fbRafge2UBA2qOyh4jPpkwJCxfgKFkKFJipIVhxzinzxhySLm4roRXJrV8eLOdJLfeZmaymJlJzCTkpVft2fj9awHmIA70N25QzNmz3nQn7dKHVyvWs2Rv87iGdQ9jtp4l6IPJqNGD1g1G0Dm+wj0CgzptJ2L69mtoXOLJRNuyo3dF1COymCaL21btkmtSyDGLB/pkSX0ldJW1vrYYPnLAK656jrckdFzvvxfv87SnCLu3XFEruY8lEjJzMEJn6tP0O+5yPPMWZsE8TLFnfd0xJ1e7Fp/Mc+0Qdk7SI6ggGpX6Zv25M6AR6UtFjQi0h7dz2lNOZ78SWI+wc+kUhQF8LSoWf6jvlSTyl7g/qmnhbtyB/YJ5Unbqu7rqnHAgDNTYigiqoPzgu46ldck3i53wKpc5lixYkSU2AohziAEZ5xIPRXEgmbMatgMmaxJioae/v+LDx6j5xbGugLqeac/BI2z++RMKLpePR0oilq0359STs4SUEFEvOsdOdJCZDikjEMqgo6bWdKSLxPd7JczwSUa9BsXRNTBiDeqKFTibbTWgXSkm1FHNU31usKR6oxH3SHq11jZ6wqUxnidrBhH1kaKy75Gdyu5isR8iouoWMtpU58YycqqhBlSBsY1+rwkEkgVjSHUnHBUg5r1GDW+p7ePw0SVMStoZCBbtghJBQnvbA2MCrJGGnzmVq+yFxCNv8qI0J8o29oecgRc6l6vUtTAQlq5/lnfLGE9bR89urukHGYseEnQ7pYsI+JthmfS9RkqVsPBBlF8nw4pSDlmK1nOdBPXaIGQh7Srruj5kg0hJpJraFx05nwHiSkIKvcXAkHXrLE9K/qO35q36bsp7e5np0fd2UvXqtVD0HyQlLgc1+6M94GdHcRK9kgjTExo/Dl6hzyGGOCrwUa+zSpf0gnd28lcqC6xgTNomBnHmlLsaxwGQNYCzWRqrv/mUOs5QFL/kK6SG2n5H2Wk9N56SAcerQsXNQKrQnOK7QeYrKclZkv7bXSvaS9UDUldgZmz6hsjWEDSjrEKchAd1AuXGwV4sX2qbnehbfkgq8RO09aG3ysczDMqlp7crU/S8tG/ct/p6rI+kNXrhUKdmPgjjvoIzv65FdB1shHdIFsA37yGAHs/uCdUSdjPmN1ob239lkKBM1w8dYTed7WucIOzgW+UFyCr792o8Hdv0ukbB9khaoLFHL/m98Hq5+JTooNL06XEChQzO2MxLndAyp+lgYJdghTAgUnOkAKFOHr3lOtvKKT2GVND3VN++dmyJs0xcQhA1wprPjQitwutH7h/VETWFQI3vetXo86A7lu9CGFuxESOIJwvb7V+qgFB/v6iAd17z+xx7iPsZStSFqac0HldihgiFr6Q2wCh20UFDozD5doRON2acbSEiojGkglshOwVFsDuJrHiMm6qRn9wzVd5FpZEaWstJnbmvNSClRkx7ihYQ9kDElYYe9YemasQXuOVlnMpGUEnVmQIpJudir1zpQ9e1MOJiVOA5RbdODcXpQp+zppr/jenP2WykkJHY9UQ9kPRxToTb2aRjC7kRjB02jDpfYAbqBRrkD3algjfQZRhrkidrtK7ILTxB10m7saVq2lqq3uq8loYXUdaK0WUNZXyUOcEjXLYLV4PcNnlSdIeoUSU+uSQ1Hwu7lRk/kqTWo4wEkDjMJ6i5Uic8FtcKV73AWnp+MyneBRsH00wtls/HCFWpU2EHq8BoqMZC1RIVGXEPqHSCADo21dRsHPOOZWyihFZTbQhI/lBe4/6yQ4W/A4EG/AJST3qFwK975CaJOhnHBvl9WutaeVE2FdKVX6mIcAqwGv/MojaX21bcS/uIcyZp99ZpPXnG88owJd0pCChzIFobmBPUV2DN9h7OU5iBru0YFaT9Osq51bbcPxyh7Xh0aSEdQAmi9wVRCGsczlEuRh8ChbLXpazhsj4nax1yp+phEfQzkojDi7yE5VwXqcN9xExjHYLN0zdgO95isHSg7dZ6olwzMsUTdI8o2RjmDxQtj9OUjLf6aWNocKFt0Sh0e2+f7/cSg6VTftTaSdaUHspYY4qs7UaHRe0AMtmytFaTY2b4R9l19u2QzJ2Rrsi7vuYul6hxRO6R8GrKOeYmwt1m+DxNmn7me5ksmY34fykxQxPsXtyvitbBVSNgAq8OPCY11quwzuSX3lKwTtuqEbcsnainqfjsQD3Bpz+MgpMpPbegGk8IUoX6/HCHlBrHSLGNB3b4jnFeP20epw2P7tbPX+9ti+JL1Tu9Q2xAuGd2fTg+hXUp0NqzLtK3RQesOSsxThW+JQ3tAj9qLJow+gQdrTgcSc37iUkLUblvWma3Af2GOo1hqYY4c4nh1su7Eux48x364ZO9IFr4LOYczxnGgtbsHy48/B9xTsnbIq78piTo1S+9faPebII2Y5EhiFvFPWo3ne6+OsjFtJF3n7OCxOpyKvy6B1BJOxq5Q2e8iuCeOvDUUWrSo0JpwLm3aFaKC0OWq8ENoH+L6ye0zSL2kbEqqVsQkC0iTb/B7IpHMEmxtmlhr7pi6tqFpKa0ODwnZX1bT1sPS9VFwX2zWxxUJTgJlUjVJ1IT3rSNwn8hNKxGZR6pHrV2ij85814k1mEFLi2FbVIYm2X/WYDKG151PRuVsnMTSUq/oCXuwZteQ2IkKtWfdrnRtVebDVknYIPv7QaxyVnzeheXXStVz2nFx6iVISdVbrEyWO3ZzYk5IyFPX7dAmkJRZrDefjfong6MZjLm455J1Wqomy3ok4MNXRccSLoVc3HIcNzxHCqB+b4FYEs0lQAlte7Dn5alp0QHYkW1Ib4pRuXqFsNoDjcaj834Cpb20pbpcat7EGY9SpW5oq56LWNU9JVXH+3LZwo6t6vdBqcOP0m6BU2e4DYiNp+wZfgRw6NZdxPSyl+RRhAd4SMqDp3K/jbBp+0gR9hwEE4QoHekxvaL9c6AkIa27Ph5aQ5pYaWcKICCEmUJJlyoKEkprVDbIq4FEhZ0N37JJV7SEsqFxqQnM2OP3sOrwtei1Aolnz8En6T53urdvIOOxdofCbFs1kdOc6ucUcs+tn8s7VT7uo2+eCbZPOL7l+hfbrP1wLkfYfTgXmY6U1eFbg9Xgdw7lqqdYBR5ur4IXPbBjo+r3xypZCiP1d6QW98vk+jqcYXrd7DmYIzHGA/9IrY/hr7OlupAshc7EUotQRS4hIMUgZQu3DQJCG+9xaa3c0l5z10dKFe6fU3xuS80Eh5CqqeQiS++hue7TRF1ixliKlGknh6nndspxLNfWVnbuUUQIMZk3GNThgwaP1eGMZbgnkvVE/m+7bUoFnvo9mqn7OcG9/cWDCBHONSWtxGWo7yWISY2SRKfU+6PjAmnHTEoUupFUq6CDodgJ1jWkkcm1U5YPqnDAZEQTogJ02/cTnrNPrEY9tIPS3HJLoaACJzMFf4KUtlNPSdYlSEmxc47JIX6/lvYvre1ZLl3HHuOUdzhnNzse7os3+D2QrGmiHnbnZvFR1rJUucjBzJe0g98j0golylia9j1zcwMs5Wy2RDojpcUCSTTlUTz8HqQ8RyI9gfSuZ/Qb4zzBB3u29PZ5moTeWz+MQy65DltkGts6W1npPVPevxiOuN13/29QLiNV5/oxdd+3wvgZnF5NC4id7MonKCVlw6x8vgQdRpP04032frKEvQpKrP+cAe4BWfvwX7CE/bovsE49miLsof2QuPvvsQrcqsbdb/9vsu0D5WouQTxAhueiPOnPSNfKrlGtPNL2iTtF4jlypFXUeXV4cT0bStWLVdwZDU2s/vbLU9tKUdpX/3lN9TGuN57srmk/1SeqD6PJxgxVffzdJ+wxqHHnng29jNW4w0+MABKEPF5Va2yTWtVy4aCTC69K2RSXEvaawa7k2JyUNfS5G0hbKyjdokOLzi5/2aEz9K0VOq3RaQ2tQZJ40L9CSSt3Xrfpxb0EI38HHWppAIx++8cOZaal4TnStV/vEkk7R95jx7ox2cY266kJSsqmPqXBAoZnhrJfh9I1wOFch4NzMFvzOQec5ki0Gmkbdb9Qh/cSCfJ7HiWz8Bxxpj1ry9TiwXZCWpnbbsm5zL02gbOZp/Y23uGDZN2JFo1oYKhbBYTd9sQ9rMHl6iH72d9rSvpZZiIokfyOhVyWPCCcEJnf29mpk33K1LlWNb70fUxNHqeuQ6lqP5Co7W/aa59yRktp9Ziwl0BrsfpzDriDDmYTzmT+TBcIZ76F8B3G4oE/duzq/xaEd+mIiFwsc1C353w2JcG7duOUpGscrHKOP/E5xWXdeXVoIXQFJRp0aNCiQYUajXXYGQ4avnYYyDo1gEpIdB6ZUc5u/r7RvTrxUC63OhudHc/zb/CIOidZbmljnnKAdJiaKJU4lpWuVZ56T/3tqeOodlz/Us9N73CWC+fq3+1U/DU7nDFo3L6ocFBQRD04fcREPaXuo6TbGFMDkC/dTamofQe0lN0x6wRzIBtgro5R/yLp2rerDqrwBp3w/xoJu7F/3adDh044G7eRyoF5qtxgiU7Str3NK7G1BOubEsLtXUDUfvmkOnihinpOP5cgpfZeooan9lOmmbn1xf1JStiepD1yOBPRuDR65s5D0jsV3Bc1+B2TrDPLXgqKpDHaVroMpkM8aw96k5Bm52ZkijOGOQnbX086S8yEdH0bGPpv4quVbtCJBhI1WjRmoNMVBNrAC3ywWdt/YvB+XnI+t5URaw4CLUUULjcmbJ+Q8vbX0olNVrXt7UtN2OZOBueWJ+3Po+xhRN+8UC5K4vbLlfQ31J4Nmhl6/Wu/Aol0OBdjDrReR7jnErp1x8jaIXpJSaIeE/dUpiiHkphSChRhh9vSatgUYfd1LRgg5/Q76Es0WE+pw51mTwgZnIdJhmLirTs0JizLU/NfAFD2mihodKLty1M506nwpZJzi9WiW6nDp+5JyT2bMi0EZRP+DsP+9ROU3L0+NPKTkGkJe0oTtAQjNXu01nW8nGagDs8SNqvDS7HW7sw266NjfMF7ZzL3myDtgKCJbaTHp8XU4Lk1KMI2BJhOq2iO21aKXCqh+7Z5syZ1h043kFaqduhJWw8rS3WiNQ5ookWLxgv7Mur02N5/bKSSxcwhs/g+UtL1lETpq7+3Jupkv1dK0iUTlv47lb98YoI73r5wsk2sCQCMpWxH2MgmTAGYsBlzcEfIOqH+9qTo2DOzf4kJyZok8gzmDFYpdXgpASYH7ZnrYZdiroq/72c8cfByn/dEbR3NWtz0202F6B3JKuz67621a/dKcd3CZUTbGkuk62Bg9rdnno8S6RsYtBJT5bYi6mOZTJZ4ewcTkoJ7RE9+5hE2vbJduEBNrJ3pCTu5nGZfEavE10AJ6DWJTc4kKcodIWuHaaKOSZoK8cltPwTmqmBjb3FKyh7KTjvVbAHK8zrnjeuk606bxTiUL5HAkrUw61eb83CR2NYRzUrXgwNbFxAVZbstOYctrhFF2j7p5o8N1waP+zTljR9vP0Xb/BKCBsbnklpVbK4mIweSpIl7QEV79GU9D/FBRR6tf52UsFm6nsJ9STd6B8haRL/ofN++l6YpN22zHu1fQNpLCKBUquvLxMtrZpY63BLUuWVDYgjpWkNCoUGrAOlJ1k56ruwj6gK3HFE7QnYq8MHhbDpbVSmBrbVdU2pYSpKmtBCxE1RsV4+PJ79vEOM85/6W1Jf7TSFnp85tW4spkqa2GSezKtByUSFdA6HTKnEGg8IdIGsCnm2a3D2h7p5DyinpMdu9jVWMGmpE2HEfbxs56VpBWZ/vDgpNJGGHC1I4ojbq8DAE6RRjpFNS9lLCzrZ1SxL1oTVP5lmhNUQp8t66TyX1kVqQBGEPz6pEOgabpesSsIPZ2SE16w2l6mAZxYL46rCuZTiW/W/JoE4dPwdT0leJ7VpDokPjVYCesOMc6j5RK+0yuIWq8FgFfigSS11H0hlxFH41VouXEHZJX6Y0K3PuMWXemDo+1FbRauY5TmXUttzE7Fje6dS1AdBL177DWdBfISOVuIOzc/vq8KFWJmwaeqXNepW9+4g4zlN9MBAq8MCRLPT0niJqIWTyU4IcSW49gPhONsAwGOiewNKf2wKl0tTaxFtrbezXSrfmL5r+t9ItWn2DTjc9UZt9Ta8Cp7N6xdcoT2JjG2l6UpO7jqlrTYWbjfpI9DlefS13T6nsXiWqZMCQzFJynSrXv3cztVa5vh/jWY7ff2pM8H8P40u0trwnNIwFh1wO8TMfou8gfvM3fxPf+73fi6/6qq+CEAIf/ehHs+V//dd/HUKI0ef555+f1e4dkawzKm/fqcxti9XgR5yJbz3A+BJXyhv50Fhru4YGFBpI7AIJW0FBevdGOWkdgwTdWaKPpeop5yRg2WA/9xhKIvUJ292rWNIe/XZRAxmHqNI0nBSCzG5Omo/qS0mS/j6qzrnvF3XvSqXqOchdyyn7dOqYwCzm26+BXltk9snI6czZrweHM6EVS9cFOLaD2csvv4w3vOEN+IEf+AH87b/9t4uP+/znP4+HDx/2v1/zmtfMaveMyTqtusipv4O0gBl13VJ16SHsZXPbnDuYTdlTi+ooJOx+MLOErbXqx6CYsI0de6w67tOVOmKO/o7LT5NYCbFlpenENadU38CYuOeS9pq+jhy9UtJ0IqpgkrQSKV2pyfEcJ76Se0Q5jM6ZlE9pBkqOjwnbhXUB6AnbfacJGwhDtyQ4lCuNY9us3/a2t+Ftb3vb7HZe85rX4Cu+4itmH+dwZ3Qsvgp8smzBYJVVCS5ckvGQWKPqLonfXYMpVaZPtKFKfFCDu5hqsz8m6i4gav8axI5JpedDqatT5bLXL5FtbUpFPqUez/b9AOph9z7k3ol4P0XU4TFjB8/4Hqak6rXnWOrtnds+2YanEne/KU1frr0wusVsYRwGL730UvC5ubnZtP43vvGNeO1rX4u/8Tf+Bv7rf/2vs4+/A2Rd9oLFL0n8Iq3BkoGDan8L27hfZg05HQPBIh8EYatezW0WqzC2bIVOtyFRJ0gttvUm+7FEHT5B0iXHpOzaub5t7XcQP38lfhs+MVMETtl04+9zTTWHtk9vSdT0+UbXOen4yrbruVBKrP4AwDPPPIMnn3yy/7zvfe/bpH+vfe1r8cEPfhC//Mu/jF/+5V/GM888g7e85S34vd/7vVn1nLEaPAdvQEiQ9CGRUyUfyzM87g/Vj1z5QzjE5Tyf3TZfJW6IIFz2EvClq7GdmiLqVH+SfS2QgukDM/eVMFMkVd8T++M+5WzJo24Qpp8S++zcZ4KaLA/tHuYdLFWBU5OLQ/TFv2++Dbu/ln5a0oz9mj3D89jKZv2lL30psClfXl6u7JnB13/91+Prv/7r+99vfvOb8T/+x//AT/3UT+Hf/Jt/U1zPmZJ1uLScryqa7cyScNzJ2QapzEfHIuhD1OOThE8QPgkccpITD2rBWCTispFau5Coi7QRCY/yqTJFiSz8MqX2amI/dR9Kn4kcUU8R9txneGxq2vb5iftTOhk4pgmLImz/uyPstIZGwoVy9b/Zdj3CVjbrhw8fBmR9SHzHd3wHfuu3fmvWMWdK1gQK7dUpbLUi0RzCnjsALlVpTw2UATEvIOxS6ZWUplPpNZ3EEfUhPscpop66ZnPs9aOyRN3Z+F93HyLiLiHtnJRdghLJcorwKG/2ZHszEg0ZqVNirj06Zwf39+fanrqOayeqMUn7k9B4eyxd+8Q8lq4Z54zPfvazeO1rXzvrmDMk61x6UX/7YdXdqbzCgSfoCuJes5BEyfY5i034ZbZEf60iwvY6mWw35UyWkqZnX8uZRD0rPSwI4k6QNkXYfv/WEolvXy41F1HPTlxv+JuWZuPUsFuCJvJlUvVcMxJVNqUSNw1gpA53Wc90vy+WrlkV7nBsb/Avf/nL+MM//MP+9xe+8AV89rOfxatf/Wp8zdd8DX78x38cf/RHf4Sf//mfBwB84AMfwOtf/3r85b/8l3F9fY2f/dmfxSc/+Un85//8n2e1e4Zk7eC9LCn7WOzIkUBKnTaFktjXuP5DYIlzWGr5wFG5hZJFTtqmbNfU8pAlbVBEvYakqb7PJepU34Pzjq+/VoF6vISw4/0ppNTfsSQ9SKKR5I2hfbr+TLwy8Yz1kxJRQeuu+H67+krC5IrqWjHRWSNxj6RtICRs0O8AS9c0lBZQK8h67rGf/vSn8d3f/d397/e+970AgHe961348Ic/jD/+4z/GF7/4xX7/fr/H3//7fx9/9Ed/hMceewzf8i3fgv/yX/5LUEcJhNZlpnkhToXXCXu1VYEL728uvrqolblqRsKre0r1F/wlyGZO+NBSlKw6Fu+bwmxHpxm2VEpy3oKo5zicIdPG3HMHomuduAelKurcdoqoA89kS7pz35W5iO8RlSa2N294EQP+MTGS16rXGpSHas3B0vcifnYD/wsoaBvxAPvdFHZJUlw9p03cfb8PhJdeeglPPvkk/vtz78YTu4vF9fxZs8c3/6f/By+++OLRbNZLcCoMXIhUelGqaPmAtgVyKyWR3fP2JyWxQlX42tWhSrDEZr2k7hIJOxy4l09yivpWSNRLfA9IiXmhhE2hlKil3AUELfvvIXkH57DwedO6A4SEWyWtlyi97ym1c3/dJjIWUt8PhTU+BGPJOTSF+It90ElSWBUO4N7kBj8zsnaIXwzZ/w2kaMijvLAOc73EY5QQePLYDQi7ROW6FVLqcNN2wobtH5/ISraUQMf1E8cWE/VUuzQBA5aINiDsHFHH0rT7LhH+9p9lmXiG1cQ1Hq9BbZ7TCkYFrrxFW0jbt5+ic+K9mEvUo0nTjOclp/FJlY3b8c/Pbe+dzShVuKBSkDJ4PeuTQ96xLHh5Ihv1MXNlx5IhkCbsUFqoEKd0nHJCG9V3BAm7FCXSH0XYAHobNjBWX1JLJcZEPWlzthhNRgrjp9cRtStDDNzwrluCsFNIPmMTRC1F3ZO0tKYugQrSkrWAhMz5ZHivoRoRswKws987W8bLSiYAqdETNkVeo+YSxJoj6kP4lMzxc4nvTdFkXNhnAOMlNKOCYOn6fuCMyJqAtVEHmw6cAKXEMzRF2Kk+Ud6i/b6IgLck7DV2aYepVJuptvoyBMnEUnbquOAvQdRT16HoOq0g6ryXtK/iTJCtR9h+nZR0nbPT5iRqCQkpdhBCohK7gKB7NTmZjpd+zuL1x4frsINCByFMvncFk6nOJ+yUNsVPJkJfSyS3HyOuOn+fhzIpwqY8w7VWE+9yrBq/v1BY6WAGVoNvCP9ipmzRtMr70M4yU5hSg1PSdawOpwjb1U3WOeG9uxVmOXBFYUnBPmIik9sW/15C1GWdnlPHtPqdnrCFUrbRKoxt2KV26hg5oq5EDSl2kGJnidoQs5G2bRkrGcuCthXC6++IWwvzXcA8zx0AKQBgB+3y09nnPJau4+tWqgYHlhP1VM6FEkl9rukrvx+DdC3Aq3FFOHbo1m3hTMh6HtbaqkvJ1YGUDhNx2FP15gjb1EtLF7REQpN2iTfxpNp1qYPRBGmX2gJTxE32La6j5NmgHNuShDzPTj4+z4Gw/bbmaCMcfKl6dE8j1bdP1BV2PUlX2KGC/a7rweGsYHlOR9hKqJ6w3XrjCh060UBAotMNIBooDQgoVMKstta3EdmxR6tXIS9pH4qoqTK5fAtT21Io05CxdA0Ysl0jWTNZHwRu0IjDtTwEjjVpFSGwzFYV11GCUuk6R9h9Pd4AXkrcpY5IU2VLUKISBNKkPeVdO2WPDn6n7u8SG2UBUc+qL/FM5NThxV7gHlE7qVqKmiTqWlz2ZO1IusIOtTZ/K105ZTmkFXdTtn5HyQM5a3Qwi68oodChgURllkAVgNQVOmGWQ1W6gYSEhgxJOnj2x34dqfOfizUZDEvzLQD5SSYQCQNC9glSpon5fkvX9wFnQNYls5504pOs48mEOnktqHCuVNsxMffvXewpCu93YtBMkUDcdow1RF0iAZP9miDtbJuR+ruIqGciL93kB95xuZykNU+69hEQM7GKlk/UvbrbEnUtLrHDFSrU2OlLQ9SosNMX2KGGqdEQtYSAEMKeibC91vZcNNxioBoarVbo0EHhEi1adOjQoUYrWjS2tk40gcrbnK8CNCDlDko1I3V4/z4R5LqlNF3im0LVE6wBPkOSdvX377Wgx4axKpyla1aDnxw8qXoCFHHPVW1vhSl1eCA5RN9jKbuvM6Vio1TLGQIH5jmZpRzfSq9dTmoucUajyh6bqEvV8vQAOiZtirBTTmSlGNupB6Lu/8Kov3e4wg6XlpwvsNMXqFDhEjV2okIFgVoYj2QJgcqRtX0NlTaErb2/hqg1Om2U4C1qdLrDHhWEk6rdOfVpNjtoOxx1MCRO2a/7a7aR41gpUVPb6ed4nvkrBwFZmC/cEfb9lK4V1k1XzmWqc+JknSDmQAXuqcIjkl4yq70tUB7hIykbSNrtUrbe0oG+5FptdY0mVd0pabtEkp1Qk8dI2oUniTpnp566RmMpmmoruD4TqvBYaxM4lPVKbCcnG6KucYkdLnGhH2Cnd7jADmZrjVpIXMoKO2kIWgigEgJSjN/KTgMaGkqb752S6LRGoxUarVBphc62L7VH0n2HAal3UFaylpC9d7i7Fu68/InsFEr9H+aW8ctOStuZMkveJV+63mpSyjgPnDhZO0xL1XPsshQCB6UJMjkUclJ2/91Xs0XE3W+fqX6L26eQ88hdSuKlpD1ZT6Lc3MU1suWKz2/ZddjiWYudynrP796xzDqTWRt1L1HrHS5x2UvTl6LChZS4qAQupEAtBSoBVAJGwhbDm6gBdMpMIToFdFqj0wKN0tgpgRsl0ChhyNqqG/tEKmJwSNPCKdGNZCo1evu10KpXC1N+HME1mHEdU4l1tgZlf8/BaRRibROdZpXXumY1+K0jJ1UfB4cg7ZKBeYq0h8pSx1ekem8KKfWdPymIPdS3CJEqsaenyoc7DjvYDpjn/U0jJV1PS90OqQmqr/4eYqmN+rsWzi59iZ2+xIW+wgNc4RI1LkSFS1nhqpK4rASuaoGdFNgJoJLoCduXrjWAVgl0GmgV0GqBVmk0SmCvNKpOoBESe6WMnVsPB2p0gb26Q4NK7Pqr4OzXvjo8lqxL382p+3QMrVqpNiA4N7sSFzybtR/GxWtdW1PMqoU8NuzMAXGiZD2Oqw7zgE8v2uE73Uyh5EXeWsqeQ9q5/am++ykqpzB5bvHE3W32CPsYKnKq/n6yMFP9Pbc/ZVje5hptiHMso9TfLpba9/i+wAPs9CWu9GN4gCs8Ji5wJStcSonLSuKqEnhQCzxWAxfSfGoJ1EKjdtK1PVtjozYk3SjYj8BeAdcdcCMFbjqNSglUnYBQGKRGDQg0vVRdiR2g0avDU+Fc8XW7LaQ1QvPt1rGne7hPBktnhjupFKT3T7q+DzhRsg4hIvIeEXVxPWNP0jnOUYsG01xsauwhnrGfJh3EJuxhc/tMOq75RKwJSXuCsKdWScr1P1tmhfq7BHQftpCqy9tfcv8o9beE9f62EvUlrnCJSzzADo/JGg9qn6QFHq+BBxXwWKWxkxo7obGTQCU0KqEhBayN2qi5G0vWe6v6vu6ACynwyBL9rgMqISE6QHdAv9ifJ2Er+7cPBYvCuVx5c6JjzdOsa7tgDDgmYm2Cvz2QroNn/X5K16wGvzUkspVZgqYPCaXpfnPSAzo9i82pZOcMnrlECZPq7VFd6+zPq9b/pZ7jmYRNqcvneH9TSHl/z/HgjjGtHVjjVJaqr0DzU+AR7kvVAEbe34GdGpe41IaoH1SGqB+vZS9NP14DT9Qaj9caj1UKF1LhUmrspPLIWkNpgU4L7JXEjZK46WRP1o+kUaFL65hWWVu3uXaAUronXg0FpQdJusJu9Nx1uh1U5i7ygYjHdvsPjUP4s8T5Fcw22QvKLF3TMGrwdcefA06QrEPQC3aE6u+wfCrt6PKEBVNhU0vbWqsGL8X0QJ+wVeshReSItJ1EBM9Ld6FKfE54Umku8uEc5nn3pnFMiSVN4lSoXfy8+Utbulzfg0PZJXbaeH0bG3Wo9nZE/XCn8ETd4bFK4arqcFkp7KRCLRUqaZ3CtECnBG66CtddhX0lcd1VeNRJ1ELY59c9OAJaa7RSoKskOl2hj83WGq3oUKNFZ4la6Q5SGKnanIMK1OBrHB7XStVFprUNQrhS55GXrhl3FSdM1oVSNZZJZnMQE+u0F/N2sZYppAh2rRTuQzu7t1DQuhuTNuH0kyLsKWe0XBpSf/94R6bOE1RvbglKkxSEagkJKeRgr9a19fy+wCVqXMoKl9VYon6iVnhy1+GJusXjdYurusNF1eGiblFXCpV0iVCAtquwbytctzWu2woXUmEnK9TShV3J3vlHawEFE97VSIFaSyhdoYXCTu96BzPKfu3CuQK7tbeUpKk/YRbZ8DmY622+xTgQS9d5R8r7pwpnNfitIH/RRIa0cypw6oUhVU6ZGXlO3UY6QyVeVEoF3v8tUnd6GaoKEpqkFmBIkX1YyEjXJnGkHEgbyn5XUMZtNzy/DGED86XjdOG0+ntbop6qu2AisUJLQkvTmedd2JW0XEy1TSNqaHqHHSpcWK9vY58WAVE/3HX4il2DJy/2eNXlHlcXLXa7FnWtUNUaolZGqusEukZiv69wta9xvd/huq1x2daoGg2gxrCikct4JozXuDLSoVaA0tp6g18iXgzEkPcQzpUi7NE1S/pPLEsrulTLtZSwg3GBMCvBStW8wIdddWuCO6aOPwecGFk7eB7gBGIP8H77gSVsYJ5qvORFzRH2yAGNIGq33yflFKGnJjAy6ruLgzWOPzWEbqExJm2FFlLW0JpQlWeknk1WBTuaRL1RXVpNEvYihzI/D7jLVmal6sBWbdXfl9jhUtTG+7saiPrhTuPJWuEJS9T/14NrPPngGo89tsfuQYf6gYLcAXIHiNq4g+tWQ910uHzU4uqVFlfXDa5vdri4uUAlzBCoNDwpX/ROaa0y8rKChlIVOl2jtYlRnHf4IE26RUKG5TQprc1WEQnB9T2AbXoJAul68pm8X9K11s4csPz4c8AJkXVZtjKyiKAJabrFedK1j1Sig5FTGkHYSXIvIGzT7zFR+3ZKf59AFZBxv8Yx6HWK7U645Q2VDZ3RogtI21lVFVr/ZAO1uO9AlJQQsFyaHpotPT4ul3tWSlWqRxoUC/wbhrCtKrBVO6l6Z6XqK1nhynp9P6iBJyxRP7nr8HDX4NVXN/hzjz/Cq77iGpdPKlSPC8gHFcSFhLiogNoGW7cK+rqDerlD/Wct6j/rsHtZWVW5m/ABbpjRWhqi1kBTC7uvsoOtRmsdzfz4ayW6PruZFANh+x7ipf4Sc5KgzB1Hcm2VStej1bxS0rWQQQpSIfLnwrgbOCGydpgnVftSBVl+pc1oKuyIsmfPJWxSNT6hEk8Rte9wJMXgYGTKhlJYkqxhyFpDQeoKSnQ9aZu/biAeyLYn7Qk7tjvH+HoWEW4BUU+FW5Vt3wBH9EiOJ6tmoQ5zf3vC1hUurPp7J01msqtK9OFZj9cdnti1ePJij4cPrvGqhze4+r8U6q+oIB/uIC5r4KoGdjVEbe9jqyBuGoirG4hdA7FTkBcNZGVouO0kGmUIWmmBVms02oR3NcrLeCYlWi1xoSt0qNFhB6UVlOjQoYJEZTzF7XKafSSCF5VwW3bp1PjiR5xsmSscHjGH742MVOFOur77qnC1conMNcceEydI1ikc1mY054WPnaEOTdgD8XVZx7IUUTuS7r97RJ1z6FLooN3ShnY5QynsWsRAkLBCooYWCkq1oU3RnY83sBxCXenaGWO7dg4tvQzPzLQneEzUsQrcaFTMMpdSS9SoDe0JYdKIVgIXFXBVaTxWmRCtB1WLxy4aPLhqcPFEh/rJCvLVlxBPXEJc7oCrHbDbmcBpAKJVwKMbQApIKQDsAShAdVBdg6atsO8qNEqazGZaYK8ErqVJtNJURiXeSWCnpXFWQ4VKV1Bihw4tKuzMGthaBt7hbjIkRAWtu/4abBVjX4KpcWUOSScTohDStUuS4mc1u09q7xh6pc1as816DkpU4LRjjTl6LGGs6k2pKjxSV/t9OBRhU8fFiInarGFcRWQ9fAcwLLAAQAknORvC7tyACUvU0e3qr5QGpKx7wk7FYrvrNMuOvUj9fWiiPhZ5h9tKHAv7e28la7NilsROGMI20jVwKYELacKzHtu1eHDR4PJBa1TfT+wMUT/xAHhwCVzUlqztM9x2QG0Xy1QastWAbqC7DvVNh8vrFpd1h6uuQ6Mk9pXCjRK4Uka63iuBVgKNBHZKYi8G6VpDGTu2MHHX5lls4MK4JCS6A1z/0jHk0NEeQVup8UhIb9JyuwsRMQ6PEyFrh7wKvEeCtFNYI12XhByl1OJLCdvvS0zYtlN9mwoKUofbHJxE7atDHUlXujZSmB3QpX/N7XKHyqoc/TWIhTDhNkJUZinD/hBpJeumJ2wyRjuyY/vnXHK9/eseX8MQtyxRewNp/3tjjNau9uzV/uTM3nFUEKjkkO97J4ELqXHlYqnrFpeXLXaPKVRPSIgnLgxRP/GYIevdDvrCkrXSEG0Ll8pMtJ35qzSqVmN3rXB53eJB06BV0ny0wE0ncCUrXEvgyknW2uQX3ymJRkvsUPdJUkzsdWPOT1dAP9lLx14fGsck6r5N9/z0OcG9sadXjceq8PsBdjC7bUSOZaUShcNIEl1hP5obI5wiX78vFGH75VJe4j5hQ0f1ajVSkwtLso6oK+xQ67q3Y8r+nwjPwVJ1ixZG1qnRiD1aNL2U3Ws0IKHQoNMtpNz1hN33NaEWd+1QUvYcdeYhB+lVdR/Bk9h/L3x7dTBJ0xVqSOykxE4aLfZFBStZa1xJhau6w4OLBheXHarHNeSrdsDjl8DjV8DjD6AfXAGXl4aorWSt9w0gJYRLIaU0pDIretR7jd11h6t9g31b4UFniPiqUrhSElfK5BBvJNBJk0O8lgI7Xdl1sKt+gujOw52bxjZS5BJp9FBE3Tu7Tjir9Y6bo7hrygP8ftit2WZ9NEQZymZK1TkJO0WU2epTHpgFyHlxU5OHVBy234+wc/3Bw28byiItcZvzG9qRVsqKiXqnd6itvOVWOo6vvYZGixqNNpQtUKFCjcadh5fH2UFBARJQKlKXR2rx/ppgftjNMdTf0305lhSXsF/HEQJBuJ5HbtrpUoz6u7Yq8J0l6gupcVl1uKpbXFy0qK86VI9LiAc1xGMXhqgffwx4cGUl6wsjTQMQlR0+lALaDqJpoDsF0XSQNwq7x7pBuu4k9qrCg0riUadxJYW1WRvSvqyAvTKrc1XaagJ0ZSeJ/RnYd2xIlLJWFR5HgtwmppzVUu8INbkdL53JOHecAFk7eC9KwgYX/B4NVuVS9hZIh2rNI2yAfkmnkqfEhA1P+jYqxDDWWqKC1EbtbQZBs1CidUkyZC1EoA5X0Ki06vdD25hsX1IepSFtoCF7CTtF2ClvcdrTdQ6OFX+9Baa9/c2PshXlpFfGTdB8zYmEsLm6h89OKuykRl11JunJpYbYSeP5fXlhbNT2o68ugaoGpNOKAEIp6KaBuDDlxb4DLhuIywbyskN1obDbdbhoOly0HS6ky3JmJgzu79AnY1uXWgRmmv68D8A/t03SJQi8y2PpOvYDiU0wdxwaYpWTGDuYrcbgGdu/TDFBgx60tkIJedB26XLCNuVp0k5J2SNbrzASbYWKVIU7xEO3hEANCSGsLEaQdbAGMYBWS0+iN5btPoGFh0634aAynGwxYW+F7Un6dAZCIfKrzwmrORHCELW5x+aWCOHWpjaraAnhrE/C+8jhr5D2t22rqgHZmL9OPV4bPbu4qCAuWlQ7jao2sdc7m1/creLlFvjw18nuJxX2CQWMdmg4n9MnVof4vd+qznxdlDr8boMX8jgK8irwlJdrHFc62cqBZ85zCTtZj2+b8gYoivAROZRpPcSfanQwMdC7Pl4acIRu4NTePlFXwknYrk5DobZzfkehhMLOuaGJridlp943asodFBqSsO1JZgk7mCAlpIVw4DrfQWoqbKu4nijpjT8xE0L0JG1ayoxSvR1aDX+7biBqwGyL0ZO8uWVCApXUqKRCJYZVu4QlaAH7V8A+j15VHmHPxW1rUSjN35I+hQJDZs1r9ga/8zg9ybo4dKJcDX4MzFG1l5SNbdojVbjzsvaIzy16ALsAghJm1SyFDsZlJ7zdg2p0IOrKSmCAT32eJyrQq9o1OnSihnJLGnqpIQETh50ibEPGVZKwHdbGzt5VqZrMdy9czHuZE5SCSf+p7d9OCagO0I0C9h2wb4GmMZ+6Gh63ytbfdUDXmo8TT3wxRQpAKgipIaVGJbXjcdiMpZ4QLyCFtpPIAt+VE0bq3Z5LqL65g4k4DXYwu1UQ6u7IXgeED/OxMfcFIqVjTNjaI/V4rAqPCVvDOphZCVfZZCZaKEveqs/77cORdiVMaI/rUeUc1zQAZVWSWkCig3axr35olkCfGtJfKWkJYbttFI6d/MJg6/YKNEITZXLq7xxcqIoOiFpCdQJqL6D3GrrpIG4a4GYPUT8yj5hSEFoZtbf9jX1j4q3dJ9JJCuHU7RpCaKtyt1ZG0fuqlff9lidLU2GgU2MRFa64pB7GALZZHxWEujtDxGQmp8wLtHSlnSnEhD12IEvHYPvbXF0pUDbtWMJ2BObirjVkL107Vbi2cdMqMeA5qScYRJX5Yiafpo9Ka+utW6MSO9uqqd8tcTh0tIESXt8JwvavZeydG0jbMx1ntpFGjk8OpQP+8LuafE8ctPdXAcPCGp1E21boWgG9Nzm/9aMbiFcu7EQZ0E4lXlt1eGskbxETdqugW9qQaIja72P4N4ayiXnG56GSz/E5gKXl7cA265NAmbdsadq/JaQ99VJt9dLNUY3HKvHheLPPkKOyx7j0Jp39mFWNlDZ7KEgxqCGNNGQcgRR0b42utUSHyny0RC12fTt+PKxLDSmc1O/BT/QwWlDlVqRnH8drO3vfE57gKbhFWyhVuNba3neTiKSzRN0oiUZVaDuJ7qaCulHQr7TAdQs82vd2amFt17quh2QoeyN942YPNC10Y8jaEDagOwGtB1Wj0iJ46nT016Xz0AifTw3VPxvjBTm6oxJfSrqeKw2XStlFfWLiv/M4cbLOY/bL4UumS4g7QSBbqayW2r0H8u4CxzPnaKa1ilThbig0/7f2eKWNg4/QRgKKISFsjdY5DbLP2ewCvFw8rLLZzVxq0w6el7q32AeF25U6bnPQ20716TsWmt+6l0BMxrAhxrmxGcaatkLTSLSPBKpHCvLlBri4NtM252DWdhB11ZM1PLLWr+yB6wb6poXeK6hGQ7XCqNeVQKfcqluiz1w25FMx/dN99rwhg56Trs05RUTtnaNP6P22Iz1Ha8aAVLTHNnDe4Xc3MQqrwQ+Ow16gqbCJUZhUlCWItA8fAXMTNMQZztw2R9SUKtyskdVZVbZRL0ttLN7S2TH7+o1EloO5lq39bqVAbWRwbUlaEANpf7wN8aISP6Qcz+L963E3JJPAhwDozRNaa7PClQY6pa1ELXCjBK6VxCtthYtmh92NXebyzxqIiz0qCcPs+9bYsHc3Jke4I+u2A24a6Ed74JU99KMG+pUO6pFC9whoryvsb2rsvUU9GiXQWLJuFND2Ur7rozLLZXrmlT6iQZvn6NjS9BLMCd2KNWUUNDp77qqfmNy2Df8UwGrwMwD1cJO2u4SdOCgTqcrnkgaFNTGhSxK5+ETvnLy0JWml2yFNqJCoUaHTFRqEcdRQ1ilMAEIIqzo1A6lTS/rqyT4W1krYbsrjx8D7qu2eSApyOh9HFX4bg9106OHaeOJ+1TR0aNGhherJsFECewXsFXCtBF5pJa5kjZ3UJrxKalQvKkB2gNpDXncQV3vAWyazJ+uuA65bI01fd3Z9a4XuZY3mFYmb6wo3TY2bpsZ1W+G6q7BXEjduIQ9H2ApoldHymKllh0506ERrzkObhTz659mZeQjS7tXlt0jmS0K3Us/CqU9KGMfByZD13FCNNVnJcuQ9yhSE8ctSQiLUYHssD0+X4axXg4sOCtIuaWnU0g32QR+VJz0rCEgtIKwq3KkmO2vzjGHU3O67t5KXkOh0F6jCc4ivPZkUxjqZFd2DM3TiKX1GcuWcQ5by/nVQaLVCqyT2SmPXAdedwKNO4BUpcdlVkM3gpW3u/Q1U06F+RUFetRCXDcRVBXEhBw/EVkPvO+gbBXWtoG40ukcwRP3yDo+uL/DKzQ4vNzs86mpcK4nrTuC6A27sZ99p3HQajSXrFgqtXcCjRYMOjTkDT6oGplXgPg7lZArMu2eHeB7vu3TNoVsnjJGXdeGqWimMpOgo+QDpqdzbhSkp/fhhF3EYly9dK60ghBnMOtFAokIrWggXTgVAoRq+a4HaqqYdjM3TOv9oDZd32PfIlajQWXW4g7D5nF2fSuzUywe1udmb7t4gp3uithKozTJnJFZDiLWy0nUHXHfAdWWka5PMxvgrVFJBCo2rrkF3o1BfKciLDvJSQVyIQTnQaugWUDcaag90NwLdXuLmUY1H1zu8crPDo3aHV9oar7QVXmklHilpl8g00n1jpepGK0vNLTphVnvr0EDpdiRV+6rwgLQJqfqQRE0hR95LQj6B8BwCFfiZTUYPAY111vgz0YKfAll7D/YGkmdK4p5ShU8RdrrcekkohTXag8CGjQ6dVTu7JCkmYYlRXe81cDEcaKRxhF4FJtRHQ2tt7Yl6sCP2f5cNiqUD2O17iG+B1P0kts+4985PAaIObLza+ieY9aGNKrxVGnulcaOAGytd74SAFBVqAdTtDtWNiYdWWuCibbG7UZA7hfpSQ+7U0F0FqEaga2T/d7+vcLOv8crNhSVqo/6+VtJK1uiJet+ZvjSWrH2pukPTq/NzUrU7f+r7IbB25a3S551S789KqgJezOMu4ZbIOqd2mAhfimOZewerQbrOxTNne5Uh7KWS3yFV31QIV7y/Tz0qDAl32qwNLFFBCitNi8uAsLWVomV0n/xI7daRgLBbxEDYoXpynVSTI2g6n/hc6XoF4mdgYejOoeB7ULfoTJocrVApgVpp7DuBfWWk69pmrjM+DRUEdpBCo1MSl22FXd2hrjvsrhVkpSFqe+7KEHSnBLrWxGo3bYVH+xqP2h2+3NR41FV4pTNStVGBGzW4I21H1I2Opeq2J+ycVJ2zUR9aqi7xmUkdl1XbH2TCcTc9wk1iH/YGvxXM9oQuJOzStksJ27Wz9jzmYErlH6/IFacgFSJ0NnPvr7Jx0jt9AYUKLhgLMCFbvnOZTTSKThgZyHnsur99P+K+Ja7V7InUjOQoB7ETpuqznvVrUboEbN+s7oLnUYtuZPLo0GKvO7P8pBaoOoGd1HjUupWuXBy0b14SuOkqXLU1drLDRd2hrsyiHIM/gzBpSpVE00mzFGZX4bqzau+uMpK1Enilk3i5E3jUASaMW+O6U7juFPbKqsDFHg1u0ODGSNfafCi1922ov9dkL6PKF0nYS85jZgKhc4bCuun5uVylkyTrFELVrnUyIggbQJ88xMdih7SM09mcOpd4ey6Nvfa3uRSkSje2I6HN3YR2KWhdQwmNSg/LKwLoSXsgbNUTdYfGfh8cgdz59KpLz3v3OAil680Iu6SOjQjbh0uIMgXfR8FNxpzJo8EeZnVoCalsTvjWPAeVdRbr+vhn0Yd2XVcSV53CTijU0q6aRZB1pyVaJdAqiRslsVdykKaVkaZf6QReboGXW+BRq/GoU3ilVbjuOlzrBjfYo8EejTBE3eobKDRQ2nw63VpJugvIOSbqY9mog8RMKa/+iUn9pIRtzyVZ5h6R8n3HLZO1SxVKqSEy8YY9MY8JG8CItP3jfMx5wajQrlS9JW1MtedjpPqfUvkT0rW0RuiesOET9c6oHEWNCi1q1HDrCNOqcDMBaEVrbaKNp64cQoaGgcZzjtH0QLqUTNctrTlDZb7JgJi3V089K6Spw70DkY+CS4bjS9cNGkgnOStAiBpVJyAbBeUlKjF/JRotcN0pXFZmScudVKiFRi2HHN9qRPDmuJvOSNKOqPfKSNOPOuDlRuPlVlmpusO1bvEIe9yIG9yIR2hwg9aTrH2tTaz+nkPUZMjmnMl2Qqou8ZEpda5M9ZPujwzu/4D4uT6iaegWoLWAXqMGZ2/w4yAmbCAksNSyk3E5CtQLRqUuXSNtT7VX2tfcfq1NClJH2L7DGTSghUKFHQCbREOYa+oGdn+ZQichu/hXQ81mQO3QDJK0J1X3/cAwsC51FCvPE76hdH0bksvERJKC76NgtB9GujbOWjU61GiEsIuxCAjllkR1z68ZuEx2MxMDfVNVuFSOrGW/xGXtrczWabMYSKMsYWubbMWGh+0Vehu1k6h9or5Ba9Tf4gYdWkPUnvq708YjfClRl9iGtzRZ5YSAWZ7gkVQ9JRQc2rHuVMFq8GOjYHbqI6XyzRFnTtqeq8bKkfYazLFjUTN9X+vgS9fGLm0IW8GEyWobHy20tBJ4h8qupuXye5twrBB9jvFeih6IuutDbDqvPO3FuvXgck4D1ppnhZ50Es8Cut6hsLOEaDLL2VSwllTNqqbSZIKCibMHTAhfq4G9BHZCoBIaOwlUQgc6gk4Drc35bbKTmb++E9leORu1tn8Hor7BjVF/4wYNrg05e+pvyqGshKhnEeNMH5ep8SJ33OGeUQmRSeNre4C75mTGGcwOhrwn+KIQJ18FGKuMZ5D2lJQd1+nXkVO/rc1JnupjamEPatISE7buk6Qos/iDtos/iOF6dggXhPDVbX0cr/3rE3XvEIQw1CYXFxufZxHEOEFKzjN80UB5ZPJfGqPvnMz8zHVKY+RQ2D8/GlZ7IiDdymrahOy5hTcMCQM3QliSFjbL2dBL5crDXxjEkPNN54dnATedxqNW4Vp12OtuIGrR4EY8QoubPqZ6rP4OPb9Nf2miPgQZzlGBz8WS/p7TxJSxDU5Hsg4wzxs2BknQCdKeQ9h9n6h6C2MvcznJ/frn2q8pwnaOZBRhQzfGLm2zRkLYwVF3JguZdWrqrBTmk7YfnqW1R9gEUftSdcouV+oY1J/fpvHWx7LnzbdXz332B5NQGK4XZ69z82WprU+CApSubDpS2dusWwXslSNq9IRdiWHKrWGeH41wgZBWATdKY98Be6VNljKlRhL1Xlz3HuCtvkHrJGvdhurvQon6Ngnstog7h/sQa80LeRwcOeey6MEXrmza+zLOagaUkfZSwk71w2+j6PhE8pUc5hK2LeQaHKnEjSQmeynbqdD7tKGooKLMZD5BA9bz23qCK930JD2Qui9hT9vgwmuUIedjSNdzsIXts8DLOAXnn+CH60GYTHWuLrPNlFdaQeHSJExRtc0dLtEp4EYKXCpD1JUAKimwE0Cw3jngrZaFfqEQFz/tSNplKDNe341xJsM1GnHTO5SliNr8DZ+brUk6q1WLJthb2reTE9iZk/gx7KTsHkjfrAY/Q1CEDURqa8KmnYtdnrKbTx1T8nJRhF2aNCGedKQ8XikbdmdVpka+qnrSBpr+epjlL5tRncBA0K4/sTQdqMEjNWapc1AKFIGXEva8hrYPjRmej7JBv1/JbAZJuOVSjff/DhANlB7MHAD6hDbQgMYOrVbY6Qo7YRZ5qaXAXklL1sZmXUsrWQu/rUG6bpW2krXJSrZXykt40g1EHXl9t9o6lakbo/wmiJp6Vg4tSa/NVhYjpzErRfCcQ0JHpMwq8ruJkybrVHKILIF6DzG5PyE5L8kvXiJVzyXuuQjs00RIVzgpGQjbLEdZeWpxn7SHwUAH96AiB0qflGOiDsvQRD0HwTkRnuHnmJLU3b+SZCjZhDiWoH3NifP+VxrorIStxXBPnIStoVDpGh0qdLqG1hq1lmiV7pOmVFKgUtpmOhvY2mS886Rrl5FMmZzkDTrsbbS3s0/v8WggaTRo9bW1Vx+fqLdwLFuCJHEXTlzJdyE5MfW33y0nM1aDHxNCYnAuiyQIIl6xBP6ATanHU2Fevjp5CXIx09lYywXqcIcpwva/9ypuIYPBHfCdz7reu9jZiE3dvpouVG8HtumZRJ06b3owSlxHgrBT0vX4+I3s1slnZsazVKACX6Iah8CQH14PdmtTodGSVGIHpXfOpQs7XaHWEpWQxqlMGdIWYkiU42Dyxg9LqfrSdIMWLVo0ounV3j5RkzZqj6Rjoj4WSeeiLShQZjdqf277GgfU+ypRsxr8IBDRr/SMxqTCDCWOpQSasmdvOUuOkXpxcmReWgeFHGH7dSUJ2xxoVOOW+CSG+igMawqPJezQXr1eogZowg4k6QOorLeuN1aBT0nVvgrcJ48ptbizXTvpGvZvp4FK7NDpJngd+yx26NBhZ7KgoUIFS9jWGa0Shqalpwd3q7EBQKdNLQ3MGtpOmu7Q2MQs1j4dEHXjEfQ4YgDYjqhL3/lSLRs1jkz1cY6dGogmrQdxsmScA25Jss6HSFED11piTdmzgXnEPXfWnKtnpDJfIV27PuQIO2psrAmLnNE6pCUFn6CHbbQ0HZSfKUXE16mUsI++yMcBpOqSNdFD89CYYPzsdT5hu3A9P4udS6RS2XXVOlSodG3XZxOoUA0JVXSoBjdt6f7T9Glouz59aIcWDYy62zmTOYfETrdDTDVCqXoLop6l6p4g6ni8mJKo43LU79JICGq7u4/JY++4RzhL1sdCpAL3t8dETTqPzRh8Y8KeK11PxWFT9urcAFNC2HNVWyU27KSn+NCJUZ10W7T9Ovi9kKSD7sROgZ5a3ifsvv0sYQ91zlKFu3tPqd9nnkf8m5Kk/TJOqu7Li3EZCu5e+4RtrpVxJFRwqm+zelqHBlWfetZ8j1PPVrqGJN7DYZlUhc5K6T5Ju7zxLt+3r/p2hD2HqIMJ3MR7WYK5/ioUQc/RggW/F7wPea2SS44ChB7hB5ys3iLYZn1gJEO2PPU3kCdqt30NYY/3z3c06/uywGZ0KMLuj5sg7F4tDsz2OQkJeDuSLjUhUGQ8/TysHLBmaXhShEqovAmNEqXq9rf5xE5J1QBB2ObAQBqT/uBvU836+eKldTqU2mS16yMFbH97c4jNbKeEIV2XM75DY53HutVEPec9yJL4Rl7e1CSwpBywnKTj7/1+F1tfdI3ulpPZfcARyZog54SkYUrLSaJeijk27EPbtoHDELaPSQnb1Y8hz3pZvQnb24bSdFyPv8+XLkbXp0C6HuN4kkfOeWwkXbskNQVEnQpddIQNhE6ESjcIYuy1yRPv8sUrdANZC5N+lmrHz2wX54t32cjM366IqP2+U9+XXN+tw7AolPRxSycytl0bDcIaVbY+kznL7avBKRV4AnNUTqVqskOS9NzBpYSw59Y71b+A+ET44peYHbYmaQopB55YHTjX+eaw3rPzpWpSio6IemQaEjR5xvDvtU/aoS27g7DhXEIP19RkPqv6v1R7VBpaX5p2Mfk+UcdpRH2inqOFmdq2JUkX+VlQXuQbvAd9/YR2rNxuffdU4Qq8kMeBMEGA0cC1hYMZRbxT6vCl7awB6bRSsMpXUd8mQrqCNnubd7r+EulnanGFkntavGCKp+IdVLpISNe5AWuLwSzvbBb0P6P+pog6RdKpCWfqt68ad8+GU6M6BzQNZcK8YLLY+ZMFGRGScipre+2cNE2loPW/Uzg0UZdcq6V9mls2hxzpkxNNISF63wRgelGPuwFeIvOAEBCk7S8mz9xAlJOKSPsQ4RDi7wttuHQK0py0nRsYl2Bsnx0vBjLXuaWEsOf0PXXOpSsgpa4nOcCmpCMR3Z/Ybhc43hCHk8/RGsJOmXYiidj6ZvjPeAlRC1GN3ouUw1funrrfqs8T7pN2mMkulYI2bM9zNoSTogdpmspsl1N/U32eegdj00FqX+7YUR8Kn+UlWGJT9yev/fFZqVoidDJjnCtuWQ0+ToRCecPG3ynHilK7TUBSU85mR7BXTyElbVPpSdcgN5lJlaV+r80wFccST/VD2JXCwm1ugBoGsZx0vR1hp30wgu1EuRKilmLXl5c9YVOJO7qR81ew39/mhekJO7FxgVqmr8qu3uU5lgk6Ba2KJms+MeeIOuzbdAjT1PY5z9AcrJmQzzG15fqfGgN6VbiAJ11P9e9uOJmxGvyWQQ1oW4FUixdIz7n6Sn+XzvCpegJJuGBpzimkrsOc44PfM4k6p7Ic1L0F96J3mBpP5ih1eDnmEHaun2n1dwlRS7kbCNsj6eRE06UQtaStvEnKONlIaMIQkEHOeOeIZvZZaVynycevd7TqWoKoKal6DkqIOvccxe8npSGacnSj8tRTbcyRpqc0iv5Elb5+4fN7V+OtOc76NhAP3okBO3awWIO5tuu5ZE4Rd4l0QP2OHasA2qY7B0u1B3Mloin0khxxz0sJ21Y0/J5Uh09J12GZNApsp4T6O0fUUtSB2tsRtf/b1DdecrX36BchYbv240Q2Qx+rIIwvlrYH8wu9YIxfp5/dzk9NG08YtvTzSG0vmSDnzyflk5E5JkHec6XplFAR9yUVtrW9mYdxWzgtsl6BydzRM7Em3nqy7oUEGZB0dI5rSHtOf2ipYv5yhZQ0FBN1QHBkJeO+BWpxPUzGAi/32c9HanArmXT555mXqlJEXYk6IGlnR3bH+dCQRh1tNZzC2pulXdwjJJ74vnm/iWx2g410SEXrQxEkHKegjffnMHWfUpqZrPSaMKWVwD8nanuub5RQQE24pjRK5Hvv3VenSRKBHftu26011inzz0SwPj2y9qUN97vf59lpKOl6bczhFurx0gFpDWFTx1OOaEBIpPkVm8aq+lKntTWgiDom6Zzt10mQrpzWKpCqg8FNyGgQA8qka1eu4HwIos6pv6WsSaKWYtdL0+57JXbovbYxdvICAAUJIYYMZW4JzJLBOpwI2vtK+QNgSEWbrMuTpv26p97PVLjiFKaWuE1p51K/44QsVP9LpfGScyjVKCXfe+LZ3UJoOQcYNfhyj25Wg28AUvoqJOx5hLMuDelSzFGJ+8f4+1O/+3oWaAfWkPTUsVMStSOx4Xc6jtjc58qoVxGrV60N24sl9p3NQpVhKWHnkbIxDm0S6u+IqJ192knTUhgyr8QOFXZBfTRZu/jmgbQ7mKUxJaRdDjW8duZ6JVS9MTl4K7LlQNaLiPgS13lujvwpFXKplDqFZBIgYgJSOp7EaWTj/o6ee89GPWg6vDEsetbH7d1Nu/V9wEmTtcPI8SLhDZ1SR205uzwGiZe0nyPtQ86mc4PoHKKOt/lEHZO0f5yvfu2dlyI7NQAIXY0HNusti4lrteQajs8tGnCjAbg/31j1bYlaip39SEjsrFQtUWFnax87mLkBW4rKZBDTJqCqEkDnwrCcLdrTQuWIh8poVzKZTdWZJOgFZp2cQxZtF/YnHeXvC3X+fh3kMf7EkAAl6VMkHUvRqYx8o/MQzkfj7tun74sa/HaZZwUo+1TKxnkbBBur1PxPqtwcUC97at8WoPoe7J9J1NRyj3EssRS1IS4YFXAldmaf/ZhttSW6KqrDDz9KD9xDv/KD55xzI9XfCCcevlTtq759oq7ErifqCjvscGX/XmY/FUx5KWojjaPq1efuOlGgyGdYR2v4AOZ+5z7+saWgCbbqP/G2knqSY0LB0Ef6Z6Sec62GT6q+0USoGj1nAuEzHJQnzEPxuW6d6Okc4LzB13zm4Dd/8zfxvd/7vfiqr/oqCCHw0Y9+dPKYX//1X8df/at/FZeXl/i6r/s6fPjDH559nrd3Z7UCiDjLJZhjE5oCZbdK7V+KzTIcZc5pS8Ke6u9Sog6ky4ionb3WqIJ3PXn7H0c+EpIk7KANT/08ZznLqeuYk6jNuabV31LWff8ooq6wQy0uUcN8drjEpX6AS/0AO32Bnb6MPnabLe8TdmXrNepzvy8Z8iYIGKAJPEXqcX1TyD7TGZJe87yXHJsLeRsR9ISKPNUHSpr2J6epvk6ROmNbvPzyy3jDG96An/mZnykq/4UvfAHf8z3fg+/+7u/GZz/7WbznPe/BD/3QD+E//af/NKvdW1GDa2iEq25ZdY22L7kXqz/18qbs15SzWa6uXAjXGht20rabWVxjbZtb4raI2ql9e2ILpA9734VJxKHQALoBRG0dq5qBsLP2O0qlOFYbxqYG6rymzjveF6q/nYRtCLUWlz1ZO0m60oZ0a1ReCtAh7lrZf53oIHUYk93bMUUHs2jHoA5PoUSqnPN8zjE3UNd7CU7h/ZmLYZJJpUkdfAzWOtPeJSgcNynK2972NrztbW8rLv/BD34Qr3/96/H+978fAPCX/tJfwm/91m/hp37qp/Dcc88V13MyNuuesAhnL4ix2tv3fk7lyz3GAz03HCU4lggPm0PYh/b0XEPUKZWm25ciat+xKuX97O6zQgchBtuj0GbdZGntupCAUu3wPPnOZsm84fZAymlocvCPHIWmnMo8T+9aXPZEXYtLq/I2BL3Tl+abvrBUW0Haya776+TZTnfYo4JAYzUOFRrcoBHXkHrX2/EVAIHIc37m+7KVVixnx57bzpStOt6+LqRvO/h9pFZRm9PPUs/7uwKtzWfN8QDw0ksvBdsvLy9xeXm5omcGn/rUp/DWt7412Pbcc8/hPe95z6x6boGsndPD8Ft78aBBOE4iUxIQPrC55SQPQdilJFqi7naETTmIUU4kfv2HxOaOZBFR+2UdOQ9q7zqw1/pk3Q9ewpC1e2b89JfOmUrY+xQ7mw2NOyccCnMdcyYk7ciprPf8Joj6Ao+hxg47fWHJeocddriwNC0hIIToibo/b63QokKFCrWu0YgKjd735w4BtNrmDYdCJXa9I1qKAJbc67nw68hJ2iX7S9oqfX+2Ju/YWdJvx//rl/e/xx78FKbGurvoCa4hoLA8dEvbY5955plg+0/+5E/in/yTf7KmawCA559/Hk899VSw7amnnsJLL72ER48e4cGDB0X13K5krZVnP/RU4d6yfTnC9uEIe6sXjCLkWOqN903V1/c1rneDBTaOJRGsJWq/XG8zhRwRdY0rS+DD0oyuHGA9wY3idzCbxElSRDdSh1PhLSkP23LCTjg0UfbzKFSr8mzUTqLewdmgd7jADjvURiEuKlSWpIVA/xcw0kGnJVqtILV1U9JiGPCFXb5SDPZkpY30rSGLJpYU5pgGSrCWuEvqXjJGjKTwFUJA/B6E+/JOkfPbvB/S9Rb40pe+hIcPH/a/t5Cqt8QRyZoYTeMSCYk1VoU7TKnQtpCuSzKZTanCU8QeDDyFhD3Vfm7bWmxB1CP1N6H6lr5jFHaWoi1ZW3utsqTTOok6ImytOygxSNemCDHRstJ1nrCB/JKa0xidv+cYJ4UM7NM+UV9ihx0qXMgKOyFRCYFKDkN5SNYajTIkLdXwrmkoKD1I0kb9bkh68A1x+b/pkK4SUM/1GkyR6m3YpMk+xRqaW7KVB/0qund3Q8reSg3+8OHDgKy3wtNPP40XXngh2PbCCy/g4cOHxVI1cIuS9djJzG1XJstUQqrOSbfr+zSo4VPtlkq7RSrwaHDLEXbxOZyQt3oMSv3tq8EdUcc221qbxCBOqpYQUNrI1ZQjlYIyzlpQvXRtiKiipetosKVJItIAZc4vFSbk26qdViE8Z+fRvcMlLnGJGjtR4VJUuJASOylQSwEpTOw0AEjL1kprdFqgVhpVNximJUQ/iXHaCC2UvR7KajZUVhXuUKqC3pK0t1RFl5xj327GF2Yo4wkBKXNd4WRuDUbCwKRgcrek7WM7mM3Fs88+i4997GPBto9//ON49tlnZ9VzIg5mg906BuWlfUyVr7GlD3blEtve3NWnAmJ2A0Rkxy7t721hqfrbqIHHRO2Iq8IOla56z+fQocrkqG4c+ZpGAQANjFf0lHTtD7jB9+SgnifqYUMihEw4b3cXflYZ72/rSHahr/AAV7hEjQtRYSclrqoKF1LgshKopSFqKcKproZApzT2SqASGrKDkcKVkZ6V1mjRoRMNKuwsaXfQqGEywI1V4f41SE5KM5qfrSIacgRLTZ5jrZQ/hpRo4lLbKe3LEq3d0uuSC4ej7tPwXRVK2owSfPnLX8Yf/uEf9r+/8IUv4LOf/Sxe/epX42u+5mvw4z/+4/ijP/oj/PzP/zwA4O/8nb+Dn/7pn8Y/+Af/AD/wAz+AT37yk/i3//bf4ld/9VdntXsiZO1ASy39Q0eowpcifjEHsgyla4qwg/KjWe2y9ZxHg04kZR8Tc9I9Jh0ACZMFqf52hB0R9YV+gEvrBe28nyVEoI3poFDrykiXvpVFABpX0EIF0jUACB1K18OzFarDXZ+X2Dapa+QnaJE2uYuzU/uJTpxEfSVrXEiJq0riQgpc1QJXlcDOkrUAIIX5ACaxQ6ME9gqW0CWqzjqhKVjCVujQQAkFJTpIvbPaicq66Nh7Vrwesi2emVTOnXDOueY+QW2pEqf60I8FHmH3bVMe24US9dy+k0uMujzmXvt0drq7ofaOcewlMj/96U/ju7/7u/vf733vewEA73rXu/DhD38Yf/zHf4wvfvGL/f7Xv/71+NVf/VX8vb/39/Av/+W/xFd/9VfjZ3/2Z2eFbQGnQNaBk1lB8Q1Ww0pJxqWETdUxlZksN6OPy1BSdtF5zUjNuATxIFZC1IFUjcFeO6i/aaJ+oB/0jlU1JCohLUF5al9oSO0RtW+zhkIndpC6gRYmlEuh7ScKwbXOxWIXkAclVcfn7Mr5zmUysFNf4FIPRP1YVeGiEriqJC4r4LFa4KoCaqsCr6UjbG2vh0CjgBsF64RmiFwICbSAUhp7bdrq0PaTGWPPXu9oBqyTsimnq/hZm1LRU34elHTdHzMhEcdOqz5hu+NjKXvq3KZgctpXo/6m+jqSrAMp+n5I08dON/qWt7wFOmMkp7KTveUtb8FnPvOZmS2FuH2yjkCRZryPjE8uHGRKX/hRm17bc9sosQOWkPYabLXk5xIJxhFU/7tXaA9pRCtLy460dqjxADvUQqK2jlXOAxqAjSnWxpHKqnmVWVILSih0wsUZW+kRFaBb234FIdSYpAnp2j/nKZttHGozttE7m7uNlg4Sn1xghxqXosaVHIj6QS3woAIeq4GrCriQGpUlbAndS9atBloF1MrpHjwHMy3R6Qp7XaHRO3S4QGvV4R2afsJi+plfq3oKx1CLz+2HT9hTGE1ICcLu6ySk7DV9HWvX4vEvlKrLVOCMu4JbIuu0k04pDrXedMrORb1M1LEl23LtTm07daSkavc7lqoFPE9ofQkXrHSJGhfSOFbVQgxk3at9BVqrv1KqNmStFZRQqNGig5HYhWiM7RoDgU46CxG/qfMpuh6ja1DZc7eTFdSotZmkmPAsGUjUV/bzoAKuKo0rqVELR9jD7F5D4EYJVJ3pvQagtYAC0GngRgnsVGXNCkZnYRLI2NzbMxzNTHtprc9WDmZzzRAl9uvcsVNtB1J0pBZ32+L6iv1NgvFmWNms1B7el9NEqtfgPO4eiR9bDX5bODnJOiWpbFl/MYEmCLsUW8xub4uw5y5T6I7pv8eq394jepCqpZBDDuyesC96VbCTMCshrOp3kBaV1miESQ5iwpYUWuygtFnDuRUNKrFDqy0RaSJHODC2z4rBM3yO89AQz0w72vVte57h0n7MZMXEUV9IaZ3JjOr7QQU8XgOP1xqPVQqXUmMnNSqhQ7LWAhdKoBISChJKCzOIQaBTwKWU2MsKO1VjryvjaS92aMUNJCSUHu6PmQinzz3wy8hMmg85Ac05nAEgSXpp2ynPcIqw+2NWnqNP2PR+K10H94K6Jk49fiaMtABbhW6dOk6OrHMIiXz5splzQHmSzpWui9tKDCZTduIp9WwOa1V4c+BL1ULI3mbrq78vsetVwQ9qmfSA1hDYd8brWWsJhRpKabhFJBrs0doMaEo3qERtwwLT50kOugu8fSevgzMBYIdK1zam2oZoVcaZzCfqJ2qFV9Uaj9cdLqXChVTG41vofjrbaoG9kqj6AX4g7E4LNErgupP2atc2fv0GFXZo7WRCapnNF57ClJR9G5PNmFDj71PHxSghbCDj+0A5zqIDvCgYY2/230lPwiYcy/w+j/s9tUjSmTAUo8dZkTWw3FFqjpcvZTtydbj9h8ISVfuSeoCxU89I7VcoXafWFfalavfXSdXCkypj9bcj6sdqE650WYUe0MqYprGzqnGlJVqt0ekaWmsorbHDBRortXdiZ4na65dNQdoPtAnp2pxfnrDTcdVjm7X7GI3CsDjHDhVqaTy+H6tCon6iVnhi1+FVdYerqsOF7FBJjUoOfeqUxHVb9dK2kTakCenSwL4SuKgkdqrCTtdo7LUXlrCN7XrwCo8lUodsGtoEaW9F2Kn3NtWnVBpf93tuezFhA4TkPkHaFPLagPz1Dq5HgszvOhROO856K5wFWccOSiXlSzx4ZzubRaR9TiiZ5Cwh7BxRh9tDFXkfX90vVnGBC5up66qSuKpE71y1k+YzImtpeuBSbXZaQ6karVY2k3bd5xxXWo4mD0ZySTwHBGED40E4eBaiyY/7618L51zmUqi60LRaSFxKc96X1qHs8VrjVbXGq2qFJ+oWT+waPNi12FUd6kpBWOlaaQGlJC6aGrU0vzst0GqB1oZ0XVcCF1JgJyRqLVHp2kjYhCrc9Xuxw9SEP8mW0nbp5IEiRP+3vy2HPhokI2WbcsvOL1gfIeqv71gWTqBK79O50NI8sM36IAhTjo6zmFH2vrTqe1w2rwqf4+Wbw3TIx+mT+azJT4KwJ+9HRNy9Ctwjq0GyrEMHKxuqdFXBStbaekCjd5oSEGgroNECF51AKyUapUyol3aJVpwDleyTpDivcPpco6xUo8nL8nvrJixSDOlTDWkaj/faZii7rIArqfFAKlxVCo/XLV61a/H4RYOriwYXdYeqVhByIOuulZDSJBe66ST2yn0ELuQw4amE0WnU9ro4TYfDFiFcwJiwtyDokveUdAojQi5L3/20Q2JaLV5Sb86eX5LEZVRfNCaZ49To+b2LOHbo1m3hZCXrOPQl2DdzNjzZ1sY2tVOXwHPnWzI4ud+peiknQV8FDgwraEkdEtZOCvMRgyf0pRy8n51k3dl8tMZWC+ykcUTbSYm2MxKr9J3K4DmaeU5et7lEotTumgxOdDtpYqlracK0rqTCVdXhQW2I+uqywW6nUO0UIDWEMCpv1UgIodF2EldtjQtr275w19MRtTBOaEK7jHAuWYuZ0JzNyIXI0S3jcOYIG0BA2qu91RNq8aJjt9IsEPbrWIK+y85l9wmnSdYTZEL9LgmzAEICPaX1XreeMKSc8ZbY6qh6cu1N1S2d3Ra7XrKuAukSuKiAC2mI+kJq7KSJK3YLV3Ta9KjxpMZaCuyUxF5IVLqyTlTW61pIe0x4jn2/iVW4zI+xdE2f2NhOTZ+/r1kwVF1DQlhvd0fYO2G8vi+kwmWlcLlrcXXZ4PKqQ3WpUO10SNZ7c3KXTYeLpsOuUdj1YV7aLgACswhI12cn7zUcfv9ymgeHEqfGQ4ZXTvWFsrnHpJ1zGC0m3XgS6z1X1JiUwiaahxMaz44JVoOfGEpUrsbBI6GyJe2x9MtxrIf+UF6ypURNqeH8fUukhEDl7RzLPC9w16eBrIwyvLJE7STrCytZXkiNq0pjZ6VqZzTptAYg0WjguhO4qAR2ncZeClRKoNbSTgLqXmIc+jNTgiwl7FwVpImnsn1yxAlL1oaoL6XGZdXhsm5xuetwcdlh91iH6lJDXKCfuGgNyEpDa4WLfYuLmw6XlcKu1Z7625gQJMw1rGGkawhgoO4KsFnelghjU9ESKULa8j0IbLk5x60N0/nmxpyloDKubYczYadCaAhoiOmCmePPAWdB1imnpZLjFoU9RS/J5qE7ByLpXN0UoVJl5hJ2iqjJspFzlVukw61aXYuBqHcSuJQmEYgfW+xqbu14cyOrvvxFJbBX2iwlaclI2gVDKDWvczJL9ZXKbkYXnn8/fWlWwix7WUn0YWpGsjbS8UXVYbdrUV8p1I9ryEsBUcPoywFAAarSqJXC7rrDrlL2WNXb+iuBoQ0hIbSdLlnCltFkWNgFUMJ80zYH9Qx77MGk66gvU86iAKFtm5nON4eSVbr6siXarc2J+n5K3XcJZ0HWMVISISVdT9mSSgh8SThGCseQpvttzh5ZQNS5epfEccepN/vlLT3nKqmljfkdUoo6YjESoZGoL+Sg0nVz4EoLdNoRuehzZju7rAsUcxODwVYu0U2cV3Kg9c95gTQWh3T5q4g5qbcnVgHUjrArZezUl8oQ9ZWA2EkImyRGtwoSgG405E6hrjpUwhG1cRdzMeoCopfIfd1G3MdzBfWOUnkS4gmFKbeOtIuW1TzgRP2+QmOdKvtc9Ay3SNYKa1OOOpSEHE3Zi+L6qH2LMlodAUuJmjrPlK1tSkNBOpd5qu+4jETVk5V0DlbSqHadM9lAWvBsr/bVUjLI5FVL5zxl1eVCGAcu7S2k0V+TME94LpFF8n4X2v5L0C8W1pOoRm09vSuhIKWCrDTkDhA1DFHXQ+5VASP1i8rU4T7AsDKXC3tzkPB+uG1Cbma/2zoZCiVBl0jVwNjh07cpB+9AJl67FJSgkCvbf2cSX4z7YrO+hSckd2XsSzfhQFKqYip9yZwds7dnJuozdS6/ZId4IXNEnSo3dZ4pMqaOS5UfHzvYq8f7Hel4q0XBLVShIaxE7S8L6bJ3BeVdljMxWGEdYjUveW6kXbn8nvWS88T1GIXZFLeAkHH9OpQGlTbRDWRrpQ+yTSLmlyy3wUI0Rf2ZmEgHKv1CojfHLU+7uxapd245zoSZGCPcuhpcQ5sUkALw5w6j1Htxar5oVpz2/k7kLc4l+iC8SIPfGY/yLZbJW3N8Kj93ST1znMpy9SWvQaRuFZ562kmCguAirQXQLwU53r/1FGjJdRjVQfQqflYUFJQdPBWsV7fN521IVkBpAa1N4hcoQLcaYqcBZdcEVUYNjlZDt4DuBDplP/Z4V5/pw7rBeipEaksk2yoM2XIIpWdiqcyCvs9VlZcmE8q9P6P3llhnnHpW/Rjr+xC2xXHWR0GoCu8fMjc2Czkm7SgkY8qhKmmrzmVYSti7aXX7tFS5FLPtyxmintIQUAsRrCHuQfU8bjd3zXRPKo60BNyQY6RDt5616NMMUr10nD+yyYpq9ttZ4rMwV+Oi0MGsyG0IW2uNzkrAnY0j7z9KQHUCqhGQrbFNw01wtSFwtQdUA3SNRNtJNEr2WcyGerWRsjX6SULQpxn3eypcaksUOZERUnOqX/F69f4xW5L2HM1eSdtx+Vwq0vuE+6IGvyWyDjOZhRgI3CdqbQcnKm9x/LBP2Vf7XkwQeS4PcO74kra3HtRKQ7TSUu+wNm9ucpKsl1SPL3NcGqRKQ2BKm4UqgGFZSEdkThqd+74NKUfTjoeljoZzidqXuBxRuwGnVSFZN8okeWlbia4RqBoFLdGPMFoZ87u60ehuJJrGlleGsFvt6jMk7Yha22mC+UcP+scKYUxl85o6Zm0buXpLSHu1Q1pmcjuMacNCHqPjM/4WIdgT/C7g1tXgwFgV3hOzJ10PA+X0S1L64k95gqcSHsTHz0l+cBC7deJ6bGlvL+/LuA1Z2K62d9pJ0UYiNKQlIzbuvCUgnT3WfCyJrej/1IIvues4vubxohY25MijyhYKndbotIZbeKONCbutoJoW3Y05WrSuPkfWwpJ1jX1XodESjV1tq/MmAQrDNaLgzn2p1HYo6XpO+7l9g1ltLF1T9cw5lyUT+dn166lJVG7fmYiQM6HtvzXHnwNukayddD1WhQvPfp1LbD9XbeSOmyo7h7Dj+kv6sCVSC2nEv+MVoHz4tvhS6brUDl6KXv1tpUBDLgKNGvLHKwFoa9R2C1Q4KdR3otIeGeUkx2NibK/uoISy62/bhT010KhhSctGC9x0lVFrtxXavYR4ZGYlorKTEg3oVqDdC+xvKjRthX1X4aaTfR2NEr2ErbW2JgY7XRCq789tglJJ+9/7v2tzls+cTOTK5wSHJe9ESqp231ntTYPV4LeMQe1Nh1gsXipzwpZNHlNA2Lm2DoUcUU/Z8uNtax2qlsDZbDtt7LadUwEbfyncKBM3bcqasKzOU4M3WuBGiUF1rPQgZVsCHLXpSbYlmHNtSu+3ss+2hoIWyqr7da+ubhSwVwI3ncBeSlx3NW6aGrvrnel7pyKyNurv/b7Go32N67bCXsn+2jTKXlOl7VKi1k4u3FU6DTXpbTyDt4kl0jVjDHYwOyoUNGSvCtcag3StVeBoBixTM0/FDqfKzW3nGCgJzQr+FoQUzSXsEue+5LGDXIcWnSFqrdEo87lRwE1nlnSUwqi6W61R2/hp0z5wo8wiHu5jpFJHfKqXHLWVrmOpbO69Li6fUJNr3fWSrOuTgu7V4K3S2HcaeyVwrYBXOhNLftHUuJA7CKGhlMCu6SD6dasFus44lV3vd3i52eGVrsIrncR1J3Btr9GNcpMBc/U727bpQwe3dc55Hhpb262X4GAasYSJKj1+mbWrqYxyQ5nTuG+Mw+CWyTrnaDZI14g8ws2R3WhbCYnnwh6yzlQToRi542KsUeXNtddPbR/XX574ZQlc3b4auINCoxV2WqJRGvtO4Eaa1JuAyVTW2sxk/kIejTJ5wa879ETfKKCxdXai9Yjatasm1Ympa5Aj7VmqVZiBVwmFDg06dOi0MueugOsOdu1p2JXIKtTNBTQE2k6i3itU0qmHBTrrTPaoqfFnzQ5/1tZ4pZN4xV6b686YFfZKo1Hm2rTo0KGBimO+D6BqXeuINUcFvpWW6GgknRi/fBU4Bd88oLVCmJsi/n23wWrwW4DvaLZGugbG0s2U09BWs9IlIR1Fg9DMFa9SoVPDfn+t4bKFT9ZgFFpj1cCd6NBqowpvtELVCeykka6rzjBzo8wSj1IM3uDKLo153VnJugNabclIG2nVSK2N5zS10T1eOZAbonZhWwotWjTo0GiFGysRX0iXOlWiagavjkZJ7KRzqzPmAWPrlrjuKvxZW+PLrcTLrcAjS9Z7Bdx0jqjdtXFR3u7TnY1kdi79dJiaQE85fsbvZy6ULRA+AgXvmTDSArgohzXHnwNOgKxpR7N+rx4vyD6FQ3k9H2JBAlfv0lSHSxy9pmbtW0NBobLnaPKO7Tw1sDKSJSpD1tYh6qYTfbKuWhgbdiXQO5tpOMka2HeDVO1U4J2VHHuinqkyXKphmFZnOm3RoF3QUGi1QqNUrwq/liLIdV4JCaC2ExdNk7WSeKWVRgWuLFHb69Na80CjzZXphNVsnChJr+nT2onmVlL1mnq07gKnMqcCj7UMGr7fAXXOp3dvGctwAmRNIJCuB6LuBzq9zQs1J+XgXEyFhY3Kx7HdB5oYlOKQA7izkxqqaNBihz0qVFqi0xL7zq6y1QGAUQcPZG37B6MKH2zWRnJslOqlaucJrrXzuh6nx8xnmVpnEkh7EQ8DrDl/Q6CtVYXfdBoXlUDV+YuTmLpaDWu7H8QBF+J1Ham+na36ptN2QjOowBVC57I5CVGOiVIVeOm92lK9veYdLYv5HvsQ+M/vKU6ybgOpacqc488BJ0LWg+3aRJr6dmwrcVtVeH+EVl7CADpD0bExqZovcFSatWjAUgcvO2vfCnMSTigoCDEQlRKX1mbboYGEVEaSvO5cpjKgk2aVLbcmM2Ceik4bqfGmM85ZxttZoUGHFm1ks+6S6sQc5iyXOueZ096ERaNDiw6ttd3vlcB1a94Dl+9cCkDDxE7XAr1kDTiyNg53r3QCL7cCr7SGsG9684C5Nq215ZtJQtNPnNy5xTbsudhyEr1FP0oTJM0pU3LcknczFb4GIJhoOpNFL2lrfxJqf99htXcMtlnfCjxVeK8uHEK35uQLTw2qJQNBPItf6tSVK3/bs+KUI9FokQnCAW8q1p3OSNXZe2hU4VqbsKUOjSFpIQBtVoMSyi3CIfvY41qax6Fyj4cNc7rpNPZK40Yp3OjeXQuNaJwlGAqNJ6GF6sT4HKdA5YUfpVQtUYHb89dQaIWh6kZ32CuFXSeH1cVcjLk1D1xKsxqZn0N9CPkyUvUj+7luNa47hetW40Z19mpYSd5OZHybdXxe5PkL7308wDO8tQf4sfIbbNWPUOsTStBKtYNzom4Hwp5UgZ8JGzGyODGyNvCla8pmHZN2UC5DINTvEpXokuQHU304BcL2MUVe8f7iwcc5BXqTLQWFTjQQkJaoKzR6DwhAakPaUIBChU5JVFKgUnbpS2vLVjYMqbVEfa06Q3bosLe01FmSdirwEk/wOSBX6SoJW9MdIIBON+hEgw41GuzR4AI3ukKlBGRn7NEask/8slfGS7wScZ1+fDbwSqtx06En6r1S2Ntr40vV7vosjU4YJRDamBjnJkJZu4Rt4P+xkYq7SOOUESwCr2+CqKH9ief9k6oBACsdzM7lcp0QWecczQZVuPZs1ikbNiUBxvtG35eEZWW8sFPlb4uw5xCsf0xpveSkylO3KyjI3gehsyTUQAiJBjf9fez/WsJuoFAr2Ydt1Z5I2Tqy1oaob9CiQYNGNGjEzUDYVgXeqxFnkkAp5gzSwyIenVHai9YQtq6xV2adb6UllFZ25S2gqQSu7Rrf/kqZSg/S9b7TuLYfE7OtcK3aXqpusA/MAwqdmczokpC2sSPklgS9hVQ918/gFBOSjMKygDRRW6l6fG1ORxA4NNhmfcsYwrholWvJUncpafpQWJRi8AQI21cRp44rrWuyD5507duuBeQw8epNIUAFiQ4alRaQQqD1/BkUNBqtRkTdWau1H5LkHKpG57zRtZ+TqMdMXBS08AnbqKcbtKghIbSw2oXQHLCTRgUuvUmL0mYt60bBxmob9fdeGRt442Rp0aAV7SBVR9fklLDKm3utY+BGzp3kZD6zglv8TPrSNACSqKel6jMRGxmTOFmyHjBI2hoKQgOaSEMa29HmkHaJhBWrwksdSibrPTBRx5OanNbB/722XyP/AW1ye0srPcfSdYem75dA048xFSrUWqKChNQCwiMpbeOpnZ3aqXg70RoC1E3kBd6N7NUlWOKTAJhBf9rZz5gEKuyM7CsqVLoyZA1AKXefZK/qliPJWvfStTMLOKLe63ZkGuisY1lPBl7SmNS1id+tqWu4lcS6WEW/YM35pZiTcCiFHFH7zmRjor7fUjVgVOAcZ310lC3uMWwbL50J0IMJKR1urAY9FZRIvZS3aVm5YdIyUoET7aY8z520K4RR7QUBAF58Vt23WJmVsHUoWRuyNgTdYI9G7NHgppcclW6K7NTU87FUU1IC0zeTHEVA2knGDq026mrTJw2t3bUGdtrP4jaMMGbZS1PepGvteq/4vVN/W6l60DioQAUOYLUn+Bps7Vjm4xCRIZOJTuLJfCFRx2rvMqK+31I1q8FPAJQqHNHSmf3zKNyfKvtinxJBH039TZgM4v1+f6ac76YWUomXH3QOVbF0rbQhJfO9672jAdjvHTrsICFQ6RoS0nqJGxiybtGJzhC0p+JVukWnm8GxzA99idSLFNYSder6GOleetoF2x+h0MJoGSQkoC2hQ0MpjU5UaCxZAxhdB7dOtTMLOG3D3rPhu0lMh8ZoHeykwPUrd179RNjardeab5aaUfzj1+KQoVxBHZnwP8qJzN9eStQDzoV6toOZ1C6fmKw59pg4YbL2pWsFl37UT5TiO5YNDk7zyTh+8Q/ldJJybjtUW7GNH8irB/1+zXW+y0nXQXsJwpZiICiXilSJDhVaSC0hhUQFs+qUIyoztDmJsbHRyp6q10rVxtaXVoGX3P/tPZ1tCJvooKxkLex0pLFkrrRLmWrSsQotjD3bj9uCGWyGACyFxkxhem3DcH0GonaSva8CH52zm1QkNFa3Fb61BQ7xjuds1BTi8SBnnwbgTTjtYubBPuuROdS46lwYp4cTJOv84h6mBB3GRdmr+2NmDAC5erZAqi8UKWYl2AXSyZr487CertdikNebkK6FqJKEPUibbn8HhZ1xxhISEhVaS9wOfapOz2HKkZEjpM4RtaaJOsYWRO3uWclxSltHO2u7lqggRWWvwQ4dOlS6QocaFSpD6Tp8P5QdmN1KWo3Y9wt1+NqGIVRL2etTnmo0Rcy3EYK4tUZkC5RGAlB5DGKi9qXqWJoeEBN12Mp9AidFORHEqvA4jGsIE0pn/1nc9kg6Hdtr59SV3JeRXNes3V2KnH1/jqaCkuZLCRsaVrpUYRgXAI2qt6c64nZZt1zK0p6srXrXSdVDbvDwnEod6ebc4zn3SUFBwEm2EspY5HsJ20xWFGpde/nDJdy/uC73txNdoG3wr4/S7RDK1juYGcc7NcPxbmoFukPFW0/hmER9yIl8/6wGRD3sp4nalTkT5tkQLtJzzfHngBMla8rRzKnCx17gILzDN+tJYZ1Z9W9C/V1KhCnCLnEiWzOhoBa/KMlYFtwfyn4NkDZsd9sdeQ2/BxW+RoUORg3olr40jlLGA1z5al7CFj9HKhs7Cc1PBZvTOrg+OUc73zNeu0mhUKi0dbMTVU/UvinA/HVKcDV4w8P9bfrro61T2RxMSc9b2K/PBavs7ISt2peqRzZq+NdmbKMeEzXjLuNEyXoao5CkCSeqVW0F5JOWrud4tB7L0a10ElFST1xfTupPETaAgbR7idpImUZmNCFMgUocrQ3pMkQm4RKtDHZoJzF22qp6e9V3PlyrPOxmec729PUfUrB2GqjErveMd4O1srZ7icrYtPW4LiUG7YFP0s4sYCYEndU4eOTgSdXp8w5jlk9JHR63P7V9bf+WOKTl/ELibZQzmUEpUZ+LjLgtWA1+64gW97ASdLzOtZ/RDAhJYssXtZSwc8cPZ9aR2x1uQ40Yq4WzqnninCkHs5iwAQSkbTfYCoyUDd30hN2rxAEILT2SDiVrJy1qDM5knW6tKrxM3R1jyrN7rSex1gpKwJyjp1HoAEAbiVrBI2n7FwKBGrzXHvTkPiwN6kjaXaN+ItNLbePnMOlsFpmYplTht40S34PYWe4Qbabgv2cponYEHdurmahDMFmfDHxV+PB7iiQp6eaQhE216R/Xf/dCM3Lt5Oqb2rc1ZjufEf4DRaTtEbZABSlMog/A2KoVmoA4fZL2JcVOt8HgF6jBCwgmRdSz7NdZou686+NrGBpUAugwTE60UCONApUFy185KyZp6hqZY6bt9eegoi7VZKSOo8w/8QSt1IFsCpRJLJgopYiaJep7jzMgawPf0ayX4qCGmGvPdp3CFg5oOTXwpBq8gKiXYosBpESqHh1bqGWIpZjBiWasFnd2bAUjeQoh0WnXzthRLHaScr/92OpSzCHqJdfcf3Ypm70jbMCaBzyNghb0ebjzcwQNjCcyZhtN1JNOdivTd67FnJCxpDo8Wi++5Jg59ZeAes9GE0qWqGfDOJitiLPerisHxYmTdcHiHhhsoqVYS9qldtvgGIKo4/ZvI+xkjp19Sb25MDqftHqkHM9g73N075JEHUtLBedRYvecHcZVWD6WsAVMyJqLgxaQAJG4JPZyB0CStGtjDknnzklrtZkqfEn8dmxySZaLJeQMaaeO2QIjok5I1fTBp6/duE2wGvwE4duunSp8UpqO1KbBvhVqvpiw47aoNqfsp6X2761I/RgqzpxdMNCQeCt1kYRN1OdLiTFRL5GqfVCx0lPS9dzJVyBdR3kDIIZBRIrxM+Y/V753d4qkAc++XTRxqXp1/ciUcSLkMZeo431zJxqh1iWhQco4lMXmGIqwWf3NSOGMyDqUrs1AD6QczYD0rBrwBr+NCNuvM1XW/7sl5pL3IQfbOapwv3yckc7BLa1JISZpVx8dmxreH4qEciiRrkvIOfSKH8wH/cQE6CcrTpPQaZ+ox5J0eI40QftllqBEFb7mXUpJ10vTm+Zszr6THBA9FwVSdV44CK9TiqiH72GWMgpM1HnwQh4ng1C0SkrXXr7wEpW4P7NeS9iujtz+0Xc/HGaiv0udnHJ9mbNvaTtT3tCuTEzYZgf674pIaBerwd3AF+ybkKyn7JtUVEFOuppr13VaBEfg7jx9bYIhrLFDGdXmFEHn+ufa8MmGkq4Pial2SvaXbKc826f61X9Hvi4fUxJ1sC8hVa+xxd4XaDvlWXP8OeAMyBoYCJuWrv184RTxpUI2tgw/mWP7HaUcnGlz91FuEy0bbLcMx5njBBSYFdxxHpGFfQylx5GTzgotRlZ1KgZCo49Nk6oDmSTGC18TQo60CS4hyqguz4ZdStD+MVSdazHXE5vaPprIzZxYl05u55D/lMd50rM8o/r2peqQqCmwVJ0CS9YnCwUNGUjXgB0ErerQqQ9T+aupl3SO9DBHfeqXd/1ciqUhKltLRWtt5jExh6aEcbKb8fGh+ncLok6pTSnpOi7j98XUlVcZT6VhBcbahNgUkFNv5wg63k4R9il4fwOeiWQBYafqdJiyw0+R/BxQBB57gIeIpWomasZZkTW9wEcsXZuSylOtDoQNEIQXDQJ+vZM9ythn59RDYZDi1kk/pecx1Y+tQV27EXFHZJJK4rE1UefL5sm7lLAdKMKOtQlO2qaO7b8XEvRcOFV4KUreicV9makJmyNVF/ktzNB+xSaY3PrUBtRiHYwSrL1q53LFz4isQ1BZzbS7bToavArXuvaRmom7fT4x5MKTwj4n1MKgpbicOq9kcJlzrktNA1uEueSImyKKFFHPUY9myyWk6vF9Gk8kSgk7zg8enKs3J6XMAH4dwe/MtUqd4ykgZxbx/wLD8zmHbEts2Ye6Tr76myJqMhVu8Jul6inwetYnCcJ27RG22d4CoobSLYRdyahHMAhOx1v6SA0OJSQxByXq1hIcwyEo2/5MD9vwWNqGn5IkKSlmadvlNtdYbV6RhO1jNOGDR+qBndZ7JvOrxY4Ieqv7vkQVPjWBJdspJNIUYY/7Pe1zEPc31494X4nTZE/MxUSdC9ViMAzOjKwBn7A1JARET9j9ylwYQrkMkYeDv1OPx0guSpGxlW3lKUup56ntqW0Om/SlULqeul7x79mkHdm2TV1jlfewj0oaUpa0JqW5SEnVfTnfZ4IgbEd6lN3SN9u4usYanPx9KJW6R57rM9XVQshAY7XmuU/5AVDe9ZT9ekkbc3wOyHoKVeD+87meqFmqLgEnRTkTaOiesM3LjV6qBqQldAkXChR7mzpsnUChFPEkwu/LXNX31ECWc4wbaQ4WqrcPcZ0ooou3L2l3qYrTSdXBvbOkPMfGm5K0+/pmkuGUU1vc5y1Q8m6UTJjIKI4gjCxUiVPtFps5El77a1TewbMYfS9SfbvyTNSz4fzq1xx/Djgdw9UshA+0/4D7s9X45YizWs0Z7JdKumEdaZXfsYj63JBzoor3U78p+Nc7tX8JBul7/kRnrMr3QntyxxWW2wqHiECIJ6zh++B9z/h0LGrLTuRH5rKZIH0oME5568oMGEvVRO2L+8W4WzhjyTq0X5MOZ7aI0AjU4v3hFoM6skxlSiEVAjIqV5AadYkz2RbYSqV/SJRMtHIq3qlruXbiFdc1NZkosaMutR2n2ihFSltwqOeEeu5DW/5YFd7brzdIHFRip58yOwT2akKqHuooXaiDMQUN57O0/PhzwBmTtQ/a4WxE2L2dWyZJe6l91a+DQs7bVAhJEvWWHt+ldZZOOrZE0XkWhkGN6l5gr56D+HrNdczKEWsJ6W7jo7DNZHAJiedC4frtXupVyo9hy4U3Zt8/It599hrVmmOq1+C+qMHPnKx9UXnscDZ4iFvbtZA0aUf3ypF2bhGQuSiR9HLOMDG2Juk55Q/hUDdqI+ER7re/5hqUXI+5DlihVDwvPtkh5+S09LpTmobc83iIiVq5o19Ypl/720tDG0vb8/oxfq7G0vqYsLP+AN7zqHQTmdzayCmSiZqxDGdO1kDwFlsSBuBlOAulbDdDT6nGA3WbS6iyIhFDTqLuv6OatKOW1HksLCGP1EBdQqCxB/gSzCXqQ2AN2cY4xDlsUWcpgdK+H0Q2tSg/ua8V2+K6+BPkOfeHcigjVeAjaRpgz+9tYcb1dcefA+4AWTskbNhAIGWblxIjW3bvjUsQdtBKAfmWbO/3W6KmypMOKkewFR56MlA8KVlou9uKyJaENzlCWWr7zNVdevyxJyPxhHYu8ZHEXbigSGlbgfo8o7XJHZ/aRqq/neo7skszUW8PVoOfNTzCDuKwHQknnM+AXi0e28eA+USdA7WyU87B5lAo8W4/hiQ/9xoWSW+ENL+UyNYstrIFfNVvSdnb0h6U9iPrBDixApizYY+OW3DOVD8oOzT12ydp97uXpJ2tOijLRH0IKL2SrM9EtL5jZB3asGnCTjifiaiKje5f0SpOt6T+LvZ8LlR7FyVSWamOXnINtiSuvn2x3DErdw5TyXoob/FUG2vPO17dbFEdiQnvFFHHv6mkM379SS9tyiFNhCv0xSFzo20Jwg5I2leFk+pvJmrGOtwxsgbGhI3Iji3RZzrzCBtOfWmJus+elHiH0oPNPO/jmLBLllncAou0AUcK2aGw1A66FqkB3t+/hQd9bqKTilAoJe4Uxs5ceSeqOE9Bvu5x/vxVoZGFpoWcx/yI4L1EMVRmPP9vX5d3nwKSjj2/ZxE1Yw20/bfm+HPAHSRrINRpA1N2bNhsZ460zfcwZAR9jeELPwcpqTImavf9UIR9EFKb6S1/CmraOZhS6a6qu/Da5cIKKYk11eeiydEMqZpM8zrTWatcy5NfjCcnuVNe5D5px45iwJicg79+bLXn3a10C5aojweNdVOfc7kDd5SsHWi1OICEt7gjakI17tCv4DW2aaeQC5uhiJo8kyPbjW/b4/xU4a85vUQVnpPUUmWA6N5MkPbWE6EtHB2XOOuNju+TtExreHLXkCJtf7sh3260DQg1XjFJs+qbcUjccbIGKLV4mPFsUIsDGdU4QoKmljikkCJqSlo/BkolrTkq7zUx6FupwEf1ZtLDLmmPcjgEsMrxbI5ESredVo+XXddwDe2h3rFUfSgtzzEd92KzQc4WHav7R9K3R9LuNxP17YC9we8UfMIGUmrx2Fu8J2V7L7UYE7bzHB9aipKpTAyc1EIepyLVTjnuLKlvKxzKfj6LLD2yA9abLkqkRVK9O1PSnvKANmVCydJ9L9EMlGKptzhVdknbo22RwxjpbJYkaSCMpXblmagPDa1X2qzPxBv8vAyHqxDfkGjW62bD0QzZ9+6E933kWJKwcwHHJV8XAkaFgs3tS0ryCLYR5zrVxhzSLnJmmuG0lCPjErL0Ha3MtsNJnyn1c6qvlH01V5fZTkn3kQp4ROT0cz7dVrhsJFkmVjMX3Jfw+G75Ryso1Qb32H//tW7NODEaB/xxA0A/VjBR32X8zM/8DL72a78WV1dXeNOb3oTf+Z3fSZb98Ic/DCFE8Lm6uprV3j2RrB3SjmcA7S0OoFeRO1s2gMnY7LU5xoNeFw5Wa9XZpX3IOd/4v1eHDhUSDoAgMcfUNS9W70flgmsZ3HM/IUgVkF3Kg3pqQhfbSku8nF297vz7Z3FCunYaAvc9aN93ospMSFPnMQWtVX8tzbmM1fLDuUWq6In+LJmYBl7dQCBBj+uNJW/ffs1EfSzchhr8F3/xF/He974XH/zgB/GmN70JH/jAB/Dcc8/h85//PF7zmteQxzx8+BCf//zn+99CCLJcCvdIsvbh3xw3M7ZStna2JpvXN5o5p6TsWMKOB5FTUW1ToIiCImf/OxVjukXbU9t9jPwBUGXDm+b0KSXVBr8JIs55Ty+dJPh9yu0b+pULA0vd144kaqrdlDZlTr/T5btRP8n9if6Q/dXjd9N9lG6hbP5u/7sJv2rDd9xu88cGf3wYpG4m6mPCG4kXf+biX/yLf4Ef/uEfxrvf/W584zd+Iz74wQ/isccew8/93M8ljxFC4Omnn+4/Tz311Kw27ylZA+bliUkbkVrcbPdfSl81FtutkrGaM4htKfGtVS2XqBznTkCWkFMpUW6NLOFYNSlVjpI8105glvaVIjB6okVf7/i8Rs+xViQx5vpJ9p0ie9LeGw6p/nay3gRBp+rzSZgk6RKCjkiaifruY7/f43d/93fx1re+td8mpcRb3/pWfOpTn0oe9+Uvfxl/4S/8BTzzzDP4W3/rb+H3f//3Z7V7z9TgFCaynsHtHpzPAMDPfgbP2Sy39CbleDaVSSnoqU6rROfCtV9SLvibIYGpvvXnTLS9VMJeGsZTghGpRCrmoRPuz9hpK0X+qf7N7Wv8TKTua7x95AEekxpB1HP6mCLqkmxmQ5ay8Qpc8cTC2ZpT7ebeo2A/MUkZq7mHbQjqoP1h3F7GYRFP5pYcDwAvvfRSsP3y8hKXl5ej8n/6p3+KrutGkvFTTz2FP/iDPyDb+Pqv/3r83M/9HL7lW74FL774Iv75P//nePOb34zf//3fx1d/9VcX9fMeS9Y+aAk77XwG+GrxkYQdSdlTkmNa5bdeypwr+VD7S4m65Pectg8tTWfbT9lk0Y2vi3//KUkwkFjneU+XDkSxhOy3lbpflMQaE3VAihExxvXF/cifF6G6RvyudNGHIk/i/KN7QGk8fHMW5ThK26OH8sN50CY1t5dxeGylBn/mmWfw5JNP9p/3ve99m/Xx2WefxTvf+U688Y1vxHd913fh3//7f48//+f/PP71v/7XxXWwZN1jLGGbrX5u8fKYbFLKFiAlMz+LUrJ3BeSWk1bj44vzfScG/pJJQC6kraTtEmwlNZfW3zsPalqqjyXHkBznRQlMqX77NiMJOXZ0nHI4o+zY1AQjde+zz4KtO9YexdcxKJ/guBH5Rp7lsRmKOnbYkDqG/l0uSQ+lGMfBUruzfzwAfOlLX8LDhw/77ZRUDQBf+ZVfiaqq8MILLwTbX3jhBTz99NNFbe52O/yVv/JX8Id/+IfF/WTJOkDa8Sx2PgtDvPxwj3aYqUcSQ6y2IyUYSgoolK6pY6fKJq9ENFgn605Ik/H55drNnf/c6zIF3/6cLZe4htQ1yZ1Dv6/A1htf1/JzGkuSVF/jfdn+umc1ce8p+YTqU+p8cs9WUv7RfmhVNy6fOb5/d/VYknbvceyfEhP1OAd1LEkPZ8s4Pzx8+DD4pMj64uIC3/qt34pPfOIT/TalFD7xiU/g2WefLWqr6zr89//+3/Ha1762uH8sWY9AhXcBKVv2SMq227MreSVW9aLseXMlWGq/6UJIUKTtNT6+IOXi+JixJNVL8Smbb6btVBjWXKJeGvuenfREmpNUVqwtIwKoOihbNCVF+9vi4/y6k5PKKe1PwX6qb07Cdn1KXad4wutPhIP2U8cTGoT5knTuHJmkbwPK/ltz/Fy8973vxbve9S5827d9G77jO74DH/jAB/Dyyy/j3e9+NwDgne98J173utf1qvR/+k//Kb7zO78TX/d1X4f/83/+D/7ZP/tn+J//83/ih37oh4rbZLJOIp31DEBy2U33Pel8Fs8FJgb6op4WEO7Udj9OmarbDYy5dqccmOL48yIJ1yuzJGZ9E5L2zjtJNqZgsv2SSRfVXkmf/f2+02JAztYEE5M3VRe5eEVGxey3nerbMGEb+hb2uSOv37gfY6KmSDrnB2JAnc8SomaCvm1ooaHFcrJekv3sHe94B/7kT/4EP/ETP4Hnn38eb3zjG/Frv/ZrvdPZF7/4RUg5vA//+3//b/zwD/8wnn/+efy5P/fn8K3f+q347d/+bXzjN35jcZtCF+ZaE+I+83ocvO7n9RbWlj1sNwOTzSLmVvCCGejdtjjDWFKCzBBajrxyCUv64yc8t/3vJYM1Va8Irkl+QC+pz9RZRtpbeF1Tg358TqPvTnsww86fam+JNB4vDuNf+/jakZPEiBCpfsQalLl9i/s39Dn/HkwRNS09x0hI00E9U3ZpJukpmJC3w+Gll17Ck08+iW978v+HWlwsrqfVe3z6xQ/hxRdfDGzWp4b7zMAzQEvZZk8+vziZ9cytme1JOaP5gN96QoKdkkxz5DQn5IkasHMk5ktTpJQ30Vd3bLwvPu/Sycqo7wuIOj4nW2DYFuxPL6tY2l7qmBJHs0HiD6+9OT68dpSpJEXUZBhUtK109bi4f14F2eNGRD0p+ZdJ2OXSNJP0qUFjnYMZr2d950ARNjCyZY+W3bSqPyEHO7a3TWtaLWjqoO2Tc9SiQ+/HqtzpM15msxzFjVNLjc7sH2VrXeLRXXou1LH+b4q0c8fPchibSdJxmdFiM1Ef4wnFqK8RUQdtp54bT82dw6JngyBkR9Qpb27q+BBTJE3VdR6D+n2DgoI4ss36NsBkPQuxwRmYlV/cOZ7BlOmJGnSyCLIHOr2k4LSdjuh+pp1RvSUSaRSGRhLHVB+p/iVs+6WqfP93SgpL3QNKeqRIm5pk5NSyJSYB8prl7kNK8nd+ExNIEvXE9SD7RJ1ffD4JH464fHDvvO95G7SPuN05JE10jME4MpisF4FeEMTsCdXisbe4L2Wb8gNxm0MKbIAEmSXVy+RgP91ESpqaVMMGpFsmfY7bi+olVKaU9iGldYiJt1QVn5z8pMh5QoXrjqUy1qXaGDYUzv77Z4+49oKeUFDtptTMU30Nnt8S7U0/mR0mrsn6i0k60y6T9J0DFTo49/hzAJP1KsSJVIAptbj7HkjTnrTsyHuy2dGmxAM3kqpnrLlcSNJxmX7ApojD7p8z+McqU4q0gYT6OEHUS65BgIQGgawitkWnyC3Vr8JJRFBfdO37/omwbMpOXmYP9upKqMFz5xdP7nRK40P2Y4Kkp0xFwUuUKstEfQ5QQkGs8AZnNfi9QTrEK5f5DBikbrPRrzKt6k6S7VzJy6svd3yZd61XHUFiaWl7oj5CxZzLp54NSfLPIx7UU5OjXP88MvTrDMLVcuRRYuOdQdJxmYC0k85m9L2d57g1dios7XP2OSHKpaXpaZJOOxGxNM04DzBZb4ICwsZQJCZtf9tQY2KhDYLIl6pxRpJNkXMO1RYt3U4Sd6LNUdmESjcoS/SVkhrJicIcdXTUZnIiQtQdEClF3jNMGVNSbnDd/WgEU2G+HVLdDND3nibtUX8mJ2bpXUtIetrDl6XpuwJ2MGPMBG3HHiVR8YqYQUf25emMUgVOOiP4+zMq9WLb95z2hnZLHNuSg7hPxv2lTdjAC8PeKDIaN5uwm2YmIe44imyzNljQpF/sKJgA1Ve/f1PqZteHlIc17eCXfuZK7Mn5bA+UdD+edC0naHM04zzBZM1YiLQdG4hIGwiIhibnsqUswzrKjqcHyHlSHAVh1f5L2xuVdbZqirSByXE2VkeXqmZLHNFSttqpOoYK8ruXTZ6AkaTrL+NZ6GCGhddiqn+09F1SlpakyxzFkr2ZUZZximAHM8YK5GOyfYh44EpKn6WEHT94tAoyRahDeyV1h+3kjk8T+HRZHSktREDkhQ5zM5ylQqTKTXt0l9Sxvh8pe3Eo6Y7imifqSodC0VoUqux0OKL3XpQ+c0zSjHsKJuuDIRWTDQTE5g0aI6kb8Jyfls7+ClXUs44vb8dhur1hQjKW1hxxuI3uT8L22zca2zFjAhpPbOjtqf5SmJ4ElUqgOUx5vtPtlN73jNo5Uz6v/qbKl9SNAnt06XVkgr6LUOggZq4THx9/DmCyPjhypO0wxGj7IMl7KTJSV7pfEXJ9IVXtOQKfUu9SXsaG0FNSYondOJfpao4JIG2aiIhlFqmmMdcUER8Taium71X2OkVOdak6ss8WdT7xdVhF0kzM9wVm4dI1avDzeFaYrI8GirQdEpJJLHVPtkCQfb9zpsQ+lxySjmJjG/zksSOVuR9yREuJszKkzVaFT9VXgnL18fL6KKQmPTRmxy0niTvTx5K6/M1M0gwGk/XxUULaDrS6vLylGeS9oL5k3cFBC4hoRACRGhxjh6SyteO2slf7KJOwl5dZi4TtuvSYOfdvRtll0gyTNGMMTorCODBKYo/WPERUyFeGvDPlSjB1DNVWEfl7sdY522vKmzhor9gGGx+YsYuvMR1QmMxet+CZWGIjJ80JCyZrM+pII9VvJmiGs1mXR8xQx58DmKxPAoVBw7OQsx+6Vo832M1tKwh1AxKk7bDCYW6t5D/nsJLJ0lY+Cqk6Z04GSu9bqpw7R3r/mnNlombcLzBZnxwOQdxAcaKUyWNLMaeNafX/iLSBCYlxwonqEKSYQIrIqO2HkFCXmD+2msjNS/NZViODEWJdnPVxzFHrwWR90ogHpkNJ3YfAFi9PhrSB5erpDQgrR6priO4Q2g7y2q3qw7TWZt7xJWCSZtBQusP8ZzA+/vTBZH1WmDtgbSmZ3wbSpA0Q5JMLB1qo3k1hewep5YNNKZaf8xS5LtBsFIEJmsFwYLK+0ziUSv3YmA5tcygh8PmkdQwTQqpcrr154VvzsFb7Uno8EzJjHTjdKOOOYY1KfemAWtIGVXeJd3wi3/kqIp5bfgvpcU57W9VRGA99UDBJM7aBRge9QjOl2Ruccdo4xmC5tI1cLLpDqep1zvFzcR4z8jFuu99M1IztYOKkOc6awbhFzFXjl6irb/PFLCGpUzZXbOHwyETNYCwBkzXjTBCvZDaFU5otzyGoQ0UALMFUv5l4GbcPzg3OYJwcbsNh7jZs/VN1iIn9ufJL2mMwThdad9ArxgPNoVsMxiExxzFtSV1TZXNtHZrw5tbPBMxgnDuYrBl3CHMk77UEFpM2EyKDcRtgBzMG46yRsnFvTapM0gzGbcKEbq1Qg3PoFoNx22AiZTAYdwNM1gwGg8E4W2i9MoPZERf1WQMmawaDwWCcLe6LzfrwqwcwGAwGg8FYBZasGQwGg3G24DhrBoPBYDBOHJzBjMFgMBiME4dxMFsjWbPNmsFgMBgMxgZgyZrBYDAYZ4xupSKbbdYMBoPBYBwURo3NanAGg8FgMBi3DJasGQwGg3G2uC+SNZM1g8FgMM4WCgpi1UIe50HWrAZnMBgMBuPEwZI1g8FgMM4WrAZnMBgMBuPEsTZd6LmkG2U1OIPBYDAYJw6WrBkMBoNxtjC5vTk3OIPBYDAYJ4u1Nme2WTMYDAaDcWDcF7JmmzWDwWAwGCcOlqwZDAaDcbZYm9TkXJKiMFkzGAwG42zBanAGg8FgMBgnAZasGQwGg3G2uC+SNZM1g8FgMM4Ya8n2PMia1eAMBoPBYJw4WLJmMBgMxtmC1eARtG4P2Q8Gg8FgMGbjvoRusRqcwWAwGIwTB6vBGQwGg3G20HrlQh6aF/JgMBgMBuPA6ACIFcczWTMYDAaDcVAYB7HlZH0ukjXbrBkMBoPBOHGwZM1gMBiMM8Y6yZrV4AwGg8FgHBor1eBgNTiDwWAwGIwtwJI1g8FgMM4WeqUae+3xxwJL1gwGg8E4Y6gNPvPxMz/zM/jar/1aXF1d4U1vehN+53d+J1v+l37pl/AN3/ANuLq6wjd/8zfjYx/72Kz2mKwZDAaDwZiBX/zFX8R73/te/ORP/iR+7/d+D294wxvw3HPP4X/9r/9Flv/t3/5tfP/3fz9+8Ad/EJ/5zGfw9re/HW9/+9vxuc99rrhNoc8lyIzBYDAYDIuXXnoJTz75JIAaYk2cNTSAFi+++CIePnxYdMyb3vQmfPu3fzt++qd/GgCglMIzzzyDv/t3/y7+4T/8h6Py73jHO/Dyyy/jV37lV/pt3/md34k3vvGN+OAHP1jUJkvWDAaDwThj6FX/5oZu7fd7/O7v/i7e+ta39tuklHjrW9+KT33qU+Qxn/rUp4LyAPDcc88ly1NgBzMGg8FgnDnWK4hfeuml4Pfl5SUuLy9H5f70T/8UXdfhqaeeCrY/9dRT+IM/+AOy7ueff54s//zzzxf3jyVrBoPBYJwdLi4u8PTTT8PkBl/3edWrXoVnnnkGTz75ZP953/ved/RzyoElawaDwWCcHa6urvCFL3wB+/1+dV1aawgR2r0pqRoAvvIrvxJVVeGFF14Itr/wwgt28jDG008/Pas8BSZrBoPBYJwlrq6ucHV1ddQ2Ly4u8K3f+q34xCc+gbe//e0AjIPZJz7xCfzoj/4oecyzzz6LT3ziE3jPe97Tb/v4xz+OZ599trhdJmsGg8FgMGbgve99L971rnfh277t2/Ad3/Ed+MAHPoCXX34Z7373uwEA73znO/G6172uV6X/2I/9GL7ru74L73//+/E93/M9+MhHPoJPf/rT+NCHPlTcJpM1g8FgMBgz8I53vAN/8id/gp/4iZ/A888/jze+8Y34tV/7td6J7Itf/CKkHFzC3vzmN+MXfuEX8I//8T/GP/pH/wh/8S/+RXz0ox/FN33TNxW3yXHWDAaDwWCcONgbnMFgMBiMEweTNYPBYDAYJw4mawaDwWAwThxM1gwGg8FgnDiYrBkMBoPBOHEwWTMYDAaDceJgsmYwGAwG48TBZM1gMBgMxomDyZrBYDAYjBMHkzWDwWAwGCcOJmsGg8FgME4cTNYMBoPBYJw4/j8f6vFb1wJpJgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(recon_OSEM[:,:,61].cpu().T, cmap='magma', interpolation='gaussian')\n",
    "plt.xlim(30,100)\n",
    "plt.ylim(100,40)\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4dbfa37",
   "metadata": {
    "papermill": {
     "duration": 0.004738,
     "end_time": "2024-11-19T22:25:56.953094",
     "exception": false,
     "start_time": "2024-11-19T22:25:56.948356",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "You may wonder what the units are above. Since we scaled the system matrix by factors with units of Counts / second / MBq, the image has units of MBq $\\cdot$ seconds. Why?:\n",
    "* Image $\\times$ system matrix = projections and projecitons have units of counts\n",
    "\n",
    "So we can convert our image to units of MBq by diving by time, which in this case is the measured time per projection (15s):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "714b16f7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:25:56.962900Z",
     "iopub.status.busy": "2024-11-19T22:25:56.962434Z",
     "iopub.status.idle": "2024-11-19T22:25:56.977389Z",
     "shell.execute_reply": "2024-11-19T22:25:56.976745Z"
    },
    "papermill": {
     "duration": 0.021786,
     "end_time": "2024-11-19T22:25:56.978579",
     "exception": false,
     "start_time": "2024-11-19T22:25:56.956793",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.002"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pydicom.dcmread(file_NM).RotationInformationSequence[0].ActualFrameDuration / 1000"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pytomo",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 39.972697,
   "end_time": "2024-11-19T22:25:58.202095",
   "environment_variables": {},
   "exception": null,
   "input_path": "t_dualpeak.ipynb",
   "output_path": "t_dualpeak.ipynb",
   "parameters": {},
   "start_time": "2024-11-19T22:25:18.229398",
   "version": "2.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}