{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 3: Operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "from spectpsftoolbox.operator2d import GaussianOperator, Kernel2DOperator\n",
    "from spectpsftoolbox.kernel2d import NGonKernel2D\n",
    "from spectpsftoolbox.utils import get_kernel_meshgrid\n",
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Operator Intro"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operators in `spectpsftoolbox` are used to convolve kernels with various inputs. Let's start by recreating the ngon kernel we made in tutorial 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "collimator_length = 2.405 \n",
    "collimator_width = 0.254 #flat side to flat side\n",
    "sigma_fn = lambda a, bs: (bs[0]+a) / bs[0] \n",
    "sigma_params = torch.tensor([collimator_length], requires_grad=True, dtype=torch.float32, device=device)\n",
    "# Set amplitude to 1\n",
    "amplitude_fn = lambda a, bs: torch.ones_like(a)\n",
    "amplitude_params = torch.tensor([1.], requires_grad=True, dtype=torch.float32, device=device)\n",
    "\n",
    "ngon_kernel = NGonKernel2D(\n",
    "    N_sides = 6, # sides of polygon\n",
    "    Nx = 255, # resolution to make polygon before grid sampling\n",
    "    collimator_width=collimator_width, # width of polygon\n",
    "    amplitude_fn=amplitude_fn,\n",
    "    sigma_fn=sigma_fn,\n",
    "    amplitude_params=amplitude_params,\n",
    "    sigma_params=sigma_params,\n",
    "    rot=90\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this operator takes in `xv`, `yv`, and `a` and returns the correpsonding kernel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nx0 = 255\n",
    "dx0 = 0.048\n",
    "x = y = torch.arange(-(Nx0-1)/2, (Nx0+1)/2, 1).to(device) * dx0\n",
    "xv, yv = torch.meshgrid(x, y, indexing='xy')\n",
    "distances = torch.tensor([1,5,10,15,20,25], dtype=torch.float32, device=device)\n",
    "kernel = ngon_kernel(xv, yv, distances, normalize=True).cpu().detach()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operators can be constructed using kernels. For example, we can create an operator using a 2D kernel as followS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngon_operator = Kernel2DOperator(ngon_kernel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operators apply convolution with the kernel to some input `input`. \n",
    "\n",
    "* We'll create an input of shape (12,255,255) corresponding to 12 different `a` positions. The input will consist of two squares in each plane. We'll see convolution with the `ngon` kernel in each plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nx0 = 255\n",
    "dx0 = 0.48\n",
    "x = y = torch.arange(-(Nx0-1)/2, (Nx0+1)/2, 1).to(device) * dx0\n",
    "xv, yv = torch.meshgrid(x, y, indexing='xy')\n",
    "distances = torch.tensor([1,25,40,55.]).to(device)\n",
    "input = torch.zeros((distances.shape[0],Nx0,Nx0)).to(device)\n",
    "input[:,120:130,120:130] = 1\n",
    "input[:,120:130,170:180] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The meshgrid `xv` and `yv`above correspond to the dimensions of the input. Typically, the kernel need not be as large as the input. We can manually define the meshgrid of the kernel `xv_k` `yv_k` or use built in functions to obtain them if we know the spatial extent of the kernel. **The thing you need to enforce is (i) The dimensions of `xv_k` should be odd and (ii) that the spacing in `xv` and `xv_k` is the same**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_width = 24 #cm, need to manually determine this\n",
    "# Automatically makes (i) same spacing as xv and (ii) odd kernel size\n",
    "xv_k, yv_k = get_kernel_meshgrid(xv, yv, k_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use our operator on the input. Again, note that `xv_k` and `yv_k` refer to the kernel used and not the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngon_output = ngon_operator(input, xv_k, yv_k, distances, normalize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAIfCAYAAABtmb+NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABArElEQVR4nO39ebDleUHf/78+n7Pc/d7eZ9+CMzBDgKGionEZzYAIAwYxGrcSBhWjUMRQlkYUl1SxWJGKfA0W+jNsUmAEg0aRSSTiVqCioBFZAgoMzEzv3Xc/6+fz++PcvtPNzPT0DP3pjcej6tLd537OOe97mPc95/n5vM/nFHVd1wEAAADOuvJ8DwAAAAAuVaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCFfUtH9pje9KUVRbH9NT0/n8ssvzzd+4zfmVa96VQ4ePPiA6/zcz/1ciqJ4RPezsbGRn/u5n8sf//Efn6WRn39vectb8p3f+Z157GMfm7Isc/311z/i2/jlX/7lPO5xj8vU1FRuuOGG/PzP/3yGw+HZHyyXLHP40bnvvvvy0z/90/nqr/7q7NmzJ4uLi/kX/+Jf5Nd+7dcyHo9P2faP//iPT3mMT/76i7/4izO6v4MHD+b5z39+9uzZk9nZ2Xz1V391/s//+T9N/GhcAszrs+PAgQPZvXt3iqLIO9/5zgd8f21tLT/6oz+aK6+8MtPT07n11lvzm7/5m2d8++Y1j4R5/ehdf/31D/oc/O/+3b87ZTvP1xeX9vkewPnwxje+MY973OMyHA5z8ODB/Pmf/3l+4Rd+Ib/4i7+Y//7f/3ue+tSnbm/7Az/wA/nmb/7mR3T7Gxsb+fmf//kkyTd8wzeczaGfN7/xG7+R/fv35yu/8itTVdUjjuVXvOIVefnLX57/+B//Y77pm74pH/zgB/PTP/3Tueeee/Jrv/ZrDY2aS5U5/Mj8zd/8Td7ylrfk+77v+/Lyl788nU4n73nPe/LDP/zD+Yu/+Iu84Q1veMB1XvnKV+Ybv/EbT7nsn//zf/6w99Xv93P77bfn+PHjee1rX5t9+/blda97Xb75m785733ve3PbbbedtZ+LS4t5/cV50YtelOnp6Yf8/nOf+9x88IMfzKtf/ercdNNNedvb3pbv+q7vSlVV+e7v/u7T3rZ5zaNlXj86X/M1X5Nf/MVfPOWyyy677EG39Xx9kai/hLzxjW+sk9Qf/OAHH/C9z372s/U111xTLyws1Pv37/+i7ufQoUN1kvpnf/Znv6jbuZCMx+Ptv99xxx31ddddd8bXPXz4cD09PV2/8IUvPOXyV7ziFXVRFPU//MM/nK1hcokzhx+do0eP1oPB4AGXv+hFL6qT1Hfffff2Ze973/vqJPU73vGOR3Vfr3vd6+ok9fvf//7ty4bDYX3LLbfUX/mVX/mobpNLm3n9xXvnO99Zz8/P129+85sfdP6++93vrpPUb3vb2065/GlPe1p95ZVX1qPR6LS3b17zSJnXj951111X33HHHQ+7nefri8uX1PLy07n22mvzmte8Jqurq/nVX/3V7csfbKnLH/3RH+UbvuEbsnv37szMzOTaa6/Nt33bt2VjYyOf+cxnsnfv3iTJz//8z28v83j+85+fJPnUpz6VO++8MzfeeGNmZ2dz1VVX5dnPfnb+/u///pT7OLFk5O1vf3t+6qd+KldeeWUWFxfz1Kc+NZ/4xCceMP677rort99+e5aWljI7O5ubb745r3rVq07Z5q//+q/zLd/yLdm1a1emp6fz5Cc/Ob/1W791Ro9PWT76/1Tuuuuu9Hq93Hnnnadcfuedd6au6/zO7/zOw97GPffckxe+8IW55ppr0u12c+WVV+bf/Jt/kwMHDiS5//F629velp/4iZ/IFVdckfn5+Tz72c/OgQMHsrq6mhe+8IXZs2dP9uzZkzvvvDNra2uP+mfiwmMOP7SdO3em0+k84PKv/MqvTJJ8/vOff9jbOFPvete78tjHPjZf/dVfvX1Zu93O937v9+av/uqvcs899zzsbTzcY/H85z8/8/Pz+fjHP56nP/3pmZubyxVXXJFXv/rVSZK/+Iu/yNd+7ddmbm4uN910U9785jeftZ+Pc8u8fnhHjx7Ni170orziFa/Itdde+6DbvOtd78r8/Hy+/du//ZTL77zzztx77735y7/8y9Peh3nN2WReXxjM63NLdJ/kmc98ZlqtVv70T//0Ibf5zGc+kzvuuCPdbjdveMMbctddd+XVr3515ubmMhgMcsUVV+Suu+5Kknz/939/PvCBD+QDH/hAXv7ylydJ7r333uzevTuvfvWrc9ddd+V1r3td2u12nvKUpzzoxH7Zy16Wz372s/n1X//1/Nqv/Vo++clP5tnPfvYp78P8b//tv+WZz3xmqqrK61//+vze7/1eXvKSl5zyQvp973tfvuZrvibHjx/P61//+vzu7/5ubr311vzbf/tv86Y3veksPYIP7iMf+UiS5AlPeMIpl19xxRXZs2fP9vcfyj333JOv+IqvyLve9a689KUvzXve85780i/9UpaWlnLs2LFTtn3Zy16WgwcP5k1velNe85rX5I//+I/zXd/1Xfm2b/u2LC0t5e1vf3t+/Md/PL/xG7+Rl73sZWf3B+W8M4cfmT/6oz9Ku93OTTfd9IDvvehFL0q73c7i4mKe/vSn58///M/P6DY/8pGP5IlPfOIDLj9x2T/8wz+c9vpn8lgkyXA4zHOf+9zccccd+d3f/d084xnPyE/+5E/mZS97WZ73vOflBS94wfYLiuc///n5m7/5mzMaPxce8/r0XvKSl+SGG27Ii1/84ofc5iMf+UhuvvnmtNunvqvwxLx8uOdh85qzzbw+vT/90z/NwsJCOp1ObrnllrzmNa95wDlYTvB8fZE434faz6XTLXU54bLLLqtvvvnm7X//7M/+bH3yw/TOd76zTlL/7d/+7UPexiNZ6jIajerBYFDfeOON9X/4D/9h+/ITS0ae+cxnnrL9b/3Wb9VJ6g984AN1Xdf16upqvbi4WH/t135tXVXVQ97P4x73uPrJT35yPRwOT7n8Wc96Vn3FFVecsnz84TzS5eU/+IM/WE9NTT3o92666ab6m77pm057/Re84AV1p9OpP/rRjz7kNicer2c/+9mnXP6jP/qjdZL6JS95ySmXP+c5z6l37dp1hj8BFwpz+OzM4bqu6//1v/5XXZblKWOu67r+0Ic+VP/7f//v63e96131n/7pn9ZveMMb6ptvvrlutVr1XXfd9bC32+l06h/6oR96wOXvf//7H3R568nO9LF43vOeVyepf/u3f3v7suFwWO/du7dOUn/oQx/avvzIkSN1q9WqX/rSlz7s2Dk/zOtHP69///d/v+50OvXf//3fnzK+L1xueuONN9ZPf/rTH3D9e++9t05Sv/KVrzzt/ZjXPFLm9aOf1z/yIz9Sv+ENb6j/5E/+pP6d3/md+nu+53vqJPX3fu/3nrKd5+uLiyPdX6Cu69N+/9Zbb023280LX/jCvPnNb84//dM/PaLbH41GeeUrX5lbbrkl3W437XY73W43n/zkJ/Oxj33sAdt/y7d8yyn/PrH36bOf/WyS5P3vf39WVlbyIz/yIw95xsdPfepT+fjHP57v+Z7v2R7Dia9nPvOZue+++x50j9/ZdLqzUT7cmSrf85735Bu/8Rtz8803P+z9POtZzzrl3yeuc8cddzzg8qNHj1pifgkyhx/ehz70oXzHd3xHvuqrvuoBS+Ke/OQn55d+6ZfynOc8J1/3dV+XO++8M+9///tzxRVX5Md//MfP6PYf7Xw/k8fi5Nt55jOfuf3vdrudL/uyL8sVV1yRJz/5yduX79q1K/v27dt+vLk4mdcPtLy8nB/6oR/KT/zET5zRSZO+mOfhL+b65jUPxbx+cK973ety55135uu//uvzr//1v85b3/rWvPjFL85b3/rWfPjDH97ezvP1xUV0n2R9fT1HjhzJlVde+ZDbPOYxj8l73/ve7Nu3Ly960YvymMc8Jo95zGPy2te+9ozu46UvfWle/vKX5znPeU5+7/d+L3/5l3+ZD37wg3nSk56Uzc3NB2y/e/fuU/49NTWVJNvbHjp0KEly9dVXP+R9nnjf84/92I+l0+mc8vUjP/IjSZLDhw+f0fgfjd27d6fX62VjY+MB3zt69Gh27dp12usfOnTotD/fyb7wtrrd7mkv7/V6Z3S7XBzM4Yf34Q9/OE972tNy44035g/+4A+2x3M6O3bsyLOe9az83//7fx/0ZzzZ7t27c+TIkQdcfvTo0SQPnIsnO5PH4oTZ2dkHnKm52+0+6O13u11z/SJmXj+4n/qpn0qn08mLX/ziHD9+PMePH9/ekbyxsZHjx49vR80XMy+/2Oub1zwY8/qR+d7v/d4kediPAvN8feH6kvzIsIfy7ne/O+Px+GE/cuDrvu7r8nVf93UZj8f567/+6/zyL/9yfvRHfzSXXXZZvvM7v/O0133rW9+a7/u+78srX/nKUy4/fPhwduzY8YjHfOIEEqc7EdKePXuSJD/5kz+Z5z73uQ+6zWMf+9hHfN9n6sR7uf/+7/8+T3nKU7Yv379/fw4fPvywe+j37t17Vk/0xKXLHD69D3/4w3nqU5+a6667Lv/7f//vLC0tnfE4T7x4f7g92k94whMecJKaJNuXnW6+n8ljwZce8/rBfeQjH8lnPvOZXH755Q/43vOe97wkybFjx7Jjx4484QlPyNvf/vaMRqNT3td9JvMyMa85+8zrR+bEc/CZnNjY8/WFyZHuLXfffXd+7Md+LEtLS/mhH/qhM7pOq9XKU57ylLzuda9LMlmymTxwz9jJiqJ4wJGld7/73Wd0hsAH8y//5b/M0tJSXv/61z/kMp3HPvaxufHGG/N3f/d3+fIv//IH/VpYWHhU938mvvmbvznT09MPOHnEm970phRFkec85zmnvf4znvGMvO9972t8CTwXN3P49HP4b//2b/PUpz41V199df7wD/8wO3fuPOMxHjt2LL//+7+fW2+99bSfA5wk3/qt35qPf/zjp5wNeTQa5a1vfWue8pSnnPaoxpk8FnxpMa8fel7/0i/9Ut73vved8vVf/st/STI5C/T73ve+zM/PJ5nMy7W1tfz2b//2Kbfx5je/OVdeeeUpO8QfjHnN2WReP/LX3G95y1uSJF/1VV912u08X1+4viSPdH/kIx/Zfn/FwYMH82d/9md54xvfmFarlXe9613be28ezOtf//r80R/9Ue64445ce+216fV6ecMb3pAkeepTn5okWVhYyHXXXZff/d3fze23355du3Zlz549uf766/OsZz0rb3rTm/K4xz0uT3ziE/M3f/M3+c//+T+f8fLpLzQ/P5/XvOY1+YEf+IE89alPzQ/+4A/msssuy6c+9an83d/9Xf7rf/2vSZJf/dVfzTOe8Yw8/elPz/Of//xcddVVOXr0aD72sY/lQx/6UN7xjnec9n4++tGP5qMf/WiSyRHqjY2NvPOd70yS3HLLLbnllluSJH/yJ3+S22+/PT/zMz+Tn/mZn0kyWZ7y0z/903n5y1+eXbt25Zu+6ZvywQ9+MD/3cz+XH/iBH9i+7kP5T//pP+U973lPvv7rvz4ve9nL8oQnPCHHjx/PXXfdlZe+9KV53OMe96geOy5e5vAjm8Of+MQntn+2V7ziFfnkJz+ZT37yk9vff8xjHrP9mH33d393rr322nz5l3959uzZk09+8pN5zWtekwMHDjxgx9n3f//3581vfnP+8R//Mdddd12S5AUveEFe97rX5du//dvz6le/Ovv27cuv/Mqv5BOf+ETe+973npXHgkuTef3I5vWtt976kN97/OMff8oRxGc84xl52tOelh/+4R/OyspKvuzLvixvf/vbc9ddd+Wtb31rWq3W9rbmNWeTef3I5vXb3va2/I//8T9yxx135Lrrrsvx48fzjne8I7/5m7+Z5z//+XnSk560va3n64vMeTh523lz4kyKJ7663W69b9+++rbbbqtf+cpX1gcPHnzAdb7wTIof+MAH6m/91m+tr7vuunpqaqrevXt3fdttt9X/83/+z1Ou9973vrd+8pOfXE9NTdVJ6uc973l1Xdf1sWPH6u///u+v9+3bV8/OztZf+7VfW//Zn/1Zfdttt9W33Xbb9vUf6gykn/70p+sk9Rvf+MZTLv+DP/iD+rbbbqvn5ubq2dnZ+pZbbql/4Rd+4ZRt/u7v/q7+ju/4jnrfvn11p9OpL7/88vpf/at/Vb/+9a9/2MfuxOPwYF8nnzHyxLgf7CySr33ta+ubbrqp7na79bXXXlv/7M/+bD0YDB72vuu6rj/3uc/VL3jBC+rLL7+87nQ69ZVXXll/x3d8R33gwIHTPl4PdfbMEz/PoUOHzuj+uTCYw49uDn/h4/aFXyeP5VWvelV966231ktLS3Wr1ar37t1bf+u3fmv9V3/1Vw+43RNnJf30pz99yuX79++vv+/7vq/etWtXPT09XX/VV31V/Yd/+IenHeMjeSye97zn1XNzcw+43m233VY//vGPf8Dl1113XX3HHXec8f1zbpnXj/65+Qs91PjqenK24Ze85CX15ZdfXne73fqJT3xi/fa3v/0B25nXnA3m9aOb1x/4wAfq22+/ffv17uzsbP0VX/EV9a/8yq884Kznnq8vLkVdWxMAAAAATfCebgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGhI+0w3LIoz3hQ4x+p69KiuZ17DhevRzuvE3IYLmedsuPQ83Lx2pBsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABrSPt8DoDnjt7w4efwNp93mK27/WD50/NfP0Yge3nN3/WTe8Yd7TrvN+N0fTvdn3nqORgQXFvMaLk3mNlx6zGtOEN2XssffkOpJTzrtJjP1p8/RYM7MXKd82DG3P33PORoNXIDMa7g0mdtw6TGv2WJ5OQAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADTE53Rfwp5y+8czXX/2tNv85frbztFozsw7jr0xn969etptlotj52g0cOExr+HSZG7Dpce85oSiruv6jDYs9DlcqOp69KiuZ17DhevRzuvE3IYLmedsuPQ83Ly2vBwAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCFFXdf1+R4EAAAAXIoc6QYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaEj7TDcsijPeFDjH6nr0qK5nXsOF69HO68TchguZ52y49DzcvHakGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABrSPt8DAAAAuLQUW3+WJ11SpE590jbV1p8nX3a+XazjvrCJbgAAgLOiSFKmSJEU5eTvW3+e+G5Spa6rJFVSV1tBez5D9v7Qnoy7vTXmyWX3b/FQ4xbfD0d0AwAAfFFOju12iqKd8pQ/7z9yXNdV6lSp6lHqk74mITvOuY3YIkVaSVGm2Brv9pi3dxjcP+4TY67qUbL1Jb4fnugGADivJi/WJ38rHnSL838k7At9wZGxB3H/clQvxrnUnQjXrWAtu2mX02mVU2mV3ZRlJ2XRTrk1z8f1KFU9zLgaZFz1M64GqapBqnowidh6nObnzdZOgqK1FdndlGU3rbJ7/7iLTlrFJBerrZ0EVTUZ92i8mXE9GfdkzMNzMOaLl+gGADgvHuQIU8rJv7denNepto6CVVtHws73UaUHP5r3hWPO1pG883sED86F+4O7VU6nVU6n055LpzWXbmsunWI23WI2rXTSSjtVqowzzCj99Ku1DKq1jMabGY7XMxr3Mq56SQap66S5eX5ycHfTKqfTbk1vj7nbmk+nmE07U2mlkzJlxhllnGEG9UYG1VoG1XqGo/VJfFe9jKts/X4yzx+M6Aa46Dz4UaX7XYhPdhfjmKFJrVNe8J58dKlVTqUoyu0lqOOqv/XCdpBx1bv/yNJ5W4Y6iYvJ0byZtFpbY95ailrX1fbRu3HVz2h7zIOtI3jjczhmaNKpc6Ldms1UZylTrcXMtHZmLjszVy1mpp7eSu4yk91PVTYzyEaxlrX28Wy2jmVzdDyDYjXDcZnROGkuvE8N7nZrNp3WXKY6S5luLWam2Jn5ekdmq/nMpJtWypQpMkqVfkbpFZvZaK1lvXUsm+Wx9EbHMxhtHcGveinqCO8HIboBLgqnLl0rTjqqlJw4GjY5uUm9vTTtQnhh29peurZ9FO8kJ8Z9/xE8T9R8KWilKDrbL9K77YVMtRczVS5mqpjPVGa358ow/fSzll61nP54NYPhaobj9YyrjaQanMM5M/kdVJTdtMr7X6RPtRYzVc5nKvPpZCpJMs4ww/TTq1cyqNbSGx3PcLSe4XgtVdXbCokL4fcTfDG23hbyBcE919mb+WJfdlZ7syNzWWx3M9MuM90q0iqTuk6GVZ3N0WyWh3NZrpayXCxmuT2dja2l3JMdV0nqXjPHuU8K7m57IdPtpcy192Wx3pud1a4sldNZ6LQz0y7SbU12mo+rpDeusz6ay9poKUezmOOt2ZRFZ/t933WqVFUmP6Sl5qcQ3QAXvK0X6K3ZdFsL6bRn0ylnUhad7S3qVBlWmxmNNzMYrWY03khd9VJndJ7GfOIF+v1R0WnNpb11NOyEUd3fXlY3HG9kPN5IXffjiZpLV5GiaG0tQZ3PdGdX5tp7s1Rcnh3VzixkJtNlK51ickSsNx5nvR7keLGSlc6hrJb7szk8mv5wK1vPSXifGtzT3Z2Zae/KfGtfdtSXZUe1kLmym245OSI2rKtsjkdZyWaWW0ezUh7MenkoGSTDZCu8vSDn4ndip3KrnE63vZC5zt7sKK7Knmpf9rZns3uqnZ1TRRY7yWy7TquY9Gi/KrM2So72yxztdXJk0E2r6KRol6lO7IxOlXE9OssBuzWXt9533mnNZbq9lPn25dlVX5m99Y7smZrK7ulWdnST+XYy3Zrc77BKNsdlVoZljvXbmet1MjWeSqtsJ+1M3u+99VaYyU50Tia6AS5oRYqik057MbNT+7LYvjJLuSyL9VKm091eqjbMOGvt9RzrHMzK+N6s9e9Lf3g8dbWZ83NEqUzZms1Ue2fmpy7PUuuq7Kj2Zi4z6aS1vVStVwyy0lrOsdyb1eG92egfymgc4c0l6sQL3skRpunOrix1rs7e6prsy1L2znSzc6rMQifplsm4TjZGnSwPpnKoN5NDo4Xsb012XNV1lf6oyriukrrpOT45mleW05nqLGW2syc7y2tyeXVlLuvMZtd0Kzu6RWZaSVkkgypZHXZzrD+dg725HMx8DnU6k5Ut/SrD7fd4n6+dgnA2TD4GrCy7abdmMt3ekfliX/ZU+3Jldz5XzLZy+Uyyb2qcnZ1xFjqjdIoqVZKNUTvHh+0c6rayv9NKd2M6RW93Uifj1jBVZ3KitboepapOrGQ7S7bOIdHeeu/5bGtPdtVX5vLsyhWz3Vw2U+bKmTq7u+Ps7A4z3RqnTDKsy6wO2zk6aOVAt5WZdiet9flkdHnG5Sjjdj9VNQnv+3cWWNFygugGuIAVaaXVms3s1L7s7dyU66sbctXMTPbNlJlvJ1OtZFQlm+PkSH8h923szmfK3bl3autF+XB0Ho4obe0oaM1ncfqaXFk8LtcXe3PFQju7p+9/Yd4fJyvD5ODmjtzT25vPdOdzOMl6b5TReBhP1lySismL9G57IXPtvdlbXZNrWrtyzXwn18zVuWxqnB2dUaZb44zrIqujdg7227lno5PZ9fmUvatTt6qMOv2Mq/5J7+9uLmCLFCmKdjqt2Uy3l7JUXpUrq6ty7cxcrp0vc8V0lb1Twyy0xymLOr1xK8eH7dzXa2VxYzrdtd2pqjqj9mTM43qw9aK8iJ1rXJyK7XlxIl5nWjuzs9qbve3ZXDHbyrVzda6bHebq2c3sntvI7Mwgne44dVVkc7OTY+sz2b8+m5nWVFpFmaqeyrC3I4NyM8PWRkatzYyqXk7E/dl5TtzaUVBs7Sho7chCsTe7q6Xsm+nmqrky18+Oc81sP5fPbmTn/GampkcpyjrDQSvrG1M5vD6TxY2ZdMtOknaGa7PpVbvTL9cybE9W3FXFid9L5vgJovsSNn7Li5PH33Dabb7i9o/lQ8d//RyN6OE9d9dP5h1/uOe024zf/eF0f+at52hEcJ5tvdBdbF+Z66sb8s93zOVxi1Wun93MnuleZjqjDKsyK/1uPr8xk/+31k3n6OUZFJvpd1YyHK9lfM4DtkyrnM50d1f2Fjfkps5l+ec7W3nswjBXzvSy2O2nVdZZH3RypD+VT69P5WMr86mX/1l6nZX0h8sZVxtb70uHS8n9R5im2otZKi7P5cWOXDPfyeMWq3zZfC9XL6xl59JGujOjVMMy62tTuW95PkvtubTLdobVXDYHl2WztZxhez2jqpdxPWgwYIvtHQWd1uSo2O7qslwxNZsbFsrcND/Kly2s5bKltczND1K0q/TXOzm2MpvPrc5ntjWdqu5ksLorG8VK+q3lDIrVVEX7Ajr3BDwKRZmyaKdVTmWqtZC57MzOYj57p9u5Yia5bnaYGxdXcsW+lcxdMUp7dyfFdDv1qMriai9L925mfv8gZbGUcT2d3rjM5mg6a6MdWS+X0m+tpjVaT1WcOA/CWRjy1o6C1oml5eViFqsd2dWZyr6ZMlfNVLl+rpd/tut4dl+5ke7lRcqFToqySLUxytLRlSzdu5npQ0tJ5jOoutkcdbK+vpi1ckd65XIG5erkbOZFeQ5W4Vw8RPel7PE3pHrSk067yUz96XM0mDMz1ykfdsztT99zjkYD59uJo0tzWcpluWpmJjcvVvny3cu5/pojmbu6TrnQTj2oMjw4zlWfmc/UwT1ZHXZz4Ni+HGnvyEZxKOOsn+tRpyy7mW7tyL56b25YaOVJS/088bLD2XP9ejr72inaRcbLo2x8vszuz+1OlaUc689m/+DyrLQ/n8HoeOp6cE7HDedCceJFermYHdXO7JmZHOH+svlebt53JLsfs5nONdMp5uaTUZXFg5tZ+KfNdO6u0quWsjxoZWW4kGPFzmy0jqRfLGf8BScoPPvKtLaOjM0WO7Irc7l8ppVrZ8e5ecdKrr/uaGb/WZnW3pmkLFKvDrL0+WOZ/8d+6uzO+ngmK4NujvR2Z6W1P61yKsNx02OGJk2OGG8HbDGbuWoxC+1OlrpF9k6Nc/XsZq7Yt5KFG6u0r9+RYs9iMjuVjMZpLW+ktfN4yvZaBqNW1kftHB92c6xfZnE0m+UsZqOcSb9sp6jaKTI4C7vUJidEK4oyZdFJp5zJbJYyn+ksdlvZPZVcPj3M1Yur2XPNeqZvmk559c5kaS4pi5Qb/bQOLqc1t5q6Xs7GcLJE/ki/zEKvm9lqMavlbNqt6QzH6zm7R+gvfqIb4AJWFGXarZks1kvZM13mutnNXH/V0Sx+xXTKG/YlOxeSwTCtzx/OFVOHs7I5nX9a35mdyzOZqhdSFu0UKc754q6yaGeqnM9SMZXLp+s8Zmkllz1+I90n7Uuu3JO0W2kdW037k/tzw/hoDvam8/+mpjPXW0yrnEoajwg49yYLUsu0ym6mivksZCa7pspcNjXOtYur2f2YzXRv3Z3iustSL8ylGAzT3n8kC1P7c9Xmcg5szuTeqekc6HQyO1xMu5iazPGiPGtHwh503EWZsmynU85kvl7KUqeTPdPJ1TP9XLV3OXO3dNK6+Ypk367U7VbK5dVM7T6QPTmWqzemc1+vm3umWlnszaZbzKfV6qYctVOfh99NcLYUJx3p7hazmamnM9cps9hJdnVH2T23MTnCff2OFDdcnnrfnmRmOhmPUxxbTtkuM7VxMLuPr2fX+mwWO53Md8pMl+1M1TNpFVOTM4OnTIry7Kxm2Xo/d1m20y6n0q1nMlt2Mt8psqNTZc/UILt2bmTq2k7KG/YmV+9LvXNpcv/rGynmZtJOMr+6nN3Lm9nVm85ip5u5TpmZ3vTkJK/l5GzmTf9eutiIboALWJEyraKT6XSz2En2TPcye3WV8jGXpb7x+tS7diX9XsqZ6XSOb2TPP65nx5EdmWm30h5NbX+Mxzkfd1Gmm9nMdcrs6FTZtWM9nevnkxuvSXXVlUm3m+Lw4ZRVlbl7Ppfd9/Sz0JnOTKa3IyLeC8Yl6MSR7k6mMl22MtdJdnRGWVrcOsJ9zd7U112demkp9WCQYnoq5epm5u8+nF0H+lnqTGWmVWRmMJN2OXXKpxg048R7QDtpF9OZqqcy1ymz0K6za6qfucuHaV17eerrrkx9+WVJu50sL6eo6nQObWTXZ9ezc3kx8+1Wpst2uplNWbQnL+LholZux2UrnXTSSrcsMt2qM98eZ3ZmMFlSvmcx9Z5dqfftTWZmkvEo6XRSbPbS2r2Sqfn1zHeHmW9VmW610i3LdMaTs5mXZTNzpSzaKdNJp+6k2yoz3UpmW3XmOsNMLY1S7l5M9u5IvW9v6h07krJM1ier5orVjbR3rWV+tp/Z5XGmW0m3LNJJO61MdhJ84ceDIroBLnhFyrRTplMmM51RWkvtZGk+9a5dqXfsTIaD1KtrKZZmMjW9kqmyTqfMeX/SK1KmXRSZbo3TnRml2Lkr9c6l1Lt3J612UlUpdh5Oa6HMdHucbpm0U6b0ZM2l6qQXz6100inKtIvcP0cWFlIvLUyCe3EpGQ5SrK+nWJhJa67IdHucdpl0yiKttLbOhH5iCec5GH7KtNJKq0i6ZZ1ua5zWXJEszCRLi6mXliZzO0mxMJdyoZvu1Ga6ZZVOmbSKImVdpiha52S8cC6UW5FZZvI53J0y6RRVOt1xiun2ZEn57HQyM5N6ZnYS3b1+ipnpZLaT1lSdTjlOq5h8pFi7KFLUZWM7zU+8NijTSpkyrbLYut/JnC47SWa6qWemU58Yc1lOFqfPribTnRTTrbTaVTplvTXu4sStbe+IsGLtVB4NgItAlTpVnQyrMvWgSgbDpN+bvCjv9yb/HowyroqM6qS6AA4Q16lSJxnXRepRmXowSgbDFL3JuNPvJ4NR6nGdcVVkXGfrGnCJOuljf+rJp/Fuz5FqWCaDcYrBMBkMJnN7MEiGw2Q83p4ndT1Z/1Fncltn9aOEzsCJ30WT+VqkHiUZjE4a72Ay/vEo9XCcalxMfgdk8glCVc7teOFcqusT/50n1bhIPaqS0XjypDweTXY2j0aTv49Hyaja2r7YfrtFVSd1UW3P8ebHPJnTk/suJoMfj1OMxsmJsVbV1njHk49MGdWpq2L7tcbkd5Ln79NxpBvgAlanSlUP0y9GWR8lK/1uBvvHaX3+cMqZ6dSra0m/n+LzBzK6bz0rq7NZH5XpjUcZpX9exz7OML1xndVRmdWVqSzeu5LW5+5LRuMU3U6Ko8dS33M4g/11lgfdbI6TYSaf73nuP+YMzo0Tc3qUYXrjcTZGnayO2llfm8rioc209x9JMT2VYn09GQ5T7D+U+sBKekfKrA472RgX6Y3HGWSYcT1MnXFyDl6c13WVKsMMt+b1xrjISr+b/uGkc3A5xb0HkqpO0WolK6vJ/iMZH+plbX02a6NWNkZJvx5nXAxTVcNTdkDAxWny+dnjepRxhhlmnGFVp1+V2Ri10+t1Uq320lreSHFsOel0kl4/GY8m/z62mnq5n+FaK5vDdjbHRQbjZFhXGWWcOuPJc+FZniuTHeLVZMzFMMOqzqBKeuMim8N2BmtlppZ7KY6vpFiYm1ypLCerbo4tp15ez3h1lH5/Nr1xK71xkWFVZbR1m3V94nPFzfGTiW6AC1hdVxnVg2wUGzk2WMznN2Zy1d3zuWLmcDrHN1LsmE36w0lwfzT5/MpCDvaLrI2GGRYbqermPrv3dKpqlH69lpXhMAd60/nc8cUs/sOBLIw+m/KKw0mnner4Rgb/tJEDn1vMPZvTOdavs1asZlz148maS1ZdpapH6Wct6/Ugy4OpHOy3c9/yfBb+cTML3f0pVzdTLM0mw1HqAyvpf2Ithw4s5b7N6RzpJ6vDcTbKtYxGG6mqUcNHuycvnqt6mGG1mY1yLWvDHTnS72R/bzpX7p9P9+Mr6VZ3pzh8PGm1Ui9vZHz3ctY+Vea+lYUcGrSyPKizVvfTz1qqerT1wt+ONS5u9+9E66efYTZHddZGyfKwnWPrM9mxf3PrLOVlis3eZEn5eJQsr6X+/JEM7+ll5fhcjg+6WR0VWR/V6VXjDIrNjMb9yfzOifA+S/OlrlJVo4zqXoZFP71qnM1RO6ujMkf7U1k9Np2ZezbTWTicIkmxvDp5T/dmLzl4LNW9K+kdKHJ0fTbLw1bWR8nmqMpmehnV/VT18Jwdpb+YiG6AC1adpMpwtJ5jnYO5b2N3/t9ad/KxYL2p7P7HjUzPHs9oVGZldTafW1nMP6zM5HNrVQ4Vx9Mbr6SqBuf8hW2dOuN6kP54OYfK5XxufSofWZ5L/Y+X56rDq1mYX0/ZqtPb6OTw6p58amU+n1xr5d6NXpaLQxmNe+d8ySycC3XqVPUoo/FmetVyjhcrOdKbzT0b7Sy159K5u8qV6ytZuPtwWnNF6mGd3rEyhw4s5VPHl3L3ZjsHN6scHW9mvTyWYbW5tWOt2flS11XG1SCD8VrWW8dzdLw7h3rtfHajnYWjO1J/rMjuI2uZ2rWaolVktFZn9eBUPn9kR/5pbTb3bBQ50htmtVjOoFrLaLyZ+jztEISzY3I0t6pHGVeD9Ku1bBRrWR3N5/iglUODVvavz2b+3kH2lWuZ2jiY1u6VZLYzWVK+3M/wnl6Of3Yq9y4v5ECvm2ODIqvDcdbqfnrF2tb8Hp7F58M6qavtHQXDajMb7ZWsVbuyMuzkyKCVg/1udi4vZPozwyxmOe3lXoql6aQsko1hxoc20/tclcP7F7J/Y2Z7h9rqaJResZF+vZpxNZjsWHO0+xSi+xL2lNs/nun6s6fd5i/X33aORnNm3nHsjfn07tXTbrNcHDtHo4Hzr65HGY7XszK+N58pd6dz9PKsDrv5p/Vd2dHZmemyyrgusjoqc6hf5O61Kv+4uZqDxafTGxxLVQ9y7p/0qlRVL5vDYzk489l8cn0myUIO9uezZ2Uu860qraJOrypzfFjkvo0in1kb5u7cl7XR/snne3pBziWpSn3iRfp4NSudQzkwWsj0+nzaZTu9aikHNmey61A/0+1xxlWR1WEn921O5+7Ndj6zVudAb5Aj5eFsjI9kONqYzPGm345Rj1LVgwxHG1lvH86RYkfmNrrpllNJprM6aueylfnMd4cpizq9UStH+lO5Z7ObuzfKfH59nAPj1RwvD6Q/Xsm4GnhBziXgxHzuZ1CtZa19PMeqpSz22ltn659Ku1zKYNTK7uPrmV5cT2uqTjVOhmutrByfy73LC/ns+mw+v9nK4V6dY4PJzqlevZxRtbk1V0Znbed5nXr7d9BovJleayUrxWqO9qcy3ykz126lW84mSS7bXMvC/s1053opymTcL7Kx3MnR44v53Op8PrfZzf7NIkf64yzXG1kvj2cwWj9pp5r5fTLRfQn76+P/v/M9hEesN7g3fz74/873MOCCUdfjjMYbWevfl3unygyKzRw4ti87l2cy026lVRRJkt54lLXRMIeK4zlYfDrLg7szGC2nrgY59++NrlNXg/RHyznW/3TqbpW1jevyuY0dWWh10m0VSYqM6zqbo1GOVZs5WN6XI9Wnsz44mHG1sfU+VbgE1aOMq14Gw9WslvuzvzWVsnd1htVclget3Ds1naXOVLplnXFdZH1c5Gg/ObhZ5b5eP/flYI7X96Y3Op7heP0crGapt1+oj6rNbA6P5lj33rTSSdb2pDfu5ki/k53dTubadcrU6VVFVoZFjvSTAxuj3Ddcz8Hy3qyPDmYwWt3aUWDHGhe5ehLdo3Evw/F61ssjWS4Wc7jfTbecTlmUGdYzWR62s2d9NvOHhumU41SZvHf6aH8qB/qd3LvZyr0bdfZvDnOkXs1yeSi90UpG482TdpyfvYCt61GqapDRuJf+eDkr7UOZqWYytdFKq2inqjvpjedztN/NjuVBZjqjlKkzGLeyOuzkcL+b/b1O7tksct9GlUODzRwrj2S9OpzheP3+HQW15/GTiW6AC9o4ddVLf3g8dV2l31nJkfaOTNULk8/h3voQinGGGRRr6Y1X0h8upz9aPq/xWmec8XgtG4NkXA2y0Tmc+1pL6VQzaVWdrW2qjOp++llNb3A8vcGxDMdrqap+nESNS9MkYKtqkMF4NeWwnSJl6laVzcFlOT6cz1K7k9l2mU45OZNxbzzOymCcY9VmjpSHc7y+N+ujQ+kPlzOuelvx2vQRpSp1Nci46GUwWs1aeSBpT37vrG/szpH+dBY6rUy3yhQpMqzqbIzGWR4Nc6ReydFyf1bH+7M5PJbheOOkHQXmOReryVy+fxXIevrlSpbbB9IqOil6uzOup7I5KnN80M2Obmd7lVedIpvjIqujIkf6RQ736hzsDXNovJaj5YFsjI9kMFrdeqvV6CyfSK1K6iJVPcio2kx/vJr18nCOlFMpx0WyvpBR1c76qJXDg8kOwOlycv/DusjGuMzxwdaOwN44B3v9HC6OZKU+mP54dXtHwf1vHzHHTxDdABe4OqPU1Wb6w1GG47VsFIdSFu2TPgsz2+8tq6rB5AmvGmwF9/l6wqtT18OMRyvZrHrpj5ZTFu3tcZ8Y8/a468HkhXg9TBzl5pJWTV6oV730h8uTS+phNlvLOVbszOxoMTPDmbTSSp0qgwyzUa5lvTyWjfGR9EbHJy/Iq43U9bmK1xM7C3oZJslgsgpn2J4sKT1S7ch0bzZTmUqSjDPOZrGZjWIlGzmWjdGR9IfLk6Pc52xHATRtErCTVSCTHVIbRTtFu0zqZNjbkfXRVI4PWpnrlJlpTT7jvk4yGCfrozorg3GODydHuI+WB7JaHZysYhmtb8fr2Z3j9+8sGFe9DEft9Ip2inYrKZNqXKe3PpfVYSdHumXm2mU6W5/RPa6TzXGyOqiyMhzn2KiXI8XRHMt9kx0Fw9WMqpN3FJjjJxPdABeFyUeHjMfDVNl4yK3qrZOvXRh7l+utHQbjVFU/RYqH3GriQhk3NKme7BCrBhkn6Q2rjKt++q3VbLSOpFPOpFV0UqYz+VifepjRaCPDajPD0cbW+zx7qare1vLNc7WTajz5POGt8K62ltX2Wsez2ppPq5iaHOFLmSrDjKp+htVmRtVmhqP1jE6M+bzvEISz5aSj3dUgw/H6/TuVW1UG5WbWRjtyfDST2bKTblmmvfWWsH5VpVdNTpq2WixntTya9erw9k61yRwfbM3xsx2vW0e7q0FG2bh/zO1xRmU/m/XerPUXMjvoZrZsp1NOxjyq6wyqKhvVMCtZz2p5LGv1kWyMJzvVtt/uYo4/KNENcNGok4wvwqexi3Xc0JSTwruutl6wb6RfnFgR0klRlJPVIBmnqkanrAjJ9tGvc70qZCu8xxtbJ2PqpRx10yuOpSzbKdLaHndVD09afTOaHJXfDgi/DbhUVKnrJBlktDUdT3zUZ7+1mvVyKTNZzFQ9k864m6LeektYMU6vWE+vWEuvXs5gtJb+eHWyg2q8cdJ5D5qYL5PfP0WdraXxa/efib3VT7+1ltVyMdOZz1Q1k9a4NblWUWWYQTbL9fSzll61nP54NYPh6tZ7uXsnrb7hCxV1XZ/RI1MU+hwuVI/2o1fMa7hwfTEfqWRuX0yKJGWKorX1Z7l1+Yk/q/vP9F2f+Gzr8x2uW2NOkRRl7h/3/WNOsj1usX0qz9mXmvvncFF0U5bdtMvptFszabdmtlav3L8SJMkDVoOMxpMzlZ84wj3Zsdb00eIiRVpJ0U5ZdtMqp9MuZ9JubY29mEq7nEqZ+8/DMq6HGVUbGdWTs5+Pxpv3r2LZPnnal+Zcf7h5LbrhEuAJHC49ovtLzYO//eJUF+IL2Yt13OeP5+xL0ck7oiYRWxTttIpuynKyeqUs7///r66rrRUsw63PtT5xxu9zHa4n7zBopyy2xl12t8f8wPOwDLc/i7vaOlP5+Vt9c+F4uHlt9gIAnHcXa5herOOGs+nE26iKpK5TjUdJ0U5V9FJUk08pyPYqlmytWqm2QnZrGfnW20nuv71zNe6t1TR1lXExSlJmXG2drPVBxl1tjXdyndEFsvrmwie6AQAAvminxnddJ0UGk2+dHK/JSW8XSc5vtNZb/ztO6nGScium86Bjnmx7IrTvvz6nJ7oBAADOmvuXWtcn3oJRP9jS6wspWE+M5cSJT4uHGPPJ23KmRDcAAEAjLtZAvVjHfWEqH34TAAAA4NEQ3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0JCiruv6fA8CAAAALkWOdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ/7/kOTOlrvbOOIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x600 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2,4,figsize=(10,6))\n",
    "for i in range(4):\n",
    "    ax[0,i].imshow(input[i].cpu().detach().numpy(), cmap='magma')\n",
    "    ax[1,i].imshow(ngon_output[i].cpu().detach().numpy(), cmap='magma')\n",
    "    ax[0,i].set_title(f\"Distance {distances[i]} cm\")\n",
    "[a.axis('off') for a in ax.ravel()]\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chaining Operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to chain operators in `spectpsftoolbox`. For example, in SPECT acquisition, the component above models the PSF component from the collimator structure only, but there is also an PSF component resulting from uncertainty in position when detected in the scintillator. The PSF is best modeled as\n",
    "\n",
    "$$PSF_{tot}[x] = PSF_{scint}[PSF_{coll}[x]] $$\n",
    "\n",
    "$PSF_{scint}$ is best modeled using a Gaussian function with width equal to the intrinsic resolution of the scintillator crystal. We can define the operator as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/anaconda/envs/pytomo_install_test/lib/python3.11/site-packages/torch/functional.py:504: UserWarning: torch.meshgrid: in an upcoming release, it will be required to pass the indexing argument. (Triggered internally at ../aten/src/ATen/native/TensorShape.cpp:3526.)\n",
      "  return _VF.meshgrid(tensors, **kwargs)  # type: ignore[attr-defined]\n",
      "/data/anaconda/envs/pytomo_install_test/lib/python3.11/site-packages/torchquad/integration/utils.py:255: UserWarning: DEPRECATION WARNING: In future versions of torchquad, an array-like object will be returned.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "intrinsic_sigma = 0.1614 # typical for NaI 140keV detection\n",
    "gauss_amplitude_fn = lambda a, bs: torch.ones_like(a)\n",
    "gauss_sigma_fn = lambda a, bs: bs[0]*torch.ones_like(a)\n",
    "gauss_amplitude_params = torch.tensor([1.], requires_grad=True, dtype=torch.float32, device=device)\n",
    "gauss_sigma_params = torch.tensor([intrinsic_sigma], requires_grad=True, device=device, dtype=torch.float32)\n",
    "scint_operator = GaussianOperator(\n",
    "    gauss_amplitude_fn,\n",
    "    gauss_sigma_fn,\n",
    "    gauss_amplitude_params,\n",
    "    gauss_sigma_params,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could use this operator directly on the input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_tot = scint_operator(ngon_output, xv_k, yv_k, distances, normalize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OR we can chain together the operators to get a single operator. We can use the `*` operator to make an operator that first applies the operator on the right, then the operator on the left:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "spect_psf_operator = scint_operator * ngon_operator\n",
    "# Now apply directly to inpuy\n",
    "output_tot = spect_psf_operator(input, xv_k, yv_k, distances, normalize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAIfCAYAAABtmb+NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABA50lEQVR4nO39eZDteUHf/78+n7P03n33WZk7BAaYUWD4RUHjMhpQhAGDGI1bCYOIUShiKEsjikuqWKxIBctgod+ETQqMYNAoMolEcClQUdCICAEEBmbm7rf3Pn2Wz+f3x+nuuXeWO3eG+7kbj0dVM3NPf8457+7hfT/n+fm8z+cUdV3XAQAAAM658kIPAAAAAC5XohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIV9W0f3mN785RVHsfE1OTubKK6/MN3/zN+fVr351jhw5cp/7/MIv/EKKonhIz7O+vp5f+IVfyAc+8IFzNPIL761vfWu+53u+J4997GNTlmWuv/76h/wYv/qrv5rHPe5xmZiYyCMf+cj84i/+YgaDwbkfLJctc/jhufvuu/OzP/uz+dqv/drs27cv8/Pz+ef//J/nN37jNzIajU7b9gMf+MBpv+NTv/7iL/7irJ7vyJEjef7zn599+/Zleno6X/u1X5v/83/+TxM/GpcB8/rcOHz4cPbu3ZuiKPKud73rPt9fXV3Nj//4j+fqq6/O5ORkbr755vzWb/3WWT++ec1DYV4/fNdff/397oP/7b/9t6dtZ399aWlf6AFcCG9605vyuMc9LoPBIEeOHMmf//mf55d+6Zfyy7/8y/nv//2/52lPe9rOti984Qvzbd/2bQ/p8dfX1/OLv/iLSZJv+qZvOpdDv2B+8zd/M4cOHcqTn/zkVFX1kGP5la98ZV7xilfkP/yH/5Bv/dZvzYc//OH87M/+bO688878xm/8RkOj5nJlDj80f/M3f5O3vvWt+cEf/MG84hWvSKfTyXvf+9786I/+aP7iL/4ib3zjG+9zn1e96lX55m/+5tNu+8qv/MoHfa7Nzc089alPzeLiYn7lV34lBw4cyOtf//p827d9W973vvfllltuOWc/F5cX8/pL8+IXvziTk5MP+P3nPve5+fCHP5zXvOY1ecxjHpO3v/3t+d7v/d5UVZXv+77vO+Njm9c8XOb1w/N1X/d1+eVf/uXTbrviiivud1v760tE/WXkTW96U52k/vCHP3yf733+85+vH/GIR9Rzc3P1oUOHvqTnOXr0aJ2k/vmf//kv6XEuJqPRaOffb7311vrgwYNnfd9jx47Vk5OT9Yte9KLTbn/lK19ZF0VR/8M//MO5GiaXOXP44Tlx4kTd7/fvc/uLX/ziOkl9xx137Nz2/ve/v05Sv/Od73xYz/X617++TlJ/8IMf3LltMBjUN910U/3kJz/5YT0mlzfz+kv3rne9q56dna3f8pa33O/8fc973lMnqd/+9refdvu3fMu31FdffXU9HA7P+PjmNQ+Vef3wHTx4sL711lsfdDv760vLl9Xy8jO57rrr8trXvjYrKyv59V//9Z3b72+pyx//8R/nm77pm7J3795MTU3luuuuy3d+53dmfX09n/vc57J///4kyS/+4i/uLPN4/vOfnyT59Kc/ndtuuy033HBDpqenc8011+TZz352/v7v//6059heMvKOd7wjP/MzP5Orr7468/PzedrTnpZPfvKT9xn/7bffnqc+9alZWFjI9PR0brzxxrz61a8+bZu//uu/zrd/+7dnz549mZyczJOe9KT89m//9ln9fsry4f9f5fbbb0+v18ttt9122u233XZb6rrO7/7u7z7oY9x555150YtelEc84hHpdru5+uqr86//9b/O4cOHk9zz+3r729+en/qpn8pVV12V2dnZPPvZz87hw4ezsrKSF73oRdm3b1/27duX2267Laurqw/7Z+LiYw4/sN27d6fT6dzn9ic/+clJki9+8YsP+hhn693vfnce+9jH5mu/9mt3bmu32/mBH/iB/NVf/VXuvPPOB32MB/tdPP/5z8/s7Gw+8YlP5OlPf3pmZmZy1VVX5TWveU2S5C/+4i/y9V//9ZmZmcljHvOYvOUtbzlnPx/nl3n94E6cOJEXv/jFeeUrX5nrrrvufrd597vfndnZ2XzXd33Xabffdtttueuuu/KXf/mXZ3wO85pzyby+OJjX55foPsUzn/nMtFqt/Omf/ukDbvO5z30ut956a7rdbt74xjfm9ttvz2te85rMzMyk3+/nqquuyu23354k+aEf+qF86EMfyoc+9KG84hWvSJLcdddd2bt3b17zmtfk9ttvz+tf//q02+085SlPud+J/fKXvzyf//zn81//63/Nb/zGb+RTn/pUnv3sZ5/2Psz/9t/+W575zGemqqq84Q1vyO///u/npS996WkvpN///vfn677u67K4uJg3vOEN+b3f+73cfPPN+Tf/5t/kzW9+8zn6Dd6/j33sY0mSxz/+8afdftVVV2Xfvn07338gd955Z776q7867373u/Oyl70s733ve/O6170uCwsLOXny5GnbvvzlL8+RI0fy5je/Oa997WvzgQ98IN/7vd+b7/zO78zCwkLe8Y535Cd/8ifzm7/5m3n5y19+bn9QLjhz+KH54z/+47Tb7TzmMY+5z/de/OIXp91uZ35+Pk9/+tPz53/+52f1mB/72MfyhCc84T63b9/2D//wD2e8/9n8LpJkMBjkuc99bm699db83u/9Xp7xjGfkp3/6p/Pyl788z3ve8/KCF7xg5wXF85///PzN3/zNWY2fi495fWYvfelL88hHPjIveclLHnCbj33sY7nxxhvTbp/+rsLteflg+2HzmnPNvD6zP/3TP83c3Fw6nU5uuummvPa1r73PNVi22V9fIi70qfbz6UxLXbZdccUV9Y033rjz55//+Z+vT/01vetd76qT1H/7t3/7gI/xUJa6DIfDut/v1zfccEP97//9v9+5fXvJyDOf+czTtv/t3/7tOkn9oQ99qK7rul5ZWann5+frr//6r6+rqnrA53nc4x5XP+lJT6oHg8Fptz/rWc+qr7rqqtOWjz+Yh7q8/Id/+IfriYmJ+/3eYx7zmPpbv/Vbz3j/F7zgBXWn06k//vGPP+A227+vZz/72afd/uM//uN1kvqlL33pabc/5znPqffs2XOWPwEXC3P43Mzhuq7r//W//lddluVpY67ruv7IRz5S/7t/9+/qd7/73fWf/umf1m984xvrG2+8sW61WvXtt9/+oI/b6XTqH/mRH7nP7R/84Afvd3nrqc72d/G85z2vTlL/zu/8zs5tg8Gg3r9/f52k/shHPrJz+/Hjx+tWq1W/7GUve9Cxc2GY1w9/Xv/BH/xB3el06r//+78/bXz3Xm56ww031E9/+tPvc/+77rqrTlK/6lWvOuPzmNc8VOb1w5/XP/ZjP1a/8Y1vrP/kT/6k/t3f/d36+7//++sk9Q/8wA+ctp399aXFme57qev6jN+/+eab0+1286IXvShvectb8k//9E8P6fGHw2Fe9apX5aabbkq320273U63282nPvWp/OM//uN9tv/2b//20/68ffTp85//fJLkgx/8YJaXl/NjP/ZjD3jFx09/+tP5xCc+ke///u/fGcP21zOf+czcfffd93vE71w609UoH+xKle9973vzzd/8zbnxxhsf9Hme9axnnfbn7fvceuut97n9xIkTlphfhszhB/eRj3wk3/3d352v+Zqvuc+SuCc96Ul53etel+c85zn5hm/4htx222354Ac/mKuuuio/+ZM/eVaP/3Dn+9n8Lk59nGc+85k7f26323n0ox+dq666Kk960pN2bt+zZ08OHDiw8/vm0mRe39fS0lJ+5Ed+JD/1Uz91VhdN+lL2w1/K/c1rHoh5ff9e//rX57bbbss3fuM35l/9q3+Vt73tbXnJS16St73tbfnoRz+6s5399aVFdJ9ibW0tx48fz9VXX/2A2zzqUY/K+973vhw4cCAvfvGL86hHPSqPetSj8iu/8itn9Rwve9nL8opXvCLPec5z8vu///v5y7/8y3z4wx/OE5/4xGxsbNxn+717957254mJiSTZ2fbo0aNJkmuvvfYBn3P7fc8/8RM/kU6nc9rXj/3YjyVJjh07dlbjfzj27t2bXq+X9fX1+3zvxIkT2bNnzxnvf/To0TP+fKe692N1u90z3t7r9c7qcbk0mMMP7qMf/Wi+5Vu+JTfccEP+8A//cGc8Z7Jr164861nPyv/9v//3fn/GU+3duzfHjx+/z+0nTpxIct+5eKqz+V1sm56evs+Vmrvd7v0+frfbNdcvYeb1/fuZn/mZdDqdvOQlL8ni4mIWFxd3DiSvr69ncXFxJ2q+lHn5pd7fvOb+mNcPzQ/8wA8kyYN+FJj99cXry/Ijwx7Ie97znoxGowf9yIFv+IZvyDd8wzdkNBrlr//6r/Orv/qr+fEf//FcccUV+Z7v+Z4z3vdtb3tbfvAHfzCvetWrTrv92LFj2bVr10Me8/YFJM50IaR9+/YlSX76p386z33uc+93m8c+9rEP+bnP1vZ7uf/+7/8+T3nKU3ZuP3ToUI4dO/agR+j3799/Ti/0xOXLHD6zj370o3na056WgwcP5n//7/+dhYWFsx7n9ov3Bzui/fjHP/4+F6lJsnPbmeb72fwu+PJjXt+/j33sY/nc5z6XK6+88j7fe97znpckOXnyZHbt2pXHP/7xecc73pHhcHja+7rPZl4m5jXnnnn90Gzvg8/mwsb21xcnZ7q33HHHHfmJn/iJLCws5Ed+5EfO6j6tVitPecpT8vrXvz7JeMlmct8jY6cqiuI+Z5be8573nNUVAu/Pv/gX/yILCwt5wxve8IDLdB772MfmhhtuyN/93d/lq77qq+73a25u7mE9/9n4tm/7tkxOTt7n4hFvfvObUxRFnvOc55zx/s94xjPy/ve/v/El8FzazOEzz+G//du/zdOe9rRce+21+aM/+qPs3r37rMd48uTJ/MEf/EFuvvnmM34OcJJ8x3d8Rz7xiU+cdjXk4XCYt73tbXnKU55yxrMaZ/O74MuLef3A8/p1r3td3v/+95/29Z//839OMr4K9Pvf//7Mzs4mGc/L1dXV/M7v/M5pj/GWt7wlV1999WkHxO+Pec25ZF4/9Nfcb33rW5MkX/M1X3PG7eyvL15flme6P/axj+28v+LIkSP5sz/7s7zpTW9Kq9XKu9/97p2jN/fnDW94Q/74j/84t956a6677rr0er288Y1vTJI87WlPS5LMzc3l4MGD+b3f+7089alPzZ49e7Jv375cf/31edaznpU3v/nNedzjHpcnPOEJ+Zu/+Zv8p//0n856+fS9zc7O5rWvfW1e+MIX5mlPe1p++Id/OFdccUU+/elP5+/+7u/yX/7Lf0mS/Pqv/3qe8Yxn5OlPf3qe//zn55prrsmJEyfyj//4j/nIRz6Sd77znWd8no9//OP5+Mc/nmR8hnp9fT3vete7kiQ33XRTbrrppiTJn/zJn+SpT31qfu7nfi4/93M/l2S8POVnf/Zn84pXvCJ79uzJt37rt+bDH/5wfuEXfiEvfOELd+77QP7jf/yPee9735tv/MZvzMtf/vI8/vGPz+LiYm6//fa87GUvy+Me97iH9bvj0mUOP7Q5/MlPfnLnZ3vlK1+ZT33qU/nUpz618/1HPepRO7+z7/u+78t1112Xr/qqr8q+ffvyqU99Kq997Wtz+PDh+xw4+6Ef+qG85S1vyWc+85kcPHgwSfKCF7wgr3/96/Nd3/Vdec1rXpMDBw7k137t1/LJT34y73vf+87J74LLk3n90Ob1zTff/IDf+4qv+IrTziA+4xnPyLd8y7fkR3/0R7O8vJxHP/rRecc73pHbb789b3vb29JqtXa2Na85l8zrhzav3/72t+d//I//kVtvvTUHDx7M4uJi3vnOd+a3fuu38vznPz9PfOITd7a1v77EXICLt10w21dS3P7qdrv1gQMH6ltuuaV+1ateVR85cuQ+97n3lRQ/9KEP1d/xHd9RHzx4sJ6YmKj37t1b33LLLfX//J//87T7ve9976uf9KQn1RMTE3WS+nnPe15d13V98uTJ+od+6IfqAwcO1NPT0/XXf/3X13/2Z39W33LLLfUtt9yyc/8HugLpZz/72TpJ/aY3vem02//wD/+wvuWWW+qZmZl6enq6vummm+pf+qVfOm2bv/u7v6u/+7u/uz5w4EDd6XTqK6+8sv6X//Jf1m94wxse9He3/Xu4v69Trxi5Pe77u4rkr/zKr9SPecxj6m63W1933XX1z//8z9f9fv9Bn7uu6/oLX/hC/YIXvKC+8sor606nU1999dX1d3/3d9eHDx8+4+/rga6euf3zHD169Kyen4uDOfzw5vC9f2/3/jp1LK9+9avrm2++uV5YWKhbrVa9f//++ju+4zvqv/qrv7rP425flfSzn/3sabcfOnSo/sEf/MF6z5499eTkZP01X/M19R/90R+dcYwP5XfxvOc9r56ZmbnP/W655Zb6K77iK+5z+8GDB+tbb731rJ+f88u8fvj75nt7oPHV9fhqwy996UvrK6+8su52u/UTnvCE+h3veMd9tjOvORfM64c3rz/0oQ/VT33qU3de705PT9df/dVfXf/ar/3afa56bn99aSnq2poAAAAAaIL3dAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBD2me7YVGc9abAeVbXw4d1P/MaLl4Pd14n5jZczOyz4fLzYPPamW4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGhI+0IPgOaM3vqS5CseecZtvvqp/5iPLP7X8zSiB/fcPT+dd/7RvjNuM3rPR9P9ubedpxHBxcW8hsuTuQ2XH/OabaL7cvYVj0z1xCeecZOp+rPnaTBnZ6ZTPuiY25+98zyNBi5C5jVcnsxtuPyY12yxvBwAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAa4nO6L2NPeeonMll//ozb/OXa28/TaM7OO0++KZ/du3LGbZaKk+dpNHDxMa/h8mRuw+XHvGZbUdd1fVYbFvocLlZ1PXxY9zOv4eL1cOd1Ym7Dxcw+Gy4/DzavLS8HAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGiI6AYAAICGiG4AAABoiOgGAACAhohuAAAAaIjoBgAAgIaIbgAAAGhIUdd1faEHAQAAAJcjZ7oBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABoiugEAAKAhohsAAAAaIroBAACgIaIbAAAAGiK6AQAAoCGiGwAAABrSPtsNi+KsNwXOs7oePqz7mddw8Xq48zoxt+FiZp8Nl58Hm9fOdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBD2hd6AAAAAJeXYuuf5Sm3FKlTn7JNtfXPU2+70C7VcV/cRDcAAMA5USQpU6RIinL871v/3P5uUqWuqyRVUldbQXshQ/ae0B6Pu7015vFt92zxQOMW3w9GdAMAAHxJTo3tdoqinfK0f95z5riuq9SpUtXD1Kd8jUN2lPMbsUWKtJKiTLE13p0x7xwwuGfc22Ou6mGy9SW+H5zoBgC4oMYv1sf/VtzvFhf+TNi93evM2P24ZzmqF+Nc7rbDdStYy27a5WRa5URaZTdl2UlZtFNuzfNRPUxVDzKq+hlVmxlV/VRVP1XdH0dsPUrz82brIEHR2orsbsqym1bZ3Rl3q5zYGXO1dZCgqsbjHo42MqrH474nvs/3AYNLh+gGALgg7ucMU8rxn7de6Napts6CVVtnwi70WaX7P5t37zFn60zehT2DB+fDeB4XZTdl0U2rnEynPZNOaybd1kw6xXS6xXRa6aSVdqpUGWWQYTazWa1mWK2nP1rLYLSW4aiXUdVL0k9dJ83N81ODezzmdmtyZ8zd1mw6xXTamUgrnZQpM8owowzSr9fTr1bTr9YyGK6N47saj7uoY54/ANENcEm5/zNKp7sYd3YPNu6LcczQpNZpL3jvfXapKMqdJaijanPrhW0/o6q3c2bpwi1DbadVTm6dzZtKq7U15q2lqHVd7Zy9G1WbGe6Mub91Bm90HscMTbpnTpRFN+3WdCY6C5lozWeqtTsz2Z2Zaj5T9eRWcpcZH36qspF+1ovVrLYXs9E6mY3hYvrFSgajMsNR0lx4nx7c7dZ0Oq2ZTHQWMtmaz1SxO7P1rkxXs5lKN520kiTDVBlklPViPeut1ay1TmajPJnecDH94dYZfOH9gEQ3wCXh9KVrxSlnlbadelbp4nlh29pZurZzFu8U9dZFWWpL0/iy0kpRdLbOLk2n257LRHs+E+V8JorZTGQ6rbqdqqgyyGY2s5petZTN0Ur6g5UMRmsZVetJ1T+Pc2b7bN5kWuVkuu258bhb85ks59PNdDqZSFmXGRXDbGY9m/VqNqvl9IaLGQzXMhitpqp6WyFxMfz9BF+KrbeFbB2E2g7umc7+zBYHsrvan12ZyXy7m6l2mclWkdbWLrA/qrMxnM7KcDYnq4UsFfNZak9mvRin2fjAVZK618x57lOCu9uey2R7ITPtA5mv92d3tScL5WTmOu1MtYt0W+OD5qMq6Y3qrA2nszpcyInMZ7E1nbLo7Lzvu06VqkrGk9zbSk4lugEuelsv0FvT6bbm0mlPp1NO7ezots+GDaqNDEcb6Q9XMhytp656qTO8QGO+5wX69k6905pJ+15nw6p6sLOsbjBaz2i0nrrejB01l68iRdHaWoI6m8nOnsy092ehuDK7qt2Zy1Qmy1Y65fiMWG80ylrdz2KxnOXO0ayUh7IxOJHNwVa2npfwvmf57HZYTLX3ZLZ1ILvqK7KrmstM2U23LFOmyKCusjEaZjkbWWqdyHJ5JGvl0RSDMv1htsLbC3IufdsHlbcPRM109mdXcU32VQeyvz2dvRPt7J4oMt9Jptt1WsW4RzerMqvD5MRmmbleO8f73XSKiRxvl6m2D0anyqgenuOA3ZrLW+8777RmMtleyGz7yuypr87+elf2TUxk72Qru7rJbDuZbI2fd1AlG6Myy4MyJzfbmel1MjGaSKtsJ+2M3++99VaY8UF0TiW6AS5qRYqik057PtMTBzLfvjoLuSLz9UIm0x0vVSvqDDLKanstJztHsjy6K2ubR9IbHE9dbeTCnFEqU7amM9HendmJK7PQuia7qv2ZyVQ6aaVMkWFRpVf0s9xaysnclZXBXdnoH89guCy8uUxtv+Adx+tkZ08WOtdmf/WIXFnsyt6pTnZPlJnrJN0yGdXJ+rCTpf5Ejvemc3g4l0Ot8YGruq6yOawyqqutF+VNzvPx2bxyKyymO/uyu3xErqyuzhWd6eyZbGVXt8hUKymLpF8lK4NuTm5O5mhvJoczm6OdztZZsGEGO+/xvlAHBeFcGH8MWFl2025NZbK9K7PFgeyrDuTq7myumm7lyqnkwMQouzujzHWG6RRVqiTrw3YWB+0c7bZyqNNKd30yRW93qrrKqDVI1RlfaK2uh6mq7Y/pOke2riHR3nrv+XRrX/bUV+fK7MlV091cMVXm6qk6e7uj7O4OMtkapUwyqMusDNo50W/lcLeVqXYnrbXZZHhlRuUwo/Zmqmoc3vccLLCiZZvovqy1zmKbi3EyPNi4HR3ny0eRVlqt6UxPHMiBzuNysDqYa6amcmCqzGw7mWglwyrZGCXHN+dy9/refK7cm7sn/zFVPcjmYHgBzihtHShozWZ+8hG5trgxB8t9uWq2nT0T97ww3xwly4PkWG9XvrCxP5/rzuZY/l+qqp/haJCL8+8n+BIV4xfp3fZcZtr7s796RK5r78m1M508YqbOFROj7OkO0i2rjOoiK8N2jmy2c+d6O5Nrsyl716ZuVak6g4yqzVPe393gkFOkKNrptKYz2V7IQnlNrq6uyXVTM7lutsw1U1X2dgeZa49SFnV6o1YWB+3c3Wvli+uT6azuTapk2N4cv9e77m+9KC9if86lqdiZF9vxOtXand3V/uxvT+eq6Vaum6lzcHqQa6c3sndmPdNT/XS6o9RVkY2NTk6uTeXQ2nSmWhNpFWWSiQw2dqVfbmTQWs+wtZFh1ct23J+bfeLWgYJi60BBa1fmiv3ZWy3kiqlurpkpc/30KI+Y3sxVM+vZNbORiclhirLOoN/K2vpEjq1NZX59Kt2yk6Sdwep0Nkf7s1muZtAer7iriu2/l8zxbaL7stXK6P97YYonPOqBNxkN8+RnfjZ/vfj/nb9hPYjn7vnpvOsP5pLWA/9fc/A//zYTr3xHTGK+LGy90J1vX52D1cF85a6ZPG6+yvXTG9k32ctUZ5hBVWZ5s5svrk/l/6120zlxZfrVRnqdxQxGqxmd94At0yonM9ndk/3FI/PozoF85e5WHjs3yNVTvcx3N9Mq66z1Ozm+OZHPrk1kfnk29dI/S6+znM3BUkbV+tb70uFycs8Zpon2fBaKK3NlsSvXznTyuPkqj57t5dq51exeWE93aphqUGZtdSJ3L81moT2TdtnOqJrJRv+KbLSWstleybDqZVT3GwzYYudAQac1Piu2t7oiV01M55FzZR4zO8yj51ZzxcJqZmb7KdpVNtc6Obk8nS+szGa6NZmq7qS/sjurxf5stpbSL1ZSFe2L6NoT8DAUZcqinVY5kYnWXGayO7uL2eyfbOeqqeTg9CA3zC/nqgPLmblqmPbeTorJduphlfmVXhbu2sjsoX7KYiGjejK9UZm1wWRWh7uyVi5ks7WS1nAtVbF9HYRzMOStAwWt7aXl5Xzmq13Z05nI/qky105VuX6ml3+2ZzF7r15P98oi5VwnRVmkWh9m4cRyFu7ayOTRhSSz6VfdbAw7WVubzUqxK71yKf1yZXw186JM7Md3iO7LWPGER2X0z/9/D7zBaJjp+vD5G9BZmO+2Mvqqr0rK8gG36Xz2rpy7I35wMRvvHLvtuSzkilwzNZUb56t81d6lXP+I45m5tk45107drzI4Mso1n5vNxJF9WRl0c/TkFTne3pX14miqrJ/36xuXZTeTrV05UO/PI+daeeLCZp5wxbHsu34tnQPtFO0io6Vh1r9YZu8X9qbKQk5uTudQ/8ost7+Y/nAxdd0/j6OG5u2cGWtNZaKcz65qd/ZNjc9wP3q2lxsPHM/eR22k84jJFDOzybDK/JGNzP3TRjp3VOlVC1nqt7I4mMvJYnfWW8ezWSxllAfeZ54bZVpbZ8ami13Zk5lcOdXKddOj3LhrOdcfPJHpf1amtX8qKYvUK/0sfPFkZj+zmTp7szaaynK/m+O9vVluHUq7NZXBaLXhMUOTxmeMdwK2mM5MNZ+5dicL3SL7J0a5dnojVx1YztwNVdrX70qxbz6ZnkiGo7SW1tPavZiyvZr+sJW1YTuLg25ObpaZH05nKfNZL6eyWbZTVO0U6Z+D/fj4gmhFUaYsOumUU5nOQuYzlYVuK3snkismB7l2fiX7HrGWycdMprx2d7Iwk5RFyvXNtI4spTWzkrpeyvpgvET++GaZuV4309V8VsrptFuTGYzWcm7P0F/6RDfARawoyrTKiczXC9k3Webg9Eauv+ZE5r96MuUjDyS755L+IK0vHstVE8eyvDGZf1rbnYWlyUzUc+PPz70AyqKdyXI+u8qJXDlZ51ELy7niK9bTfeKB5Op9SbuV1smVtD91KI8cnciR3mT+38Rk5nq70ionMt5RW5bG5Wf7c60nitnMZSp7JspcMTHKdfMr2fuojXRv3pvi4BWp52ZS9AdpHzqeuYlDuWZjKYc3pnJocjKHNzqZHezKYjGRcuvTDM7VmbD7H3OZsmynU05ltl7IQqeTfZPJtVObuWb/UmZu6qR141XJgT2p262USyuZ2Hs4+3Iy165P5u5eN3dOtDLfm063mE1Zjn8HdQoznEtWccqZ7m4xnal6MjOdMvOdZE93mL0z6+Mz3NfvSvHIK1Mf2JdMTSajUYqTSynbZSbWj2Tv4lr2rE1nvtPJbKfMZNnORD2VVjExvmBqyqQoz81qlq33c5dlO+1yIt16KpNlOzOdIrs6VfZN9LNn93omruukfOT+5NoDqXcvjJ9/bT3FzFTaSWZXlrJveT17epOZ73Qz0ykz1ZscX+S1HF/ktem/ly41ohvgIlakTKvoZDLdzHeSfZO9TF9bpXzUFalvuD71nj3JZi/l1GQ6i+vZ95m17Dq+K1PtVtrDiZ2P8Ti/gx7vbDuZzHS7zK5OlT271tK5fja54RGprrk66XZTHDuWsqoyc+cXsvfOzcx1JjOR7k5EwOWoyPhAWicTmSxbmekkuzrDLMxvneF+xP7UB69NvbCQut9PMTmRcmUjs3ccy57Dm5lrT2SqVWSiP5F2OX5RnkbPdG+/B7STdjGZiXoiM50yc+06eyY2M3PlIK3rrkx98OrUV16RtNvJ0lKKqk7n6Hr2fH4tu5fmM9tuZbJsp5vp8cFAc5xLXrkTl6100kkr3bLIdLvObHuU6an+eEn5vvnU+/akPrA/mZpKRsOk00mx0Utr73ImZtcy2x1ktlVlstVKtyzTGXXTKjopy2bmSlm0U6aTTt1Jt1VmspVMt+rMdvuZWBim3Duf7N+V+sD+1Lt2jVegrq0lSYqV9bT3rGZmqp/p1iiTraRbFumknVY6431446tvLj2iG+AiV6RMO2U6ZTLVGaa10E4WZlPv2ZN61+5k0E+9sppiYSoTk8uZKOt0yux8NNcFHXdRZLI1SndqmGL3ntS7F1Lv3Tu+bkNVpdh9LK25MpPtUbpl0k6Z0s6ay9XWfNx5kV6U6Za5Z47MzaVemBsH9/xCMuinWFtLMTeV1kyRyfYo7TLplEVaaaVM55Q53vwyziJlWmmlVSTdss5Ee5TWTJHMTSUL86kXFnauyVIszKac66Y7sZFuWaVTJq2iSFmXKYqzudArXBrGH5Q3/t9WmbSKpFNU6XRHKSbb4yXl05PJ1FTqqelxdPc2U0xNJtOdtCbqdMpRWsX4I8XaRZGibm7/vR3EZVopU6ZVFlvPW6dTVik7Saa6qacmU2+PuSzHi9OnV5LJTorJVlrtKp2y3hp3sf1o4+cotpeWs81vA+ASUKVOVSeDqkzdr5L+INnsjV+Ub/bGf+4PM6qKVHVSbS3pOqcfM/IQ1alSJxnVRephmbo/TPqDFL3xuLO5mfSHqUd1RtV4mWllsSmXs635uP0ZvFXqjOp75kj6oxT9QdLvj+d2v58MBslodM88qccLTMfz69SPEjo/c33776J7xp2kPzxlvP2t8Q9SD0apRsV4u8RSUy57dT3+qpJUoyL1sEqGo/FOeTQcH2weDsf/Phomw2pr+2Ln7RZVndTFeH6fnzGP5/T4uYvx4EejFMNRsj3Wqtoa72j8kSnDOvXW641k++8kE/xMnOkGuIjVqcYf/VUMszZMlje76R8apfXFYymnJlOvrCabmym+eDjDu9eyvDKdlWGZ3miYYTYv0KDHLxRGGaQ3qrMyLLOyPJH5u5bT+sLdyXCUottJceJk6juPpX+ozlK/m7VhMsj48z0v5MECaNJ2LA8zSL+q0hslK8N2VlcmMnd0I+1Dx1NMT6ZYW0sGgxSHjqY+vJze8TIrg07WR0V6o1H6GWRUD1JnlPMR3HVdpcogg615vT4qsrzZzeaxpHNkKcVdh5OqTtFqJcsryaHjGR3tZXVtOqvDVtaHyWY9yqDYTFUNdv6egEvX+KDXqB5mlEEGGWVQ1dmsyqwP2+n1OqlWemktrac4uZR0OklvMxkNx38+uZJ6aTOD1VY2Bu1sVkX6o2RQVxlmlDqj8b7wHM+V7b+DRhlkVAwzqOoM66Q3KrIxaKe/WmZiqZdicTnF3Mz4TmU5XnVzcin10lpGK8Nsbk6nN2qlNyoyqKoMtx6zrrcPBprjpxLdABexuq4yrPtZL9Zzsj+fL65P5Zo7ZnPV1LF0FtdT7JpONgfj4P548sXluRzZLLI6HGRQrKeqhxdk3FU1TL9ez/JgkMO9yXxhcT7z/3A4c8PPp7zqWNJpp1pcT/+f1nP4C/O5c2MyJzfrrBYrGVWbGe+sHTXnMlRXGVWb2cxqVqrNLPa7ObLZzt1Ls5n7zEbmuodSrmykWJhOBsPUh5ez+cnVHD28kLs3JnOyn6wMRlkvVzMcrqeqmj5INX7xXNWDDKqNrJerWR3sysl+J4d7Ezl2aDbdTyynW92R4thi0mqlXlrP6I6lrH66zN3Lcznab2WpX2e1Hv/cVT3ceuFvjnNp2z4wPsxmNjPIxrDO6jBZGrSztDaZXYc2tq5SXqbY6I2XlI+GydJq6i8ez+DOXpYXZ7LY72ZpUGR9VKdXjQ9ODUeb4/md7fA+R/OlrlJVwwzrXjaLXnrVKGuDdlaGZU5sTmTl5GSm7txIZ+5YiiTF0sr4Pd0bveTIyVR3Lad3uMiJteksDVpZGyYbwyob6WVYb6aqB+ftLP2lRHRfzraXrjyAYji86JZyVnXGY64f+J0P9dBHD/DlYrzgejBcy8nOkdy9vjf/b7U7/liw3kT2fmY9k9OLGQ7LLK9M5wvL8/mH5al8YbXK0WIxvdFyqqp/3l/Y1qkzqvvpjU7maLmUL6xN5GNLM6k/c2WuObaSudm1lK06vfVOjq3sy6eXZ/Op1VYObWxmqTia4ajnTDeXpTp1qnqY4WgjvWopi8Vyjm1M585uOwvtmXTuqHL12nLm7jiW1kyRelCnd7LM0cML+fTiQu7YaOfQepUTo42slSczqDa2Dqw1O1/qusqo6qc/Ws1aazEnRntzeKOdz3U7mT2xK/U/Ftl7fDUTe1ZStIoMV+usHJnIF4/vyj+tTufO9SLHe4OsFEvpV6sZVZupL9ABQTg3xmdzq3qYUdXPZrWa9WI1K8PZLPZbOdpv5a61mUzfNciBcjUT60fS2rucTHfGS8qXNjO4s5fFz0/krqW5HO51c7JfZLk/ymq9mY1i+Z6APWf7wzrZfmvL9kG09nJWqz1ZHnRyvN/Kkc1udi/NZfJzg8xnKe2lXoqFyaQskvVBRkc30vtClWOH5nJofWrngNrKcJhesZ7NeiWjqj8+sOZs92lE92VrlKc84/OZyZEH3KJKnb9ce/t5HNODe+fim3PHlWf+7M6TxXJ85h9fLup6mMFoLcuju/K5cm86J67MyqCbz67tyUJndybLKqO6yMqwzNHNInesVvnMxkqOFJ9Nr38yVd3P+d/pVamqfjYGJ3Nk6vP5zNp0ktkc2ZzNvuWZzLaqtIo6varM4qDI3etFPrc6yB31oayODo0/39MLci5LVertF+mjlSx3jubwcC6Ta7Npl+30qoUc3pjKnqObmWyPMqqKrAw6uXtjMl/YaOdzq3UO9/o5Xh7L+uh4BsP18Rw/l2fB7k89TFX3MxiuZ619LMeLXZlZ76ZbTiSZzMqwnSuWZzPbHaQs6vSGrRzfnMidG93csV7mi2ujHB6tZLE8nM3R8ikH1rwg51K2PZ83069Ws9pezFK1K8d77a2r9U+kXS6kP2xl7+JaJufX0pqoU42SwWory4szuWtpLp9fm84XN1o51qtzsj8+ONWrlzIYrWVUjffh5+rgeZ165++g4WgjvdZylouVnNicyGynzEy7lW45nSS5YmM1c4c20p3ppSiT0WaR9aVOTizO5wsrs/nCRjeHNooc3xxlqV7PWrmY/nAtw9HG1kE18/tUovsy9uGlX7/QQ3jINja/mA9svu5CDwMuGnU9ynC0ntXNu3PXRJl+sZGjJ6/IwtJkptqtdMrxhUx6o2FWh4McLRZzpPhslvp3pD9cSl31c/6Xadepq142h0s5ufnZfLpbZXn9YL6wvitzrU66rSJlUWRQ1dkYDnOy2siR8u4crz6btf6RjKr1rfepwmWoHmZU9dIfrGSlPJRDrYmUvWszqGay1G/lronJLHQm0i3rjOoia6MiJzaTIxtV7u5t5u4cyWJ9V3rDxQxGa+dhNUu980J9WG1kY3AiJ7t3pZOJZHVPeqNujm92srvbyUy7Tpk6/arI0qDIsc3k8Powdw/WcqS8K2vDI+kPV7YOFDiwxiWuHkf3cNTLYLSWtfJ4ThbzmdrspFtOpizKDOqpLA3a2bc2ndmjg3TKUaqM3zt9YnMihzc7uWujlbvW6xzaGOR4vZKl8mh6w+UMRxup6v45D9i6Hqaq+hmOetkcLWW5fTRT1VQm1ltpFe1UdSf9ajYnNrvZtdTPVGc4ntejVlYGnRzb7OZQr5M7N4rcvV7laH8jJ8vjWa9O7hwoqOuh6zbci+gGuKiNxgE7WExdV9nsLOd4e1cm6rm0hxNppZMqo/FFXIqNbIxOZnOwlM3h0gWN1zqjjEarWe8no6qf9c6x3N1aSKeaSqvqpEwrowwyrDezmZX0+ovp9U9mMFpNVW3G+7m5PI0Dtqr66Y9WUg7Gn2dbt6ps9K/I8mAuh9rtTLfLdMrxlYx7o1GW+6OcrDZyvDyWxfqurA2PZnOwlFHV24rXpl/cVqmrfkZFL/3hSlbLw0k7GWQza+t7c3xzMnOdViZbZYqMD6itD0dZGg5yvF7OifJQVkaHsjE4mcFo/ZQDBeY5l6rxXL5nFchaNsvlLLUPp1V0UvT2ZlRPZGNYZrHfza5uZ2eVV50iG6MiK8MixzeLHOvVOdIb5OhoNSfKw9moTqY/XNlaEXKu47VK6iJV3c+w2sjmaCVr5bEcLyfSGrWStZkMq3bWhq0c2RwfAJwsx88/qIusj8os9rcOBPZGOdLbzLHieJbrI+mNFk87UGCOn050A1zk6gxTVxvZHAwzGK1mvTiasminKO75HM/t95ZVVX+8w6v6W8F9oXZ4dep6kNFwORtbZ73Lor0z7u0x74y77o9fiNeDePsIl7dq/EK96mVzsDS+pR5ko7WUk8XuTA/nMzWYSiut1KnSzyDr5WrWypNZHx1Pb7g4fkFeraeuz1e8bh8s6GWQJP3xKpxBe7yk9Hi1K5O96UxkIkkyyigbxUbWi+Ws52TWh8ezOVgan+U+bwcKoGnjgB2vAhkfkFov2inaZVIng96urA0nsthvZaZTZqo1/oz7Okl/lKwN6yz3R1kcjM9wnygPZ6U6ko3BiQyGaw3F6z0HC0ZVL4NhO72inaLdSspkNLoivbWZrAw6Od4tM9Mu09n6jO5RnWyMkpV+leXBKCeHvRwvTuRk7s766Hj6g5UMq1MPFJjjpxLdAJeE8UeHjEaDVFl/wK3qrYuvXRxHl+utAwajVNVmihRn3PLiGTc0qR4fEKv6GSXpDbauZt5ayXrreDrlVNrF5PgMeKqM6kGGw/UMqo0MhusZVhsZVb1UVS91Pcr5O0i19XnCW+FdbS2r7bUWs9KaTauYSLuY2NpykGG1mUG1kWG1kcFwLcPtMV/wA4JwrpxytrvqZzBau+egcqtKv9zI6nBXFodTmS476ZZl2sV4P7hZVelV44umrRRLWSlPZK06tnNQbTzH+w3F69bZ7qqfYdbvGXN7lGG5mY16f1Y35zLd72a6bKdTjsc8rOv0qyrr1SCr2chSeTyr9fGsj8YH1Xbe7mIly/0S3QCXjDrJ6BLcjV2q44amnBLedbX1gn09m8X2ipBOiqIcrwbJKFU13FkRUtfDU17Uns9VIVvzuE6q0frWxZh6KYfd9IqTKcvxuJPtlTeDrTEPt1ax9LcOEji4xuVkfDAq6Wf7w3W2P+pzs7WStXIhU5nPRD2VzqibVt1OlSqjYpResZZesZpevZT+cDWbo5XxAarR+s51D5o5QDX++6eos7U0fvWeK7G3NrPZWs1KOZ/JzGaimkpr1EqZcvyZ3ulno1zLZlbTq5ayOVpJf7Cy9V7u3imrb5zlvreiruuz+i9ZFPocLlYP96NXzGu4eH0pH6lkbl9KiiRliqK19c/tj8zc/md1z5W+6+oiWRWyNeYUSVHmnnHfM+YkO+MW26ezz77c3DOHi6KbsuymXU6m3ZpKuzWVTjmVVjExfq/31hy592qQ4Wgjo6p/yhnupoL79HEXaSVFO2XZTaucTLucSru1NfZiIu1yImW2DqZtrbwZ1VvjHm19bY35y/3A2oPNa9ENlwE7cLj8iO4vN/d++8XpATt2Mb6QPXXc9zfm5OIc94Vjn305OvVA1Dhii6KdVtHdWQVSlvf896vramsFy2Drc623VrHUw/McrqceMGinLLbGXXZ3xnzPgcDca8zbq1iGWwcJzvfqm4vLg81rsxcA4IK79wvsS+XF66njvlTGDOfa9tuoiqSuU42GSdFOVfRSVONPKcgp8TpetVJtXVB068KCW28nuefxzte4t1bT1FVGxTBJmVG1dbHW+xl3tTXe8X2GF8nqm4uf6AYAAPiSnR7fdZ0U6Y+/dWq8JjsfBXbho7Xe+t9RUo+SlFsxnbMY8z3358xENwAAwDlzz1LrevstGPX9rQS5mIJ1eyzbFz4tHmDMp27L2RLdAAAAjbhUA/VSHffFqXzwTQAAAICHQ3QDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDRDcAAAA0RHQDAABAQ0Q3AAAANER0AwAAQENENwAAADREdAMAAEBDirqu6ws9CAAAALgcOdMNAAAADRHdAAAA0BDRDQAAAA0R3QAAANAQ0Q0AAAANEd0AAADQENENAAAADRHdAAAA0BDRDQAAAA35/wPs3QkYNFO2MwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x600 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2,4,figsize=(10,6))\n",
    "for i in range(4):\n",
    "    ax[0,i].imshow(input[i].cpu().detach().numpy(), cmap='magma')\n",
    "    ax[1,i].imshow(output_tot[i].cpu().detach().numpy(), cmap='magma')\n",
    "    ax[0,i].set_title(f\"Distance {distances[i]} cm\")\n",
    "[a.axis('off') for a in ax.ravel()]\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, we can combine operators together in complicated ways"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_operator = (ngon_operator + scint_operator) * ngon_operator + scint_operator\n",
    "\n",
    "output_weird= spect_psf_operator(input, xv_k, yv_k, distances, normalize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check normalization. Since there were 2 squares each with $10 \\times 10 $ pixels of magnitude 1, the blurred image should sum to 200 at each plane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([200.0000, 200.0001, 200.0000, 200.0000], device='cuda:0',\n",
       "       grad_fn=<SumBackward1>)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output_weird.sum(dim=(1,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Normalization is carefully tracked for combined operators"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python39_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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}