{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Modeling of the thickness of the sensor\n", "\n", "In this notebook we will re-use the experiment done at ID28 and previously calibrated and model in 3D the detector.\n", "\n", "This detector is a Pilatus 1M with a 450µm thick silicon sensor. Let's first have a look at the absorption coefficients of this sensor material: https://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z14.html\n", "\n", "First we retieve the results of the previous step, then calculate the absorption efficiency:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "# use `widget` instead of `inline` for better user-exeperience. `inline` allows to store plots into notebooks." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "wavelength: 6.968e-11m,\t dist: 2.845e-01m,\t poni1: 8.865e-02m,\t poni2: 8.931e-02m,\t energy: 17.793keV\n" ] } ], "source": [ "import time\n", "start_time = time.perf_counter()\n", "from matplotlib.pyplot import subplots\n", "import numpy\n", "import fabio, pyFAI, pyFAI.units, pyFAI.detectors, pyFAI.azimuthalIntegrator\n", "import json\n", "with open(\"id28.json\") as f:\n", " calib = json.load(f)\n", "\n", "thickness = 450e-6\n", "wavelength = calib[\"wavelength\"]\n", "dist = calib[\"param\"][calib['param_names'].index(\"dist\")]\n", "poni1 = calib[\"param\"][calib['param_names'].index(\"poni1\")]\n", "poni2 = calib[\"param\"][calib['param_names'].index(\"poni2\")]\n", "energy = pyFAI.units.hc/(wavelength*1e10)\n", "print(\"wavelength: %.3em,\\t dist: %.3em,\\t poni1: %.3em,\\t poni2: %.3em,\\t energy: %.3fkeV\" % \n", " (wavelength, dist, poni1, poni2, energy))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Absorption coeficient at 17.8 keV" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "µ = 1537.024174 m^-1 hence absorption efficiency for 450µm: 49.9 %\n" ] } ], "source": [ "# density from https://en.wikipedia.org/wiki/Silicon\n", "rho = 2.3290 # g/cm^3\n", "\n", "#Absorption from https://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z14.html\n", "# Nota: enegies are in MeV !\n", "Si_abs = \"\"\"\n", " 2.00000E-03 2.777E+03 2.669E+03 \n", " 3.00000E-03 9.784E+02 9.516E+02 \n", " 4.00000E-03 4.529E+02 4.427E+02 \n", " 5.00000E-03 2.450E+02 2.400E+02 \n", " 6.00000E-03 1.470E+02 1.439E+02 \n", " 8.00000E-03 6.468E+01 6.313E+01 \n", " 1.00000E-02 3.389E+01 3.289E+01 \n", " 1.50000E-02 1.034E+01 9.794E+00 \n", " 2.00000E-02 4.464E+00 4.076E+00 \n", " 3.00000E-02 1.436E+00 1.164E+00 \n", " 4.00000E-02 7.012E-01 4.782E-01 \n", " 5.00000E-02 4.385E-01 2.430E-01 \n", " 6.00000E-02 3.207E-01 1.434E-01 \n", " 8.00000E-02 2.228E-01 6.896E-02 \n", " 1.00000E-01 1.835E-01 4.513E-02 \n", " 1.50000E-01 1.448E-01 3.086E-02 \n", " 2.00000E-01 1.275E-01 2.905E-02 \n", " 3.00000E-01 1.082E-01 2.932E-02 \n", " 4.00000E-01 9.614E-02 2.968E-02 \n", " 5.00000E-01 8.748E-02 2.971E-02 \n", " 6.00000E-01 8.077E-02 2.951E-02 \n", " 8.00000E-01 7.082E-02 2.875E-02 \n", " 1.00000E+00 6.361E-02 2.778E-02 \n", " 1.25000E+00 5.688E-02 2.652E-02 \n", " 1.50000E+00 5.183E-02 2.535E-02 \n", " 2.00000E+00 4.480E-02 2.345E-02 \n", " 3.00000E+00 3.678E-02 2.101E-02 \n", " 4.00000E+00 3.240E-02 1.963E-02 \n", " 5.00000E+00 2.967E-02 1.878E-02 \n", " 6.00000E+00 2.788E-02 1.827E-02 \n", " 8.00000E+00 2.574E-02 1.773E-02 \n", " 1.00000E+01 2.462E-02 1.753E-02 \n", " 1.50000E+01 2.352E-02 1.746E-02 \n", " 2.00000E+01 2.338E-02 1.757E-02 \"\"\"\n", "data = numpy.array([[float(i) for i in line.split()] for line in Si_abs.split(\"\\n\") if line])\n", "energy_tab, mu_over_rho, mu_en_over_rho = data.T\n", "abs_18 = numpy.interp(energy, energy_tab*1e3, mu_en_over_rho) \n", "mu = abs_18*rho*1e+2\n", "eff = 1.0-numpy.exp(-mu*thickness)\n", "\n", "print(\"µ = %f m^-1 hence absorption efficiency for 450µm: %.1f %%\"%(mu, eff*100))\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAAAx/klEQVR4nO3dd3wVVfrH8c+ThITeQ+9NBASEgCA2VBQbrB107cpa2KauZd1dd93i/lz76qqIfRV7QUVRUbGCBKTX0DsR6TXl+f0xEzfGlEvIzU1yv+/X675yZ+bMuc9kIM+dOWfOMXdHRETiV0KsAxARkdhSIhARiXNKBCIicU6JQEQkzikRiIjEOSUCEZE4p0QgIhLnlAikQjGzC83sg3zLbmadwvePmtkfYxdd5WFmfzaz/8Y6DqkclAik3JnZUWb2lZltM7PvzexLM+sH4O7Pu/tJhe3n7le7+1/LN1qwwPlmNsnMNpnZBjP7wMzOimDfMWa2yMxyzezSAtseNbOd+V77zGxHMXUdb2YzzGy7mS0zs1FlcHj56x9hZivMzAqsTwqP+/Sy/DypOJQIpFyZWV3gHeDfQEOgJfAXYF8s4yqKmSUCLwBXAf8AOgCtgT8Do8zs8YJ/OAuYBVwLzCi4IUxstfNewDjglSLiqAa8ATwG1APOB+41s16lPbZCvAnUB44tsH4o4MD7ZfhZUoEoEUh56wLg7uPcPcfd97j7B+4+G8DMLjWzLwrb0cyeNrO/5VsebmYzw2/IS81saLi+hZmND682Mszsqnz7/NnMXjazZ81sh5nNM7O0YuK9jSBJDXH3Se6+092z3P0r4FSgLnBRUTu7+8PuPgnYW9wvxcxqAWcDzxRRpGH4Wc95YBqwAOhWXL1h3dXMbJyZvWZmyeHv5zUzyzSz5Wb2qzDWvcDLwMUFqrgYeMHds0v6LKmclAikvC0GcszsGTM7xcwalKYSM+sPPAv8juBb7DHAinDzi8AaoAVwDvAPMzs+3+7DwjL1gfHAQ0V8Ri3gSuA6oJqZPRneIvnUzJ4CjgKuB35VmmMo4GwgE/issI3uvpHgiuEyM0s0s4FAW6DQpJnvGGoQfNPfB5wHZANvE1yptAROAH5jZieHuzwDnBPuh5nVA86g6AQlVYASgZQrd99O8AfUgceBzPDbe9MDrOoK4El3/9Ddc919rbsvNLPWwCDgZnff6+4zgbH8+FvuF+4+wd1zgOeAom6vDAQ+cfdd4ee1AjoBlwMnAwnuvhZodICxF+YS4FkvfhTIccCfCP6ofw7c5u6riylfl+B2zlLgsvB4+wGp7n6Hu+9392UE52EEgLt/CWwEzgzrOA9YHP4epYpSIpBy5+4L3P1Sd28F9CD45n7/AVbTmuAPXEEtgO/dPX+j60qCb795NuR7vxuobmZJhdTVBFgbvj8MeNPdt4d/PL8AMLM6wK4DjP1HzKwNcBzBFU5RZboSXMVcDCQD3YGbzOy0YqoeAPQE/pkvwbQFWpjZ1rwX8HsgfyJ+lv8lzouKi0uqBiUCiSl3Xwg8TZAQDsRqoGMh69cBDcM/0Hna8L8/6AfiO6B5+H4O8DMzq2Nm7QmuahoA/wGeLEXd+V0EfBkmmKL0IPhmPjG8AloEvAucUsw+HwB3ApPyXXGtBpa7e/18rzrufmq+/Z4DTghvPw0Ani/tgUnloEQg5crMuprZDWbWKlxuDYwEphxgVU8Q3C8/wcwSzKylmXUNb5V8BdxpZtXNrCfBbZ3S9Kn/Ghgc3i9/gqDdISN8/wHBbZqvKeZqJmycrQ4YQTtDdTMr+P/uYoJkWJxvgc5hF1Izs47A6cDs4nZy97sIej1NMrPGwDfADjO72cxqhO0NPSzsvhvus4Lgimcc8KG7byisbqk6lAikvO0AjgCmmtkuggQwF7jhQCpx92+Ay4D7gG3AZILbHhAklnYEVwdvALe7+0cHGmh4e+kF4P7wfvrl7t7U3Y9390uBfu7+H3fPLaaaD4A9wJHAmPD9MXkbw2/drSik26iZvWdmvw9jWUrQNvEgsD083tcI2j9KOo6/EjQYf0TQ9fR0oDewnOCqZ2y4Pr9nCH6fui0UB0wzlIkULWw7eIXgS9PfgJlATWA4cBMwwN13xixAkTKgRCBSgvBWzqUE38gPBfYDnwD/cPe5MQxNpEwoEYiIxDm1EYiIxLnC+k5XaI0bN/Z27drFOgwRkUpl+vTp37l7amHbKl0iaNeuHenp6bEOQ0SkUjGzlUVt060hEZE4p0QgIhLnlAhEROKcEoGISJxTIhARiXNRSwT5JvEo9MnLcOCsB8MZpGabWZ9oxSIiIkWL5hXB0wRznRblFKBz+BoFPBLFWEREpAhRSwTu/hnwfTFFhhPOyOTuU4D6Zta8mPIHZVnmTv7v/YVoSA0RkR+LZRtBS4JJMvKs4cezSP3AzEaZWbqZpWdmZpbqwyYt2MQjny5l7OfLS7W/iEhVVSkai919jLunuXtaamqhT0iX6Mqj23PqYc24870FfJXxXRlHKCJSecUyEawlmHc2TytKN51gRMyMu87pRcfU2owe9y1rt+6J1keJiFQqsUwE44GLw95DA4Bt7r4+mh9YOyWJxy7qS1Z2Llc/N529WTnR/DgRkUohmt1HxxHM53qIma0xsyvM7GozuzosMgFYRjAH7OPAtdGKJb8OqbW57/zezFm7jd+/PkeNxyIS96I2+qi7jyxhuwPXRevzi3Nit6ZcP6QL9364mO4t63HFUe1jEYaISIVQKRqLo2H04E6c3L0p/5iwgC/VeCwicSxuE0FCgnHPeb3p0LgW170wg9Xf7451SCIiMRG3iQCCxuPHL04jN9e56tl0du3LjnVIIiLlLq4TAUC7xrV46II+LN64g+tfnklurhqPRSS+xH0iADimSyq3ndaNifM28sCkJbEOR0SkXFW6OYuj5fJB7ViwfjsPTFrCIc3qcOphURv2SESkQtEVQcjM+PuZPejTpj7XvzyTuWu3xTokEZFyoUSQT0pSIo9dlEbDmslc9Ww6m7bvjXVIIiJRp0RQQGqdFB6/JI2tu7MYpWEoRCQOKBEUonuLetx3fm9mrt7KTa/O1jAUIlKlKREUYWiPZtw09BDGz1rHg5MyYh2OiEjUqNdQMa45tiMZm3Zy30eLaZ9ai2G9WsQ6JBGRMqcrgmKYGXeedRj92jXgxldmMWPVlliHJCJS5pQISpDXk6hZ3epc9Uy6xiQSkSpHiSACDWsl8+Sl/cjKyeWyp6exbU9WrEMSESkzSgQR6tSkNo9e1JcV3+3i2uenk5WTG+uQRETKRFQTgZkNNbNFZpZhZrcUsr2tmU0ys9lm9qmZtYpmPAfryI6NufOsw/gyYzN/eGOuupWKSJUQzakqE4GHgVOAbsBIM+tWoNjdwLPu3hO4A7gzWvGUlXPTWvPL4zvxUvpqHv5E3UpFpPKL5hVBfyDD3Ze5+37gRWB4gTLdgI/D958Usr1Cun5IF848vCV3f7CYN79dG+twREQOSjQTQUtgdb7lNeG6/GYBZ4XvzwTqmFmjghWZ2SgzSzez9MzMzKgEeyDMjH+efRgDOjTkpldnM2XZ5liHJCJSarFuLL4RONbMvgWOBdYCPxncx93HuHuau6elpqaWd4yFSklK5LGfp9GmUU1GPZvO4o07Yh2SiEipRDMRrAVa51tuFa77gbuvc/ez3P1w4LZw3dYoxlSm6tWsxtOX9SOlWiKXPvkNG7ZptFIRqXyimQimAZ3NrL2ZJQMjgPH5C5hZYzPLi+FW4MkoxhMVrRrU5KlL+7FtTxaXPvUNO/bqGQMRqVyilgjcPRsYDUwEFgAvu/s8M7vDzIaFxY4DFpnZYqAp8PdoxRNNPVrW45Gf9yVj006u/u909mfrGQMRqTyssvWFT0tL8/T09FiHUahXp6/hxldmMaxXC+4/vzcJCRbrkEREADCz6e6eVtg2jT5ahs7p24pNO/Zy1/uLaFInhT+cXvCxCRGRikeJoIxdc2xHNm3fx9gvltOkbgqjjukY65BERIqlRFDGzIw/nt6NzB37+MeEhTSqlcLZfSv0yBkiEueUCKIgMcG49/xebNm9n5tem03DWskM7tok1mGJiBQq1g+UVVnBPAZ9ObR5Ha55fjrTV2pSGxGpmJQIoqhO9Wo8dWl/mtWtzuVPT2PRBj19LCIVjxJBlKXWSeG5K46gerUELnpiKqs2a4YzEalYlAjKQeuGNXn28iPYl53LRU9OZdMODUUhIhWHEkE5OaRZHZ66rB+ZO/Zx8RPfsG23hqIQkYpBiaAc9WnTgMcu6suyzF1c+vQ37NqXHeuQRESUCMrb0Z1TeXDk4cxavZVRz6WzN+sno26LiJQrJYIYGNqjGXed04svMzbzy3HfkpWjQepEJHaUCGLknL6t+Muw7nw4fyM3vDyLnNzKNfifiFQderI4hi45sh279mdz1/uLqFEtkTvPOkwjlopIuVMiiLFrj+vEnv05/PvjDGokJ3L7Gd0wUzIQkfKjRFABXD+kC7v25fDkl8upXi2Rm4ceomQgIuUmqm0EZjbUzBaZWYaZ3VLI9jZm9omZfWtms83s1GjGU1EFI5YeyoVHtOHRyUu5/6MlsQ5JROJI1K4IzCwReBgYAqwBppnZeHefn6/YHwimsHzEzLoBE4B20YqpIjMz/jq8B/uzc3lg0hKSkxK4bnCnWIclInEgmreG+gMZ7r4MwMxeBIYD+ROBA3XD9/WAdVGMp8JLSDD+eXZP9ufk8q+Ji0hOTOCqYzrEOiwRqeKimQhaAqvzLa8BjihQ5s/AB2b2S6AWcGIU46kUEhOMe87tRXaO8/cJC0hMMC4/qn2swxKRKizWzxGMBJ5291bAqcBzZvaTmMxslJmlm1l6ZmZmuQdZ3pISE7h/RG+Gdm/GHe/M55mvVsQ6JBGpwqKZCNYCrfMttwrX5XcF8DKAu38NVAcaF6zI3ce4e5q7p6WmpkYp3IqlWmICD448nCHdmnL7+Hk8+/WKWIckIlVUNBPBNKCzmbU3s2RgBDC+QJlVwAkAZnYoQSKo+l/5I5SclMDDF/ThxEOb8Ke3lAxEJDqilgjcPRsYDUwEFhD0DppnZneY2bCw2A3AVWY2CxgHXOruGmshn+SkBP5zYV8lAxGJGqtsf3fT0tI8PT091mGUu/3ZuVz7/HQ+WrCJvwzrziVHtot1SCJSiZjZdHdPK2xbrBuLJUJ5VwYnhW0GT3yxPNYhiUgVoURQiSQnJfDwhX04pUcz/vrOfMZ8tjTWIYlIFaBEUMnk9SY6rWdz/jFhIQ9/khHrkESkktOgc5VQtcQEHji/N9USjH9NXMS+rBx+O6SLBqoTkVJRIqikkhITuOe83iQnJfDgxxnsy8nllqFdlQxE5IApEVRiiQnGP8/qSXJSAo9NXsa+rFz+dHo3TW4jIgdEiaCSS0gIRi2tnpTI2C+Ws3t/Nnee1ZNEJQMRiZASQRVgZtx22qHUTEniwUlL2JOVy73n9aJaovoCiEjJikwEZnZ9cTu6+71lH46Ulplx/ZAu1ExO5J/vLWTP/mweuqAP1aslxjo0EangivvKWKeEl1RAVx/bkb/+rAeTFm7isqemsXNfdqxDEpEKrsgrAnf/S3kGImXnogFtqZ2SyI2vzObCsVN55rJ+1K+ZHOuwRKSCKrGNwMyqEwwX3Z1gdFAA3P3yKMYlB+nMw1tRO6Ua170wg/Mfm8KzV/Snad3qJe8oInEnktbE54BmwMnAZIJ5BXZEMygpG0O6NeXpy/qxZstuzn7kK1Z8tyvWIYlIBRRJIujk7n8Edrn7M8Bp/HTKSamgjuzYmHGjBrB7fw7nPPo189dtj3VIIlLBRJIIssKfW82sB8Ek802iF5KUtZ6t6vPyLwaSnGic/9jXTFm2OdYhiUgFEkkiGGNmDYA/EswwNh+4K6pRSZnr1KQ2r15zJE3rVefiJ7/h/bkbYh2SiFQQJSYCdx/r7lvcfbK7d3D3Ju7+aHkEJ2WrRf0avPKLgXRvUZdrn5/OC1NXxTokEakAIuk1lAKcDbTLX97d74hg36HAA0AiMNbd/1lg+33A4HCxJtDE3etHGLuUQoNayTx/5RFc+/wMfv/GHDbt2MuvT+iswepE4lgkQ0y8BWwDpgP7Iq3YzBKBh4EhwBpgmpmNd/f5eWXc/bf5yv8SODzS+qX0aiYn8fjFadzy2hzu/2gJG7fv5a/De5CkISlE4lIkiaCVuw8tRd39gQx3XwZgZi8CwwnaGAozEri9FJ8jpVAtMYG7z+1Js3opPPzJUjJ37OPfI/tQI1lDUojEm0i+An5lZoeVou6WwOp8y2vCdT9hZm2B9sDHRWwfZWbpZpaemZlZilCkMGbG707uyh3DuzNp4SZGPj6FzTsjvugTkSoikkRwFDDdzBaZ2Wwzm2Nms8s4jhHAq+6eU9hGdx/j7mnunpaamlrGHy0XD2zHIxf2ZcH67XrwTCQORZIITgE6AycBZwCnhz9LshZonW+5VbiuMCOAcRHUKVEytEczXrhqANv2ZHHWI18xY9WWWIckIuUkkkSwo5DXugj2mwZ0NrP2ZpZM8Md+fMFCZtYVaAB8HWnQEh192zbg9WsHUad6EiPHTOG9OetjHZKIlINIEsEMIBNYDCwJ368wsxlm1reondw9GxgNTAQWAC+7+zwzu8PMhuUrOgJ40d29tAchZad941q8fs2RwbMGL8xg7OfL0KkRqdqspP/kZvY4wf37ieHySQTPFTwFPODu5TruUFpamqenp5fnR8alvVk5XP/yTCbM2cDPB7Thz2d0V/dSkUrMzKa7e1ph2yL5nz0gLwkAuPsHwEB3nwKklFGMUsFUr5bIQyP78ItjO/DfKau44pl0duzNKnlHEal0IkkE683sZjNrG75uAjaGD4zlRjk+iaGEBOPWUw7lzrMO44uM7zj30a9Zs2V3rMMSkTIWSSK4gKDHz5vhq024LhE4L1qBScUxsn8bnrmsP2u37uFnD3+pHkUiVUyJbQQVjdoIYidj004uf3oaG7bv5V/n9GR470KfDxSRCqhUbQRmdn/4820zG1/wFaVYpQLr1KQ2b143iN6t6vPrF2dyzweLyM2tXF8kROSnihtr6Lnw593lEYhUDg1rJfPfK4/gD2/O4d8fZ5CxaSf3nNeLmsmRDFslIhVRkf973X16+HNy3rpwgprW7l7WQ0xIJZKclMD/nd2TLk3r8I8JC1j5yG4evySNlvVrxDo0ESmFEhuLzexTM6trZg0JHi573MzujX5oUpGZGVce3YEnLu3H6u93M/yhL5i24vtYhyUipRBJr6F67r4dOAt4NnyA7MTohiWVxeBDmvDGdYOoW70aFzw+hXHfaNYzkcomkkSQZGbNCbqKvhPleKQS6tSkNm9cN4iBHRtz6+tz+OObc9mfrUdMRCqLSBLBHQTjBWW4+zQz60Aw5pDID+rVqMZTl/bjF8d24LkpK/n52Klk7tDcBiKVgZ4jkDI3ftY6bnp1Fg1qJvPIz/vSu3X9WIckEvcOdqwhkQMyrFcLXrvmSBITjPMe/ZqXpqndQKQiUyKQqOjeoh5vjz6KIzo05ObX5vD7N+awL7vQCehEJMaUCCRqGtRK5unL+nPNcR15YeoqzntsCuu27ol1WCJSQJEPlJnZ9cXt6O56lkBKlJhg3Dy0K71a1ePGV2Zz+r+/4N8jD2dQp8axDk1EQsVdEdQp4VUiMxsaTnqfYWa3FFHmPDObb2bzzOyFAwtfKouhPZrz1uhBNKqVzEVPTOXhTzI0TpFIBRG1XkPhfAWLgSHAGoI5jEe6+/x8ZToDLwPHu/sWM2vi7puKq1e9hiq3XfuyueX1Obw9ax3Hd23Cvef1on7N5FiHJVLlFddrqMSRwsysOnAF0B2onrfe3S8vYdf+BM8eLAvreREYDszPV+Yq4GF33xLWWWwSkMqvVkoSD47oTVrbBvzt3fmc/u8v+M+FfejZqn6sQxOJW5E0Fj8HNANOBiYTTFKzI4L9WgKr8y2vCdfl1wXoYmZfmtkUMxtaWEVmNsrM0s0sPTMzM4KPlorMzLjkyHa8/IuB5OY65zzyNc9+vYLK9kyLSFURSSLo5O5/BHa5+zPAaUBZTVifBHQGjgNGEgxoV79gIXcf4+5p7p6WmppaRh8tsXZ4mwa8+6ujGdSpEX96ax6jx32reZFFYiCSRJD3P3OrmfUA6gFNIthvLdA633KrcF1+a4Dx7p7l7ssJ2hQ6R1C3VBENaiXzxCX9uHloV96fu4Ez/v0Fc9dui3VYInElkkQwJpyH4I/AeIJ7/HdFsN80oLOZtTezZGBEuH9+bxJcDWBmjQluFS2LKHKpMhISjGuO68i4qwawNyuXs/7zlW4ViZSjEhOBu4919y3uPtndO7h7E3d/NIL9soHRBAPWLQBedvd5ZnaHmQ0Li00ENpvZfOAT4Hfuvrn0hyOVWf/2DZnw6//dKrrmvzPYtlu3ikSircTuo2b2p8LWu/sdUYmoBOo+WvXl5jpjv1jGXe8vomnd6jwwojdp7RrGOiyRSu1gB53ble+VA5wCtCuz6EQKSEgwRh3TkVfDgevOHzOFhz5eQo4eQBOJigN+oMzMUoCJ7n5cVCIqga4I4suOvVnc9sZcxs9ax4AODbnv/N40r6e5kUUOVFkPQ12ToAeQSNTVqV6NB0b05u5zezF7zTZOeeBz3p+7IdZhiVQpkUxeP8fMZoevecAi4P6oRyYSMjPO6duKd391NK0b1OTq/07n1tdns3t/dqxDE6kSShxiAjg93/tsYGPYI0ikXLVvXIvXrjmS+z5azKOTlzJ12ffcP6K3hqcQOUhFXhGYWUMza0gwnETeaw9QN1wvUu6SkxK4eWhXXrhyAHuycjjrP1/x0MdLyM7JjXVoIpVWcbeGpgPp4c9Mgqd+l4Tvp0c/NJGiDezYiPd/fQynHtacuz9YzPljprBq8+5YhyVSKRWZCNy9vbt3AD4CznD3xu7eiOBW0QflFaBIUerVrMaDIw/ngRG9WbxxB0Mf+Ixx36zSE8kiByiSXkMD3H1C3oK7vwccGb2QRA7M8N4tmfibY+jduj63vj6HK55JZ9OOvbEOS6TSiCQRrDOzP5hZu/B1G7Au2oGJHIgW9Wvw3yuO4PYzuvFlxnecfN9nvDNb/0xFIhFJIhgJpAJvhK8m4TqRCiUhwbhsUHve/dXRtGlYk9EvfMvoF2awZdf+WIcmUqFFbarKaNGTxRKJ7JxcHvtsGfd/tJh6NZL5x5k9OKl7s1iHJRIzpXqy2MzuD3++bWbjC76iFKtImUhKTOC6wZ0YP/oomtRJYdRz0/nNi9+ydbeuDkQKKu6BsufCn3eXRyAi0XBo87q8NXoQD3+SwUMfZ/Dl0s387Wc9OFlXByI/OKBbQ+EENa3dfXb0Qiqebg1Jac1bt43fvTKb+eu3c3rP5vxlWHca1U6JdVgi5eKgBp0zs0/NLO9p4hkE8wrfW9ZBikRb9xb1eGv0IG4Y0oWJ8zYw5L7PeGvmWj13IHEvkl5D9dx9O3AW8Ky7HwGcGEnlZjbUzBaZWYaZ3VLI9kvNLNPMZoavKw8sfJEDUy0xgV+e0DkYwK5hTX794kyufCad9dv2xDo0kZiJJBEkmVlz4DzgnUgrNrNE4GGCiWy6ASPNrFshRV9y997ha2yk9YscjC5N6/D6NUfyh9MO5aulmxly72c8N2UluZr8RuJQJIngDoK5hZe6+zQz60Aw5lBJ+gMZ7r7M3fcDLwLDSx+qSNlKTDCuPLrDD08l//HNuZz32Ncs2bgj1qGJlKtIJq9/xd17uvs14fIydz87grpbAqvzLa8J1xV0djjXwatm1rqwisxslJmlm1l6ZmZmBB8tErk2jWry3BX9ufvcXmRk7uTUBz/n3g8XszcrJ9ahiZSLSBqLu5jZJDObGy73NLM/lNHnvw20c/eewIfAM4UVcvcx7p7m7mmpqall9NEi/5M3+c1H1x/LaYc158FJSzj1gc/5aul3sQ5NJOoiuTX0OHArkAUQdh0dEcF+a4H83/Bbhet+4O6b3X1fuDgW6BtBvSJR07h2CvePOJxnL+9Pdq5zweNTueHlWWzeua/knUUqqUgSQU13/6bAukhmKJsGdDaz9maWTJA8fvREctgInWcYsCCCekWi7pguqUz8zTFce1xH3pq5luPvmcyL36xSY7JUSZEkgu/MrCPgAGZ2DrC+pJ3C6SxHEzQ0LwBedvd5ZnaHmQ0Li/3KzOaZ2SzgV8ClpTgGkaiokZzITUO7MuHXR3NI0zrc8vocznn0K+av2x7r0ETKVIlPFoe9hMYQzEGwBVgOXOjuK6Mf3k/pyWKJBXfn1elruPO9hWzbk8UlA9vx2yGdqVO9WqxDE4nIQT1ZHPYSOpFgKOquwLHAUWUbokjFZmacm9aaj284lhH9WvPUV8s5/p7JvPmtnkyWyq+40UfrmtmtZvaQmQ0BdgOXABkED5eJxJ36NZP5+5mH8ea1g2herzq/eWkm54+ZwsINul0klVeRt4bM7C2CW0FfAycQTEhjwK/dfWZ5BViQbg1JRZGT67w0bTV3TVzIjr3ZXDSgLb8d0oV6NXS7SCqe4m4NFZcI5rj7YeH7RIIG4jbuHtPJYJUIpKLZsms/93y4iBemrqJBzWR+d/IhnJvWmsQEi3VoIj8obRtBVt4bd88B1sQ6CYhURA1qJfO3nx3G+NFH0SG1Fre8PofhD39B+orvYx2aSESKuyLIAXblLQI1CNoJDHB3r1suERagKwKpyNyd8bPWceeEhWzYvpdhvVpw8yldaVm/RqxDkzhX3BVBkTOUuXti9EISqZrMjOG9WzKkW1Me+XQpYz5bxgfzNzDqmI5cfWwHaiYXNymgSGxE8kCZiBygmslJ3HDSIUy64VhOPLQpD05awuC7P+XV6Wv0dLJUOEoEIlHUqkFNHrqgD69ePZBmdatz4yuzGPbwF0xZtjnWoYn8QIlApByktWvIG9cO4oERvfl+535GjJnCVc+mszRzZ6xDE1EiECkvCQlB+8HHNx7H704+hK+Xbuak+z7jj2/OJXOHRjeV2ClxrKGKRr2GpKr4buc+7v9oMeO+WU31pAR+cWxHrjy6vRqUJSpK9UBZRaVEIFXN0syd3PX+QibO20hqnRR+c2JnzktrTbVEXbBL2TmoQedEJLo6ptbmsYvSePXqgbRtWJPb3pjLSfd9xruz12tAOykXSgQiFURau4a8cvVAHr84jaQE47oXZjD84S/5fEmmEoJElRKBSAViZgzp1pT3f3MM/zqnJ5t37ueiJ77hwrFTmbFqS6zDkyoqqonAzIaa2SIzyzCzW4opd7aZuZkVev9KJN4kJoTzH9x4LLef0Y1FG3Zw1n++4spnpmmGNClzUWssDkcsXQwMAdYQzGE80t3nFyhXB3gXSAZGu3uxLcFqLJZ4tGtfNk99uZzHPlvGjr3ZnNazOb89sTOdmtSJdWhSScSqsbg/kBHOcLYfeBEYXki5vwL/B2hkU5Ei1EpJYvTxnfnipuO5bnBHPlm4iZPu+4zfvjSTFd/tKrkCkWJEMxG0BFbnW14TrvuBmfUBWrv7u8VVZGajzCzdzNIzMzPLPlKRSqJezWr87uSufH7TYK48ugPvzV3PCfdO5sZXZrFq8+5YhyeVVMwai80sAbgXuKGksu4+xt3T3D0tNTU1+sGJVHCNaqfw+1MP5bObBnPJwHa8PWsdx9/zKTe9qoQgBy6aiWAt0DrfcqtwXZ46QA/gUzNbAQwAxqvBWCRyTepU509ndOPzmwbz8wFteXPmOgbf8ym/e2UWKzfrlpFEJpqNxUkEjcUnECSAacAF7j6viPKfAjeqsVik9DZu38ujk5fywtRVZOXk8rPeLbl2cCc6Nakd69AkxmLSWOzu2cBoYCKwAHjZ3eeZ2R1mNixanysSz5rWrc7tZ3Tn85sHc8VR7Xlv7gaG3DeZ0S/MYMF6dTuVwmmsIZEqbPPOfYz9YjnPfb2SnfuyOfHQJlw3uBOHt2kQ69CknGnQOZE4t213Fk9/tYInv1zOtj1ZHNmxEdcN7sSRHRthZrEOT8qBEoGIALBzXzbjpq7i8c+XsWnHPnq1qsc1x3VkSLdmJCYoIVRlSgQi8iN7s3J4bcYaxny2jJWbd9OhcS1GHdOBM/u0JCUpMdbhSRQoEYhIoXJynffmrufRyUuZu3Y7qXVSuGxQOy48oi31alSLdXhShpQIRKRY7s6XGZt57LOlfL7kO2olJzKifxsuG9SOVg1qxjo8KQNKBCISsXnrtjH28+W8PWsdDpx6WHOuOro9PVvVj3VochCUCETkgK3buoenvlzOuG9Ws3NfNv3bNeTyo9ozpFtTNSxXQkoEIlJqO/Zm8dK01Tz15QrWbt1Dm4Y1ufTIdpzXrzW1U5JiHZ5ESIlARA5adk4uE+dt5MkvlzN95RbqpCRxblprLjmyLW0b1Yp1eFICJQIRKVMzV2/lyS+WM2HOenLcOaFrEy45sh1HdWqsB9QqKCUCEYmKjdv38vyUlTw/dRWbd+2nY2otLjmyHWf1aaXbRhWMEoGIRNXerBzenb2eZ75ewew126idksRZfVpy8cC2mk6zglAiEJFy4e58u3orz329kndnr2d/Ti4DOzTiooFtGdKtKdUSYzYXVtxTIhCRcvfdzn28NG01L0xdxdqte2hSJ4UR/dswol9rWtSvEevw4o4SgYjETE6u8+miTTw3ZSWTF2diwPFdm3DBEW04tksTPZNQTopLBGrNEZGoSkwwTji0KScc2pTV3+/mxWmreGnaGj5akE6LetU5v18bzuvXiub1dJUQK1G9IjCzocADQCIw1t3/WWD71cB1QA6wExjl7vOLq1NXBCKV3/7sXD5asJFx36zi8yXfkWAw+JAmjOjfhsGHpJKktoQyF5NbQ2aWSDBn8RBgDcGcxSPz/6E3s7ruvj18Pwy41t2HFlevEoFI1bJy8y5emraaV6avIXPHPprUSeGcvq04L6017RrrQbWyEqtbQ/2BDHdfFgbxIjAc+CER5CWBUC2gcjVYiMhBa9uoFjcN7cpvh3Th44WbeHnaah6dvJT/fLqUI9o35Ly01pxyWDNqJutOdrRE8zfbElidb3kNcETBQmZ2HXA9kAwcH8V4RKQCq5aYwMndm3Fy92Zs3L6XV6ev4ZX01dzwyixuHz+P03s255y+rejbtoGeXi5j0bw1dA4w1N2vDJcvAo5w99FFlL8AONndLylk2yhgFECbNm36rly5Mioxi0jF4u5MW7GFl9NXM2HOenbvz6F941qc3aclZ/ZpRUt1Q41YrNoIBgJ/dveTw+VbAdz9ziLKJwBb3L1ecfWqjUAkPu3al817czfwSvpqpi7/HjMY0L4RZ/dtxdAezTSkRQlilQiSCBqLTwDWEjQWX+Du8/KV6ezuS8L3ZwC3FxVoHiUCEVn9/W5en7GW12asYdX3u6lRLZGTuzflzD6tGNSxkXodFSJmD5SZ2anA/QTdR59097+b2R1AuruPN7MHgBOBLGALMDp/oiiMEoGI5HF3pq/cwuvfruWdWevYvjeb1DopnNGzBWce3pIeLeuqPSGkJ4tFpMrbm5XDJws38ebMtXy8cBNZOU6H1FoM79WS4b1bxH1XVCUCEYkrW3fvZ8KcDYyftZapy7/HHXq1qscZvVpwRq8WNK1bPdYhljslAhGJW+u27uGd2esYP2sdc9duxwz6t2vIGb1acEqPZjSqnRLrEMuFEoGICLA0cydvz1rH27PWsTRzF4kJxpEdG3F6z+ac3L0Z9WsmxzrEqFEiEBHJx91ZuGEHb89ax7tz1rNy826SEoxBnRpz2mHNOal70yqXFJQIRESK4O7MW7edt2evY8Kc9az+fg9JCcbAjo049bDmnNStaZW4faREICISgbyk8O6c9UwIrxQSDI5o34ihPYLhL5rVq5wNzUoEIiIHyN1ZsH4H781dz/tzN7Bk004AereuH46J1JQOqbVjHGXklAhERA5SxqadTJy3gYnzNjB7zTYAOjepzUndm3JSt2Yc1rIeCRV4tjUlAhGRMrR26x4+mLeBD+Zt5JsV35OT6zStm8KJhzZlSLemDOzYiJSkxFiH+SNKBCIiUbJ1934mLdjEh/M38tmSTHbvz6FWciJHd07lxG5NGXxIaoVobFYiEBEpB3uzcvh66WY+XLCRSQs2snH7PsygT5sGHN+1CScc2oRDmtaJyfhHSgQiIuXM3Zm7djuTFm5k0oJNzFkbtCu0rF+D4w5J5fiuTRjYsVG5zbymRCAiEmMbt+/lk4Wb+HjhJr7I+I7d+3NITkpgQIdGDD4kleMOaUL7KA6Mp0QgIlKB7MvO4Zvl3/Ppokw+WbSJZZm7AGjbqCbHdUnl2ENSGdChbK8WlAhERCqwVZt38+niTUxelMlXSzezJyuH5MQE+rVvwDGdUzm6cyqHNj+4tgUlAhGRSmJfdg7pK7YweXEmkxdlsmjjDgBS66Twh9MOZXjvlqWqt7hEoEk+RUQqkJSkRAZ1asygTo35/amHsmHbXj5fkslnS76L2jwKUZ3Y08yGmtkiM8sws1sK2X69mc03s9lmNsnM2kYzHhGRyqZZveqcm9aaf488nAEdGkXlM6KWCMwsEXgYOAXoBow0s24Fin0LpLl7T+BV4K5oxSMiIoWL5hVBfyDD3Ze5+37gRWB4/gLu/om77w4XpwCtohiPiIgUIpqJoCWwOt/ymnBdUa4A3itsg5mNMrN0M0vPzMwswxBFRCSqbQSRMrOfA2nAvwrb7u5j3D3N3dNSU1PLNzgRkSoumr2G1gKt8y23Ctf9iJmdCNwGHOvu+6IYj4iIFCKaVwTTgM5m1t7MkoERwPj8BczscOAxYJi7b4piLCIiUoSoJQJ3zwZGAxOBBcDL7j7PzO4ws2FhsX8BtYFXzGymmY0vojoREYmSqD5Q5u4TgAkF1v0p3/sTo/n5IiJSsko3xISZZQIrS7l7Y+C7MgynMtAxxwcdc3w4mGNu6+6F9rapdIngYJhZelFjbVRVOub4oGOOD9E65grRfVRERGJHiUBEJM7FWyIYE+sAYkDHHB90zPEhKsccV20EIiLyU/F2RSAiIgUoEYiIxLm4SQQlTZJTGZlZazP7JJzcZ56Z/Tpc39DMPjSzJeHPBuF6M7MHw9/BbDPrE9sjKD0zSzSzb83snXC5vZlNDY/tpXBYE8wsJVzOCLe3i2ngpWRm9c3sVTNbaGYLzGxgVT/PZvbb8N/1XDMbZ2bVq9p5NrMnzWyTmc3Nt+6Az6uZXRKWX2JmlxxoHHGRCCKcJKcyygZucPduwADguvC4bgEmuXtnYFK4DMHxdw5fo4BHyj/kMvNrgqFL8vwfcJ+7dwK2EAxrTvhzS7j+vrBcZfQA8L67dwV6ERx7lT3PZtYS+BXBxFU9gESC8cqq2nl+GhhaYN0BnVczawjcDhxBMA/M7XnJI2LuXuVfwEBgYr7lW4FbYx1XFI7zLWAIsAhoHq5rDiwK3z8GjMxX/odylelFMJLtJOB44B3ACJ62TCp4vgnGuhoYvk8Ky1msj+EAj7cesLxg3FX5PPO/+UwahuftHeDkqniegXbA3NKeV2Ak8Fi+9T8qF8krLq4IOPBJciqd8FL4cGAq0NTd14ebNgBNw/dV5fdwP3ATkBsuNwK2ejDQIfz4uH445nD7trB8ZdIeyASeCm+HjTWzWlTh8+zua4G7gVXAeoLzNp2qfZ7zHOh5PejzHS+JoEozs9rAa8Bv3H17/m0efEWoMn2Ezex0YJO7T491LOUoCegDPOLuhwO7+N/tAqBKnucGBFPbtgdaALX46S2UKq+8zmu8JIKIJsmpjMysGkESeN7dXw9XbzSz5uH25kDeXA9V4fcwCBhmZisI5sE+nuD+eX0zyxtNN/9x/XDM4fZ6wObyDLgMrAHWuPvUcPlVgsRQlc/zicByd8909yzgdYJzX5XPc54DPa8Hfb7jJRGUOElOZWRmBjwBLHD3e/NtGg/k9Ry4hKDtIG/9xWHvgwHAtnyXoJWCu9/q7q3cvR3BefzY3S8EPgHOCYsVPOa838U5YflK9c3Z3TcAq83skHDVCcB8qvB5JrglNMDMaob/zvOOucqe53wO9LxOBE4yswbhldRJ4brIxbqhpBwbZE4FFgNLgdtiHU8ZHdNRBJeNs4GZ4etUgnujk4AlwEdAw7C8EfSeWgrMIeiREfPjOIjjPw54J3zfAfgGyABeAVLC9dXD5Yxwe4dYx13KY+0NpIfn+k2gQVU/z8BfgIXAXOA5IKWqnWdgHEEbSBbBld8VpTmvwOXhsWcAlx1oHBpiQkQkzsXLrSERESmCEoGISJxTIhARiXNKBCIicU6JQEQkzikRSJVmZjlmNjMcxXKWmd1gZqX+d29mv8/3vl3+USNL2O83ZnZxaT+3QF13m9nxZVGXCGiGMqnizGynu9cO3zcBXgC+dPfby6C+dgTPMfQoYZ8kYAbQx/83Tk6pmVlb4HF3P+lg6xIBXRFIHHH3TQTD944On85MNLN/mdm0cHz3XwCY2XFm9pmZvWvBHBaPmlmCmf0TqBFeYTwfVptoZo+HVxwfmFmNQj76eGBGXhIws0/NLC183zgcLgMzu9TM3gzHoF9hZqPN7PpwoLkp4XDDuPtKoJGZNYvm70vihxKBxBV3X0Ywtn0Tgqc4t7l7P6AfcJWZtQ+L9gd+STB/RUfgLHe/Bdjj7r09GNYCgrHhH3b37sBW4OxCPnYQwciZkegBnBXG83dgtwcDzX0N5L+1NCOsV+SgKRFIPDuJYOyWmQTDdzci+MMO8I27L3P3HIJhAI4qoo7l7j4zfD+dYGz5gpoTDCMdiU/cfYe7ZxIMpfx2uH5Ogbo3EYzKKXLQkkouIlJ1mFkHIIfgD6kBv3T3iQXKHMdPh/4tqjFtX773OUBht4b2EIyF86OPCX9WK6a+3HzLufz4/2v1sF6Rg6YrAokbZpYKPAo85EEviYnANeFQ3phZl3DCF4D+4Wi1CcD5wBfh+qy88gdgAdCpwLp+4c/jCG5VHaguBIOxiRw0JQKp6vIad+cRjOT4AcGolgBjCYY2nhF2A32M/33rngY8RPBHfDnwRrh+DDA7X2NxJN4Djimw7kQzm0Yw7v73ZvarSCsLE1EngtFIRQ6auo+KFBDeGrrR3U8vwzrfAG5y9yVm9mlYf6n+kJvZmQRdUf9YVvFJfNMVgUj5uIWg0bgsJAH3lFFdIroiEBGJd7oiEBGJc0oEIiJxTolARCTOKRGIiMQ5JQIRkTj3//CbCXSUErsYAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "depth = numpy.linspace(0, 1000, 100)\n", "res = numpy.exp(-mu*depth*1e-6)\n", "fig, ax = subplots()\n", "ax.plot(depth, res, \"-\")\n", "ax.set_xlabel(\"Depth (µm)\")\n", "ax.set_ylabel(\"Residual signal\")\n", "ax.set_title(\"Silicon @ 17.8 keV\")\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is consistent with:\n", "http://henke.lbl.gov/optical_constants/filter2.html\n", "\n", "Now we can model the detector\n", "\n", "## Modeling of the detector:\n", "\n", "The detector is seen as a 2D array of voxel. Let vox, voy and voz be the dimention of the detector in the three dimentions.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Detector Pilatus 1M\t PixelSize= 1.720e-04, 1.720e-04 m\n", "0.000172 0.000172 0.00045\n" ] } ], "source": [ "detector= pyFAI.detector_factory(calib[\"detector\"])\n", "print(detector)\n", "\n", "vox = detector.pixel2 # this is not a typo\n", "voy = detector.pixel1 # x <--> axis 2\n", "voz = thickness\n", "\n", "print(vox, voy, voz)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The intensity grabbed in this voxel is the triple integral of the absorbed signal coming from this pixel or from the neighboring ones.\n", "\n", "There are 3 ways to perform this intergral:\n", "* Volumetric analytic integral. Looks feasible with a change of variable in the depth\n", "* Slice per slice, the remaining intensity depand on the incidence angle + pixel splitting between neighbooring pixels\n", "* raytracing: the decay can be solved analytically for each ray, one has to throw many ray to average out the signal.\n", "\n", "For sake of simplicity, this integral will be calculated numerically using this raytracing algorithm.\n", "http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.42.3443&rep=rep1&type=pdf\n", "\n", "Knowing the input position for a X-ray on the detector and its propagation vector, this algorithm allows us to calculate the length of the path in all voxel it crosses in a fairly efficient way.\n", "\n", "To speed up the calculation, we will use a few tricks:\n", "* One ray never crosses more than 16 pixels, which is reasonable considering the incidance angle \n", "* we use numba to speed-up the calculation of loops in python\n", "* We will allocate the needed memory by chuncks of 1 million elements\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from numba import jit \n", "\n", "BLOCK_SIZE = 1<<20 # 1 million\n", "BUFFER_SIZE = 16 \n", "BIG = numpy.finfo(numpy.float32).max\n", "\n", "mask = numpy.load(\"mask.npy\").astype(numpy.int8)\n", "from scipy.sparse import csr_matrix, csc_matrix, linalg" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array([0, 0, 1, 1], dtype=int32), array([0, 1, 1, 2], dtype=int32), array([0.00029791, 0.00029791, 0.00059583, 0.00059583], dtype=float32))\n", "The slowest run took 10.11 times longer than the fastest. This could mean that an intermediate result is being cached.\n", "6.04 µs ± 7.92 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "5.92 µs ± 46.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "@jit\n", "def calc_one_ray(entx, enty, \n", " kx, ky, kz,\n", " vox, voy, voz):\n", " \"\"\"For a ray, entering at position (entx, enty), with a propagation vector (kx, ky,kz),\n", " calculate the length spent in every voxel where energy is deposited from a bunch of photons comming in the detector \n", " at a given position and and how much energy they deposit in each voxel. \n", " \n", " Direct implementation of http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.42.3443&rep=rep1&type=pdf\n", " \n", " :param entx, enty: coordinate of the entry point in meter (2 components, x,y)\n", " :param kx, ky, kz: vector with the direction of the photon (3 components, x,y,z)\n", " :param vox, voy, voz: size of the voxel in meter (3 components, x,y,z)\n", " :return: coordinates voxels in x, y and length crossed when leaving the associated voxel\n", " \"\"\"\n", " array_x = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n", " array_x[:] = -1\n", " array_y = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n", " array_y[:] = -1\n", " array_len = numpy.empty(BUFFER_SIZE, dtype=numpy.float32)\n", " \n", " #normalize the input propagation vector\n", " n = numpy.sqrt(kx*kx + ky*ky + kz*kz)\n", " kx /= n\n", " ky /= n\n", " kz /= n\n", " \n", " #assert kz>0\n", " step_X = -1 if kx<0.0 else 1\n", " step_Y = -1 if ky<0.0 else 1\n", " \n", " #assert vox>0\n", " #assert voy>0\n", " #assert voz>0\n", " \n", " X = int(entx//vox)\n", " Y = int(enty//voy)\n", " \n", " if kx>0.0:\n", " t_max_x = ((entx//vox+1)*(vox)-entx)/ kx\n", " elif kx<0.0:\n", " t_max_x = ((entx//vox)*(vox)-entx)/ kx\n", " else:\n", " t_max_x = BIG\n", "\n", " if ky>0.0:\n", " t_max_y = ((enty//voy+1)*(voy)-enty)/ ky\n", " elif ky<0.0:\n", " t_max_y = ((enty//voy)*(voy)-enty)/ ky\n", " else:\n", " t_max_y = BIG\n", " \n", " #Only one case for z as the ray is travelling in one direction only\n", " t_max_z = voz / kz\n", " \n", " t_delta_x = abs(vox/kx) if kx!=0 else BIG\n", " t_delta_y = abs(voy/ky) if ky!=0 else BIG\n", " t_delta_z = voz/kz\n", " \n", " finished = False\n", " last_id = 0\n", " array_x[last_id] = X\n", " array_y[last_id] = Y\n", " \n", " while not finished:\n", " if t_max_x < t_max_y:\n", " if t_max_x < t_max_z:\n", " array_len[last_id] = t_max_x\n", " last_id+=1\n", " X += step_X\n", " array_x[last_id] = X\n", " array_y[last_id] = Y\n", " t_max_x += t_delta_x\n", " else:\n", " array_len[last_id] = t_max_z\n", " finished = True\n", " else:\n", " if t_max_y < t_max_z:\n", " array_len[last_id] = t_max_y\n", " last_id+=1\n", " Y += step_Y\n", " array_x[last_id] = X\n", " array_y[last_id] = Y \n", " t_max_y += t_delta_y\n", " else:\n", " array_len[last_id] = t_max_z\n", " finished = True\n", " if last_id>=array_len.size-1:\n", " print(\"resize arrays\")\n", " old_size = len(array_len)\n", " new_size = (old_size//BUFFER_SIZE+1)*BUFFER_SIZE\n", " new_array_x = numpy.empty(new_size, dtype=numpy.int32)\n", " new_array_x[:] = -1\n", " new_array_y = numpy.empty(new_size, dtype=numpy.int32)\n", " new_array_y[:] = -1\n", " new_array_len = numpy.empty(new_size, dtype=numpy.float32)\n", " new_array_x[:old_size] = array_x\n", " new_array_y[:old_size] = array_y\n", " new_array_len[:old_size] = array_len\n", " array_x = new_array_x\n", " array_y = new_array_y\n", " array_len = new_array_len\n", " return array_x[:last_id], array_y[:last_id], array_len[:last_id]\n", "\n", "print(calc_one_ray(0.0,0.0, 1,1,1, 172e-6, 172e-6, 450e-6))\n", "import random\n", "%timeit calc_one_ray(10+random.random(),11+random.random(),\\\n", " random.random()-0.5,random.random()-0.5,0.5+random.random(), \\\n", " vox, voy, voz)\n", "%timeit calc_one_ray.py_func(10+random.random(),11+random.random(),\\\n", " random.random()-0.5,random.random()-0.5,0.5+random.random(), \\\n", " vox, voy, voz)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we are able to perform raytracing for any ray comming in the detector, we can calculate the contribution to the neighboring pixels, using the absorption law (the length travelled is already known). \n", "To average-out the signal, we will sample a few dozens of rays per pixel to get an approximatation of the volumic integrale. \n", "\n", "Now we need to store the results so that this transformation can be represented as a sparse matrix multiplication:\n", "\n", "b = M.a\n", "\n", "Where b is the recorded image (blurred) and a is the \"perfect\" signal. \n", "M being the sparse matrix where every pixel of a gives a limited number of contribution to b.\n", "\n", "Each pixel in *b* is represented by one line in *M* and we store the indices of *a* of interest with the coefficients of the matrix.\n", "So if a pixel i,j contributes to (i,j), (i+1,j), (i+1,j+1), there are only 3 elements in the line. \n", "This is advantagous for storage.\n", "\n", "We will use the CSR sparse matrix representation:\n", "https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_.28CSR.2C_CRS_or_Yale_format.29\n", "where there are 3 arrays:\n", "* data: containing the actual non zero values\n", "* indices: for a given line, it contains the column number of the assocated data (at the same indice)\n", "* idptr: this array contains the index of the start of every line.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from numba.experimental import jitclass\n", "from numba import int8, int32, int64, float32, float64\n", "spec = [(\"vox\",float64),(\"voy\",float64),(\"voz\",float64),(\"mu\",float64),\n", " (\"dist\",float64),(\"poni1\",float64),(\"poni2\",float64),\n", " (\"width\", int64),(\"height\", int64),(\"mask\", int8[:,:]),\n", " (\"sampled\", int64), (\"data\", float32[:]),(\"indices\", int32[:]),(\"idptr\", int32[:]),\n", " ]\n", "@jitclass(spec)\n", "class ThickDetector(object):\n", " \"Calculate the point spread function as function of the geometry of the experiment\"\n", " \n", " def __init__(self, vox, voy, thickness, mask, mu, \n", " dist, poni1, poni2):\n", " \"\"\"Constructor of the class:\n", " \n", " :param vox, voy: detector pixel size in the plane\n", " :param thickness: thickness of the sensor in meters\n", " :param mask: \n", " :param mu: absorption coefficient of the sensor material\n", " :param dist: sample detector distance as defined in the geometry-file\n", " :param poni1, poni2: coordinates of the PONI as defined in the geometry \n", " \"\"\"\n", " self.vox = vox\n", " self.voy = voy\n", " self.voz = thickness\n", " self.mu = mu\n", " self.dist=dist\n", " self.poni1 = poni1\n", " self.poni2 = poni2\n", " self.width = mask.shape[-1]\n", " self.height = mask.shape[0]\n", " self.mask = mask\n", " self.sampled = 0\n", " self.data = numpy.zeros(BLOCK_SIZE, dtype=numpy.float32)\n", " self.indices = numpy.zeros(BLOCK_SIZE,dtype=numpy.int32)\n", " self.idptr = numpy.zeros(self.width*self.height+1, dtype=numpy.int32)\n", " \n", " def calc_one_ray(self, entx, enty):\n", " \"\"\"For a ray, entering at position (entx, enty), with a propagation vector (kx, ky,kz),\n", " calculate the length spent in every voxel where energy is deposited from a bunch of photons comming in the detector \n", " at a given position and and how much energy they deposit in each voxel. \n", "\n", " Direct implementation of http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.42.3443&rep=rep1&type=pdf\n", "\n", " :param entx, enty: coordinate of the entry point in meter (2 components, x,y)\n", " :return: coordinates voxels in x, y and length crossed when leaving the associated voxel\n", " \"\"\"\n", " array_x = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n", " array_x[:] = -1\n", " array_y = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n", " array_y[:] = -1\n", " array_len = numpy.empty(BUFFER_SIZE, dtype=numpy.float32)\n", "\n", " #normalize the input propagation vector\n", " kx = entx - self.poni2\n", " ky = enty - self.poni1\n", " kz = self.dist\n", " n = numpy.sqrt(kx*kx + ky*ky + kz*kz)\n", " kx /= n\n", " ky /= n\n", " kz /= n\n", "\n", " step_X = -1 if kx<0.0 else 1\n", " step_Y = -1 if ky<0.0 else 1\n", "\n", " X = int(entx/self.vox)\n", " Y = int(enty/self.voy)\n", "\n", " if kx>0.0:\n", " t_max_x = ((entx//self.vox+1)*(self.vox)-entx)/ kx\n", " elif kx<0.0:\n", " t_max_x = ((entx//self.vox)*(self.vox)-entx)/ kx\n", " else:\n", " t_max_x = BIG\n", "\n", " if ky>0.0:\n", " t_max_y = ((enty//self.voy+1)*(self.voy)-enty)/ ky\n", " elif ky<0.0:\n", " t_max_y = ((enty//self.voy)*(self.voy)-enty)/ ky\n", " else:\n", " t_max_y = BIG\n", "\n", " #Only one case for z as the ray is travelling in one direction only\n", " t_max_z = self.voz / kz\n", "\n", " t_delta_x = abs(self.vox/kx) if kx!=0 else BIG\n", " t_delta_y = abs(self.voy/ky) if ky!=0 else BIG\n", " t_delta_z = self.voz/kz\n", "\n", " finished = False\n", " last_id = 0\n", " array_x[last_id] = X\n", " array_y[last_id] = Y\n", "\n", " while not finished:\n", " if t_max_x < t_max_y:\n", " if t_max_x < t_max_z:\n", " array_len[last_id] = t_max_x\n", " last_id+=1\n", " X += step_X\n", " array_x[last_id] = X\n", " array_y[last_id] = Y\n", " t_max_x += t_delta_x\n", " else:\n", " array_len[last_id] = t_max_z\n", " last_id+=1\n", " finished = True\n", " else:\n", " if t_max_y < t_max_z:\n", " array_len[last_id] = t_max_y\n", " last_id+=1\n", " Y += step_Y\n", " array_x[last_id] = X\n", " array_y[last_id] = Y \n", " t_max_y += t_delta_y\n", " else:\n", " array_len[last_id] = t_max_z\n", " last_id+=1\n", " finished = True\n", " if last_id>=array_len.size-1:\n", " print(\"resize arrays\")\n", " old_size = len(array_len)\n", " new_size = (old_size//BUFFER_SIZE+1)*BUFFER_SIZE\n", " new_array_x = numpy.empty(new_size, dtype=numpy.int32)\n", " new_array_x[:] = -1\n", " new_array_y = numpy.empty(new_size, dtype=numpy.int32)\n", " new_array_y[:] = -1\n", " new_array_len = numpy.empty(new_size, dtype=numpy.float32)\n", " new_array_x[:old_size] = array_x\n", " new_array_y[:old_size] = array_y\n", " new_array_len[:old_size] = array_len\n", " array_x = new_array_x\n", " array_y = new_array_y\n", " array_len = new_array_len\n", " return array_x[:last_id], array_y[:last_id], array_len[:last_id]\n", "\n", " def one_pixel(self, row, col, sample):\n", " \"\"\"calculate the contribution of one pixel to the sparse matrix and populate it.\n", "\n", " :param row: row index of the pixel of interest\n", " :param col: column index of the pixel of interest\n", " :param sample: Oversampling rate, 10 will thow 10x10 ray per pixel\n", "\n", " :return: the extra number of pixel allocated\n", " \"\"\"\n", " if self.mask[row, col]:\n", " return (numpy.empty(0, dtype=numpy.int32),\n", " numpy.empty(0, dtype=numpy.float32))\n", "\n", " counter = 0\n", " tmp_size = 0\n", " last_buffer_size = BUFFER_SIZE\n", " tmp_idx = numpy.empty(last_buffer_size, dtype=numpy.int32)\n", " tmp_idx[:] = -1\n", " tmp_coef = numpy.zeros(last_buffer_size, dtype=numpy.float32)\n", "\n", " pos = row * self.width + col\n", " start = self.idptr[pos]\n", " for i in range(sample):\n", " posx = (col+1.0*i/sample)*vox\n", " for j in range(sample):\n", " posy = (row+1.0*j/sample)*voy\n", " array_x, array_y, array_len = self.calc_one_ray(posx, posy)\n", "\n", " rem = 1.0\n", " for i in range(array_x.size):\n", " x = array_x[i]\n", " y = array_y[i]\n", " l = array_len[i]\n", " if (x<0) or (y<0) or (y>=self.height) or (x>=self.width):\n", " break\n", " elif (self.mask[y, x]):\n", " continue\n", " idx = x + y*self.width\n", " dos = numpy.exp(-self.mu*l)\n", " value = rem - dos\n", " rem = dos\n", " for j in range(last_buffer_size):\n", " if tmp_size >= last_buffer_size:\n", " #Increase buffer size\n", " new_buffer_size = last_buffer_size + BUFFER_SIZE\n", " new_idx = numpy.empty(new_buffer_size, dtype=numpy.int32)\n", " new_coef = numpy.zeros(new_buffer_size, dtype=numpy.float32)\n", " new_idx[:last_buffer_size] = tmp_idx\n", " new_idx[last_buffer_size:] = -1\n", " new_coef[:last_buffer_size] = tmp_coef\n", " last_buffer_size = new_buffer_size\n", " tmp_idx = new_idx\n", " tmp_coef = new_coef\n", "\n", " if tmp_idx[j] == idx:\n", " tmp_coef[j] += value\n", " break\n", " elif tmp_idx[j] < 0:\n", " tmp_idx[j] = idx\n", " tmp_coef[j] = value\n", " tmp_size +=1\n", " break \n", " return tmp_idx[:tmp_size], tmp_coef[:tmp_size]\n", "\n", " def calc_csr(self, sample):\n", " \"\"\"Calculate the CSR matrix for the whole image\n", " :param sample: Oversampling factor\n", " :return: CSR matrix\n", " \"\"\"\n", " size = self.width * self.height\n", " allocated_size = BLOCK_SIZE\n", " idptr = numpy.zeros(size+1, dtype=numpy.int32) \n", " indices = numpy.zeros(allocated_size, dtype=numpy.int32)\n", " data = numpy.zeros(allocated_size, dtype=numpy.float32)\n", " self.sampled = sample*sample\n", " pos = 0\n", " start = 0\n", " for row in range(self.height):\n", " for col in range(self.width): \n", " line_idx, line_coef = self.one_pixel(row, col, sample)\n", " line_size = line_idx.size\n", " if line_size == 0:\n", " new_size = 0\n", " pos+=1\n", " idptr[pos] = start\n", " continue\n", "\n", " stop = start + line_size\n", " \n", " if stop >= allocated_size:\n", " new_buffer_size = allocated_size + BLOCK_SIZE\n", " new_idx = numpy.zeros(new_buffer_size, dtype=numpy.int32)\n", " new_coef = numpy.zeros(new_buffer_size, dtype=numpy.float32)\n", " new_idx[:allocated_size] = indices\n", " new_coef[:allocated_size] = data\n", " allocated_size = new_buffer_size\n", " indices = new_idx\n", " data = new_coef\n", "\n", " indices[start:stop] = line_idx\n", " data[start:stop] = line_coef\n", " pos+=1\n", " idptr[pos] = stop\n", " start = stop\n", " \n", " last = idptr[-1]\n", " self.data = data\n", " self.indices = indices\n", " self.idptr = idptr\n", " return (self.data[:last]/self.sampled, indices[:last], idptr)\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 4.69 s, sys: 7.11 ms, total: 4.7 s\n", "Wall time: 4.7 s\n" ] }, { "data": { "text/plain": [ "(array([0., 0., 0., ..., 0., 0., 0.], dtype=float32),\n", " array([ 2, 2, 4, ..., 1023180, 1023181, 1023182],\n", " dtype=int32),\n", " array([ 0, 0, 0, ..., 1902581, 1902582, 1902583],\n", " dtype=int32))" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "thick = ThickDetector(vox,voy, thickness=thickness, mu=mu, dist=dist, poni1=poni1, poni2=poni2, mask=mask)\n", "%time thick.calc_csr(1)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 32.1 s, sys: 183 ms, total: 32.3 s\n", "Wall time: 32.1 s\n" ] } ], "source": [ "thick = ThickDetector(vox,voy, thickness=thickness, mu=mu, dist=dist, poni1=poni1, poni2=poni2, mask=mask)\n", "%time pre_csr = thick.calc_csr(8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Validation of the CSR matrix obtained:\n", "\n", "For this we will build a simple 2D image with one pixel in a regular grid and calculate the effect of the transformation calculated previously on it. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAHWCAYAAAA/5CPqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABr7UlEQVR4nO29f8wlR3nn+33OOzOe2ASPjR3HjE3siElY3+hmcUZgxCqXGxNinAgjLWHtJIsX+Wa0eyEhgWwwm6uwIf+E3Vx+ScjJLHYwUcKPOCiMWG+4XgcUbRS8HgJrsA3xxPzwGBsbMDZgjD1znvvH6X7fOnWqqqu6q7urz/l+pFdvd3XVU09Xn66n69dToqoghBBCSH/MxlaAEEIIWXdobAkhhJCeobElhBBCeobGlhBCCOkZGltCCCGkZ2hsCSGEkJ4Z3NiKyGUi8gUROSYi1w6dPyEkD3yXCYlHhlxnKyJbAP4RwM8COA7gdgBXqepdgylBCOkM32VC0hi6Zfs8AMdU9V5VfRLABwBcMbAOhJDu8F0mJIGhje1+APcZ58erMELItOC7TEgCu8ZWwEZEDgE4BABbsuunTtvaN65ChPSJCB576uGvq+rZY6vSB8vv8+6fOm3vWSNrtEE89dTi/+7d4+qxQagA3378Aef7PLSxvR/A+cb5eVXYNqp6GMBhADh999n6gtNeNpx2hAzN7l342NcPf3lsNVrQ+C4D1vt86jP1kuf8qnkREFn+78J3LZQmF7VuoeNYGaaubeTU8V04ZMjxrwEi0P0/1Jw+Ql6yPqn35pLZJv2I6GyGW/7h95zv89DdyLcDOCAiF4rIHgBXAjgysA6EkO7wXR6aKW0aM7auY+fvYNCWraqeEJHXAvgYgC0AN6jqnUPqQAjpDt9lEmSIXoem/Atj8DFbVb0ZwM1D50sIyUurd9mugOtzV8Xc1A3pq8zrcF/aVLnmeUy3spm/7/5smSZtDFWoLOxu7Dma+zS7GkrXfee6r6bnGspjRCNc3AQpQsgaY1Z2TWO2TRVjaCw3VoeU67aRjj1uM2ab8jEQiiuyc73WYYbl81hZMfo0jdXGyG/6WPJh3mdM3Fwt70g918LY1o45pONXSy45Jeq0rnJK1Cnnva0drhZeqNVnV4quFmRMi8aHywiG5IRapCFSWrapNOlrX69btr7ehBy65DBkMTK6GM3QR0HMxC/H7y70xk/e2JoesFQ1W+XdRU4ur1x9yCmljHLJ2YQyWitcLds63Nei8V2zW1Fty9qVn0vHsWcjt+keNVu2tYyYabFTmo2c67mnXnPE0UC6yRtbm7YVnF1pr6ucEnUqTU4JOq0tvhZem4+mPlo0LtmhlvUUx2y7yEtNH9tr4LuW0s3tegaxaQdgsrv+qKq3VZPa2tk0OSXqVJqcMXUiFW0Nsdl1Wh/nNDS2oczZdZpTjk/mPL/4ZB1SrtnxUrvNY9M25duRSRrbmMorNk5TvNiKckg5Oe8/Jg7LqDlOrnvbSFytvlCF6Quz04XkuuLlfkZ85v3h+qBKSZuSR8LvQwJxJmdsUyqtocZOhzY2XfMZI0+WEQ0uIZvM5MZs6/GvpoorZpxMRKIqwCZZJcoBWEYx10sqo43DHGMzw3xxm2S1yd8+7vpBlJo+Zs1r7nS+87ZNr7a65KDL8+r6mwGSxoEnZ2xrQhVcSqU2hJwUWU1GoLR7YxnllbNRuGZ1usqqqfxSy9dl5H3xfPnV13yziuvzpgrdZaTapgtRy7TT5TCSsTLaTJByxfNhPgOXrNCs80ys7WxkVwXXpmKLbZm0lZNLp3WVU6JOY8pZa3wGyVdJxhjf1Hc39pmkLA3xncfMem7zgRH7wWAu/fF9JKTk20WXJvmu30YKTR9tdtwcwzqb0LKtMSu4LhWbXVG2lZVLjktul7R9lFEuOSyjDcI1ackV7gsLhdvXXcalvt7UQgrhq6htAxZbobf9vcTKtuPNAcxalm0XXVINXCiu6/nFys5laF15Bn5Xk5sg1Sd1Jdu18s8tp0RKK6MSy6pEnQiZPBN9rybfsgXyVmoltLI2QU5OWesqZ2PwjXe2keM7D8lMzS/U4rbvg+zQR2uyTSs15XkndrOHlv6shbElhEwY13inWSGGum3tijdX/kD7CVJdDW7Obk5XGY7Vn2k+17b3FxrfTp0g1Ta/AKEJUuxGJoSUh8/L0VitxaZJTDl6M+zJTGP2kLjuLXbSV0iefX8psmImPfmOC+htYsuWEDIcodnIpiHd8lTsrkp7qNnIbYxr11nFoa7vWD3azkaONbiuWeUp+CY7dTWQTeXK2ciEkLUlZTZyHR6qjNs6JoipJLuss22SEYqTcl8xM7NdZdV2pneTHk0ziH3PzRWW2hWc0o08wmxkGltCyLi0NagxlbNJ7ASp2KVFruMcLaY+uspNmWN5fGqzzKspXuxHhH2tpwlSIWhsCSHjkms8zW7VxBq+vsbzck50Isu4urhT1tm2yacjazNBqjRn8SXKKVGn0uSUptNG4GqhpnRzmi2mUCvHvO46zvHMYltZpBu+ZxmbNiUP3+/DIWetdv1xUVoFWcspTa8csIyaKUkXQkhPJLZ6J29szYqtSyWXS05OSrs3ltFwcjaC1BZt7HhcVx2GpO1G7rnSza2/IXXJQZfn5Wq5NvVyxMTxMGlj66rM2lRwueS40uWS01ZWaXJc6VhGG0poaU1T3L502GQmbQ3KZ7LFG6rEcm0MnlpR+uKnjgWWdm9DlVEOOamySpOz8bQ1uLknWcXStVXXtgbOlW5m/LVlTCvS5bmbjjVcfzFpEpjkbOSYyktV0eSjNlYOgKCs2Mq0SadccmJllVZGQ5Z1ik5DldHGoRq/5KapjMfqRu5qaLpuAp+aNrR5fBt5tYzcBjdWZtdu5K5pEt7nybZsCSFrBj9E0ljn2nsN722StxTTOsgZpyleTJyY/HLJGSNOjnsbsqzHiMNWrQXLoywmaQ2mw2SLN1RxpVRqueSE4qdWtKXd21BllENOqqzS5Gw8TWNlqelS8xv6WXHMthsTGrOdrLEF3JVYm4otlxxXulxy2soqTY4rHcuIEFIUMeO5Qy39EZHzReTjInKXiNwpIq+rws8UkVtE5J7q/xlVuIjIu0TkmIjcISIXt83b0sN5PJacnJR2byyj4eQMzSjvc8w62/rcvGYfu2R10cHH3HPcxNz4b69nbTOb2bc2NkbW3PG/61pbO23XGdopenR57qlrbO10ifl3admeAPAGVb0IwCUAXiMiFwG4FsCtqnoAwK3VOQC8FMCB6u8QgOs65L1Erkott5zS9MoBy6iZknRJYPj32dWF6+vWNa/Zxy5ZXXTwMfMcNzEz/ttdtm1q4Bnay5o5/vvktdEnVo8YeTF0ee45upET8m+99EdVHwDwQHX8bRG5G8B+AFcAeFEV7UYAnwDwxir8fbpYA/FJEdknIudWcjpTWqW9rnJyylpXObllDcEg77Mq5Imn0pWbCTAfaa2yuTzJd9yUHljENVtBbeWEsPN46gQwE8j3T8TLsPUz8zfDYmWlvAd1Hh2W14zOzK9rlnW2InIBgOcCuA3AOcYL9yCAc6rj/QDuM5Idr8KyGFtCSB76ep/1ie/j5D1fzK4vdA7IbOe/HVZjn/vCXMTGK4naaD362Lh6lIjveUrXZrmfzsZWRJ4G4C8B/IaqPmZ+0auqikjSJ6mIHMKiWwp7Z0/rqh4hJIFe32ecCsxP5lR3Bz25/N8+dp37wkLypwg9mK3ie549PudOZlxEdmPxYv6Zqn64Cv6aiJxbXT8XwENV+P0AzjeSn1eFLaGqh1X1oKoe3DPb20U9QkgCfb/Pu3FKf8oTUjitW7ay+OS9HsDdqvo249IRAFcD+IPq/0eM8NeKyAcAPB/Ao7nGawkh3RjifZZT9mDX/h+xM25uefnixKTdYOZf/yYAYHbWmTuB5hgyyY8I8E/uS126kV8I4F8D+KyIfKYK+w9YvJQfEpFrAHwZwCurazcDuBzAMQCPA3h1h7wJIXnp/32ezTB/+qnpmuU2tl0nJZn5+iY6+Sb7xJDpI0Ie2w2IYP6DP1DGh0kOHULPzZ7AlpI2EzrzdxZ3mY38PwD4tL/UEV8BvKZtfoSQ/hj1fQ5VknYF7YrbxaF8GwNg6+M6jtXL5XC/bToXdhnZGxiYc4T6dnGUw9h32ZDCng2eQx/bQU9A5qQ9SNWkbmHXJKs0OTlklSanllWanNLKiFQ0ORvoIpd0Y2KTtMdiLYxtLuoKsmtFmVtOiZRWRiWWVYk6FUVMS7ApbKwybuNBqk3akLy2Hp98OnTRq4suPnlrxuSNrVmhdanc7LRtZeWS0yS3bdqcZZRLDstog/C5xOvSldsnXd0Q5nLVGHKnaOsYihc6t2U1pXXF7eL6MeRKMje5fj+xbh4xcWPrqtjaVHZ9V/a5dFpXOSXqNKactcbn5q7k2bE5dsRxyUh1+RiKH3IJ6cuvjcvJVH1SaNiRSE4W/i41/IYna2xDlVjKmNkQcmKux+aZUnmzjIaR0xSXY7g9MJRxDvkeTm1x5drOztYph2FrKy9khDNvH6hbGZ95rt9Pgo/kLO4ahyS1Agz5qE2plNdRTok6lSgHCPs6zqXTRuDrPm7Tjdy2LGNn8tozde2ZvCkyQobZvjZzhMWmXbqmlV9ph0xfulzjyS7M8oq5Z1ecOtxV9nZ4aKZ1zm5kk8BvcnLGVkSiK7ehHNTH6pRLTkw+LKPmfEoro43AtxZ1yHWgsS2u2gB02aXHZzibDEebVrQts3aKb34cxBrT1I+RVDmhcJ+RTZVZWL/t5IwtEFdRxlRsdZyQrNgKskmn0uTEymIZDVtGG4Vr3DampWCG9W2kfYa2S1drVz2AcOs8Zau9JsMWKyPnPQ5hJHtaZxtiksYWCFdwqRWbr9JdVzkl6lSanDF1Wmt8FVwo3LXdWxdijENTd6jLOYRtwGJbkm27kZviDNmNHNvK7dpNHXpurmfgS8tu5O60rdjsinJd5ZSoU2lyStBpbcnVjRxyy9dESssvNk7o2GeIgfhuZJPY1mhTN3JMl20MKV6oUrqn23TZx+aTs0ck8h0vrFc7HbMyK2VsrRQ9XHJK1q0UPUq5N4K4buQp0XU2ci3DddxWRg5yypuSVdqEbmSTEivH0nRaVzk5ZZUmZ+2JGbMtlaHHbIeU25ZNNbiRrOEtEUImgW8ZUExY2yVDXeji/an+n8MjUl/uGkPhMfK6eNryyeqLvsZsA6xFy5YQMkFSZiObY7Q5ZyOnpO86GzlX08a1JKmNjJiwNvJydpP3BWcjE0I2mqYJU6Ex3KbK047nMtpDGm9g+DFX1ySprozZP9rleTUZytjfn5lkk2YjE0ImTKh72GUkfXFDsl3nubqlU9O3baF2Sec7b2s0c62zbUOX59UmbcNSn7Xfz5YQQggpGbZsCSHj4lvKE+swwNfStbsYfeO+Nrm6kWO6ONmN3J7QbyU1fQhfj4orKruRCSHFEePUomn8NlTJxnQbp2A7p3AZGVceLg9OXVwdNnlHiknr0ilVlksfc/JW167uprRjdyHXmGv0A3LXwtjG7MwypJxaVi45QDn3xjIaTs5akuqusQ/a7vqTktaWERorDRmXlN1/fMbTdy20/Md2SWkb0KalQ7ZRb6N7Du9WQ+DqPXFQqvrRmK7x+t4EPjV9LjldYRmlySmljNYSkZ0/O3woYms915Kdti4EXTvSzBzHrvQpXqdCOrfR3U4bo7NLn5g4Lpmp9z8mDb/hKdxCEm0rODvdusiJkd023brK6UKfz20tqB1SxDi0KIUuDhZCjh5SWq2p+9t2yc8V39XqzO14wiVzCAcXuWj4DU/W2KqqtyJLreDWWU6JOk1Fzlg6kYLJVfHbhqyvzdy7yojteu6Sd644hTNJYxtTecXGaYoXW1EOKSfn/cfEYRk1x8l1bxvH0JNcgPjKvck9o+9629aY7fKwq2Gbdyhb1721dcfoui/bdWUOg+t6BrlcZGZgchOkUiqtnBNwQnJSjEQOOTH5pMRlGTXHLWlC2KTJtcVeLavN7yFlHDE0Azl1i72UCVdtJkvZMpu22IuRFSJGlmsMNiSvy7Z/OV1G9kCBKoURkegKqyneOsspUaepyhlKp41g6mO2bVubobHImBZblzHbLi1S13GKrBzdzgW0SnMwuZZtjWvDb/NaDjkpsnLKAfwtr1z3VpqcFFnrXNYbR8jI1i2d1OVC9pIVGOem3FC+obCQ8bKXzMTSx9jsXHdat8Y1UQ06X1iK3+VDo2u6Jhmzhus+RmpiTtbYAu4Krk3F5qsoU2XlkuOTVZqcNrCsN5yUbuSZkcb12/OFN3UnNlW2diUe06WcIt+ME7MOtY6XIhNY7UauwhsNrZkml04xdOlGnkAf7QRUDNPXht9t5eaS0yS3bdpccrrIKk2OnTZnGZENp+0aV5+cLvFz1/ZDyJuShep7na2IbInIp0Xko9X5hSJym4gcE5EPisieKvyU6vxYdf2Crnk7dMmSvjQ5JVJaGZVYViXq1ERJ7zOxCNXWbSYHdWmZt5U9tIyhiHjXc9zO6wDcbZy/FcDbVfXZAB4BcE0Vfg2AR6rwt1fxspAyIShGVmlycsgqTU4tqzQ5pZXRCIz+PhfFlCr8MUiZbbyuRL7rnYpGRM4D8PMA3lOdC4CfAXBTFeVGAC+vjq+ozlFdv1QmXCMRsm4U+T67RIrh7nHIKiRHV3BX14N2+rbyhmjVtpU1hMHO9btJkNP1tt4B4LexM6T9DADfUtUT1flxAPur4/0A7gOA6vqjVXxCSBm8A6W9z64JUGZYyUuGXORwsFDKDOG+KFm3DrQ2tiLyCwAeUtVPZdQHInJIRI6KyNEn50/kFE0I8TDI+3zi8dV1tvb/WLoY2VSfxF08Efk8QoW8KoXy9i1LCsULrfWNlRWrT1ud4YibWu4hXWxyfaSZv+cGmV2W/rwQwMtE5HIAewE8HcA7AewTkV3V1+55AO6v4t8P4HwAx0VkF4DTAXxjVXc9DOAwAJy+++yJfbYSMln6f593na3ypa/2fiNkwfx7i8bK7KkTDTFJNmb+buXWxlZV3wTgTQAgIi8C8Fuq+ssi8hcAXgHgAwCuBvCRKsmR6vzvq+t/o3QYS0gRDPE+68mTOPnYY73oT/ycfOrJsVUg6Gco+o0AXi8ix7AYw7m+Cr8ewDOq8NcDuLaHvAkheeH7TEgGsniQUtVPAPhEdXwvgOc54jwB4Bdz5NcnpTmdL9F5fWn3ts5lNAZ9vc+yezd2nfNMKzBiM4FUD1IEADB/5FuACGb7TrcuzIHZJq7RGYDZDPiK+9Kk3TXW1L1XXSu33DvKlGRMWEbxOpVSRmvHri3MzzIq/thdf8Y0tm3cNbpk2PTtrnEGyBOLMdv5M56eJseUlUunGLq4ayyFgLGd0m04MSu2LpVcLjk5Ke3eWEbDySFrRt/uGsfYXi7Wo1VbeV3XJA+Jufbbw6Rbtq7KrE0LJ5ccl6xcctrKKk2OSxbLaIPwba3X5sMkZtcf3/WUXX/s5ShNLV3XrkNNMrviXF5T7fpjXUva9aerPqaclP10Qzo0PddQHn0Z7/q3GCjXqXw3rBBqNaRuDJ5DTii+qo6iU2lyQvFzyUmVVZockoCvYsvtpSkUr6/8c2PIjNr1p0+6bDjQ9rkU0EqepLGNqbxyxmmKF2tMh5IzRpwc9zZkWY8Rh0Y3M32VZ1snCqnpR0LG/h36nFy0TZtTfo9M0tgSQgjZQCZssSapeswYWM44TfFix+Ri5Ax9b7nisIya43DsNjN9lefY3cg9M+lu5Fj5Be5GNNkJUiLi7ZZLqdRyygHcXYqplWyJ98YyGkbO2mOWRezSn9ykjA/GLvcJhbsm7jTJdIU1dYXaMmvXgeYuPDmX/qToZJ+H7tlVTjEUaGBNClMnDVcl1qZiyyWnS7oYOWPeW59llLOs16WM1hqXUR17HDFEl/E+Y1buylhp6oYIifl5r8USmjHdx7h0aGODNWDSxhZYrsy6VGy55ITkdkm7TnJyUtq9lVhGRWEv9RnDyHaZjNOh4l/qvm3a3SZlt5sm3dro7NqJp61Osbv+xOiTEqeNvj0y2W5kk1AX3phy+mzljgXLqJlcZbSWTLkbOSWtS8bQ3ch2eAndyKHwqXcjN9RBk2/Z1pRWaZcop0SdSpNTmk5kZEqqIWM9NnVds9rFwLWRkSq7NKTZexQwndshhBDS90zerozhNrINI3wMl1wchBDiZoyeA9aWpANrMWZLCJkgtsGsz82x2zrMvGam6zrW2zZ9m67WXMbaNZbcRkZMWBt5fW+6UBIJH31Tui1CyDrh2pQgtFGBawZz10lVbdO3mSWca3asa6ZwWxl2WBd5oWVHbWWVTsLvhy1bQshwuIzpHMDMCLdno6p1zRfWpoJ2zRLO6afXFzfnrj/eJT8Nu/7kXCpkp+1rJyOb1s89wki65Lqap5EGly1bQsi4pOxj2tT9GerONOXaeeRyscgadX1w/UaaCHQrr0XLtl7TmGO5Ra61nyXKAVhGTXKA7mWUs6zXjpR1tjMjjav1YIe7xkZjxiZD5z53jSmuG9uss3XR1V1jdS3KN3Lf62x9snzllCuPNmSSz+8wg7qS7OqQoDQ5OSnt3koso5oSdVobhvqQKWnij6/13aYF5ouX2sIP+T5u4wAkdJ6TnL+fep3tuju1MCu0LpWbnbatrFxymuS2TZuzjHLJYRltEPVEJ99EqBQ5Q9B14o9vMlOK3CbXhk0uFX35xMjyxfXpk9O1Y5/uFnP+fny/aYtJG1tXxdamsuu7ss+l07rKKVGnMeVsFFPxjdxkGHPNwM3pG3neoWxz+Ua2ZdhhKb6RYz5KfHkUwGSNbagSS6ngmuTEymqKlyJnqHvLJYdllK+MSAK+Mu1awaYYgL7yz40hc2X3oSEJtc5j07fJowCjO7kJUqkVe2iSSoqsdZUDsIya5ADDlNHGYk8SCi39CaWP6Sa1r7WZkOPrXk3JfxPJVR6hpVqxrd/U52zmFyLwW52csU3ZUWVqTuXH2L2IZRQXNwc0tB6aJsY0eXmyZyE3Vai2cW+7k4+dpyt/GtwdcpWH7/fi++gKpU3JJ4bAOz7JbuSYSis2TlO8nHJy6J1LTkocllFznBxlRDLRdc1srjW3XfLvUWbU0p9S6bIUaGRrN7mWbU2oZZJasflkraucEnUqTc6YOq01Ofez7duvse2DONQiTm0pN62zda3VTZEJuNfZxtKUJraFmrKcx9cybbOUqcBmZIEqxeOqxNpWbHa6dZVTok6lySlBp7Ul19KfIelrwlLK0p8UuV3S+9LlkNE2/pp0xU/a2ALLlVkpY2ul6OGSU7JupehRyr2tJT4HACWXVV/dui6PV13zD7WSuziZaOtwoqtjizYyCqXTbYjIPhG5SUQ+LyJ3i8gLRORMEblFRO6p/p9RxRUReZeIHBORO0Tk4jy3ED9GN5ScnLIoZzhZpckZmlLe57WiraHzyekrfhty5zF1Q9vwzne9lXcC+GtVfQ6AnwRwN4BrAdyqqgcA3FqdA8BLARyo/g4BuK5j3oSQvPB9dpFjzWzX7tHUdah9d8f2sW516l3IfXmQEpHTAfw0gOsX+eiTqvotAFcAuLGKdiOAl1fHVwB4ny74JIB9InJu2/wJIfmY3Ps8VM9BSnevL70tx76WokeTP+OmeE35x87E9uWT2zdyXy3bmN9PrKH3DY1YdLmVCwE8DOBPROTTIvIeETkNwDmq+kAV50EA51TH+wHcZ6Q/XoURQsZnWu9zyROqfEyhlTYFHYcis6HvIm4XgIsBXKeqzwXwXex0MQEAdLEOIumtEJFDInJURI4+OX+ig3qEkAT6f59PPJ5N2aInVPmYwvhjDh2ncJ+5SPgddimW4wCOq+pt1flNWLysX6u7k6r/D1XX7wdwvpH+vCpsCVU9rKoHVfXgntneDuoRQhLo/33edWpvyhNSOq2Nrao+COA+EfnxKuhSAHcBOALg6irsagAfqY6PAHhVNYvxEgCPGt1ThJAR4fs8AFPoos2h4xTuMxcJwxldPUj9GoA/E5E9AO4F8GosDPiHROQaAF8G8Moq7s0ALgdwDMDjVVxCSDn0+z4/dQLy1a/3oHbh6BwQT7smdK2tzIr5d74LAJg98I3ldMBOWjUsYxs9InUhHY2tqn4GwEHHpUsdcRXAa7rkF9ADQHcnArnk1LJyyQHKuTeW0XByhqbv91lPnMDJhx9upxxpDcu8DCb/OWL6os21Z2hXOXX6XHK6wjJKk1NKGRFC1ofJbkTgo21rya4Y10WOT/Y63FtfcrrQ53ObOiKC2V5OehyK+ZNPAQBme3aPrMmG8T138KSNra+STK3g1lVOSFYqpd1baXJyy1pL9p4C/NgFO+fmLi+hfUh91wbYL1ZUt7ek8x3HyhCz96SFnDq+C5eMrfu/BswEeu4POdPbOoVkpepUy26znV8tb3JbAc6wWAr0GfflSRrbGAMSU8HlkhMjK6ccoHk8cCg5dRyWURgaXAe2ZybfoJYIMPOUcZuBMN+2dnYcLFf428czQOfW9nWBzePrdC7j4TQogY+I1gbI0qWzvAgZK+GRH0dBnWK26KuPU9LmIqD75MZsU1pqQ429xeaTS07XfMbIk2XEMdzJsI5LV8xJx1P+HcY8G1+ckZ/r5IxtSuugKW6srHWVU6JOpckZUifiIdU42HvF9lnJDuWQn0yeSXYj15WWr6WQWpGGWhwpFWkuOUD3exuqjFjWzTrRyAYQaTamvji+8NA+rL4w33WrSzhKTtsmjCtdiuG1dZyJ81pU17Epyx5XbzNWHhoqiJUZ+9wKbUIWqlYcrkqsTcXmS5MqK5ccX5rS5LSBZb3hiEx78/hcNaZrp56meLFyY8JiZcWOq8fICMVxHafImACTv4W+jEBbubnkNMltmzaXnC6ySpNjp81ZRmTDKWnz+Jy1fR8GcM03j59kN7KPrhVd3RVYmpwSKa2MSjRyJepUFK7Wrf17d5WhGRbTDd0FX6urTWszp8GdO8Jc8ZpkxXaNh+TEzO5Oldk3uX43Ce/4lL4bvIhIES3IvuTkkFWanFpWaXJKK6O1xq7sXJWfL6wOH/JjdO45jk03R54JVW318MmICY+RZ95nF/qe1JaThN/fWhhbQsiGMcbHzNBdwWNQwn62Uyinmk1r2RJC1pSmymxqPQg5xjpzTNhizT84LHJCCCGkZ2hsCSFlENtKLaE123U5zrqxCffYERYRIYQQ0jM0toSQcWha1mOH99GiHdu5Rg6HE13ym1l/bXVxye7C0JYp9ffV4ndCY2tRmtP5EtfZlnZv61xGa41r6Y+v3Ppa7hOz/MhFm6U/XeXY6dqktdPM0X1JUq4lTbbMIQn99nzxE1kLY6uqWSq33JV/ScaEZRQno6QyWntSWpVmyyNn67OtrBLGbIfw4JSSdsot2wGY/C2ZFVuXSi6XnJyUdm8so+HkbAQpTi3slm2uVm5K+tKcWnRtjbrCNsWpRa53M0HOpN01uiqzXK772sqxdcolp62s0uS4ZLGMNgifgW1T+fnShFwZ1tebNiEPXfNtGO/ZPN4pM2S4U3fV8XqD0sXOP9Z1UV3s/NPHvq8hg5vrvpqeayiPPpqX5u8w8K5P1tiGWg0pFVwuOSFZueSkyipNTkgWy2hDMMtBdcdHbRtftTFb7LlI3WLPZVDteKFjlxHwyQzpGLsFXS2z3mLPup/oLfZiSN0Wr+kjx1dOMXCLvfzEdM/ljNMUL1ecFJ1KizNUGa1zOW40JZdPXy29XJOs+sivKe1Q/otz5FGIr+VJGtt1hRXyesPni53xVpZFPGMZijHyzf2RABSzH24BKqQT0yWXM05TvNguwhg5Q99brjgso+Y47EqGf/P4TSS29h1r5nMOH8xd8hw6756Z7G2EKq6USi2nHF/81Iq2xHvLJYdlRJyUXEZdl8R0ldEm/RCbx5MkJl1krkqsTcWWS06XdDFyxry3PssoZ1mvSxkRksQYNXkozxzerdrKGYuGd31Kt+LErMy6VGx9VYq5dFonOX1Rwr2VXkYkkZJqSN/YY8qM35h4XbqtN7H1HDksMpXbCZK7NdlVXi45trwSWPeyzkFJuqwtY5Rxjq7g0vezzbl5fMl77eb6/STI6XRbIvKbInKniHxORN4vIntF5EIRuU1EjonIB0VkTxX3lOr8WHX9gi55O3ShnAY5JepUmpzSdBqSwd9nu4xcZdYUNlY5013jatpNbNUm0PqWRGQ/gF8HcFBVfwLAFoArAbwVwNtV9dkAHgFwTZXkGgCPVOFvr+IRQgqgmPfZN1s5NIs5dnazGc+Wl2uWdGr6oVt/oV1/2jKmYezyvMxn7vtrSmMRchjStZh2AfgBEdkF4FQADwD4GQA3VddvBPDy6viK6hzV9Utlip//hKwvZbzPoXW4Pv/JMet2ff6VU2TE5EGGo0t5u9Z8m2FNvzXHdQno09rYqur9AP4QwFeweCkfBfApAN9S1RNVtOMA9lfH+wHcV6U9UcV/Rtv8CSH5mPz7zO920ga7hdrUu9Hhd9alG/kMLL5uLwTwTACnAbistSY7cg+JyFEROfrk/Imu4gghEQzyPp94vDlB25bKCLu4ZCHnrj1t0uXYiagAV4itsVuyTT0cHTbO6NKN/GIAX1TVh1X1KQAfBvBCAPuqbigAOA/A/dXx/QDOB4Dq+ukAvmELVdXDqnpQVQ/ume3toB4hJIH+3+ddpzZr0bblMMLs0ixwzLYbfY7ZxqRJoEsxfQXAJSJyajVWcymAuwB8HMArqjhXA/hIdXykOkd1/W+UzmIJKQW+z4TUxM4BSKDLmO1tWEyM+AcAn61kHQbwRgCvF5FjWIzhXF8luR7AM6rw1wO4tm3eDl2yOXlfZzkl6lSanByycpb1UJT0PhMSjec9k+892TptX+k67Werqm8G8GYr+F4Az3PEfQLAL3bJr2/qCrLr/qOlyclJafdWYhlNlXV7n8kGs2d3XLzacJp1R08fypPdPL7GbEF0qXD7amXlMgK57o1l1KxTCXLWkvkc8t0Wkx5F3BWgLzwjMldotQm77zhWhsyN30YLOXV8F04ZTy0mkcvj329MHyUvUZ/Ue7Nltkm7JKdT6pbM/J3Fkza2rsq/TQXXp5wSdSpNTok6jSlnXdHvP4mTx744thqbx2OPja0BwYSdYoVaWSktsCY5sbKa4qXIGerecslhGeUrI0LIejK5lm1qxR5qVaTIWlc5AMuoSQ4wTBltBLOtsTXYHOYnF/9Z5sNy0h08OWMrItGV29ScyqfcWy45LKO4uDnYdEMre0/B1rOfvROgujPuGhp/dV2ry3LqPQb1vbdJB4TTPvDw4v+5Z/vT+co1VZfcv+0+ZJrkHO839FSRxVx+B5MztkBcRRlTsdVxQrKGlFPHG0JOrCyW0bD3ttaIQHcZI1fb5SXhSVBLcbFaCU/Z4M7RbjCv9toUSCvVBKOlMp8v/PeqyCKt7f2prS65ByT7kGmSy9jav8UeNyIYjVDllVqx+eKvq5wSdSpNzpg6bQy2F54Yzz0hWW11yCGnLV08QbVJO8OOoTXldPEi1YcVGcoy2R6hYn9nLX4nkzW2gLsSa1ux2elKk9OWTSqjnGXNruOe6KsVOuXW7RBM2X9xn6gu/3bsc1f8lkza2ALLldm6VZC57m1TyiiXnNLKeq3oq0xY1mEmX9N76PLcXS3ZptateS0x70mO2drkrHBzTr7JoVdpxrHEMqrllSAjp5yNpmkMtynMJxNwT77J9cxyjQWS8WnxLPvcPH7tYIXbTGlltM5lvbGMYbBoJIlJ5t/DWrRsCSETxtfKTK3s7HRmy6SnSS+N+tCA94Pd9QvEl3VqT0gm2LIlhJAUQuN5xE/OrvoJQmNLCCmb0ivX0vXbJJqexYjPisaWELLe9FnB9iGbxnstobElhKw3Uxs3nZq+JAoaW0JI2ZRsfJqcIJBhaXoWIz6rtTC2ubYwy7kVWk45Jd0by2g4OYQQBxN9t9bC2NaUYgTq9LnkdMWUwzLqV07fMklLcj4Lnyw+7zAbXj6TN7Z2hda2gltXOTGy26ajnH5lEUI8mN33qV35qe9kpqGCSRtbX0WWWsGtq5y2aVLkTL2MhihrGtyOdJ2d27STSxt5LtmcRRwm5zpbu9xjn0Gb30Km388kPUjFVF4xfndzyYmRlVMO0OyicCg5dRyWUZicfqDXDttXsV2e5rnLr3Ed5kvnej4uL1P1eYy+Tce+fDeZXOVhP6/QM/ClTX3OZn4BJJD/5Fq2Ka2EocbzYvPJJadrPmPkyTJiC5eQogi9jzHXEt/nyRnblNZBU9xYWesqp0SdSpMzpE4E3SvAPvPvg9R9ZufGX678przXbdfnVY/H2n8Ws+88sZxfi3wn2Y1cV1q+lkJqRRpqcaRUpKXJAVhGTXKAsspo4wl1Jzf5JG6qAGO32OtSgaduPpDa3OnaPJph1bhOrsll0GWzh4Su5PnT9kalWdst9lwVWJtKzZcmVVZpcnxpWEbNacYso40j98SZVLl9PR8+935xTZBKSZcat+PznLSxBfJVZLactnJLk5OT0u6trzLqIqfE5zZ5cpRhzudgV8CxM2HbyF8n+vqoSsmjjQ6pxtzD5I2tSQmVrZk+l5yumHJYRv3K6VvmRtG24rQNYOpziE0Tk38pmPqMWev7yiZnL0RbQ9wzkxyztSmxwi1Np9Lk5JS1rnI2gi5jbn1j6takZ9dnnmMM2aVjW71qWb5x7pRn1sW4pqTJ2OXbB43fOCJyg4g8JCKfM8LOFJFbROSe6v8ZVbiIyLtE5JiI3CEiFxtprq7i3yMiV/dzO4SQEKO/z7m7XNvQJe8+ukJzdn26ZLYt85Dx6ms83GeYCzSeqcR0KLwXwGVW2LUAblXVAwBurc4B4KUADlR/hwBcByxeZgBvBvB8AM8D8Ob6hSaEDMp7Mfb7nKulMxQ5jHPulm8JDKFTiffdkkZjq6p/C+CbVvAVAG6sjm8E8HIj/H264JMA9onIuQB+DsAtqvpNVX0EwC1YfeEJIT3D93lNWQejtA73EKDtUPk5qvpAdfwggHOq4/0A7jPiHa/CfOGEkPHh+zwFcnfd9jmRaF0NZ4dx+87z0nSxkj/bDAcROSQiR0Xk6JPzJ3KJJYRE0Ov7fOLxXGKHI+ekm76MWF/LmtrEz7kUqk+D3WXimCssQl5bY/u1qjsJ1f+HqvD7AZxvxDuvCvOFr6Cqh1X1oKoe3DPb21I9QkgCw7zPu07NrnjvpDi5H0JWaFOFHORyf5hDVp+MUGZtje0RAPUMxKsBfMQIf1U1i/ESAI9W3VMfA/ASETmjmkjxkiqsOEpzOl+i8/rS7m2dy2gg1vZ9LpbSulnb+kfu83c+pXfI41PZpHGdrYi8H8CLAJwlIsexmIX4BwA+JCLXAPgygFdW0W8GcDmAYwAeB/DqhR76TRH5fQC3V/Heoqr2JI3WxG6FFiunK7WcXNuq5ZDDMorXqZQy6oMpvM/EQ9NGDL41sTHyZlZ413ey4HdgLBqNrape5bl0qSOuAniNR84NAG5I0i4Cs2LrUuHmkpOT0u6NZTScnL4Y/X327fk6ZOWcaozaGDKXDFe4fWzmk4pPZq1zqkyXg4xUWaZTEDs8xmmIq2xCefkcbxTwHk7aXaOrBTF216SdLpectrJKk+NKxzIig9KXQ4Yu2JNs6vM2Di1iw9vq1UaGLyz2HtvEaVOGPTJZYxuqxFI3BvfFT60o+5aTKqs0OaH4LCMyGKmtsikTMZY4OjH6lX4PEUzS2MZUXjnjNMXLFSdFp9LiDFVG61yOhGSnkFZdkCnomAEpuRIQkW8D+MLYejg4C8DXx1bCAfVKoxS9fkRVzx5bib7h+5wM9UqjFL2c73Ppu/58QVUPjq2EjYgcpV7xUC9Swfc5AeqVRql61UyyG5kQQgiZEjS2hBBCSM+UbmwPj62AB+qVBvUiQLnlTb3SoF4tKHqCFCGEELIOlN6yJYQQQiZPscZWRC4TkS+IyDERuXbgvM8XkY+LyF0icqeIvK4KP1NEbhGRe6r/Z1ThIiLvqnS9Q0Qu7lG3LRH5tIh8tDq/UERuq/L+oIjsqcJPqc6PVdcv6FGnfSJyk4h8XkTuFpEXFFJWv1k9v8+JyPtFZG8J5bVp8F0O6sf3OV6vSb/PRRpbEdkC8G4ALwVwEYCrROSiAVU4AeANqnoRgEsAvKbK/1oAt6rqAQC3Vueo9DxQ/R0CcF2Pur0OwN3G+VsBvF1Vnw3gEQDXVOHXAHikCn97Fa8v3gngr1X1OQB+stJv1LISkf0Afh3AQVX9CQBbAK5EGeW1MfBdboTvcwRr8T7XHnlK+gPwAgAfM87fBOBNI+rzEQA/i8WC/HOrsHOxWDcIAH8M4Coj/na8zHqch8UP/WcAfBSAYLGIe5ddblhsefaC6nhXFU960Ol0AF+0ZRdQVvsB3AfgzOr+Pwrg58Yur03747sc1IXvc7xek3+fi2zZYqdga45XYYNTdT88F8BtAM7RxX6eAPAggHOq46H0fQeA38bO7pPPAPAtVT3hyHdbp+r6o1X83FwI4GEAf1J1h71HRE7DyGWlqvcD+EMAXwHwABb3/ymMX16bBt9lP+8A3+co1uF9LtXYFoGIPA3AXwL4DVV9zLymi0+mwaZyi8gvAHhIVT81VJ6R7AJwMYDrVPW5AL6LnS4mAMOXFQBUY0pXYFF5PBPAaQAuG1IHUg4lvcuVPnyfE1iH97lUY3s/gPON8/OqsMEQkd1YvJx/pqofroK/JiLnVtfPBfBQFT6Evi8E8DIR+RKAD2DR9fROAPtEpHa7aea7rVN1/XQA38isE7D4mjyuqrdV5zdh8bKOWVYA8GIAX1TVh1X1KQAfxqIMxy6vTYPvshu+z2lM/n0u1djeDuBANdNsDxYD4UeGylxEBMD1AO5W1bcZl44AuLo6vhqL8Z86/FXVzLxLADxqdLlkQVXfpKrnqeoFWJTH36jqLwP4OIBXeHSqdX1FFT/716iqPgjgPhH58SroUgB3YcSyqvgKgEtE5NTqedZ6jVpeGwjfZQd8n5OZ/vs85oBx6A/A5QD+EcA/AfidgfP+F1h0k9wB4DPV3+VY9PnfCuAeAP8dwJlVfMFixuU/AfgsFjPm+tTvRQA+Wh3/KID/CeAYgL8AcEoVvrc6P1Zd/9Ee9fnnAI5W5fVXAM4ooawA/B6AzwP4HIA/BXBKCeW1aX98lxt15Pscp9ek32d6kCKEEEJ6ptRuZEIIIWRtoLElhBBCeobGlhBCCOkZGltCCCGkZ2hsCSGEkJ6hsSWEEEJ6hsaWEEII6RkaW0IIIaRnaGwJIYSQnqGxJYQQQnqGxpYQQgjpGRpbQgghpGdobAkhhJCeobElhBBCeobGlhBCCOkZGltCCCGkZ2hsCSGEkJ6hsSWEEEJ6hsaWEEII6ZkoYysiN4jIQyLyOc91EZF3icgxEblDRC7OqyYhpDRE5DIR+UL13l87tj6ElExsy/a9AC4LXH8pgAPV3yEA13VTixBSMiKyBeDdWLz7FwG4SkQuGlcrQsolytiq6t8C+GYgyhUA3qcLPglgn4icm0NBQkiRPA/AMVW9V1WfBPABLOoBQoiDXGO2+wHcZ5wfr8IIIesJ33lCEtg1dIYicgiLrmacdtppP/Wc5zxnaBUImRyf+tSnvq6qZ4+tRyrm+74lu3/qtL1njaxRJlQBkfZpa3LI8PHUCUAA7N69mq5tvk065JA7cR773gPOdzWXsb0fwPnG+XlV2AqqehjAYQA4ePCgHj16NJMKhKwvIvLlsXWwiHrnzff99FOfqc//3w4tLoi4K+s5IKrQutKu+t7kZBXm6oub7xy60poypcpTTaNgy5zvyFqhNrIeY6tWmFfGdgRZSdOEqEYZWzn+tYX8/T+0o2t9/7PVglzRNcagrwhZvpeVcp4jClsXXxk5yxdo/BhKLfMgZlHOgVv+4fec72qubuQjAF5VzUq+BMCjqvpAJtmEkPK4HcABEblQRPYAuBKLeqAz3oowpbay4matXPtmHRdkzoz/ifeX/OwiP0aykHAvUS1bEXk/gBcBOEtEjgN4M4DdAKCqfwTgZgCXAzgG4HEAr07RlxAyLVT1hIi8FsDHAGwBuEFV72xMWFecdgVat0Rmumj9zIz4ZmvVlXa2U7Hq3IqnutSiWqm4Z5a8Or4tqxZpyHYaAdvIh2R40njDDJmS0h1s6Lod2yHf1jXaxBmtSGf5umgwUq5yc6Uz4y2ViadcnPqZvzfA3/qOMayBOFHGVlWvariuAF4TI4sQsh6o6s1YfGinJHIf2+d1hbcVIc9XOXpaN0tdzXMsGeulvK24YukucFTeZsVtdEebXditx2xd3dsxLbhK16X4ro8A1/Mwu5/tjxI7LpYNdOveBHsowboW1LmBlWGEWl5M97YZp0Xvw+ATpAghJJqmMU8LZyVtGc9UlmRaRseU5zIO9fVgvq5KvDbakWOcvZFq2A1EddHybBhnB+C8V+f4einUuiYY3XUcHSCETJW52zBJqAXrkWOnN+V6jZ/HMMe0JrfjzFfzj8ZMO7ahzYV9T677Ctyr/ey8mD0H5p9HZhZsvQP3QWNLCCmbuefYQ2xFmq3CJXF0/HhIMrhDEvlxxW5kQkhROCexnER008C1RMgedw3OeG6qOH3LlnJRQvdxH7S8r+hu5ITu5mxd07lnIxNCSBa2x7ock6EqVsZd7W5Vx8xRu8t3ZWatY3LT9mkorovUMWM7XWjiEdBskFzdo9Y62pVzz7rTle5xO10Msc4tWq6x7aNF29bUhuYDNEFjSwgZHsfY6PaMXcMYbi9zscLstIsLhsG18/PNfDZnCYfiN4VvZyxR47tLuuZoJbvub67AzHDAYRhp58xo+7yLXo7ZyW3SB7F7GCKeTZJ8l4hWqRbQ2BJChsez9lVFliu0WWVwrXh1+qUlNZYck6Uzq5Xn9Prkanm6Wn6Oa0td1mIZOxvP+lS7hdzKIM88a5od+W/fQ+iaa210iEjvWCtLo1z5x5ASv2M3cq1vSnc0J0gRQkiIkFGZwiSrtjpO4d4mBFu2hJDhMD1DbcFZoS+tzay8Sq2EAdvpFeJet2nLrJPD0Ur2xF1aU2vrafoXdsir5dT5OdtAItCt1StqxTbP5aRLkAcRwOEHeXsSmdmdb1wTuyXrasG5Wut2j0FEc66+t+17nDuekQ/fMEBTnm1btsb92M+oCRpbQsjw2BWk4cRhaYJT5b5xJcxOC8vBxHzZmPgm3aiI0+uTK42tt+v6tjwrTtBhvs81YROeSU/OeED0RCzveG6iHq6JaimEjKjPq1ezUEc6Q2aT4e5yPzS2hJDBkJPLS2+WKre6xdZmgpRtTE159nXDOGy34HxxbRoq45Xx3oi0yS2z+dwwlLo6O3v7w6G65jBIMp+7Jxi1mY1sury0DHsr0xSRd9KkNlecppnZvnybc/BCY0sIGQy7a9HVBbzUheiaIGWmn7sr9eAEKSB+MlWIJmf3oZnOgfybs5XV1n2tz8wwJvYEKTPPqms5VC7RmLpYecVOkDLjJmsQMKD+TNuZzS7rc2lsCSHD4dvhpqMvYGcXoE+ey9A5ZkdnHwvsw8evLbPp3GJl5rSZLuX+HQY9FrvMorpzXXondCOPAY0tIWR47ArS2O90ZYKSb4IUsNOqmnsmPTm2y1tZsmHFNY3uyhpgS+/gNnvG0iSvAYmcRLTNHO0MRoTRbfy48C1jCrXwY+7Ndlbi2w6xiYT4nTxItVzDw6U/hJDikJOeKbcDuzHUwFrRUXajaVNj+/RskiWy82fLyn3vLS3R4M+gg8WksSWEFIfu3vJ3OZ/k+s9BGGqd7YZYIXYjE0LGwVWZRy39sdK6XD+a3c6Obfu8GxIMvfTHyrMXfEuEYvLt4L6xj6U/Tm9f2xfj3TWmeH9aiusrM8fSMRsaW0LIcNgbEXQ1NJ51tEHZ1tKf6I0IfBW7pxJfWpZkjnN2cEtYL4OKNix2niHjn9vTVMsu3u17jF2ek7Ie2NAtZQJc425R3GKPEFIsjgrKbgk6W0aOXYK2cVTSyY4PQnFdRiBUiUc4iOhlX92GMkmW0VIHM9eoSViu/EO6D9XVXWHrH2xpW2xIbzkhZDJ4jJxUfzGsxHOdO2RxQ/n+aF22sa15V7o21xJI+U2yZUsIGQy7xRrltSemVehoyXkNrrWEZ6Xd1KRTUyWe4kGqTatTdbl1H/QgZehi5uNawpO6ZtWVZxsHEzZtjWRsF3Obe7TT2sERstiyJYQMhr2UxjzvtIyjrfcgV7pcergcTORYqlKvX63/TGJr9Kb7TtHTzDPHPcamb7MOd6jfmAO2bAkhw2IbhNohhcc1Y4y7RgAr44MmS2dGC85l4Fe8KflwxGnM15E22V1jTGtNZNGqNQ2MI09vzikt7oA+bT+golO1NYAt09FdIyFkGpheoOpKuqk1NsPORgGulpzl8WnFKNsTsWxHFZ4t9qLG4lKMUq5WdCiN71qEYw6xPwZiu1ljW8Qd3HECgQ8jW5eI2d5ZnWHY3fgeaGwJIcNjV4jmuJ/NPC6eOZO5aV1kzBZ7UXSttNuO2XpayVnzSdEnRnbHZV6Ny29qXUJh9rKsOoon3BVnhcj7orElhAyHxw/uynGivE7rbD2OKJzpHDJMktfZ5uhCjp00ZBjpbOtsQ0ubelhnu33dpUMMnglcsRP12I1MCJkEUcZNM+1na8l0HjfFTZCzIs9j6Gx5yUtiUmcjuwxM7tnI27oYeblmescQMuAx6ZpI6SK3k7ZKtYCzkQkhgxFqGRQzG7kLQ8xGBsqbjVz/DTkbuTTZDbBlSwgZjJWNwo0WRieHEqldiWYLrsus2ybZrvSh8elYHewucrtbvkl+05rfri1boL1h66tFa8YdweiyZUsIgYicLyIfF5G7ROROEXldFX6miNwiIvdU/8+owkVE3iUix0TkDhG5ODnP0Lhh3XVsndfpxL7mSm+GhbomfV3DdTozfW1Y5o7rDl293o/mhiyXHN+fjzkg5m5IrrHwuvs5NDZu32vMX5Os1D9fWl8eQ9KkV0AfGltCCACcAPAGVb0IwCUAXiMiFwG4FsCtqnoAwK3VOQC8FMCB6u8QgOtiMrGdWNgOLpbWhJrXjWt1Otc1536rDrnOa7505p+50bzretNfndbsAm4jA1iRpVuyLHtm5FfnE+ruDekZ+jPjtimTNvft07mJrt3cTXoFZEd3I4vIZQDeCWALwHtU9Q+s688CcCOAfVWca1X15lj5hJDxUNUHADxQHX9bRO4GsB/AFQBeVEW7EcAnALyxCn+fqiqAT4rIPhE5t5LjywSz7z+1OHZ0V4oRb3u2cD2ZxQyr04daFXalF2oBpcSdG9dmAQMQkjO3wrccbR5X97Ypz5bh4qkTwEwg3z+xo2udzqd7SM8YbLltDFvb1qpP35h7jaWDoY4ytiKyBeDdAH4WwHEAt4vIEVW9y4j2/wD4kKpeV30R3wzggtaaEUJGQUQuAPBcALcBOMcwoA8COKc63g/gPiPZ8SrMa2z1ie/j5Of/Kbu+xMO88gTy2HfG1YMAiG/ZPg/AMVW9FwBE5ANYfNmaxlYBPL06Ph3AV3MpSQgZBhF5GoC/BPAbqvqYmN28qioiSc0OETmERTcz9uLUHQNAhoNlXgSxY7a+r1iT/wjgV0TkOBat2l/rrB0hZDBEZDcWhvbPVPXDVfDXROTc6vq5AB6qwu8HcL6R/LwqbAlVPayqB1X14G6c0p/yhBROzqU/VwF4r6r+vyLyAgB/KiI/oapL89XML91nPetZGbMnhLRFFk3Y6wHcrapvMy4dAXA1gD+o/n/ECH9t1cv1fACPBsdrAcie3dh13o/UGYZnmLrGK2PGy3xLO3x5hWQ2ObdoWqvqk2GPLaeOAzbNTK6Yf/2bwGyG2Zn7VseRZzO/fmY+viVSseXZ5/Kf0G8olZglT7H3cq87ONbYxnzFXgPgMgBQ1b8Xkb0AzsLOlzCqa4cBHAaAgwcPDjxvmxDi4YUA/jWAz4rIZ6qw/4CFkf2QiFwD4MsAXllduxnA5QCOAXgcwKsbc9jawskzn7Y49lWUc6xuJuAKs9LULK3jNdae+tbwLsXzyFpOEF5HurLrT5MMADqL7WA0ZM5NRd1GQL7zXUAE832rZe7Ks7UbRNMQW7q03vUn9Lx88WINdKRuPv/JXho2Wog1trcDOCAiF2JhZK8E8EtWnK8AuBTAe0XknwHYC+DhSPmEkBFR1f8BeL3RXeqIrwBek5zRttMFdVZM5mYCAHZcM6rubLM3d6ep40Gq3X/M3Xtqg2C1fM14K/J8hjLQolopwAiD3cpdY1Me5jVTZ6t8W8ttSlthloem3Ksdz3heMfk2ym1R7q74SwY48FEHRI7ZquoJAK8F8DEAd2Mx6/hOEXmLiLysivYGAL8qIv8LwPsB/JvqhSSEELLhDOYhrEuaRFLuKXrMtloze7MV9rvG8V1YdEURQoiT0EYEao6bGS0G26Wjdw9Ws0Xqarm4WqyheD7ajOX60sWORdtpQu4a7Y0IXPm4Wuj2GG0svo0I2tLWW1RMax3oNtbruK9Yg0sPUoSQQTG7ir2+kQOTgLJtDxeTzjcmOEd8Wjue7eowtks3NDHKlOdzo1gb6aZ82xpaW04srnJKTZ+hy7tvaGwJIaMQnHASmKXbaXegVHx5da05zfQxs5pT79mnX+4aP4c8z+SqlWgnOu4+PzI0toSQQVnxhzwVbF1j6/5cBju0TMi13Z5Pj5l1HptHCl1luGZ67/LcYC6de4Zb7BFCBqUedw2OddVdzZ70rrjbx6G1rl27G+sx0fq/a9ZsTH7mVng5ujVtw2/v+lPnU4dvefK1u8rb0nWs1BpXX7nWZmw5B55Z0t7rBjS2hJDhsI1ozKQhe0zO5VyhycDZRmTWbOhXmFv/fWltQ2AaaJfMrQY5tl6OJSb2/sArPQZ23r4PCQR09WDqsrSMym6Itl3208dHQej5h0jZtMKCxpYQMixmhWWvfe0ySzQ2bdN4aShdSsVe62RuQ9el63k7TL3d78FueTt/u8xmgbgBlvIcY2DSnIkdE28kaGwJIcNgLc0JtmAiupij4rqum8tjmpbeeFp7LmcYSzLt9AZLLcGWHxjR6zvny7pst3zhyLdFC3GlhW23bO0eiWi9DRkx8VJkxsj1YX8YshuZEFIcxkQWFfG6qzLR2WynQk9phQZbeLraug7FreMYXY86l9UK255Vq+ocG10xSK6ZxyGjNHMYegcCrGwgv53Odc921+rSml0s329luFZa0r5eA1fZwIq71MqOXH+cqzu4bfqEWeOcjUwImRwaayxjrjfRcolSsg62rFyzbNvW8qF0sTJ99zCRGcROIpcq2dDYEkJGoes629GXDaUY/BRZrmtj36vLUuRaY9ulHIculw750dgSQkahcelPgwep4IYBOQm5D8ylQ8uyiKLtEh57+ZB9LUZuzHh626VZbeJ3pcMz55gtIWQU7OUqSwRaEGbLVnyTrXIRciIR2yJzjo9ix1g1eZBqs9zEnlnsOnaVmzmm6ktnn4eMbkoXv5l3bAuyz2fvy891HAGNLSFkMOx9apfC4Da8pgMMu+t4yeACSxWv15hbFeZK/rGTeFz4ljWZ1+pJPa7JUTFYk7W8zDyGzDz35T/EUiofobzbOrho+qDpSPDDsYLGlhAyLDPrv7HMQ+eys29tHWeOZaNotaqWvFE5xnNXZj4HWk4uWd50EWtdV/JtiB9D1LIf05g6DGvjh0hK6zJgGNuMq0vblm2u2eoBmu4ndJ1jtoQQQkjP0NgSQobFN/GmxE1dYiZHeehjAle0Mwszb/tYdWV7w6V0MS4wDVm9MPTEpwFgNzIhZDBCm8dvo7rd/arzqlvXCltJa8lY6cwLGI5gXPO8yUevLc/Osz63PFhFG1BTluvDxJx0NUOezeNjZgq7No+vaNVZm/qRMqBh7jLKy5YtIWRQfJvHLxGx9GcQQgYmdsmK61ruzeNrmaHj2kjHbh4fs/QnR2+ErxWeknYoOuRHY0sIGYXJbx7fxRlDm83jU/KIWfrTVkZ9HrOHbgyuyVypaYeiQ37sRiaEjML097P1rA925WfLMluEfbTOXPvZ2uEpMtrQZumUK33IwI3dsk1YA82WLSGkHDpMSArKKbVbcpNoW7Zr8kzYsiWEDINhMKU+98XrEt40wSm0ebyZ1teyq7aqc+6849o5x2Jls/XUnWscm8f74+qSrtvOQVzJLV19eSytRbavrWw72NJQ2uXmk9O29d1WrybYsiWETJ4BnBYspUsd6yy0Nm07xq2VUwrTPabrOJo2euQaFy4AtmwJIeXQNI6b4tWoTeXua9nGtKDsPV8LYXu5VYt0ofNk2qSP3UR+AtDYEkKGIbSdmr3G05yhak+U8aWrz+049sSc0EzgpZatZ4IUqtaiywDY8hwyVrpak2cyp20eb7q5XNk83lcuLj2jdLOVaNnD0NTNa/uZBpo/iEIG2/UbCsWJkWNBY0sIGQzXRgQAlgzp0pjmlgCVY4ul1tn27Fpr/NDhG3nF1651brf6lnwj28czwxmEx1A3OqWfWfp6dA7KSdmIwFymNGv4cEnNw8b6yGlqUZsbTKTu/rTCVnzU5F2aLGp9U3oM1qBxTgiZJG1rH7sF1sE5fDIxa2JjrxmEdGw1PjoR7HtTy1iXiDo+jmKgsSWEbCMiWyLyaRH5aHV+oYjcJiLHROSDIrKnCj+lOj9WXb8gpx5FGZcYY9qkb4MxjQkbndh7bUmR92zRRUcaW0KIyesA3G2cvxXA21X12QAeAXBNFX4NgEeq8LdX8bJSfOU7sH6tyqOrjqFx9hIoTZ8ANLaEEACAiJwH4OcBvKc6FwA/A+CmKsqNAF5eHV9RnaO6fmkVv5m27gPHIEdlntEgdDa4gbHqlTghI9vGtWIKEzKisUT/1EXkMhH5QtVtdK0nzitF5C4RuVNE/jyfmoSQAXgHgN/GztSlZwD4lqqeqM6PA9hfHe8HcB8AVNcfreI3E7GkJrhJwZCkbjDQJV4DrcvE464yuAlE/d/1Z17rizX0+BU1G1lEtgC8G8DPYvHC3S4iR1T1LiPOAQBvAvBCVX1ERH6oD4UJIfkRkV8A8JCqfkpEXpRR7iEAhwBg7+6nL19s6/0n5OnIMUu43qKvtW/ekPxYX7mBfJdmXxszkEOemqLpwxglGEKBuwVt35fzXmN8Iw+MS+9YYpf+PA/AMVW9FwBE5ANYdCPdZcT5VQDvVtVHAEBVH4rWghAyNi8E8DIRuRzAXgBPB/BOAPtEZFfVej0PwP1V/PsBnA/guIjsAnA6gG/YQlX1MIDDAHD6rrN19qWvNmuic0AWnW5ihIkEOuJ02QKvxFW/ZU+J24jMljdPiJFl5G+m9VbjtUxzv1oH8+9+DwAwO3kymOeK3JzY5eGL1kYXmfWjc1Oe5mlC0lhju91lVHEcwPOtOD8GACLyd1isePqPqvrXCboQQkZCVd+ERc8Uqpbtb6nqL4vIXwB4BYAPALgawEeqJEeq87+vrv+NarjJoydP4uQjj/SiP/Fz8ltPjq0CQd7pCbsAHADwIgBXAfgvIrLPjiQih0TkqIgcffjhhzNmTwjpgTcCeL2IHMNiTPb6Kvx6AM+owl8PwDmPgxCyILZlW3cZ1ZjdSTXHAdymqk8B+KKI/CMWxvd2M5LZrXTw4MHpjG4TsiGo6icAfKI6vheLYSQ7zhMAfjFFruzejV0/XM2vCo2futw1Nvk6brvPaIpMAJjPgdls53+TPJ+MmtksfUxSdVmGh/mjjy2yOP3pO7rW6bYc7pZW9txt0UVrl0nb8dbY5zX0BKmY+7nPHRxrbG8HcEBELsTCyF4J4JesOH+FRYv2T0TkLCy6le+NlE8IWXd2beHk2fsWx4HNzM3JQQt3jcvuC5f64+Y7aUy8LhiBeHeNdlxXeMAhxbYce3KW5V5SXUbbwZK8CCMjTy66j+dn71u9B1eeTVsVmm4qY4xhS89XEvNx5dMxOpN2HwFR99PF2KrqCRF5LYCPYTEee4Oq3ikibwFwVFWPVNdeIiJ3ATgJ4N+r6sqECULI5rLtC9f2aWxiLnGp41UzWwEAJz0GzczHIdN5nBLX3hFoZs1yrnWy49v/zQ+MmULqFqTPf7PLCEbtQlRNoDL3763TiaPsV1q2WP4omhn34ct/1lC2vjxtA17o0p8VP9sJRG9EoKo3A7jZCvtd41ixGLt5fStNCCHrjdEik/rcF69i2xDZBjKlYnYakcj1s3MA0JWNE5Y3j9ed/zNXejOOtdTH3Dy+qWVpyBxs83hjIrMaG0LEbR6v8ZsZuJ5RXWahhv8Ym8dbPRPea3aW7XMkhJDMhIxIjIGpK19fa9JVOYdar9ZxEY42RqK+98HLYI72RrVvEsqCxpYQMjypFbZtRLvSZHRJeZRscCN+OzS2hBBCSM/Q2BJChid1kkldU+Vy3eeq+QpyC0gclGqtIjdlKFV9Qsgm0nUjdtso28uF2hrZKl3x2/71SNtN07NQqqVKKItSb4EQss5EzEReYo64sbHYCVLm+F/CxKsuk4OkjZOIghhtghRQ7iSpyPFagMaWEEIGIdZ5BVlPotfZEkJINnxeiGwvRfX/GfxddqYsbzeycd101mDGccmyHDusdKWGuqVrvWsZRsssuH1e6laADlRk26FE227f2HTe+wiNs/u8U7nSh7DKtXeayiRwnZ9ahJBhsCuiDhVX7KSUzoRqyNy1p31POe7R2oKvqDHnknQZABpbQggZAbUNq49Io1SUIa0J9Ua4rpd4D5mgsSWETIs1rpA3kg15nhyzJYQMR1NLJnWs0toEYMVvrUueb7zWDnONITaND7rGf61lR0v+g31j0aFymqnlg3gHb+vWHne2Nz0w42AnbgxLeXZpvpmt3Sbfxfa4OBDWt2lsvWmrxQwfBDS2hJDB8U4Qsio1c6KPa1s8MSdRmXF83ZTAauXpM/i2XF/Fbht8uxKP/YCI3fUnN035uCaVxRjimDFn364/pl6xpBj62C581/W+d/0hhJBOtNn1p07j2hav7c4/qbv+1Dv5uHb9Ma/ZMs01vkar2Lvrj08Hx444yWtdrbxjdv2x87Nb0y4dllvt6m8Vuj5EUp5RjN4huuz60xIaW0LIMFQtHbVagNsGwGztGkts7JatmT64fMbHrKFla7JlHde6wmq1mfHM+1PduVZV8NtGujbCrvW3Qf3m3m7kbWaA1DORjSU428VVtzrN8jPLxSpXdexMux22tD+vdQ8pE6TMvLdWk3hpYzjbdgt36E6msSWEDMbK+lRj7Wp9fak6m1X7qLrSYqclZVeB9tjlShXp6JJ2xrW7OUNjw9bHwFIXtzXmun0f5sdHA0sfIqkGxtcVGjKGsa3nwJrlNjOkt3stYtK2HU9taTS7zPimsSWEFEcxy1hil+QMoG/QEYaPrnr5JhGVsh1hykfByNDYEkLGITDWJmZ3rSudZwzU1dJZuZ6CPfEqNHHJNzbZkO/SGK7RKl5qIbelqyHqWHZLPRIG9f05hxDs/JpmCg+IrWPKRyHX2RJCyiUwaSeJNpVyqqExHPW3NZBmui5ySsK8D/vYjONlDcoAoLElhBSEt9JNNDy9GqmeK/91MLDZSfjY2Z697pjF3jehZ8duZELIYGxXRic9EeoKsu5CNcPq9L51sVb6peVBdsVrdFGurNV1ya1ZmnmbYAB8MpA487aWGdOyn+vCN7K5jKZOJ4n3GkugTJI+Ilxx7UlpdtgcAFx56Gqz0jfxzZen53rKPdHYEkIGp3HXG995zBiZHd+uWH0yY51I5NxpppS+xdInGsU89z6cWqQ6vAiUYSmPmhBCulV8obilzG7eZDb8GdDYEkLKITBmm9TyiukSTsmfdGfMWdV9wW5kQkhxmGOv9bkvXkhGU5g9Bmv/T3XXCCx3HXvcF267bjRxdDdHuWsMkeKucb48VtvKXSOWlyKZYTarZeJw/hGDWe7bwjONK9t6tcEc642ctc6WLSFkvenafemoJaN3uvFcC67PLLC71aWvNnmIajuGmpK+NAsWKI/SVCWEjISI7BORm0Tk8yJyt4i8QETOFJFbROSe6v8ZVVwRkXeJyDERuUNELs6szPL/mLhjYG2fB8DtE3jgmrYYD1zAcnmE/CVPAfteEqCxJYTUvBPAX6vqcwD8JIC7AVwL4FZVPQDg1uocAF4K4ED1dwjAdVk1ifC8tBK37fUmQl3Bc+PPzK+hO7nvtbS55fvkReVjlodv7L1tF3GuWeGx2PeSAI0tIQQicjqAnwZwPQCo6pOq+i0AVwC4sYp2I4CXV8dXAHifLvgkgH0icm5DJuFz+5qvFWRfS12ekcIM/m7kGZb/bP1sOVZ6tVt7djpfuAO7S3dn0wYJxnPS1G0eEb6N61k1/dU6jN2N7DLkbNkSQjpyIYCHAfyJiHxaRN4jIqcBOEdVH6jiPAjgnOp4P4D7jPTHqzAyAvQ61QOZrSONLSEEWKxMuBjAdar6XADfxU6XMQBAVRVuFz1eROSQiBwVkaNPnvhuNmUnhaNV20hiq8nZqm3DFC1CbCt4ZKJVFJHLROQL1YSIawPx/qWIqIgczKMiIWQAjgM4rqq3Vec3YWF8v1Z3D1f/H6qu3w/gfCP9eVXYEqp6WFUPqurBPbtO6035jWDKE4tInLEVkS0A78ZiUsRFAK4SkYsc8X4QwOsA3GZfI4SUi6o+COA+EfnxKuhSAHcBOALg6irsagAfqY6PAHhVNSv5EgCPGt3NpA8iu4o7dSkPPeEoB/YEtUKJdWrxPADHVPVeABCRD2AxQeIuK97vA3grgH+fTUNCyFD8GoA/E5E9AO4F8GosPsg/JCLXAPgygFdWcW8GcDmAYwAer+KGeeopyFe/3qyFGjWnzHbCxNM2UNd04Vn4uit+TLxCsdu8AmD++OMAgNmDD4+mx9pj/yYDv6FYY+uaDPH8pTwX6+zOV9X/KiI0toRMDFX9DADX8M+ljrgK4DVJ8k+cxMmHh6v4yYKT3//+2CoQZBpWFpEZgLcBeENE3O0JEw/zxSOEELIBxLZsmyZD/CCAnwDwCVkM4v8wgCMi8jJVPWoKUtXDAA4DwMGDBzlfnZANQUQw27t3bDU2hvmTTwEAZnt2j6zJhvE9d3Cssb0dwAERuRALI3slgF+qL6rqowDOqs9F5BMAfss2tISQDeaUPcCPX7g4Djml920S3rTBt4kZ1+XNyXZC4ZNnxjWHkk3H/DMrnn1vAe9R284xbN1Dk5wiNyLY+upi4rjuP2c7bHsjgi1HWbYYsrb18PqMbtr31VVmMf2urTciaJmuJrTp/GfcSaKMraqeEJHXAvgYgC0AN6jqnSLyFgBHVfVIS5UJIZuC4X1HRdyTaewKzN7QvWn5i5neNrh1BWkaRV+8+tislA2jutitx7rmkmnLgLUzzsxzT0EH/44dhxwIsONBanueWcD7kb0TToQhW1nT6/pwcCrXcM++XXlW8mvZOZpjGZXv9+Mheos9Vb0ZixmIZtjveuK+KFYuIYREE2rhliSTlE3OZx4pi/vZEkLKJ9UfbZMfZlfruT6m68PyyP0xlEteghwaW0LIYJjjnM6uUJHlMc0tWYxRwjEeWG+I7pBhxhXrmt0SWdkU3Y5rHs+qtDNHa8boIgcC46r2hupWujptUE5M96nZhbxd7pHd8bF52DjKxOU+sr6/7TFk67l75eWio9xa35V7C8idgEdJQsha0rb2sfaQXarwHJVd1r1dU8dXI/Ju2oknaqeeEgnorcYHhtevc6H3rI6PoxhobAkhxVGUgYkxpk26dthYIBTWO/aWcuZ5Zlrd38Bl0uUZ0NgSQobFNYvXPh6IqP1YfccNaTVkiFtW2tGVvW0gzePqQ8Ypy47v07dnw1tqq7YLNLaEkGHxLSkZwTVxktN+M25Eum3Zrrh9T8Ly5Vn9ee87Va++7mMNJ6nR2BJChsVX60ypNooci42Nm52U9a0lMhU9E5jSz5sQQgiZJDS2hBCyCfTRWlzDFmhf0NgSQsgmEZrpTePZGzS2hBCySYQmSJHeoLElhAyHq8bpUgu5HFskyF3auSckC0ibFNVEH8tmEvJ2OmbIoU+GSWGxZbg0AS0lv5Fa7zS2hJDhcC3v6bLkp3bZ6GqpRcjdTjdvkAWkLfeJIXfrsoU8SVzOFK1DB1mxZbi0tGrMco+ExpYQMhyF1ThJHoE4nhnPAC3bHHkNCTciIIQMR93adGzIbm/ODgDq3vU2G05n8ib2RgSu8Prc3tDAbHnZeXRthTZht/pq/avwqNa7b4N3uyx8cloaQTF1bmLgVmrj7yUAjS0hZDC2jejcs8sLsFSBysmdsO0q7qQ//na60HXrfKXq9Bkc83iOxa44vrzt7lQ7bc1MIXOrv7vJmKnGdb3PdbHzT62rmbdEeLVayUO3dfbmP1suL/GN4foMuXk9xuCO0CW89LtNMLw0toSQUfBuqwYsbTtnx7G3oms0TjEtsVBaF7PANTNPM46pg7FFYLLv5Jh76eJ/utY9FNe85jO8oYlLoXtOaRV3ebZt4EYEhBBSKDmNQQ5ZI/igJjS2hBBCSO/Q2BJCSJ+EuqJzyYqRZ64p9smayMzeKUJjSwgZFHPM1XdcJCUvSSmtJh/y3vvcVzcjpT0iQsga05tBbSt3DC9OuWtdW159PmuYwBUTlkoOw9e34aQHKULIJtLFAIfSprhdXHFb2NcOObnl+gxtLLk2JDDcXQ5C27IcsQXMpT+EkOExl74gU4u34zKQ4FKkvmhrMGaJejYtUwLal93QTbZcvRipa3g7/kbZsiWEAABE5DdF5E4R+ZyIvF9E9orIhSJym4gcE5EPisieKu4p1fmx6voFfenVZVw3Jn7xY8VkhyGHC+rWc0QrWkPbFlbQ2BJCICL7Afw6gIOq+hMAtgBcCeCtAN6uqs8G8AiAa6ok1wB4pAp/exWvOR+1PBlVYaEWpXnNjhdypO+V64nX2rF9QPbKtQZvVtF51BsnGBsoAPZxwPWifa+2bil6mXp0KbuYvH36xsh03WPqn4em3zBAY0sI2WEXgB8QkV0ATgXwAICfAXBTdf1GAC+vjq+ozlFdv1QkrvlQV0p2BbViOD0VWLJxbHvNlm27XAwZqlBFbRvJHBW+Kc9ngE0j7bt3Mw/boLvy8DnI8Okde3+xBrWt4W1L07MIyKaxJYRAVe8H8IcAvoKFkX0UwKcAfEtVT1TRjgPYXx3vB3BflfZEFf8ZQ+q8tpTWrT0z/k9lfLZAaGwJIRCRM7BorV4I4JkATgNwWQa5h0TkqIgcffLEd5sT+FoGsa0XVys0FB7K33UccnXYtfs0t5xcxLZmm2jbZZ6ja7ovEvSisSWEAMCLAXxRVR9W1acAfBjACwHsq7qVAeA8APdXx/cDOB8AquunA/iGLVRVD6vqQVU9uGfXad21TK10Q8a7rUwyDqU+p0i9oo2tiFwmIl+oZh9e67j+ehG5S0TuEJFbReRHEtQlhIzLVwBcIiKnVmOvlwK4C8DHAbyiinM1gI9Ux0eqc1TX/0a11NqQkPGJMrYisgXg3QBeCuAiAFeJyEVWtE9jMZPxf8diwsR/yqkoIaQ/VPU2LN7bfwDwWSzqhsMA3gjg9SJyDIsx2eurJNcDeEYV/noAKx/ghIyFfO/J4nY3inVq8TwAx1T1XgAQkQ9gMb5zVx1BVT9uxP8kgF/JpSQhpH9U9c0A3mwF34vF+2/HfQLAL7bNa3DnEWSz2LN72PxUGydzxRrb7ZmHFccBPD8Q/xoA/y1SNiFkEzg5x+yxxwEAOhPIPDCealdcTZWZbbxd6V24ZDYtFXJtDO+T2bQ8pI3bwcgJQ/MTi0nks+88vuolqs63y7IoF6lemnLl61IlixRbaHup2d01isivADgI4P/wXD8E4BAAPOtZz8qdPSGkUPTJJ3Hi3i+NrcbGMf/2t8dWgSB+gtT2zMMKc1biNiLyYgC/A+Blqvp9lyBzduLZZ5+dqi8hhBAyOWJbtrcDOCAiF2JhZK8E8EtmBBF5LoA/BnCZqj6UVUtCyHow2xpbg81hfnLxn2U+LCfdwVHGVlVPiMhrAXwMC5+pN6jqnSLyFgBHVfUIgP8M4GkA/qLy2vYVVX1ZBtUJIWuA7D0FswMHwnGs8brtnXiMMVLb4fuKm0cj3pJLR3OM0pJRywz5Wl6RE7Mn7IhjtnjgIUBmwA+ftTpmO3N0avocgaSQY8zWl3fq+HqIHjxTbf/ePue+Hj1mq6o3A7jZCvtd4/jFLXUkhGwCIsCuRSWvIpD5fCe8YmWprggUAObznXi2ofQYW4jsXGswttsGusnYzrGzPaBrEC7G2JpLUmaONE3UPo4bEJkBM4Fube3oOgdkPofuijC2bZbO2GI3xdi6fj8W9CBFCBkHV6suNMM3phVoG2Tff1+eTcehGjOmpRtDzxvMq6tV60tn+kOeOcJCuJ6t7xn6wmKf+9CYv6dI3bh5PCFkOFS3K3uzO1hUl7t9XV3BtQgrnW/z8+209hIXn0G1w+yuYlOOr0Uas2TJ3vi9jZGO2Tx+ZsiYA9gyNp13LU8KtcpdxjWmZR+zPKopro/Y5V1d8gilTbwHGltCyHB4KqO2G7ivGOiW3YrOMdshGKrFFtMin4qjkbZl5jOOvvH3zOVBY0sIGYy2RpV0ZF3KPbehDcnMXGYcsyWEDA9rHrJhsGVLCCkLx7Kcuns3qmVszm72dTNbcZL1m0qX6zoROzkuh6wEVuYQeOD3JSGEVLCbuz86lW2hzyXlnmhsCSFrB40mKQ0aW0IIIaRnaGwJIYSQnqGxJYQQQnqGxpYQsnYM7pyCkAZobAkhhPTOOn4ApdwTjS0hhFSso0FYC9bgudCpBSGkGFY2Hqi3sxNZXs5jbBnn28BgSWZ90uCQwukEw7cxgZGfLQPA8iYIgXxdy5TqjRnssFbkdHFoE9LJfmZN2Xk2owgnGtbBSNP9hK6zZUsImR5zbO+32mSEOrdWQ+nN/XJbyhDVoI5N1xtpmza0t2/kBvY+vetw+96S73PgFq+pb6quNLaEkOmRUHMV4eCigw5F6N8Sn+7q6IFYd2hsCSHlMk8Md7Q2so7DumRZYUuttdiWr53OOBe7hTlka87Mrz6Ozd+lf4V9b85nNPS9BpDvfG/5vKHF7oJjtoSQcajHYx1st3hcccwmgmfM1paVbbzT3JatodUWHLM1zl2tO1+Lbzu06X7sTe8d+QbT+vZ4deGTFxizDbVoJSXvgcZs9Wk/sPjPMVtCyCRZ4xpITcPsokUXarCyj90VJ9SFm6tbd4O6h2NZ4586IaRIfLVOh9qo7dhfUrqIuKlbAGbBbmW31TOHXoEWfyyttjwcaXu9FGhsCSHD4apx+qqFfHJdla0RVxOXrESTwRA5ZaacW9T3unK/fRu8jrR+RiO2uDlmSwgZn7YG1xjT7TQ2a8kajKEMRp9GZmAD1vpDyEpnjq371jpnybeCLVtCNggRuUFEHhKRzxlhZ4rILSJyT/X/jCpcRORdInJMRO4QkYuNNFdX8e8RkatbKBIft0stFZOWteAk6NTjEEgXmsTlbfnbzND4O+LPjJDN4r0ALrPCrgVwq6oeAHBrdQ4ALwVwoPo7BOA6YGGcAbwZwPMBPA/Am2sD3TsRNVbMDFhvBcoasT+GKFtXd3bCbO9WRN4Xf1qEbBCq+rcAvmkFXwHgxur4RgAvN8Lfpws+CWCfiJwL4OcA3KKq31TVRwDcglUDvr6UMtO2kDWoQ5A88axpJvgI0NgSQs5R1Qeq4wcBnFMd7wdwnxHveBXmC8/CJnkVchJTK2+QoW1ND7+jlTkBCRaUxpYQso2qKoBsNbmIHBKRoyJy9MkTjzcn8NRIKpJWW1lx2xhwu9s5On7E+J2XmefYZuDZv50wy8NXNoF7LenjK2rsNv0SIWRD+FrVPYzq/0NV+P0AzjfinVeF+cJXUNXDqnpQVQ/u2XVqsyae2cCi2mmmcJtZyqkO8rfjGJskePHJm3uOp4xZHr6yCdzrpLY9DNwHjS0h5AiAekbx1QA+YoS/qpqVfAmAR6vu5o8BeImInFFNjHpJFUZiKai1BgAyXxfLXi5cZ0vIBiEi7wfwIgBnichxLGYV/wGAD4nINQC+DOCVVfSbAVwO4BiAxwG8GgBU9Zsi8vsAbq/ivUVV7UlX4zLHdJsSI9g9nc2m1YIshYTfWbSxFZHLALwTwBaA96jqH1jXTwHwPgA/BeAbAP6Vqn4pVj4hpH9U9SrPpUsdcRXAazxybgBwQwdFls+NSmvFyYDdtdpQuS2lny+H28cqsiwzxtDRKLUjwzBAcMw0YXMG8zfS2alF5H1FGVsR2QLwbgA/i8XMw9tF5Iiq3mVEuwbAI6r6bBG5EsBbAfyrODUIIRuBuSTDrMxm1U4vM4XOjclQvjCT2c5Y6dIkJZFF00Cr9C7qeLYswJnGDIlZp9skw5WmCZ1bnrIiNiDYNiyBPG1dozu6jV16VsqkZe9C8Hl54sXo69Jv6bdVM893L6nJnwfgmKreq6pPAvgAFmvwTMy1ejcBuFSksIEJQkiZ5KgqJtht3NkjEqvY9sys/65rGYntRnatq3u+L46qnhCRRwE8A8DXuypJCFkP5MSizy3ox9iYebw9lqgKaN16MuT5NmdXwdIKJnsDdrOb2TZYIb1q5oBsOQydKc/UzWzJm13bMwCaWLPP54YMj661SFVgPl9u9c0BUVntXTD3BI7tKrf0EPNWRJYfViz183Z9SLi6bBO69bcl1sMHKbh+KwkMPkFKRA5h4foNAL5v+mgtkLNQ/scCdexO6foBwI+PrUBXHvveA9/5/z79li+MrYeDUp9/Hr3y39l6l1d3fsQVGGtsY9bV1XGOi8guAKdjMVFqCVU9DOAwAIjIUVU9GKnD4JSuH0Adc1C6fsBCx7F1yMAXSiznUp8/9UqjVL1qYvsvbgdwQEQuFJE9AK7EYg2eiblW7xUA/qaazUgIIYRsNFEt22oM9rVYLFzfAnCDqt4pIm8BcFRVjwC4HsCfisgxLBydX9mX0oQQQsiUiB6zVdWbsVjkbob9rnH8BIBfTMz/cGL8oSldP4A65qB0/YBp6NhEqfdAvdKgXi0Q9vQSQggh/TLBlWmEEELItBjE2IrIZSLyBRE5JiLXOq6fIiIfrK7fJiIXDKFXgn6vF5G7ROQOEblVRJxTu8fU0Yj3L0VERWTQWXkx+onIK6tyvFNE/nxI/WJ0FJFnicjHReTT1bO+fGD9bhCRh3zL4aoNAd5V6X+HiFw8pH5tif3t9pT3+dUzrX93r6vCzxSRW0Tknur/GVX4oGUsIlvV7+2j1fmFVR14rKoT91Thg9WRIrJPRG4Skc+LyN0i8oISyktEfrN6hp8TkfeLyN4SyisaVe31D4sJVf8E4EcB7AHwvwBcZMX5vwH8UXV8JYAP9q1Xon7/J4BTq+N/N6R+sTpW8X4QwN8C+CSAgyXpB+AAgE8DOKM6/6HSyhCLMZ9/Vx1fBOBLA+v40wAuBvA5z/XLAfw3LNbmXwLgtiH166vce87/XAAXV8c/COAfq2f7nwBcW4VfC+CtY5QxgNcD+HMAH63OPwTgyur4j4zf42B1JBaeAP+v6ngPgH1jlxcWTpO+COAHjHL6NyWUV+zfEC3b0l09Nuqnqh9X1Xrn609isc54SGLKEAB+Hwuf1E8MqRzi9PtVAO9W1UcAQFUfwrDE6KgAnl4dnw7gqwPqB1X9Wyxm8vu4AsD7dMEnAeyTah/agon97faCqj6gqv9QHX8bwN1YVNxmnXMjgJdXx4OVsYicB+DnAbynOhcAP4NFHejSq/c6UkROx+Kj73oAUNUnVfVbKKC8sJjQ+wOy8ONwKoAHMHJ5pTCEsXW5etzvi6OqJwDUrh6HIEY/k2uw+JIbkkYdq+6b81X1vw6pWEVMGf4YgB8Tkb8TkU/KYhepIYnR8T8C+BVZbD13M4BfG0a1aFJ/qyVQjM5VV+JzAdwG4Bxd7M0LAA8COKc6HlLfdwD4bew4DnwGgG9VdaCd91B15IUAHgbwJ1X39ntE5DSMXF6qej+APwTwFSyM7KMAPoXxyysaTpBKQER+BcBBAP95bF1MRGQG4G0A3jC2LgF2YdGV/CIAVwH4LyKyb0yFHFwF4L2qeh4W3WN/WpUtmTgi8jQAfwngN1T1MfOaLvoaB12WISK/AOAhVf3UkPlGsAuLoYzrVPW5AL6LRbfxNiOV1xlYtFYvBPBMAKcBGPqDvRNDVCQprh4hAVePPRGjH0TkxQB+B8DLVPX7A+lW06TjDwL4CQCfEJEvYTF2ckSGmyQVU4bHARxR1adU9YtYjJ0dGEg/IE7Ha7AYA4Kq/j2AvVj4Wy2FqN9qYYyus4jsxsLQ/pmqfrgK/lrd3Vn9r4c1htL3hQBeVr2vH8CiO/SdWHTD1v4PzLyHqiOPAziuqrdV5zdhYXzHLq8XA/iiqj6sqk8B+DAWZTh2eUUzhLEt3dVjo34i8lwAf4yFoR16rLFRR1V9VFXPUtULVPUCLMaVX6aqQ/nTjXnGf4VFqxYichYW3cr3DqRfrI5fQbWJuoj8MyyM7cMD6tjEEQCvqmaAXgLgUaNrr1Riyr03qnG66wHcrapvMy6Zdc7VAD5ihPdexqr6JlU9r3pfr8SizvtlAB/Hog506dV7HamqDwK4T0TqjS8uBXAXRi4vLN7NS0Tk1OqZ1nqNWl5JDDELC4suuX/EYlbi71Rhb8HCIACLSu0vABwD8D8B/OgQeiXo998BfA3AZ6q/I0PqF6OjFfcTGHA2cmQZChZd3XcB+CyqGYSF6XgRgL/DYsbsZwC8ZGD93o/FeNRTWLQwrgHwbwH8W6MM313p/9mhn3HOch8w73+BRZfnHcb7ezkW43e3Ariner/PHKuMsfgIrWcj/2hVBx6r6sRTqvDB6kgA/xzA0arM/grAGSWUF4DfA/B5AJ8D8KcATimhvGL/6EGKEEII6RlO/iCEEEJ6hsaWEEII6RkaW0IIIaRnaGwJIYSQnqGxJYQQQnqGxpYQQgjpGRpbQgghpGdobAkhhJCe+f8BSiKccNwJ0s4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dummy_image = numpy.ones(mask.shape, dtype=\"float32\")\n", "dummy_image[::5,::5] = 10\n", "#dummy_image[mask] = -1\n", "csr = csr_matrix(pre_csr)\n", "dummy_blurred = csr.T.dot(dummy_image.ravel()).reshape(mask.shape)\n", "fix, ax = subplots(2,2, figsize=(8,8))\n", "ax[0,0].imshow(dummy_image)\n", "ax[0,1].imshow(dummy_blurred)\n", "ax[1,1].imshow(csr.dot(dummy_blurred.ravel()).reshape(mask.shape))\n", "pass" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "ax[0,0].set_xlim(964,981)\n", "ax[0,0].set_ylim(0,16)\n", "ax[0,1].set_xlim(964,981)\n", "ax[0,1].set_ylim(0,16)\n", "ax[1,1].set_xlim(964,981)\n", "ax[1,1].set_ylim(0,16)\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Least squares refinement of the pseudo-inverse" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 5.58 s, sys: 4.58 s, total: 10.2 s\n", "Wall time: 1 s\n", "(1, 30, 0.0046368041028449205, 0.0005109819474116924, 2.1354864217638663, 4.833787403420624, 2175.5691041334367)\n" ] } ], "source": [ "blured = dummy_blurred.ravel()\n", "\n", "# Invert this matrix: see https://arxiv.org/abs/1006.0758\n", "\n", "%time res = linalg.lsmr(csr.T, blured)\n", "\n", "restored = res[0].reshape(mask.shape)\n", "ax[1,0].imshow(restored)\n", "ax[1,0].set_xlim(964,981)\n", "ax[1,0].set_ylim(0,16)\n", "\n", "print(res[1:])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Pseudo inverse with positivitiy constrain and poissonian noise (MLEM)\n", "\n", "The MLEM algorithm was initially developed within the framework of reconstruction of\n", "images in emission tomography [Shepp and Vardi, 1982], [Vardi et al., 1985], [Lange and\n", "Carson, 1984]. Nowadays, this algorithm is employed in numerous tomographic reconstruction\n", "problems and often associated to regularization techniques. It is based on the iterative\n", "maximization of the log-likelihood function." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAHeCAYAAACok2NLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABAAUlEQVR4nO3debwcVZn/8c83NwsGAgHCmoRFWRSZUTAC6qAoqIAIzIyjoCgwKCMjuIsgjji4DO7gTxwmIzEuEEREBhUEXJDRASRAWEJAAiIJJCSsIpGsz++POhcqfbvv7bq3uqtv9/f9euWV7qrTp55zu08/VedUVSsiMDMzs+qMqToAMzOzXudkbGZmVjEnYzMzs4o5GZuZmVXMydjMzKxiTsZmZmYVczK2pkn6hKRvVR2HWdkk7SdpcQvq3UFSSBrbYP18SfuVvd2ySdpC0l2SnpeeXyHp6BZsZ7akzw6yPiTtVPZ2myXpfkkHpMcnSfpCWXU7GbdYevP+KukpSU9I+j9J75U06v72EfH5iHh31XHY6Cfp7ZLmSvqLpCXpy/3vqo6r3SLixRFxTdVxNOEUYHZE/BUgIg6KiO9UHFPV/ht4h6Qty6hs1CWEUerNETEJ2B44E/g4cF61IZlVQ9KHgbOAzwNbAdsB3wQOqzAsa0DSBOBo4PtVxzJc9UYmGo1WNCsingGuAN41knr6ORm3UUQ8GRGXAW8Djpa0O4CkayQ9e8Qp6RhJv809D0n/KumedIT9GUkvSEfZf5Z0kaTxqex+khZLOlnSsnTUcbikgyX9QdJjkj6Rym4taYWkzXPb2lPScknjauOX9GlJ30+P+4ffjpW0SNLj6Yj/5ZJuS6MA38i99gWSfiXpUUmPSDpf0uSa7d6S2vdDST/ID1dJOkTSvNzowt+W8qZYW0naBDgDeF9EXBIRT0fE6oj4SUR8LJWZIOksSQ+lf2elhJD/fH8k9/k+Nq3bW9JSSX257f29pNuGqrcmxo9Lurhm2dmSvt7fBknnpW0/KOmz/duU1Cfpy+kzfh/wpiH+Hvlhz0+nz/73Uz+4XdIukk5NbV0k6Q251x4raUEqe5+kf6mp++QU40OS3q3cEG/6W3xZ0gOSHpZ0rtIQdB17A09ExOJc3et9Z9Vst0jd9UyRdHVq128kbd9gO818b75P0j3APbnPzsclLQW+LWmMpFMk3Zu+my6StFmujndK+lNad1qdMK5hiPe4WU7GFYiI3wOLgX0LvOyNwMuAfYCTgZnAUcB0YHfgyFzZrYENgKnAp8iGU45Kr98X+DdJO0bEUrIP01tzr30ncGFErG4yrr2Bncl2MM4CTgMOAF4MvFXSa1I5Af8BbAu8KMX9aQBlOxI/BmYDmwFzgL/v34CkPYBZwL8AmwP/BVxW74vUOt4ryD6bPx6kzGlkn/OXAi8B9gI+mVu/NbAJ2ef7OOAcSZtGxA3A08DrcmXfDlzQZL39LgQOljQJsgRL1kf665kNrAF2AvYA3gD0J4X3AIek5TOAtwzSznreDHwP2BS4BbiS7Ht6KtlOzH/lyi5L29oYOBb4mqQ9U8wHAh8m64s7AfvVbOdMYBeyv8VOPPddUc/fAHcXaEORuut5B/AZYAowDzi/wGtrHU72HbVber412XfM9sDxwEmpzGvIvpseB84BkLQb8J9k34nbkn33TKupfwHZZ2nkIsL/WvgPuB84oM7y64HT0uNrgHfn1h0D/Db3PIBX5Z7fBHw89/wrwFnp8X7AX4G+9HxSev3eNa8/PD1+G/C79LgPWArs1aAtnwa+nx7vkOqdmlv/KPC23PMfAR9sUNfhwC3p8auBBwHl1v8W+Gx6/J/AZ2pefzfwmqrfX/8r3B/eASwdosy9wMG5528E7k+P+z/fY3PrlwH7pMefBWalx5PIkvP2Tda7OLfut8C70uPXA/emx1sBK4Hn5coeCfw6Pf4V8N7cujekfjK2QVvvJ30/pP51dW7dm4G/1OnLkxvUdSnwgfR4FvAfuXU7pdfuRLZj/DTwgtz6VwB/bFDvaWQ76Pll15D7zsotL1R3ndfPzm8L2AhYC0xPzwPYqV4M1P/efF3u+X7AKmCD3LIFwP6559sAq4GxZDsQ+Vg2TK8/ILdsZ2BtGX1jRGPmNiJTgccKlH849/ivdZ5vnXv+aESsza2r9/qN0uP/Ac6VtCOwK/BkZEfuZcW1EYCkrYCzyY7MJ5Ht7T+eym0LPBjp050syj3enmxY/6TcsvHpdTa6PEo2DDk2ItY0KLMt8Kfc8z+x/nv9aM1rV/Dc5/kC4P8knQD8A3BzRPTXNVS9eReQJdnvsv7R9fbAOGCJpP6yY3ju87ot639289trRm3/eaROX94IeELSQcDpZEehY4CJwO25OObm6srHtEUqe1OuDSLbGa/ncbI+24yiddfzbKwR8RdJjzHw71q4rmR5ZHO9/bYHfixpXW7ZWrKdrvW2GRFPS3q0pr5JwJPDiGsAD1NXQNLLyZJx//zG02Qf4H5bD3hRi6QP5kVkw9jvJBsia4XPk+2p/k1EbJy2199blwBTleu9ZMPY/RYBn4uIybl/EyNiTotitda5juzI8vBByjxE9iXZb7u0bEgRcSdZAjyI9ZNo0Xp/COwnaRrZlEl/PYtS/FNyn8WNI+LFaf0S1v/sbtdM3EWlKZofAV8GtoqIycDlrN+n8kOq+ZgeIUvsL861YZOI2Ij6biNL+M0oWnc9z8YqaSOyYeV671Mz35u1P0tY+3wRcFDNd8sGEfEgNe+lpIlkQ9V5LwJuHawxzXIybiNJG0s6hGxO6vsR0b8XOw/4B0kT0wkWx7U5tO+SDfEcSuuS8SSyIbcnJU0FPpZbdx3Z3uiJksZKOoxsPq/ffwPvVXaCjiRtKOlN/XN6NnpExJNkw3/nKDuxcKKkcZIOkvTFVGwO8Ell17ZOSeWLnMl7AfABsumPH+aWN11vRCwnGwb9NtkQ64K0fAlwFfCV1J/HKDs5sf/ciIuA90uaJmlTskuCWmE8MAFYDqxJR8lvyK2/CDhW0otSEvm3XNvWkfWpryldliNpqqQ3NtjW74HJqd8Oahh113OwpL9L55J8Brg+IuodFc9j5N+b5wKf6z9JLH02+s/qvxg4JBfLGQzMma8hO6N6xJyM2+Mnkp4i2ws7Dfgq2QkX/b5GNhfxMPAdRnbCQmER8TtgHesP6ZXt34E9yYZ0fgZcktv+KrIhxeOAJ8iOmn9KdgRCRMwlOzHmG2RDZgvJdh5sFIqIr5CdXPRJsmSyCDiRbM4TsnnfuWRHZLcDN6dlzZpD9iX5q4h4JLe8aL0XkJ0AdUHN8neRJcM7yT6PF5PNNUKWiK4kO1q6mdznvEwR8RTwfrKk+zjZKMBlufVXAF8Hfk3WX65Pq1am/z/ev1zSn4FfkE1T1dvWKrK53KPqrZe0r6S/5BY1rFvSdsquLR9sxOACsuH3x8hOOq27Xcr53jyb7O92VfqOvp7shC8iYj7wvhTPErK/c/6M8g2Ag9O2R0zrT9NZr5L0K+CCiOiIO2xJugE4NyK+XXUsZqOdpBcBdwATBpmrH+z1WwD/C+wR6cYfvS6dwzI9Ik4upT4nY0tz2FeTfbCeqiiG15CdIf0I2Rm35wLPT8OCZlaQpL8nm0eeSHb0ti4iDq80KGuo0DC1pFnKLj6/o2b5ScruWzo/N+9jo4Ck75ANI32wqkSc7Eo2tPcE8BHgLU7EreX+3PX+heyyr3vJzsk4odpwbDCFjowlvZrsJJzvRkT/3aNeSzYP+qaIWClpy4hY1pJozaw07s9mnaPQkXFEXMvAa2NPAM6MiP6TbdxxzUYB92ezzlHGTT92AfaV9DngGeCjEXFjvYKSjie7BRl99L1sIhuXsHmz7rDLy55fd/lNN930SERs0a4wGE5/1riXbTh+s4GFWn1OStHqNXSR9csXfEEntbeX2god196d/2b6gGWD9eUykvFYsouy9wFeDlwk6flRZ/w7ImaS3VOZjbVZ7K39S9i8WXe4eu4P6y6X1KrLzeoZVn/eZIOt45XT3jmwtrVrBy4r09p1Q5fJ6yt4NWdfkRtH0Vnt7aW2Qse194q5Xx2wbLC+XMZ1xouBSyLze7LrVaeUUK+ZtZ/7s1kFykjGlwKvBZC0C9nF8I8M9gIz61iX4v5s1naFhqklzSH75YspkhaT3SVlFjArXR6xCji63pCWmXUW92ezzlEoGUfEkQ1WNbpdmZl1KPdns87RcT+heOVD8wqVf+O2L21JHO3SS+3tpbZC8faOahIxbuDXyYJTip0E/oI5xU6SmXDf8kLl68U4qIIn+XRSezutrc+/sNgJWRvcW+yquk5rb1H+oQgzM7OKORmbmZlVzMnYzMysYk7GZmZmFXMyNjMzq5iTsZmZWcWcjM3MzCrmZGxmZlYxJ2MzM7OKORmbmZlVzMnYzMysYh13b+rRfj/ionqpvb3UVije3qsL/pZ6R2lwb+pNbh9XqJqVk4sdH0wo+gPxBcsXvd/x5Ns6qL0d1tZVk4vdh3uDTntvby3W3qJ8ZGxmZlYxJ2MzM7OKFU7GkmZJWpZ+fLx23UckhaQp5YRnZq3ivmzWOYZzZDwbOLB2oaTpwBuAB0YYk5m1x2zcl806QuFkHBHXAo/VWfU14GQgRhqUmbWe+7JZ5yhlzljSYcCDEXHrEOWOlzRX0tzVrCxj02ZWomb7cir7bH9etebpNkRn1r1GfGmTpInAJ8iGtQYVETOBmQAbazPvdZt1kCJ9Gdbvz5tM3Nb92WwEyjgyfgGwI3CrpPuBacDNkrYuoW4zax/3ZbOKjPjIOCJuB7bsf5468YyIeGSkdZtZ+7gvm1VnOJc2zQGuA3aVtFjSceWHZWat5r5s1jkKHxlHxJFDrN9h2NGYWdu4L5t1jo67N7WZjULr1jFmxTMDFk+9bFFrt7u22A29tbbY/ZG1ek2h8tv+pHPa20tthTa096cF2/v/ihX37TDNzMwq5mRsZmZWMSdjMzOzijkZm5mZVczJ2MzMrGJOxmZmZhVzMjYzM6uYk7GZmVnFnIzNzMwq5mRsZmZWMSdjMzOzivne1GY2chIxrs7XyRi1drMUu79w3RgH09dXrHwHtbeX2gqjv70+MjYzM6uYk7GZmVnFCiVjSbMkLZN0R27ZlyTdJek2ST+WNLn0KM2sdO7PZp2j6JHxbODAmmVXA7tHxN8CfwBOLSEuM2u92bg/m3WEQsk4Iq4FHqtZdlVE9M+0Xw9MKyk2M2sh92ezzlH2nPE/A1c0WinpeElzJc1dzcqSN21mJWu6P69as6KNYZl1n9KSsaTTgDXA+Y3KRMTMiJgRETPGMaGsTZtZyYr25/FjJ7YvOLMuVMp1xpKOAQ4B9o+IKKNOM6uG+7NZ+404GUs6EDgZeE1EeKzKbBRzfzarRtFLm+YA1wG7Slos6TjgG8Ak4GpJ8ySd24I4zaxk7s9mnaPQkXFEHFln8XklxWJmbeT+bNY5fG9qMxu5Bvembu3dfIF1Bae0C96PuOj9jjuqvb3UVhj17fXtMM3MzCrmZGxmZlYxJ2MzM7OKORmbmZlVzMnYzMysYk7GZmZmFXMyNjMzq5iTsZmZWcWcjM3MzCrmZGxmZlYxJ2MzM7OKddy9qa98aF6h8m/c9qUtiaNdeqm9vdRWKN7eUW3dOsaseGbA4jtP3rJQNTudv7pQ+fF/XFaovNauLVZ+9ZpC5e/8WMH2XtC69vZSW6Hz2luUj4zNzMwq5mRsZmZWscLJWNIsScsk3ZFbtpmkqyXdk/7ftNwwzaxs7stmnWM4R8azgQNrlp0C/DIidgZ+mZ6bWWebjfuyWUconIwj4lrgsZrFhwHfSY+/Axw+srDMrNXcl806R1lnU28VEUvS46XAVvUKSToeOB5gAyaWtGkzK1FTfRlq+vPYSW0Izax7lX4CV0QEEA3WzYyIGRExYxwTyt60mZVosL6c1j/bn8f3eefabCTKSsYPS9oGIP1f7AIxM+sU7stmFSgrGV8GHJ0eHw38T0n1mll7uS+bVWA4lzbNAa4DdpW0WNJxwJnA6yXdAxyQnptZB3NfNuschU/giogjG6zaf4SxmFkbuS+bdY6Ouzf1aL8fcVG91N5eaisUb+/V61oTR1tIxLiBXyeTbyv2FbNqcrHNju8rNrhXL8ZB9fUVKj759s5pb6vbuultxcqv3KxQccY/0FnvbdH2FuXbYZqZmVXMydjMzKxiTsZmZmYVczI2MzOrmJOxmZlZxZyMzczMKuZkbGZmVjEnYzMzs4o5GZuZmVXMydjMzKxiTsZmZmYV67h7U5tZ99jmqiVVh9BWvdTera9eWnUIbdXq9vrI2MzMrGJOxmZmZhUrLRlL+pCk+ZLukDRH0gZl1W1m7eX+bNZepSRjSVOB9wMzImJ3oA84ooy6zay93J/N2q/MYeqxwPMkjQUmAg+VWLeZtZf7s1kblZKMI+JB4MvAA8AS4MmIuKq2nKTjJc2VNHc1K8vYtJmVbDj9edWaFe0O06yrlDVMvSlwGLAjsC2woaSjastFxMyImBERM8YxoYxNm1nJhtOfx4+d2O4wzbpKWcPUBwB/jIjlEbEauAR4ZUl1m1l7uT+btVlZyfgBYB9JEyUJ2B9YUFLdZtZe7s9mbVbWnPENwMXAzcDtqd6ZZdRtZu3l/mzWfqXdDjMiTgdOL6s+M6tO4f4sEeMGfp2oxJjqWhfFyvf1FSper02D6aj29lJbYdS313fgMjMzq5iTsZmZWcWcjM3MzCrmZGxmZlYxJ2MzM7OKORmbmZlVzMnYzMysYk7GZmZmFXMyNjMzq5iTsZmZWcWcjM3MzCpW2r2pzayHrVvHmBXPDFy+dm1rt7t2XaHiKhiPVq8pVL6T2ttLbYXR314fGZuZmVXMydjMzKxiTsZmZmYVKy0ZS5os6WJJd0laIOkVZdVtZu3l/mzWXmWewHU28POIeIuk8cDEEus2s/ZyfzZro1KSsaRNgFcDxwBExCpgVRl1m1l7uT+btV9Zw9Q7AsuBb0u6RdK3JG1YW0jS8ZLmSpq7mpUlbdrMSla4P69au6L9UZp1kbKS8VhgT+A/I2IP4GnglNpCETEzImZExIxxTChp02ZWssL9eXyfR7HNRqKsZLwYWBwRN6TnF5N1ZjMbfdyfzdqslGQcEUuBRZJ2TYv2B+4so24zay/3Z7P2K/Ns6pOA89OZl/cBx5ZYt5m1l/uzWRuVlowjYh4wY6T1XPnQvELl37jtS0e6yUr1Unt7qa1QvL2dpHB/lohxA79OFpyyRaHtvmBOsfv/TrhveaHy9WIcVF9foeILPl6wvRe2rr291FbovPYW5TtwmZmZVczJ2MzMrGJOxmZmZhVzMjYzM6uYk7GZmVnFnIzNzMwq5mRsZmZWMSdjMzOzijkZm5mZVczJ2MzMrGJOxmZmZhUr84ciSjHa70dcVC+1t5faCsXbe/W61sRRpcnzxhUqv3ojFSrfab+KPvnWYu1dtfHobW8vtRVg03mtTZc+MjYzM6uYk7GZmVnFSk3Gkvok3SLpp2XWa2bt5/5s1j5lHxl/AFhQcp1mVg33Z7M2KS0ZS5oGvAn4Vll1mlk13J/N2qvMI+OzgJOBhueESjpe0lxJc1ezssRNm1nJzqJAf161ZkXbAjPrRqUkY0mHAMsi4qbBykXEzIiYEREzxnXcietmBsPrz+PHTmxTdGbdqawj41cBh0q6H7gQeJ2k75dUt5m1l/uzWZuVkowj4tSImBYROwBHAL+KiKPKqNvM2sv92az9fJ2xmZlZxUq/v1dEXANcU3a9ZtZ+7s9m7dFx96Y2s+6xzVVLqg6hrXqpvS1vq4rdy7rVtr56aUvr9zC1mZlZxZyMzczMKuZkbGZmVjEnYzMzs4o5GZuZmVXMydjMzKxiTsZmZmYVczI2MzOrmJOxmZlZxZyMzczMKuZkbGZmVjEnYzMzs4o5GZuZmVXMydjMzKxipSRjSdMl/VrSnZLmS/pAGfWaWfu5P5u1X1m/Z7wG+EhE3CxpEnCTpKsj4s6S6jez9nF/NmuzUo6MI2JJRNycHj8FLACmllG3mbWX+7NZ+5V1ZPwsSTsAewA31Fl3PHA8wAZMLHvTZlaypvvz2I3bG5hZlyn1BC5JGwE/Aj4YEX+uXR8RMyNiRkTMGMeEMjdtZiUr0p/Hj/XOtdlIlJaMJY0j67jnR8QlZdVrZu3n/mzWXmWdTS3gPGBBRHy1jDrNrBruz2btV9aR8auAdwKvkzQv/Tu4pLrNrL3cn83arJQTuCLit4DKqMvMquX+bNZ+vgOXmZlZxZyMzczMKuZkbGZmVjEnYzMzs4o5GZuZmVXMydjMzKxiTsZmZmYVczI2MzOrmJOxmZlZxZyMzczMKuZkbGZmVjEnYzMzs4qV8kMRZbryoXmFyr9x25e2JI526aX29lJboXh7u9GCj00pVP4Fc9YWKj/hj8sLlW+1wu29sGB77+uc9ra8rX98pFD5Vlvw0S1aWr+PjM3MzCrmZGxmZlax0pKxpAMl3S1poaRTyqrXzNrP/dmsvUpJxpL6gHOAg4DdgCMl7VZG3WbWXu7PZu1X1pHxXsDCiLgvIlYBFwKHlVS3mbWX+7NZmykiRl6J9BbgwIh4d3r+TmDviDixptzxwPHp6a7A3QU2MwXorNPrWquX2ttLbYXi7d0+Ilp7KmeO+3Ppeqmt4PYOpmFfbuulTRExE5g5nNdKmhsRM0oOqWP1Unt7qa3QPe11f25OL7UV3N7hKmuY+kFgeu75tLTMzEYf92ezNisrGd8I7CxpR0njgSOAy0qq28zay/3ZrM1KGaaOiDWSTgSuBPqAWRExv4y6c4Y1HDaK9VJ7e6mt0OHtdX8uXS+1FdzeYSnlBC4zMzMbPt+By8zMrGJOxmZmZhXrmGQs6QOS7pA0X9IHc8tPknRXWv7FtGwHSX+VNC/9O7eywIepXnsl/SDXpvslzcuVPzXdmvBuSW+sKu7hKNLWLn5vXyrp+tSmuZL2Sssl6evpvb1N0p6VBl+SXurPvdSXwf05LSu/P0dE5f+A3YE7gIlkJ5X9AtgJeG16PCGV2zL9vwNwR9Vxl93emjJfAT6VHu8G3ApMAHYE7gX6qm5Hi9rale8tcBVwUCpzMHBN7vEVgIB9gBuqbkML/wZd1597qS8Ps72j9r0drL2t6M+dcmT8IrKgV0TEGuA3wD8AJwBnRsRKgIhYVmGMZWrUXiDbuwLeCsxJiw4DLoyIlRHxR2Ah2S0LR4OibR3tGrU3gI1TmU2Ah9Ljw4DvRuZ6YLKkbdoddMl6qT/3Ul8G9+eW9edOScZ3APtK2lzSRLK9i+nALmn5DZJ+I+nludfsKOmWtHzfKoIegUbt7bcv8HBE3JOeTwUW5dYvTstGg6Jthe58bz8IfEnSIuDLwKmp/Gh+bxvppf7cS30Z3J9b1p/bejvMRiJigaQvkB36Pw3MA9aSxbcZ2eH+y4GLJD0fWAJsFxGPSnoZcKmkF0fEnytpQEGDtLffkXTJnuUw2tqt7+0JwIci4keS3gqcBxxQWaAt1Ev9uZf6Mrg/08L+3ClHxkTEeRHxsoh4NfA48AeyvYpL0iH/74F1wJQ0xPNoet1NZPMuu1QV+3A0aC+SxpINg/wgV3xU356wSFu7+L09GrgkFfkhzw1Njur3tpFe6s+91JfB/ZlW9efhTmyX/Y/nTubYDrgLmAy8FzgjLd+F7PBfwBakkx6A56fGblZ1G0ba3vT8QOA3NWVfzPonfdzH6Drpo0hbu/K9BRYA+6Xl+wM3pcdvYv0TPn5fdfwt/Bt0ZX/upb48jPaO6ve2UXtb0Z87Ypg6+ZGkzYHVwPsi4glJs4BZku4AVgFHR0RIejVwhqTVZHvX742Ix6oLfVgGtDctP4KaYa2ImC/pIuBOYE0qnx8a6nRNtxXoyvdW0nuAs9PRwzM899ODl5PNQy0EVgDHVhFwC/RSf+6lvgzuzy3pz74dppmZWcU6Zs7YzMysVzkZm5mZVczJ2MzMrGJOxmZmZhVzMjYzM6uYk7GZmVnFnIzNzMwq5mRsZmZWMSdjMzOzijkZm5mZVczJ2MzMrGJOxjYkSZ+W9P0G6/aVdHe7YzIz6yZOxiWS9HZJcyX9RdISSVdI+rsOiOt+SS35IfuI+N+I2LUVdZuZ9Qon45JI+jBwFvB5YCuy3778JnDYMOoa8NOW9ZaZmVl3cDIugaRNgDPIfuvykoh4OiJWR8RPIuJjqcwESWdJeij9O0vShLRuP0mLJX1c0lLg22lo+GJJ35f0Z+AYSZtIOi8ddT8o6bOS+nJxvEfSAklPSbpT0p6Svke2Y/CTdMR+ciq7j6T/k/SEpFsl7ZerZ0dJv0n1XA1MGaTt+0lanHt+v6SPSbpN0tMp3q3SKMFTkn4hadNc+R9KWirpSUnXSnpxbt3mkn4i6c+Sbkzt/W1u/QslXS3pMUl3S3rrsN9EsxKlz/vRFcewg6QYzo68pOelvvekpB82+ZprJL27eKR162o4mlf7ndNuI/m7DsbJuByvADYAfjxImdOAfYCXAi8B9gI+mVu/NbAZsD3P/VD1YcDFwGTgfGA22Q+S7wTsAbwBeDeApH8CPg28C9gYOBR4NCLeCTwAvDkiNoqIL0qaCvwM+Gza5kfJfkB7i7TdC4CbyJLwZ4CiXyr/CLwe2AV4M3AF8AlgC7LP3PtzZa8Adga2BG5O7ex3DvB0+tscnY9D0obA1SnWLcl+2PybknYrGKtVIH3Z/jXtIC6VNFvSRk2+trQv/VaJiIMi4jsjqWOwczXa4C1kI3ybR8Q/1a6sOLau5GRcjs2BRyJizSBl3gGcERHLImI58O/AO3Pr1wGnR8TKiPhrWnZdRFwaEevIEuzBwAfTkfcy4GtkSQiypPzFiLgxMgsj4k8NYjkKuDwiLo+IdRFxNTAXOFjSdsDLgX9LsVwL/KTg3+P/RcTDEfEg8L/ADRFxS0Q8Q7bDskd/wYiYFRFPRcRKsp2Jl6QRgD6ypH56RKyIiDuB/JfbIcD9EfHtiFgTEbcAPwIGfHFYx3pzRGxEtoO6B3BqGZW2ekpHmW7/7twe+MMQ32mjTidPAXb7B6pdHgWmDPGmbgvkk+Of0rJ+y1OyyluUe7w9MA5YkoaWnwD+i+yoEGA6cG+T8W4P/FN/PamuvwO2STE9HhFP18RaxMO5x3+t83wjAEl9ks6UdG8air8/lZlCdhQ9lvX/BrV/j71r2vAOsqNoG0UiYilwJVlSBhpPo0j6HLAv8I10VP2NtDwkvU/SPcA9adl7JC1M0xiXSdqWBiQdI+l3kr6RhmbvkrR/bv01kj4n6XfACuD5kl6Zpk+eTP+/sqb8u3PP/1nZFNLjkq6UtH1u3Ytz0y0PS/qEpAPJRpPeltp5ayrbcKoq9acvS3pE0n3Amwb7u0t6UYrzCUnzJR2alv878Kncto+reV3d2JLt09/xKUlXSZqSe13DqbEGXq5suu1xSd+WtEGDdoSknXLPZ0v6bHo84inAon/X4XIyLsd1wErg8EHKPESWQPptl5b1izqvyS9blLYxJSImp38bR8SLc+tf0GDbtXUvAr6Xq2dyRGwYEWcCS4BNlQ0D52NthbeTDcUfAGwC7JCWC1hONiQ/LVd+ek0bflPTho0i4oQWxWotImkacBCwMD1vOI0SEaeRjbacmN7vE3NVHQ7sDewm6XXAfwBvJdvJ/BNw4RCh7E22QzsFOB24RNJmufXvJJtCmgQ8lWL8OtnI2FeBn0navE77DiNLXv9AtpP5v8CctG4S8Avg52Q7wjsBv4yIn5OdDPqD1M6XpOpm02CqCngP2YjRHsAMsqHmuiSNIxvxuopsh/4k4HxJu0bE6TXbPi//2kFig6xPH5vqHE/23g36njaKkWzn+o1k32u7sP60XhEjmgKkwN91JJyMSxART5LtSZ4j6XBJEyWNk3SQpC+mYnOAT0raIu0tfgpoes4lIpaQdZyvSNpY0hhJL5D0mlTkW8BHJb1MmZ1ye98PA8/PVfd94M2S3pj2+jZIe5DT0tD2XODfJY1XdmnWm4f5pxnKJLIdjEeBiWQdvL+9a4FLgE+nv+cLyebD+/0U2EXSO9Pfepykl0t6UYtitfJdKukpsh2rZWQJEAaZRhmivv+IiMfSNM87gFkRcXOaAjkVeIWkHQZ5/TLgrHTy5Q+Au1n/KGh2RMxPQ7dvAO6JiO+laZI5wF3U7yvvTbEtSK/9PPDS1D8PAZZGxFci4pk0ZXNDveAkbcXgU1VvTfEviojHyHZGGtmHbITqzIhYFRG/IutTRw7ymmZ8OyL+kN6Di3hutGM47+k3cm353AhiG+kUYJG/67A5GZckIr4CfJhs72052RfMicClqchnyT58twG3k52s9NmCm3kX2d7mncDjZHt226Tt/5DsA3sB2V77pWR7g5B9eD6Zhoc+GhGLyPYMP5GL9WM893l4O9lRwmNkX5DfLRhns75LdsTyYGrT9TXrTyQ7Yl4KfI9sh2YlQEQ8RfaFeATZCMNS4AvAhBbFauU7PCImAfsBL+S5s/YHm0YZTH4aY71poYj4C9lO31RlN6r5S/o3P/eaByMiP4pUO5XUsP5c+al14toeODvXlsfIRn+mUnx6abCpqm1rYhxsemlbYFFKRkPFX8TS3OMVpCkphvee1ral4TTDEEY6BVjk7zpsHTFx3S0i4nzWPxs4v+4ZsrOI319n3TWsPxxLRHy6TrkngRPSv3rbOBc4t87y/wH+p2bZDcBrasumdfeRzcsNqTb2iNihZv1RNc+/RXYU3/8FWXsd9ndzZZeTOzKR9AVgcW597ZGLjUIR8RtJs4Evkw0190+jvKfRS5pYvt60UJp22Zws4d7Pc0kib6ok5RLydsBlzdSfK//zOvUuAj6Xvh/Wk46Ojxj4kgHb66+nf6qq3olVS1h/Kmew6aWHgOmSxuQS8nbAHwZ5zWCxDWWo97Se2rY81KDcCrKRtX5bk/ueoNgU4Ej/rsPmI2PrWMquI/7bNOy+F3Acg18+ZqPXWcDrJb2EQaZRUtnaaZd65gDHSnqpsuv5P092Vv/9g7xmS+D9acrjn4AXAZc3KHs52TTJ2yWNlfQ2YDeyod5a5wKnKl1Dn04W6j/r/6fANpI+qOxeBJMk7Z1r5w5KZ243MVV1UYp/mrJr+U8ZpK03kCWxk1N79yMbYh9qXr3ferE1Yaj3tJ73pbZsRnZp6A8alJsHvD3VeyANDjIaKfnvOmxD/iElzZK0TNIdDdZL0teVnbV4m6Q9yw/TetQksnnjp8k64leoOcK37pBGQb4LfKqJaZSzgbcoO8v26w3q+wXwb2SXuy0hOwmo0RFovxvIrnl/hGzK5y0R8WiD+h8lm+/9CNnw98nAIRHxSJ2yPyabQrlQ2dm7d5CdsNY/3fJ6skS4lOxM8Neml/bfbONRSTenxw2nqoD/Jjsr/VayabBLGjU0IlalbR6U2vtN4F0RcVej19SoF1tDTbyn9VxAliTvIxvKbzSt9wGytjxBdq7ApUNGP1Apf9eR0PpTJHUKSK8G/gJ8NyJ2r7P+YLIz8Q4mm2c8OyL2ri1nZtapJB0DvDsiSrmXvKRrgW9FRKvOt7AuM+SRcWQ3fXhskCKHkSXqiIjrgcmShjrRwsysK0maSDaM/seqY7HRo4wTuKay/plmi9OyJbUFJR1Pus5rww03fNkLX/jCEjZv1t1uuummRyJisOsxh0XSLLKh1mX5US9JJwHvA9YCP4uIk8vedreStCXZ9dI/AX47RHGzZ7X1bOqImAnMBJgxY0bMnTu3nZs3G5UkteRSCrIbHXyD3Bnskl5LNtr1kohYmZJL14uI2WR/j5HWs4zsulWzQso4m/pB1j/te1paZmYdrMEU1AlkN4Lov557WdsDM+tBZRwZXwacKOlCshO4nkynipvZ6LMLsK+ye0A/A3w0Im6sVzA/7dQ3ZvzLJk4c+Eub0adCGx+zcm2xaNcVvNy1WDjEuL6hC+XLd1J7e6mt0PL2omIb2GWXgbfJH2zKachkLGkO2R1ypij7DcnTye5W0n+TicvJzqReSHbd2rGFIjazTjKW7M5t+5D9etdFkp4fdS67yE87bTxpauz10n8dUNnKzcYX2viG9z5RqLxW1N5YaXAxttgX8JqtNilUvpPa2/K2blqwrfc9Uah8q9/b1VsXa2+MLTaQ/KtfDfwRssGmnIZMxhEx6P1AUyd9XzPBmVnHWwxckvr17yWtI7tN5fJqwzLrbr4Dl5nlXUq66YSkXchuhDDgRhZmVi7fm9qsRzWYgpoFzEp33FsFHF1viNrMyuVkbNajBpmCOqrBcjNrESdjMxuxkFjzvIEn0Dw8o9hXzLS/1PsxpcbGL1pdqDzjisWzZsNi5Ze9rFj5qa1sb6vbWvC9nfp0Z723aycWK//obsVOWCvKc8ZmZmYVczI2MzOrmJOxmZlZxZyMzczMKuZkbGZmVjEnYzMzs4o5GZuZmVXMydjMzKxiTsZmZmYVczI2MzOrmJOxmZlZxXxvajMbMa0Nxj+5asDyTe8aV6yiDvt9qPGPryxUftO7Cn6ldlB7e6mtMIz23qMWRZLxkbGZmVnFnIzNzMwq5mRs1qMkzZK0TNIdddZ9RFJImlJFbGa9xsnYrHfNBg6sXShpOvAG4IF2B2TWq5yMzXpURFwLPFZn1deAk+m4U27MupeTsZk9S9JhwIMRcWsTZY+XNFfS3NVrnm5DdGbdy5c2mRkAkiYCnyAboh5SRMwEZgJsvOFUH0WbjYCPjM2s3wuAHYFbJd0PTANulrR1pVGZ9QAfGZsZABFxO7Bl//OUkGdExCOVBWXWI3xkbNajJM0BrgN2lbRY0nFVx2TWq3xkbNajIuLIIdbv0KZQzHqek7GZjZjWrmXMEwPPqN7shqeKVRQFzwMrWL7oWWb12jSYTX/fOe0t3NYnVxQq30lthda3d8PFywtuoRgPU5uZmVXMydjMzKxiTsZmZmYVayoZSzpQ0t2SFko6pc767ST9WtItkm6TdHD5oZqZmXWnIZOxpD7gHOAgYDfgSEm71RT7JHBRROwBHAF8s+xAzczMulUzR8Z7AQsj4r6IWAVcCBxWUyaAjdPjTYCHygvRzMysuzWTjKcCi3LPF6dleZ8GjpK0GLgcOKleRfkbyy9f3trTxM3MzEaLsk7gOhKYHRHTgIOB70kaUHdEzIyIGRExY4sttihp02ZmZqNbM8n4QWB67vm0tCzvOOAigIi4DtgAmFJGgGZmZt2umWR8I7CzpB0ljSc7QeuymjIPAPsDSHoRWTL2OLSZmVkThkzGEbEGOBG4ElhAdtb0fElnSDo0FfsI8B5JtwJzgGMiit77zMzMrDc1dW/qiLic7MSs/LJP5R7fCbyq3NDMbNSQYNzAr5MYs65YPUX34dcVrL9OjIOSChXvqPa2uq0Fy4/295axrf0pB9+By8zMrGJOxmZmZhVzMjbrUZJmSVom6Y7csi9Juivd1vbHkiZXGKJZz3AyNutds4EDa5ZdDeweEX8L/AE4td1BmfUiJ2OzHhUR1wKP1Sy7Kl1BAXA92X0FzKzFnIzNrJF/Bq5otDJ/e9tVa1e0MSyz7uNkbGYDSDoNWAOc36hM/va24/smti84sy7U2gunzGzUkXQMcAiwv2/eY9YeTsZm9ixJBwInA6+JCI89m7WJh6nNepSkOcB1wK6SFks6DvgGMAm4WtI8SedWGqRZj/CRsVmPiogj6yw+r+2BmJmTsZmNXIwdw5pNB57EtXaDvkL1TFj6l2IbXrmqWPmCVm9W7MS0dRNGb3t7qa1QvL1jNt6gRZGk+ltau5mZmQ3JydjMzKxiTsZmZmYVczI2MzOrmJOxmZlZxZyMzczMKuZkbGZmVjEnYzMzs4o5GZuZmVXMydjMzKxiTsZmZmYV872pzWzkJNaNH7hvv+SVEwpVM/XXawuVH7/4r4XKF/1x5hhX7HilcHuvWVeo/PhFzbe3l9oKrW/vw3sWu5d1UT4yNjMzq5iTsZmZWcWcjM16lKRZkpZJuiO3bDNJV0u6J/2/aZUxmvUKJ2Oz3jUbOLBm2SnALyNiZ+CX6bmZtZiTsVmPiohrgcdqFh8GfCc9/g5weDtjMutVPpvazPK2iogl6fFSYKtGBSUdDxwPMGHC5NZHZtbFmjoylnSgpLslLZRUd9hK0lsl3SlpvqQLyg3TzNotIoJBrhiJiJkRMSMiZowfv2EbIzPrPkMeGUvqA84BXg8sBm6UdFlE3JkrszNwKvCqiHhc0patCtjMWuphSdtExBJJ2wDLqg7IrBc0c2S8F7AwIu6LiFXAhWTzSnnvAc6JiMcBIsId2Gx0ugw4Oj0+GvifCmMx6xnNJOOpwKLc88VpWd4uwC6Sfifpekm1Z2gC2RyTpLmS5i5fvnx4EZtZKSTNAa4DdpW0WNJxwJnA6yXdAxyQnptZi5V1AtdYYGdgP2AacK2kv4mIJ/KFImImMBNgxowZRe9eZmYliogjG6zav62BmFlTyfhBYHru+bS0LG8xcENErAb+KOkPZMn5xlKiNLOOpjXrGPf4MwOWT7m92P6+1hXcRx9T8OrMguXHPrGyUPnC7V1b7H7NheLv6ytU9ahuK7S8vZvfoULli2qmtTcCO0vaUdJ44AiyeaW8S8mOipE0hWzY+r7ywjQzM+teQybjiFgDnAhcCSwALoqI+ZLOkHRoKnYl8KikO4FfAx+LiEdbFbSZmVk3aWqcISIuBy6vWfap3OMAPpz+mZmZWQG+HaaZmVnFnIzNzMwq5mRsZmZWMSdjMzOzijkZm5mZVczJ2MzMrGJOxmZmZhVzMjYzM6tYWT8UYWa9LAKtWjNg8UYLHitcT6HiYwreL7iv2PHHmJWrC5XvqPYW/NuM6rZCy9u7wT1PFipflI+MzczMKuZkbGZmVjEnYzMbQNKHJM2XdIekOZI2qDoms27mZGxm65E0FXg/MCMidgf6yH461cxaxMnYzOoZCzxP0lhgIvBQxfGYdTUnYzNbT0Q8CHwZeABYAjwZEVfVlpN0vKS5kuauWrui3WGadRUnYzNbj6RNgcOAHYFtgQ0lHVVbLiJmRsSMiJgxvm9iu8M06ypOxmZW6wDgjxGxPCJWA5cAr6w4JrOu5mRsZrUeAPaRNFGSgP2BBRXHZNbVnIzNbD0RcQNwMXAzcDvZ98TMSoMy63K+HaaZDRARpwOnVx2HWa/wkbGZmVnFnIzNzMwq5mRsZmZWMSdjMzOzijkZm5mZVczJ2MzMrGJOxmZmZhVzMjYzM6uYk7GZmVnFnIzNzMwq5mRsZmZWsabuTS3pQOBsoA/4VkSc2aDcP5LdYP7lETG3tCjNrKPF2DGs2XTgbxqv3Gx8oXo2vPeJQuW1YnWh8lGoNKzeepNC5Vdt2jnt7aW2QuvbG1ttVHALxQx5ZCypDzgHOAjYDThS0m51yk0CPgDcUHaQZmZm3ayZYeq9gIURcV9ErAIuBA6rU+4zwBeAZ0qMz8zMrOs1k4ynAotyzxenZc+StCcwPSJ+NlhFko6XNFfS3OXLlxcO1szMrBuN+AQuSWOArwIfGapsRMyMiBkRMWOLLbYY6abNrEUkTZZ0saS7JC2Q9IqqYzLrZs2cwPUgMD33fFpa1m8SsDtwjSSArYHLJB3qk7jMRq2zgZ9HxFskjQcGnp1lZqVpJhnfCOwsaUeyJHwE8Pb+lRHxJDCl/7mka4CPOhGbjU6SNgFeDRwDkM4VWVVlTGbdbshh6ohYA5wIXAksAC6KiPmSzpB0aKsDNLO22xFYDnxb0i2SviVpw9pC+XNAVq9+uv1RmnWRpq4zjojLgctrln2qQdn9Rh6WmVVoLLAncFJE3CDpbOAU4N/yhSJiJjATYONJU4te5mlmOb4Dl5nVWgwsjoj+ewZcTJaczaxFnIzNbD0RsRRYJGnXtGh/4M4KQzLrek0NU5tZzzkJOD+dSX0fcGzF8Zh1NSdjMxsgIuYBM5ouL7F2g74By5ftWewrZtqfB5wnNqjxiwqe5D2m2GDgugkD2zSYou2d+lSx+x1PeKBAe/uKxV64rXt0UFuh5e19ZPcJhcoX5WFqMzOzijkZm5mZVczJ2MzMrGJOxmZmZhVzMjYzM6uYk7GZmVnFnIzNzMwq5mRsZmZWMSdjMzOzijkZm5mZVczJ2MzMrGK+N7WZjdiYNesY/+hfByzf4tZi9//V6nVlhdRgAypUfNxjzxQqX7S9Y1atLVS+lQq39baCbW31e1tQ0fZOub1FgSQ+MjYzM6uYk7GZmVnFnIzNzMwq5mRsZnVJ6pN0i6SfVh2LWbdzMjazRj4ALKg6CLNe4GRsZgNImga8CfhW1bGY9QInYzOr5yzgZKDh9SiSjpc0V9LcVWtWtC0ws27kZGxm65F0CLAsIm4arFxEzIyIGRExY/zYiW2Kzqw7ORmbWa1XAYdKuh+4EHidpO9XG5JZd3MyNrP1RMSpETEtInYAjgB+FRFHVRyWWVdzMjYzM6uY701tZg1FxDXANU0URCtXD1i80fxHim6wWPkWG/PMqkLlN5o/8P7cg+qg9vZSW6F4ezf4wxOtCSTxkbGZmVnFnIzNzMwq1lQylnSgpLslLZR0Sp31H5Z0p6TbJP1S0vblh2pmZtadhkzGkvqAc4CDgN2AIyXtVlPsFmBGRPwtcDHwxbIDNTMz61bNHBnvBSyMiPsiYhXZdYeH5QtExK8jov8WPNcD08oN08zMrHs1k4ynAotyzxenZY0cB1xRb0X+9nnLly9vPkozM7MuVuoJXJKOAmYAX6q3Pn/7vC222KLMTZuZmY1azVxn/CAwPfd8Wlq2HkkHAKcBr4mIleWEZ2Zm1v2aOTK+EdhZ0o6SxpPdHu+yfAFJewD/BRwaEcvKD9PMzKx7DZmMI2INcCJwJdkPjV8UEfMlnSHp0FTsS8BGwA8lzZN0WYPqzMzMrEZTt8OMiMuBy2uWfSr3+ICS4zIzM+sZvgOXmZlZxZyMzczMKuZkbGZmVjEnYzMzs4o5GZvZeiRNl/Tr9OMv8yV9oOqYzLpdU2dTm1lPWQN8JCJuljQJuEnS1RFxZ9WBmXUrHxmb2XoiYklE3JweP0V2f4HB7kdvZiPkZGxmDUnaAdgDuKHOumd/+GXV2hUDXmtmzXMyNrO6JG0E/Aj4YET8uXZ9/odfxvdNbH+AZl3EydjMBpA0jiwRnx8Rl1Qdj1m3czI2s/VIEnAesCAivlp1PGa9wMnYzGq9Cngn8Lr0wy/zJB1cdVBm3cyXNpnZeiLit4AKvWbsGFZvudGA5Ss3HVdo2xsuHDA1PSiteKZQ+aLqtWkwo7m9vdRWKN7eddtu3KJIMj4yNjMzq5iTsZmZWcWcjM3MzCrmZGxmZlYxJ2MzM7OKORmbmZlVzMnYzMysYk7GZmZmFXMyNjMzq5iTsZmZWcWcjM3MzCrme1Ob2YhpzTrGLX1qwPKnpk8pVM9GK1cVKh9/HrjNQfX1FSpe7O7LHdbeXmortLy9T76kWHuL8pGxmZlZxZyMzczMKuZkbGYDSDpQ0t2SFko6pep4zLqdk7GZrUdSH3AOcBCwG3CkpN2qjcqsuzkZm1mtvYCFEXFfRKwCLgQOqzgms67ms6nNrNZUYFHu+WJg79pCko4Hjk9P/3LlXWfePaCmuxpuYwrwyEiCbIulBcvXb28vtRXcXgD0g4/WK7t9o0qaSsaSDgTOBvqAb0XEmTXrJwDfBV4GPAq8LSLub6ZuMxudImImMHM4r5U0NyJmlBxSR+qltoLbO1xDDlM3OX90HPB4ROwEfA34wkgDM7PKPAhMzz2flpaZWYs0M2fczPzRYcB30uOLgf0lqbwwzayNbgR2lrSjpPHAEcBlFcdk1tWaGaZuZv7o2TIRsUbSk8Dm1Iyj18wxrZR0x3CCbqNOn/vo9Pig82Ps9PgAdm3nxlIfPhG4kmxqalZEzC95M8Ma3h6leqmt4PYOS1tP4MrPMY2GeYVOj7HT44POj7HT44MsxnZvMyIuBy5vYf0984XdS20Ft3e4mhmmbmb+6NkyksYCm5CdyGVmZmZDaCYZNzN/dBlwdHr8FuBXERHlhWlmZta9hkzGEbEG6J8/WgBcFBHzJZ0h6dBU7Dxgc0kLgQ8Dzdw+bzQMZXR6jJ0eH3R+jJ0eH4yOGNcj6QOS7pA0X9IHc8tPknRXWv7FtGwHSX+VNC/9O7eywIehXlsl/SDXnvslzcuVPzXdZvRuSW+sKu7hKtLe0f7eQsP2vlTS9alNcyXtlZZL0tfT+3ubpD2b3lBE+J//+Z//lfYP2B24A5hIdl7KL4CdgNemxxNSuS3T/zsAd1Qdd5ltrSnzFeBT6fFuwK3ABGBH4F6gr+p2tLC9o/a9Hay9wFXAQanMwcA1ucdXAAL2AW5odlu+HaaZle1FZF9CKyIbWfsN8A/ACcCZEbESICKWVRhjWRq1FciOlIC3AnPSosOACyNiZUT8EVhIdvnoaFG0vaNdo/YGsHEqswnwUHp8GPDdyFwPTJa0TTMbcjI2s7LdAewraXNJE8mOFqYDu6TlN0j6jaSX516zo6Rb0vJ9qwh6mBq1td++wMMRcU96Xu9S0altibQcRdsLo/e9hcbt/SDwJUmLgC8Dp6byw35/W56Mh/opNkkT0nzDwtRJd2h1TAXj+7CkO9P4/y8lNby3aFUx5sr9o6SQ1PZLdZqJUdJb099yvqQLOik+SdtJ+nX60rhN0sFtjm+WpGVqcO39iOai2iwiFpDdhe8q4OfAPGAt2TDfZmTDdx8DLkpHUkuA7SJiD7JzTi6QtHGdqjvOIG3tdyTdc5Q4nPaO2vcWBm3vCcCHImI68CGy86ZGvLFWjrf3kc2JPB8YTzZXsltNmX8Fzk2PjwB+0Mb5gGbiey0wMT0+oZ3xNRtjKjcJuBa4HpjRaTECOwO3AJum51t2WHwzgRPS492A+9v8N3w1sCcN5tcYwVxU1f+Az6d+/nPgtbnl9wJb1Cl/Tbs/w2W3NT0eCzwMTMutPxU4Nff8SuAVVcfdqvZ203ubby/wJKC0TMCf0+P/Ao7Mlb8b2KaZult9ZNzpt9IcMr6I+HVErEhPrye7zrqdmv05u8+Q7cE9087gkmZifA9wTkQ8Dm2fL2wmvkZzQG0REdcCjw1SZNhzUVWQtGX6fzuyObYLgEvJdm6RtAvZjtEjkrZQdg98JD2fbMftvgrCHpYGbQU4ALgrIhbnil8GHJFGBHcka+vv2xnvSBVp72h/b6Fhex8CXpOKvA7oH5a/DHhXGsnaB3gyIpY0s51W34GrtFtpVhhf3nFkRyftNGSMachyekT8TNLH2hlc0szfcRcASb8jO1L9dET8vD3hNRXfp4GrJJ0EbEj2xdJJGs1FNdXRK/AjSZsDq4H3RcQTkmYBs9JQ/Crg6IgISa8GzpC0GlgHvDciBtsx6TQD2pqWH0HNEHVkl4VeBNwJrEnl88O8o0HT7SUb8RnN7y3U/yy/Bzhb2U2unuG52zxfTjaKtRBYARzb7Eb8e8ZNknQUMIPn9oY6gqQxwFeBYyoOZShjyfaK9yMbXbhW0t/kOnLVjgRmR8RXJL0C+J6k3SNiXdWBjUYRMeBEnTQqcVSd5T8CftSOuFqhXlvT8mMaLP8c8LlWxtRKRdo72t9baPhZ/i3ZTwbXLg/gfcPZTquHqTv9VppN/VScpAOA04BDI12W0UZDxTiJ7Fq4ayTdTzafeFmbT+Jq5u+4GLgsIlZHdknHH8iSc6fEdxxwEUBEXAdsQPYjEp3CP2to1sVanYw7/VaaQ8YnaQ+ySflD2zzP2VSMEfFkREyJiB0iYgeyee1DI6KdPy7QzPt8KdlRMZKmkA1bt2vuqJn4HgD2T/G9iCwZL29TfM0Y9lyUmXW+lg5TR4OfYpN0BjA3Ii4jOyX8e8pupfkY2RdlWzQZ35eAjYAfpvPKHoiIQxtWWk2MlWoyxiuBN0i6k+zSgI9FRFtGQJqM7yPAf0v6ENnJXMe0cacQSXPIdlamSFoMnA6MS/Gfywjmosys86mN3zdmZmZWh+/AZWZmVjEnYzMzs4o5GZuZmVXMydjMzKxiTsZmZmYVczI2MzOrmJOxmZlZxf4/FTY/7XIaIsEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fix, ax = subplots(2,2, figsize=(8,8))\n", "ax[0,0].imshow(dummy_image)\n", "ax[0,1].imshow(dummy_blurred)\n", "ax[1,1].imshow(csr.dot(dummy_blurred.ravel()).reshape(mask.shape))\n", "ax[0,0].set_xlim(964,981)\n", "ax[0,0].set_ylim(0,16)\n", "ax[0,0].set_title(\"Dummy image\")\n", "ax[0,1].set_xlim(964,981)\n", "ax[0,1].set_ylim(0,16)\n", "ax[0,1].set_title(\"Convolved image (i.e. blurred)\")\n", "ax[1,1].set_xlim(964,981)\n", "ax[1,1].set_ylim(0,16)\n", "ax[1,1].set_title(\"Retro-projected of the blurred\")\n", "ax[1,0].set_title(\"Corrected image\")\n", "pass" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def iterMLEM_scipy(F, M, R):\n", " \"Implement one step of MLEM\"\n", " #res = F * (R.T.dot(M))/R.dot(F)# / M.sum(axis=-1)\n", " norm = 1/R.T.dot(numpy.ones_like(F)) \n", " cor = R.T.dot(M/R.dot(F))\n", " res = norm * F * cor\n", " res[numpy.isnan(res)] = 1.0\n", " return res\n", "\n", "def deconv_MLEM(csr, data, thres=0.2, maxiter=1000):\n", " R = csr.T\n", " msk = data<0\n", " img = data.astype(\"float32\")\n", " img[msk] = 0.0 # set masked values to 0, negative values could induce errors\n", " M = img.ravel()\n", " #F0 = numpy.random.random(data.size)#M#\n", " F0 = R.T.dot(M)\n", " F1 = iterMLEM_scipy(F0, M, R)\n", " delta = abs(F1-F0).max()\n", " for i in range(maxiter):\n", " if delta:4: RuntimeWarning: divide by zero encountered in divide\n", " norm = 1/R.T.dot(numpy.ones_like(F))\n", ":5: RuntimeWarning: invalid value encountered in divide\n", " cor = R.T.dot(M/R.dot(F))\n", ":6: RuntimeWarning: invalid value encountered in multiply\n", " res = norm * F * cor\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0 1.7867398\n", "100 0.0016492605\n", "200 0.00027656555\n", "262 9.918213e-05\n", "CPU times: user 4.39 s, sys: 10.6 ms, total: 4.4 s\n", "Wall time: 4.4 s\n" ] }, { "data": { "text/plain": [ "(0.0, 16.0)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time res = deconv_MLEM(csr, dummy_blurred, 1e-4)\n", "\n", "ax[1,0].imshow(res)\n", "ax[1,0].set_xlim(964,981)\n", "ax[1,0].set_ylim(0,16)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion of the raytracing part:\n", "\n", "We are able to simulate the path and the absorption of the photon in the thickness of the detector. \n", "Numba helped substentially to make the raytracing calculation much faster. \n", "The signal of each pixel is indeed spread on the neighboors, depending on the position of the PONI and this effect can be inverted using sparse-matrix pseudo-inversion. \n", "The MLEM can garanteee that the total signal is conserved and that no pixel gets negative value.\n", "\n", "We will now save this sparse matrix to file in order to be able to re-use it in next notebook. But before saving it, it makes sense to spend some time in generating a high quality sparse matrix in throwing thousand rays per pixel in a grid of 64x64." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 34min 16s, sys: 7.83 s, total: 34min 24s\n", "Wall time: 34min 16s\n" ] } ], "source": [ "%time pre_csr = thick.calc_csr(64)\n", "hq_csr = csr_matrix(pre_csr)\n", "from scipy.sparse import save_npz\n", "save_npz(\"csr.npz\",hq_csr)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total execution time: 2109.755 s\n" ] } ], "source": [ "print(f\"Total execution time: {time.perf_counter()-start_time:.3f} s\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.9.5" } }, "nbformat": 4, "nbformat_minor": 4 }