diff --git a/examples/component_example.ipynb b/examples/components.ipynb similarity index 100% rename from examples/component_example.ipynb rename to examples/components.ipynb diff --git a/examples/convolution.ipynb b/examples/convolution.ipynb new file mode 100644 index 0000000..4fa9ea6 --- /dev/null +++ b/examples/convolution.ipynb @@ -0,0 +1,144 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "f42e34d0", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from easydynamics.sample_model import DeltaFunction, Gaussian, Lorentzian, SampleModel, DampedHarmonicOscillator\n", + "from easydynamics.utils import convolution \n", + "from easydynamics.utils import _detailed_balance_factor as detailed_balance_factor\n", + "\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "600c0850", + "metadata": {}, + "outputs": [], + "source": [ + "# Standard example of convolution of a sample model with a resolution model\n", + "sample_model=SampleModel(name='sample_model')\n", + "gaussian=Gaussian(name='Gaussian',width=0.5,area=1)\n", + "dho = DampedHarmonicOscillator(name='DHO',center=1.0,width=0.3,area=2.0)\n", + "lorentzian = Lorentzian(name='Lorentzian',center=-1.0,width=0.2,area=1.0)\n", + "delta = DeltaFunction(name='Delta',center=0.4,area=0.5)\n", + "sample_model.add_component(gaussian)\n", + "sample_model.add_component(dho)\n", + "sample_model.add_component(lorentzian)\n", + "sample_model.add_component(delta)\n", + "\n", + "resolution_model = SampleModel(name='resolution_model')\n", + "resolution_gaussian=Gaussian(name='Resolution Gaussian',width=0.2,area=0.8)\n", + "resolution_lorentzian = Lorentzian(name='Resolution Lorentzian',width=0.2,area=0.2)\n", + "resolution_model.add_component(resolution_gaussian)\n", + "resolution_model.add_component(resolution_lorentzian)\n", + "\n", + "energy=np.linspace(-2, 2, 100)\n", + "\n", + "\n", + "y = convolution(sample_model=sample_model,\n", + " resolution_model=resolution_model,\n", + " energy=energy,\n", + " )\n", + "\n", + "plt.plot(energy, y, label='Convoluted Model')\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('Intensity (arb. units)')\n", + "plt.title('Convolution of Sample Model with Resolution Model')\n", + "\n", + "plt.plot(energy, sample_model.evaluate(energy), label='Sample Model', linestyle='--')\n", + "plt.plot(energy, resolution_model.evaluate(energy), label='Resolution Model', linestyle=':')\n", + "\n", + "\n", + "plt.legend()\n", + "# set the limit on the y axis\n", + "plt.ylim(0,6)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fede1a58", + "metadata": {}, + "outputs": [], + "source": [ + "# Use some of the extra settings for the numerical convolution\n", + "sample_model=SampleModel(name='sample_model')\n", + "gaussian=Gaussian(name='Gaussian',width=0.5,area=1)\n", + "dho = DampedHarmonicOscillator(name='DHO',center=1.0,width=0.3,area=2.0)\n", + "lorentzian = Lorentzian(name='Lorentzian',center=-1.0,width=0.2,area=1.0)\n", + "delta = DeltaFunction(name='Delta',center=0.4,area=0.5)\n", + "sample_model.add_component(gaussian)\n", + "sample_model.add_component(dho)\n", + "sample_model.add_component(lorentzian)\n", + "\n", + "resolution_model = SampleModel(name='resolution_model')\n", + "resolution_gaussian=Gaussian(name='Resolution Gaussian',width=0.2,area=0.8)\n", + "resolution_lorentzian = Lorentzian(name='Resolution Lorentzian',width=0.2,area=0.2)\n", + "resolution_model.add_component(resolution_gaussian)\n", + "resolution_model.add_component(resolution_lorentzian)\n", + "\n", + "energy=np.linspace(-2, 2, 100)\n", + "\n", + "\n", + "temperature = 15.0 # Temperature in Kelvin\n", + "offset = 0.3\n", + "upsample_factor = 5\n", + "extension_factor = 0.2\n", + "\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('Intensity (arb. units)')\n", + "\n", + "y = convolution(sample_model=sample_model,\n", + " resolution_model=resolution_model,\n", + " energy=energy,\n", + " offset=offset,\n", + " method=\"auto\",\n", + " upsample_factor=upsample_factor,\n", + " extension_factor=extension_factor,\n", + " temperature=temperature,\n", + " normalize_detailed_balance=True,\n", + " )\n", + "\n", + "plt.plot(energy, y, label='Convoluted Model')\n", + "\n", + "plt.plot(energy, sample_model.evaluate(energy-offset)*detailed_balance_factor(energy-offset, temperature), label='Sample Model with DB', linestyle='--')\n", + "\n", + "plt.plot(energy, resolution_model.evaluate(energy ), label='Resolution Model', linestyle=':')\n", + "plt.title('Convolution of Sample Model with Resolution Model with detailed balancing')\n", + "\n", + "plt.legend()\n", + "plt.ylim(0,2.5)\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "newdynamics", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/sample_model.ipynb b/examples/sample_model.ipynb new file mode 100644 index 0000000..732fa86 --- /dev/null +++ b/examples/sample_model.ipynb @@ -0,0 +1,144 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "64deaa41", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from easydynamics.sample_model import Gaussian\n", + "from easydynamics.sample_model import Lorentzian\n", + "from easydynamics.sample_model import DampedHarmonicOscillator\n", + "from easydynamics.sample_model import Polynomial\n", + "\n", + "from easydynamics.sample_model import SampleModel\n", + "\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "%matplotlib widget" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "784d9e82", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d473593b3aa14baf8f8c4dd432169d44", + "version_major": 2, + "version_minor": 0 + }, + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAwG1JREFUeJzsnQV4U1cbx/91dzcqUEopFIq7uztsDAbb2AYT5hvfjClzNjZgYwIzGDB0wGC4W9FiRQp1t9T9e95zm9JCWypJY+/vee6TmzRNTpJ7z/nfV/XKysrKwDAMwzAMw+gM+qoeAMMwDMMwDNO0sABkGIZhGIbRMVgAMgzDMAzD6BgsABmGYRiGYXQMFoAMwzAMwzA6BgtAhmEYhmEYHYMFIMMwDMMwjI7BApBhGIZhGEbHYAHIMAzDMAyjY7AAZBiGYRiG0TFYADIMwzAMw+gYLAAZhmEYhmF0DBaADMMwDMMwOgYLQIZhGIZhGB2DBSDDMAzDMIyOwQKQYRiGYRhGx2AByDAMwzAMo2OwAGQYhmEYhtExWAAyDMMwDMPoGCwAGYZhGIZhdAwWgAzDMAzDMDoGC0CGYRiGYRgdgwUgwzAMwzCMjsECkGEYhmEYRsdgAcgwDMMwDKNjsABkGIZhGIbRMVgAMgzDMAzD6BgsABmGYRiGYXQMFoAMwzAMwzA6BgtAhmEYhmEYHYMFIMMwDMMwjI7BApBhGIZhGEbHYAHIMAzDMAyjY7AAZBiGYRiG0TFYADIMwzAMw+gYLAAZhmEYhmF0DBaADMMwDMMwOgYLQIZhGIZhGB2DBSDDMAzDMIyOwQKQYRiGYRhGx2AByDAMwzAMo2OwAGQYhmEYhtExWAAyDMMwDMPoGCwAGYZhGIZhdAwWgAzDMAzDMDoGC0CGYRiGYRgdgwUgwzAMwzCMjsECkGEYhmEYRsdgAcgwDMMwDKNjsABkGIZhGIbRMVgAMgzDMAzD6BgsABmGYRiGYXQMFoAMwzAMwzA6BgtAhmEYhmEYHYMFIMMwDMMwjI7BApBhGIZhGEbHMIQWsHz5crHduXNH3A8KCsI777yD4cOHV/v8VatWYfbs2VUeMzExQX5+fr3et7S0FHFxcbCysoKenl4jPgHDMAzDME1FWVkZsrKy4O7uDn193bSFaYUA9PT0xCeffAJ/f3/xo/76668YO3Yszp07J8RgdVhbWyM8PLzifkMEHIk/Ly+vRo2dYRiGYRjVEB0dLTSELqIVAnD06NFV7n/00UfCInjixIkaBSAJPldX10a9L1n+5AcQCUqGYRiGYdQfmUwmDDjydVwX0QoBWJmSkhKsX78eOTk56N69e43Py87Ohre3t3DjdujQAR9//HGNYrEm5FZDEn8sABmGYRhGs9DT4fAtrRGAYWFhQvBRHJ+lpSU2bdqE1q1bV/vcgIAA/PLLLwgODkZmZia++OIL9OjRA5cvX67VFFxQUCC2ylcQDMMwDMMwmoZeGQXNaQGFhYWIiooSgu7vv//GTz/9hIMHD9YoAitTVFSEwMBAPPTQQ/jggw9qfN7ChQvx3nvv3fc4vSdbABmGYRhGM5DJZLCxsdHp9VtrBOC9DBo0CM2bN8cPP/xQp+dPnjwZhoaGWLNmTb0sgBRDoMsHEMMwDMNoGiwAtcgFfC8U21dZrD0obpBcyCNGjKj1eVQqhjaGYRhGeZBdori4WMzNDNMQDAwMhFFHl2P8dEIALliwQNT8a9asmajrs3r1ahw4cAC7du0Sf585cyY8PDywaNEicf/9999Ht27d0KJFC2RkZODzzz9HZGQknnjiCRV/EoZhGN2Gwnni4+ORm5ur6qEwGo65uTnc3NxgbGys6qGoJVohAJOSkoTIo0mDTLqU3EHib/DgweLvFBtYudBjeno65syZg4SEBNjZ2aFjx444duxYneIFGYZhGOV5bm7fvi2sN1SglxZutuAwDbEg04VEcnKyOJ6oRrCuFnvWyRjApoBjCBiGYRQHVXGgBZtKdJH1hmEaA1mRybvn6+sLU1PTKn+T8frNvYAZhmEY9YKtNYwi4OOodvjbYRiGYRiG0TFYADIMwzCMFkNxlJs3b4a6069fP7zwwgt1fv6qVatga2ur1DFpMywAGYZhGKYRULLB3LlzRSUKKhVGfeaHDh2Ko0ePQhu4c+eOEJGUnBMbG1vlb5R8KS+3Qs9jNAcWgAzDMAzTCCZOnIhz587h119/xfXr17F161ZhzUpNTYU2QeXUfvvttyqP0WemxxnNgwUgwyiJ8IQsfLDtCuIz81Q9FIZhlATVkj18+DA+/fRT9O/fX2Qwd+nSRdSnHTNmTMXzvvrqK7Rt2xYWFhaig9S8efOQnZ19nztz27Ztol89ZUFPmjRJZLKSyPLx8RFly55//vkqBbLpcWphSq1M6bVJjC1durTWMUdHR2PKlCni/ezt7TF27Ng6We8effRRrFy5sspjdJ8evxdqxUrfA1lEqRbfG2+8IYp7y8nJyRHl2ywtLcXfv/zyy/teg5o5vPLKK+Iz0Wfr2rWrqPHLKAYWgAyjBKi60msbLuLnI7fxyE8nkZ5TqOohMYxGnke5hcUq2epaIY0EDG0UY1db9ynKSF2yZAkuX74sBN2+ffvw2muvVXkOiT16zl9//YWdO3cKsTN+/Hjs2LFDbL///rtob0r97itDzQzatWsnrJAktObPn4/du3dXO46ioiLhnrayshLCldzUNP5hw4aJ2nm1QYKW6ugeOXJE3Kdbuj969OgqzyM3MXXW6ty5My5cuIDly5fj559/xocffljxnFdffVWIxC1btuC///4Tn/Xs2bNVXufZZ5/F8ePHxfdx8eJF0bKVxnnjxo1ax8noUCFohlE3Tt1Ow4XoDLF/KzkHs1edxuo5XWFuzKccw9SVvKIStH5H6ujU1Fx5f2idzleKfyPrHTUX+P7779GhQwf07dsX06ZNE00J5FRObiCrHYmhp59+GsuWLasizkgsUR97giyAJPoSExOFSKNmBWRl3L9/P6ZOnVrxfz179hTCj2jZsqUQdYsXL65ohlCZtWvXioLbP/30U0WRbbLikTWQRNiQIUNq/KxGRkZ45JFH8Msvv6BXr17ilu7T45Whz0RWzu+++068R6tWrRAXF4fXX38d77zzjhC6JAj/+OMPDBw4UPwPiWJPT8+K16AGDjQuuqWi4ARZA0kY0+Mff/zxA38bpnbYAsgwSmDFoQhx2z/ACbbmRjgfnYF5f55FUUmpqofGMIwSYgBJ4FDsH1moSEiRECRhKGfPnj1C7JA7k6xvM2bMEDGClVvekdtXLv4IFxcXIRZJ/FV+jLpfVaZ79+733b969Wq1YyWL3M2bN8UY5NZLcgNTEe5bt2498LM+9thjWL9+veikRbd0/17ovWkMlbu4kEgll3dMTIx4H7I2kktXDo2BXN9ywsLChKubBK18nLSR1bAu42QeDJsjGEbB3EjMwt5rSaC5753RQUjLKcT0n07gQHgyXvv7Ir6c3A76+tzeimEehJmRgbDEqeq96wN1miCLG21vv/226C3/7rvvYtasWSK+btSoUSJT+KOPPhJih9ynjz/+uBBC8q4n91rSSEBV9xhZ8BoKiTBqf/rnn3/e9zcnJ6cH/j/FMZJFj2IOAwMD0aZNG5w/f77B46ltnJR1fObMGXFbmcqCmGk4LAAZRknWv2FBrvB1tBDb8ukd8cRvodh0LhbO1iZYMDxQ1cNkGLWHxI6mhk2Qu1Zee49EDIk2SnSQd6dYt26dwt7rxIkT990ncVYdZJkkN7Czs3ODW6CR1Y+SWMhdXR303hs2bBBxlHIrILmlyepIbl4SwCRsT548KUrnEBRLSBnU5D4nQkJChAWQrJ29e/du0DiZ2mEXMMMokITMfGw+L9XJerKPX8Xj/Vs54/NJwRUCkayCDMNoPuTGHTBggIhno0QF6mVMrtHPPvtMZNcSLVq0EPF93377LSIiIkRcH8ULKgoSV/R+JKAoA5jenxJBqmP69OlwdHQUY6MkEBovuawpu5jcs3WB4h2p9iFZOauDxCFlGj/33HO4du2aSPQga+hLL70kBDBZ8Mj6SYkglAxz6dIlYSmt3LqNXL80VsoU3rhxoxjnqVOnsGjRImzfvr2B3xRTGRaADKNAVh67jaKSMnTxtUdIM7sqf5vQwROtXK1AyYXHbqWobIwMwygOEjMUy0ZJF3369BEuUXIBk0iiJAiCMnSpDAyViqG/k/uVhIyiePnllxEaGiqsZpRcQu9Fmb7VQe7mQ4cOCcvbhAkThLWOxBjFANbVIkiJLyQi6bY6KM6RspZJsNFnp2QXeo+33nqrSuYyWfYog3jQoEEiqYRc05WhZA8SgPT5KD5w3LhxOH36dIXVkGkcemV1zXVn7kMmk8HGxgaZmZkNNqUz2kNWfhF6LNqHrIJi/PxoJwwMdLnvOR9uu4KfjtzGtM5e+GTi3QxBhmEgRAhZenx9fUVMHfNgKEmEMozr00JNV6jteJLx+s0WQIZRFGtORQnx18LZEv0DnKt9Tk9/R3F7+EZKneuMMQzDMIyiYQHIMAqgsLgUvxy5UxH7V1OWb1dfexgZ6CE2Iw+RqXfLPzAMwzBMU6KZ6VUMo2acjUpHgiwfDhbGGNteKlpaHZTR2KGZHU7eTsPhmynwcbRo0nEyDKNd1KWFG8NUB1sAGUYBhMVkitvOPvYwMay9fljvcjfwkRvJTTI2hmEYhrkXFoAMowAuxkoCsK2nzQOf28tfKrZ67FYqSko5DpBhGIZpelgAMowCuBgj9f0NroMAbOthA2tTQ2TlF1f8H8MwDMM0JSwAGaaRZOYWVSR0BHvYPvD5Bvp66NFc7gbmeoAMwzBM08MCkGEaSVi5+9fbwRw25lX7dtZEL3kc4E0WgAzDMEzTwwKQYRrJhXI3Lrl260qvFo4V2cM5BcVKGxvDMAzDVAcLQIZRUAZwO88Hu3/lkLXQ085MtI07dTtNiaNjGIaRWLVqFWxt6z5PMdoNC0CGUZALuC4ZwHL09PTuloNhNzDDaAUJCQmYP38+WrRoIVqPubi4oGfPnli+fDlyc1Vf+H3q1Km4fv26qofBqAlcCJphGkFKdoHo6qGnBwS516+fZM8WjlhzKpoTQRhGC4iIiBBijyxsH3/8Mdq2bQsTExOEhYVhxYoV8PDwwJgxY1Q6RjMzM7ExDMEWQIZRgPvXz9ECVqZ1SwCRQ5nAJBzDE7OQJMtX0ggZhmkK5s2bB0NDQ4SGhmLKlCkIDAyEn58fxo4di+3bt2P06NHieV999ZUQhxYWFvDy8hL/l52dXfE6CxcuRPv27au89tdffw0fH5+K+wcOHECXLl3Ea5DgJOEZGRkp/nbhwgX0798fVlZWsLa2RseOHcWYqnMB37p1S4yPLJWWlpbo3Lkz9uzZU+W96X1J0D722GPiNZs1ayYELaP5sABkmEZwsQHxf3LsLYwrrIZHb7EVkGHuo6wMKMxRzUbvXUdSU1Px33//4ZlnnhGirKawD0JfXx9LlizB5cuX8euvv2Lfvn147bXX6vxexcXFGDduHPr27YuLFy/i+PHjePLJJytef/r06fD09MTp06dx5swZvPHGGzAyqv7ilITniBEjsHfvXpw7dw7Dhg0TQjUqKqrK87788kt06tRJPIcE69y5cxEeHl7nMTPqCbuAGaYRyAs51yf+71438KVYGU5GpGF8iKeCR8cwGk5RLvBxzb21lcr/4gDjuvXqvnnzJsrKyhAQEFDlcUdHR+TnS9Z9EoeffvopXnjhhSrWtQ8//BBPP/00li1bVqf3kslkyMzMxKhRo9C8eXPxGFkb5ZB4e/XVV9GqVStx39/fv8bXateundjkfPDBB9i0aRO2bt2KZ599tuJxEokk/IjXX38dixcvxv79++/7vIxmwRZAhmkgNOHLW8DVpQNIdbQvtxxejpMpdGwMw6ieU6dO4fz58wgKCkJBQYF4jFysAwcOFDGB5FKdMWOGsCDWNUnE3t4es2bNwtChQ4W17ptvvkF8fHzF31966SU88cQTGDRoED755BPh5q0JsgC+8sorQkCSa5jcwFevXr3PAhgcHFyxT5ZGV1dXJCUlNeAbYdQJtgAyTANJlBUgOatAdPZo7dYwAdimvHZgeEIWCotLYWzI12QMU4GRuWSJU9V71xHK+iVhdK9blGIACXnixZ07d4TljlyoH330kRBzR44cweOPP47CwkKYm5sLFzFdXFamqKioyv2VK1fi+eefx86dO7F27Vq89dZb2L17N7p16yZiCB9++GERd/jvv//i3XffxV9//YXx48ffN24Sf/R/X3zxhfgMNM5JkyaJsVT5Ku5xIdNnLS0trfP3w6gnvNowTCMLQPs7W8LM2KBBr0G1AKkvcGFJKW4kZSl4hAyj4VBcG7lhVbGVx9TVBQcHBwwePBjfffcdcnJyanwexeSRcKKYOhJrLVu2RFxcVYHr5OQkyslUFoFkRbyXkJAQLFiwAMeOHUObNm2wevXqir/R67744osiLnHChAlCMFbH0aNHhTWRxCElppBlj0QqoxuwAGSYRmYAN9T9K7+SllsBL8eyG5hhNBWK4aMEDUqWIKscuVLJIvjHH3/g2rVrMDAwEFY2suZ9++23omzM77//ju+//77K6/Tr1w/Jycn47LPPhPt26dKlwpIn5/bt20L4UfIHZf6SyLtx44Zw4+bl5YnYPcoSpr+RwKNkkMoxgpWh+MCNGzcKgUnZw2Q5ZMue7sACkGEayN34v8ZV1pcLwEtx0usxDKN5UEIGZclS7B0JNEquIDFIYo9crZRgQY9RGRhKBiGr3Z9//olFixZVeR0SayQmSfjR8ymOkP5fDrmJSVBOnDhRWPooA5gSTJ566ikhMimecObMmeJvVI5m+PDheO+996odM43Fzs4OPXr0EPGEFFfYoUMHpX9XjHqgV3ZvsAFTZygby8bGRmRkUb0lRneg0ybkg93IyC3C1md7NkoEbjkfi/l/nUeHZrbYOK+nQsfJMJoEZcyShcvX11d00mAYZR1PMl6/2QLIMA0hJj1PiD8jAz0EuFo16rWC3CUL4JV4GUpK+XqMYRiGUT4sABmmEQWgA92sYWLYsAQQOb6OFjA3NkB+USkiku92BGAYhmEYZcECkGEaUwC6PH6vMUhlZCQXBMcBMgzDME0BC0CGaQByoaYIAVglEYQzgRmGYZgmQCsE4PLly0WlcgrkpK179+5V0uarY/369aJVDgWGUv2jHTt2NNl4Gc3neqLkqm1VbrlrLPKewJfKM4sZhmEYRplohQCkxtfU8oaKbIaGhmLAgAEYO3asaLZdHVQ486GHHhLV1yltnxpr03bp0qUmHzujeWTmFokOIERzp7r1Cq2rBfBKnAylnAjCMAzDKBmtEIBUv4iaVVNRS6p9RC12qKfhiRMnqn0+9U4cNmyYaJhNNZeoPhPVPqIq7gzzIG4mSx07XK1NYWVatUVSQ2nhbCnawGUVFCMqrW49QRmGYRhGpwVgZUpKSkTfQ2rHQ67g6qAK6lSsszJUAJMerw1q5k21gypvjO5xM0ly//q7WCrsNY0M9BFYXk6GE0EYhmEYZaM1AjAsLExY/UxMTPD0009j06ZNaN26dbXPpT6LLi4uVR6j+/R4bVDFdiocKd+8vLwU+hkYzRKAzZ0UJwCJIE4EYRiGYZoIrRGAAQEBop/hyZMnMXfuXDz66KO4cuWKQt+D2vtQ1XD5Fh0drdDXZzSDG0qwABJtygtCX2YLIMMwDKNktEYAGhsbi0bbHTt2FJY66qFIsX7V4erqisTExCqP0X16vDbIuijPNJZvjO5aAFso2ALYxuNuJjB3aGQYzWLWrFnQ09MTm5GRkfAqDR48GL/88gtKS0srnufj44Ovv/76vv9fuHAh2rdvX+WxtLQ0vPDCC/D29hZrnLu7Ox577DFERUU1yWditButEYD3QiccxexVB8UG7t27t8pju3fvrjFmkGHk5BYWizZw8sQNRdLSxQqG+npIzy1CXGa+Ql+bYRjlQ8mF8fHxuHPnjihF1r9/f8yfPx+jRo1CcXFxvV6LxF+3bt2wZ88efP/997h586aIb6fbzp07IyIiQmmfg9ENDKEFkGt2+PDhaNasGbKysrB69WocOHAAu3btEn+fOXMmPDw8hGWQoBOyb9+++PLLLzFy5EhxUlH5mBUrVqj4kzDqTkRyjri1tzCGg6WJQl/b1MgA/i5WuBovE1ZAD1szhb4+wzDKhbxEck8SrTlUXYJE3MCBA7Fq1So88cQTdX6tN998E3FxcULwyV+T1jha16jixTPPPPPAercMo/UCMCkpSYg8uvKi5AwqCk0nCZnfCTKX6+vfNXb26NFDiMS33noL//vf/8TJtHnzZrRp00aFn4LRZfevnDbu1kIAXo7NxNCg2kMSGEbboVCIvGLJ4t7UmBmaCXduY6G6tBSStHHjxjoLQPJgkWFi+vTp94UmmZmZYd68eWL9Iiuhvb19o8fI6CZaIQB//vnnWv9O1sB7mTx5stgYpj7cSJJqALZQcAJI5YLQ68/E4FIcZwIzDIm/rqu7quS9Tz58EuZG5gp5Leo6dfHixYr7r7/+uhBwlSksLKyoXJGcnIyMjAxRp7Y66HESx2Qd7NKli0LGyOgeWiEAGUZrLICVEkEYhtEOSKxVtiZSEwJKGqnMkiVLcOjQofv+j2GUBQtAhmlACRhFJ4DICXSzhr4ekJRVgKSsfDhbmSrlfRhGEyA3LFniVPXeiuLq1avw9fWtuO/o6CiqVlSmsivXyckJtra24v9qej0SlPe+BsPUBxaADFNHCotLEZmaq5QagHLMjQ3h42ghkk2uxmexAGR0GhI5inLDqop9+/aJRgUvvvhinf+HYtanTJmCP//8E++//36VOMC8vDwsW7ZMdK/i+D+mMWhtGRiGUTSRqTkoKS2DpYmh6AOsLMgKSFAyCMMwmgOVHqOOUrGxsTh79iw+/vhjjB07VpSBoUTF+kD/S8KPkhkp25caD5CLmIRfUVERli5dqrTPwegGLAAZpp7u3+bOlgrJDqyJ1uUC8BoLQIbRKHbu3Ak3NzdR7JlqAu7fv1/E9m3ZsgUGBgb1ei0HBwecOHFC1BJ86qmn0Lx5c2EVpNvTp0/Dz89PaZ+D0Q30yjjKtMHIZDJRdobawnFXEO1nyd4b+Gr3dUzs4Ikvp7RT2vvsvZqIx38NRYCLFXa92Edp78Mw6kZ+fj5u374t4uVMTTn8gVHe8STj9ZstgAyjLgkg97qAbyVno6C4RKnvxTAMw+gmLAAZpp4lYPyVLADdbExhY2aE4tIy3EiU3pNhGIZhFAkLQIapA5T8EZHcNBZAii8MdLMS+5wIwjAMwygDFoAMUwdi0nNRUFwKY0N9eNkrvyzF3UxgqfMIwzAMwygSFoAMUw/3r5+jBQyoUnMTCcBrCWwBZBiGYRQPC0CGUaMEEDmBrndrAXKiPsMwDKNoWAAyTL0SQKTYPGVDnUbI0pieW4REWUGTvCfDMAyjO7AAZBg1tACaGhkIdzPBiSAMwzCMomEByDAPgFywt+QWQCX1AK4tDvAKC0CGYRhGwbAAZJgHQC7Y7IJi4ZL1cZCsck0B9wRmGKapuXPnjihFdf78eVUPhVEyLAAZ5gHcSJJKsXjbm4syME2FvBbgtQQuBcMw6s6sWbMwbtw4qBsk5jZv3lzn53t5eSE+Ph5t2rRR6rgY1cMCkGHqmADSVPF/91oAqQB1fhG3hGMYXaGwsFBl721gYABXV1cYGhqqbAxM08ACkGHqmADSlPF/hLOVCewtjFFaBlxPZCsgw2gqBw8eRJcuXWBiYgI3Nze88cYbKC4urvh7v3798Oyzz+KFF16Ao6Mjhg4dKh6/dOkShg8fDktLS7i4uGDGjBlISUmp8n/PP/88XnvtNdjb2wvhtnDhwoq/+/j4iNvx48cLS6D8Pt3S/Xu36lzAJSUlePzxx+Hr6wszMzMEBATgm2++qdb6+cUXX4jP5+DggGeeeQZFRUVK/V6ZxsESn2EewM1E1VgA5S3hjt5MFXGAwZ62Tfr+DKMOCVhleXkqeW89M7MKUdQYYmNjMWLECCGSfvvtN1y7dg1z5syBqalpFbH266+/Yu7cuTh69Ki4n5GRgQEDBuCJJ57A4sWLkZeXh9dffx1TpkzBvn37qvzfSy+9hJMnT+L48ePifXr27InBgwfj9OnTcHZ2xsqVKzFs2DBh3SPocRJ2BN1OmjQJRkZG1Y6/tLQUnp6eWL9+vRB2x44dw5NPPimEHo1Fzv79+8VjdHvz5k1MnToV7du3F5+VUU9YADLMA7iZ3LQ1AO8tCC0JQLYAMroHib/wDh1V8t4BZ89Az7zxbR+XLVsm4uq+++47IShbtWqFuLg4Iebeeecd6OtLjjh/f3989tlnFf/34YcfIiQkBB9//HHFY7/88ot4revXr6Nly5biseDgYLz77rsVr0Hvs3fvXiEAnZycxOO2trbCOihH/jgxf/58EfNHorA6SBi+9957FffJEkhCc926dVUEoJ2dnXhvEpn0GUeOHCnGwQJQfWEByDC1kJpdgLScQpAhoLlT01oACc4EZhjN5urVq+jevXsVayJZ6LKzsxETE4NmzZqJxzp2rCp0L1y4IKxp5P69l1u3blURgJUhK1xSUlKdxrZixQr8/PPPwqpXWRTey9KlS4X4jIqKEpZIilEk615lgoKCKiyM8nGEhYXVaRyMamAByDB1SADxsDWDmfHdya2paFWeCSxvCacIlxTDaArkhiVLnKreuymxsKhaYooE4ujRo/Hpp5/e91wSV3Ludd3SHEFu2wdB4vK5557DmjVr7hORlfnrr7/wyiuv4MsvvxRC1srKCp9//rlwOVemoeNgVAcLQIapSwJIE8f/yaG4Q0N9PcjyixGXmS+EKMPoCiI5QQFuWFUSGBiIDRs2VLmAozg/ElIUW1cTHTp0EP9HCRuNycglYSaP95NDMXoU9/e///0PEyZMqPX/aaw9evTAvHnzqlggGc2Hs4AZpi49gF2aPv6PMDE0qEg+uRrHbmCGUWcyMzNF9mzljRImoqOjhbWNEkC2bNkiYvYocUMe/1cdlEWblpaGhx56SMTnkejatWsXZs+efZ+gqw0SkBSLl5CQgPT0dOHCJcsixRfS2Ohx+VYdFFcYGhoq3ptiD99+++0a4wUZzYIFIMPUpQagCuL/5HAcIMNoBgcOHBDCqvL2wQcfYMeOHTh16hTatWuHp59+WpRVeeutt2p9LXd3d2F9I7E3ZMgQtG3bVpSJoYSO2oTjvZDrdvfu3SJ5hMaTmJgohCiJQnoPcifLt+p46qmnhJWQsnq7du2K1NTUKtZARnPRKyO7NNMgZDIZbGxsxFWftbW0SDPaRdeP94hWcBvn9UCHZnYqGcOKQ7fw8Y5rGNHWFcumqyYjkmGagvz8fNy+fVtkmlKZFIZR1vEk4/WbLYAMUxOy/CIh/lRRA7Ayrd1sxO1ldgEzDMMwCoIFIMM8wP3ram0Ka9Pqi6Q2BUHu0tVpZGquEKUMwzAM01hYADKMmnUAuRc7C+OK7N8rbAVkGIZhFAALQIapgRtJWWohACtbAdkNzDAMwygCFoAM88ASMOogAMvjAGMzVT0UhmEYRgtgAcgwDygCrcoSMHLYAsgwDMMoEhaADFMNuYXFiEnPU2kR6Mq08ZAsgDeTs5FfVPcisAzDMAxTHSwAGaYaIpJzxK2DhTHsLYxVPRy4WJuIsZSUluFaghSbyDAMwzANhQUgw9SSANJcDRJACOoh2rrCDcxxgAzDMEzjYAHIMLUlgKiJAKzsBuY4QIbRLlatWiVavGkCCxcuRPv27et9Abt582aljYlpGCwAGaYabiSqnwCsSAThTGCGUTtmzZolhA5txsbGaNGiBd5//30UFxdDm3jllVdEH2FG8zFU9QAYRp0tgC2cVZ8Acm8pmKsJWSgqKYWRAV+/MYw6MWzYMKxcuRIFBQXYsWMHnnnmGRgZGWHBggXQFiwtLcXGaD5asYIsWrQInTt3hpWVFZydnTFu3DiEh4c/0OQuv1qTb9x8nCEKiksQmZarNjUA5Xjbm8PSxBCFxaW4lSwJVIZh1AcTExO4urrC29sbc+fOxaBBg7B161akp6dj5syZsLOzg7m5OYYPH44bN25U+xp37tyBvr4+QkNDqzz+9ddfi9ctLS3FgQMHxJpFlrhOnTqJ1+zRo8d9697y5cvRvHlzYZEMCAjA77//XuXv9Bo//PADRo0aJV4jMDAQx48fx82bN9GvXz9YWFiI171161aNLuDTp09j8ODBcHR0hI2NDfr27YuzZ88q6BtllIlWCMCDBw+KK60TJ05g9+7dKCoqwpAhQ5CTI2Vy1oS1tTXi4+MrtsjIyCYbM6O+3EnJFdm2VqaGcLYygbqgr6+H1m5yNzDHATLaT1lZGYoKSlSy0Xs3FjMzMxQWFgr3MAk6EoMksOi1R4wYIdaqe/Hx8RHCkSyJlaH79DokDuW8+eab+PLLL8VrGxoa4rHHHqv426ZNmzB//ny8/PLLuHTpEp566inMnj0b+/fvr/K6H3zwgRCn58+fR6tWrfDwww+L55LVkl6Xxvrss8/W+BmzsrLw6KOP4siRI2IN9vf3F5+NHmfUG61wAe/cufM+6x5ZAs+cOYM+ffrU+H909UNXawxTUws4OkbUCcoEPnUnTSSCTOyo6tEwjHIpLizFivkHVfLeT37TF0YmBg36XxJNZJ3btWuXsPZRAsTRo0eFNY34888/4eXlJR6fPHnyff//xBNP4Omnn8ZXX30lrIpkUQsLC8OWLVuqPO+jjz4SFjfijTfewMiRI5Gfny+8WV988YUQjPPmzRN/f+mll4RAo8f79+9f8RokCqdMmSL2X3/9dXTv3h1vv/02hg4dKh4jEUnPqYkBAwZUub9ixQqR0EKGGbIsMuqLVlgA7yUzUwqSt7e3r/V52dnZwqROJ+LYsWNx+fLlWp9PcR0ymazKxmgf6pgBfG8m8CUuBcMwase2bdtEfBwJMBJ+U6dOFSKMrHNdu3ateJ6Dg4NwyV69erXa16EwJgMDA2HFkxs1SLSRdbAywcHBFftubm7iNikpSdzSa/fs2bPK8+n+ve9Z+TVcXFzEbdu2bas8RqKypvUuMTERc+bMEZY/cgGTZ43W1qioqAd+X4xq0QoLYGUoPuKFF14QB3qbNm1qfB6dfL/88os4+Ekw0lURXZ2RCPT09Kwx1vC9995T4ugZdWoB569GCSD3ZgJfjZOhtLRMuIUZRlsxNNYXljhVvXd9IZFGcXcUc+fu7i6EH7l96wv9P7llye07YcIErF69Gt988819z6MEEzlybwWtgfWhuteoz+uS+zc1NVWMjwwqZLEkKyK5vhn1RusEIMUCUrwDxSPUBh2gtMkh8UcBsBQQSzER1UExEWRGl0NXRGQ9ZLSLm+UlYMgFrG7QmIwN9ZFVUIzo9Fx4O1ioekgMozRIfDTUDasKKGmCyr9UhtYVKgVz8uTJChcwCSZK2GjdunWNr0VuYDJiLFu2TPw/CcH6QO9LbmcSaHLofm3v2RDoNWmMFPdHREdHIyUlRaHvwSgHrRKAFKhKJvhDhw7VaMWrCbriCQkJEdlPNUFXNrQx2gv12ZVn2LZyUz8LIJV+aeVqhYsxmbgUK2MByDBqDrlGKcSI3KRkYKBqFRSv5+HhIR6vTcB169ZNxOVRcgcllNSHV199VcT20bpGSSX//PMPNm7ciD179kDRn4+yiykbmYwi9L71HSujGrQiBlCepUTxEvv27YOvr2+9X6OkpEQE2crjKBjdjf8rLi2DrbkRXK3VsyxQRUFojgNkGI2AXLkdO3YUSRHkeaI1i+oEVna1Vsfjjz8uXKmVs3vrCsURkluWwpuCgoKE+KRxUHkXRfLzzz+LMjcdOnTAjBkz8Pzzz4skTEb90StTRK67iqEsJ4qRoAwpiu2TQwGp8isRiqegKy6K4yOoQjtdXZG5PiMjA59//rnIyKLM4bqayOlqh96DYggp8JXRfNadjsZrGy6iR3MHrJ7TDerI7yci8fbmS+jT0gm/PdZF1cNhGIVByQa3b98WF/Fcl1Uq0bJ+/XpcvHhR1UPRuuNJxuu3driAKeiWuPfKRl43iaCMpMr1k+iKhUzyCQkJojgnXZ0dO3ZM4fERjGZxJV7KdJPX21NH2pRbAK/EZQpLgrqVqmEYpnFQFi0VhP7uu+/w4Ycfqno4jJaiFQKwLkZMqpxemcWLF4uNYaoVgOUiSx1p5WoNSv5NyS5EUlYBXNTUVc0wTMOgkKY1a9YIN25D3L8MozMxgAyjqAsJKq+i7gLQzNigIkP5QnSGqofDMIyCobp/VHd27dq1oh4gwygDFoAMU05Mep4or2JsoI/mTupXAqYy7b1sxe15FoAMwzBMA2AByDDlUHs1wt/FUpRbUWdCmtmJWxaADMMwTENQ71WOYZoQTUgAudcCSC7gklKNT+RnmCrUt5sFw1QHH0c6kATCMIrgigbE/8lp6WIFc2MD5BSWiNqFAa7qV7SaYRrSAo2qNcTFxcHJyUnc5yx3piHx3FQ/MTk5WRxPdBwx98MCkGHKuapBFkADfT0Ee9rgREQazkWlswBktAJarKlmW3x8vBCBDNMYzM3N0axZsyol4Ji7sABkGACZuUWIzcgT+4EaYAEk2nvZCQFIcYDTujRT9XAYRiGQtYYWbep/Sx2aGKYhUPa0oaEhW5BrgQUgw1SK//OyN4O1ae3tmdSFkGacCcxoJ7RoU5u0B7VKYxim4bBdlGEqCcBAV82w/hEh5Ykg4YlZyC4oVvVwGIZhGA2CBSDDaFgCiBxna1N42JqBGuFcjGErIMMwDFN3WAAyjIaVgKmuHMy5KBaADMMwTN1hAcjoPIXFpbiZlKVxFkCC4wAZhmGYhsACkNF5qI5eUUkZrE0NhUtVk6jcEo5qXzEMwzBMXWAByOg8Fe5fd2uNKxnQxsMGhvp6SM4qqChjwzAMwzAPggUgo/NUJIC42UDTMDUyQGB53CK7gRmGYZi6wgKQ0XmuxGeK20A3zeymwYkgDMMwTH1hAcjoNBQ3p4klYCrDiSAMwzBMfWEByOg0FDcnyy+GkYEe/J012wJ4KTZTZDQzDMMwzINgAcjoNFfjpfIvLZytYGyomaeDr6MFbMyMUFBcimsJkjWTYRiGYWpDM1c8hlEQ56PTxW2Qhrp/CcpcrlwOhmEYhmEeBAtARqcJvSMJwM4+dtBkOBGEYRiGqQ8sANWUklIu6qtsKF5ObjHr6G0PTaaDtyRgQyPTVD0UhmEYRgNgAaiGhN5Jw9CvD+FERKqqh6LVXI7LFHFzduZGaO5kAU2mo7cdDPT1EJ2Wh5j0XFUPh2EYhlFzWACqIRvOxor2ZAs2hiG/qETVw9F69y9Z/zStA8i9WJoYIthTKmR9IoKtgAzDMEztsABUQxaMaAVnKxPcTsnBN3tvqHo4WovcXdpJw+P/5HT3cxC3x2+x5ZhhGIapHRaAaoi1qRE+GNdG7K84FCHquzGKLwB9JlKyAHYqj5/TdLqVC0AKHaDPxzAMwzA1wQJQTRka5IoRbV1FMsgbGy+iuIQL/CqSO6m5SMkuFLX/2pa7TjUdsmQa6uuJ4tYx6XmqHg7DMAyjxrAAVGMWjgkSBX4vxcrw85Hbqh6O1iXaEMEeNjAxNIA2YG5siHbl5WDYDcwwjKaTml2AX47cRhEbQJQCC0A1xtnKFG+ODBT7X+2+jjspOaoektZQ4f710ezyLzXFAXIGOcMwms6nO6/h/W1X8Mr6C6oeilbCAlDNmdzREz1bOIhyJZQVzLFdiuF0uQVQW+L/5HRvXp4IwnGADMNoMGci07AuNEbsz+zuo+rhaCUsANUcKk+yaHwwTI30xaL+58koVQ9J40nPKcSt5JyK+nnaRIdmdjAy0EN8Zj4iU7keIMMwmgfFvL+9+bLYn9rJS+vmaXWBBaCakpNZgNLybiDNHMzx2tBWYv/jHVcRncYLuyLcvy2cLWFnYQxtwszYACFe0mTJbmCGYTQRMnRciZfB1tQIz/XwVfVwtBYWgGrIhb3R+OPt47h2LL7isVk9fNDFxx65hSV47e+LFeKQqT+n5fX/tPSqslslNzDDMIwmkZxVgC/+Cxf7z/u5YfunZ8WayCgeFoBqCMVuFReW4sSWWyjMKxaP6evr4bNJwTAzMhAL+x8nI1U9TI3lTEUHEC0VgH72FZnAHAfIMIwmsejfq8jKL0aIqzXKLmaguIC7YSkLFoBqSNt+nrB1MUdeVhHO7LxT8biPowXeGC65ghftuIYojvGqN9Ra72KMVFi7s5ZlAFeOA6T6hklZBaKbDMMwjCZw6nYaNp6NBXXmfNTOHvnZRWItbNPPQ9VD00pYAKohBob66Dmxhdg/vzcamcl3hd6Mbt7CwpNXVIJX/r7AruB6Ql1VCktK4WhpDG8Hc2gjpkYUB1heD5DdwAzDaEjixztbLon96UHuiD+dJPZ7TmoBAwOWKsqAv1U1xbutA7wC7VBaXIZjG29VPC5cwRPbwdzYQFwtrTp210LIPJjQyLvuX8qw1lbk5WBOREjxjgzDMOrMyqN3cC0hC7bmRuicqS/WPq/W9vBuI81ljOJhAaimkDjpOdlfmMIjziUjNlwSLvKs4AUjAisKZd5IzFLhSDWL0PL4P211/95bEJrjABmGUXeofSU1OyBebtcMMWGp0NPXE9Y/bb5QVzUsANUYB3dLBPWRYh+O/H2jirv3ka7N0LelkygQPf+v8ygs5lY5D4L6KoeWZwBrawKInPbNbGFiqI+U7ALcSs5W9XAYhmFq5N0tl0VYU2dvWxiFycRjbXq7izWQUR5aIQAXLVqEzp07w8rKCs7Ozhg3bhzCw6U08tpYv349WrVqBVNTU7Rt2xY7duyAutFltC+MzQyREp1dpSwMXRV9PikYduZGol7S4j3S1RNTM+ej05GRWwRrU0O08bCBNkP9jTv5SCL34PUUVQ+HYRimWv67nIA9VxNhqK+Huc3ckBqTDRNzQ3QezfX/lI1WCMCDBw/imWeewYkTJ7B7924UFRVhyJAhyMmpOQPy2LFjeOihh/D444/j3LlzQjTSdumSFISqLphZGqPzSKkNTuWyMISztSkWTWgr9r8/eAsnOeC/VnZfkYKK+7dyhpEOBBX3D3AWt3uvJqp6KAzDMPeRU1CMhVuljh9PdvfBnQOxYr/zSF+x9jHKRStWwZ07d2LWrFkICgpCu3btsGrVKkRFReHMmTM1/s8333yDYcOG4dVXX0VgYCA++OADdOjQAd999x3UuSzM6R1Vkz6GtXET/YIpzOuldRcgyy9S2TjVHbrKJAYGukAXGFT+OSlZiI8LhmHUjcW7ryMuMx9e9mbommco1jgu+9J0aIUAvJfMTKnOm719zYH+x48fx6BBg6o8NnToUPF4TRQUFEAmk1XZmqwszCSpLMzFfdHISKxa/+/dMUFoZm8uAmkXbpGuppiq3EnJwc2kbOFmoNhJXYDqRlK7u+LSMhwMT1b1cBiGYSq4HJeJleVVLN7s7Y8r5da/XpP9uexLE6F133JpaSleeOEF9OzZE23atKnxeQkJCXBxqWoJovv0eG2xhjY2NhWbl5cXmgqfto6iNExpSRmOrL9R5W+WJoZYPLUd9PWAjediseW8dCIx91v/uvjaw8bMCLrCwEDnKp+fYRhGHRLy/rcxTNyObOOKotBUkeTo09aBy740IVonACkWkOL4/vrrL4W/9oIFC4R1Ub5FRzdtf8Jek/yhb6CHyEupuBNWNbC/o7c9nhvgL/bpxOIOEFWRCyC5W1RXkH/eA+HJotAqwzCMqqH6tRdiMmFlaognWrgh6nKaWNt6TpLWMKZp0CoB+Oyzz2Lbtm3Yv38/PD09a32uq6srEhOrWkXoPj1eEyYmJrC2tq6yNSUUG9FugGR1JCtgyT2lX54b0EJYuHIKS/Ds6rMoKOYeikRmbhFOl9f/0zUBSG3hKFM8M6+oogg2wzCMqohOy8UXu6QqHQuGtsLl7VJf+3YDvcQaxzQdWiEAqdAtib9NmzZh37598PV9cPp49+7dsXfv3iqPUQYxPa7OdBrhAzNrY2Qm5eHCvqoWSEMDfSyZFiIW/MtxMtEvmAEOXE8SroaWLpaiiLYuYaCvJ7KeiT1X2A3MMIxq1+o3N18SNf+6+tojILMMmcl5MLc2Fmsb07Toa4vb948//sDq1atFLUCK46MtLy+v4jkzZ84ULlw58+fPF9nDX375Ja5du4aFCxciNDRUCEl1hmoC9hjfXOyHbr+DnMyCKn93tTHFl1PaVZjZd12uOaZRV9hdLnx0zfonR/65yQ3OXUEYhlEVm87F4tD1ZBgb6uPdwQE4869k/es+oTmMTQ1VPTydQysE4PLly0VMXr9+/eDm5laxrV27tuI5VBYmPv5uIeUePXoIwbhixQpROubvv//G5s2ba00cURcCurrC2ccaRQUlOL7pbp9gOQNauWBOb8kK+trfF0V2sK5CHVIOXk/WqfIv99KnpROMDfRxJzUXt5I5NpRhmKaHuhK9v+2K2J8/0B9xhxLEGubia42ALjWHXjHKQysEIFk1qtuoNqCcAwcOiPqAlZk8ebLoGELlXShxZMSIEdAEqEdin6ktxX74iQTE35LK3lTm1aGt0M7LVsR+6XI84Ok7acjKL4ajpTHae9lCF6Es8a5+UkkkLgrNMIwq+GDbFdGJqZWrFcZ4OIi1C3pA7yktxZrGND1aIQB1EbpqCuzpJvYP/RVepU8wQSb27x4KEW3PzkVliJNPl92/1BWD4uF0lcGtJevn3qtSNxSGYZimYt+1RGw5HydKlX0yvi2OlZcya93DTaxljGpgAajBdB/XXPRMpD7BVw7fX/vPy94cX09rL/b/OBGFv8/EQJcgK/Dea+Xxf+UCSFcZUJ4IEhqZhvScQlUPh2EYHYG8UAs2hon9x3r6wuB2jlizaO3qNk6KZ2dUAwtADcbMyhhdx/iJ/RNbIpCXXVhtPOALg6TaSm9uCsOl2PvdxdrK9cRsRKflCWtob39H6DKedubC9UKG4v3hbAVkGKZp+HDbFSTKCuDraIFnevji5NYI8TitXbSGMaqDBaCGE9TbHQ6elijILcaJzdKJdS/PD/AXFqCC4lI8/ccZnbEA7bwkZUD3bO4Ac2POMJO7gbkrCMMwTQFdbK4/EwM9PeCzScE4v/2OWKscvSwR1If7/aoaFoAajr6BPvpOkxJCrhyNQ+Kd+/sT6+vrYfGU9qJfcEx6HuavPS/q4mkz1PVi7ekosT+6nbuqh6N2XUFyCopVPRyGYbQYWX4RFmyQXL+ze/iiWZkhrhyTKnFQEiOtS4xqYQGoBbi1sBWlYVAGHFoTjrJqxJ2NuRF+mNERpkb6og7Tpzu1u0j0vmtJiMvMF0WxR7SVkmV0nWBPG/g4mCO3sKTCOsowDKMMPtp2FQmyfHg7mOOVwS1FsiKtUQHdXMWaxageFoBaAhXSNDI1QFJklrAEVkegmzU+myQViV5xKKLCQqaN/H5CKjA6pbMXTI0MVD0ctUBPTw8TOkgtEjee062EIIZhmg4yMqwNjRau388ntcPt04libTI2NUD38kYGjOphAaglWNiYoOtoKSGEikPnZVUf5zemnbsowkm8uekSjt9KhbZxOyUHh2+kiMlnehdvVQ9HrRgfIsXdHLuVijgdLhDOaCa5hcX48r9w9Pt8PyYtP4YFGy/ilyO3cfhGsnA5MuqR9fv6hoti/9HuPmjraFnRsKDLaD+xVjHqAQtALaJtPw84eEgJIdV1CJFDWcEUF1dcWiaSQkgwaRN/llv/+rV00rnevw+CSgNRD07qCEdtmRhGU0o6bT4XiwFfHMS3+26KrjahkelYcypadJeY8fMp9P50P07dTlP1UHWehVsvIz4zX4SbvDYsAMc23RJrEiUr0hrFqA8sALUtIeThALF/9Vg84m9m1OgK/HxSsOiMQVdrj686jYxc7cgMzi8qEVlnxIzubP2rjolyN/DZGO4NzKg9F2MyMHH5Mbyw9ryIKfOyN8Piqe3wzbT2eG5ACwwLcoW7jamYy2b8fBJ7you/M03PjrB4cWFJ+R1fTmmPjKhsXCtP/Oj3cIBYoxj1gX8NLcOtuU1Fh5CDa8JRWlJa7fMoLm7FzI7wsDVDREoOnvr9jBBPms4/F+LEQuBpZ4a+LaXix0xVhrd1FclA1Bf4Qozu1IVkNA8KUSHxdzYqA+bGBnh1aAB2v9gX40M8Mba9B14eEoDvZ3TE3pf7YWB5qaun/jiD9aHRqh66zpGUlS9qzRJz+zVHe08bHFwdLu637ukGVz8bFY+QuRcWgFoIBdmaWBgiNTYHF/fXHOzvbGWKnx7tJHrFnrydhvl/nRPlUzSZP8rdv9O7eut067fasDI1wtAg1worIMOoIzeTsvHU76EoKilD/wAn7H+lH57p36LapC4zYwMhBMm6TSWuXv37In44WHMYDKNYyJPwxoYwpOcWobWbNeYPbImw/TFIi8sRa1E3TvxQS1gAaiFmlsaiTRxx6p/byE4vqPG5lBlMlkBjA33supwoEkM01S14ITpDWLTos0zpJLk5meqRZwNvvRCHwmLNFv2M9pGSXYDZq05Bll+Mjt52WP5IR7hYm9b6P0YG+vhicjCe7CMlwy369xqLwCZi7eloUXqL5t7FU9ujMKtIrD1Ej/EtxJrEqB8sALWU1j3dRZPtooISHP1barxdEz2aO2LJQ+1F3Aal7n+6UzLba6r1b2SwGxwsOdOsNnq1cISzlQkycovExM0w6gKFosz5LVS0caTi9StE/dK6lXKi+Ob/jQgUyQfEF/+F41rC/cXxGcURlZqLD7ZdEfvkog9wtcKR9TfE2kNrUGAPrsOqrrAA1FL09PXQ96EAUQrl5pkkRF6uvdzLsDZuWDShrdj//uAt/Hio+rZy6kp8Zp6wZhGPdGum6uGoPeQel5eEYTcwoy6UlpbhpXXncS4qAzZmRlg5u3ODLubm9m0uOt+Q+/iV9RdQpOGhLeoKfa/P/3UOOYUl6OJrj8d6+SLyUipunU0Sa49YgzgUR21hAajFODWzQnB/L7FPHUKKCmtP8pjauRleH9ZK7H+04yr+PClZ1DSBj3dcEwHgnbzt0KGZnaqHo1FuYOrXmaYj/aEZ9eabvTewIywBRgZ6onNRcyfLBr0OWQI/Ht9GiMhLsTJ2BSuJJXtv4Hx0BqxMDYXrt7S4VCQfEsEDvMQaxKgvLAC1nC5jfGFpZwJZSj5Cd9x54POf7uuHp8pjaCge8NdjD/4fVXMiIlVk/9IV58IxQWLyZx4MuWraeFgLKwlbARlVExaTie/23xT7iyYEo5ufQ6Nez9naFAvHtK4QluwKVvy8e/f3aisqSoRuv4Os1Hyx5nQZ7avqITIPgAWglmNsaojeU1uK/fP/RSE1LrvW55N4emN4qwoR+O7Wy2rtDqasZSo8Skzv2gxtPLjUQH14qIvkLl959I7GZ4AzmgslIr369wWRwTuyrRsmdVRMEte49h7sClYCmblFeHHteVFQnhLuRgW7IzU2G+d3S+1Fac2htYdRb1gA6gB+7Z3gE+wo4muoLlNZaVmdRCAVWZW7g5eWX+mpY+LHtYQs2Job4eXBUuA3U3eobIaDhTFiM/KwPUwq2MowTQ1Zkug8trcwxntjgxT2uuwKVjxUJWLBpoui24evowXeHR0k1hRaW2iN8W3nKNYcRv1hAagj9JnWEoYmBoi/mYmrx+PrNHFSkdWXBkvWw893heOr/8LVqkQMlYr4cvd1sf/KkADYWXCpgfpC2ZWP9vAR+ysORajV78voBpfjMrGs/ALzvTFBcFRwBv+9ruCI5Nq9IEztrAuNrojTXDItBBYmhlLnqVuZYo2Re5wY9YcFoI5gZW+KruUxGcc23ESurG5B/88P9K9IDFmy7yZeXn8BBcXq0THk853hyMovRpC7dYUrk6k/M7p5w8zIAJfjZDh2q/ZscYZRJOSSfXX9RdGXnFq6jQpWTskQcgX3C3ASruBPd15TynvoAuEJWSIsSH7R3dbTRqwlxzZKAp7WGFprGM2ABaAOEdzfE45elqIx94NqA1aG2vp8OK6NKB2y8WwsZvx0SuVZo5R5tu6M1O7p/bFB3PWjEZDlVF44+wc1jvdktI/lB27hSrxMhHB8MK6N0hK46HXfHBEoap1SwfvTd9KU8j7aTE5BMeb9SS1DS9Hb3xFzektx4lTzj9YUWltojWE0BxaAOgQ14u43vZXIlr1+KhFRD6gNWJlHunlj5azOsDIxxKk7aRi/7Kho1aSqmn/z/jgjApAndPBAR297lYxDm3iit59YHA9dT8bVeM6WZJTPjcQsfLtPuhBdODoITlbKLd7u72KFqZ2lslgf77jK4Q71gL6rtzZfEv3DXa1N8fXU9tDX1xM1/26cThRrSv9HWok1htEc+NfSMVx8rCtqAx5YHS6qtdeVPi2dsHFeD3jZmyEyNVeIwL1XE9GUZOQWYubPpxCXmQ8/Rwu8PVKK7WEah5e9OYa3ldxv6pz1zWiXoCCX7MBWzhjb3r1J3vfFQS1hbmwgCk1THBtT91Zvm87FCk/Ltw+HiOLchfnFIvFDXvPP2dta1cNk6gkLQB2tDUhxGlSv6eQ/EfW+it48r6foz0nxd4//Gor/bQoT7gFlk1dYgid+DcWNpGy4WJvgt8e7cOKHApGX/qGOKnEZeaoeDqPFUCjJydtpMDXSb9LanZQQInddfrbrGvfBrgNX4mR4p1LcX2cfyeNCvX6z0vLFWsI1/zQTFoA6CNVn6vuwVDLl4t5oJEXWz+VHV3+r53TFE72kk371ySiMXHIYZ6PSoSyoRt1za84iNDId1qaG+O2xrvC0M1fa++kiwZ626OZnLwLyVx6VGrkzjDJqyJELVp5kRtbnpuTJPn7C3UxeDHn/cKZ6svKL8Mzqs0Io9w9wqrhITLwjw8V9Ugx23+kBXPNPQ2EBqKN4t3GAf2cXEUe3/49rKKlngVQTQwO8Nao1Vj/RFe42priTmovJ3x/HF7vCFW4NpMnnjY1h2HM1CSaG+vjp0c6iiwWjeJ7q07xC1Cdl5at6OIwWQpa31JxC+Dtb4olekqBoSqhsCbmCiSX7biAzr6jJx6AJUE0/Kp59OyVHzPFfTZHi/mitoDWD1g5aQ7yDGtexhVEdLAB1mF6T/WFiboiU6Gxc2CNdzdWXHi0c8e8LfTCuvbuo4k8FXft8th8/H7mN/KLGl4s5E5mGUd8ext9nYkSSwrcPhYim44xyoFIZwZ42orn7l7ukGosMo8js/dWnpG4RlPVrbKiaJYiy3ls4WyIjtwjLDqhnkXtVs/zgLZExbWygj6XTO1SE29BakRqTDRMLQ7GGMJoLC0AdxtzaGD0nSd0+Tm27jYzE3Aa9DlXZ/3paCJZP7wAfB3Nxdf/Btivo/8UBYUlqSN1AWX4R3tochknfH8f1xGzRIWDZ9A4YEuTaoDEydYNisd4ZJSXWUJmdS7GZqh4SoyXQBeKbm8Kk7P0Qj0b3+m0Mhgb6WDC8VUUbRI55rcr+8CR88V94RZmtkGZ2Yp/WCForiJ4T/cUawmguLAB1nFbd3eDZyg4lReVm/Qe0iasNyiLd/VJffDKhLdxsTEWrIEoQ6fjBHhFHsuV8rBB2tRWFPROZLtrODfryIP44ESUWi8kdPbH3pb4Y1kY5RWKZqnTyscfodu7iu39/2xUul8EohN+P3xHFximG938jA1U9HAxo5Sy8CRRisri8oxAD3EnJwfw158T5/3DXZphWXmSf1gYRLlRUKtaMVt35YlzT0Svj2b3ByGQy2NjYIDMzE9bWmpsCL0vJw5r3T6K4sFQkh7Tp49Ho1yT3L1n/qL1YguxuLBm1D2rrYSMselamRrAyNRRdKKgYLIm/3MK71kLqM/nR+Dbo0dyx0eNh6gf1Bh745QFR9JUsryPKS8QwTENIyMzHoK8OIrugWBSVp7qi6sC5qHSMX3ZMhJf8O7+PzscWU/w2lfcirwtVelgzp1uFm/7SwRgcXHNdtHt76O0usHY0gyYj05L1uzFw6g4jTuRuY5uLiu7U0ocSRBrbzod6zD7WyxezevjgYmwm/rucgP+uJIri0WejMmr8PztzI3T1dUBPf0dh+aPXYZoeD1szPNmnOZbsvSEyNslawr8F01De3XpJiL+QZrZ4WI3aNpJrc3gbV/x7KQGf77omEsx0Oenj1b8vCPHnbGUiQnrk4o/KvRzbeEvsdxvrp/Hij5FgAcgI2vb3xM0ziUiIkOHAn+EY9WywQmpzUdZYey9bsb02rJVoxE7WPqohSCUG6JYWhmb25uje3AEtna3E/zCq5+m+flh3Ohox6XkiqeeZ/lK8KMPUB7r4o2QCQ309LJrQVu3O71eGBoiLU6oycOp2ms4mmS3ec10UxyYvzfJHOoiaiQQ5CWlNoKYBrn42CO7H7d60BY4BZAQ0KfefEQh9Qz3RIo5axSkDPydLjAp2x0NdmgkL08tDAvDu6CDM7umLVq7Warc4NJSy0lKk/fYboh57HAW3pCtnTcPc2BBvlAfKU1xmYiVXPsPUBbrIe2fL5Yr6e3SOqxvNnSwrWsR98q9utojbcCYG3+6TsqE/Ht+2SnvN6ycTxJpgYKiPATNbQU9L5miGBSBTCXs3C3QeKRV3PrzuOnJlhaoekkZSFBeHqNmPIfHjRcg5dgyJH30MTYVadJHbjmIz39x0SScXR6bhfPnfdRED7O1gLoo+qysvDPQXscgUnkLWSl3iZEQq3th4Uew/0785JneSxDBBa8Dh9VK/5s6jfGDnaqGycTKKhwUgU4WQIc3g6GWJgpxiHPpLKgPA1A0SR5lbtiBizFjknjwJPTMzwNBQiMDcs2ehiVAYAAXtUy2wPVcT8dtx7pzA1L3m36/H74j9j8a1VesYUnJ3PtHbt6JQNXUe0gWoyPNTf5wRPZlHtnXDy4OlDlHy+ezgmnCxFtCa0H6w+sRuMoqBBSBTBQMDMvMHClfsrbPJuBGqW1fDjXH5xr36GuJefwOl2dkwa9cOfps2wnb8OPH3lO+WQlMJcrfBghGSK/ij7VdxOY5rAzK1QyWdFmy8W/Ovl7/6Z/KTi5qqE0Qk5+DPk1Kxam0mI7cQj606LYphU4z2l1PaVQnBuXkmCRHnksVjtCbQ2sBoF/yLMvfh5GWFDsOlMg2H/mJXcF3IOXIEsm3bhMXPaf7z8P7zDxj7+MDhqac03gpIUDb3oEBnFIqezOcU3u6P0S6W7b+Fq/EykdX/phrU/KsLVJbqpcFSizgqgpySXQBtJbewGI//GiosgJTx/+PMTlUstDTnH1oj1UbsONxbrAmM9qE1AvDQoUMYPXo03N3dhdtq8+bNtT7/wIED4nn3bgkJCU02ZnWm03AfOHhYIj+7iF3BdUD233/i1m7KZDjOnQs9QynB3tjTUyusgHRufD6pHVytTYWF5N2tUmA/w1Tn+qUeuwQleDlYmkBToOS0Nh7WojrBJ/9egzZCha/n/nFW1F2loty/zOoMJyuT+1y/+TlFcPC0RMfhPiodL6M8tEYA5uTkoF27dli6tH6LbHh4OOLj4ys2Z2dnpY1Rk6CMr4GPsiu4LpQVFyN7z16xbzVkyH1/1xYrIPUC/Xpae1E0l3ozbz4Xq+ohMWpoWXpx7XnR9m1UsJtIItIkDPT18P7YNmKfjnESSdpW6++V9Rdw8HoyTI30sXJ25/uKX1d2/dIaQGsBo51ozS87fPhwfPjhhxg/fny9/o8En6ura8Wmr681X0mjcWrGruC6kBsaipKMDBjY2sK8U6f7/q4tVkCC+rfKszmpzR9ZexhGDhUNJ7ciWYop8UMRtUSbmg7N7DClk1Tr7p0tl4SY1QbIsvfeP5ex9UKcqMm4/JGOVcq9EOz61S10Xu20b98ebm5uGDx4MI4eParq4ai3K3hNOJcBqYascvev5aCBFa5fbbUCEs8N8EevFo6iNMyjv5wSsV4Ms/9akujfTVBCgY25ETSV14e1Eu5R6l28+qR2ZL4v2XsTvx6PBGly+n36B1T1drHrV/fQWQFIou/777/Hhg0bxObl5YV+/frhbC2Lc0FBgegfWHnTKVfwOXYFV5f9m7V7j9i3rsb9q41WQHKT/TCjIzo0s0VmXhFm/HxSdHhhdJfU7AK8+rdUS+6xnr7o2UL9s35rg+IWXx0qlUT5fFe4+HyaDBVyp04fxMLRQRjb/v5+7zdOJ7LrV8fQ2V84ICAATz31FDp27IgePXrgl19+EbeLFy+u8X8WLVokmkfLNxKNuuIK7jhCuhok90BOhmZPhook7/x5FCcnQ9/SEhbdutX6XGEF1NMTVsCipCRoMhYmhlg5uwtau1kjJbsQ0386iei0XFUPi1FRXNnrG8JE1qy/syVeG3a3lpwm83BXbwS5W0OWX4yPd2huQsg3e24IEUuQqH20x/2Wvez0AhHmQ3Qa6cOuXx1BZwVgdXTp0gU3b0rtcKpjwYIFyMzMrNiio6OhK4h4kGZWKMgtxr7fdbNdUnVk7Sp3//bvDz1j41qfS1ZAk0Cpnl7uqdPQdGzMjPD7413Q3MkC8Zn5eOTnk9wuTgf5Zu8NUSScioUvntperQs+NyQhhFymG87GYPvFeGgSNEd/9V94heWPhHl1/bzpefv/uCrmdmdvK3QYJsV9M9oPC8BKnD9/XriGa8LExATW1tZVNl2BioAOmtVauAWiLqfhypE46Do0cWbt3i32rYYMrtP/WHSVrIQ5J45DGyBX2Z9PdIOXvRkiU3MxbulRXODEEJ3h37B4IQCJD8e3QRsPG2gTHb3tMLdvc7FP7dI0xcpNcxPVMlxS3t/3fyNaYV6/+8UfcflwnJjTRbgPzfFc8Fln0JpfOjs7Wwg42ojbt2+L/aioqArr3cyZMyue//XXX2PLli3C4nfp0iW88MIL2LdvH5555hmVfQZ1x97dAl3H+on9I3/fRGZyHnSZ/EuXRd9favlm2atXnf7HoltXcZt74iS0BVcbU6x+ohv8yi2Bk384jnWndcc6rqtciZPhpXUXKuL+plTqIatNvDi4peiHTbUBn//rnOhyou4u+Q+2XcXS/bfE/bdGBuLJPpKIvReaw49ukERit3F+oh88oztojQAMDQ1FSEiI2IiXXnpJ7L/zzjviPtX4k4tBorCwEC+//DLatm2Lvn374sKFC9izZw8GDhyoss+gCbQb6AW3FjYoLijBvt+uokxLSiQ0Kvu3Tx/oU9/fOmDWsROZU1EUE4PCmBhoC1725tjyTE8Mbu0iCs2+tuEi3tocJvYZ7YOSIub8Foq8ohL09ncUFiZtxchAH0umhcDK1BDnojKweLfkUlVH8otKMO/Ps/jl6G1x/93RrfFEb+mivTqhuPfXK2Iud/e3RbsB2ingmZrRK+NgrgZDWcCUDELxgLrkDqarxr8+PCUmjh4TWyBEB5uE02kTMWw4CiMj4f7lF7AZObLO/3tn2kMiecTtow9hO3EitAlaVL4rzzikmYUyhT+dGAx/Fw4q1xZI1FO856nbafBxIOHfS6NLvtQVigF8ZvVZERP4+2Nd1a6/MYnyJ34LFSKV4jE/nxxcbbavnHO7o3Bsw00YmRhg2ttdYO1Yt4tYbUGmo+u3VloAmabDxskMvSZJ8SQnttxCaqzulQApuH5DiD9K/LDs269e/2te7gbO0SI3sBwqIUGFon9+tJOwmJyNysCwbw7j/X+uQJZfpOrhMY2E3J8vrD0nxJ+liSF+erSTTog/YmSwm2gVRxc2L647jyQ1SniiMkwTlh8T4k+enFWb+EuJyRZzN9FzUgudE3+MBAtApkG07uUOn7YOKC0uw+5fLqOkqFQn3b8WPXvCwLJ+cTPycjG5J05obTb1gFYu2PF8bwxp7SI6KZBLasAXB7AuNFpYCRnNg8Tf82vOYUdYgrAwffdwCFo465Zl951RrdHSxRLJWQWi9BGVvlGHAtwk/igJy9PODBvm9kBXP4can19cVII9Ky+Ludsn2FHM5YxuwgKQaRDU4qn/jECYWRkhNTYHJ7ZGQJfIPnRI3FoNGlTv/zVr315YDql+YOFtKVZHG6G4wBUzO+G3x7qIBBGqF/ja3xcx5OtD+P34HWQXFKt6iEw9xd+/lyTxR4XA+93TSUIXMDM2wI8zO4lWdzeSsjH9x5NIyylUmSv+o+1XMHvVaWTkFqGdly02zeuJFs6Wtf7fyS0RYs6mubv/I600sl0foxhYADINxtzaWIhA4vyeKMSEa1fj9JoozctD/tWrYt+8q+TOrQ/6pqYwK09WyjlxAtpOn5ZO2Dm/D94cEQgrE0PcTMrG21suo/vHe7Fw62Xc4i4iai/+nltdVfz1b6V74k+Ot4MF1jzZDc5WJghPzMLDP55AehOLQCpHQ9n2Px6WLiBn9fDBuqe6wcnKpNb/i7mWhvN7pAz9ATMCxRzO6C4sAJlG4St3IZQBe1ddQUGu9sd55V0MA4qLYejsDCOPhrlPLLrL3cDaFwdYHcaG+pjTxw/HFgzAe2OChEUwq6AYq47dwcAvD2LI4oP45N9rOH0nDcVqXmZDlyAr7dw/zmLn5QTxG/4wU7fFnxxfR0kEkuC6lpAl3MEZucoXgRQysulcDEYsOSzqbVK8HwnyhWOCYGJYewFu6vG791fpwrV1b3fh/mV0m+o71zNMPaAg4tjwdJEdfHDNdQx5PAjaTN45qV+0WYcODXafyC2HuSdPin7Cevq6cS1mZWokWlHN6OaNo7dS8OuxO9gfnozridli+/7gLdiaG6Grrz2CPW3R1sNGbHYWbKloam4kZuHpP87gVnKOEH8rdNTtWxPNnSyxZk5XTFtxAlfiZZjyw3F8NaW90ophX4rNFBbz0Mj0iiLVSx4KgYdt3RI4qNUbtXyjJL6eE6svCs3oFiwAmUZjbGqIQbNbY+MXZ0VDce82Dgjo6gptJfesJADNO0hu3IZg1qYN9M3NUZKZiYLwcJgGSq50XYGyhXv7O4mNLCcHryeLYHYSgxTPtOtyotjk0CJHVpdmDubwtjeHt4M5PGzNhQXGwdJY1GpjFMc/F+Lw+oaLyC0sEfFuS6d3EIKDqQolwaye0w0P/3hSXMCMXXoU8/o1x7MDWjzQIlef8i5f/Hcdf52OEhnIZkYG4vWf7ONX5+M+/GSCmJv19PUw6LHWYs5mGD4KGIXg6meDziN9cOqf2zi4JlzcpytNbYOsdXnnpe4HZiEdGvw6ekZGMOvcCTkHD4lyMLomACtja24sSlbQRu7f89EZopxFWGym2G6n5CA2I09sqKFVt72FMZwsTYT1UGxmxuLW2sxIuMno1trUUNza0mZuLO4bsnC8L7Fg0b9XsfLoHXG/R3MHYWVytKw9tkyXaelihZ0v9Ma7Wy5je1g8vt13E/9dThR1+MiK3VAoVnb1ySisPxMtupAQY9u7443hreBmU/e5VfLMhIt9mqNdfbWrXR/TcFgAMgqj4zBvRF9JQ/ytTFEaZsIrHaCvZQtswc2bKJXJRPs301YBjXot6gtMApDKwTjMnqWwMWoyJMg6+diLTU5mXhGuxcsQmZaLqNRcRKXliv2EzDyRWUxlZigTsyHZmFSr0M7cWAhIR0vp1t7CRAT4U4s7F2va6L6pcINqKxRbtu9aEj7afhURKTniMbJkvTwkAAb6nCX6IEggk5V0ZFg83t58SSSHUF9sSoAaFeyOIUEusDY1qlMnj91XEvHnyUiciEireLy1mzXeGxuEzpXOi7pQUlIq5uKi/BLRwanjcJ8GfT5GO2EByCgMEnvkXlj74Wkk3pbh9PY76Dqm+jZEmkreOanXtFlwsLDiNYaKvsCnT6OsqKjRr6etkAWP6ppVV9uMagqm5xYiObsAKVmFyMgrFC5kciun5xZBllckBCQVoc7MK664Ly9BQ5YV2khU1gaFepIrlOqsedmZi9tmDhbwdTSHr6Ml7MyNNLacBsX6vb/tCg7fSBH3SQh/PL4thgRpbxiHshjR1g3d/Bzw7tbLwo1+IDxZbMYb9YUY7OZnL4SghYkhLE0NYWSgJyx9FN93KVaG64lZKC6vk0m6m+ppTu/WDH38nRokxEO33xFzsbGZFKZDoRcMI4cFIKNQrB3M0O/hAPz382Wc+fcOvALtRZ9JbSGvPP7PrBHxf3JMWrWCvo0NSjMzkX/5sqgPyNQPWtAcLCkO0ARwrV9pExKDGXlFooRHarkFkeKtyKpIhX4TZPlIyMxHUlY+ikrKEJ+ZL7bTd+4vd0TuZF8nSzR3soC/s5WoxUZbM3tztbWgXUuQ4bfjkVh7OlpYUanEy+xePni2fwuRrMM0DLIif/tQCOYP9Bft47ZdjBM1A/dcTRTbg6ALjSmdvTCtsxfc65jgUR1xN9LFHEz0mx4g5maGqQwLQEbh+Hd2QdSVVFw7niDcD1Pf6gJTC+1YUHLPnRO35h0aHv8nhzJ/Lbp0Qdbu3SIOkAVg00HB8xXC0an255KVkQRiTHouYtLzxBadnovI1BzcTs5BXGY+ZPnFoiwHbZUhtzFli1L3CIoVo83f2VIUyVaFMMwtLMa2C/FYczpKxFnKoY4tb44MFDXuGMVAFwDzB/mLLTwhS8QHUjxrTkExsvOLhRU6v7gEPg4WaONujSAPG5FB7G5j2mhrMpV82f3LFZE00qqHG/w7uSjsczHaAwtARin0ntoScTczIUvOw4E/r2HonDYa6yKTU5ySgqKoKOEPVJRYM+/cWQjA3PLSMox6Whkp25i2kGb3Z8LmFZYgMk0Sg+TOu5mcjRuJ2YhIyUZ+USmuxsvEVhmTcmHo70JWQ0v4OFrAz9FC3FKPXUVBlj2y9FHv3pMRaThyM6XC/W1IGaGBLpjV00e4LRnlEeBqJbamiuekOVde8qX3FP8meV9G82AByCgFKjNA9QA3fnYGt84m4/LhOLTpU3Nzck0q/2Li7w8DK8VM5mYhkpDMP39BTNyaLpJ1EWoP1srVWmz3Wg7JWkgJARTbRbF24SQMk7NRUFwqasfRVl1CgYedGdysTeFmawo3G1MhPq1MjETSCsWO0T5RVFqK4pIy4dImy1JcZh7iMvLLLZW5wiJJ1snK+DiYY2rnZpjU0fOBnSMYzYPmWppz9Q30MPjxIC75wtQIHxmM0nDxsUa38c1xbMNNHFl/A27NbeDgUXufSnUm7+w5hcX/yTENCICeiYmoB1h45w5MfH0V9tqM6i2HVLeQtsGtXapY5UickZXwelIWIpJzcCclB3dSc0T8YYqIQyyAVGyo8VgYG4is6i6+9sLSF+Jly8kAWkpKTDaOrLsh9ruPby7mYIapCRaAjFJpP9ALMdfSEXU5Fbt+vITJCzrDyEQxBVKbGrmbVhHxf3L0jI1hGhQkkkvyLlxgAagDUOwfxdrRNqiSMCQoWzkyJVdY8hLKk07k5W6odV5WflFF/JjcjUvxjPSaZIl0tzET1kNKHvCwNUWgm7UoIcL1DrWfooIS/PfTJZQUl4pi/O0GeKl6SIyawwKQUSqi8vysQPz14SmkJ+Ti8Lrrogm5plGan4/8K1crWsApErN27SQBeP48bMeNU+hrM5oFlQhp62kjNoapD4fXXhdzrIWNMQY+GijmXoapDb4sZJSOmZUxBj8WBOgBV4/G4/rpBGga+WFhQFERDJ2cYOSh2FhGeUJJ3oWLCn1dhmF0g+unEnD1WLyoV0lzLc25DPMgWAAyTYJngB06jZCq0B/4IxwZibUX3lU3civi/zooPFHDrH07cUs9gUtzpC4MDMMwdYHm0gN/Sq3eaI71COCezUzdYAHINBmdR/iIotAUq7Lrp0soLiyBphWANldgAogcIxcXGLq5Udoo8i5dVvjrMwyjndAcuvPHS2JOpblVfpHNMHWBBSDTpK3iqDSMmZURUqKzcXi9lK2m7pSVliL3/HmlxP9VjgMkKA6QYRimLhxedwOpMdliTqW5Vdt6rzPKhY8WpkmxsDXB4NlSPOCVw3EIP6n+8YCFERGiXZueqSlMW7VSynvI3cAsABmGqQs0d145EifmUor7o7mVYeoDC0CmyfFqbS/cwcSB1eFIi8/RiPZvZm3bQs/ISLkWwAtSQWiGYZiaSIvLEd0+iM4jfUXPdYapLywAGZXQaaQvPFvZoZjiActjWNQVuVVOmb16qRYgicuStDQURUcr7X0YhtFsaK6kuL/iwlIxh3LcH9NQWAAyKoE6EZDbwtzaWFzNHlwdrraWr7zzF6q0bVMG+sbGMGkdWGEFZBiGuReaI2muTI/PgbmNVF6Lu7owDYUFIKMySPwNeSJI1K6ieBbqYaluiBZtt25VcdMqC3N5PcBzHAfIMMz9XD4UK+ZKKvJMSR80hzJMQ+FOIIxK8Whph27jmuP4pluiS4iTlxVcfNWnf2XeRak4s5F3Mxg6OCj0tQtKCvD39b9xIPoAWtq1xFAfW1CEIVsAGeZ+SstKcTH5InZH7sb19Ovo79UfE1tOhImBbiQ/JNzOFFm/RLdxfmLuZJjGwAKQUTkhQ5oh8bYMEeeTsXNFGKa82RlmlupxZSu3xsmtc4qgsKQQG25swE9hPyEpN0k8diL+BLZnlmE5JZ1cu4LLMWcR5KmckjMMo0lcTrmMLbe2YG/kXiTlSeeL/Jz5+dLPmNN2Dib4T4CxgXrMGcogL6sQu1ZcQmlJGfxCnBAyuJmqh8RoAewCZlQOddYY8GggbJzNkJ1egP9+uozS0jKtTAD559Y/GLFxBD4++bEQfy7mLpjfYT5G+o1EvoMF0iwB/ZIyvL/qUWyP2K6Q92QYTWVbxDY8vONhrLm2Rog/CyMLca7QOUPnDp1DH538CCM3jRTnljZCc+F/P18Wc6OtizkGzgxUeDciRjdhCyCjFpiYGWL4U23x96ehiLmWjlP/RKDb2OYqLwAtdwErQgDuidyD/x35n9h3Nne+z3JBlsEr+2cBR86hRWypeG4ZyjDKb1Sj35thNA0SdG8dfUu4fvt59cPklpPRza1bxfkys/VMbLyxET9e/BEJOQnifDE3NMdA74HQJk5tjRBzoqGxPoY91QbGZrxsM4qBLYCM2uDgYYn+j0iFls/8G4nbF5JVOp6CmzdRmp0NPXNzmPj7N+q1YrJi8M7Rd8T+1ICp2DFhB6a1mlbFbUX7Xt0Hif1+mW5i4XvzyJtaa9lgmLqIv0ktJ+Gb/t+gj2ef+84XOod2TNwhzini7WNvIzY7FtoChcWc2Rkp9vvPaAUHd0tVD4nRIlgAMmpFyy6uaNvfU+zvXnkF6Qk5qnf/UgFow4ZfdReVFOHVg68iqygL7Zza4fUur9cYuC7vCOIdVYBJ/hPFAkgLIYtARlegY50ufOjYJ6vf293ehr5ezUsVnUt0TgU7BSOrMEuca3TOaTo09+1ZdUXsB/f3RMvOrqoeEqNlsABk1I6ek1qIxuZF+SXYsTwMBXnFqq3/10j37+Kzi3Ep9RKsja3xeZ/PYaRvVGtBaBgaoiQlBW80e1wsgHJL4M7bOxs1DoZRd+gYp2OdQh+mtJyCt7q9Vav4k0PnFJ1bdI6FpYTh67NfQ5OhOY/mPpoDaS7sMamFqofEaCEsABm1w8BAH0PntIGlnQkyEnOxZ+UVlKkgKeRuAkjD6//ti9qH36/8LvY/6vUR3Czdan2+PvUbbt1a7OefOy8WQBKBtCC+e+xd4UpmGG0kOitaHON0rNMx/2a3N+sk/uS4W7rjw54fiv3frvyG/VH7oYnQXEdzHs19NAfSXEhzIsMoGj6qGLWECpwOe6otDAz1cediCk7vuNOk71+SkYHCiIhGWQDjsuOE+1YesE6B7HXBvINU/iX3zBmxAL7Z9U10cO6A3OJcEeheUqq+bfMYpiHQMU2WPzrG6VinY74+4k9O/2b9MaP1DLFP5x6dg5rG6e23xZxHc9/wp9tysWdGabAAZNQWFx9r9H04QOyf3na7SZNC5Nm/xt7eMLRrWMHVT059ImKS2jq2xQsdXqjz/5l1lARg3pmz4tZA30BYD6kExrmkc1h5eWWDxsMw6sovl34RxzYd4x/3/lgc8w3lxQ4vinNOVigT56CmJX2c3i5d7NLc5+ytPkXxGe2DBSCj1gT2cEPbvh4VSSGpcdkaUf/vaupV7I/eL6wYH/b6EEYGNcf91WQBpCxkakVHeFp54o0ub4j9peeWitdnGG3gSuoVLDu/TOwv6LIAHpbS+d5Q6Fyjc04PeuIcvJZ2DZoAzW3k+iXa9vMUcx/DKBMWgIza03OK/92kkGUXkZ9d1HQCMKRhAvD7C9+L2+G+w+Fn41ev/6WWc2R5RFlZxTiIsc3HYlCzQSguK8Ybh99AfnF+g8bGMOoCHcMLDi8Qx/Rg78EY03yMQl6Xzjk69yqfi+oMzWk0txUVSEkfPSdz0gejfFgAMmoPBUBTAVQrB1PIUvKx88cwlJSUKu39ykpKkHeh4QWgyeKwL3qfsEA8Gfxkg8Zg1rGjuM0tdwMTVP3/ne7vwNHMERGZERqf6cgwi88sFseyk5kT3un2jkI7XDwV/JQ4B/dG7UV4WjjUFZrLaE6juc3a0VTMdZz0wTQFWnOUHTp0CKNHj4a7u7uYRDZv3vzA/zlw4AA6dOgAExMTtGjRAqtWrWqSsTL1h3oDj5wXDCMTA8SGZ+BoeVN0ZVBw8xZKc3Kg38AC0HKLwzDfYfW2/skxL48DzD17psrjdqZ2eL/H+2L/z6t/IjQhtEGvzzCqho7d1ddWi/33e74PW1Nbhb6+n60fhvkMU3sr4JF1N8ScRnPbiLnBatMHndF+tEYA5uTkoF27dli6dGmdnn/79m2MHDkS/fv3x/nz5/HCCy/giSeewK5du5Q+VqbhnUIGzW4N6AFhB2Nx6ZByKv7L3a6mwcHQM6hfMDpZGsjiQJaHp4OfbvAYzMrjAPPDLqG0sLDK33p79sZE/4lin/qgFpVqftFbRregY/bDE1LJFur00cujl1Le56l2khVwT9QetbQC0hx26WCsmNMGP9ZazHEM01RojQAcPnw4PvzwQ4wfP75Oz//+++/h6+uLL7/8EoGBgXj22WcxadIkLF68WOljZRqOX3sndB0jWdUO/3UdseHpCn+PvHPnGlz/r8L65zNMWCAairGPDwzs7VFWUID8y5fv+/uLHV+EnYkdbmbcxJ9X/mzw+zCMKvjjyh+4lXkL9qb29cqQry/NbZtjqM9Qsf/DxR+gTsSEp4s5jOg21g++7ZxUPSRGx9AaAVhfjh8/jkGDpL6rcoYOHSoeZ9SbjsO84d/ZBaWlZfj3hzBRMFU5ArB+8X9kYSBLA1kcyPLQGCiMwaxDiDSes3fjAOXYmNgIEUgsu7AMCTkJjXo/hmkq6FhdfmG52KdjmI5lZSKPBdwduRvX0yXBpWpoztr5Q5iYw2gu6zDUW9VDYnQQnRWACQkJcHFxqfIY3ZfJZMjLy6v2fwoKCsTfK29M00PiaMCMVnDxtUZBbjG2Lb2A/BzFuEGLEpNQGBlJb1JRjqWuyC0MZHEgy0NjMe9wfyJIZca2GIsQ5xDkFefhs9OfNfr9GKYp+PTUp+KYpYLPisr6rY0Wdi0wxGeI2sQCUsbvtu8uiLmL5rABM1spNPmFYeqKzgrAhrBo0SLY2NhUbF5eXqoeks5iaCwFTFvZmyIzKQ//fh+GkuLGZwbnnj4tbk0DA2FgXfcirJTJSBYGYf0Lbpz1T455JQtgWdn9rfDkXUIM9AzEex+OOayQ92UYZXEo5pCwktMxW99Wb41Bfk7SeULnqqqgOYq8FpnJeWLuojnM0KjhRa8ZpjHorAB0dXVFYmJilcfovrW1NczMzKr9nwULFiAzM7Nii46ObqLRMtVBLZJGPhMMI1MDxN3IwIHV4dUKpfqQe+qU9NpdutTr/9ZeWytuqd0bWRwUAfUE1jMxkdrS3b5d7XMC7AMwPXC62F90ahEKSgoU8t4Mo4yaf4tOLhL7jwQ+gpZ2LZvsvf3t/CtaMa4LXwdVQHPTgT+vibmK5iyau7jNG6NKdFYAdu/eHXv37q3y2O7du8XjNUHlYkggVt4Y1UJZc9QsnTwo147F4+yuyCYXgLlFudh6a6vYn9ZqGhSFnrExzIKDpfc4U7UcTGXmtZ8HZzNnRGdF4+ewnxX2/gyj6HZvMdkxcDZ3xtz2c5v8/R8KeEjcbrm5RZyzTQ3NTdeOJ4i5iuYszvhlVI3WCMDs7GxRzoU2eZkX2o+Kiqqw3s2cObPi+U8//TQiIiLw2muv4dq1a1i2bBnWrVuHF1+UAusZzcE7yAG9pkjWhBObI3AjtKplt17xf3fuSPF/naT4u7qwLWIbsouy4WPtg25u3aBIKvoCn5USU6qD+qe+1uU1sU8CkIQgw6gT0bK7FyevdX5NHLNNTTf3bvC29hbn6vbb25v0vW+cThRzE9F7aksxZzGMqtEaARgaGoqQkBCxES+99JLYf+edd8T9+Pj4CjFIUAmY7du3C6sf1Q+kcjA//fSTyARmNI/g/p5iI/asuiLcLE0R/0dunb/C/xL7UwKmKDymSZ6Icm9B6HsZ4j0EXd26orC0EJ+d4oQQRr349PSn4tikCyQ6VlUBnZtTWk4R+39d+6vR4SJ1Je5GOvb8KvX4DR7gKfr8Mow6oDUCsF+/fuKEvneTd/egW+r8ce//nDt3TmT33rp1C7NmzVLR6BlF0HOyP3zbOaK0uAw7ll9EekKO0t2/55LO4Ub6DZgamColo1GUotHTQ1FkFIqTk2t8HmUR/q/L/2CoZ4gDMQdEsD3DqAMHow/iYMxBcWwu6LpApRmvlDlP5yqVgzmffLfPtrJIi8/BjuVhYk6iGqY9J9W/sxDDKAutEYAMo6+vh8GPB1WUh/nn2wvIlVXtoqFoAUiWBGKk30il1DMjS6RJS8m9nVuLG5igwtOPtH5E7H9y6hNOCGFUDh2DdCwSM1rPaHBrREVB5+gIvxFif821NUp9r5zMgirlXgY91lrMUQyjLrAAZLQKI2MD0TPY2skMWan52L70AooKSpQS/5eSl4LdUbvF/tSAqVAW8r7AeQ9wAxNPt3u6IiHk18u/Km1MDFMXVl1aJSV+mDk3uji6opCfq1QShs5hZUBzzvalF8UcRHOR6GNuzOVeGPWCBSCjdZhZGWP0s+1gamGEpMgs7PrxEkpKSutk/atP/N+G6xtQXFqMdk7tEOgQCGVh1lESpDknTj7wuRRc/3Knl8X+jxd/RFx2nNLGxTC1QcfeT2E/if1XOr+iksSP6mjt0BrBTsHi3N14Y6PCX5/mGppzkqOyxBxEcxHNSQyjbrAAZLQSWxdzjJgXDAMjfUReSsWB36/VGvRdX/cvLR7rr69XuvWPsOjRQ1gmC8LDa40DlDPcdzg6uXRCfkk+vgj9QqljY5ia+Pz05+IYpGORemOrE9MCpHJNdA7TuazQWn+/XxNzjqGRvqj1R3MRw6gjLAAZrcWtuY1UI1BfD9dOJFSUYVCEAKTA9sTcRNiZ2FU0m1cWhnZ2oig0kVOHXtUiIaTr/yo6hByNParU8THMvRyJPVLR8YOORXVrdUat4ejcpb7ElKCiKE5sviXmGppzhsxpA1c/5fY5ZpjGwAKQ0Wp8gx3Rb3pARSHWC/vur5FXlJgo9f/V169z/N/acKnzxwT/CTA2UL57x6JnT3Gbc/RonTsfPBz4sNj/4MQHovcqwzQFdKx9eOJDsU/HIB2L6oaJgQnG+4+v0sWnsVzYG42zu6RSY/0fCRBzD8OoMywAGa2ndU93dB0jZR8eWX/jvkLRuafqV/+PEiyOxx8XfX8ntZyEpkAuALOPHqtz/bJn2j8DF3MXxGbHYsXFFUoeIcNI/HDhB3HMuVq44tn2z0JdmdxysrilczkmK6bRhZ5pbiG6jvVDYA93hYyRYZQJC0BGJ+g43Btt+3oAZcCelVcQdTm1we5feeB4D/ce8LRqmqKuZiHtoWdujpKUFBRcv16n/6Gge3K/ybMxb6bfVPIoGV2HamLKs8+pLqW5kfrGv9G5S+cw0ZhkEJpLqPg8QUWeOw7zVtgYGUaZsABkdAKKQeo1tSVadHRGaUkZ/v0hDPG3Mu8RgJ0f+DpFpUXYdGOT2J/YciKaCn1jY5h37iT2c47UPaZvQLMBGOA1AMVlxXj/xPsoLas9G5phGgodW+8ff18cawObDUT/Zv2h7kz0l87hTTc3iXO7vsTfzMC/34eJOYXmll5T/NUu3pFhaoIFIKMzUBHWQbNbo1mQPYoLS0WNwIQLkXfj/8rLrdTGoehDSM1PhYOpA/p59UNTYlnPOEA51H3B3NBcdC1RRtkLhiE23NggumvQsfZGlzegCfT36g97U3tRD7C+3XNSYrKwbelFFBeVijmF5hYu9MxoEiwAGZ3CwFAfw55sKzKEqUL/9p9vINfMqc7xf+tvSKVfxrUYByN9IzQl8jjA3NBQlObn1/n/RCxWiBSL9dWZr5RW/JbRXeiYWnxmsdh/LuQ5ccxpAkYGRuJcJv6+/ned/y8jMRdbvzmPwrxiMZcMe6qtmFsYRpPgI5bROYxMDER9LgdPS+QX6uN8u+eg17nPA/+PAtuPxR6r4jpqSoz9/GDo6oqywkLkhj64K0hlHmr1EALtA5FVmIXPTn2mtDEyugkdU3RsUZFlOtY0Cfm5TOWS6lI4PTs9X4i/vKwiOHpZirmEu3wwmggLQEYnMTE3wui5QTDPT0a+qQOOpAc/sG8wxf6VoQxd3brCy9oLTQ3FFln07NEgN7ChviHe7f4u9PX08e+df0V9QKYGCrKA3LSqGz3GVMt/d/4TxxQdW+90fwcG+polhppZN0NX167i3KZYwNqgOWLL1+eRlZYPG2czjH6uvZhLGEYTMVT1ABhGZdwIQ7tzS3Cu48vIzLTFlq/PYdxLITCzvL+uH3ULkCd/NFXpl5riADM3bKy3ACSCHIPwWJvHRHuuD45/gBDnEDia6WitsuJCICEMiDkFJF0BMmMBWax0W1iD2DO2Amw8AGsP6da5NeDZBXBtCxga66zrl+pMEo+3eRxBDkHQROicPplwUsTIPhX8lLhgupe8bBJ/54T719LeBGPmt4e5tW7+7ox2wAKQ0Vmy9++HWUEaejtexlG9/kiLyxGunXEvhtx3VX845jCS8pJEwPhAr4EqG7N59+5SW7jr11GUlAQjZ+d6/f/cdnNFsPv19OtCBH7d/2vdyFosLQViQ4Fr24Hok0DcOaC47nGUAhKGydekrTKGpoB7CODVFWg1EvDoJJKKtB2qR/ne8feQUZCBlnYtxbGlqVC2PHUGScpNEl1M7k3wys8pEnMDzRHmNsYY+0IIrB3MVDZehlEELAAZnYQWr6x9+8W+66CuGNs2BJu/OouU6Gz88+0FcXVvbGpYJcORGNN8jAgcVxWiLVxQEPIvXULOsWOwHScFsNcV6lryca+PMW37NOyL3od/Iv4Rn0lrRV/MaeDKZuDKFsnCVxkze8CrC+DWHrD1KrfseQLW7oDhPYs7dVKRxQGZMdLrZEQD8eclMZmXDkQdl7ajX0uv03os0Hoc4NlZa8Xg1ltbcSD6gLCW0TGlyvOisdB5QefBr1d+xYbrG6oIwML8Ymz77oKYG8ysjIT4s3VW3/qGDFNXWAAyOknhrVsoio6GnpGRcKvqW1gI0bf5q3NIvC3D9qUXMerZdiJhhPqFHo49rLLkj+qygYUAPFp/AUgE2AdgXrt5WHJuCT45+Qm6uHbRmKzNOpGXAZz/Ezi1Aki/U9WFGzAM8OsvWescmgtrap0wtgAc/aWtMtSVJfWWJAQj9gPhOyWBeGKZtNn5AF2eBEIeAUy1py8snROfnPqkouMMHVOaDtX1JAF4KPaQ+Hx0ThQVlAjxR3OCiYUhxswPgb2bhaqHyjAKQTsvTRnmAWTtl6x/5t26CfFHOHpalVv+DBB3I0PUCaQFgKx/VOS2s2tn+Nj4qIEALE8EOXYMZWTlagCz28xGsFMwsoqy8M7Rd+rcXk6tSbkBbH8F+Ko1sOt/kvgj0dd2CjBtDfDqTWDiT0DIdMCxRd3FX23Qa9Br0WvSa9N7TFstvSe9N42BxvJloDS2FM3vxkLnwttH30Z2UbY4hmYFzYI24Gvji04uncTno1hAufiLv5kp5oQxz7eHo6elqofJMAqDBSCjk2SXu3+tBlTtVuDsbY3Rz7eHkakBYq9n4J/vzmPT1S3ib1MDpkIdMG9f3hYuNRX5V6426DXIbfdRz49gamAqeqH+cfUPaLTwWz8L+K4TcPpHoCgHcAoERn8DvHIdmPgj0GoEYGSq/LHQe1AcIL0nvfeor6Wx0JhobN91lMZKY9ZQ/rz6J07EnxDHDh1D1SVMaCpTW0nnOJ3zdO7ThSDNBTQn0NzAMNoEC0BG5yhOS0Pe+fNi37Lf/d08XP1sxNU+TfzxNzLR5dxEuBi7iUBxdUDP2BiWvXqJ/axduxr8OmTNfKnTS2L/q9CvcD5J+k40BorD2/IMsLQLcJkytPWAgBHAzC3AvONAx1mAsQpjtei9O82WxkJjorHRGGmsS7sCW56VPoMGQccIHSvEy51eVguLuCKhBC8617ucmyDOfZoDaC6gOYFhtA0WgIzOkX3wkIjdMmkdCCM3t2qfIxeBJYaF8JD5Y/yN54Ai9TldrIcPE7eynTsb5b6dFjANQ32Giv6trxx8Ben56VB78jOBnf8Dvu0AnPsDoP7GJK7mHgUeWgP49VOMe1dR0FhoTDS2p49IYy0rAc79Ln0G+iz0mdSctPw0cYzQsTLMZ5jaWMQVSrE+xl1/Du4yf3Hus/hjtBn1WdEYponI3rdP3Fr1q71ZfbZdMra2WopCg3wgzkLEA1FGoDpg2bcv9MzMRCJL/qXLDX4dKgHzXo/34GPtg8TcRLxx+A2UlJZALSGhe3E98F1n4MRSoKQQ8OkNPL5HElcuGlCDzrWNNNbHd0tjp89An4U+U9jf0mdUQ+iYWHB4gThG6FhZ2GOh1pUPkmf76sVbiHN+S6vvkGPPbRMZ7YUFIKNTlBYUILu8iLJl/9oF4NrwtUi0uoOkfqHCFUTxQP8sOY+C3CKoGn1zc1j26yv2ZTv/bdRrWRhZ4Kt+X4mYrmNxx7AibAXUjuRw4NfRwMYngOxEwKEF8MgG4NF/AK/O0Dio/AyNnT4DfRb6TBselz4jfVY1Y8XFFeLYMDM0w+J+i8Uxo03QOU11/ugcp4SPxH6nkWQVKeYAhtFWWAAyOkXuqVMoy82FobMzTINa1/i8nKIcUSOPGNtriKj9ZWJuiIQImWgFlZ+tehFoPWy4uM36t3FuYMLfzl+08SKWn18uFnu16dix/2NgeU/gzmGp6PKAt4G5x4AWg9TL1VtfaOz0GeizDHhL+mz0Gemz7l8ElKj+GCOo//XyC8vF/tvd3kYLuxbQJuhcpnNalHoxN8TYF0MwtueQilqHuUW5qh4iwygFFoCMznX/kCd/6NVSoHd7xHYhAsndRX1CXXysRZs4U0sjJEdlYfPisw/sHaxsLPv0FtnARXFxyL94sdGvN7r5aFHnkHqivn7odUTKIqFSEi8DPw0EDn4KlBYBLYcDz5wE+rwCGJpAa6DP0udV6bO1HCZ91oOfAD8OABKvqHRodAy8fvh1cUxQuzQ6RrQJOoc3fXVWnNNU5JnOccr2pX7fdO7THLAtYpuqh8kwSoEFIKMzUM08efcPy3vKv1R5XlkZ/gr/S+xToLs81onqBI5/qYNoBZUam4NNX55Fdno924kpEH0zM1iVu7Fl/+5UyGsu6LoAbRzaiPZe8/bME4H/TQ7FIB5ZDKzoByRcBMzsgEm/AA//JRVW1lbosz28Vvqs9Jnps6/oK30XKojLpN9+7p654ligY+KNLm9Am6Bzl85heXu3cS91EOc4oa+njykBU8Q+uYG1ok4mw9wDC0BGp9y/xQkJ0LeygkW3bjU+71zSOdxIvyFi4sa0qNomzd7dQohASzsT0RR+w+dnxK3Ks4F37WpwUejKmBiY4NuB38LD0gNRWVF4bt9zyK9vz9zGQIWTfxkG7FkoJUiQ1W/eSaCN6juwNBn0WeeVWwPpO6Dvgr6T9KazyNJvTr99dFa0OBbomKBjQ1uofO5a2ptg/Msd7uvwQa3haA6gvtnnkzWsRBLD1AEWgIzOkLFxo7i1HjEC+qY1FwX+65pk/RvhNwLWxvcXf7V1Mcf4VzrAxtkM2WkF2PjFGSRHZ0EVWPTuLTqZFMfHI+/8BYW8pqOZI5YNXCY++8XkiyL7s0kyg6lf7/d9gJhTgIk1MG65lDFr5QKdgz7zQ38BY5dJ3wV9J9/3lr4jJUO/NWWD029Px8CyQcvEMaEtkLuXzlk6d8W5/HKHanv72pjYYLivFGe75toaFYyUYZQLC0BGJyjJykLWf7vFvu2E8TU+Lz47Hv9F/if2a6tzZu1ghgmvdISjlyXysoqw+cuziLuZoYSR146+iQksBw5QSDZwZfxs/fBN/29gpG+EPVF78OWZL6E0ivKBbS8B62YCBZmAZxeppl/7hzU7yaOx0GenFnP0XdB3Qt8NfUfbX5a+MyXxRegX2Bu1V/z2SwYsgZ+NH7QFyvLd/NVZcc7SuUvij87lmpjWapq43X1nt+gPzDDaBAtARicQBZPz82Hs5wfT4OAan/f71d9RUlaCLq5d0Nqh5ixhwtxaihtya2GDwvwS/PPNedwJS1FdNvBOxbiB5XRy7YQPe34o9n+/8jt+ufQLFA61RKNEj9Cfpfu9XgRm7wBsmyn+vTQV+i7oO+n5gnT/9E/AT4OU0k6OfmN5W8CPen2Eji4doS3Qubl1yXlxrtI5S+cuncO1QXMA9QCn4td0DjCMNsECkNEJMjduqrD+1VTANrMgE39f/1vsz24zu06va2JmKPqEerd1QHFRKXYsD8PVY/FoSix69RRxjcVJScg7e1ahr01u8Bc6SMJj8ZnFihWBV/8BVvQHEi8B5o5STbxBCwEDI8W9h7ZA38ng94DpG6TvKjFM+u6uKi5D9eewn8VvTLzY8cUK96c2cPVYnDg3S4pK4dPWQXT4oHO3LswOkuYCmhtkhTIlj5Rhmg4WgIzWUxBxG3nnzgEGBrAeUzWpozLrr69HXnGeqInX071nnV/fyNgAw59ui4CurigrLcO+364idMedJssc1Dc2htXAgWJftkNxbmA5j7d9HPPazRP7JBB+CvupcS9I8YR73gPWPgIUZgHePSU3J9XEY2rHf5DUTq5ZD+m7Wztd+i4bGaNJv+nXZ78W+/Paz8NjbR6DNkDnYOiO29j32zVxbtI5OuzptjA0Nqjza/Ty6IUWti2QW5yLdeHrlDpehmlKWAAyWk/mJsn6Z9mrF4ycnat9TkFJAf68+mfFFX9921wZGOhj4KxAdBjqLe6f3BqBQ39dR2lp04hA6xGStSZz+3aU5uUp/PXntp8rhAHxzdlvGi4Cc1KBPyYCR76S7nd7Bpi5BbByVeBotRxrN+DRrUA36fcQ3yV9p7kNK9nz48UfxW9KPNP+GcxtNxfaAJ17h9Zcx8mtt8X9DsO8xTlK52p9oLlA7hGgOaKQMrMZRgtgAchoNWUlJcjcImVO2oyvOflj261tSMlLgYu5C4b5SqVV6gstFN3HN0fvqS0BPeDSwVjs/CEMxYXKz6C16NkTRh4eKM3MROY25RSuJWHwbPtnxT4Jhh8u/FA/K2c81bXrB0TsB4zMgYk/A8M+ZpdvQ6DvbNgi6Tuk75K+0x/6St9xHaHf7vsL32PJuSXi/nMhz+Hpdk9DG6Bzjs69S4dixblI52T3cc0b3L94uM9wOJs7izmCC0Mz2gILQEaryTl2TMTGGdjY1Fj8ubSsFKsurxL7M1rPENmPjSG4vyeGPtEGBob6uH0hBZsXn1N61xA9AwPYPfyw2E//40+luZ+faveUEArEd+e/w0cnP0JxafGD/5HKl/wyFMiMAuz9gCf2AG0nKWWMOgV9h/Rd2vlK3y19x3UoFUO/2YcnPsTS80vF/edDnseTwU9CG6Bzjc45OvfoHBw2p404JxuDkYERZgTOEPs0V9CcwTCaDgtARjdq/40eLWLlquNA9AHckd2BlZGVaIWmCFp0dMaY+e1Eb1HqMfr3p6Gi44AysZ04AXqmpigID0deaKjS3oeEwqudXoUe9ESXhGf3Povswuzqn0xZydTXlsqXUE9Vv/7AnH2AS5DSxqdz0Hf55H7pu6XvmL7rA59I33010G9Fv9m66+vEb0i/5ZzgOdAG6Byjc0309bUwxJj57dG8Q/VhH/WFWuFZGlniduZtHIw+qJDXZBhVwgKQ0VpKMjKQvWfvA2v/ya1/kwMmw9LYUmHv7+5vh0mvd4KNkxmyUvNF54Hoa8prrWZgawub0VKv1rQ/pHhGZTEzaCYW918MM0MzHI07ihn/zkBcdlzVJxVkA+tnSn1t5fF+0/+W2pwxioW+U/pu5XGBBxYB6x8FCqtedNBvRL8V/Wb029FvSL+lNhB9NU2cY3Su0Tk36bVOcPe3Vdjr09xAc0TlOYNhNBkWgIzWkvnPNpQVFcGkVSuYtm5dY9s32sjtOz1wusLHQJ0GJr7eUaoVmFeMbUsu4MrRe4SSArF75BFxm7VnD4oSlFu4dmCzgVg5bCWczJxwM+MmHt7+MM4nlbfMyoyR2pdRqRdyqY/5rjzer26lN5gGQN8txQXSd03f+dWt5W73WPFn+m3oN6Lfin4z+u3oN9QG6Jza9u0FcY7RuUbnHJ17iuaRwEdgqG+Is0ln7x7rDKOhsABktBISfmkrV4p928nVx5pRnNx3574T+6ObjxZB3srAzNIYY+eHwL+zi8hM3P/7NRxZfwOlJYqPIzINaAnzzp2BkhKk/yW1tFMmQQ5BWD1yNVratURqfipm75yNn4++j9IfB0i16iycgFnbgA5S/BTTBNB3Td851QtMCBO/xU9H38esnbPEbxRgFyB+M/rtNB06h46suyHOKTq3WnZxEecanXPKgOYI6hFMyOcOhtFUtEoALl26FD4+PjA1NUXXrl1x6tSpGp+7atUqkRFWeaP/Y7SDzG3bURQXBwMHB9hOrD6u73j8cZxKOCWsf08FP6XU8RgY6WPwY63ReZSvuH9hbzS2Lb2I/JwipVkBM9atR2lBAZSNq4Urfhv+G4b5DBMdE76+uR5PW5YixaWVFO/XrJvSx8DcA33nc/aJ3+ApyxJ8c3O96HBD2ay/Dv9V/GaaDp072767gAv7osX9LqN9MWh2a3GuKROKgaU542TCSRyPO67U92IYZaI1AnDt2rV46aWX8O677+Ls2bNo164dhg4diqSkpBr/x9raGvHx8RVbZGRkk46ZUV7pl9QVK8S+/axHoV+NsCfr35KzSyp6/rpbuit9XHSR0WWUL4Y92QaGxvqIvpKGDZ+dQXqCYpNDrAYOgKGrK0rS0iD7V/GFoavDwtAcn8EZ7yWnwrS0FMfNzDDRzgRHc6TFmWl6jubGiN/ghJkZzEpL8X5yKj7Vcxa/laZD5wwle0RfTRfnEp1TnUf6NrjMS33wsPTAlIApYp/mkKYq+M4wikZrBOBXX32FOXPmYPbs2WjdujW+//57mJub45dfam5dRZOFq6trxebi4tKkY2aUQ9buPSi8fRv61tawe+ihap+zJ2oPLqdeFoHwT7R9oknHR1mJE17tCEt7E2Qk5uLvT0IV2kNYz9AQdtOmKb0kTAXFhcDmedDb/yEmZOdgrWM/tLT1R1pBOp7e8zTeOvIW0vPTlTsGpgL6ruk7p++efgP6Lf5y7Ifx2TnQ2/eB+K3Eb6ah0LlC50xmUp44hya+1lFhmb51ZU7bOWLuuJR6CXujpEQzhtE0tEIAFhYW4syZMxg06G4rKX19fXH/+PGaTfTZ2dnw9vaGl5cXxo4di8uXLzfRiBllQWInZcUPYt/+kekwsLSstgbat+e+FfszW8+Eg5lDk4/TycsKk9/oDLfmNqI5/fZlF3Fq223RrkoR2E6ZDD1jY+RfuiS1wVMW1H3i9/HAhdVUjBAY+SX8Ri/F6lFr8HArqS7hlltbMGbzGGy5uYWtJUqEvtvNNzdL3/WtLaLEC/0G9FvQb4IRX0i/Ef1W9Js1sHOIqqBz49Q/Edi+9KI4ZyjZg84hR0+rJh8LzRlUM5SguaSkka34GEYVaIUATElJQUlJyX0WPLqfUEMmZEBAgLAObtmyBX/88QdKS0vRo0cPxMTE1Pg+BQUFkMlkVTZGvcg5fBgFV65Cz9wcdjOqTzz459Y/opaXjYkNHg16FKrC3NoYY18MQZu+HkAZcHrbbWxfrpi4QEN7e9iMlYLVk79arBzhlXoL+GkQEHkEMLYCpq8DOkvWVBMDEyzougC/D/9d9FbOKMjAW0ffwuP/PY6b6TcVPxYdh77Tx3Y9hrePvi2+a/rOKS6TfgP6LQRd5gAPr5N+K/rNfh4s/YYaAJ0TdJF0evsdcZ/OmbEvhIhzSFXMCpol5pCIzAj8E/GPysbBMDotABtC9+7dMXPmTLRv3x59+/bFxo0b4eTkhB9+kKxH1bFo0SLY2NhUbGQ5ZNSLlB+k2D+7qVNhaGdXbc/fZReWVbhxrGgxVCHUqaDvQwEY+GigCF6PDEvF+kWnkRJTQ2HleuA4bx70TEyQGxqKnEOHoFAijwE/DQTSbgE2XsDju4AWdy3wcto7t8faUWvxYscXYWpgitMJpzHxn4l488ibiM2WypMwDYe+Q/ou6TsNTQwV3/FLHV8S3zl99/fhP0j6reg3S70p/Yb0W6oxKTFZ4pyIvJQqzhHq50vnDJ07qoTmjifaSBc8y84v4x7BjMahFQLQ0dERBgYGSExMrPI43afYvrpgZGSEkJAQ3LxZs3ViwYIFyKReq+VbdDQHuKsTuadPI+/MGegZGcF+1qxqn7MufB0SchJEOQdK/lAXWnV3w8RXO8LKwRSylHxs+DRU1DZrjOXOyM0N9jOkjOCkL78SyTEK4cJa4LexQF464N4BeGJvrZ09KGPysTaPYdPYTaLuHLXR2nprK0ZtGoVFJxeJ/qpM/aDvjL47+g7pu6TvlL7bzeM2Y3ab2bW3M6Tfin4z+u3oN6Tfkn5TNYOOfToHNnx6RpwTdG7QOdKqmxvUhWmtpom5JD4nXswtDKNJaIUANDY2RseOHbF3791gXHLp0n2y9NUFciGHhYXBza3mycXExERkDlfeGPWz/tlMnAAjl/uDwjPyM7DiovScue3mwtRQvcr+ODWzwpQFndEsyB7FRaWittneVVdRmF+HXrs14DBnjkiGKbh+HZn/NNJNRWJ030fApicBsnYEjgFmbQes6pY85Wnlia/7f43VI1ajq1tXEYu5+tpqjNg4Ah+f/BhRsqjGjU8HiJRF4qMTH4nvjL47+g7pu6TvlL5bylCtE/Sb0W8XOFr6Lek33f+x9BurAXTM07FP5wCdC3ROTPlfZ3GOqBM0hzzd7mmxT3NLZkGmqofEMHVGr0xLorKpDMyjjz4qXLhdunTB119/jXXr1uHatWsiFpDcvR4eHsKNS7z//vvo1q0bWrRogYyMDHz++efYvHmzSCahLOK6QDGA5AomayCLQdWSffQooh9/AjAwQPOd/8K4Gvf8u8fexcYbG9HCtgXWjV5Xu5VExcHuZ/+LxMktEWI9tnM1x9A5beDg0bA2dak//4ykz7+Aobsbmv/7L/RNymPC6kNRPrDlGeDS39L9Xi8CA96hbCs0lBPxJ/DNmW9EJiVBSQv9vPqJxJyOLh2bpKSHJkBTNLl3f7/yu+hbXUYBoxQH59AG8zvORze3RtRZpH7Be98Djn4t3W8zCRi7FDBS3cVRamw2dv14CekJudDT10PXMb7oMMRb7KsjRaVFmPLPFNFhhXqJL+yxUNVDYuqAjNdvaE1fpqlTpyI5ORnvvPOOSPyg2L6dO3dWJIZERUWJzGA56enpomwMPdfOzk5YEI8dO1Zn8ceoD6WFhUh8/wOxbzf94WrF35nEM0L8Ee92f1dtxR9BC13HYT4iQ/i/ny6LhZDKXvSe2hKBPd3qLYzspk9H2u9/oDguHumr18BhdvXu8RrJSQH+ehiIPgnoGwKjFgMdGt8/loRL15FdRUFdEjeHYg5hf/R+sZFIH9t8LEb6jYSTuRN0keTcZGyP2C4yeklcyOnj2UeI5C6uXRovkmlOHPweYO8HbH9JEvjUxm/an4CFI5rc5XskTnT2IKufha0JhjwepNB+vsqA5pJ3ur+Dmf/OxIYbG0RXIbqAYRh1R2ssgKqAryDUg5Tvv0fy19/A0MkJfv/uuK/0S1FJESb9M0lk601qOUkIQE0hL6sQe1ZeQdQVqWSHX4gT+j/SCqYW9ROwGRs2Iv7NN2FgY4Pmu/+DQV2P16SrwOqpQEYkYGoDTPkd8OsLZUC/z59X/hQxbfkl+eIxfT19dHfvjjF+Y9DXqy8sjCygzWQXZgshvDViq+gyQbF9BCV3UAuy6a2nw8/GTzlvHnEAWDsTIDemrbeUMezcCk1BfnYR9v9xDRHnk8V9cvkOmtUaZlaqy/KtLwuPLRQCsLlNc6wfvR5GBup7kcnw+k2wAGwELABVT2FMDCJGjkJZQQHcv/gCNqNG3vecHy/+iCXnlsDe1B5bx20VpRs0CXIJn9sTJVzCpSVlwjJCLa88A+zq/holJbg9bhwKbtwUcYHOL7/04H+6sQf4ezZQIAPsfCRB4BQAZUNxVLvu7BLles4nn69iaens2lm4ift59oObpfokAzSG+Ox4HIg5INy71JqQ4vrktHdqjzEtxmCI95CmOW6Tw4HVU4D0O4CJNTB5ZbXZ3Yok5loa9qy6ipyMAugb6KHrWD+EDGqmti7f2o5bqsGYlp+G+R3mN3mBeaZ+yHj9ZgHYGPgAUj3Rc+che/9+mHftimarVt7nEouWRWP81vGi/Mui3oswym8UNJWkSBl2/3JFdA+BHtBhSDN0GeVX596nWfv2I2bePEp5h+/f62EaUIuYO7kC2Pk6qU/Au6dk+bNwUEnSA1kESRDSfmXI0tLJtRM6uXQSLjdNcRWTa5di+igsITQhFLcyq9bi87b2xlCfocLiR/tNTk4qsPYRIOoYxSMAwz4Fuj6p8LcpKSrFyX8icG53lKiDaetiLly+6pboUR/oouV/R/4nai9uGrMJXtZcKkxdkfH6zQKwMfABpFqy9u1DzLxnhKDx27wJJs2bV/k7HdrUDutY3DERb7Zi8AqNTywoKijBkfU3RKwUQYkhVBeNOos8CPo+Yp59Dtl798KkVSv4rlsruoVUoaQI2PkGcPon6X77R6SYP0PVu+KoeDdZyWgjy6DcPSqnmVUzBDkEIcA+AIH2geJWFV1eKpOal4rwtHBcTbsqbqn9YFRW1WxncnOTpU9YNr36wdfGFyqnuADY9iJw/k/pPhX4HvYJoCC3ZnJ0FvauuoLUWKkPduve7ug1yR9GJgbQZOgcm7N7Dk7Gn0RP955YPmi5xs852oqM128WgI2BDyDVUZqXJ1y/RXFxNbo05VfjxvrGogZdM+tm0BYiziXjwOpryMsqgr6+HjqP8kGHod7QN6jdGlickoKI0WNQkp4Oh6efgvMLL9z9I7UGWzcTuHNY5ORi0EKg53xqmg11g0r6CAtauSXtWtq1iuzYypDb38vKS2wkEKkUDdVtczB1gKOZI6xNrIUAawgkQGUFMlGTLzU/FUm5SYjJihECLzorWmzkDrwXer8AuwBhtZRbL21N1TDRgZaGo98AeyirtQzw7QNM/hUwt2/wS5aUlOLszkiEbr+D0tIymFkZod/0VvBrrxnW27pAluoJWyagsLQQH/f6WCSFMOqHjNdvFoCNgQ8g1ZH46WdIW7lSKm2ybRv0zc3vm4SpNENucS6eC3kOTwYr3oWlanJlhTi4JlyIQcLZ2woDH20Ne/faEyVkO3chloSfvj581qyGWbt2QOIVYM00KdnD2BKYsAJodX88pbqSVZiFi8kXKyxtJAjpGKhOFFbGUM9QxNaZG5nDzNAM5obSrYF+VUsU9XrNK84Tx5O4LcoVMV/FZbXXaKTSNuTGbWXfqsIyGewUrPIONPXi2nZg45NAYbYUC/rQX4BzYL1fJi0uB3t/vYKkyKyKhKZ+DwdoVKJHXfnhwg/47vx34niiklMqceUztSLj9ZsFYGPgA0g1ZO3di5hnnhX7nsuWwWpA/yp/p3i/R3Y8IkQAWVh+GvLTfQu6tkCn7/VTiTi89joKcouhb6iHTsMla2BtrbJiX3kVsm3bYOzrC99PnoL+9nnSAk/Zn7TAu2h+OSQSaSQCK1vkyEJHFjvaZIWK6eVNAlJuUSQLY2WLIy38JC41nsTLwJqH7l4gTPwJCBhep38tKS7FmZ2ROPPvHZHEZGJuiD7TWsK/s4vWukcpkeeJ/54Q1mkS/b+P+P1uT2ZGLZDx+s0CsDHwAdT0FEZH4/aEiSjNyoL9o4/CZcEb9z3nwxMfYm34WtiZ2OHvMX8Ll5+2k51eIFzC1EuYICsglYtx9as+c7QkM1O4gouTkmDXMgeuHTIBn97AlN8a5eLTJKg8ELluyZJXYd0rkm7vjS8kt62wDhrdtRLKhZ/OlPug5JD1j94NEej/JtD75VqLgSdEZIryLmT9I3zaOgiXL2WyazuJOYmY/M9kpBekY1rANLzZ7U1VD4mphIzXbxaAjYEPoKaltKAAdx56CAVXrsIsJATev/0q+v5W5r87/+Hlgy+LfQrA7uXRC7oCnco3Q5NweN11ERtIa3RwP09RVsPY9J6a7/mZyP78YUT/GSHues3pDMsXflZYkD+jpdybJBQwEhi/XKoRWYnCvGKc3BqBiwdiRPggxfpRIfMWHZ211upXHUdij2Dunrli/6t+X2Gw92BVD4kpR8brNwvAxsAHUNMS/867yFi3DgZ2dvDdtBFGrq5V/k4uPor7yy7KxuNtHscLHSslOOgQVFT3yN83EH4iQdy3sDFGz8n+dxffpGvA2ulA6k3En7FHxg1T6FtainhAE39/VQ+f0QTO/gZsf1nqI+zQApj6pygaLS5CziSJTPXczELx1IBuriLD19RSNy8uFp9ZjF8u/QIrIyusHb1WhAcwqkfG6zcLwMbAB1DTkbF5M+LfWCAyUr1++hGWPXtW+Tu58GbtnIUrqVcQ4hyCX4b+AkNqW6bDRF1OxcG/rkOWnCfue7ayQ99O0bA9+DRQlANYe6J0/M+IfnspckNDYeTuDp91a2Ho2LQtwBgNJfYMsHYGIIsVcYHpfZfj0GkvxFxLF3+2djJD34daollr1ZbiUYdewY/tfEyULqIyRSuHrRQhBIxqkfH6zQKwMfAB1DTknj2HqMceQ1l+PhyfexZOzzxzX8D1i/tfFN0UKC7r79F/w9WiqnVQVykuKsHZXVE4u/MOSorLoI8ihFhsQYegBBhPXSH6vRanpyNy2kMojIyEaXAwvH9dBX0zXqCYOpCdjMK1T+LMFTeczxmLUhjBwFAPHYZRIlIzGBppZ/JVQ7q9TN42WcSbUq3Hxf0W6/wFqqqR8fqNhhXAYpgmIu/SZUQ/+aQQfxZ9+8BxrhRPI4euXyjpg8QfZdl9O+BbFn+VoAW4Sy8DPNTqOzQzPiMW6DM5k/DntZdw9UKRaDNnaGcHrx++F32C8y9eRNzrb6CstGoSBMPcCx07Vy8W48/wl3E2Z5I4trxNQsWxRscci7+7UNvCJf2XiDmKCpnTnMW2F0bVsAWwEfAVhHLJv34dUTNmioxV806d4PXjivssU0vPL8X3F74XWZoUZD2w2UCVjVctCf8X2PQ0kJ+BMhMb3A5egaMnbCrcwtR2q9cUf7i3sBVu4KjZj6GsqAj2jz0G51df0amAfabuxN3MwJF1N5AclVXh7u3ZLQO+F5+CXkEmQIWtx/8ABAxT9VDVir1Re/HSgZdElvncdnMxr/08VQ9JZ5Hx+s0CsDHwAaQ8Cm7fRiSJv5QU4ZZs9ssvMLCsWuB4Xfg6fHDiA7H/dre3MSVgiopGq4ZQK6+97wPHv5Pue3QEJq0E7LxFD9aL+2MQuuM2CvNLxJ+pEwNlCxuc3ou4114Xj1GZHec3XmcRyFSQFp+DE5tv4faFFHHf2NQAnUb6imxz0ZM6PRJYPwuIOyv9Q/dngYHvAIbaX/alrvC8pR7IeP1mAdgY+ABSDoUxsYh85BEUJySInrUUk0buycrsurMLrx16TVxJP93uaTzTvmpcoE5DWb4bngASw6T7XecCg9+/r58vdRI5+U8Erh6JE12/SOe16uGGlgXnkPWFtEDZTJwAt/ffh54Bu/N0may0fJzedhvXjsdXHCuBvdzRdbQfzK3v6eRRXAjsfgc4uVy679oWmPgz4BSgkrGrI9+d+w4/XPxBeC4+6/MZhvoMVfWQdA4Zr98sABsDH0CKJy/sEmLmzUNxcjKM/fzg/ftvMHSomkW44foGvH/ifSH+JvpPxLvd32UrFUGnMtVn++8toDgfMHcAxi59YMcGKtJ7Ystdqw5Zclp65MFhzUIYF8hgNXQo3D//DPrG2teyi6kdukg4918kwg7Eio4ehG87R3Qb1xz2brW3HMS1HcDWZ4HcVMDQFBj6EdDpcbXsLd3U0LK78PhCbLyxUYhAmsMm+E9Q9bB0Chmv3ywAGwMfQIpFtus/xL3+ukj4oHp0VO7FyMWlynN+DvsZX5/9WuyT+CMXira2easXWYnAP/OB6/9K95sPAMYtB6zqnhBDXRuObbyJ+JuZ4r6BQRk8IvehWeR/sOvaHp5Lvrmv5zKjvcLv/O4ohB2MQXGhJPzc/W3RfXzzGrvLVEtWArB5LnBrn3S/5XBg9DeAVdXzWheh/tJ0IUsikHix44t4rM1jqh6WziDj9ZsFYGPgA0gx0CGY+sMKJH8tCTuLPr3h8dVXMLC0rPIcKqi68vJKcZ8KPc/vMJ8tf3T6hv0N/PsqkJcOGBhL7t4uT9XaoqvmlytD5KVU4e5LipQC/PVLCuEZexAtjG7Db/EimPj5KuGDMOpATmYBLuyJriL8nL2t0HmUL7zbODTsfKOM8pPfA3velQpHm9kBI74A2kzUeWsgnW90QUuFoonZQbOFENT5ea0JkPH6zQKwMfAB1HhK8/KQsHAhMrdsFfftZs6Ay2uvQc/QsErPVrpS3nxzs7j/SqdX8GjQoyobs1pZ/ba/BFzbJt13DZasfq5tGv3S1QrB0iK4pYSi07QO8JxSu1uZ0SwyEnNxbk8Uwo8nVLh6Gy387iXhErD5aSChPDa11Shg5FdsDQSw6tIqfHnmS7E/vsV44dnQmR7TKkLG6zcLwMbAB1DjyLtwQdScK7xzh/yNcH37LdhNm1blOXHZcXjl4CsISwmDgZ4BFvZYiHEtxkGnoVP24lqpJytZ/fSNgL6vAb1eVHgv3wohuPUGkqKl0jEoK4WHWRq6PTMQrv663eVB00m8IxMxfrfOJYuevYSLrzU6jfBRnPC7t5fwkcXAwc+A0iLJGjjsUyB4is5bAzfd2CTiAim2OdgxGF/0/ULUD2SUg4zXbxaAjYEPoIZRVliI5OXLhduX3EOGzs5w/+xTWHTrVuV5h2MOY8GRBaJ6vrWxNRb1XoQ+nn2g06TcALa9CNw5rHCrX23QNBEXnooTPx5BQs7dY93Z1QjtR7aEXwcnGBhwXXlNoKSkFBFnk3FxfzQSImQVj3u3dUCHId5wa2GjfBeksAbOBRIuSvd9+0jWQEfd7kV9KOYQFhxeAFmhTHQ1WtRrEXp79lb1sLQSGa/fLAAbAx9A9Sc//DriFyxA/pUr4r716NFwfevNKmVeqLXbsvPL8GPYj+J+G4c2+KLfF/Cw9IDOUpQPHPlKsp5QHBVlVfZ5Feg5X+FWvwcRvfUgTv8RigTbtigrb2dlYW2ENv08EdjTHRY2XPNNXeP7rh6Nw6WDscjJLBSP6Rvowb+TC0KGNIODx92Y2yaBrIFHvwEOfS5lrVP8Klmxe70EGJlCV4nNjsXLB17G5dTL4v6ctnNEmStOdlMsMl6/WQA2Bj6A6k5xWhqSlyxBxrr1wupnYGsL14XvwnpY1U4BN9Jv4L3j7+FC8gVxf1rANLza+VUY0+Kgi9DpeX0nsOt/QFqE9FiLwcCIzwF71SVjFCUlIfLjr3D9ahFi3Xuh0EQS8Pr6evAJdkRgTzc0C3IQ9xnVUVpahqjLqbhyJA53wlJF+zaCave16euB1r3UQLCn3QZ2vALc3CPdt/cDhn4MtByms27hwpJCfHb6M6wNXyvut3Nqh4XdF6KFXQtVD01rkPH6zQKwMfAB9GBKCwuR/vvvSFn+PUqzs8VjVkOGwOWtN2Hk7FzxvIKSAvxw4QesvLQSxWXFsDCyELWxhvvqcLIBuclI+N0+KN23cgOGfwoEjlGbhTH70CHEvfchYotcEOPRFzIbv4q/WdqZoFV3NwR0dYWtC5ePaeqkjvCTCaJwc3Z6QcXjVMKlbT8PNO/gDANDNXLZ0zJ0ZYsU15oVLz3m108Sgi5B0FV2ROwQCXA5RTkw1DcUZWKeDH5S9BRmGoeM128WgI2BD6DahZ9s61ak/LACRdHR4jHT1q3hsuANmHfuXOW5pxNOC6tfpCxS3O/v1R//6/o/uFrUvYadVpGdDOz/EDj7m0i4EK6xbvOA3i8DptZqmcmdsmw5UleuRLaxE+LdeiDBqyeKcHeRcvaxRssuLsLdeF/nCEZhtftuhCbi+smEisxtwtTCCAHdXIVV1sG9id289SVfBhz+EjixTAp10NMHOswE+r8FWDpBF0nIScBHJz/CgegD4r63tbe4OO7sWnUeZeqHjNdvFoCNgQ+g+ynNyUH6uvVIW7kSxUlJ4jFK8nB68UXYjB0DvUq16cjdu+TckoqJzcnMSQi/gc0G6mYdLMroPfYdcGI5UJQjPdZ6HDD4PcDOB+pOYXQ0UpYuQ+bWrSgt00eyUzsktxmJ5DIXYeAh9PT14NHSFs1DnODb3kn17kctiOu7fT5ZZPHGXs+ocPHS9+wVaI9W3VxFn2fRp1eTILcw1Q0kqyBhZAF0mwv0eFbKHNYxaJneE7UHi04uQnJesnisn1c/PB/yPPztdDtxpqHIeP1mAdgY+ACquvhnrP8b6WvXojQzs0L42c+eDbupU6p0kIjOihZJHtsjtqMMZaIV0uSWk0VhZytjK+gcBVnAie+BY98CBdJ3B/cQyf3l3QOaRkFEBFK+WwrZjh3ifqGRFdLaj0Kie3ekZlYKZNcD3PxshBD0aesg3MQ6KfzrAU3X5N6leL6Ic8lIuJ1ZUb5FbmkN6OqCFh21xNIaeUwKg4g7J903tQF6PAd0fRow0b25grKDl5xdgvXX14tyMXrQwyi/UZjXfh48rTxVPTyNQsbrNwvAxqDrBxC5ebP37kXG+vXIOXa84nFjb284zHkC1mPGVOkfeyvjFn6/8ju23NoiMn2Jwd6D8WzIs/CrFDumUxY/6t1LFj/ql0o4twb6vwm0Gqk2cX6NyfhOW7UKsm3bUFZUJB4rcG8JWffJSDD2QXK8lIkqx8rBFN5BDmjWxkFYCY1N7xYD12UK84uFdS/qUioiL6ciKzW/yt+pbp9fiJOw9Nk6a2GsJS1R17YD+z8CkqTqAaLPNVkEOz+hkxbBiMwIfHfuO+yO3C3uU3zg2OZjMaP1DDS3ba7q4WkEMh1fvwkWgI1AFw+gspIS5J4OhWznv8ja9R9K0tOlP+jpwaJHD9hOnQKrgQOhZyBZeujwOhp3VAi/Y3HHKl6nh3sPPN/heQQ56GCAtywOOL4UOLMKKJQSY+DQAui3AAia0KAWbupMcWqquEhIX72mIiyAKGvTGbJOY5Fo4IW4yByUFt+diih72NnHCu4t7YQYpOQFXRGEJPgSbmUK0Rd7PR3JkVkim1eOvqEePPxtRbY1iT5LOx0pmVJaAlzeBOz/GEi7JT1mbAl0nAV0fwawdoeuQaViyCJYeW7t6d5TCEGaY9miXjMyHVy/74UFYCPQlQOILH15oaGQ7d6NrP92oyS13FpV7ua1mTgBthMnwdjzbp2+5Nxk7Li9QzQ6p6tVgly9A7wGYGbQTIQ4h0DnIDfWqR+Bi+ukLgiEcxDQ6wVJ+Blot8AhK2DW3n3I3LIF2YcPA8WSFZgEr1FIZ+R2GIIUS3/ERhdBllLVykUxbQ4eFnDxsRYWL3J12rtaiMc1GYrZS0vIQdIdGRJvy0RnjtTYnIpYPjnWjqairA5ZSD0C7GBkosM14UqKgUsbpBqCSVKtPNENh7qJdJkjhU/oGGcTz4qL7H3R+4RrmCCvygT/CRjpNxKOZo6qHqLaIdOR9bs2WAA2Am0+gApjYpBz+DCyDx1GzsmTKMvNrfgbFW22HDwI1kOHwaJ7t4q+vfnF+SKhg1y8dEUqn4iopAv1t5weOF334lSogPOVzZLwiw29+7h3T6nobYtBGu/qbWhdSNm//4qEkfwL5d0gyjH284Ne90GQeYUgpcQecRHZ97k9CUNjfdi7W8LRwwIOnpYiw9XW1VzEvqmb5YOmWcrSzUjIRWpcNlJjspESm4O0uGwUF0rnyb3ucLJ8erS0g7u/LawdzVQybrWGlq4bu4GjXwORR+8+7tFJEoKUQKVjBaUpvnr11dXYdHOTKB1DUAtNsgaOaT5GJI6YUhF5Btq8ftcVFoCNQFsOIDoEiqKikBsaKty7uadPoyg2tspzDJwcYdm3ryT6unWFnpHUfYLatB2OPYx9UftwJPYI8orL+8WWFy+lSWeE7whYkqtGl6CG9+fXABf/uhvfR1aKoHFAl6cALy7hIKcoLg5Z+/Yje98+5Jw6ddcySBgawqxdO+h17IVstyCk6zkiOTYPSZGyaoWT+BcTA9g4mcHWyUwIJwtbE1GTkG5pM7MygqGRYi1oxUUlyMsqQk5Ggdio9h7dylLykJGch8zkPBQXlNQ4XudmVlWsm1b2vEjXi+jTwKkfgMub71rXKU4weBrQ/iHAtS10iazCLPx7+19svbW1oqg+YWZohl4evTCg2QD09ugt2s3pKjItWb8bAwtAHTyASmQy5F+6hLyLYcgLC0PexQsoSU6p+iQDA5iFtIdl7z7/b+9MgKM4zjb8SXvpPkAHCCQOgcE2hzjMYTvGKQjY4BgnKQeTVEyo+IzjwnEu7IrtcqpSxFfsMiHluJLgP1UxOFTZ5DeJ7cLYxj/3IQi3AhhJCHQggW5pd7U7f709u6vdRSskkNDszvsUTc/09Ejbmpmdt7+v+2tJueNr4hg/XllVYNX776X/yq7zu2Tb+W2yv2q/CtzsZ2jyUDUrDcJvZLrxQ5f0Kc01Ioc36MKv+nBnedpwkenLRaYuM20ss57iaWrSLc87dkjrjp1KHIYQHy+OceMksahI3IVF0pJRII0dKVJX2Sp153RLYU++0ay2eElIsYkjCcmq9hEqxWq3qO1w1zJcsh1ur3S4POJB7vaKs7VDnK1uaW92q/0rAaMkLHtYcq0zJUt6ThJXTOnLZ7D4f0T2rRVpDOrE5k7UheDE+0VSOgPQm4HShlIlBDd9tUkqWyo7n4E4q0wbMk1uz7tdZuXNkhsyb1DDdMxCY5S+v/sSCsAYvoE0r1dZ8pwlJdJ+okScJSdU7g/MHILNJokTJ6ogzUnTp0vilCliSUkWj9cjpxtOy8Gag7Knao/sqdwjl5y+iR8+xmSMUcGbEb/vpsE3Gc791q80Vooc/1Dk+P/qbiif21sFbx53t8jkpfrSbTE+vq/fLNNnz0rLzl3Sunu3tB08eLkgxJdYQoIkjBsnjptuFNvYceLOHS1tyUOkqSVOmi62B6xxzfXt0lrvCplQ0ZdAxCVl2CUlIyFgdYQlLz0nUc3Ohfgz1OobsT5O8NRmkYPv6kspIqg0gMDB8IubFouMv0ckbaiY6Xk6VndMtpRvkc/Pfi6n6k+FHM90ZMrMoTNVgOminCIpTC+M6fWHGw3+/r4eUADGwA0Eq4mrrFzc5WXiLC0V1+mvVCw215kzorVfPnYK2PLzleBLmDRREidNUqt0xCckqMkbxy8el6O1R+XghYNy6MIhaXb7Zqr6SLImyfQh02XW0Flyx/A7VGR604DHpfqo/nIp+Ujk7O7Q4xh/BEsDJnUkDRqoTxmzuKurpe3AQWk7cEDajh4R57Hj4g0anxqMNTdXHIWFYh85UuyjRqncNnKEaJk54mr3SnuLW5wtHeJs61AuXLiUdeueJ6Dj/UA3wG2sWwnj1bYj0SqOZKtaacORbBN7gsVcnZ9oofWiPmnkP+tEzu0POhAnkj9D76iNna+HYDLR9cPKS19WfCm7Knep1ZiCh++AFFuKTMqeJEXZRXJz1s1y46AbJTspdjwYjQZ5fw8kFIBRcAN5mluko7pK3OcrlUXPn1znKsRdfrYzFEsXYKyevbBQt5CMHy8J4/Xck5YkZxrOqF4g4vOduHhCjtcdl7r2zhm+wYJvYvZEmZozVWbnzZYJWRPEhvFsZqGlTqRsm75Y/clPRZrCrFDDZ4jcdK++Rm+micSwQazcrrIyaT96TNqPHxPnqVPiPHlSOs53urouw2oV27A8sQ/PF9vw4WIbPkxseXliG5ontryhYs3ODoQxIjHGpVLdYo8VRir2hh5LzRMZO0+fmDXidpHkwWIW3F63HKk9IjvP75TimmLV8Q8XhGBwwmC5cfCNMn7QeBVvEN6fUemjonJt4kYKQApAI95ACJPR8OEmXfRVVom3OdQC1xWWrCyxFxSoIMyOwtFqJqVl5Aipy7TK2bZzUt5Yrnp85U3lakxIRXNFYJZuMBgDgvABeMAxiQOuADzkCDRqKrfu2V0ipdtFSreJXDgeetyaKDLqa7prFwGb0zvD3xBj4GluVkLQ9dUZcZUilfpSWSAodUQsFiUCrbk5YsvJEWtOrgp3ZM3KEmt2lsotg7PEOigzMBmKRCEN5/QA07Dmn0FYojDBA4sg3MUjbxMpmC2Sap61yRGoH8YBDP3BJBIYCBDOK9I7Iz81X3mCClIL9DytQG1jPXejvjsaKQApAI14A1344x+l9s3VIWXxaWliGzJEbMOGqWTNGyru7AxpyE6UC4MsUiUNUt1SrQb5nm8+r1J1a7V4tK5nHoI0e5oSd+jJYQAwenbIMVPMVO6h6iO6a0il4tDB48Evg1FzdAsBXgo2E/2NYiyQOYJRuysqxHW2QtwVZ1XIo45KdLYqxV1VFToL+QrEp6eLNTNTLIMHK0FoycjwJf92ugqbhOcXuSUtTY1ZpKvYYLjb9DG8sPCf2dq54kgwacNEhk3Vh3kgz51gqmEesAhiAiA8RcjhOYJIxPJ0kUAImtykXMlLyVMJkwRzk3NVmT9hJvJAPA+NFIAUgEa8gSqKt0n1gR3SkGGXulSRmhSPXJAm5Z6ta6tTi4HXttUGllPrDrhqEXtvRKreK/P3ziD8YM43zYsI6+3WnRKpPSlSc1wXfRjL15XYw4CvQO//dn093mQGUjWNQKytlY7qal0oIq+uUdtY0QTHPDiOYOjeK8/87RKbTSypqSopYZiaIvHJKRKfkiLxqSliQZ6cHJqSkCepNbXjExNVHodks5nnGb6etNTq6xDDAwBhCEHYhfVLiUIIwdybRXJuFMkaq6/qY5J1iiEf8F6CEAx4mZA3lUlFU4VyLV8JWAgRqDo7MVsGJw5W76VBCYMCKTMhU7mZYU3sSxopAGNLAK5Zs0ZeeeUVqaqqksmTJ8vq1atlxowZEetv2LBBnnvuOSktLZWxY8fKSy+9JAsXLhzwG+iN/W/IX478pUd18YCgF5WTlKP3qJJzZVjKMJXQ48KDZYqp/VgmCiEgGs6KXCoTqS/Vc4z5gfBr6mZMWEaBvnqA6tlPExk6WcRhsriFpNdCEeGUsCoOglp7Ll6Sjot14qmv70yX6sXT0CDehga9bmOjiCeyRf6qsFh0QZiYKHFJiRKfgJQgcYm+PCFB4h0OPffvJzgkzo4yh37MkSBxDrvE2e2+fYfE2fz7eh6SIDrNNkbS2SxS+R89mDs8BVjVp748cv3UoboQzBypjwvO8OXp+XoYmhieXesH7mIYKuCNOtd8TiV4qeCZqmmtUfnF9os9+lkPTXxIVkxd0aefr5ECUIzpnL8K3nvvPXn66aflrbfekpkzZ8obb7whCxYskJKSEsnJuTzu044dO2Tp0qWyatUqueeee+Tdd9+V++67T4qLi2XChAkykAxLHaYsdP7ejz9Hz8jfU0KO3pId4UZiWdS11Yu0XdRdta21ushruaDnzdW6sMPauk1VeCt3//OSs0WybtB76ei1Izgseu0J5g2GSq4OCCC4fpF6OvwdfW1vS4t4m5rE09gk3iaIwibxtjSrmfzeZv0Y9lHPg7oqtep5W6to2G5t7RzH6PGoMcI9GSfcp1gsnWKwu2S16uMkbXoeZ/WVqXKrmpDTWWbx7VslzqIfV6sMYdt/XG1bdAHq21afxb8db/Hl8fq5vjwO62ujjiUoR1217/t5wfuqftCsbnQIMRYQyU97g+5NQNB3eBTgXaj9r/79hO8lpNL/u/xvF2fRxxNi7WIIxZRcXRTi+wl5UpbuWk5EyohasQjDAwwTSBhL3hUuj0t5tSAU/Z4tWBQvtV9S4tCfw5hB+p6YsQBC9N1yyy3yhz/8Qe17vV7Jz8+XJ598UlauXHlZ/SVLlkhLS4ts2rQpUDZr1iwpKipSIrInmL4HgVsHbmjE2FLJLdLh9KX2zhzja9ytQXmriKtFd8siV9uNIu2NIs4G/YvVn3oDLJ2Yyad63CM6c79bBl+mhMQAEIDetjY9QRD6t9vbVegnbxty7Ds7c2e7b79dNJezc9vtEq/TJZoTdZzidTlFc7n1fZcrkEwJBGC4IPTvRyqLwz+PWpEkDkHy4QbFttelb8dpqKJyva4ekUZpTV+Z/1erfQssrjbEIdK3fbleDhFt9+1D2Oo5ygLbKvnq+pMSuzb980Jg4vdjGxuqLE5vF75TVR6nB0fXP1RnO8PqdJ7nr+c7N7xe4NzQupeV+eoinBNSX9Jo9vd3rFgAXS6X7N+/X5555plAWXx8vMybN0927tzZ5Tkoh8UwGFgMN27cGPH3OJ1OlYJvoH4BIQoQqiCgzX252g/bDsl95Wqsiha27e3cD0+wtMF6pnKtcxviTiX/tlsPsKpy/UvtuoDliiDe0ENWvWTkOXpvGT1ojMNBQFeUMeAyMQGwplmQrtOLS9kJIDohDN0+Ueju8G2jzK2Xdfi2fQkTavRtX479Dl85Es71eDrLUA85vmc6UK7v4ztIw74q99XDGExVF3EbPSLhZXC3d5V7vYHtK47jRLvx8/y71/yXvNpZ4/jNEOHmFOJZS+ZJ9ouhEyPJtRMTb8taDMz2eCQ3NzekHPsnTpzo8hyME+yqPsojAXfxiy++KP1O9TF9SbFoBPGgsNi41Z/bRWxJvpTYmcOlgvWB7cl6cqSJJKTpYg8uWWwrF0gmRR0hA4yyzNjtYrFjyEmyxApK2AYJQs2jd4jVPo75BONluf8cfx2cpwWVqeNBZVh9Jvy47+dgX9/21/GVezpEw9jD9iblPdFc8JzAgwILbqtIh0v3sridormdPg8MBLWvg97h0oWx+v0dgVzZH4PtCGoD/+nlajPc1uA3S3Z5HvaDJiJdVkc/FuJrDPu9gd8ZfH7Qts3RdbB3cm3wzdoLYGEMthrCAgg3c58zZq4ugEJ9ARG2g8oCeXyYCR0mfn+5b0KIMvvHByV9DIzaVjnG1sBFoLsL9OO6K6Ez97kiIPYwFhHlnJFICIkSAu5bs0xqCR62A/Ho9+SovKMzD/YCoSzgIQryGIV7kdTP95f5PEnBHqpgz1S410rfiLyNYPukz4kJAZiFwKwWi1RXV4eUY3/IkK6njqO8N/WBw+FQqd/B8kRIhBBCSF+hBK+v4w7PCzE1MREfxG63y7Rp02TLli2BMkwCwf7s2bO7PAflwfXB5s2bI9YnhBBCCIkVYsICCOCaXbZsmUyfPl3F/kMYGMzyXb58uTr+4IMPyrBhw9Q4PrBixQqZM2eOvPbaa7Jo0SJZv3697Nu3T95+++0BbgkhhBBCSP8SMwIQYV0uXLggzz//vJrIgXAuH3/8cWCiR3l5uZoZ7OfWW29Vsf9+/etfy7PPPqsCQWMG8EDHACSEEEII6W9iJg7gQMA4QoQQQkj00cj3d2yMASSEEEIIIT2HApAQQgghxGRQABJCCCGEmAwKQEIIIYQQk0EBSAghhBBiMigACSGEEEJMBgUgIYQQQojJoAAkhBBCCDEZFICEEEIIISYjZpaCGwj8i6ggojghhBBCooNG33vbzIuhUQBeA01NTSrPz88f6I9CCCGEkKt4j6enp4sZ4VrA14DX65Xz589LamqqxMXF9XnvBMLy7NmzMblOIdsX/cR6G9m+6CfW28j2XT2apinxl5eXJ/Hx5hwNRwvgNYCbZvjw4f36O3DTx+KD7Yfti35ivY1sX/QT621k+66OdJNa/vyYU/YSQgghhJgYCkBCCCGEEJNBAWhQHA6HvPDCCyqPRdi+6CfW28j2RT+x3ka2j1wLnARCCCGEEGIyaAEkhBBCCDEZFICEEEIIISaDApAQQgghxGRQABJCCCGEmAwKQANQWloqP/rRj2TUqFGSmJgohYWFauaTy+Xq9rz29nZ54oknZPDgwZKSkiLf+c53pLq6WozKb3/7W7n11lslKSlJMjIyenTOD3/4Q7XKSnC66667JFbahzlYzz//vAwdOlRd+3nz5snJkyfFiFy8eFG+//3vq4CsaB/u2ebm5m7PufPOOy+7fo899pgYhTVr1sjIkSMlISFBZs6cKXv27Om2/oYNG2T8+PGq/sSJE+Xf//63GJnetO+dd9657FrhPKPy5Zdfyje/+U21kgM+68aNG694zhdffCFTp05Vs0rHjBmj2mxkettGtC/8GiJVVVWJEVm1apXccsstajWtnJwcue+++6SkpOSK50Xbc2hUKAANwIkTJ9Sycn/605/k6NGj8vrrr8tbb70lzz77bLfn/fSnP5UPP/xQPQxbt25Vy9J9+9vfFqMCQXv//ffL448/3qvzIPgqKysDad26dRIr7Xv55ZflzTffVNd79+7dkpycLAsWLFDi3mhA/OH+3Lx5s2zatEm9nB555JErnvfwww+HXD+02Qi899578vTTT6vOVnFxsUyePFn97Wtqarqsv2PHDlm6dKkSvgcOHFAvK6QjR46IEelt+wDEffC1KisrE6PS0tKi2gSR2xPOnDkjixYtkq9//ety8OBBeeqpp+Shhx6STz75RGKljX4gooKvI8SVEcF7C0aMXbt2qe8Vt9st8+fPV+2ORLQ9h4YGYWCI8Xj55Ze1UaNGRTxeX1+v2Ww2bcOGDYGy48ePI6SPtnPnTs3IrF27VktPT+9R3WXLlmmLFy/Woomets/r9WpDhgzRXnnllZDr6nA4tHXr1mlG4tixY+re2rt3b6Dso48+0uLi4rRz585FPG/OnDnaihUrNCMyY8YM7YknngjsezweLS8vT1u1alWX9b/73e9qixYtCimbOXOm9uijj2qx0L7ePJdGA/fmBx980G2dX/7yl9rNN98cUrZkyRJtwYIFWqy08fPPP1f1Ll26pEUjNTU16vNv3bo1Yp1oew6NDC2ABqWhoUEGDRoU8fj+/ftVbwkuQz8wiRcUFMjOnTslloBbAz3YcePGKetaXV2dxAKwSMA1E3wNsTYlXHVGu4b4PHD7Tp8+PVCGz431sGG57I6///3vkpWVJRMmTJBnnnlGWltbxQjWWjxDwX97tAX7kf72KA+uD2BRM9q1utr2Abj0R4wYIfn5+bJ48WJl8Y0Voun6XStFRUVqWMk3vvEN2b59u0TTew909+4z03Xsb6z9/htIrzl16pSsXr1aXn311Yh1IBzsdvtlY81yc3MNO97jaoD7F25tjI88ffq0covffffd6mG3WCwSzfivE66Z0a8hPk+4G8lqtaov6u4+6/e+9z0lKDCG6dChQ/KrX/1Kuafef/99GUhqa2vF4/F0+bfHkIyuQDuj4VpdbfvQwfrrX/8qkyZNUi9ifP9gTCtE4PDhwyXaiXT9Ghsbpa2tTY3BjXYg+jCcBB01p9Mpf/7zn9U4XHTSMPbRyGAYFNzyt912m+osRiKankOjQwtgP7Jy5couB+QGp/Av43PnzinRg7FkGDsVi23sDQ888IDce++9aqAvxnlg7NnevXuVVTAW2jfQ9Hf7MEYQvXNcP4wh/Nvf/iYffPCBEvPEWMyePVsefPBBZT2aM2eOEunZ2dlqbDKJDiDiH330UZk2bZoS7xD0yDGu3OhgLCDG8a1fv36gP4ppoAWwH/nZz36mZrF2x+jRowPbmMSBAcp4YN9+++1uzxsyZIhy89TX14dYATELGMeM2sZrBT8L7kRYSefOnSvR3D7/dcI1Q8/dD/bxEr4e9LR9+Kzhkwc6OjrUzODe3G9wbwNcP8x2HyhwD8GCHD5rvrvnB+W9qT+QXE37wrHZbDJlyhR1rWKBSNcPE19iwfoXiRkzZsi2bdvEyPzkJz8JTCy7krU5mp5Do0MB2I+g94zUE2D5g/hDz23t2rVqvE53oB6+oLds2aLCvwC41srLy1VP3oht7AsqKirUGMBgwRSt7YNbG19auIZ+wQd3FNw1vZ0p3d/twz2FzgbGleHeA5999ply2/hFXU/A7Etwva5fJDB8Au3A3x6WZYC2YB8vo0h/AxyHm8oPZi5ez+etP9sXDlzIhw8floULF0osgOsUHi7EqNevL8EzN9DPWyQwt+XJJ59UXgF4dfCdeCWi6Tk0PAM9C4VoWkVFhTZmzBht7ty5aruysjKQguuMGzdO2717d6Dsscce0woKCrTPPvtM27dvnzZ79myVjEpZWZl24MAB7cUXX9RSUlLUNlJTU1OgDtr4/vvvq22U//znP1ezms+cOaN9+umn2tSpU7WxY8dq7e3tWrS3D/zud7/TMjIytH/+85/aoUOH1IxnzP5ua2vTjMZdd92lTZkyRd2D27ZtU9dh6dKlEe/RU6dOab/5zW/UvYnrhzaOHj1au+OOOzQjsH79ejXj+p133lGznB955BF1LaqqqtTxH/zgB9rKlSsD9bdv365ZrVbt1VdfVTPuX3jhBTUT//Dhw5oR6W37cN9+8skn2unTp7X9+/drDzzwgJaQkKAdPXpUMyJ4rvzPGF5lv//979U2nkOAtqGNfr766istKSlJ+8UvfqGu35o1azSLxaJ9/PHHmlHpbRtff/11bePGjdrJkyfVfYkZ+PHx8eq704g8/vjjaub5F198EfLea21tDdSJ9ufQyFAAGgCEX8DD3VXygxco9jHN3w9Ewo9//GMtMzNTfbF961vfChGNRgMhXbpqY3CbsI+/B8CXwPz587Xs7Gz1gI8YMUJ7+OGHAy+waG+fPxTMc889p+Xm5qqXNToBJSUlmhGpq6tTgg/iNi0tTVu+fHmIuA2/R8vLy5XYGzRokGobOjl4+TY0NGhGYfXq1aoTZbfbVdiUXbt2hYSwwTUN5h//+Id2ww03qPoIKfKvf/1LMzK9ad9TTz0VqIv7ceHChVpxcbFmVPwhT8KTv03I0cbwc4qKilQb0RkJfhaNSG/b+NJLL2mFhYVKuOO5u/POO5WBwKhEeu8FX5dYeA6NShz+G2grJCGEEEIIuX5wFjAhhBBCiMmgACSEEEIIMRkUgIQQQgghJoMCkBBCCCHEZFAAEkIIIYSYDApAQgghhBCTQQFICCGEEGIyKAAJIYQQQkwGBSAhhBBCiMmgACSEEEIIMRkUgIQQQgghJoMCkBBCCCHEZFAAEkIIIYSYDApAQgghhBCTQQFICCGEEGIyKAAJIYQQQkwGBSAhhBBCiMmgACSEEEIIMRkUgIQQQgghJoMCkBBCCCHEZFAAEkIIIYSYDApAQgghhBCTQQFICCGEEGIyKAAJIYQQQkwGBSAhhBBCiMmgACSEEEIIMRkUgIQQQgghJoMCkBBCCCHEZFAAEkIIIYSYDApAQgghhBCTQQFICCGEECLm4v8BcKdXHPT0xB0AAAAASUVORK5CYII=", + "text/html": [ + "\n", + "
\n", + "
\n", + " Figure\n", + "
\n", + " \n", + "
\n", + " " + ], + "text/plain": [ + "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sample_model=SampleModel(name='sample_model')\n", + "\n", + "# Creating components\n", + "gaussian=Gaussian(name='Gaussian',width=0.5,area=1)\n", + "dho = DampedHarmonicOscillator(name='DHO',center=1.0,width=0.3,area=2.0)\n", + "lorentzian = Lorentzian(name='Lorentzian',center=-1.0,width=0.2,area=1.0)\n", + "polynomial = Polynomial(name='Polynomial',coefficients=[0.1, 0, 0.5]) # y=0.1+0.5*x^2\n", + "\n", + "sample_model.add_component(gaussian)\n", + "sample_model.add_component(dho)\n", + "sample_model.add_component(lorentzian)\n", + "sample_model.add_component(polynomial)\n", + "\n", + "\n", + "x=np.linspace(-2, 2, 100)\n", + "\n", + "plt.figure()\n", + "y=sample_model.evaluate(x)\n", + "plt.plot(x, y, label='Sample Model')\n", + "\n", + "for component in list(sample_model):\n", + " y = component.evaluate(x)\n", + " plt.plot(x, y, label=component.name)\n", + "\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "d35179d8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[Gaussian(name = Gaussian, unit = meV,\n", + " area = ,\n", + " center = ,\n", + " width = ),\n", + " DampedHarmonicOscillator(name = DHO, unit = meV,\n", + " area = ,\n", + " center = ,\n", + " width = ),\n", + " Lorentzian(name = Lorentzian, unit = meV,\n", + " area = ,\n", + " center = ,\n", + " width = ),\n", + " Polynomial(name = Polynomial, unit = meV,\n", + " coefficients = [Polynomial_c0=0.1, Polynomial_c1=0.0, Polynomial_c2=0.5])]" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# sample_model=SampleModel(name='sample_model')\n", + "sample_model.components" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "newdynamics", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/easydynamics/sample_model/__init__.py b/src/easydynamics/sample_model/__init__.py index a64ffd2..2b1274f 100644 --- a/src/easydynamics/sample_model/__init__.py +++ b/src/easydynamics/sample_model/__init__.py @@ -6,6 +6,7 @@ Polynomial, Voigt, ) +from .sample_model import SampleModel __all__ = [ "SampleModel", diff --git a/src/easydynamics/sample_model/sample_model.py b/src/easydynamics/sample_model/sample_model.py new file mode 100644 index 0000000..ff23f7c --- /dev/null +++ b/src/easydynamics/sample_model/sample_model.py @@ -0,0 +1,306 @@ +import warnings +from typing import List, Optional, Union + +import numpy as np +import scipp as sc +from easyscience.base_classes import CollectionBase +from easyscience.global_object.undo_redo import NotarizedDict +from easyscience.job.theoreticalmodel import TheoreticalModelBase + +from easydynamics.sample_model.components import DeltaFunction + +from .components.model_component import ModelComponent + +Numeric = Union[float, int] + + +class SampleModel(CollectionBase, TheoreticalModelBase): + """ + A model of the scattering from a sample, combining multiple model components. + + Attributes + ---------- + name : str + Name of the SampleModel. + unit : str or sc.Unit + Unit of the SampleModel. + components : List[ModelComponent] + List of model components in the SampleModel. + + """ + + def __init__( + self, + name: str = "MySampleModel", + unit: Optional[Union[str, sc.Unit]] = "meV", + data: Optional[List] = None, + ): + """ + Initialize a new SampleModel. + + Parameters + ---------- + name : str + Name of the sample model. + unit : str or sc.Unit, optional + Unit of the sample model. Defaults to "meV". + data : List[ModelComponent], optional + Initial list of model components to include in the sample model. + """ + + CollectionBase.__init__(self, name=name) + TheoreticalModelBase.__init__(self, name=name) + if not isinstance(self._kwargs, NotarizedDict): + self._kwargs = NotarizedDict() + + self._unit = unit + + # Add initial components if provided. Mostly used for serialization. + if data: + # Just to be safe + self.clear_components() + for item in data: + # ensure item is a ModelComponent + if not isinstance(item, ModelComponent): + raise TypeError("Data items must be instances of ModelComponent.") + self.insert(index=len(self), value=item) + + def add_component( + self, component: ModelComponent, name: Optional[str] = None + ) -> None: + """ + Add a model component to the SampleModel. Component names must be unique. + Parameters + ---------- + component : ModelComponent + The model component to add. + name : str, optional + Name to assign to the component. If None, uses the component's own name. + """ + if not isinstance(component, ModelComponent): + raise TypeError("component must be an instance of ModelComponent.") + + if name is None: + name = component.name + if name in self.list_component_names(): + raise ValueError(f"Component with name '{name}' already exists.") + + component.name = name + + self.insert(index=len(self), value=component) + + def remove_component(self, name: str): + """ + Remove a model component by name. + """ + # Find index where item.name == name + indices = [i for i, item in enumerate(list(self)) if item.name == name] + if not indices: + raise KeyError(f"No component named '{name}' exists in the model.") + del self[indices[0]] + + def list_component_names(self) -> List[str]: + """ + List the names of all components in the model. + + Returns + ------- + List[str] + Component names. + """ + + return [item.name for item in list(self)] + + def clear_components(self): + """ + Remove all components from the model. + """ + + for _ in range(len(self)): + del self[0] + + def normalize_area(self) -> None: + # Useful for convolutions. + """ + Normalize the areas of all components so they sum to 1. + """ + if not self.components: + raise ValueError("No components in the model to normalize.") + + area_params = [] + total_area = 0.0 + + for component in list(self): + if hasattr(component, "area"): + area_params.append(component.area) + total_area += component.area.value + else: + warnings.warn( + f"Component '{component.name}' does not have an 'area' attribute and will be skipped in normalization." + ) + + if total_area == 0: + raise ValueError("Total area is zero; cannot normalize.") + + if not np.isfinite(total_area): + raise ValueError("Total area is not finite; cannot normalize.") + + for param in area_params: + param.value /= total_area + + @property + def components(self) -> List[ModelComponent]: + """ + Get the list of components in the SampleModel. + + Returns + ------- + List[ModelComponent] + """ + return list(self) + + @property + def unit(self) -> Optional[Union[str, sc.Unit]]: + """ + Get the unit of the SampleModel. + + Returns + ------- + str or sc.Unit or None + """ + return self._unit + + @unit.setter + def unit(self, unit_str: str) -> None: + raise AttributeError( + ( + f"Unit is read-only. Use convert_unit to change the unit between allowed types " + f"or create a new {self.__class__.__name__} with the desired unit." + ) + ) # noqa: E501 + + def convert_unit(self, unit: Union[str, sc.Unit]) -> None: + """ + Convert the unit of the SampleModel and all its components. + """ + self._unit = unit + # for component in self.components.values(): + for component in list(self): + component.convert_unit(unit) + + def evaluate( + self, x: Union[Numeric, list, np.ndarray, sc.Variable, sc.DataArray] + ) -> np.ndarray: + """ + Evaluate the sum of all components. + + Parameters + ---------- + x : Number, list, np.ndarray, sc.Variable, or sc.DataArray + Energy axis. + + Returns + ------- + np.ndarray + Evaluated model values. + """ + + if not self.components: + raise ValueError("No components in the model to evaluate.") + result = None + for component in list(self): + value = component.evaluate(x) + result = value if result is None else result + value + + return result + + def evaluate_without_delta( + self, x: Union[Numeric, list, np.ndarray, sc.Variable, sc.DataArray] + ) -> np.ndarray: + """ + Evaluate the sum of all components except delta functions. + + Parameters + ---------- + x : Number, list, np.ndarray, sc.Variable, or sc.DataArray + Energy axis. + + Returns + ------- + np.ndarray + Evaluated model values. + """ + + if not self.components: + raise ValueError("No components in the model to evaluate.") + result = None + for component in list(self): + if not isinstance(component, DeltaFunction): + value = component.evaluate(x) + result = value if result is None else result + value + + return result + + def evaluate_component( + self, + x: Union[Numeric, list, np.ndarray, sc.Variable, sc.DataArray], + name: str, + ) -> np.ndarray: + """ + Evaluate a single component by name. + + Parameters + ---------- + x : Number, list, np.ndarray, sc.Variable, or sc.DataArray + Energy axis. + name : str + Component name. + + Returns + ------- + np.ndarray + Evaluated values for the specified component. + """ + if not self.components: + raise ValueError("No components in the model to evaluate.") + + if not isinstance(name, str): + raise TypeError( + (f"Component name must be a string, got {type(name)} instead.") + ) + + matches = [comp for comp in list(self) if comp.name == name] + if not matches: + raise KeyError(f"No component named '{name}' exists.") + + component = matches[0] + + result = component.evaluate(x) + + return result + + def fix_all_parameters(self) -> None: + """ + Fix all free parameters in the model. + """ + for param in self.get_parameters(): + param.fixed = True + + def free_all_parameters(self) -> None: + """ + Free all fixed parameters in the model. + """ + for param in self.get_parameters(): + param.fixed = False + + def __repr__(self) -> str: + """ + Return a string representation of the SampleModel. + + Returns + ------- + str + """ + comp_names = ", ".join(c.name for c in self) or "No components" + + return f"" diff --git a/src/easydynamics/utils/__init__.py b/src/easydynamics/utils/__init__.py index 9cf350f..a6bd0bf 100644 --- a/src/easydynamics/utils/__init__.py +++ b/src/easydynamics/utils/__init__.py @@ -1,3 +1,4 @@ +from .convolution import convolution from .detailed_balance import _detailed_balance_factor -__all__ = ["_detailed_balance_factor"] +__all__ = ["_detailed_balance_factor", "convolution"] diff --git a/src/easydynamics/utils/convolution.py b/src/easydynamics/utils/convolution.py new file mode 100644 index 0000000..a282be8 --- /dev/null +++ b/src/easydynamics/utils/convolution.py @@ -0,0 +1,749 @@ +import warnings +from typing import List, Optional, Tuple, Union + +import numpy as np +import scipp as sc +from easyscience.variable import Parameter +from scipy.signal import fftconvolve +from scipy.special import voigt_profile + +from easydynamics.sample_model import DeltaFunction, Gaussian, Lorentzian, SampleModel +from easydynamics.sample_model.components.model_component import ModelComponent +from easydynamics.utils.detailed_balance import ( + _detailed_balance_factor as detailed_balance_factor, +) + +Numerical = Union[float, int] + + +def convolution( + energy: np.ndarray, + sample_model: Union[SampleModel, ModelComponent], + resolution_model: Union[SampleModel, ModelComponent], + offset: Optional[Union[Parameter, float, None]] = None, + method: Optional[str] = "auto", + upsample_factor: Optional[int] = 0, + extension_factor: Optional[float] = 0.2, + temperature: Optional[Union[Parameter, float, None]] = None, + temperature_unit: Union[str, sc.Unit] = "K", + energy_unit: Optional[Union[str, sc.Unit]] = "meV", + normalize_detailed_balance: Optional[bool] = True, +) -> np.ndarray: + """ + Calculate the convolution of a sample model with a resolution model using analytical expressions or numerical FFT. + Accepts SampleModel or ModelComponent for both sample and resolution. + If method is 'auto', analytical convolution is preferred when possible, otherwise numerical convolution is used. + Detailed balancing is included if temperature is provided. This requires numerical convolution and that the units + of energy and temperature are provided. An error will be raised if the units are not compatible. + The calculated model is shifted by the specified offset. + + Examples: + energy = np.linspace(-10, 10, 100) + sample = SampleModel() + sample.add_component(Gaussian(name="SampleGaussian", area=1.0, center=0.1, width=1.0)) + resolution = Gaussian(name="ResolutionGaussian", area=1.0, center=0.0, width=0.5) + result = convolution(energy, sample, resolution, offset=0.2) + + energy = np.linspace(-10, 10, 100) + sample = SampleModel() + sample.add_component(Gaussian(name="Gaussian", area=1.0, center=0.1, width=1.0)) + sample.add_component(DampedHarmonicOscillator(name="DHO", area=2.0, center=1.5, width=0.2)) + sample.add_component(DeltaFunction(name="Delta", area=0.5, center=0.0)) + + resolution = SampleModel() + resolution.add_component(Gaussian(name="ResolutionGaussian", area=0.8, center=0.0, width=0.5)) + resolution.add_component(Lorentzian(name="ResolutionLorentzian", area=0.2, center=0.1, width=0.3)) + + result_auto = convolution(energy, sample, resolution, offset=0.2, method='auto', upsample_factor=5, extension_factor=0.2) + result_numerical = convolution(energy, sample, resolution, offset=0.2, method='numerical', upsample_factor=5, extension_factor=0.2) + + + Args: + energy : np.ndarray + 1D array of energy transfer where the convolution is evaluated. + sample_model : SampleModel or ModelComponent + The sample model to be convolved. + resolution_model : SampleModel or ModelComponent + The resolution model to convolve with. + offset : Parameter, float, or None, optional + The offset to apply to the x values before convolution. + method : str, optional + The convolution method to use: 'auto', 'analytical' or 'numerical'. Default is 'auto'. + upsample_factor : int, optional + The factor by which to upsample the input data before numerical convolution. Default is 0 (no upsampling). + extension_factor : float, optional + The factor by which to extend the input data range before numerical convolution. Default is 0.2. + temperature : Parameter, float, or None, optional + The temperature to use for detailed balance calculations. Default is None. + temperature_unit : str or sc.Unit, optional + The unit of the temperature parameter. Default is 'K'. + energy_unit : str or sc.Unit, optional + The unit of the energy. Default is 'meV'. + normalize_detailed_balance : bool, optional + Whether to normalize the detailed balance factor. Default is True. + """ + + # Input validation + if not isinstance(energy, np.ndarray): + raise TypeError( + f"`energy` is an instance of {type(energy).__name__}, but must be a numpy array." + ) + + energy = np.asarray(energy, dtype=float) + if energy.ndim != 1 or not np.all(np.isfinite(energy)): + raise ValueError("`energy` must be a 1D finite array.") + + if not isinstance(sample_model, (SampleModel, ModelComponent)): + raise TypeError( + f"`sample_model` is an instance of {type(sample_model).__name__}, but must be SampleModel or ModelComponent." + ) + + if not isinstance(resolution_model, (SampleModel, ModelComponent)): + raise TypeError( + f"`resolution_model` is an instance of {type(resolution_model).__name__}, but must be SampleModel or ModelComponent." + ) + + if isinstance(sample_model, SampleModel): + if not sample_model.components: + raise ValueError("SampleModel must have at least one component.") + + if isinstance(resolution_model, SampleModel): + if not resolution_model.components: + raise ValueError("ResolutionModel must have at least one component.") + + # Handle offset + if offset is None: + offset_float = 0.0 + elif isinstance(offset, Parameter): + offset_float = offset.value + elif isinstance(offset, Numerical): + offset_float = float(offset) + else: + raise TypeError( + f"Expected offset to be Parameter, number, or None, got {type(offset)}" + ) + + if not isinstance(upsample_factor, int) or upsample_factor < 0: + raise ValueError("upsample_factor must be a non-negative integer.") + + if not isinstance(extension_factor, float) or extension_factor < 0.0: + raise ValueError("extension_factor must be a non-negative float.") + + if temperature is not None: + if energy_unit is None: + raise ValueError( + "energy_unit must be provided when temperature is specified." + ) + if not isinstance(energy_unit, (str, sc.Unit)): + raise TypeError( + f"Expected energy_unit to be str or sc.Unit, got {type(energy_unit)}" + ) + + use_numerical_convolution_as_fallback = False + if method == "auto": + if temperature is not None: + method = "numerical" + else: + method = "analytical" + use_numerical_convolution_as_fallback = True + + if method == "analytical": + if temperature is not None: + raise ValueError( + "Analytical convolution is not supported with detailed balance. Set method to 'numerical' instead or set the temperature to None." + ) + return _analytical_convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + offset_float=offset_float, + use_numerical_convolution_as_fallback=use_numerical_convolution_as_fallback, + upsample_factor=upsample_factor, + extension_factor=extension_factor, + ) + elif method == "numerical": + return _numerical_convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + offset_float=offset_float, + upsample_factor=upsample_factor, + extension_factor=extension_factor, + temperature=temperature, + temperature_unit=temperature_unit, + energy_unit=energy_unit, + normalize_detailed_balance=normalize_detailed_balance, + ) + else: + raise ValueError( + f"Unknown convolution method: {method}. Choose from 'auto', 'analytical', or 'numerical'." + ) + + +def _numerical_convolution( + energy: np.ndarray, + sample_model: Union[SampleModel, ModelComponent], + resolution_model: Union[SampleModel, ModelComponent], + offset_float: Optional[float] = 0.0, + upsample_factor: Optional[int] = 5, + extension_factor: Optional[float] = 0.2, + temperature: Optional[Union[Parameter, float]] = None, + temperature_unit: Optional[Union[str, sc.Unit]] = "K", + energy_unit: Optional[Union[str, sc.Unit]] = "meV", + normalize_detailed_balance: Optional[bool] = True, +) -> np.ndarray: + """ + Numerical convolution using FFT with optional upsampling + extended range. + Includes detailed balance correction if temperature is provided. + + + Args: + energy : np.ndarray + 1D array of energy values where the convolution is evaluated. + sample_model : SampleModel or ModelComponent + The sample model to be convolved. + resolution_model : SampleModel or ModelComponent + The resolution model to convolve with. + offset_float : float, or None, optional + The offset to apply to the input array. + upsample_factor : int, optional + The factor by which to upsample the input data before convolution. Default is 5. + extension_factor : float, optional + The factor by which to extend the input data range before convolution. Default is 0.2. + temperature : Parameter, float, or None, optional + The temperature to use for detailed balance correction. Default is None. + temperature_unit : str or sc.Unit, optional + The unit of the temperature parameter. Default is 'K'. + energy_unit : str or sc.Unit, optional + The unit of the energy. Default is 'meV'. + normalize_detailed_balance : bool, optional + Whether to normalize the detailed balance factor. Default is True. + Returns: + np.ndarray + The convolved values evaluated at energy. + """ + + # Create a dense grid to improve accuracy. We evaluate on this grid and interpolate back to the original values at the end + energy_dense = _create_dense_grid( + energy, upsample_factor=upsample_factor, extension_factor=extension_factor + ) + + energy_step = energy_dense[1] - energy_dense[0] + span = energy_dense.max() - energy_dense.min() + # Handle offset for even length of x in convolution. + # The convolution of two arrays of length N is of length 2N-1. When using 'same' mode, only the central N points are kept, + # so the output has the same length as the input. + # However, if N is even, the center falls between two points, leading to a half-bin offset. + # For example, if N=4, the convolution has length 7, and when we select the 4 central points we either get + # indices [2,3,4,5] or [1,2,3,4], both of which are offset by 0.5*dx from the true center at index 3.5. + if len(energy_dense) % 2 == 0: + x_even_length_offset = -0.5 * energy_step + else: + x_even_length_offset = 0.0 + + # Handle the case when x is not symmetric around zero. The resolution is still centered around zero (or close to it), so it needs to be evaluated there. + if not np.isclose(energy_dense.mean(), 0.0): + energy_dense_centered = np.linspace(-0.5 * span, 0.5 * span, len(energy_dense)) + else: + energy_dense_centered = energy_dense + + # Give warnings if peaks are very wide or very narrow + _check_width_thresholds( + model=sample_model, + span=span, + energy_step=energy_step, + model_name="sample model", + ) + _check_width_thresholds( + model=resolution_model, + span=span, + energy_step=energy_step, + model_name="resolution model", + ) + + # Evaluate sample model. Delta functions are handled separately for accuracy. + if isinstance(sample_model, SampleModel): + sample_vals = sample_model.evaluate_without_delta( + energy_dense - offset_float - x_even_length_offset + ) + elif isinstance(sample_model, DeltaFunction): + sample_vals = np.zeros_like(energy_dense) + else: + sample_vals = sample_model.evaluate( + energy_dense - offset_float - x_even_length_offset + ) + + # Detailed balance correction + if temperature is not None: + detailed_balance_factor_correction = detailed_balance_factor( + energy=energy_dense, + temperature=temperature, + energy_unit=energy_unit, + temperature_unit=temperature_unit, + divide_by_temperature=normalize_detailed_balance, + ) + sample_vals *= detailed_balance_factor_correction + + # Evaluate resolution model + if isinstance(resolution_model, SampleModel): + resolution_vals = resolution_model.evaluate_without_delta(energy_dense_centered) + elif isinstance(resolution_model, DeltaFunction): + resolution_vals = np.zeros_like(energy_dense_centered) + else: + resolution_vals = resolution_model.evaluate(energy_dense_centered) + + # Convolution + convolved = fftconvolve(sample_vals, resolution_vals, mode="same") + convolved *= energy_step # normalize + + if upsample_factor > 0: + # interpolate back to original energy grid + convolved = np.interp(energy, energy_dense, convolved, left=0.0, right=0.0) + + # Add delta function contributions + delta_contributions = _calculate_delta_contributions( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + offset_float=offset_float, + ) + convolved += delta_contributions + + return convolved + + +def _analytical_convolution( + energy: np.ndarray, + sample_model: Union[SampleModel, ModelComponent], + resolution_model: Union[SampleModel, ModelComponent], + offset_float: float = 0.0, + use_numerical_convolution_as_fallback: bool = False, + upsample_factor: int = 5, + extension_factor: float = 0.2, +) -> np.ndarray: + """ + Convolve sample with resolution analytically if possible. Accepts SampleModel or single ModelComponent for each. + Possible analytical convolutions are any combination of delta functions, Gaussians, and Lorentzians. + Falls back to numerical convolution for other pairs of functions + + Most validation happens in the main `convolution` function. + + Args: + x : np.ndarray + 1D array of x values where the convolution is evaluated. + sample_model : SampleModel or ModelComponent + The sample model to be convolved. + resolution_model : SampleModel or ModelComponent + The resolution model to convolve with. + offset_float : float + The offset to apply to the convolution. + use_numerical_convolution_as_fallback : bool + Whether to use numerical convolution as a fallback if analytical convolution is not possible. Default is False. Is True when method='auto'. + upsample_factor : int, optional + The factor by which to upsample the input data before numerical convolution. Improves accuracy at the cost of speed. Default is 5 + extension_factor : float, optional + The factor by which to extend the input data range before numerical convolution. Improves accuracy at the edges of the data. Default is 0.2 + Returns: + np.ndarray + The convolved values evaluated at x. + + Raises: + ValueError + If both sample_model and resolution_model contain delta functions. + + """ + + # prepare list of components + if isinstance(sample_model, SampleModel): + sample_components = sample_model.components + else: + sample_components = [sample_model] + + if isinstance(resolution_model, SampleModel): + resolution_components = resolution_model.components + else: + resolution_components = [resolution_model] + + total = np.zeros_like(energy, dtype=float) + + # loop over sample components. Try to convolve each with all resolution components analytically + for sample_component in sample_components: + not_analytical_components = SampleModel(name="not_analytical") + + # Go through resolution components, adding analytical contributions where possible, making a list of those that cannot be handled analytically + for resolution_component in resolution_components: + handled, contrib = _try_analytic_pair( + energy=energy, + sample_component=sample_component, + resolution_component=resolution_component, + offset_float=offset_float, + ) + if handled: + total += contrib + else: + not_analytical_components.add_component(resolution_component) + + if not_analytical_components: + if use_numerical_convolution_as_fallback: + total += _numerical_convolution( + energy=energy, + sample_model=sample_component, + resolution_model=not_analytical_components, + offset_float=offset_float, + upsample_factor=upsample_factor, + extension_factor=extension_factor, + ) + else: + raise ValueError( + f"Could not find analytical convolution for sample component '{sample_component.name}' with resolution model '{not_analytical_components.name}'. " + "Set method to 'auto' or 'numerical'." + ) + + return total + + +# ---------------------- helpers & evals ----------------------- + + +def _create_dense_grid( + energy: np.ndarray, upsample_factor: int = 5, extension_factor: float = 0.2 +) -> np.ndarray: + """ + Create a dense grid by upsampling and extending the input energy array. + + Args: + energy : np.ndarray + 1D array of energy values. + upsample_factor : int, optional + The factor by which to upsample the input data. Default is 5. + extension_factor : float, optional + The factor by which to extend the input data range. Default is 0.2. + Returns: + np.ndarray + The dense grid created by upsampling and extending x. + """ + if upsample_factor == 0: + # Check if the array is uniformly spaced. + energy_diff = np.diff(energy) + is_uniform = np.allclose(energy_diff, energy_diff[0]) + if not is_uniform: + raise ValueError( + "Input array `energy` must be uniformly spaced if upsample_factor = 0." + ) + energy_dense = energy + else: + # Create an extended and upsampled energy grid + energy_min, energy_max = energy.min(), energy.max() + span = energy_max - energy_min + extra = extension_factor * span + extended_min = energy_min - extra + extended_max = energy_max + extra + num_points = len(energy) * upsample_factor + energy_dense = np.linspace(extended_min, extended_max, num_points) + + return energy_dense + + +def _try_analytic_pair( + energy: np.ndarray, + sample_component: Union[ModelComponent, SampleModel], + resolution_component: Union[ModelComponent, SampleModel], + offset_float: float, +) -> Tuple[bool, np.ndarray]: + """ + Attempt an analytic convolution for component pair (sample_component, resolution_component). + Returns (True, contribution) if handled, else (False, zeros). + The convolution of two gaussian components results in another gaussian component with width sqrt(w1^2 + w2^2). + The convolution of two lorentzian components results in another lorentzian component with width w1 + w2. + The convolution of a gaussian and a lorentzian results in a voigt profile. + The convolution of a delta function with any component or SampleModel results in the same component or SampleModel shifted by the delta center. + All areas are multiplied. + + + Args: + energy: np.ndarray + 1D array of energy values where the convolution is evaluated. + sample_component : Union[ModelComponent, SampleModel] + The sample component to be convolved. + resolution_component : Union[ModelComponent, SampleModel] + The resolution component to convolve with. + offset_float : float + The offset in energyto apply to the convolution. + + Returns: + Tuple[bool, np.ndarray]: + - bool: True if analytical convolution was computed, False otherwise + - np.ndarray: The convolution result if computed, or zeros if not handled + """ + # Two delta functions is not meaningful + if isinstance(sample_component, DeltaFunction) and isinstance( + resolution_component, DeltaFunction + ): + raise ValueError("Convolution of two delta functions is not defined.") + + # Delta function + anything --> anything, shifted by delta center with area A1 * A2 + if isinstance(sample_component, DeltaFunction): + return True, sample_component.area.value * resolution_component.evaluate( + energy - sample_component.center.value - offset_float + ) + + if isinstance(resolution_component, DeltaFunction): + return True, resolution_component.area.value * sample_component.evaluate( + energy - resolution_component.center.value - offset_float + ) + + # Gaussian + Gaussian --> Gaussian with width sqrt(w1^2 + w2^2) and area A1 * A2 + if isinstance(sample_component, Gaussian) and isinstance( + resolution_component, Gaussian + ): + width = np.sqrt( + sample_component.width.value**2 + resolution_component.width.value**2 + ) + area = sample_component.area.value * resolution_component.area.value + center = ( + sample_component.center.value + resolution_component.center.value + ) + offset_float + return True, _gaussian_eval(energy, center, width, area) + + # Lorentzian + Lorentzian --> Lorentzian with width w1 + w2 and area A1 * A2 + if isinstance(sample_component, Lorentzian) and isinstance( + resolution_component, Lorentzian + ): + width = sample_component.width.value + resolution_component.width.value + area = sample_component.area.value * resolution_component.area.value + center = ( + sample_component.center.value + resolution_component.center.value + ) + offset_float + return True, _lorentzian_eval(energy, center, width, area) + + # Gaussian + Lorentzian --> Voigt with area A1 * A2 + if ( + isinstance(sample_component, Gaussian) + and isinstance(resolution_component, Lorentzian) + ) or ( + isinstance(sample_component, Lorentzian) + and isinstance(resolution_component, Gaussian) + ): + if isinstance(sample_component, Gaussian): + gaussian, lorentzian = sample_component, resolution_component + else: + gaussian, lorentzian = resolution_component, sample_component + center = (gaussian.center.value + lorentzian.center.value) + offset_float + area = gaussian.area.value * lorentzian.area.value + return True, _voigt_eval( + energy, center, gaussian.width.value, lorentzian.width.value, area + ) + + return False, np.zeros_like(energy, dtype=float) + + +def _gaussian_eval( + energy: np.ndarray, center: float, width: float, area: float +) -> np.ndarray: + """ + Evaluate a Gaussian function. y = (area / (sqrt(2pi) * width)) * exp(-0.5 * ((x - center) / width)^2) + All checks are handled in the calling function. + + Args: + energy : np.ndarray + 1D array of energy values where the Gaussian is evaluated. + center : float + The center of the Gaussian. + width : float + The width (sigma) of the Gaussian. + area : float + The area under the Gaussian curve. + Returns: + np.ndarray + The evaluated Gaussian values at x. + """ + return ( + area + * 1 + / (np.sqrt(2 * np.pi) * width) + * np.exp(-0.5 * ((energy - center) / width) ** 2) + ) + + +def _lorentzian_eval( + energy: np.ndarray, center: float, width: float, area: float +) -> np.ndarray: + """ + Evaluate a Lorentzian function. y = (area * width / pi) / ((x - center)^2 + width^2). + All checks are handled in the calling function. + + Args: + energy : np.ndarray + 1D array of energy values where the Lorentzian is evaluated. + center : float + The center of the Lorentzian. + width : float + The width (HWHM) of the Lorentzian. + area : float + The area under the Lorentzian. + Returns: + np.ndarray + The evaluated Lorentzian values at x. + """ + return area * width / np.pi / ((energy - center) ** 2 + width**2) + + +def _voigt_eval( + energy: np.ndarray, center: float, g_width: float, l_width: float, area: float +) -> np.ndarray: + """ + Evaluate a Voigt profile function using scipy's voigt_profile. + Args: + energy : np.ndarray + 1D array of energy values where the Voigt profile is evaluated. + center : float + The center of the Voigt profile. + g_width : float + The Gaussian width (sigma) of the Voigt profile. + l_width : float + The Lorentzian width (HWHM) of the Voigt profile. + area : float + The area under the Voigt profile. + Returns: + np.ndarray + The evaluated Voigt profile values at x. + """ + + return area * voigt_profile(energy - center, g_width, l_width) + + +def _check_width_thresholds( + model: Union[SampleModel, ModelComponent], + span: float, + energy_step: float, + model_name: str, +) -> None: + """ + Helper function to check and warn if components are wide compared to the span of the data, or narrow compared to the spacing. + In both cases, the convolution accuracy may be compromised. + Args: + model : SampleModel or ModelComponent + The model to check. + energy_step : float + The bin spacing of the energy array. + span : float + The total span of the energy array. + model_name : str + A string indicating whether the model is a 'sample model' or 'resolution model' for warning messages. + returns: + None + warns: + UserWarning + If the component widths are not appropriate for the data span or bin spacing. + + """ + + # The thresholds are illustrated in performance_tests/convolution/convolution_width_thresholds.ipynb + LARGE_WIDTH_THRESHOLD = ( + 0.1 # Threshold for large widths compared to span - warn if width > 10% of span + ) + SMALL_WIDTH_THRESHOLD = ( + 1.0 # Threshold for small widths compared to bin spacing - warn if width < dx + ) + + # Handle SampleModel or ModelComponent + if isinstance(model, SampleModel): + components = model.components + else: + components = [model] # Treat single ModelComponent as a list + + for comp in components: + if hasattr(comp, "width"): + if comp.width.value > LARGE_WIDTH_THRESHOLD * span: + warnings.warn( + f"The width of the {model_name} component '{comp.name}' ({comp.width.value}) is large compared to the span of the input " + f"array ({span}). This may lead to inaccuracies in the convolution. Increase extension_factor to improve accuracy.", + UserWarning, + ) + if comp.width.value < SMALL_WIDTH_THRESHOLD * energy_step: + warnings.warn( + f"The width of the {model_name} component '{comp.name}' ({comp.width.value}) is small compared to the spacing of the input " + f"array ({energy_step}). This may lead to inaccuracies in the convolution. Increase upsample_factor to improve accuracy.", + UserWarning, + ) + + +def _find_delta_components( + model: Union[SampleModel, ModelComponent], +) -> List[DeltaFunction]: + """Return a list of DeltaFunction instances contained in `model`. + + Args: + model : SampleModel or ModelComponent + The model to search for DeltaFunction components. + Returns: + List[DeltaFunction] + A list of DeltaFunction components found in the model. + """ + if isinstance(model, DeltaFunction): + return [model] + if isinstance(model, SampleModel): + return [c for c in model.components if isinstance(c, DeltaFunction)] + return [] + + +def _calculate_delta_contributions( + energy: np.ndarray, + sample_model: Union[SampleModel, ModelComponent], + resolution_model: Union[SampleModel, ModelComponent], + offset_float: float, +) -> np.ndarray: + """ + Calculate the contributions of delta functions in the convolution. + Args: + energy : np.ndarray + 1D array of energy values where the convolution is evaluated. + sample_model : SampleModel or ModelComponent + The sample model to be convolved. + resolution_model : SampleModel or ModelComponent + The resolution model to convolve with. + offset_float : float + The offset to apply to the convolution. + Returns: + np.ndarray + The delta function contributions evaluated at energy. + + Raises: + ValueError + If both sample_model and resolution_model contain delta functions. + """ + delta_contributions = np.zeros_like(energy) + + # Add delta contributions on original grid + # collect deltas + sample_deltas = _find_delta_components(sample_model) + resolution_deltas = _find_delta_components(resolution_model) + + # error if both contain delta(s) + if sample_deltas and resolution_deltas: + raise ValueError( + "Both sample_model and resolution_model contain delta functions. " + "Their convolution is not defined." + ) + + # if sample has deltas, convolve each delta with the resolution_model + for delta in sample_deltas: + (_, delta_contribution) = _try_analytic_pair( + energy=energy, + sample_component=delta, + resolution_component=resolution_model, + offset_float=offset_float, + ) + delta_contributions += delta_contribution + + # if resolution has deltas, convolve each delta with the sample_model + for delta in resolution_deltas: + (_, delta_contribution) = _try_analytic_pair( + energy=energy, + sample_component=sample_model, + resolution_component=delta, + offset_float=offset_float, + ) + delta_contributions += delta_contribution + + return delta_contributions diff --git a/tests/performance_tests/convolution/convolution_width_thresholds.ipynb b/tests/performance_tests/convolution/convolution_width_thresholds.ipynb new file mode 100644 index 0000000..925b2ab --- /dev/null +++ b/tests/performance_tests/convolution/convolution_width_thresholds.ipynb @@ -0,0 +1,126 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "018fa173", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from easyscience.variable import Parameter\n", + "from scipy.special import voigt_profile\n", + "\n", + "from easydynamics.sample_model import DeltaFunction, Gaussian, Lorentzian, SampleModel, DampedHarmonicOscillator\n", + "from easydynamics.utils import convolution \n", + "\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf69a4b9", + "metadata": {}, + "outputs": [], + "source": [ + "# When the width of the Gaussian is >~20% of the span, numerical issues arise. We set the limit to 10% to be safe.\n", + "gaussian_widths=[30, 20, 10, 5]\n", + "sample_model=SampleModel(name='sample_model')\n", + "gaussian=Gaussian(name='Gaussian',width=30,area=1)\n", + "sample_model.add_component(gaussian)\n", + "\n", + "resolution_model = SampleModel(name='resolution_model')\n", + "resolution_gaussian=Gaussian(name='Resolution Gaussian',width=5,area=1.0)\n", + "resolution_model.add_component(resolution_gaussian)\n", + "x=np.linspace(-50, 50, 101)\n", + "\n", + "for gwidth in gaussian_widths:\n", + " sample_model['Gaussian'].width=gwidth\n", + " y_analytical = convolution(sample_model=sample_model,\n", + " resolution_model=resolution_model,\n", + " x=x,\n", + " )\n", + "\n", + " y_numerical = convolution(sample_model=sample_model,\n", + " resolution_model=resolution_model,\n", + " x=x,\n", + " method='numerical',\n", + " upsample_factor=0\n", + " )\n", + "\n", + " plt.plot(x, y_analytical, label='Analytical, width = {}'.format(gwidth), color = plt.cm.viridis(gaussian_widths.index(gwidth)/len(gaussian_widths)))\n", + " plt.plot(x, y_numerical, label='Numerical, width = {}'.format(gwidth), color = plt.cm.viridis(gaussian_widths.index(gwidth)/len(gaussian_widths)), linestyle='--')\n", + " plt.xlabel('Energy (meV)')\n", + " plt.ylabel('Intensity (arb. units)')\n", + " plt.title('Convolution of Sample Model with Resolution Model with various widths')\n", + "\n", + "plt.legend()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cb777e0", + "metadata": {}, + "outputs": [], + "source": [ + "# When the width of the Gaussian is <~50% of the bin spacing, numerical issues arise. We set the limit to 100% to be safe.\n", + "gaussian_widths=[5, 2, 1, 0.5, 0.25]\n", + "gaussian_centers=[-50, -25, 0, 25, 50]\n", + "sample_model=SampleModel(name='sample_model')\n", + "gaussian=Gaussian(name='Gaussian',width=0.5,area=1)\n", + "sample_model.add_component(gaussian)\n", + "\n", + "resolution_model = SampleModel(name='resolution_model')\n", + "resolution_gaussian=Gaussian(name='Resolution Gaussian',width=10,area=1.0)\n", + "resolution_model.add_component(resolution_gaussian)\n", + "x=np.linspace(-100, 100, 201)\n", + "\n", + "for gwidth, gcenter in zip(gaussian_widths, gaussian_centers):\n", + " sample_model['Gaussian'].width=gwidth\n", + " sample_model['Gaussian'].center=gcenter\n", + " y_analytical = convolution(sample_model=sample_model,\n", + " resolution_model=resolution_model,\n", + " x=x,\n", + " )\n", + "\n", + " y_numerical = convolution(sample_model=sample_model,\n", + " resolution_model=resolution_model,\n", + " x=x,\n", + " method='numerical',\n", + " upsample_factor=0\n", + " )\n", + "\n", + " plt.plot(x, y_analytical, label='Analytical, width = {}'.format(gwidth), color = plt.cm.viridis(gaussian_widths.index(gwidth)/len(gaussian_widths)))\n", + " plt.plot(x, y_numerical, label='Numerical, width = {}'.format(gwidth), color = plt.cm.viridis(gaussian_widths.index(gwidth)/len(gaussian_widths)), linestyle='--')\n", + " plt.xlabel('Energy (meV)')\n", + " plt.ylabel('Intensity (arb. units)')\n", + " plt.title('Convolution of Sample Model with Resolution Model with various widths')\n", + "\n", + "plt.legend()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "newdynamics", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/unit_tests/sample_model/test_sample_model.py b/tests/unit_tests/sample_model/test_sample_model.py new file mode 100644 index 0000000..19d1042 --- /dev/null +++ b/tests/unit_tests/sample_model/test_sample_model.py @@ -0,0 +1,318 @@ +from copy import copy + +import numpy as np +import pytest +from easyscience.variable import Parameter +from scipy.integrate import simpson + +from easydynamics.sample_model import Gaussian, Lorentzian, Polynomial, SampleModel + + +class TestSampleModel: + @pytest.fixture + def sample_model(self): + model = SampleModel(name="TestSampleModel") + component1 = Gaussian( + name="TestGaussian1", area=1.0, center=0.0, width=1.0, unit="meV" + ) + component2 = Lorentzian( + name="TestLorentzian1", area=2.0, center=1.0, width=0.5, unit="meV" + ) + model.add_component(component1) + model.add_component(component2) + return model + + def test_init(self): + # WHEN THEN + sample_model = SampleModel(name="InitModel") + + # EXPECT + assert sample_model.name == "InitModel" + assert len(sample_model.components) == 0 + + def test_initialization_with_components(self): + # WHEN THEN + component1 = Gaussian( + name="InitGaussian", area=1.0, center=0.0, width=1.0, unit="meV" + ) + component2 = Lorentzian( + name="InitLorentzian", area=2.0, center=1.0, width=0.5, unit="meV" + ) + sample_model = SampleModel( + name="InitModelWithComponents", data=[component1, component2] + ) + + # EXPECT + assert sample_model.name == "InitModelWithComponents" + assert len(sample_model.components) == 2 + assert sample_model["InitGaussian"] is component1 + assert sample_model["InitLorentzian"] is component2 + + # ───── Component Management ───── + + def test_add_component(self, sample_model): + # WHEN + component = Gaussian( + name="TestComponent", area=1.0, center=0.0, width=1.0, unit="meV" + ) + # THEN + sample_model.add_component(component) + # EXPECT + assert sample_model["TestComponent"] is component + + def test_add_duplicate_component_raises(self, sample_model): + # WHEN THEN + component = Gaussian( + name="TestGaussian1", area=1.0, center=0.0, width=1.0, unit="meV" + ) + # EXPECT + with pytest.raises(ValueError, match="already exists"): + sample_model.add_component(component) + + def test_add_invalid_component_raises(self, sample_model): + # WHEN THEN EXPECT + with pytest.raises( + TypeError, match="component must be an instance of ModelComponent." + ): + sample_model.add_component("NotAComponent") + + def test_remove_component(self, sample_model): + # WHEN THEN + sample_model.remove_component("TestGaussian1") + # EXPECT + assert "TestGaussian1" not in sample_model.components + + def test_remove_nonexistent_component_raises(self, sample_model): + # WHEN THEN EXPECT + with pytest.raises( + KeyError, match="No component named 'NonExistentComponent' exists" + ): + sample_model.remove_component("NonExistentComponent") + + def test_getitem(self, sample_model): + # WHEN + component = Gaussian( + name="TestComponent", area=1.0, center=0.0, width=1.0, unit="meV" + ) + # THEN + sample_model.add_component(component) + # EXPECT + assert sample_model["TestComponent"] is component + + def test_list_component_names(self, sample_model): + # WHEN THEN + components = sample_model.list_component_names() + # EXPECT + assert len(components) == 2 + assert components[0] == "TestGaussian1" + assert components[1] == "TestLorentzian1" + + def test_clear_components(self, sample_model): + # WHEN THEN + sample_model.clear_components() + # EXPECT + assert len(sample_model.components) == 0 + + def test_convert_unit(self, sample_model): + # WHEN THEN + sample_model.convert_unit("eV") + # EXPECT + for component in list(sample_model): + assert component.unit == "eV" + + def test_evaluate(self, sample_model): + # WHEN + x = np.linspace(-5, 5, 100) + result = sample_model.evaluate(x) + # EXPECT + expected_result = sample_model["TestGaussian1"].evaluate(x) + sample_model[ + "TestLorentzian1" + ].evaluate(x) + np.testing.assert_allclose(result, expected_result, rtol=1e-5) + + def test_evaluate_no_components_raises(self): + # WHEN THEN + sample_model = SampleModel(name="EmptyModel") + x = np.linspace(-5, 5, 100) + # EXPECT + with pytest.raises(ValueError, match="No components in the model to evaluate."): + sample_model.evaluate(x) + + def test_evaluate_component(self, sample_model): + # WHEN THEN + x = np.linspace(-5, 5, 100) + result1 = sample_model.evaluate_component(x, "TestGaussian1") + result2 = sample_model.evaluate_component(x, "TestLorentzian1") + + # EXPECT + expected_result1 = sample_model["TestGaussian1"].evaluate(x) + expected_result2 = sample_model["TestLorentzian1"].evaluate(x) + np.testing.assert_allclose(result1, expected_result1, rtol=1e-5) + np.testing.assert_allclose(result2, expected_result2, rtol=1e-5) + + def test_evaluate_nonexistent_component_raises(self, sample_model): + # WHEN + x = np.linspace(-5, 5, 100) + + # THEN EXPECT + with pytest.raises( + KeyError, match="No component named 'NonExistentComponent' exists" + ): + sample_model.evaluate_component(x, "NonExistentComponent") + + def test_evaluate_component_no_components_raises(self): + # WHEN THEN + sample_model = SampleModel(name="EmptyModel") + x = np.linspace(-5, 5, 100) + # EXPECT + with pytest.raises(ValueError, match="No components in the model to evaluate."): + sample_model.evaluate_component(x, "AnyComponent") + + def test_evaluate_component_invalid_name_type_raises(self, sample_model): + # WHEN + x = np.linspace(-5, 5, 100) + + # THEN EXPECT + with pytest.raises( + TypeError, + match="Component name must be a string, got instead.", + ): + sample_model.evaluate_component(x, 123) + + # ───── Utilities ───── + + def test_normalize_area(self, sample_model): + # WHEN THEN + sample_model.normalize_area() + # EXPECT + x = np.linspace(-10000, 10000, 1000000) # Lorentzians have long tails + result = sample_model.evaluate(x) + numerical_area = simpson(result, x) + assert np.isclose(numerical_area, 1.0, rtol=1e-4) + + def test_normalize_area_no_components_raises(self): + # WHEN THEN + sample_model = SampleModel(name="EmptyModel") + # EXPECT + with pytest.raises( + ValueError, match="No components in the model to normalize." + ): + sample_model.normalize_area() + + @pytest.mark.parametrize( + "area_value", + [np.nan, 0.0, np.inf], + ids=["NaN area", "Zero area", "Infinite area"], + ) + def test_normalize_area_not_finite_area_raises(self, sample_model, area_value): + # WHEN THEN + sample_model["TestGaussian1"].area = area_value + sample_model["TestLorentzian1"].area = area_value + + # EXPECT + with pytest.raises(ValueError, match="cannot normalize."): + sample_model.normalize_area() + + def test_normalize_area_non_area_component_warns(self, sample_model): + # WHEN + component1 = Polynomial( + name="TestPolynomial", coefficients=[1, 2, 3], unit="meV" + ) + sample_model.add_component(component1) + + # THEN EXPECT + with pytest.warns(UserWarning, match="does not have an 'area' "): + sample_model.normalize_area() + + def test_get_parameters(self, sample_model): + # WHEN THEN + parameters = sample_model.get_parameters() + # EXPECT + assert len(parameters) == 6 + + expected_names = { + "TestGaussian1 area", + "TestGaussian1 center", + "TestGaussian1 width", + "TestLorentzian1 area", + "TestLorentzian1 center", + "TestLorentzian1 width", + } + actual_names = {param.name for param in parameters} + assert actual_names == expected_names + assert all(isinstance(param, Parameter) for param in parameters) + + def test_get_parameters_no_components(self): + sample_model = SampleModel(name="EmptyModel") + # WHEN THEN + parameters = sample_model.get_parameters() + # EXPECT + assert len(parameters) == 0 + + def test_get_fit_parameters(self, sample_model): + # WHEN + + # Fix one parameter and make another dependent + sample_model["TestGaussian1"].area.fixed = True + sample_model["TestLorentzian1"].width.make_dependent_on( + "comp1_width", + {"comp1_width": sample_model["TestGaussian1"].width}, + ) + + # THEN + fit_parameters = sample_model.get_fit_parameters() + + # EXPECT + assert len(fit_parameters) == 4 + + expected_names = { + "TestGaussian1 center", + "TestGaussian1 width", + "TestLorentzian1 area", + "TestLorentzian1 center", + } + actual_names = {param.name for param in fit_parameters} + assert actual_names == expected_names + assert all(isinstance(param, Parameter) for param in fit_parameters) + + def test_fix_and_free_all_parameters(self, sample_model): + # WHEN THEN + sample_model.fix_all_parameters() + + # EXPECT + for param in sample_model.get_parameters(): + assert param.fixed is True + + # WHEN + sample_model.free_all_parameters() + + # THEN + for param in sample_model.get_parameters(): + assert param.fixed is False + + def test_repr_contains_name_and_components(self, sample_model): + # WHEN THEN + rep = repr(sample_model) + # EXPECT + assert "SampleModel" in rep + assert "TestGaussian" in rep + + def test_copy(self, sample_model): + # WHEN THEN + sample_model.temperature = 300 + model_copy = copy(sample_model) + # EXPECT + assert model_copy is not sample_model + assert model_copy.name == sample_model.name + assert len(list(model_copy)) == len(list(sample_model)) + for comp in list(sample_model): + copied_comp = model_copy[comp.name] + assert copied_comp is not comp + assert copied_comp.name == comp.name + for param_orig, param_copy in zip( + comp.get_parameters(), copied_comp.get_parameters() + ): + assert param_copy is not param_orig + assert param_copy.name == param_orig.name + assert param_copy.value == param_orig.value + assert param_copy.fixed == param_orig.fixed diff --git a/tests/unit_tests/utils/test_convolution.py b/tests/unit_tests/utils/test_convolution.py new file mode 100644 index 0000000..70fb820 --- /dev/null +++ b/tests/unit_tests/utils/test_convolution.py @@ -0,0 +1,852 @@ +import numpy as np +import pytest +from easyscience.variable import Parameter +from scipy.signal import fftconvolve +from scipy.special import voigt_profile + +from easydynamics.sample_model import ( + DampedHarmonicOscillator, + DeltaFunction, + Gaussian, + Lorentzian, + SampleModel, +) +from easydynamics.utils import convolution +from easydynamics.utils.detailed_balance import ( + _detailed_balance_factor as detailed_balance_factor, +) + +# Numerical convolutions are not very accurate +NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE = 1e-6 +NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE = 1e-5 + + +class TestConvolution: + @pytest.fixture + def sample_model(self): + test_sample_model = SampleModel(name="TestSampleModel") + test_sample_model.add_component(Gaussian(center=0.1, width=0.3, area=2.0)) + return test_sample_model + + @pytest.fixture + def resolution_model(self): + test_resolution_model = SampleModel(name="TestResolutionModel") + test_resolution_model.add_component(Gaussian(center=0.2, width=0.4, area=3.0)) + return test_resolution_model + + @pytest.fixture + def gaussian_component(self): + return Gaussian(center=0.1, width=0.3, area=2.0) + + @pytest.fixture + def other_gaussian_component(self): + return Gaussian(name="other Gaussian", center=0.2, width=0.4, area=3.0) + + @pytest.fixture + def lorentzian_component(self): + return Lorentzian(center=0.1, width=0.3, area=2.0) + + @pytest.fixture + def other_lorentzian_component(self): + return Lorentzian(center=0.2, width=0.4, area=3.0) + + @pytest.fixture + def energy(self): + return np.linspace(-50, 50, 50001) + + # Test convolution of components + @pytest.mark.parametrize( + "offset_obj, expected_shift", + [ + (None, 0.0), + (0.4, 0.4), + (Parameter("off", 0.4), 0.4), + ], + ids=["none", "float", "parameter"], + ) + @pytest.mark.parametrize( + "method", ["analytical", "numerical"], ids=["analytical", "numerical"] + ) + def test_components_gauss_gauss( + self, + energy, + gaussian_component, + other_gaussian_component, + offset_obj, + expected_shift, + method, + ): + "Test convolution of Gaussian sample and Gaussian resolution components without SampleModel." + "Test with different offset types and methods." + # WHEN + sample_gauss = gaussian_component + resolution_gauss = other_gaussian_component + + # THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample_gauss, + resolution_model=resolution_gauss, + offset=offset_obj, + method=method, + ) + + # EXPECT + # Convolution of two Gaussians is another Gaussian with width = sqrt(w1^2 + w2^2) + expected_width = np.sqrt( + sample_gauss.width.value**2 + resolution_gauss.width.value**2 + ) + expected_area = sample_gauss.area.value * resolution_gauss.area.value + expected_center = ( + sample_gauss.center.value + resolution_gauss.center.value + expected_shift + ) + expected_result = ( + expected_area + * np.exp(-0.5 * ((energy - expected_center) / expected_width) ** 2) + / (np.sqrt(2 * np.pi) * expected_width) + ) + + np.testing.assert_allclose(calculated_convolution, expected_result, atol=1e-10) + + @pytest.mark.parametrize( + "offset_obj, expected_shift", + [ + (None, 0.0), + (0.4, 0.4), + (Parameter("off", 0.4), 0.4), + ], + ids=["none", "float", "parameter"], + ) + @pytest.mark.parametrize("method", ["auto", "numerical"], ids=["auto", "numerical"]) + def test_components_DHO_gauss( + self, energy, gaussian_component, offset_obj, expected_shift, method + ): + "Test convolution of DHO sample and Gaussian resolution components without SampleModel." + "Test with different offset types and methods." + # WHEN + sample_dho = DampedHarmonicOscillator(center=1.5, width=0.3, area=2) + resolution_gauss = gaussian_component + + # THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample_dho, + resolution_model=resolution_gauss, + offset=offset_obj, + method=method, + ) + + # EXPECT + # no simple analytical form, so compute expected result via direct convolution + sample_values = sample_dho.evaluate(energy - expected_shift) + resolution_values = resolution_gauss.evaluate(energy) + expected_result = fftconvolve(sample_values, resolution_values, mode="same") + expected_result *= energy[1] - energy[0] # normalize + + np.testing.assert_allclose(calculated_convolution, expected_result, atol=1e-10) + + @pytest.mark.parametrize( + "offset_obj, expected_shift", + [ + (None, 0.0), + (0.4, 0.4), + (Parameter("off", 0.4), 0.4), + ], + ids=["none", "float", "parameter"], + ) + @pytest.mark.parametrize( + "method", ["analytical", "numerical"], ids=["analytical", "numerical"] + ) + def test_components_lorentzian_lorentzian( + self, + energy, + lorentzian_component, + other_lorentzian_component, + offset_obj, + expected_shift, + method, + ): + "Test convolution of Lorentzian sample and Lorentzian resolution components without SampleModel." + "Test with different offset types and methods." + # WHEN + sample_lorentzian = lorentzian_component + resolution_lorentzian = other_lorentzian_component + + # THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample_lorentzian, + resolution_model=resolution_lorentzian, + offset=offset_obj, + method=method, + upsample_factor=5, + ) + + # EXPECT + # Convolution of two Lorentzians is another Lorentzian with width = w1 + w2 + expected_width = ( + sample_lorentzian.width.value + resolution_lorentzian.width.value + ) + expected_area = sample_lorentzian.area.value * resolution_lorentzian.area.value + expected_center = ( + sample_lorentzian.center.value + + resolution_lorentzian.center.value + + expected_shift + ) + expected_result = ( + expected_area + * expected_width + / np.pi + / ((energy - expected_center) ** 2 + expected_width**2) + ) + + np.testing.assert_allclose( + calculated_convolution, + expected_result, + atol=NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE, + rtol=NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE, + ) + + @pytest.mark.parametrize( + "offset_obj, expected_shift", + [ + (None, 0.0), + (0.4, 0.4), + (Parameter("off", 0.4), 0.4), + ], + ids=["none", "float", "parameter"], + ) + @pytest.mark.parametrize( + "method", ["analytical", "numerical"], ids=["analytical", "numerical"] + ) + @pytest.mark.parametrize( + "sample_is_gauss", + [True, False], + ids=["gauss_sample__lorentz_resolution", "lorentz_sample__gauss_resolution"], + ) + def test_components_gauss_lorentzian( + self, + energy, + gaussian_component, + lorentzian_component, + offset_obj, + expected_shift, + method, + sample_is_gauss, + ): + "Test convolution of Gaussian and Lorentzian components without SampleModel." + "Test with different offset types and methods." + # WHEN + if sample_is_gauss: + sample = gaussian_component + resolution = lorentzian_component + else: + sample = lorentzian_component + resolution = gaussian_component + + # THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample, + resolution_model=resolution, + offset=offset_obj, + method=method, + upsample_factor=5, + ) + + # EXPECT + expected_center = sample.center.value + resolution.center.value + expected_shift + expected_area = sample.area.value * resolution.area.value + + gaussian_width = ( + sample.width.value if sample_is_gauss else resolution.width.value + ) + lorentzian_width = ( + resolution.width.value if sample_is_gauss else sample.width.value + ) + + expected_result = expected_area * voigt_profile( + energy - expected_center, + gaussian_width, + lorentzian_width, + ) + + np.testing.assert_allclose( + calculated_convolution, + expected_result, + atol=NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE, + rtol=NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE, + ) + + @pytest.mark.parametrize( + "offset_obj, expected_shift", + [ + (None, 0.0), + (0.4, 0.4), + (Parameter("off", 0.4), 0.4), + ], + ids=["none", "float", "parameter"], + ) + @pytest.mark.parametrize( + "method", ["analytical", "numerical"], ids=["analytical", "numerical"] + ) + @pytest.mark.parametrize( + "sample_is_gauss", + [True, False], + ids=["gauss_sample__delta_resolution", "delta_sample__gauss_resolution"], + ) + def test_components_delta_gauss( + self, + energy, + gaussian_component, + offset_obj, + expected_shift, + method, + sample_is_gauss, + ): + "Test convolution of Delta function sample and Gaussian resolution components without SampleModel." + "Test with different offset types and methods." + # WHEN + if sample_is_gauss: + sample = gaussian_component + resolution = DeltaFunction(name="Delta", center=0.1, area=2) + else: + sample = DeltaFunction(name="Delta", center=0.1, area=2) + resolution = gaussian_component + + # THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample, + resolution_model=resolution, + offset=offset_obj, + method=method, + ) + + # EXPECT + expected_center = sample.center.value + resolution.center.value + expected_shift + expected_area = sample.area.value * resolution.area.value + width = sample.width.value if sample_is_gauss else resolution.width.value + expected_result = ( + expected_area + * np.exp(-0.5 * ((energy - expected_center) / width) ** 2) + / (np.sqrt(2 * np.pi) * width) + ) + + np.testing.assert_allclose(calculated_convolution, expected_result, atol=1e-10) + + # Test convolution of SampleModel + @pytest.mark.parametrize( + "offset_obj, expected_shift", + [ + (None, 0.0), + (0.4, 0.4), + (Parameter("off", 0.4), 0.4), + ], + ids=["none", "float", "parameter"], + ) + @pytest.mark.parametrize( + "method", ["analytical", "numerical"], ids=["analytical", "numerical"] + ) + def test_model_gauss_gauss_resolution_gauss( + self, + energy, + sample_model, + resolution_model, + offset_obj, + expected_shift, + method, + ): + "Test convolution of Gaussian sample components in SampleModel and Gaussian resolution components in SampleModel." + "Test with different offset types and methods." + + # WHEN + sample_G2 = Gaussian(name="another Gaussian", center=0.3, width=0.5, area=4) + sample_model.add_component(sample_G2) + + # THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + offset=offset_obj, + method=method, + ) + + # EXPECT + sample_G1 = sample_model["Gaussian"] + resolution_G1 = resolution_model["Gaussian"] + expected_width1 = np.sqrt( + sample_G1.width.value**2 + resolution_G1.width.value**2 + ) + expected_width2 = np.sqrt( + sample_G2.width.value**2 + resolution_G1.width.value**2 + ) + expected_area1 = sample_G1.area.value * resolution_G1.area.value + expected_area2 = sample_G2.area.value * resolution_G1.area.value + expected_center1 = ( + sample_G1.center.value + resolution_G1.center.value + expected_shift + ) + expected_center2 = ( + sample_G2.center.value + resolution_G1.center.value + expected_shift + ) + + expected_result = expected_area1 * np.exp( + -0.5 * ((energy - expected_center1) / expected_width1) ** 2 + ) / (np.sqrt(2 * np.pi) * expected_width1) + expected_area2 * np.exp( + -0.5 * ((energy - expected_center2) / expected_width2) ** 2 + ) / (np.sqrt(2 * np.pi) * expected_width2) + np.testing.assert_allclose(calculated_convolution, expected_result, atol=1e-10) + + @pytest.mark.parametrize( + "offset_obj, expected_shift", + [ + (None, 0.0), + (0.4, 0.4), + (Parameter("off", 0.4), 0.4), + ], + ids=["none", "float", "parameter"], + ) + @pytest.mark.parametrize( + "method", ["analytical", "numerical"], ids=["analytical", "numerical"] + ) + def test_model_lorentzian_delta_resolution_gauss( + self, + energy, + method, + lorentzian_component, + resolution_model, + offset_obj, + expected_shift, + ): + "Test convolution of Lorentzian and Delta function sample components in SampleModel and Gaussian resolution components in SampleModel." + " Result is a combination of Voigt profile and Gaussian." + # WHEN + + sample = SampleModel(name="SampleModel") + sample.add_component(lorentzian_component) + sample_delta = DeltaFunction(center=0.5, area=4, name="SampleDelta") + sample.add_component(sample_delta) + + # THEN + energy = np.linspace(-10, 10, 20001) + calculated_convolution = convolution( + energy=energy, + sample_model=sample, + resolution_model=resolution_model, + offset=offset_obj, + method=method, + upsample_factor=5, + ) + + # EXPECT: Combine Gaussian, Lorentzian, and Delta functions contributions + # + gaussian_component = resolution_model["Gaussian"] + + expected_voigt_area = ( + lorentzian_component.area.value * gaussian_component.area.value + ) + expected_voigt_center = ( + lorentzian_component.center.value + + gaussian_component.center.value + + expected_shift + ) + expected_voigt = expected_voigt_area * voigt_profile( + energy - expected_voigt_center, + gaussian_component.width.value, + lorentzian_component.width.value, + ) + expected_gauss_area = sample_delta.area.value * gaussian_component.area.value + expected_gauss_center = ( + sample_delta.center.value + gaussian_component.center.value + expected_shift + ) + expected_gauss_width = gaussian_component.width.value + expected_gauss = ( + expected_gauss_area + * np.exp( + -0.5 * ((energy - (expected_gauss_center)) / expected_gauss_width) ** 2 + ) + / (np.sqrt(2 * np.pi) * expected_gauss_width) + ) + expected_result = expected_voigt + expected_gauss + np.testing.assert_allclose( + calculated_convolution, + expected_result, + atol=NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE, + rtol=NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE, + ) + + def test_numerical_convolve_with_temperature( + self, energy, sample_model, resolution_model + ): + "Test numerical convolution with detailed balance correction." + # WHEN + temperature = 300.0 # Kelvin + + # THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + method="numerical", + upsample_factor=5, + temperature=temperature, + ) + + sample_with_db = sample_model.evaluate(energy) * detailed_balance_factor( + energy=energy, temperature=temperature + ) + resolution = resolution_model.evaluate(energy) + + expected_convolution = fftconvolve(sample_with_db, resolution, mode="same") + expected_convolution *= [energy[1] - energy[0]] # normalize + + np.testing.assert_allclose( + calculated_convolution, + expected_convolution, + atol=NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE, + rtol=NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE, + ) + + @pytest.mark.parametrize( + "x", + [ + np.linspace(-10, 10, 5001), # Odd length + np.linspace(-10, 10, 5000), # Even length + ], + ids=["odd_length", "even_length"], + ) + def test_numerical_convolve_x_length_even_and_odd( + self, x, sample_model, resolution_model + ): + "Test numerical convolution with both even and odd length x arrays. With even length the FFT shifts the signal by half a bin." + + # WHEN THEN + calculated_convolution = convolution( + energy=x, + sample_model=sample_model, + resolution_model=resolution_model, + method="numerical", + upsample_factor=0, + ) + + # EXPECT + expected_convolution = convolution( + energy=x, + sample_model=sample_model, + resolution_model=resolution_model, + method="analytical", + upsample_factor=0, + ) + + np.testing.assert_allclose( + calculated_convolution, expected_convolution, atol=1e-10 + ) + + @pytest.mark.parametrize( + "upsample_factor", + [0, 2, 5, 10], + ids=["no_upsample", "upsample_2", "upsample_5", "upsample_10"], + ) + def test_numerical_convolve_upsample_factor( + self, energy, upsample_factor, sample_model, resolution_model + ): + "Test numerical convolution with different upsample factors." + # WHEN THEN + calculated_convolution = convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + method="numerical", + upsample_factor=upsample_factor, + ) + + # EXPECT + expected_convolution = convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + method="analytical", + upsample_factor=0, + ) + + np.testing.assert_allclose( + calculated_convolution, + expected_convolution, + atol=NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE, + rtol=NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE, + ) + + @pytest.mark.parametrize( + "x", + [np.linspace(-5, 15, 20000), np.linspace(5, 15, 20000)], + ids=["asymmetric", "only_positive"], + ) + @pytest.mark.parametrize( + "upsample_factor", [0, 2, 5], ids=["no_upsample", "upsample_2", "upsample_5"] + ) + def test_numerical_convolve_x_not_symmetric( + self, x, upsample_factor, resolution_model + ): + "Test numerical convolution with asymmetric and only positive x arrays." + # WHEN + sample_model = SampleModel(name="SampleModel") + sample_model.add_component(Gaussian(center=9, width=0.3, area=2)) + + # THEN + calculated_convolution = convolution( + energy=x, + sample_model=sample_model, + resolution_model=resolution_model, + method="numerical", + upsample_factor=upsample_factor, + ) + + # EXPECT + expected_convolution = convolution( + energy=x, + sample_model=sample_model, + resolution_model=resolution_model, + method="analytical", + ) + + np.testing.assert_allclose( + calculated_convolution, + expected_convolution, + atol=NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE, + rtol=NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE, + ) + + def test_numerical_convolve_x_not_uniform(self, sample_model, resolution_model): + "Test numerical convolution with non-uniform x arrays." + # WHEN + x_1 = np.linspace(-2, 0, 1000) + x_2 = np.linspace(0.001, 2, 2000) + x_non_uniform = np.concatenate([x_1, x_2]) + + # THEN + calculated_convolution = convolution( + energy=x_non_uniform, + sample_model=sample_model, + resolution_model=resolution_model, + method="numerical", + upsample_factor=5, + ) + + # EXPECT + expected_convolution = convolution( + energy=x_non_uniform, + sample_model=sample_model, + resolution_model=resolution_model, + method="analytical", + ) + + np.testing.assert_allclose( + calculated_convolution, + expected_convolution, + atol=NUMERICAL_CONVOLUTION_ABSOLUTE_TOLERANCE, + rtol=NUMERICAL_CONVOLUTION_RELATIVE_TOLERANCE, + ) + + # Test error handling + def test_analytical_convolution_fails_with_detailed_balance( + self, energy, sample_model, resolution_model + ): + # WHEN + temperature = 300.0 + # THEN EXPECT + with pytest.raises( + ValueError, + match="Analytical convolution is not supported with detailed balance.", + ): + convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + method="analytical", + temperature=temperature, + ) + + def test_convolution_only_accepts_auto_analytical_and_numerical_methods( + self, energy, sample_model, resolution_model + ): + # WHEN THEN EXPECT + with pytest.raises( + ValueError, + match="Unknown convolution method: unknown_method. Choose from 'auto', 'analytical', or 'numerical'.", + ): + convolution( + energy=energy, + sample_model=sample_model, + resolution_model=resolution_model, + method="unknown_method", + ) + + def test_energy_must_be_1d_finite_array(self, sample_model, resolution_model): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match="`energy` must be a 1D finite array."): + convolution( + energy=np.array([[1, 2], [3, 4]]), + sample_model=sample_model, + resolution_model=resolution_model, + ) + + with pytest.raises(ValueError, match="`energy` must be a 1D finite array."): + convolution( + energy=np.array([1, 2, np.nan]), + sample_model=sample_model, + resolution_model=resolution_model, + ) + + with pytest.raises(ValueError, match="`energy` must be a 1D finite array."): + convolution( + energy=np.array([1, 2, np.inf]), + sample_model=sample_model, + resolution_model=resolution_model, + ) + + def test_numerical_convolve_requires_uniform_grid_if_no_upsample( + self, sample_model, resolution_model + ): + # WHEN + x = np.array([0, 1, 2, 4, 5]) # Non-uniform grid + + # THEN EXPECT + with pytest.raises( + ValueError, + match="Input array `energy` must be uniformly spaced if upsample_factor = 0.", + ): + convolution( + energy=x, + sample_model=sample_model, + resolution_model=resolution_model, + method="numerical", + upsample_factor=0, + ) + + def test_sample_model_must_have_components(self, resolution_model): + # WHEN + sample_model = SampleModel(name="SampleModel") + + # THEN EXPECT + with pytest.raises( + ValueError, match="SampleModel must have at least one component." + ): + convolution( + energy=np.array([0, 1, 2]), + sample_model=sample_model, + resolution_model=resolution_model, + ) + + def test_resolution_model_must_have_components(self, sample_model): + # WHEN + resolution_model = SampleModel(name="ResolutionModel") + + # THEN EXPECT + with pytest.raises( + ValueError, match="ResolutionModel must have at least one component." + ): + convolution( + energy=np.array([0, 1, 2]), + sample_model=sample_model, + resolution_model=resolution_model, + ) + + def test_numerical_convolution_wide_sample_peak_gives_warning( + self, resolution_model + ): + # WHEN + x = np.linspace(-2, 2, 20001) + + sample_gauss = Gaussian(center=0.1, width=1.9, area=2, name="SampleGauss") + sample = SampleModel(name="SampleModel") + sample.add_component(sample_gauss) + + # #THEN EXPECT + with pytest.warns( + UserWarning, + match=r"The width of the sample model component ", + ): + convolution( + energy=x, + sample_model=sample, + resolution_model=resolution_model, + method="numerical", + upsample_factor=0, + ) + + def test_numerical_convolution_wide_resolution_peak_gives_warning( + self, sample_model + ): + # WHEN + x = np.linspace(-2, 2, 20001) + + resolution_gauss = Gaussian( + center=0.3, width=1.9, area=4, name="ResolutionGauss" + ) + + resolution = SampleModel(name="ResolutionModel") + resolution.add_component(resolution_gauss) + + # #THEN EXPECT + with pytest.warns( + UserWarning, + match=r"The width of the resolution model component 'ResolutionGauss' \(1.9\) is large", + ): + convolution( + energy=x, + sample_model=sample_model, + resolution_model=resolution, + method="numerical", + upsample_factor=0, + ) + + def test_numerical_convolution_narrow_sample_peak_gives_warning( + self, resolution_model + ): + # WHEN + x = np.linspace(-2, 2, 201) + + sample_gauss1 = Gaussian(center=0.1, width=1e-3, area=2, name="SampleGauss") + + sample = SampleModel(name="SampleModel") + sample.add_component(sample_gauss1) + + # #THEN EXPECT + with pytest.warns( + UserWarning, + match=r"The width of the sample model component 'SampleGauss' \(0.001\) is small", + ): + convolution( + energy=x, + sample_model=sample, + resolution_model=resolution_model, + method="numerical", + upsample_factor=0, + ) + + def test_numerical_convolution_narrow_resolution_peak_gives_warning( + self, sample_model + ): + # WHEN + x = np.linspace(-2, 2, 201) + + resolution_gauss = Gaussian( + center=0.3, width=1e-3, area=4, name="ResolutionGauss" + ) + + resolution = SampleModel(name="ResolutionModel") + resolution.add_component(resolution_gauss) + + # #THEN EXPECT + with pytest.warns( + UserWarning, + match=r"The width of the resolution model component 'ResolutionGauss' \(0.001\) is small", + ): + convolution( + energy=x, + sample_model=sample_model, + resolution_model=resolution, + method="numerical", + upsample_factor=0, + )