{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial for Stacked Trapezoid model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This tutorial shows how to do simulations for stacked trapezoid model and how to fit the experimental data using this model." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from cdsaxs.simulations.stacked_trapezoid import StackedTrapezoidSimulation\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Prepare the data" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "pitch = 100 #nm distance between two lines\n", "qzs = np.linspace(-0.5, 0.5, 121)\n", "qxs = 2 * np.pi / pitch * np.ones_like(qzs)\n", "\n", "#parameters\n", "dwx = 0.1\n", "dwz = 0.1\n", "i0 = 10\n", "bkg = 0.1\n", "y1 = 0.\n", "height = 20. #same as [20., 20.]\n", "bot_cd = 40.\n", "swa = [80., 80.0]\n", "\n", "langle = np.deg2rad(np.asarray(swa))\n", "rangle = np.deg2rad(np.asarray(swa))\n", "\n", "#simulation data\n", "i_params = {'heights': height,\n", " 'langles': langle,\n", " 'rangles': rangle,\n", " 'y_start': y1,\n", " 'bot_cd': bot_cd,\n", " 'dwx': dwx,\n", " 'dwz': dwz,\n", " 'i0': i0,\n", " 'bkg_cste': bkg,\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Description of the parameters:\n", "\n", "dwx, dwz, i0 and bkg are the parameters necessary to calculate the Debye-waller factor which accounts for the real world imperfections in the model. You can read more about it [here](https://en.wikipedia.org/wiki/Debye%E2%80%93Waller_factor).\n", "\n", "The other parameters are the geometrical parameters of the model. The parameters are as follows:\n", "\n", "*y_start* - starting y-coordinate of the base of the nano structure\n", "\n", "*bot_cd* - bottom width (CD for critical dimension) of the nano structure\n", "\n", "*heights* - height of each individual trapezoid can be a single value or a list of values \n", "in case of list of values, each value corresponds to the height of each individual trapezoid.\n", "\n", "*langles* - list of all the left bottom angles of each individual trapezoid. The dictionary which is passed should have it in radians.\n", "\n", "*rangles* - list of all the right bottom angles of each individual trapezoid. The dictionary which is passed should have it in radians.\n", "\n", "*weight* - weight of each individual trapezoid to account for the fact that they could be made of different materials. here we will assume that all the trapezoids are made of the same material hence the weight is not necessary.\n", "\n", "Note: In symmetric case either left or right angle can be passed and the other will be calculated using the symmetry.\n", "\n", "\n", "So we are constructing a line space pattern and the cross section of each line looks like following:\n", "\n", "![double trapezoid](../../../assets/images/double_stack.png)\n", "\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create instance of the Simulation class and call `simulate_diffraction` method" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'Stacked Trapezoid diffraction simulation')" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHICAYAAACyBMv/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABk0klEQVR4nO3deXhTZdoG8PskbZLu+97Ssm8iIAgWRUAKKIgyijLIBwVUHAUXGHVEZ1hcQHHDUZRBRdQZBUFxZVBAEB2K7LLv3ei+722a5P3+aE9o6EK35CTp/buuXtrTc5Inp6W5+66SEEKAiIiIyEmolC6AiIiIqD0x3BAREZFTYbghIiIip8JwQ0RERE6F4YaIiIicCsMNERERORWGGyIiInIqDDdERETkVBhuiIiIyKkw3JBT2bVrFyRJwqZNm6z6PDExMZg5c6ZVn8MZzJw5EzExMVc9LykpCZIkYd26da16nnXr1kGSJCQlJZmPjRw5EiNHjrQ4LysrC5MnT0ZAQAAkScLKlSsBAOfOncPYsWPh4+MDSZLw9ddft6qO9tbQa7AXS5YsgSRJijx3W39eGsN/186D4Yba7NixY5g8eTKio6Oh0+kQERGBMWPG4O2337Y4b9myZXbzpmFtM2fOhCRJV/3gL1Lbmj9/Pn788UcsXLgQn376KW699VYAQHx8PI4dO4aXXnoJn376KQYPHmyzmk6ePIklS5ZYBDOynj179mDJkiUoLCxUuhSyIhelCyDHtmfPHowaNQqdOnXCgw8+iNDQUKSmpmLv3r1466238Oijj5rPXbZsGSZPnoxJkyYpV7CNPPTQQ4iLizN/npiYiEWLFmHOnDkYPny4+XjXrl2VKM9m3n//fZhMJkWe+6effqp37Oeff8add96JJ5980nysoqICCQkJeO655zBv3jxblgigJtwsXboUI0eOrNfK1dBrsBd///vf8cwzzyhdRovt2bMHS5cuxcyZM+Hr62vxtTNnzkCl4t/8zoDhhtrkpZdego+PD/bv31/vF0V2drYyRdmB2NhYxMbGmj8/cOAAFi1ahNjYWPzf//1fo9eVlZXBw8PDFiXahKurq2LPrdFo6h3Lzs6u93Oak5MDAPWON8TW35+GXoO9cHFxgYuLc72FaLVapUugdsKISm1y4cIF9O3bt8E3huDgYPP/S5KEsrIyfPzxx/W6ZJKTk/HII4+gZ8+ecHNzQ0BAAO65554Gm+kLCwsxf/58xMTEQKvVIjIyEjNmzEBubm6jNVZVVeH222+Hj48P9uzZAwAwmUxYuXIl+vbtC51Oh5CQEDz00EMoKCiwuFYIgRdffBGRkZFwd3fHqFGjcOLEiZbfqAbI40R++eUXPPLIIwgODkZkZGSL7on8GLt378ZDDz2EgIAAeHt7Y8aMGfVeCwD897//xfDhw+Hh4QEvLy9MmDDB4vXIY5Ya+riyVeHdd99F3759odVqER4ejrlz59Zr6m9ozE1hYSFmzpwJHx8f+Pr6Ij4+vkVdBCdOnMAtt9wCNzc3REZG4sUXX2ywdajueBX5PgkhsGrVKvNrWrJkCaKjowEATz31lMXrlMeUnDx5Evfddx/8/Pxw0003AQCOHj2KmTNnokuXLtDpdAgNDcXs2bORl5dXr460tDTcf//9CA8Ph1arRefOnfHwww9Dr9dj3bp1uOeeewAAo0aNMte1a9eueq9Blp2djfvvvx8hISHQ6XTo378/Pv74Y4tz5DEpr732GtasWYOuXbtCq9Xi+uuvx/79+696j6urq7F06VJ0794dOp0OAQEBuOmmm7Bt2zbzOQ2NuZEkCfPmzcPGjRvRp08fuLm5ITY2FseOHQMA/Otf/0K3bt2g0+kwcuTIej/PjY15ac7Yo+Z8T5YsWYKnnnoKANC5c2fz/ZbraOj5L168iHvuuQf+/v5wd3fHDTfcgB9++MHiHPnfzRdffIGXXnoJkZGR0Ol0GD16NM6fP99k3WQdzhW7yeaio6ORkJCA48eP45prrmn0vE8//RQPPPAAhgwZgjlz5gC43CWzf/9+7NmzB3/+858RGRmJpKQkvPfeexg5ciROnjwJd3d3AEBpaSmGDx+OU6dOYfbs2bjuuuuQm5uLb7/9FpcuXUJgYGC9562oqMCdd96JAwcOYPv27bj++usB1HQbrVu3DrNmzcJjjz2GxMREvPPOOzh8+DD+97//mVscFi1ahBdffBHjx4/H+PHjcejQIYwdOxZ6vb7d7uEjjzyCoKAgLFq0CGVlZS26J7J58+bB19cXS5YswZkzZ/Dee+8hOTnZ/EtX/h7Ex8dj3LhxeOWVV1BeXo733nsPN910Ew4fPoyYmBj07t0bn376qcVjFxYWYsGCBRZhdcmSJVi6dCni4uLw8MMPm59z//79FvfvSkII3Hnnnfjtt9/wl7/8Bb1798bmzZsRHx/frHuVmZmJUaNGwWAw4JlnnoGHhwfWrFkDNze3Jq+7+eab8emnn2L69OkYM2YMZsyYAQC49tpr4evri/nz52Pq1KkYP348PD09La6955570L17dyxbtgxCCADAtm3bcPHiRcyaNQuhoaE4ceIE1qxZgxMnTmDv3r3me56eno4hQ4agsLAQc+bMQa9evZCWloZNmzahvLwcN998Mx577DH885//xLPPPovevXsDgPm/V6qoqMDIkSNx/vx5zJs3D507d8bGjRsxc+ZMFBYW4vHHH7c4/7PPPkNJSQkeeughSJKEFStW4K677sLFixebbFVbsmQJli9fbv43W1xcjAMHDuDQoUMYM2ZMk/f6119/xbfffou5c+cCAJYvX47bb78dTz/9NN5991088sgjKCgowIoVKzB79mz8/PPPTT5eczXne3LXXXfh7Nmz+Pzzz/Hmm2+af2cEBQU1+JhZWVkYNmwYysvL8dhjjyEgIAAff/wx7rjjDmzatAl/+tOfLM5/+eWXoVKp8OSTT6KoqAgrVqzAtGnT8Pvvv7fLa6QWEERt8NNPPwm1Wi3UarWIjY0VTz/9tPjxxx+FXq+vd66Hh4eIj4+vd7y8vLzesYSEBAFAfPLJJ+ZjixYtEgDEV199Ve98k8kkhBBi586dAoDYuHGjKCkpESNGjBCBgYHi8OHD5nN//fVXAUD85z//sXiMrVu3WhzPzs4WGo1GTJgwwfz4Qgjx7LPPCgANvpbG7N+/XwAQH330kfnYRx99JACIm266SRgMhlbdE/kxBg0aZHHPV6xYIQCIb775RgghRElJifD19RUPPvigxWNmZmYKHx+fesdlJpNJ3H777cLT01OcOHFCCHH5vowdO1YYjUbzue+8844AINauXWs+Fh8fL6Kjo82ff/311wKAWLFihfmYwWAQw4cPr3d/GvLEE08IAOL33383H8vOzhY+Pj4CgEhMTDQfHzFihBgxYoTF9QDE3LlzLY4lJiYKAOLVV1+1OL548WIBQEydOrVeHQ19fz7//HMBQOzevdt8bMaMGUKlUon9+/fXO1/+mdq4caMAIHbu3FnvnCtfw8qVKwUA8e9//9t8TK/Xi9jYWOHp6SmKi4stXlNAQIDIz883n/vNN98IAOK7776r91x19e/fX0yYMKHJc+T7UxcAodVqLb4P//rXvwQAERoaaq5PCCEWLlxY73sWHR3d4L+rK++D/Prq/rw093vy6quv1nvexp5f/nn79ddfzcdKSkpE586dRUxMjPnnX/6907t3b1FVVWU+96233hIAxLFjx+o9F1kXu6WoTcaMGYOEhATccccd+OOPP7BixQqMGzcOERER+Pbbb5v1GHX/6q6urkZeXh66desGX19fHDp0yPy1L7/8Ev3796/31xKAes3jRUVFGDt2LE6fPo1du3ZhwIAB5q9t3LgRPj4+GDNmDHJzc80fgwYNgqenJ3bu3AkA2L59O/R6PR599FGLx3/iiSea9bqa68EHH4RarbY41tx7IpszZ47FX+IPP/wwXFxcsGXLFgA1f9UWFhZi6tSpFq9ZrVZj6NCh5td8pRdeeAHff/891q1bhz59+gC4fF+eeOIJi8GXDz74ILy9ves12de1ZcsWuLi44OGHHzYfU6vVFgPPm7JlyxbccMMNGDJkiPlYUFAQpk2b1qzrW+Mvf/lLvWN1vz+VlZXIzc3FDTfcAADm74/JZMLXX3+NiRMnNjj7qjXTqLds2YLQ0FBMnTrVfMzV1RWPPfYYSktL8csvv1icP2XKFPj5+Zk/lwezX7x4scnn8fX1xYkTJ3Du3LkW1zh69GiLrsihQ4cCAO6++254eXnVO361WpqrOd+TltqyZQuGDBli7o4EAE9PT8yZMwdJSUk4efKkxfmzZs2yGCfV3PtN7a9Dh5vdu3dj4sSJCA8Pb9XaFnKf85UfzjQgtDmuv/56fPXVVygoKMC+ffuwcOFClJSUYPLkyfX+8TekoqICixYtQlRUFLRaLQIDAxEUFITCwkIUFRWZz7tw4UKTXV91PfHEE9i/fz+2b9+Ovn37Wnzt3LlzKCoqQnBwMIKCgiw+SktLzQOhk5OTAQDdu3e3uD4oKMjiDaOtOnfuXO9Yc++J7MoaPT09ERYWZh5LIL9J3XLLLfVe808//dTg4O+tW7di6dKlWLhwIe6++27zcfm+9OzZ0+J8jUaDLl26mL/ekOTkZISFhdXr+rnysZq6/srX2pLrW6Oh709+fj4ef/xxhISEwM3NDUFBQebz5O9PTk4OiouLm/0z2xzy679yRo/cjXXlve/UqZPF5/LPbUPjsep6/vnnUVhYiB49eqBfv3546qmncPTo0WbVeOVz+vj4AACioqIaPH61WpqrOd+TlkpOTm7wZ6u97ze1vw495qasrAz9+/fH7Nmzcdddd7X4+ieffLLeX3WjR482j+voaDQaDa6//npcf/316NGjB2bNmoWNGzdi8eLFTV736KOP4qOPPsITTzyB2NhY80Jqf/7zn1s9jfjOO+/E+vXr8fLLL+OTTz6xeDMwmUwIDg7Gf/7znwavbaz/3VoaGi/S3vdEvubTTz9FaGhova9fOeslMTER06ZNw5gxY/Diiy+2+PmcSUPfn3vvvRd79uzBU089hQEDBsDT0xMmkwm33nqrYlPfG3Jli6BM1I4daszNN9+MCxcu4JtvvsFPP/2EDz74AG+++SZWr16NBx54oFXP2ZxaGmvNMhqNjV4vs4fvSWvvN7W/Dh1ubrvtNtx2222Nfr2qqgrPPfccPv/8cxQWFuKaa67BK6+8Yh617+npafEX6B9//IGTJ09i9erV1i7d7snN8BkZGeZjjf3i2rRpE+Lj4/H666+bj1VWVtabQdO1a1ccP368Wc8/adIkjB07FjNnzoSXlxfee+89i8fZvn07brzxxiYHosqzaM6dO4cuXbqYj+fk5Fj9L7Hm3hPZuXPnMGrUKPPnpaWlyMjIwPjx4wFcHrwdHBxssf5OQyoqKnDXXXfB19cXn3/+eb1WAvm+nDlzxuK+6PV6JCYmNvn40dHR2LFjB0pLSy3+7Zw5c6bJmupe31BXSXOvbw8FBQXYsWMHli5dikWLFpmPX1lXUFAQvL29r/oz25LuqejoaBw9ehQmk8ni+3L69Gnz19uLv78/Zs2ahVmzZqG0tBQ333wzlixZctVw0xZ+fn4N/ownJydb/KxdqbnfE6Dl97uhny1r3G9qXx26W+pq5s2bh4SEBKxfvx5Hjx7FPffcg1tvvbXRfugPPvgAPXr0sFikzdnt3Lmzwb9K5LEedZt0PTw8GvzFpVar6z3G22+/DaPRaHHs7rvvxh9//IHNmzfXe4yGapgxYwb++c9/YvXq1fjb3/5mPn7vvffCaDTihRdeqHeNwWAw1xgXFwdXV1e8/fbbFo8vL9lvTc29J7I1a9agurra/Pl7770Hg8FgDu/jxo2Dt7c3li1bZnGeTF7rBagZY3L27Fls3ry5we63uLg4aDQa/POf/7So8cMPP0RRUREmTJjQ6OsaP348DAaDRdg0Go31VrNu6vq9e/di3759FrU31gpnDfJf51d+f678uVCpVJg0aRK+++47HDhwoN7jyNfL3djNmQ4/fvx4ZGZmYsOGDeZjBoMBb7/9Njw9PTFixIiWvJRGXTml3dPTE926dUNVVVW7PH5junbtir1791rMRvz++++Rmpra5HXN/Z4ALb/f+/btQ0JCgvlYWVkZ1qxZg5iYGPM4NLI/HbrlpikpKSn46KOPkJKSgvDwcAA13VBbt27FRx99hGXLllmcX1lZif/85z8OuWJnWzz66KMoLy/Hn/70J/Tq1Qt6vR579uzBhg0bEBMTg1mzZpnPHTRoELZv34433ngD4eHh6Ny5M4YOHYrbb78dn376KXx8fNCnTx8kJCRg+/btCAgIsHiup556Cps2bcI999yD2bNnY9CgQcjPz8e3336L1atXo3///vXqmzdvHoqLi/Hcc8/Bx8cHzz77LEaMGIGHHnoIy5cvx5EjRzB27Fi4urri3Llz2LhxI9566y1MnjwZQUFBePLJJ81TWcePH4/Dhw/jv//9b4PTzttTc++JTK/XY/To0bj33ntx5swZvPvuu7jppptwxx13AAC8vb3x3nvvYfr06bjuuuvw5z//GUFBQUhJScEPP/yAG2+8Ee+88w5++OEHfPLJJ7j77rtx9OhRi3EWnp6emDRpEoKCgrBw4UIsXboUt956K+644w7zc15//fVNLlI4ceJE3HjjjXjmmWeQlJSEPn364Kuvvmr2mIinn37avG3C448/bp4KLrdo2IK3tzduvvlmrFixAtXV1YiIiMBPP/2ExMTEeucuW7YMP/30E0aMGIE5c+agd+/eyMjIwMaNG/Hbb7/B19cXAwYMgFqtxiuvvIKioiJotVrccsstFlPvZXPmzMG//vUvzJw5EwcPHkRMTAw2bdqE//3vf1i5cqXFgN226NOnD0aOHIlBgwbB398fBw4cwKZNm6y+gvMDDzyATZs24dZbb8W9996LCxcu4N///vdVV/Juyfdk0KBBAIDnnnsOf/7zn+Hq6oqJEyc2OFbymWeeweeff47bbrsNjz32GPz9/fHxxx8jMTERX375JVcztmfKTNKyPwDE5s2bzZ9///33AoDw8PCw+HBxcRH33ntvves/++wz4eLiIjIzM21YtfL++9//itmzZ4tevXoJT09PodFoRLdu3cSjjz4qsrKyLM49ffq0uPnmm4Wbm5vFVOqCggIxa9YsERgYKDw9PcW4cePE6dOnG5wWmpeXJ+bNmyciIiKERqMRkZGRIj4+XuTm5gohLKeC1/X0008LAOKdd94xH1uzZo0YNGiQcHNzE15eXqJfv37i6aefFunp6eZzjEajWLp0qQgLCxNubm5i5MiR4vjx441OWW1MU1PBG5om3Nx7Ij/GL7/8IubMmSP8/PyEp6enmDZtmsjLy6v3uDt37hTjxo0TPj4+QqfTia5du4qZM2eKAwcOWDxeQx91p3QLUTP1u1evXsLV1VWEhISIhx9+WBQUFFicc+VUcCFqvofTp08X3t7ewsfHR0yfPl0cPny4WVPBhRDi6NGjYsSIEUKn04mIiAjxwgsviA8//NBqU8FzcnLq1XDp0iXxpz/9Sfj6+gofHx9xzz33iPT0dAFALF682OLc5ORkMWPGDBEUFCS0Wq3o0qWLmDt3rsWU4ffff1906dJFqNVqi2nhDb2GrKws88+GRqMR/fr1q3ffGntN8j24ssYrvfjii2LIkCHC19dXuLm5iV69eomXXnrJYrmBxqaCN/f+NvZv9fXXXxcRERFCq9WKG2+8URw4cKBZU8Fb8j154YUXREREhFCpVBY/Nw39u75w4YKYPHmy8PX1FTqdTgwZMkR8//33zXotDdVJtiEJwZFOQE0/7ObNm837Hm3YsAHTpk3DiRMn6g0S8/T0rDcoc/To0fD29m6wy4TIWuSFCPfv32/TzR6JiOwZu6UaMXDgQBiNRmRnZ191DE1iYiJ27tzZ7HVdiIiIyHo6dLgpLS212PcjMTERR44cgb+/P3r06IFp06ZhxowZeP311zFw4EDk5ORgx44duPbaay0GTa5duxZhYWFNzrwiIiIi2+jQ4ebAgQMW02cXLFgAAIiPj8e6devw0Ucf4cUXX8Rf//pXpKWlITAwEDfccANuv/128zUmkwnr1q3DzJkzr7oOAxEREVkfx9wQERGRU+E8NiIiInIqDDdERETkVDrcmBuTyYT09HR4eXm1aldeIiIisj0hBEpKShAeHn7VBRQ7XLhJT0+vtzstEREROYbU1FRERkY2eU6HCzfy8uSpqanw9vZWuBoiIiJqjuLiYkRFRTVrm5EOF27krihvb2+GGyIiIgfTnCElHFBMREREToXhhoiIiJwKww0RERE5FYYbIiIiciqKhpvdu3dj4sSJCA8PhyRJ+Prrr5s8/6uvvsKYMWMQFBQEb29vxMbG4scff7RNsUREROQQFA03ZWVl6N+/P1atWtWs83fv3o0xY8Zgy5YtOHjwIEaNGoWJEyfi8OHDVq6UiIiIHIXdbJwpSRI2b96MSZMmtei6vn37YsqUKVi0aFGzzi8uLoaPjw+Kioo4FZyIiMhBtOT926HXuTGZTCgpKYG/v3+j51RVVaGqqsr8eXFxsS1KIyIiIoU49IDi1157DaWlpbj33nsbPWf58uXw8fExf3DrBSIiIufmsOHms88+w9KlS/HFF18gODi40fMWLlyIoqIi80dqaqoNqyQiIiJbc8huqfXr1+OBBx7Axo0bERcX1+S5Wq0WWq3WRpURERGR0hyu5ebzzz/HrFmz8Pnnn2PChAlKl0NERER2RtGWm9LSUpw/f978eWJiIo4cOQJ/f3906tQJCxcuRFpaGj755BMANV1R8fHxeOuttzB06FBkZmYCANzc3ODj46PIayAi+yGEQEW1Ee4ah2yUJqJ2omjLzYEDBzBw4EAMHDgQALBgwQIMHDjQPK07IyMDKSkp5vPXrFkDg8GAuXPnIiwszPzx+OOPK1I/EdmPnJIqTPnXXvRf+hM++l8i7GSVCyJSgN2sc2MrXOeGyPkcu1SEOZ8eQEZRpfnYvYMj8cKka6B1UStYGRG1l5a8fzvcmBsiorq+OZKGyav3IKOoEl2CPPDoLd2gkoAvDlzC1DV7kV1SefUHISKnwnBDRA7rmyNpeHz9EVQZTBjVMwhfz70Rfx3bEx/NGgJvnQsOpRRi8nsJqKw2Kl0qEdkQww0ROSQhBN7deQEA8H83dMIH8dfDW+cKABjRIwjfzLsJId5apOSX47s/0pUslYhsjOGGiBzS/qQCnMkqgZurGk+N6wW1SrL4eudAD8wc1hkA8O+9yUqUSEQKYbghIof0aW1gmTQwHD5urg2ec+/gSGjUKvxxqQh/pBbasDoiUhLDDRE5nOySSmw9ngEA+L8bohs9L8BTi9uvDQNwOQwRkfNjuCEih7NhXyqqjQLXdfJF3/CmF/D8v9ia8PPdH+koKNPbojwiUhjDDRE5FIPRhM/21SzuOT228VYb2cAoX/QN90aVwYRNBy9ZuzwisgMMN0TkULafykZGUSX8PTQY3y/squdLkoTptV1X//49GSZTh1q3lKhDYrghIociz3yacn1Us1cfvmNAOLx0LkjOK8fucznWLI+I7ADDDRE5jMTcMvx2PheSBNw3pFOzr3PXuGDyoEgAwL/3plzlbCJydAw3ROQwdp7OBgDc2DUQUf7uLbr23sFRAIDfzuegysAVi4mcGcMNETmMPRfyAAA3dQ9s8bW9Qr0Q6KlBZbUJR1IK27kyIrInDDdE5BCMJoHfE2vCTWyXgBZfL0kSYrvWhKKEi3ntWhsR2ReGGyJyCCfSi1BSaYCXzgV9w71b9RhyKJJbgIjIOTHcEJFDkAPJ0M4BcFG37lfXsK414eZwSgEq9Bx3Q+SsGG6IyCHI4Sa2a8u7pGTRAe4I99Gh2ihwIDm/vUojIjvDcENEdk9vMOFAUk0YGdaGcGMx7oZdU0ROi+GGiOze0UuFKNcb4e+hQc8QrzY9ltzyw3E3RM6L4YaI7J4cRG7o4g+VSmrTY8nh5lhaEUoqq9tcGxHZH4YbIrJ7CebxNi1f3+ZKEb5uiAlwh9EksD+J426InBHDDRHZtcpqIw6mFABo23ibusxdU+fZNUXkjBhuiMiuHUougN5gQrCXFl0CPdrlMeUWII67IXJODDdEZNfk1YSHdQ2AJLVtvI1MXszvVGYxCsr07fKYRGQ/GG6IyK7JrSvD2mG8jSzIS4seIZ4QAuYtHYjIeTDcEJHdqtAb8UdqIYC2Ld7XELn1Zu9FDiomcjYMN0Rkt05nFsNgEgj01CDSz61dH3tgJz8AwPG0onZ9XCJSHsMNEdmtE+nFAIA+4T7tNt5GJm++eTKjGCaTaNfHJiJlMdwQkd2Sw01rdwFvSpcgT+hcVSjXG5GYV9buj09EymG4ISK7dTK9psvIGuFGrZLQK7TmceUQRUTOgeGGiOySwWjC6cwSAEDfcB+rPIccmk6kc9wNkTNhuCEiu3QhpwxVBhM8tS6I9ne3ynPIoekkW26InArDDRHZJXkWU+8wrzZvltmYyy03xRCCg4qJnAXDDRHZpcuDia3TJQUAPUO9oFZJyC/TI7O40mrPQ0S2xXBDRHbphBUHE8t0rmp0C/Kseb40dk0ROQuGGyKyO0IInMywfstNzeNzxhSRs2G4ISK7k5pfgZJKAzRqFbqHeFr1ufpwxhSR02G4ISK7IweNHqGecFVb99eU3DLElhsi58FwQ0R2xzyYOMy6XVLA5ZabtMIKFJbrrf58RGR9DDdEZHfMg4kjrDeYWObj5ooo/5pNObneDZFzYLghIrtjzT2lGiK3EB3nuBsip8BwQ0R2JaekCtklVZAkmPd+sjbOmCJyLgw3RGRX5C6pzoEe8NC62OQ55e4vhhsi58BwQ0R2xRYrE19Jfq6LOaWo0Btt9rxEZB0MN0RkV07aeLwNAAR7aRHoqYFJAKcy2XpD5OgYbojIrpyuDRd9wmwXbiRJQu/a5zudUWKz5yUi61A03OzevRsTJ05EeHg4JEnC119/fdVrdu3aheuuuw5arRbdunXDunXrrF4nEdlGtdGE5LxyALD6ysRX6h7sBQC4kFNq0+clovanaLgpKytD//79sWrVqmadn5iYiAkTJmDUqFE4cuQInnjiCTzwwAP48ccfrVwpEdlCcl45DCYBD40aod46mz5312APAAw3RM7ANlMRGnHbbbfhtttua/b5q1evRufOnfH6668DAHr37o3ffvsNb775JsaNG2etMonIRs5n1wSLrsGekCTJps8t7w4u10BEjsuhxtwkJCQgLi7O4ti4ceOQkJDQ6DVVVVUoLi62+CAi+yS3mnQNsm2XFFATqICabRg4Y4rIsTlUuMnMzERISIjFsZCQEBQXF6OioqLBa5YvXw4fHx/zR1RUlC1KJaJWuFDbatIt2PbhJsBDA193VwgBXMxl6w2RI3OocNMaCxcuRFFRkfkjNTVV6ZKIqBHnzS03HjZ/bkmSzC1G7JoicmyKjrlpqdDQUGRlZVkcy8rKgre3N9zc3Bq8RqvVQqvV2qI8ImoDIYSiLTdAzbibg8kFuJBTpsjzE1H7cKiWm9jYWOzYscPi2LZt2xAbG6tQRUTUXjKLK1GmN8JFJSE6wPYtN8DlUHWBLTdEDk3RcFNaWoojR47gyJEjAGqmeh85cgQpKSkAarqUZsyYYT7/L3/5Cy5evIinn34ap0+fxrvvvosvvvgC8+fPV6J8ImpHF7JrWks6BbjDVa3MryZOBydyDoqGmwMHDmDgwIEYOHAgAGDBggUYOHAgFi1aBADIyMgwBx0A6Ny5M3744Qds27YN/fv3x+uvv44PPviA08CJnMD57JqVgbspMFNK1i2oZiG/i7llMJqEYnUQUdsoOuZm5MiREKLxXyANrT48cuRIHD582IpVEZES5HEuXRUabwMAEX5u0LiooDeYcKmgXLHuMSJqG4cac0NEzkueoaRky41aJaFLoIdFPUTkeBhuiMgumBfwU7Dlpu7zc9wNkeNiuCEixRVXViO7pAqAMmvc1MVtGIgcH8MNESlOnnod4q2Fl85V0Vout9xwrRsiR8VwQ0SKO6/w4n111W25aWrCAxHZL4YbIlKceaaUgoOJZV2CPCBJQFFFNfLK9EqXQ0StwHBDRIqzp5YbnasakX4127lw3A2RY2K4ISLFXTRvmKl8uAEu18EZU0SOieGGiBSlN5iQnF8OwD5abgDOmCJydAw3RKSo5LyarQ48tS4I9tIqXQ4AzpgicnQMN0SkKLl1pGuwJyRJUriaGtwdnMixMdwQkaLMKxMrvHhfXfKYm7TCCpTrDQpXQ0QtxXBDRIq6mGs/08Bl/h4a+LnXLCaYlFuucDVE1FIMN0SkqOS8mvAQHeCucCWWOtXuCJ6cx3E3RI6G4YaIFCWHm5gA++mWAoCY2rAlz+QiIsfBcENEiimtMiC3tGbDzE521nITzZYbIofFcENEipGDQ4CHBt4Kb5h5JbnlhmNuiBwPww0RKUbukrK3Vhvg8hggttwQOR6GGyJSTFJtcLC38TbA5W6pjOJKVFYbFa6GiFqC4YaIFJNipzOlgJquMk+tC4QALhWwa4rIkTDcEJFi7LnlRpIkc+jiuBsix8JwQ0SKsecxN8DlFqUkjrshcigMN0SkiMpqIzKKKgHYZ8sNcHncTQrXuiFyKAw3RKQIOTB46VzMWx3YG/N08DyGGyJHwnBDRIqouzKxvewGfiUu5EfkmBhuiEgRcmCw1/E2wOUxN5cKKlBtNClcDRE1F8MNESni8kwp+w03IV46aF1UMJoE0gsrlC6HiJqJ4YaIFHF5N3D7HEwMACpVnengHHdD5DAYbohIEfa8xk1dHHdD5HgYbojI5vQGE9IKarp57LlbCuAGmkSOiOGGiGwurbACJgG4uaoR5KVVupwmdTKvdcOWGyJHwXBDRDYnd0lFB7jb7TRwGde6IXI8DDdEZHPJuZfDjb2TxwSl5JXDaBIKV0NEzcFwQ0Q2l1RnAT97F+ajg6tagt5oQmZxpdLlEFEzMNwQkc3JWy/Y8wJ+Mhe1CpF+NXXKLU5EZN8YbojI5hxlGrhM7j5L5gaaRA6B4YaIbMpoEkjNlxfws/+WG+ByCEviWjdEDoHhhohsKr2wAtVGAY1ahTAfN6XLaRZzyw3XuiFyCAw3RGRT8rYLkf5uUKvsexq47PIWDGy5IXIEDDdEZFPyYOJof8fokgKATv413VKp+eUQgtPBiewdww0R2VRqQU24iXKgcBPpV9N9VqY3oqC8WuFqiOhqGG6IyKbkwcRRfo4TbnSuagTXbhORyhlTRHaP4YaIbCq1dsPMKH/HGEwsk1ua5JYnIrJfDDdEZFPmlhsH6pYCgE619aaw5YbI7jHcEJHNlFUZkF+mB+B44SaqdtxNan6FwpUQ0dUw3BCRzchdOj5urvDWuSpcTctE1oaxS+yWIrJ7ioebVatWISYmBjqdDkOHDsW+ffuaPH/lypXo2bMn3NzcEBUVhfnz56OykpvZETkCudXD0cbbAJcHQHNAMZH9UzTcbNiwAQsWLMDixYtx6NAh9O/fH+PGjUN2dnaD53/22Wd45plnsHjxYpw6dQoffvghNmzYgGeffdbGlRNRa5g3zHSwLing8iafaYUVMJq41g2RPVM03Lzxxht48MEHMWvWLPTp0werV6+Gu7s71q5d2+D5e/bswY033oj77rsPMTExGDt2LKZOnXrV1h4isg+OOA1cFuqtg6taQrVRILOYrcVE9kyxcKPX63Hw4EHExcVdLkalQlxcHBISEhq8ZtiwYTh48KA5zFy8eBFbtmzB+PHjG32eqqoqFBcXW3wQkTLk8SqRDthyo1ZJCPeVBxWza4rInikWbnJzc2E0GhESEmJxPCQkBJmZmQ1ec9999+H555/HTTfdBFdXV3Tt2hUjR45ssltq+fLl8PHxMX9ERUW16+sgouYzj7nxc7wxNwDH3RA5CsUHFLfErl27sGzZMrz77rs4dOgQvvrqK/zwww944YUXGr1m4cKFKCoqMn+kpqbasGIikgkhzLOlHHHMDVB3IT9OByeyZy5KPXFgYCDUajWysrIsjmdlZSE0NLTBa/7xj39g+vTpeOCBBwAA/fr1Q1lZGebMmYPnnnsOKlX9rKbVaqHVatv/BRBRi+SV6VGuN0KSgAhHbbnxZ7cUkSNQrOVGo9Fg0KBB2LFjh/mYyWTCjh07EBsb2+A15eXl9QKMWq0GAO7US2Tn5EAQ4qWD1kWtcDWtw24pIsegWMsNACxYsADx8fEYPHgwhgwZgpUrV6KsrAyzZs0CAMyYMQMRERFYvnw5AGDixIl44403MHDgQAwdOhTnz5/HP/7xD0ycONEccojIPjnqnlJ1cX8pIsegaLiZMmUKcnJysGjRImRmZmLAgAHYunWreZBxSkqKRUvN3//+d0iShL///e9IS0tDUFAQJk6ciJdeekmpl0BEzeSoe0rVJY8VyiquQmW1ETpX/lFFZI8k0cH6c4qLi+Hj44OioiJ4e3srXQ5Rh/HMl0exfn8qHh/dHfPH9FC6nFYRQuCaxT+iTG/E9gUj0C3YU+mSiDqMlrx/O9RsKSJyXHJXjiO33EiSxK4pIgfAcENENuHoa9zIImsHFV/ioGIiu8VwQ0RWZzQJpBfKA4odt+UGqDMdnGvdENkthhsisrqMogoYTAIatQoh3jqly2kTeVAxp4MT2S+GGyKyOnk38Ag/N6hVksLVtI281k0Kww2R3WK4ISKru1Q73ibSwcfbAHXWumG4IbJbDDdEZHXOMFNKJge04koDiiqqFa6GiBrCcENEVie3cjjqhpl1eWhdEOipAcDWGyJ7xXBDRFYnj0+Rx6s4ukjuMUVk1xhuiMjqnGFfqbq4kB+RfWO4ISKrqqw2IqekCoDztNzICxHKCxMSkX1huCEiq7pU22rjqXWBr7urwtW0D7nl5hJbbojsEsMNEVmVHAAi/dwgSY69xo1MboG6xFWKiewSww0RWZUcACKdpEsKuDwd/FJBBYQQCldDRFdiuCEiq0qt03LjLMJ8dZAkoKLaiLwyvdLlENEVGG6IyKout9w4T7jRuqgR4lWzRxa7pojsD8MNEVmVM3ZLAXW7pjiomMjeMNwQkVWlOWG3FGA57oaI7AvDDRFZTYXeiNzSmjEpzrCvVF2cDk5kvxhuiMhq5Dd+L50LfNycY40bWSQX8iOyWww3RGQ1zjreBrj8mthyQ2R/GG6IyGouOel4G4Br3RDZM4YbIrIaZ5wGLgvzcYMkAVUGk3lcERHZB4YbIrIaOdw4y4aZdWlcVAjzlte6YdcUkT1huCEiq3Hmbimg7rgbDiomsicMN0RkNalOPKAYqDNjii03RHalVeEmPj4eu3fvbu9aiMiJlFUZkF+771KE07bccCE/InvUqnBTVFSEuLg4dO/eHcuWLUNaWlp710VEDi6tsOYN39sJ17iRsVuKyD61Ktx8/fXXSEtLw8MPP4wNGzYgJiYGt912GzZt2oTq6ur2rpGIHNDl8TbO2SUFcH8pInvV6jE3QUFBWLBgAf744w/8/vvv6NatG6ZPn47w8HDMnz8f586da886icjBmGdK+TtnlxRweQuGNK51Q2RX2jygOCMjA9u2bcO2bdugVqsxfvx4HDt2DH369MGbb77ZHjUSkQNy5tWJZaE+Oqhq17rJKa1SuhwiqtWqcFNdXY0vv/wSt99+O6Kjo7Fx40Y88cQTSE9Px8cff4zt27fjiy++wPPPP9/e9RKRg0jNd+5p4ADgqlYhzId7TBHZG5fWXBQWFgaTyYSpU6di3759GDBgQL1zRo0aBV9f3zaWR0SOqiO03AA1M8HSCitwqaAcg6L9lC6HiNDKcPPmm2/innvugU6na/QcX19fJCYmtrowInJszr6AnyzSzw37EjljisietKpbaufOnQ3OiiorK8Ps2bPbXBQRObbSKgMKymt+Rzh7uInidHAiu9OqcPPxxx+joqL+P+SKigp88sknbS6KiBxbWu0bva+7K7x0zrnGjYzTwYnsT4u6pYqLiyGEgBACJSUlFt1SRqMRW7ZsQXBwcLsXSUSOpaN0SQGXxxSlseWGyG60KNz4+vpCkiRIkoQePXrU+7okSVi6dGm7FUdEjsk8U8rXuQcTA5ZbMJhMAiqVpHBFRNSicLNz504IIXDLLbfgyy+/hL+/v/lrGo0G0dHRCA8Pb/ciicixXJ4p5fwtN2E+OqhVEvTGmrVuQrwbn2hBRLbRonAzYsQIAEBiYiI6deoESeJfKERUX0cKNy5qFUK9debp4Aw3RMprdrg5evQorrnmGqhUKhQVFeHYsWONnnvttde2S3FE5JguFTr/vlJ1RfnXrHWTml+BQdFKV0NEzQ43AwYMQGZmJoKDgzFgwABIktTgXiqSJMFoNLZrkUTkWOTVeuW9l5xdTYjL54wpIjvR7HCTmJiIoKAg8/8TETWkuLIaRRU1a9xEdIBuKcByUDERKa/Z4SY6OrrB/yciqkueEu3n7gpPbasWQXc4kVzIj8iutHoRvx9++MH8+dNPPw1fX18MGzYMycnJ7VYcETkeeRp4R+mSAoCo2pabVHZLEdmFVoWbZcuWwc2t5h9zQkIC3nnnHaxYsQKBgYGYP39+uxZIRI6lI82UkkXWBrn0wgoYTfXHIhKRbbUq3KSmpqJbt24AgK+//hqTJ0/GnDlzsHz5cvz6668teqxVq1YhJiYGOp0OQ4cOxb59+5o8v7CwEHPnzkVYWBi0Wi169OiBLVu2tOZlEJEVdJTdwOsK9dbBRSWh2iiQXVKpdDlEHV6rwo2npyfy8vIAAD/99BPGjBkDANDpdA3uOdWYDRs2YMGCBVi8eDEOHTqE/v37Y9y4ccjOzm7wfL1ejzFjxiApKQmbNm3CmTNn8P777yMiIqI1L4OIrEDumonqQC03apWEcN/arql8jrshUlqrRvuNGTMGDzzwAAYOHIizZ89i/PjxAIATJ04gJiam2Y/zxhtv4MEHH8SsWbMAAKtXr8YPP/yAtWvX4plnnql3/tq1a5Gfn489e/bA1bVmM76WPB8RWV9HbLkBarrhUvLLcamgHEM6+1/9AiKymla13KxatQqxsbHIycnBl19+iYCAAADAwYMHMXXq1GY9hl6vx8GDBxEXF3e5GJUKcXFxSEhIaPCab7/9FrGxsZg7dy5CQkJwzTXXYNmyZU2uq1NVVYXi4mKLDyKyno60aWZdnA5OZD9a1XLj6+uLd955p97xlmyamZubC6PRiJCQEIvjISEhOH36dIPXXLx4ET///DOmTZuGLVu24Pz583jkkUdQXV2NxYsXN3jN8uXLuZknkY0UlVejpNIAoOO13ETVvl55thgRKafVi1AUFhZi3759yM7OhslkMh+XJAnTp09vl+KuZDKZEBwcjDVr1kCtVmPQoEFIS0vDq6++2mi4WbhwIRYsWGD+vLi4GFFRUVapj6ijk8fbBHpq4KZRK1yNbUX6s+WGyF60Ktx89913mDZtGkpLS+Ht7W2xgWZzw01gYCDUajWysrIsjmdlZSE0NLTBa8LCwuDq6gq1+vIvzd69eyMzMxN6vR4ajabeNVqtFlqttrkvjYjaQH5jj+hgrTZAnYX8CtlyQ6S0Vo25+etf/4rZs2ejtLQUhYWFKCgoMH/k5+c36zE0Gg0GDRqEHTt2mI+ZTCbs2LEDsbGxDV5z44034vz58xYtRWfPnkVYWFiDwYaIbOtSB5wpJZO7pdILK2Ewmq5yNhFZU6vCTVpaGh577DG4u7ftr7MFCxbg/fffx8cff4xTp07h4YcfRllZmXn21IwZM7Bw4ULz+Q8//DDy8/Px+OOP4+zZs/jhhx+wbNkyzJ07t011EFH76KgzpQAg2EsLV7UEo0kgs5hr3RApqVXdUuPGjcOBAwfQpUuXNj35lClTkJOTg0WLFiEzMxMDBgzA1q1bzYOMU1JSoFJdzl9RUVH48ccfMX/+fFx77bWIiIjA448/jr/97W9tqoOI2kdHnSkFACqVhAhfNyTlleNSQUWHDHhE9qJV4WbChAl46qmncPLkSfTr18+85ozsjjvuaPZjzZs3D/PmzWvwa7t27ap3LDY2Fnv37m1RvURkG/ICdh1pX6m6ovzdkZRXjtT8ctzQJUDpcog6rFaFmwcffBAA8Pzzz9f7miRJTa47Q0TOSQjRoVtuAK51Q2QvWhVu6g7oJSICgMLyapTpa/6wifDtqOGmdsYUww2Rolo1oLiuykoOnCOiy2vcBHtpoXPtWGvcyOSWG/leEJEyWhVujEYjXnjhBURERMDT0xMXL14EAPzjH//Ahx9+2K4FEpFjuDxTqmO22gCXW27S2HJDpKhWhZuXXnoJ69atw4oVKyzWl7nmmmvwwQcftFtxROQ4zGvcdNDBxAAQVbtKcUZRBaq51g2RYloVbj755BOsWbMG06ZNs1gtuH///o3uC0VEzk2eKdWRW26CPLXQuqhgEkBGIbvsiZTS6kX8unXrVu+4yWRCdXV1m4siIsdzeaZUx225kSQJEeYZUxx3Q6SUVoWbPn364Ndff613fNOmTRg4cGCbiyIixyOPuYnqwOEGuPz6OWOKSDmtmgq+aNEixMfHIy0tDSaTCV999RXOnDmDTz75BN9//31710hEdq5mjRt2SwGcMUVkD1rVcnPnnXfiu+++w/bt2+Hh4YFFixbh1KlT+O677zBmzJj2rpGI7FxemR4V1UZIEhDmq1O6HEVxrRsi5bWq5QYAhg8fjm3btrVnLUTkoOQ38lBvHbQuHXONG5k8Y4pjboiU06qWmy5duiAvL6/e8cLCwjZvpklEjic1v2Nvu1CX3HIjzx4jIttrVbhJSkpqcP+oqqoqpKWltbkoInIsl8fbdOzBxMDlgJdVUokqA/fZI1JCi7qlvv32W/P///jjj/Dx8TF/bjQasWPHDsTExLRbcUTkGFI7+IaZdQV4aODmqkZFtRHphZXoHOihdElEHU6Lws2kSZMA1KzlEB8fb/E1V1dXxMTE4PXXX2+34ojIMcjdUh15dWKZJEmI8nfD2axSpOaXM9wQKaBF4UbeDbxz587Yv38/AgMDrVIUETkWOdx0YrgBUHMfzmaVIiWfg4qJlNCq2VKJiYntXQcROSijSSCtsHYBP4YbAHUHFTPcECmh1VPBd+zYgR07diA7O9vcoiNbu3ZtmwsjIseQWVyJaqOAq1pCqHfHXuNGJrdgcSE/ImW0KtwsXboUzz//PAYPHoywsDBIktTedRGRg0jJu7ynlFrF3wXA5XDDbikiZbQq3KxevRrr1q3D9OnT27seInIwXOOmPrl7Tg5+RGRbrVrnRq/XY9iwYe1dCxE5ILnrhYOJL5NXKS6uNKCovFrhaog6nlaFmwceeACfffZZe9dCRA4ohdPA63HXuCDQUwOA426IlNCqbqnKykqsWbMG27dvx7XXXgtXV1eLr7/xxhvtUhwR2b8UTgNvUJS/O3JL9UjJL8c1ET5Xv4CI2k2rws3Ro0cxYMAAAMDx48fbsx4icjDyHkoMN5Y6+bvjcEohp4MTKaBV4Wbnzp3tXQcROaByvQG5pVUAgCjuK2VBvh+cMUVkey0KN3fddddVz5EkCV9++WWrCyIixyG32njrXODj7nqVszsWTgcnUk6Lwk3djTKJiMzbLgSw1eZK8gBrecd0IrKdFoWbjz76yFp1EJEDMs+UYpdUPfJ08EsF5TCaBBc4JLKhVk0FJyICOFOqKWE+bnBRSag2CmQWVypdDlGHwnBDRK12qYBr3DRGrZLMqzZzxhSRbTHcEFGrcQG/pkVxUDGRIhhuiKhVhBBc4+Yq5HDDlhsi22K4IaJWyS3Vo6LaCEkCIny5aWZDOjHcECmC4YaIWkXuagnz1kHjwl8lDeFCfkTK4G8kImqVVI63uarLC/lxrRsiW2K4IaJWYbi5Onmtm9zSKlTojQpXQ9RxMNwQUatwjZur83FzhZeuZq3U1AJ2TRHZCsMNEbUKw83VSZJ0uWsqj+GGyFYYboioVeQ9k+SuF2qYPKiYLTdEtsNwQ0QtpjeYkF4khxu23DRF3lSUM6aIbIfhhohaLL2wAkIAOlcVgjy1Spdj1y4v5McZU0S2wnBDRC1WdzdwSeJu102Jqt1fKiW/TOFKiDoOhhsiarFkDiZutugADwA1gVAIoXA1RB0Dww0RtVhybk0rREygh8KV2L8IXzeoVRIqq03ILqlSuhyiDoHhhohaLCmvNtwEsOXmajQuKvPeW4m57JoisgWGGyJqsaTaNVvkLhdqmtzClZzHcENkC3YRblatWoWYmBjodDoMHToU+/bta9Z169evhyRJmDRpknULJCIzo0mYF6TrzG6pZpFbuJK4kB+RTSgebjZs2IAFCxZg8eLFOHToEPr3749x48YhOzu7yeuSkpLw5JNPYvjw4TaqlIgAILO4EnqjCa5qCWE+OqXLcQhyCxdbbohsQ/Fw88Ybb+DBBx/ErFmz0KdPH6xevRru7u5Yu3Zto9cYjUZMmzYNS5cuRZcuXWxYLREl1Y4bifJzh4ta8V8hDkFuuUnMZcsNkS0o+ptJr9fj4MGDiIuLMx9TqVSIi4tDQkJCo9c9//zzCA4Oxv3333/V56iqqkJxcbHFBxG1njyYOJqDiZut7pgbTgcnsj5Fw01ubi6MRiNCQkIsjoeEhCAzM7PBa3777Td8+OGHeP/995v1HMuXL4ePj4/5Iyoqqs11E3VkybXjRjgNvPki/dygkoByvRE5pZwOTmRtDtWmXFJSgunTp+P9999HYGBgs65ZuHAhioqKzB+pqalWrpLIucndUjGcKdVsWhc1wmungydzUDGR1bko+eSBgYFQq9XIysqyOJ6VlYXQ0NB651+4cAFJSUmYOHGi+ZjJZAIAuLi44MyZM+jatavFNVqtFlot974hai/J5mng7JZqiZgAD1wqqEBSbhmuj/FXuhwip6Zoy41Go8GgQYOwY8cO8zGTyYQdO3YgNja23vm9evXCsWPHcOTIEfPHHXfcgVGjRuHIkSPsciKyMpNJ1FnAjy03LRETKE8H54wpImtTtOUGABYsWID4+HgMHjwYQ4YMwcqVK1FWVoZZs2YBAGbMmIGIiAgsX74cOp0O11xzjcX1vr6+AFDvOBG1v6ySSlQZTHBRSYis3RCSmkcOg1zrhsj6FA83U6ZMQU5ODhYtWoTMzEwMGDAAW7duNQ8yTklJgUrlUEODiJxWUu1U5kg/N04DbyGudUNkO4qHGwCYN28e5s2b1+DXdu3a1eS169ata/+CiKhByeZp4OySail5rZvk3JrdwSVJUrgiIufFP72IqNkSuWFmq0X5u0OSgJIqA/LK9EqXQ+TUGG6IqNmSc7lhZmvpXNUI95Gng7NrisiaGG6IqNnkmT7cMLN15OnzSdyGgciqGG6IqFmEEFzjpo3qbsNARNbDcENEzZJdUoWKaiNUEhDpx3DTGuYNNDkdnMiqGG6IqFnkbRci/NygceGvjtbgdHAi2+BvKCJqFvOGmRxM3GryvUvM5e7gRNbEcENEzcJtF9pOHqtUUmlAYXm1wtUQOS+GGyJqFg4mbjudqxphPjoA3GOKyJoYboioWRJz2XLTHszTwRluiKyG4YaIrqpmGnhtuOEaN21i3kCTa90QWQ3DDRFdVU5pFcr0RkgSEOXP3cDbQg6HbLkhsh6GGyK6qgvZNW/EUX7u0LqoFa7GsXWpDTcXckoVroTIeTHcENFVyW/EXYPYJdVWXYM9AdQERpOJ08GJrIHhhoiu6nx2TbjpVvvGTK3Xyd8drmoJFdVGZBRXKl0OkVNiuCGiq7rccsNw01auapV5peIL2eyaIrIGhhsiuqoLbLlpV91qQ+J5hhsiq2C4IaImlVUZkF5U033Clpv20TWYg4qJrInhhoiadDGnZqZUgIcGfh4ahatxDl3ZckNkVQw3RNQkjrdpf3L33oUcrnVDZA0MN0TUJHO44XibdtOlNijmllahiBtoErU7hhsiapLcdcI1btqPp9bFvIHmeY67IWp3DDdE1CS55YYzpdqX3M3HQcVE7Y/hhogaZTCazLuBc8xN+zKPu+GgYqJ2x3BDRI1KLahAtVFA56pChC83zGxPcjcfW26I2h/DDRE1Sh5v0yXQEyqVpHA1zoXTwYmsh+GGiBrFmVLWI3dLpeSXo8pgVLgaIufCcENEjTJvmMnxNu0uyEsLL60LTAJIyi1Xuhwip8JwQ0SNutxyw2ng7U2SJHOLGMfdELUvhhsiapAQ4nLLDbulrILjboisg+GGiBqUU1qFkkoDVBIQE8CWG2voxpYbIqtguCGiBl3IrlnfJsrfHTpXtcLVOCdOByeyDoYbImrQeW6YaXXmMTfZZTCZhMLVEDkPhhsiatAFjrexuk7+7nBVS6ioNiKjuFLpcoicBsMNETXIPFOKG2Zajataheja8UwcVEzUfhhuiKhB57LYcmML8hpC57JKFK6EyHkw3BBRPQVlemTWdpP0DPVWuBrn1ivMCwBwKoPhhqi9MNwQUT2nMooB1IwJ8dS6KFyNc+sdVhMe5XtORG3HcENE9ZysfaPtXduqQNbTpzbcnM8uRbXRpHA1RM6B4YaI6pG7SORWBbKeSD83eGldoDeauN4NUTthuCGiek6ZW24YbqxNkqQ6427YNUXUHhhuiMhCtdFknpbch+HGJi6Pu+GgYqL2wHBDRBYu5JRCbzTBU+uCSD83pcvpEDiomKh9MdwQkQX5DbZXqBckSVK4mo6hVyi7pYjaE8MNEVk4zcHENtcz1AuSBOSW6pFdwm0YiNqK4YaILJzkYGKbc9e4oHPtNgwcd0PUdnYRblatWoWYmBjodDoMHToU+/bta/Tc999/H8OHD4efnx/8/PwQFxfX5PlE1DKXp4FzjRtb4rgbovajeLjZsGEDFixYgMWLF+PQoUPo378/xo0bh+zs7AbP37VrF6ZOnYqdO3ciISEBUVFRGDt2LNLS0mxcOZHzySmpQm5pFSSppquEbKc3p4MTtRvFw80bb7yBBx98ELNmzUKfPn2wevVquLu7Y+3atQ2e/5///AePPPIIBgwYgF69euGDDz6AyWTCjh07bFw5kfOR31g7B3jAXcNtF2yJLTdE7UfRcKPX63Hw4EHExcWZj6lUKsTFxSEhIaFZj1FeXo7q6mr4+/s3+PWqqioUFxdbfBBRw7h4n3Lke34hpwyV1UaFqyFybIqGm9zcXBiNRoSEhFgcDwkJQWZmZrMe429/+xvCw8MtAlJdy5cvh4+Pj/kjKiqqzXUTOatT3FNKMWE+Ovi4ucJoEuZFFImodRTvlmqLl19+GevXr8fmzZuh0+kaPGfhwoUoKioyf6Smptq4SiLHwT2llCNJkjlUnmTXFFGbKNqpHhgYCLVajaysLIvjWVlZCA0NbfLa1157DS+//DK2b9+Oa6+9ttHztFottFptu9RL5MyqDEbzxo0MN8roHeaNvRfzOe6GqI0UbbnRaDQYNGiQxWBgeXBwbGxso9etWLECL7zwArZu3YrBgwfbolQip3cuqxQGk4CPmyvCfBpuCSXr4qBiovah+HSIBQsWID4+HoMHD8aQIUOwcuVKlJWVYdasWQCAGTNmICIiAsuXLwcAvPLKK1i0aBE+++wzxMTEmMfmeHp6wtPTU7HXQeTo6o634bYLypA3Kj2dWQIhBL8PRK2keLiZMmUKcnJysGjRImRmZmLAgAHYunWreZBxSkoKVKrLDUzvvfce9Ho9Jk+ebPE4ixcvxpIlS2xZOpFT4Xgb5XUL9oRaJaGwvBoZRZUI9+XGpUStoXi4AYB58+Zh3rx5DX5t165dFp8nJSVZvyCiDuiPS4UAgH4RPsoW0oHpXNXoEeKFUxnF+CO1kOGGqJUcerYUEbUPvcGEY2lFAIDrOvkpXE3Hdl0nXwDAoZQCZQshcmAMN0SEE+lF0BtM8PfQIDrAXelyOjQ5XB5KKVS2ECIHxnBDROY30oFRvhzEqrDromvCzbG0msBJRC3HcENE5i4Q+Y2VlBMT4A4/d1foDSacSC9Suhwih8RwQ0Q4nFwTbgbWjvcg5UiShIHsmiJqE4Ybog4us6gS6UWVUElA/0hfpcshcFAxUVsx3BB1cPIbaK9Qb3ho7WJ1iA5PHlQst6gRUcsw3BB1cIeS5fE2vsoWQmb9o3yhkoD0okpkFlUqXQ6Rw2G4IergzIOJub6N3fDQuqBnaM1K0YfZNUXUYgw3RB1YlcGI42k1e0oNZLixKxx3Q9R6DDdEHdiJ9GLojTWL98Vw8T67wsX8iFqP4YaoA5PH23DxPvvDxfyIWo/hhqgDO1zbKsDF++xPTIA7/D00XMyPqBUYbog6MHk8Bxfvsz+SJGFglC8Adk0RtRTDDVEHlVFUgQwu3mfX5BY1DiomahmGG6IO6mAyF++zd3KL2qHkAgghlC2GyIEw3BB1UL+ezQUADO3ir3Al1JiBUX7QuKiQUVSJCzmlSpdD5DAYbog6ICEEfjmbAwAY2TNY4WqoMW4aNYZ2rgmfu87kKFwNkeNguCHqgE5nliCzuBI6V5X5zZPskxw+GW6Imo/hhqgDkt8oY7sEQOeqVrgaasrInkEAgH2J+SirMihcDZFjYLgh6oB+OZsNABjRI0jhSuhqugR6INLPDXqjCXsv5ildDpFDYLgh6mBKKqtxIKlmphTH29g/SZLMrTfsmiJqHoYbog7mf+fzYDAJxAS4IybQQ+lyqBlG9qgdd3M2m1PCiZqB4Yaog5G7pNhq4ziGdQuARq1Can4FLuaWKV0Okd1juCHqQIQQ+KW2a2NET463cRTuGhcMqZ3V9gu7poiuiuGGqAM5l12K9KJKaFxUuKFzgNLlUAvIg793nWW4IboahhuiDmTXmZouqRu6BMBNwyngjkQeVLz3Yh4q9EaFqyGybww3RB2IPNtmJKeAO5xuwZ6I8HWD3sAp4URXw3BD1EGUVhmwPykfwOVWAHIckiSZx0ntrG2BI6KGMdwQdRDf/5GOaqNAlyAPdOYUcIc0pncIAOD7oxmoMrBriqgxDDdEHcSGA6kAgCmDoyBJksLVUGsM7x6IEG8t8sv02H6SrTdEjWG4IeoAzmSW4HBKIVxUEu66LlLpcqiVXNQq3DMoCsDlsEpE9THcEHUAG/bXvBHG9Q5BkJdW4WqoLe4dXBNufj2Xg0sF5QpXQ2SfXJQugMhRCCFQUmVAdnElsoqrkF1SCTdXNSJ83RHh5wY/d1e77O6pMhjx1eFLAIApQ6IUrobaqlOAO27sFoD/nc/DxgOXMH9MD6VLapDeYEJmUSUuFZYjs6gS7hoXhHhrEeKtQ5CXFq5q/m1N1sNwQ9SEymojEi7mYcepLPx8KhvpRZWNnuuhUWNYt0Dc0T8ccb1D7GYdmZ9OZKGwvBphPjrc3J2zpJzBlOs71YabVDw2ujvUKvsI1RdzSvHdHxn44Vg6zmWXorFtsNQqCYOj/TC6dzBG9w5B1yBP2xZKTo/hhqgBZzJL8P6vF/HD0QxUVFvOSvHWuSDEW4dgby3KqoxIK6xATkkVyvRGbDuZhW0ns+CuUWNsnxA8MqobeoR4KfQqashdUvcMirSbN0Fqm7F9QuDr7or0okr8ei5H0X3CDEYTNh28hH//nozjacUWX9O6qBDh64YwXx3K9UZk17Z4VhsFfk/Mx++J+Vi25TS6BnkgflgMJg+KhLuGb0vUdpLoYFvMFhcXw8fHB0VFRfD29la6HLIjQggkXMjDv3ZfxC91lrgP9dbhlt7BiOsdjKGdA+Chrf/Lt7LaiPPZpdhyLAPfHU1Han4FAEAl1YyRmD+mB0K8dTZ7LbLU/HIMX7ETkgTsfmoUovzdbV4DWcfS707go/8l4bZrQvHe/w2y+fMLIbDjVDZe3noa57NLAdS0yNzULRAT+4fj5u6BCPLS1uuqNZkEUgvKsfN0Nnaczsbei3moNta8Dfm6u2LGDdGYHhvDsWFUT0vevxluiAAcu1SEl7acxN6LNYvcSRJwa99Q3H9TZwyK9mvRWBohBA6nFuJfv1zAjyeyAABurmrMubkL5o7qBo2L7cYavPHTGfzz5/MY3j0Qn94/1GbPS9Z3OrMYt678FS4qCXufHY1AT9uFgdOZxVj8zQn8nljz78XX3RVzR3bDXddFIKCFdZRUVmPz4TR88GsiUvJrBki7uarx0IgueHB4lwb/mKCOieGmCQw3VFdaYQVe+/EMNh9OAwBoXFSYMjgKDwzvjOiAti90dyApH8u2nMKhlEIAQO8wb6ycMgA9Q63fVVVcWY1bXtuF3FI93rlvIG6/Ntzqz0m2NWnV/3AktRCPjOyKp2/tZfXnM5kEPvwtEa/+eAZ6owlaFxVm39QZfxnRFT5urm16bKNJ4KcTmVj9ywX8cakIABDspcVfx/bA5EFR7FIlhpumMNwQUNONtPqXC3hv1wVUGUwAgD8NjMBfx/ZApF/7dt0IIfD90Qws/vYE8sv00Lio8PS4nph9Y2eorPgL+4XvT+LD3xLRJcgDWx+/2aYtRmQbP57IxEOfHoRGrcJP829GjBVXnr5UUI4nN/5hbt0c3SsYL0y6BuG+bu36PEIIbDmWiZe3njJ37/YO88YLd/bF4Bj/dn0uciwMN01guKHtJ7Ow9PsT5l+cQzr74+8TeuPaSF+rPm92SSWe+fIYfj5ds7LsTd0C8eaUAVYZW3AmswTj//krjCaBT2YPwc3cKNMpCSEQ/9F+7D6bg1t6BWPtzOut8jxbjmXgb5uOoqTKAHeNGotu74Mp11t3pesqgxGfJiTjnzvOobjSAAC467oILLytN8fjdFAMN01guOm4knLL8Pz3J83hItRbh+cm9Mbt14bZbH0aIQQ+35eKF74/iYpqI4K9tHjnvuswpHP7/UUqhMDU9/di78V83No3FKun236wKdnOhZxS3LpyN6qNAh/GD8bo2v2n2oPeYMKyLaewbk8SAGBgJ1+8ee8Aq7YQXSmvtAortp4xr8jspXXBE2N6YEZsNNfK6WAYbprAcNPxlFUZ8M7O8/jw10TojSa4qCTcP7wzHrulu2KDFc9nl+Dhfx/CuexSqFUSnhrXE3OGd2mXbqrv/kjHo58fhtZFhR1/HdHu3Wxkf17+72ms/uUCOvm746f5N0Pn2vY1ltIKKzD3P4dwJLUQAPCXEV3x5NgecFEoUBxOKcCib07gWFrNeJxuwZ5YMrEvbuoeqEg9ZHsMN01guLFkMJqQX6ZHXpke+bUfBpMJJhNgFAISAG83V/i6ucLXXYMATw0CPDR2uRLvlYwmga8OXcJrP51BVnEVAODmHkFYdHsfdAtWftGwcr0Bf998HF/VDmYe2TMIr07u36Ym97IqA0a//gsyiyuxYEwPPDa6e3uVS3asrMqAW17fhaziqnb5vm89noG/fXkMRRXV8HFzxev39Edcn/ZrEWoto0lg44FUrPjxDPLL9ACAcX1D8PStvRxmIUB55ebCCj0Ky6tRWFGNSr0RkgSoJAlqlQR3jRoBnhr4e2gR4KmBt65tg7WdBcNNEzpquKmsNuJMZglOZhTjbFYJknLLkJRXjtT8chhMLfsR8NK6IDrQHTEBHuge7IVrI31wTYSP3fSDCyGw9XgmXt921rz+Rid/dyy6vQ9G9w62q2AmhMD6/alY8u0JVBlMCPTU4NXJ/TGqV8sXZTMYTXhq01FsPpzWrn/Bk2P49o90PFbbYvfJ7CEY2iWgxY9Rrjfg+e9OYn3two/9I33wzn3X2d36SEXl1Vi54yw+SUiG0SSgkoDJgyLxeFwPRLTzAOfWMhhNOJddimNpRTiRVoSLuWVIzivHpYJytPBXLnzcXBET4I6YQA90CfRErzAv9A33RoSvm139PrM2hpsmdIRwYzCacDqzBEdSC3EktRBHLxXiQk4ZjI38i1JJgJ+7Bv4eGvh5aKB1UUElSVBJgFHUrENRVF6Noopq5JfrG11SPcxHh+ui/TAkxh9DOvujZ4iXVWcDXanaaMJPJ7Kw+pcL5qZrX3dXPDyiK+KHxdj1G/2ZzBI89vlhnMkqAQDMHBaDJ8f1hGczu80qq42Y99lhbD+VBZUErJ15vaKr1pLtCSEwa91+7DqTA42LCm9PHYhxfUObff2R1EIs2HAEF3PLIEk13VDz43rY9Sy7s1klePXHM9h2smY9KY1ahSnXR2HWjTHoYuOWnOLKahxMKsC+pHzsT8zH8fQiVFabGjxX66KCv4cGPm6u8HV3hZurGgKASdRMty+pMiC/rAr5pXqU6Y0NPgZQE3r6hnujf5QvBkT5YmCUL4IVWCzUVhwu3KxatQqvvvoqMjMz0b9/f7z99tsYMmRIo+dv3LgR//jHP5CUlITu3bvjlVdewfjx45v1XM4WboQQyCiqNAeZwykFOJbW8D8qfw8N+oZ7o2eIFzoHeaBzgAdiAj0Q4q1r9hoSVQYjUvPLkZhbjqTcMpzKKMbRtCJcyKm/j4y3zgVDOgfghi7+uKFLAHqHeVtlrYqckiqs35eC//yegszimr2fPDRqPDC8C+4f3tlhmnQrq414+b+nzYM3Az01ePSW7pg6pFOTbzBF5dV44JP92J9UAG3tm9rYFrypkfOorDbi0c8PY9vJmpC77E/98OchnZq8JjG3DK/9dAY/HM0AUDPQ/o0p/TGsq+OMZTmUUoDXfjyDPRfyzMdu7hGE+NhojOwZbJXfO0UV1difmI+9F/OwNzEPJ9OL67XIeGpdcE2EN/pF+KB7sBdiAj0QE+De4MrNjSnXG5CSX/P7NjG3HOezS3EyoxjnskoabHWP8HWrCTqdagJP33Afu9nnrq0cKtxs2LABM2bMwOrVqzF06FCsXLkSGzduxJkzZxAcXP8vzz179uDmm2/G8uXLcfvtt+Ozzz7DK6+8gkOHDuGaa6656vNZK9wYjCaUVhng42a9naGrDEYk1f5wn8ooxrG0IhxPK0Jebd9zXV46FwyoTfPXRvqiX4QPQryb/w+qpUqrDDh2qQgHkvKxLykfB5MLUH7FXxxeOhdc18kPg6L9cF0nP1wb5dPq4HExpxQ7TmVj+6ksHEguMLdKBXpqcd+QKMQPi2nxSqn2YueZbCz59gSS82pWa+3k745HRnbFDV0CEB3gbv4eZtbuK/T+rxdxNqsUXjoXfBh/fbvOvCLHYzCa8Nzm4+bZRf93QyeM7xeGQdF+0LrUvMmV6w04kV6Mb46kYf2+VBhMApJUs9bTPyb0gZ+HRsmX0Gp7LuRi7W+J2HE62/zHVqCnBqN61mzQObx7YKsmEZhMAsn55TicUoCDyTUfZ7JK6v1BFxPgjutrW64HdvJDl0APq7VeVxmMOJdV0+11JKXmj9uz2fVrUklA92Av9K0NWT1CvNA92LNFAaulqo0mlFcZ4ePevn9YOlS4GTp0KK6//nq88847AACTyYSoqCg8+uijeOaZZ+qdP2XKFJSVleH77783H7vhhhswYMAArF69+qrPZ61wcyazBONW7oabqxphPjqE+eoQ6u2GQC8NAmsHhfm5a+ChdYG7Rg03jRoatQpCAAICQgDleiOKK6tRXFHTBZRZVIn0okqkF1YgNb8cyfnlDXYtqVUSeoZ4YWAnXwzs5IcBUb5W/UfVHAajCcfTi/H7xTzsvZiH/UkFKK0y1Dsv1FuH7iGe6BHihTAfnfk+ebu5wmAUqDIYUVltQn5ZFc5lleJsdinOZZUg44rduQdE+WLmsBjc1i/U/AvckVUbTVi/PxVvbT+H3NIq83EfN1f0i/BBdkklzmaVmo8He2nx8ewh6B3m+K2R1HZCCLz+01m8s/O8+ZibqxqDY/yQU1KFs1klFq0Mo3oG4elbeznNz09KXjn+/XsyNuxPRVFFtfm4Rq1Ct2BP8++c6AB3eGhcoHVVQeuihtEkkF+mR0G5HnmlVbiYW4ZzWaU4n11abwNdAOgS6IGhXS63Tiuxf1xd8h+Zh1MLcCSlEIdTC5FTUtXgud46F3QJ8kSErxvCfXUI83FDkJcWXjoXeLu5wlvnAhdVzRAFSQKEACqqjSjXG1BR+16VV6ZHXmnNvcoqrkJGcSUyiyqQXVKFm7q1/5YvDhNu9Ho93N3dsWnTJkyaNMl8PD4+HoWFhfjmm2/qXdOpUycsWLAATzzxhPnY4sWL8fXXX+OPP/6od35VVRWqqi5/c4uLixEVFdXu4ebXczmY/uG+dnu8xnhpXdA12BM9Q7xwTaQP+kX4oFeol12PJwEujwOS/+o5mFyAtMKKVj+eq1rC0M4BGN07GKN7haBTgH0NeGwvZVUGrNuThO2nsnAivRh6w+XuRkkCro30xc3dAzFtaDRCfZy3r51aZ8epLPxwNAO7z+VahGSgJhBf18kPM2+MwQ2tGHzsCPQGE/Yl5mPH6SzsOJVt3ruqNTRqFa6J8Mag6Mutz44wviWruBLHLhXVDGxOL8aFnFIk55W1eFBzS/UK9cLWJ25u18dsSbhRdEey3NxcGI1GhIRYTjEMCQnB6dOnG7wmMzOzwfMzMzMbPH/58uVYunRp+xTchOHdg3D6hVuRWVSJjKJKZBRVIKOoEnmleuSXVZmnWlfojSjTG1CuN0JvMEGtkiABkCQJbhp1TWrWucJL54JQbx3Ca1N1hK87ugV7WrVryZpc1CpcE1Ezqyp+WAyAmj7r89klOJtVinNZpcgprUJB7X0qqqiGxkUFrYsKOtea+9It2NPcpNorzLvZg20dmYfWBXNHdcPcUd2gN5hwJrMEx9KK4OPmimFdAxy2+4BsY3TvEIzuHQKTSdT+cZGPIC8dBkT5dogwrHFR4abugbipeyAW3d4HqfkVOJNVgnPZJTiXVYpLBeWorDahstqISoMRakmCv0ft5Ap3DaL83dEjxBPdQ7wQ7e+u2Bo/bRHirUNIH53FVP7KaiMSc8uQlFtm7h1IL6xAfpkeJZUGFFdWo6TSAKNJwCSEuavLTaOGe+2Hp9YFAZ5aBHjULBES5KlFmK9bTc+FjxsCFP7d5PTvDgsXLsSCBQvMn8stN9agc1XXDBiz4eqdjszHzRWDov0xKJpjRJpD46JCv0gf9Iv0UboUcjAqlYQ+4d7oE+4c3U6tIUkSOgW4o1OAO8bYwZo9StK5qtE7zNtpuiEbomi4CQwMhFqtRlZWlsXxrKwshIY2PNsjNDS0RedrtVpotY45sJSIiIhaTtE2No1Gg0GDBmHHjh3mYyaTCTt27EBsbGyD18TGxlqcDwDbtm1r9HwiIiLqWBTvllqwYAHi4+MxePBgDBkyBCtXrkRZWRlmzZoFAJgxYwYiIiKwfPlyAMDjjz+OESNG4PXXX8eECROwfv16HDhwAGvWrFHyZRAREZGdUDzcTJkyBTk5OVi0aBEyMzMxYMAAbN261TxoOCUlBSrV5QamYcOG4bPPPsPf//53PPvss+jevTu+/vrrZq1xQ0RERM5P8XVubM3ZVigmIiLqCFry/u1489qIiIiImsBwQ0RERE6F4YaIiIicCsMNERERORWGGyIiInIqDDdERETkVBhuiIiIyKkw3BAREZFTYbghIiIip6L49gu2Ji/IXFxcrHAlRERE1Fzy+3ZzNlbocOGmpKQEABAVFaVwJURERNRSJSUl8PHxafKcDre3lMlkQnp6Ory8vCBJktLlKK64uBhRUVFITU3lXltWxPtsG7zPtsH7bDu815cJIVBSUoLw8HCLDbUb0uFablQqFSIjI5Uuw+54e3t3+H84tsD7bBu8z7bB+2w7vNc1rtZiI+OAYiIiInIqDDdERETkVBhuOjitVovFixdDq9UqXYpT4322Dd5n2+B9th3e69bpcAOKiYiIyLmx5YaIiIicCsMNERERORWGGyIiInIqDDdERETkVBhuOqD8/HxMmzYN3t7e8PX1xf3334/S0tJmXSuEwG233QZJkvD1119bt1AH19L7nJ+fj0cffRQ9e/aEm5sbOnXqhMceewxFRUU2rNr+rVq1CjExMdDpdBg6dCj27dvX5PkbN25Er169oNPp0K9fP2zZssVGlTq2ltzn999/H8OHD4efnx/8/PwQFxd31e8LXdbSn2nZ+vXrIUkSJk2aZN0CHRDDTQc0bdo0nDhxAtu2bcP333+P3bt3Y86cOc26duXKldy2oplaep/T09ORnp6O1157DcePH8e6deuwdetW3H///Tas2r5t2LABCxYswOLFi3Ho0CH0798f48aNQ3Z2doPn79mzB1OnTsX999+Pw4cPY9KkSZg0aRKOHz9u48odS0vv865duzB16lTs3LkTCQkJiIqKwtixY5GWlmbjyh1PS++1LCkpCU8++SSGDx9uo0odjKAO5eTJkwKA2L9/v/nYf//7XyFJkkhLS2vy2sOHD4uIiAiRkZEhAIjNmzdbuVrH1Zb7XNcXX3whNBqNqK6utkaZDmfIkCFi7ty55s+NRqMIDw8Xy5cvb/D8e++9V0yYMMHi2NChQ8VDDz1k1TodXUvv85UMBoPw8vISH3/8sbVKdBqtudcGg0EMGzZMfPDBByI+Pl7ceeedNqjUsbDlpoNJSEiAr68vBg8ebD4WFxcHlUqF33//vdHrysvLcd9992HVqlUIDQ21RakOrbX3+UpFRUXw9vaGi0uH2wauHr1ej4MHDyIuLs58TKVSIS4uDgkJCQ1ek5CQYHE+AIwbN67R86l19/lK5eXlqK6uhr+/v7XKdAqtvdfPP/88goOD2arbBP7G7GAyMzMRHBxscczFxQX+/v7IzMxs9Lr58+dj2LBhuPPOO61dolNo7X2uKzc3Fy+88EKzuwydXW5uLoxGI0JCQiyOh4SE4PTp0w1ek5mZ2eD5zf0edEStuc9X+tvf/obw8PB6wZIsteZe//bbb/jwww9x5MgRG1TouNhy4ySeeeYZSJLU5EdzfzFd6dtvv8XPP/+MlStXtm/RDsia97mu4uJiTJgwAX369MGSJUvaXjiRjbz88stYv349Nm/eDJ1Op3Q5TqWkpATTp0/H+++/j8DAQKXLsWtsuXESf/3rXzFz5swmz+nSpQtCQ0PrDVQzGAzIz89vtLvp559/xoULF+Dr62tx/O6778bw4cOxa9euNlTuWKx5n2UlJSW49dZb4eXlhc2bN8PV1bWtZTuFwMBAqNVqZGVlWRzPyspq9J6Ghoa26Hxq3X2Wvfbaa3j55Zexfft2XHvttdYs0ym09F5fuHABSUlJmDhxovmYyWQCUNMyfObMGXTt2tW6RTsKpQf9kG3JA10PHDhgPvbjjz82OdA1IyNDHDt2zOIDgHjrrbfExYsXbVW6Q2nNfRZCiKKiInHDDTeIESNGiLKyMluU6lCGDBki5s2bZ/7caDSKiIiIJgcU33777RbHYmNjOaD4Klp6n4UQ4pVXXhHe3t4iISHBFiU6jZbc64qKinq/i++8805xyy23iGPHjomqqipblm7XGG46oFtvvVUMHDhQ/P777+K3334T3bt3F1OnTjV//dKlS6Jnz57i999/b/QxwNlSV9XS+1xUVCSGDh0q+vXrJ86fPy8yMjLMHwaDQamXYVfWr18vtFqtWLdunTh58qSYM2eO8PX1FZmZmUIIIaZPny6eeeYZ8/n/+9//hIuLi3jttdfEqVOnxOLFi4Wrq6s4duyYUi/BIbT0Pr/88stCo9GITZs2WfzclpSUKPUSHEZL7/WVOFuqYQw3HVBeXp6YOnWq8PT0FN7e3mLWrFkWv4QSExMFALFz585GH4Ph5upaep937twpADT4kZiYqMyLsENvv/226NSpk9BoNGLIkCFi79695q+NGDFCxMfHW5z/xRdfiB49egiNRiP69u0rfvjhBxtX7Jhacp+jo6Mb/LldvHix7Qt3QC39ma6L4aZhkhBC2LorjIiIiMhaOFuKiIiInArDDRERETkVhhsiIiJyKgw3RERE5FQYboiIiMipMNwQERGRU2G4ISIiIqfCcENEREROheGGiIiInArDDRE5rNOnT2PUqFHQ6XTo0aMHtmzZAkmScOTIEaVLIyIFMdwQkUM6ffo0hg4diuHDh+PEiRN45ZVXMGPGDLi6uqJPnz5Kl0dECuLeUkTkkEaPHo3o6GisXbvWfOzuu+/GuXPncPToUQUrIyKluShdABFRSyUnJ+Pnn3/GH3/8YXFco9Ggf//+ClVFRPaC3VJE5HCOHDnSYPfT8ePHzeFmxIgRGDBgAAYMGAC1Wo0DBw4oUSoRKYAtN0TkcFQqFYxGI4xGI1xcan6Nbd261SLc/PLLLwCAxYsXY8SIERg8eLBi9RKRbbHlhogczqBBg+Dq6opnn30WFy9exJdffom5c+cCgEW31MqVK5GUlISVK1cqVCkRKYHhhogcTnh4OD744AN88cUX6N+/PzZs2IAHH3wQoaGhCA4OBgCsW7cOu3fvxtq1ayFJksIVE5EtcbYUETmFBQsW4OTJk9i6dSs2b96Mf/3rX/jmm2+g1WqVLo2IbIwtN0TkFI4ePWrukpo9ezYuXryIoUOHYsCAAfj+++8Vro6IbIkDionIKRw7dgyzZs0CABQUFChcDREpid1SRERE5FTYLUVEREROheGGiIiInArDDRERETkVhhsiIiJyKgw3RERE5FQYboiIiMipMNwQERGRU2G4ISIiIqfCcENEREROheGGiIiInArDDRERETmV/wfL5NAdW36nxQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Simulation1 = StackedTrapezoidSimulation(qys=qxs, qzs=qzs)\n", "\n", "intensity = Simulation1.simulate_diffraction(params=i_params)\n", "\n", "#plot\n", "import matplotlib.pyplot as plt\n", "plt.plot(qzs, intensity)\n", "plt.xlabel(r'$q_{z}$')\n", "plt.ylabel('Intensity')\n", "plt.title('Stacked Trapezoid diffraction simulation')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A txt or csv file with the experimental data containing the values for $Q_{x}$, $Q_{z}$ and intensities can be in the format as shown below:\n", " \n", " ```\n", " qx, qz, intensity\n", " 0.1, 0.2, 0.3\n", " 0.2, 0.3, 0.4\n", " 0.3, 0.4, 0.5\n", " ...\n", " ```\n", "\n", "for a text file you can read the data using the following code:\n", "\n", "```python\n", "import numpy as np\n", "data = np.genfromtxt('path_to_file.txt', delimiter=',', skip_header=1)\n", "qx = data[:,0]\n", "qz = data[:,1]\n", "intensity = data[:,2]\n", "```\n", "\n", "We will suppose that we've read the data and stored it in the variables qx, qz and intensity. And we'll use the simulated intensities used in the previous section." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's introduce a bit of noise to the simulated intensities to make it more realistic." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHACAYAAABeV0mSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABzPklEQVR4nO3dd3iV9f3/8ed9TvYmCRlAIOwhUxBERcWiiErVWmvVCnW2VmqVOtu6W1Hrauvg60Rb/eHGhRNBRBEEBEE2BFnZITs5Sc65f3/c55wkkEASkpyR1+O6zpXkPvd9532OmLzymYZpmiYiIiIiQcLm6wJERERE2pPCjYiIiAQVhRsREREJKgo3IiIiElQUbkRERCSoKNyIiIhIUFG4ERERkaCicCMiIiJBReFGREREgorCjYiIiASVLh1uli5dyvTp0+nRoweGYbBgwYJWXX/33XdjGMYhj+jo6I4pWERERI6oS4ebiooKRo0axZNPPtmm62+66Says7MbPYYNG8aFF17YzpWKiIhIS3XpcDNt2jT+/ve/c/755zf5vMPh4KabbqJnz55ER0czYcIElixZ4n0+JiaGtLQ07yM3N5eNGzdy5ZVXdtIrEBERkYN16XBzJLNmzWL58uXMnz+fH374gQsvvJAzzzyTbdu2NXn+c889x6BBg5g0aVInVyoiIiIeCjfN2L17Ny+++CJvvPEGkyZNon///tx0002cdNJJvPjii4ecX11dzSuvvKJWGxERER8L8XUB/mr9+vU4nU4GDRrU6LjD4SApKemQ89955x3KysqYOXNmZ5UoIiIiTVC4aUZ5eTl2u53Vq1djt9sbPRcTE3PI+c899xznnHMOqampnVWiiIiINEHhphljxozB6XSSl5d3xDE0WVlZLF68mPfee6+TqhMREZHmdOlwU15ezvbt271fZ2VlsXbtWhITExk0aBCXXnopM2bM4JFHHmHMmDHk5+ezaNEiRo4cydlnn+297oUXXiA9PZ1p06b54mWIiIhIA4Zpmqavi/CVJUuWMHny5EOOz5w5k3nz5lFbW8vf//53Xn75Zfbt20dycjLHH38899xzDyNGjADA5XLRp08fZsyYwT/+8Y/OfgkiIiJykC4dbkRERCT4aCq4iIiIBBWFGxEREQkqXW5AscvlYv/+/cTGxmIYhq/LERERkRYwTZOysjJ69OiBzXb4tpkuF272799PRkaGr8sQERGRNtizZw+9evU67DldLtzExsYC1psTFxfn42pERESkJUpLS8nIyPD+Hj+cLhduPF1RcXFxCjciIiIBpiVDSjSgWERERIKKwo2IiIgEFYUbERERCSpdbsyNiIj4D6fTSW1tra/LED8RFhZ2xGneLaFwIyIinc40TXJyciguLvZ1KeJHbDYbffv2JSws7Kjuo3AjIiKdzhNsUlJSiIqK0qKq4l1kNzs7m969ex/VvwmFGxER6VROp9MbbJKSknxdjviR7t27s3//furq6ggNDW3zfTSgWEREOpVnjE1UVJSPKxF/4+mOcjqdR3UfhRsREfEJdUXJwdrr34TCjYiIiAQVhRsREZEOdvfddzN69OhWXXPqqadyww03+LyOQKQBxSIiIh3spptu4o9//GOrrnn77bePalBtV6ZwIyJdR20VhESAxnpIJzFNE6fTSUxMDDExMa26NjExsYOqCn7qlhKRrqE0G/45AN75va8rkQDncDi4/vrrSUlJISIigpNOOonvvvsOgCVLlmAYBh999BFjx44lPDycZcuWHdIdVFdXx/XXX09CQgJJSUnceuutzJw5k/POO897zsHdUpmZmdx///1cccUVxMbG0rt3b5555plGtd16660MGjSIqKgo+vXrxx133NElV4BWuBGRriFvI9SUw55vfV2JNME0TSpr6nzyME2zVbXecsstvPXWW7z00kusWbOGAQMGMHXqVIqKirzn3HbbbTzwwANs2rSJkSNHHnKPBx98kFdeeYUXX3yRr7/+mtLSUhYsWHDE7/3II48wbtw4vv/+e/7whz9w7bXXsmXLFu/zsbGxzJs3j40bN/Kvf/2LZ599lscee6xVry8YqFtKRLqG2krro6Pct3VIk6pqnQy78xOffO+N904lKqxlvw4rKip4+umnmTdvHtOmTQPg2Wef5bPPPuP555/nuOOOA+Dee+/l9NNPb/Y+//nPf7j99ts5//zzAXjiiSdYuHDhEb//WWedxR/+8AfAaqV57LHHWLx4MYMHDwbgb3/7m/fczMxMbrrpJubPn88tt9zSotcXLBRuRKRrqK2yPtYo3Ejb7dixg9raWk488UTvsdDQUMaPH8+mTZu84WbcuHHN3qOkpITc3FzGjx/vPWa32xk7diwul+uw379hK5BhGKSlpZGXl+c99tprr/Hvf/+bHTt2UF5eTl1dHXFxca1+nYFO4UZEuoaaCutjXTU468CuH3/+JDLUzsZ7p/rse7e36Ojodr8ncMjsKcMwvIFo+fLlXHrppdxzzz1MnTqV+Ph45s+fzyOPPNIhtfgz/d8tIl2Dp+UGoKYMIrv5rhY5hGEYLe4a8qX+/fsTFhbG119/TZ8+fQBrO4nvvvuuxWvSxMfHk5qaynfffcfJJ58MWNsNrFmz5qjWoPnmm2/o06cPf/3rX73HfvrppzbfL5D5/78kEZH2UFtR/3lNhcKNtEl0dDTXXnstN998M4mJifTu3ZuHHnqIyspKrrzyStatW9ei+/zxj39kzpw5DBgwgCFDhvCf//yHAwcOHNX2AwMHDmT37t3Mnz+f4447jg8//JB33nmnzfcLZAo3ItI1NGy50aBiOQoPPPAALpeLyy67jLKyMsaNG8cnn3xCt24tD8y33norOTk5zJgxA7vdzjXXXMPUqVOx29veRfbzn/+cG2+8kVmzZuFwODj77LO54447uPvuu9t8z0BlmK2dAxfgSktLiY+Pp6SkpEsOshLpsj66DVY8bX1+1SLo1fyAT+lY1dXVZGVl0bdvXyIiInxdjl9wuVwMHTqUX/3qV9x3332+LsdnDvdvozW/v9VyIyJdg2cqOICjzHd1iGCNhfn000855ZRTcDgcPPHEE2RlZXHJJZf4urSgoEX8RKRraBhuNB1cfMxmszFv3jyOO+44TjzxRNavX8/nn3/O0KFDfV1aUFDLjYh0DRpzI34kIyODr7/+2tdlBC213IhIYGnrMMGahrOlFG5EgpnCjYgEjtyN8Nhw+Pzu1l/bqOVGY25EgpnCjYgEBmctvPM7KN0LG95u/fUHr3MjIkFL4UZEAsOyxyDnB+vz0v1whD14DtFohWJ1S4kEM4UbEfF/2T/Alw/Wf+2qhYr81t2jpuFUcIUbkWCmcCMi/q2uBhb8AVx1MOQciO1hHS/Z27r7NJoKrjE3IsFM4UZE/NtXj0DueohMhHMeg/he1vHSowg3armRTrRr1y4Mw2Dt2rXNnrNkyRIMw6C4uPiovldmZiaPP/74Ud3jSObNm0dCQkKHfo+jpXAjIv6rzmGNtQE4+2GISYH4ntbXJftafh9nHThr6r/WmBuRNrvooovYunWrr8s4LC3iJyL+q+oAOB1g2GDY+daxOHe4KW1FuGnYagNquRE5CpGRkURGRvq6jMNSy42I+K/qEutjeBzY3D+uPN1SJXtafp+GM6VAY26kzT7++GNOOukkEhISSEpK4pxzzmHHjh2Nzlm5ciVjxowhIiKCcePG8f333x9yn4ULFzJo0CAiIyOZPHkyu3btOuScZcuWMWnSJCIjI8nIyOD666+noqJ+GYO8vDymT59OZGQkffv25ZVXXjli/b/97W8577zzePjhh0lPTycpKYnrrruO2tpa7zkHDhxgxowZdOvWjaioKKZNm8a2bdu8zx/cLbVu3TomT55MbGwscXFxjB07llWrVrX4dXQEhRsR8V9VxdbHyIT6Y95w05qWm4N+kGqdG/9jmtZ/F188WrHqdUVFBbNnz2bVqlUsWrQIm83G+eefj8u9NEF5eTnnnHMOw4YNY/Xq1dx9993cdNNNje6xZ88efvGLXzB9+nTWrl3LVVddxW233dbonB07dnDmmWdywQUX8MMPP/Daa6+xbNkyZs2a5T3nt7/9LXv27GHx4sW8+eabPPXUU+Tl5R3xNSxevJgdO3awePFiXnrpJebNm8e8efMa3XfVqlW89957LF++HNM0OeussxoFoIYuvfRSevXqxXfffcfq1au57bbbCA0NbfHr6AjqlhIR/+VpuYmIrz/Wpm6pg1pu1C3lf2or4f4evvnef9kPYdEtOvWCCy5o9PULL7xA9+7d2bhxI8OHD+fVV1/F5XLx/PPPExERwTHHHMPevXu59tprvdc8/fTT9O/fn0ceeQSAwYMHs379eh58sH65gzlz5nDppZdyww03ADBw4ED+/e9/c8opp/D000+ze/duPvroI1auXMlxxx0HwPPPP9+ijTe7devGE088gd1uZ8iQIZx99tksWrSIq6++mm3btvHee+/x9ddfc8IJJwDwyiuvkJGRwYIFC7jwwgsPud/u3bu5+eabGTJkiLfWlr6OiIiII9bbFmq5ERH/1VS48bTclOVY08RbwhNuwt33cTqsFY9FWmnbtm1cfPHF9OvXj7i4ODIzMwHrFzzApk2bGDlyZKNf2hMnTmx0j02bNjFhwoRGxw4+Z926dcybN4+YmBjvY+rUqbhcLrKysti0aRMhISGMHTvWe82QIUNaNIvpmGOOwW63e79OT0/3tvh47tuwvqSkJAYPHsymTZuavN/s2bO56qqrmDJlCg888ECjbrojvY6OopYbEfFf1cXWx4bhJioZ7GHW7KeybOjW58j38XRDxXQHhzswOcogKrFdy5WjEBpltaD46nu30PTp0+nTpw/PPvssPXr0wOVyMXz4cGpqWhi0W6i8vJzf/e53XH/99Yc817t376OareTpMvIwDMPbrdYWd999N5dccgkffvghH330EXfddRfz58/n/PPPP+Lr6CgKNyLiv7zhJqH+mM1mdU0dyLK6ploSbjwtNxHxYA+3Wm5qyhVu/IlhtLhryFcKCwvZsmULzz77LJMmTQKswbINDR06lP/+979UV1d7W2++/fbbQ8557733Gh07+Jxjjz2WjRs3MmDAgCZrGTJkCHV1daxevdrbLbVly5ajXidn6NCh1NXVsWLFCm+3lOd1Dxs2rNnrBg0axKBBg7jxxhu5+OKLefHFFzn//POP+Do6irqlRMR/NdUtBQ0GFbdwIT/PgOLQKAiPsT7XuBtppW7dupGUlMQzzzzD9u3b+eKLL5g9e3ajcy655BIMw+Dqq69m48aNLFy4kIcffrjROb///e/Ztm0bN998M1u2bOHVV19tNKAX4NZbb+Wbb75h1qxZrF27lm3btvHuu+96B+IOHjyYM888k9/97nesWLGC1atXc9VVVx31FO2BAwdy7rnncvXVV7Ns2TLWrVvHb37zG3r27Mm55557yPlVVVXMmjWLJUuW8NNPP/H111/z3Xffecf+HOl1dBSFGxHxX95wk9D4uGdQcYvDjbvlJjQKwtzhRgv5SSvZbDbmz5/P6tWrGT58ODfeeCP//Oc/G50TExPD+++/z/r16xkzZgx//etfGw0UBqs75q233mLBggWMGjWKuXPncv/99zc6Z+TIkXz55Zds3bqVSZMmMWbMGO6880569KgfdP3iiy/So0cPTjnlFH7xi19wzTXXkJKSctSv88UXX2Ts2LGcc845TJw4EdM0Wbhw4SHdWQB2u53CwkJmzJjBoEGD+NWvfsW0adO45557Wvw6OoJhmq2YA9fOli5dyj//+U9Wr15NdnY277zzDuedd16z57/99ts8/fTTrF27FofDwTHHHMPdd9/N1KlTW/w9S0tLiY+Pp6SkhLi4uHZ4FSLSYV6fARvfhbMehvFX1x9fdB989TAcdxWc/ciR77PiGfjoZhh2HhRuh9wN8Ju3YcDPOqx0aV51dTVZWVn07du3w2bLSGA63L+N1vz+9mnLTUVFBaNGjeLJJ59s0flLly7l9NNPZ+HChaxevZrJkyczffr0JhdIEpEg4Fnn5pBuqVZuweBZoTgsWi03Il2ATwcUT5s2jWnTprX4/IM3A7v//vt59913ef/99xkzZkw7VyciPtfcmJu41o65cYeb0Mj6MTdayE8kaAX0bCmXy0VZWRmJic3PeHA4HDgcDu/XpaWlnVGaiLSHZgcUexbya224iaqfkaMBxSJBK6AHFD/88MOUl5fzq1/9qtlz5syZQ3x8vPeRkZHRiRWKyFFpbkCxZ7ZU1YGWtcDUNAw3se5j2l9KJFgFbLh59dVXueeee3j99dcPOzr89ttvp6SkxPvYs6cVm+2JiO+YZvMtNxHx9SGlJeNuPLOlwjQV3J/4cD6L+Kn2+jcRkN1S8+fP56qrruKNN95gypQphz03PDyc8PDwTqpMRNpNTTmYTuvzg8MNWF1T+Zutrqnugw5/r4br3GhAsc95phRXVlYe9bosElw8Kz033B6iLQIu3Py///f/uOKKK5g/fz5nn322r8sRkY7iabWxhVoDgQ8W38sKN61pudEifn7BbreTkJDg3c8oKioKwzB8XJX4msvlIj8/n6ioKEJCji6e+DTclJeXs337du/XWVlZrF27lsTERHr37s3tt9/Ovn37ePnllwGrK2rmzJn861//YsKECeTk5AAQGRlJfHwTf9mJSODyTAOPTLCW5j9YEwv5ffhDNt/tKuJ3p/QjPb5BIKppMFvK23KjMTe+lJaWBuANOCJgLZTYu3fvow67Pg03q1atYvLkyd6vPctYz5w5k3nz5pGdne3daRXgmWeeoa6ujuuuu47rrrvOe9xzvogEkebG23h4BhWX7qWsupY73/2Rd763WnEWrN3HA78YyZnDrV+gjda5cbo3OFTLjU8ZhkF6ejopKSnU1mqHdrGEhYVhsx39cGCfhptTTz31sIOHDg4sS5Ys6diCRMR/HCncuFtuSnN/4ux/f8WeoipsBvRJiiaroILf/281l0zozR1nDyOy4To3TvcvUq1z4xfsdvtRj68QOVjAjbkRkS6iuWngHu6Wm/x9O9jjqKJnQiT/+vVoRvZK4JFPt/B/S3fy6ordrNtTzAfOSgyA0GgwXdb1GlAsErQCdiq4iAS56mLrYzMtN65Yq+UmnUKmDkth4Z8mMS4zkbAQG7efNZT/XTmBhKhQftxfiqPKHWQajrlRt5RI0FK4ERH/dIRuqcXZ1nTiKMPBP6f3IT6y8Y7FJw1M5rpTBwBgegYUh0VpQLFIF6BwIyL+6TDhxjRNnli2l0LTWsgvzpHb5C0untCb+IgQwk33FiyaCi7SJSjciIh/8oSbyIRDnlqZVcT3u4vJJtk60MwGmjHhIVw5IQ2bYU1cMBt2S7lqoc7R5HUiEtgUbkTEP3nWuWmi5eapJTsAsMUfeXfw34zt7v38612V9eEG1HojEqQUbkTEPzXTLfXj/hK+3JqPzYCMzIHWwdLmVylODLO2cKg2Q3n6qyywh0CIe4E/jbsRCUoKNyLin5oJN0+7W23OGdmD2KR062BlUfP3ca9xU0U4X28vZN2e4vpxN1rrRiQoKdyIiH/yhptu3kN7iipZuD4bgGtP7V8ffDznNsW7gF8UAHO/3KHp4CJBTuFGRPxTE+vcfLQhG5cJJ/RPYmh6XMvCjXsaeFS0NbPqi815uDQdXCSoKdyIiP9xOcFRan3eINws3pwPwBnDUhs/d9iWG2tH8LCoGHomROKoc1HqDLOeU8uNSFBSuBER/+MJNuANMGXVtXy3yxpbc+rglEbPHT7cWONqjNAoJg+xZk7lOtwL/mkLBpGgpHAjIv7HE1ZCoyDEamVZtq2AOpdJv+RoMpOjredb0XJDaBST3aFob4V7o0a13IgEJYUbEfE/Taxxs3hLHgCTh6TUn9cw3Jhm0/fyzIgKjeSE/smEhdjIr/G03GjMjUgwUrgREf9z0DRwl8tk8RZrvI2n9aXh87hq61toDuY5HhZNZJidif2SqMC9zo1abkSCksKNiPgfb7hJAGBjdin5ZQ6iw+wc17d+ajhhMWDYGl9zMO9UcCvQnDYkhXIirGNa50YkKCnciIj/OWga+OLNVpfUiQOSCQ+x159nGEced+MNN9Y4ncmDU6gwrXBTU1Xa9DUiEtAUbkTE/xzULfVFU+NtPI4Ubmoat9z0TooiMsa6pqjoMCsbi0jAUrgREf/TINwUVdSwdk8xcNB4G4+WttyERXkP9e2RBkBpyYH2qFZE/IzCjYj4H09QiUzgy615mCYMTY8jLT7i0HNb3C1VH24G97bCjaOiBJermVlWIhKwFG5ExP80mAruWZV48uDuTZ/rDTfFTT/vXecm0ntoQC9rw81wVyUb9jcIRc1NJxeRgKJwIyL+x90K4wqPY+k2K9yc1tR4G2jBmBvPOjfR3kOhUdY+U9FGtTc88eVD8M/+kL/16GoXEZ9TuBER/+MOKtmOCIora4kKszM6I6Hpc93TxZvvljq05YYwK9zEUMWqn4rA5YKVz0JlIexZcfT1i4hPKdyIiP9xB5UtxdaPqFG9EgixN/Pjqg0Digm3dgWPwsHa3Qdw5WyACmtGFlUaZCwS6BRuRMT/uMfPbHDP1B7TO6H5c9swoJgwK9yEGk5qHFUUrf/okO8tIoFL4UZE/I87qKzOtQb4NtslBa1Y5+bQcAMQTTXOrYvqn1PLjUjAU7gREf9SV+NtbVlb4A43R9VyU78ruJfN5h1gnGyUkFS0pv45z0wtEQlYCjci4l8c9VsilJlR9EyIJCW2ifVtPDzhxtHMVgpNjbkB77ibn9m+J8SsrT+ubimRgKdwIyL+xd1y4rBH48J2+PE2cPiWG2ettWM4NJ4tBRBmtdxMs1uzo1yx6e7vr24pkUCncCMi/sUdUsqxwsdhx9tA43Bz8CJ8nlYbaLTODeAddzPSlgVAbq9p1nF1S4kEPIUbEfEv7m6hAqfVjTSmd7fDn+8JN84aqKtu/JxnMLFhB3to4+fCY72f1pk2vo2c1Oj7i0jgUrgREf/ibrk54Iok1G5wTI+4w58fFgOGrdG1Xt7xNtFgGIde57bWHMA3he6vq4qtRf1EJGAp3IiIf3G3nJSaUQxLjyMi1H748w2j+XE33jVuDhpvA94BxQBfOUfwzT6n+yuz+cHJIhIQFG5ExL+4A0qJGX3kLimPZsNNE9PAPRq03CxjFPvKXbhC3CFIXVMiAU3hRkT8izuglBJ95MHEHs2FG++mmU2EG8+Ym4h4alJGAeAIdXeBacaUSEBTuBERv+KstIJFqRl15GngHkdquTl4jRuASPe9+53KqD5JAJS5Z2hpxpRIYFO4ERG/UnqgEIC6sDh6JzYRSpriDTfFjY8fbszNqEtg3JUw+a+MybC6vwrq3N9PLTciAS3E1wWIiDRUUVpIN6BbUneMg2c4NeeIA4oPWuMGIL4nnPMoAKMpB2C/I4JhNjTmRiTAqeVGRPxKXYXVapKWktryi8KbG3NzmJabBvomRRMfGcoBl6flprjl31tE/I5Pw83SpUuZPn06PXr0wDAMFixYcMRrlixZwrHHHkt4eDgDBgxg3rx5HV6niHSeWEcOAD1692/5RUdquWlqzE0DNpvByF7xlHjH3KhbSiSQ+TTcVFRUMGrUKJ588skWnZ+VlcXZZ5/N5MmTWbt2LTfccANXXXUVn3zySQdXKiKdoaqinCTTCha9+w1p+YVH7JY68tidYelxFJvu6eHqlhIJaD4dczNt2jSmTZvW4vPnzp1L3759eeSRRwAYOnQoy5Yt47HHHmPq1KkdVaaIdJKfdm5mCFBOJMnd01p+YVvWuTnIkPRY1jQ3W8rlgs/vgozxMHR6y+sSEZ8IqDE3y5cvZ8qUKY2OTZ06leXLlzd7jcPhoLS0tNFDRPxT7u6tABSFph26XcLhtGWdm4MMSYujxLTCjXlwt9TelfDNv+GzO1tek4j4TECFm5ycHFJTGw8yTE1NpbS0lKqqqiavmTNnDvHx8d5HRkZGZ5QqIm1QlrMDgKroXq27sC3r3Bykf/cYyg2rW6q24qBwU7LX+qiBxiIBIaDCTVvcfvvtlJSUeB979uzxdUki0gzXgZ8AsHXr07oLmw03npabw8+WAggLsRHbrTsAzoPDTek+66OnJUhE/FpArXOTlpZGbm5uo2O5ubnExcURGdn0D6/w8HDCw8M7ozwROQqmaRJRbrWQxKT2a93FDcONadZ3aXnH3DSxzk0TuqekQTnYHcWNnyjNtj46HeCsBXto6+oTkU4VUC03EydOZNGiRY2OffbZZ0ycONFHFYlIe8ktdZDisv54Seo1qHUXe8KNswbqquuPt3CdG49e6ekAhDkrrBDj4Wm5Aagpb11tItLpfBpuysvLWbt2LWvXrgWsqd5r165l9+7dgNWlNGPGDO/5v//979m5cye33HILmzdv5qmnnuL111/nxhtv9EX5ItKONuWUkmHkAxCWnNm6i8NiwHD/OGvYNeVd56ZlLTd9e/Wo/6Lhfcqy6z93KNyI+DufhptVq1YxZswYxowZA8Ds2bMZM2YMd95pzUjIzs72Bh2Avn378uGHH/LZZ58xatQoHnnkEZ577jlNAxcJAtv35JJklFlfJPRu3cU2G4S7d/SubjAj8nB7SzVhSM9ESk3rXEdZYf0TpQ3CjcbdiPg9n465OfXUUzFNs9nnm1p9+NRTT+X777/vwKpExBcK91nTwKtD4ojwdDO1RkS8tfheo5ablq9zA5ASG062EUMcVezNzqZ/2mBwORu33KhbSsTvBdSYGxEJXlV5WQDUxLZxuYamZky1Yp0bAMMwqAmxWoCys92BpiIfTGeDeyrciPg7hRsR8bmaOhe2EqsLOjQps2038Yab4vpjrVjnxsMVmQBAfr61x1WjwcSgMTciAUDhRkR8bkd+OT3JAyCie9+23eTglpvaKqhr3VRw69REAEoOFFgHGo63AY25EQkACjci0vkKdzQKCZtzSullWGHCSGjlAn4eEQnWR0+42bPC+hjbA6ISW3ybqLgkAKpKrZlbjcbbANSUta0+Eek0Cjci0rmyf4D/jIUF13oPbc4uI8OwWm5o7erEHge33Oz80vrY9+RW7VMVl5gCQGhNKflljkO7pdRyI+L3FG5EpHPtWw2YsOUjb1DYlFNGL/caN62eBu5xcLjJWmp97HdKq24TGt3NKsOoYHNO6aHdUhpzI+L3FG5EpHN5unmcNbBrGQB792cTb7jXpGmPcFNdAvvXWF/3Pbl194m0wk08FWzOLqtvuYlzb+aplhsRv6dwIyKdq2E3z/ZFFJY7iKiw9pQyo5JbvJrwIRqGm5++AdMFif0hvpU7jLtnS8UZFWzKKa0PY8kDrY8acyPi9xRuRKRzNezm2bGILTll3m0XjLaOt4HG4cYz3qaVXVLWfRIASKCczfsbdEt5w41abkT8ncKNiHSu0v31nxduZ0/W5qMfbwONw01Wg8HEreXpljIqyM3Pg1p3mElyhxuNuRHxez7dfkFEuqAyd7iJToGKPEKzFjcIN+3QclOyt359m8y2hJsEwBpzk+hy7y8V2Q2ik63P1XIj4vfUciMinaemon4206hfA9Cz8Jv6aeDt0XLjCTapIyA6qQ33SQAg0qgh03CvUhzbw9p5HDTmRiQAKNyISOfxjF8Ji4FjzgNgmON7+npCRHuMufFoy3gbsHYXN6wfjYMMa6AzcT0g3BNu1HIj4u8UbkSk83hmSsWmQ/poXJGJxFJFf5s79CRktv3eYTHeUAK0bbwNgM3mDUpDbNZ+V8Sl18/i0pgbEb+ncCMincczrTquB9jsFKWe0Pj51k7bbshms1pdAAw79Dnh8Ocfjrtrqr7lpmeDbim13Ij4O4UbEek8nplScT0A2BY7of652HQIjTi6+3u6pnqOhfDYtt/HPWOqn2GFMTM2vUG4KQfTPJoqRaSDKdyISOfxhJvYdACWM6r+uaMZTOzhCTdtHW/j4Z4xFWo4AThgT6ofc4MJtZVHd38R6VAKNyLSeRp2SwHrSiLZ5Mqwjh3NNHCPnsdaXVJDf3509/HsMO62q6YbhEYB7g04Ne5GxK8p3IhI5/Hu02SFm+155XzkdHdN9Tz26O9/9mNw83ZIH3l093F3S3lsroy2dhZv2DUlIn5L4UZEOk9pfctNVY2TfcVVPOk8l9JfvwvHXXX097fZICrx6O/j7pYCqDLD2HTAbn3hmTGlcCPi1xRuRKRzOGuhPNf6PLYHO/KtgBAfHUnckFPBHuq72g7WoOUmx+zG9nz3DCnPuBt1S4n4NYUbEekc5bmACbYQiO7uDTf9u7dxF/CO1GDMTS6J3lrrW240HVzEnynciEjn8HRJxaaDzcaOPCswDEiJOcxFPtKgWyrbTCSvzEFpdS2EuaeXawsGEb+mcCMinePgwcTelht/DDf13VJlod0Ba/CzWm5EAoPCjYh0jrIGLTfAjjwrIPT3x5abBt1SpjuM7cgr15gbkQChcCMinaNBy02d00VWgRVuBvhly02C99Owbj0Bd0uTWm5EAoLCjYh0jgbTwPceqKLG6SI8xEbPhEjf1tWUBt1S8SnWysk78io05kYkQCjciEjnaLD1wnb3YOJ+3WOw2QwfFtWM0Cgr4Bh2kjIGA1gzptRyIxIQQnxdgIh0EWWeTTN7smOXH8+UAms14kteh+oSMlP7ANvZXVRJXWiU9UNTY25E/JrCjYh0PNNs0C2Vzva8EsBP17jxyBgPQIppEhMeQrmjjsKaMFJBKxSL+Dl1S4lIx6ssAqfD+jw23bsont+23DRgGIZ3RldOtfvvQYUbEb+mcCMiHc/TJRWVjGkP84658cs1bprgaWHaW+7+kakxNyJ+TeFGRDqeZzBxXDoF5TWUVtdhGNA32Y+7pRrwtDDtKncPftaYGxG/pnAjIh2vtH4wsafVJqNbFBGhdh8W1XKeFqYdJe4DarkR8WsKNyLS8RpMA/cs3tfPnwcTH8QTbrYXuw9onRsRv6ZwIyIdr8E08KwC9xo3yYEx3gagd2IUNgMKakKtA2q5EfFrCjci0vEaTAPPKqgEoG9ylA8Lap2wEBu9ukVRToR1wFkDdTW+LUpEmqVwIyIdzzvmpoe35aZvALXcgDX4udITbkDTwUX8mMKNiHQ8d7dUXXQau4vcLTcBNOYGrHBTRwh1Rph1QOFGxG/5PNw8+eSTZGZmEhERwYQJE1i5cuVhz3/88ccZPHgwkZGRZGRkcOONN1JdXd1J1YpIq9VUQLU1zSjblUit0yQ8xEZ6XMQRLvQvnmnrVYa7bo27EfFbPg03r732GrNnz+auu+5izZo1jBo1iqlTp5KXl9fk+a+++iq33XYbd911F5s2beL555/ntdde4y9/+UsnVy4iLVaea30MjWJHqbVOTGZStH9umHkYnnBTbrrDjda6EfFbPg03jz76KFdffTWXX345w4YNY+7cuURFRfHCCy80ef4333zDiSeeyCWXXEJmZiZnnHEGF1988RFbe0TEh8rdf6zEpJBV6BlMHFhdUlBfc6kr3DqgbikRv+WzcFNTU8Pq1auZMmVKfTE2G1OmTGH58uVNXnPCCSewevVqb5jZuXMnCxcu5KyzzuqUmkWkDTwtNzGp7HKvcZMZgOGmR0IkYXZbfcuNwo2I3/LZruAFBQU4nU5SU1MbHU9NTWXz5s1NXnPJJZdQUFDASSedhGma1NXV8fvf//6w3VIOhwOHw+H9urS0tH1egIi0jKflJro7Oz0L+AVguLHbDPokRVFxQGNuRPydzwcUt8aSJUu4//77eeqpp1izZg1vv/02H374Iffdd1+z18yZM4f4+HjvIyMjoxMrFpH6bqlU7+rEgTZTyqNvcjQVnungDq1SLOKvfBZukpOTsdvt5ObmNjqem5tLWlpak9fccccdXHbZZVx11VWMGDGC888/n/vvv585c+bgcrmavOb222+npKTE+9izZ0+7vxYROQx3t1RdVHf2FVcB1oDiQNRorRu13Ij4LZ+Fm7CwMMaOHcuiRYu8x1wuF4sWLWLixIlNXlNZWYnN1rhku93aeM80zSavCQ8PJy4urtFDRDqRu+Wm0EjANCE2PITkmDAfF9U2fZOjNeZGJAD4bMwNwOzZs5k5cybjxo1j/PjxPP7441RUVHD55ZcDMGPGDHr27MmcOXMAmD59Oo8++ihjxoxhwoQJbN++nTvuuIPp06d7Q46I+JkKK9zsr40FrC4pwwisaeAefZOjWaOWGxG/59Nwc9FFF5Gfn8+dd95JTk4Oo0eP5uOPP/YOMt69e3ejlpq//e1vGIbB3/72N/bt20f37t2ZPn06//jHP3z1EkTkSNwtN7uqre0WArVLCqxws9SMBMBZXYb+pBLxTz4NNwCzZs1i1qxZTT63ZMmSRl+HhIRw1113cdddd3VCZSJy1EzTO+ZmS0UUUBuQa9x4dI8Np85uhZvK8hJifVyPiDQtoGZLiUiAqS62dtAGNpRYi9/1C9CZUgCGYRAVGw9AdUWJj6sRkeYo3IhIxynPtz6Gx7O1sBYI7G4pgJi4BABqqzQVXMRftSnczJw5k6VLl7Z3LSISbNxdUq7o7uSVWYtpBuLqxA11S0gEwKzWbCkRf9WmcFNSUsKUKVMYOHAg999/P/v27WvvukQkGLjDTVV4MgBJ0WHER4b6sqKjlpRohRtbrWZLifirNoWbBQsWsG/fPq699lpee+01MjMzmTZtGm+++Sa1tbXtXaOIBCr3TKliWzcgMDfMPFhqchIAoU6FGxF/1eYxN927d2f27NmsW7eOFStWMGDAAC677DJ69OjBjTfeyLZt29qzThEJRO41bvLMBCA4wk2PFKsVKsKspsJR5+NqRKQpRz2gODs7m88++4zPPvsMu93OWWedxfr16xk2bBiPPfZYe9QoIoGq3LOAn3uNmyAIN3HxVrdUjFHNrgINKhbxR20KN7W1tbz11lucc8459OnThzfeeIMbbriB/fv389JLL/H555/z+uuvc++997Z3vSISSNxjbnZUW6EmEHcDP0RY/WvYnVvgw0JEpDltWsQvPT0dl8vFxRdfzMqVKxk9evQh50yePJmEhISjLE9EAppnAb/yKCBwdwNvJDQKFwY2TLJzC4CBvq5IRA7SpnDz2GOPceGFFxIREdHsOQkJCWRlZbW5MBEJAu51bjxbL/RJDIJwYxjU2aMIc1aQU1Do62pEpAlt6pZavHhxk7OiKioquOKKK466KBEJAi4nVFjhJt+MJz0+gsiw4NiNyRVqhbSiA0U+rkREmtKmcPPSSy9RVVV1yPGqqipefvnloy5KRIJAZRGYTkwMioilT1KUrytqN0a41RJVXHzAx5WISFNa1S1VWlqKaZqYpklZWVmjbimn08nChQtJSUlp9yJFJAB5FvALTaCuOiQopoF7hETGQom1M3hpdS1xEYG9MKFIsGlVuElISMAwDAzDYNCgQYc8bxgG99xzT7sVJyIBzB1uPAv49QnwPaUasodb+4HHUMVPBZWM6BXv44pEpKFWhZvFixdjmiannXYab731FonuZcgBwsLC6NOnDz169Gj3IkUkALnH2+S54oDA3zCzEXe3VJThYFdhhcKNiJ9pVbg55ZRTAMjKyqJ3794YhtEhRYlIEHC33OyttVo5MpODZ8wNYVa4iaGKXQXahkHE37Q43Pzwww8MHz4cm81GSUkJ69evb/bckSNHtktxIhLA3KsT76uzWm6CYhq4h3shvyiq2VVY6eNiRORgLQ43o0ePJicnh5SUFEaPHo1hGJimech5hmHgdDrbtUgRCUDulpt8M4G0uOCZBg6Ae8xNtFHNrkK13Ij4mxaHm6ysLLp37+79XETksNwtN/lmfFBNAwe8LTexVPGTwo2I32lxuOnTp0+Tn4uINMkdbgqID6pp4AAkWD8Dz7V/zYsVUymrriVW08FF/EabF/H78MMPvV/fcsstJCQkcMIJJ/DTTz+1W3EiEsAadEsF0zRwAEb+CvqcRIxRzTOhj7Jnf7avKxKRBtoUbu6//34iIyMBWL58OU888QQPPfQQycnJ3Hjjje1aoIgEIGctVFlbE+Sb8fQNpplSAPZQ+NVL5NtS6GfLIemT66ztJkTEL7Qp3OzZs4cBAwYAsGDBAn75y19yzTXXMGfOHL766qt2LVBEApB7jZta7BQTE3wtNwDRyfy3zz+oNkNJzV0Ki+/3dUUi4tamcBMTE0NhobUb7qeffsrpp58OQERERJN7TolIF+Pukio04zCxBd+AYrfQXmO4tfZq64uvHoY9K31bkIgAbQw3p59+OldddRVXXXUVW7du5ayzzgLgxx9/JDMzsz3rE5FA1GCmVGpcOFFhrVovNGBkJkfzruskVoWNtw7sW+3bgkQEaGO4efLJJ5k4cSL5+fm89dZbJCUlAbB69Wouvvjidi1QRAJQMA8mbsCzpcTWWmuZDM/rFhHfatOfUwkJCTzxxBOHHNemmSICNGi5SaBvEIebPu6B0ntqYiEU7+sWEd9qc1txcXExK1euJC8vD5fL5T1uGAaXXXZZuxQnIgHKu8ZNnDcABKO4iFCSosPIr3ZvnKmWGxG/0KZw8/7773PppZdSXl5OXFxcow00FW5EpGG31IQgbrkB6JMURcFehRsRf9KmMTd//vOfueKKKygvL6e4uJgDBw54H0VFRe1do4gEmgbdUsE85gasQcX5ZoL1hbqlRPxCm8LNvn37uP7664mKCt7mZhFpO2eZ1YJRYMaTGcTdUmANKs4z3S03FQVazE/ED7Qp3EydOpVVq1a1dy0iEiTMSmsdLKKTgnYauEdmcjRFxOHCANMJlWq9FvG1Nv3UOfvss7n55pvZuHEjI0aMIDS08YZxP//5z9ulOBEJQM46QhzFAMR2S/NtLZ0gMykKJ3YOEEcSJda4m5juvi5LpEtrU7i5+mprRc577733kOcMw8DpVLOsSJdVdQAAl2mQlJLq42I6nmdMUZ4rniSbO9ww3LdFiXRxbeqWcrlczT4UbES6uMoCAIqJpndynI+L6XjxkaEkRoeR7xl3o0HFIj7XpnDTUHV1dXvUISLBwj3e5oAZS9/k4J4p5ZGZFEU+mg4u4i/aFG6cTif33XcfPXv2JCYmhp07dwJwxx138Pzzz7drgSISYCqslptC4oJ2w8yDZSY1mA7u3hFdRHynTeHmH//4B/PmzeOhhx4iLCzMe3z48OE899xz7VaciASeqhKrW+aAGevdeynYWWvdqOVGxF+0Kdy8/PLLPPPMM1x66aXY7Xbv8VGjRrF58+Z2K05EAk9xQTYAlSEJRIcH9zRwjz5JUQ0W8lO4EfG1Ni/iN2DAgEOOu1wuamtrj7ooEQlcVcXuX+5RSb4tpBP1TY4mnwTrCw0oFvG5NoWbYcOG8dVXXx1y/M0332TMmDGtuteTTz5JZmYmERERTJgwgZUrVx72/OLiYq677jrS09MJDw9n0KBBLFy4sFXfU0Q6Tm2ZNeYmJLbrrPXSJ6m+W8osU8uNiK+1qc34zjvvZObMmezbtw+Xy8Xbb7/Nli1bePnll/nggw9afJ/XXnuN2bNnM3fuXCZMmMDjjz/O1KlT2bJlCykpKYecX1NTw+mnn05KSgpvvvkmPXv25KeffiIhIaEtL0NEOoJ7KnhEwqH/Dwer+MhQnJHdwQVG9QGoq4GQsCNfKCIdok0tN+eeey7vv/8+n3/+OdHR0dx5551s2rSJ999/n9NPP73F93n00Ue5+uqrufzyyxk2bBhz584lKiqKF154ocnzX3jhBYqKiliwYAEnnngimZmZnHLKKYwaNaotL0NEOkCow1rELz4p3ceVdK7EpO7UmO4xiJoxJeJTbV7nZtKkSXz22Wfk5eVRWVnJsmXLOOOMM1p8fU1NDatXr2bKlCn1xdhsTJkyheXLlzd5zXvvvcfEiRO57rrrSE1NZfjw4dx///2HXTjQ4XBQWlra6CEiHSe6rhiA5JQevi2kk/XpHkuB1roR8QttCjf9+vWjsLDwkOPFxcX069evRfcoKCjA6XSSmtp4efbU1FRycnKavGbnzp28+eabOJ1OFi5cyB133MEjjzzC3//+92a/z5w5c4iPj/c+MjIyWlSfiLReSUUNCWYZAGnpPX1cTedqtNaNBhWL+FSbws2uXbuabC1xOBzs27fvqItqjsvlIiUlhWeeeYaxY8dy0UUX8de//pW5c+c2e83tt99OSUmJ97Fnz54Oq0+kq9udk0e4Yc2YjEoI/n2lGtJaNyL+o1UDit977z3v55988gnx8fHer51OJ4sWLSIzM7NF90pOTsZut5Ob2/iHQG5uLmlpTe8knJ6eTmhoaKO1dYYOHUpOTg41NTWNFhT0CA8PJzw8vEU1icjRycnZxwjAQTjhYV1jdWKPzKQoflTLjYhfaFW4Oe+88wBr5++ZM2c2ei40NJTMzEweeeSRFt0rLCyMsWPHsmjRIu99XS4XixYtYtasWU1ec+KJJ/Lqq6/icrmw2axGp61bt5Kent5ksBGRzlWQtx+AytAEutqfFJnJ0Sxxj7mpKclGP5FEfKdV3VKenb979+5NXl5eo93AHQ4HW7Zs4Zxzzmnx/WbPns2zzz7LSy+9xKZNm7j22mupqKjg8ssvB2DGjBncfvvt3vOvvfZaioqK+NOf/sTWrVv58MMPuf/++7nuuuta8zJEpIOUF1rj5WrDu/m4ks4XFxFKVai1cGHVgWwfVyPStbVpnZusrKx2+eYXXXQR+fn53HnnneTk5DB69Gg+/vhj7yDj3bt3e1toADIyMvjkk0+48cYbGTlyJD179uRPf/oTt956a7vUIyJHp9K9r5QRnezjSnzDFpsKpeAs1ZgbEV9q88YvixYtYtGiRd4WnIaaW6emKbNmzWq2G2rJkiWHHJs4cSLffvttq2oVkc7hLLfWdwntQqsTNxSRmA6lYKvUmBsRX2pTuLnnnnu49957GTduHOnp6RiG0d51iUiAKa6sIbK2GEIgulvXminlEZ/cC3ZBhKPA16WIdGltCjdz585l3rx5XHbZZe1dj4gEqF2FlSRirXHTVVtuktKtdbQiXFXgKIfwGB9XJNI1tWmdm5qaGk444YT2rkVEAthPhRUkGla46Uo7gjfUOzWFCtM9T6xCXVMivtKmcHPVVVfx6quvtnctIhLAsgoqSDTc25t00QHFfZKjKHAv5FepGVMiPtOmbqnq6mqeeeYZPv/8c0aOHEloaGij5x999NF2KU5EAseuggp+TtduuYmLCGWnrRt9yKMgew+9+/u6IpGuqU3h5ocffmD06NEAbNiwoT3rEZEAtauwkiRPy01U12y5AagOTwbHFkry9/q6FJEuq03hZvHixe1dh4gEuL0FJcQbldYXXbTlBqAuqjs4oFrdUiI+06pw84tf/OKI5xiGwVtvvdXmgkQk8BRX1mBUHYAIMA0bRmSCr0vymZC4NDgAzjIt5CfiK60KNw03yhQR8dhVWOkdTGxEdgOb/QhXBK/IbunwE4RoIT8Rn2lVuHnxxRc7qg4RCWC7ChpOA++6420A4rv3AiDSUejjSkS6rjZNBRcRaSiroMK7gF9XnQbu0d29kF+8WUxZda2PqxHpmhRuROSo7SqsoJu35SbRt8X4WHRiDwC6U8yu/AofVyPSNSnciMhRyyqoIAlNAwcgJgWAcKOOvTmaMSXiCwo3InJUTNMkK79hy03XnQYOQEg4VTZrT6m8bK11I+ILCjciclQKymsoc9TVL+DXxcfcANSGxgJQWJjv40pEuiaFGxE5KlkF1riS9FAt4OdhRMQBUHxAM6ZEfEHhRkSOSlZBOQApduujwg2ERCUAUF5ahGmavi1GpAtSuBGRo7LTPSOom3dAscJNWEw3AEJryygor/FxNSJdj8KNiByVnQUVgElUXbF1QGNusEdaq7nHUuntthORzqNwIyJHJauggliqsJt11gG13EC4NeYm1qj0dtuJSOdRuBGRNnO6TH5quIBfaDSERvq2KH/gHlAcR6W7ZUtEOpPCjYi02b4DVdQ6TVI9g4mj1WoDQITVLRVnVJKlVYpFOp3CjYi02U53l8vQOPegWXVJWTzdUhpzI+ITCjci0maeX9z9o6utA1196wUPT7eUUclPhZU4XZoOLtKZFG5EpM0808B7R1ZZB9RyY/F2S1VR43Sx70CVjwsS6VoUbkSkzbyrE4d4xtyo5QaAcCvcdLNZoWanZkyJdCqFGxFps6yCCkKpo2/up9aB5EG+LchfRHimglvhRuNuRDqXwo2ItEl1rZN9xVX80v4l4RX7ISYNRv7K12X5B/eA4iiXtcChwo1I51K4EZE22VVotdr8MfRd68BJN2qNGw/3mBsbTqJwKNyIdDKFGxFpk6z8Ci6wL6UHBRCTCmNn+rok/xEaCbYQwJoOvlNr3Yh0KoUbEWmTXXnFzApZYH2hVpvGDMPbNRVnVLK/pIrqWqePixLpOhRuRKRNkra/TS+jgIrQJBj7W1+X43/cXVNp4Q5M0+rGE5HOoXAjIq1XV8MpuS8BsGvo1Wq1aYp7xtTAOBeAtmEQ6UQKNyLSehsXkOrKJc9MwDbuCl9X45/c3VJ9Y6zd0rWBpkjnUbgRkVarXWW12vyvbgp90rQqcZPc3VK9omoBrXUj0pkUbkSkdYqyCN29DJdpsDT6dKLCQnxdkX/yjrmxNhVVuBHpPAo3ItI6a18FYJlrOFHdM31biz9zd0slhzoA2JGvLRhEOovCjYi0nMvpDTevO09lQEqMjwvyY+4Bxd1sVRgGFFfWUlju8HFRIl2Dwo2ItNzOJVC6lwpbLJ+5xircHI67Wyqktoxe3azZZNvy1Hoj0hn8Itw8+eSTZGZmEhERwYQJE1i5cmWLrps/fz6GYXDeeed1bIEiYvn+fwB8YpuEgzAGdFe4aZa7W4rqUu/7tF3hRqRT+DzcvPbaa8yePZu77rqLNWvWMGrUKKZOnUpeXt5hr9u1axc33XQTkyZN6qRKRbq4yiLY/AEAL1SeBKCWm8OJ8ISbEu/7pHAj0jl8Hm4effRRrr76ai6//HKGDRvG3LlziYqK4oUXXmj2GqfTyaWXXso999xDv379OrFakS5s/RvgrKE66Rg2uDKJjQihe2y4r6vyX56WG0epN9xoULFI5/BpuKmpqWH16tVMmTLFe8xmszFlyhSWL1/e7HX33nsvKSkpXHnllUf8Hg6Hg9LS0kYPEWmD7/8LwPZe5wFWq41hGD4syM+5x9xQXaqWG5FO5tNwU1BQgNPpJDU1tdHx1NRUcnJymrxm2bJlPP/88zz77LMt+h5z5swhPj7e+8jIyDjqukW6nLJcyFkPho1lEacCaLzNkXjDTQkDuscCkF1STbmjzodFiXQNPu+Wao2ysjIuu+wynn32WZKTk1t0ze23305JSYn3sWfPng6uUiQIVeRbHyMT2XDAWrRP422OwNMtVVNGfISN5BirC2+HWm9EOpxPlxZNTk7GbreTm5vb6Hhubi5paWmHnL9jxw527drF9OnTvcdcLmtTupCQELZs2UL//v0bXRMeHk54uMYFiByVqiLrY1Sit2tF4eYIPAOKARxlDEiJpqDcwfa8ckZlJPisLJGuwKctN2FhYYwdO5ZFixZ5j7lcLhYtWsTEiRMPOX/IkCGsX7+etWvXeh8///nPmTx5MmvXrlWXk0hHqbTCjRmZ6N0AUuHmCELCISTC+rzBoOLtGlQs0uF8vinM7NmzmTlzJuPGjWP8+PE8/vjjVFRUcPnllwMwY8YMevbsyZw5c4iIiGD48OGNrk9ISAA45LiItCN3y01VSBw1dS7CQmz06hbl46ICQHgc1FW7x91oULFIZ/F5uLnooovIz8/nzjvvJCcnh9GjR/Pxxx97Bxnv3r0bmy2ghgaJBJ+qAwAUYw2M7Zccjd2mmVJHFBEHFXlQXcrA1HRAY25EOoPPww3ArFmzmDVrVpPPLVmy5LDXzps3r/0LEpHG3N1S+XVWa426pFrIM2PKUcqAdOs9+6mo0tv6JSIdQ/93iciRuVtu9jmsPZIUbloovH6V4pTYcGLDQ3C6THYVVvi2LpEgp3AjIkfmbrnZVWXNPFS4aaGI+v2lDMOgvxbzE+kUCjcicmTuAcXbykIBhZsW83ZLlQBopWKRTqJwIyJH5m65ya6JwmZA3+RoHxcUIBp0S4HCjUhnUbgRkSNzt9wcMGPonRhFeIjdxwUFiAb7SwGaDi7SSRRuROTwXC7vgOIDZqy6pFqjwc7gQKPdwZ0u01dViQQ9hRsROTxHKZjWNiclRHsHxUoLHNRyk5EYxZCQbN613cyBla/5sDCR4KZwIyKH5+6SqjbCcRCm3cBbI6LxmBu7zeDy6OUMse0h5ss7wVnrw+JEgpfCjYgcXqV7dWLTWp1Y3VKtcFC3FMBI+y4AIqpyYfMHPihKJPgp3IjI4blbbgpd1gwpdUu1wkHdUpgmmTXb6p9f8Uzn1yTSBSjciMjhVdbPlOoRH0FcRKiPCwogB3VLUbKHyLoSak07ddhh9zeQ/YPv6hMJUgo3InJ4DTbNHJIe5+NiAoynW6quyhpfs38tAFvMDD5zjbeeW/l/vqlNJIgp3IjI4bm7pYrNaIakxfq4mAAT3iAMVpdC9joANtGX52tPt46vf9PbOiYi7UPhRkQOz9MtpZab1rOHQJh7jJKjBLLXAlAQO5RV5mBK44dCXTWsecl3NYoEIYUbETkss9LTchPDULXctF7DLRjc3VLO1JGAwfLuv7Se++55cNb5pDyRYKRwIyKHVV2aD0CZEac9pdrCM6g4fwtUFoBhJy5zNADv1U2EyEQo2QM7FvmuRpEgo3AjIodVU1YIQGRCd0Ls+pHRap7p4Du/tD52H8KAnt0BWJ9XAwOmWMfzNvmgOJHgpJ9UInJ47tlSickpPi4kQHm6pbKWWh97jGZImnVsd1EltZHJ1vHKAh8UJxKcFG5E5LDCaosBSE/r6dtCApWnW6p0r/UxfRSJ0WGkxoUDkOd0DziuKPRBcSLBSeFGRJpXV0OkqxKAjF4KN20SftAMs/TRAN7Wmz0OT7jJ78SiRIKbwo2INKuq1OoqcZkGAzJ6+biaAOUZcwNg2CBtOABD0q2ZZzsqrRYcdUuJtB+FGxFp1k97ra6UUiOa7vFRPq4mQEU0aLlJHgRh1oyzoe6Wm43FYdZz6pYSaTcKNyLSrP37rXBTFRJ/hDOlWQ27pdJHeT8d7F4zaE1RiHVA3VIi7UbhRkSalZefA4AzPMG3hQSyiIT6z93jbQD6d48hxGawu9rdIlZXBTUVnVqaSLBSuBGRZpUU5AFgj070cSUBLKLplpuwEBsDUmKoIAKnzdM1pXE3Iu1B4UZEmmSaJlXu1Ykj4rr7uJoA1qhbamSjp6yNSA0qQ7tZBxRuRNqFwo2INCmvzEFEbQkAsYlawK/NkgZASCT0ngjhjffm8mxEegD3mCbNmBJpFyG+LkBE/NOm7FISKAMgJDrZx9UEsJjucOOG+t3BGxjiHlScWxdNb1DLjUg7UcuNiDRpc04Z3Yxy64uobr4tJtBFJ0NoxCGHPQv57a3RQn4i7UnhRkSatDm7lARPuIlUuOkIqXHhJESFUmC6u6vULSXSLhRuRKRJm3PKSMATbjRbqiMYhsGQtFiKTPegYy3kJ9IuFG5E5BBVNU625ZU36JZSuOkox/SIpwBPuFG3lEh7ULgRkUP8uL8Ep8tVH27UctNhRvaKp0jdUiLtSuFGRA6xbm8J0VQTSp11QC03HWZ0RoK3W8rUbCmRdqFwIyKHWLenuL7Vxh4Oodo0s6P0ToyiJiIJALNc3VIi7UHhRkQO8cPeYu8aN0R2A8PwbUFBzDAMevToBYDNWa39pUTagcKNiDRSUlnLrsJKEgz3L1l1SXW4wRnpOMxQ6wt1TYkcNYUbEWnkh33FAAyIqbEOaDBxhxuZkdBgxpTCjcjRUrgRkUZ+2GvtJzUswTOYWAv4dbRRGQneGVPVpbk+rkYk8PlFuHnyySfJzMwkIiKCCRMmsHLlymbPffbZZ5k0aRLdunWjW7duTJky5bDni0jrrN1TDED/aId1QC03HS41LoJyewIA+/ft9W0xIkHA5+HmtddeY/bs2dx1112sWbOGUaNGMXXqVPLy8po8f8mSJVx88cUsXryY5cuXk5GRwRlnnMG+ffs6uXKR4PTD3mIAekZUWwe09ULniO4OQEGuwo3I0fJ5uHn00Ue5+uqrufzyyxk2bBhz584lKiqKF154ocnzX3nlFf7whz8wevRohgwZwnPPPYfL5WLRokWdXLlI8MktrSa31IHNgGSbVifuTJEJKQCUFub4uBKRwOfTcFNTU8Pq1auZMmWK95jNZmPKlCksX768RfeorKyktraWxET9ABY5WuvcXVIDU2IJcVhjb9Qt1Tnik9MBqC1rutVaRFouxJffvKCgAKfTSWpqaqPjqampbN68uUX3uPXWW+nRo0ejgNSQw+HA4XB4vy4tLW17wSJBzjOYeFRGPBwosg6q5aZTpKZlABDhKKK4soaEqDAfVyQSuHzeLXU0HnjgAebPn88777xDREREk+fMmTOH+Ph47yMjI6OTqxQJHOv2FhNFNVeWPAH7VlsH43r4tqguIqqb9UdeolHmDZki0jY+DTfJycnY7XZycxtPfczNzSUtLe2w1z788MM88MADfPrpp4wcObLZ826//XZKSkq8jz179rRL7SLBxjRNIvYs45OwWxm853Xr4IRrIX20T+vqMqKSAUgySr3dgyLSNj4NN2FhYYwdO7bRYGDP4OCJEyc2e91DDz3Efffdx8cff8y4ceMO+z3Cw8OJi4tr9BCRQxV881+e5V4ybPmY8Rkw412Y9oC2Xugs0db+UkmUsk4tNyJHxadjbgBmz57NzJkzGTduHOPHj+fxxx+noqKCyy+/HIAZM2bQs2dP5syZA8CDDz7InXfeyauvvkpmZiY5OdbMgpiYGGJiYnz2OkQCnX3VswAsCZ/MqX/4L4TH+riiLsY9FTzSqGHrHs2YEjkaPg83F110Efn5+dx5553k5OQwevRoPv74Y+8g4927d2Oz1TcwPf3009TU1PDLX/6y0X3uuusu7r777s4sXSR4VBbR7cB6AFYN+BOnKth0vrAYTHs4htOBq6KA7JIq0uMjfV2VSEDyebgBmDVrFrNmzWryuSVLljT6eteuXR1fkEhXs+MLDEw2uzIYOHCQr6vpmgwDIzoZSveRRAkrs4o4d3RPX1clEpACeraUiLSPmi2fA/ClayTH90vycTVdWLQ1qDjRKOPbnYU+LkYkcCnciHR1pom53Qo3W2MmkBrX9LIK0gkazJhavkPhRqStFG5EurrcDYRX51NphhM18CRfV9O1uVtuko1SdhVWkl1S5eOCRAKTwo1IV+dutVnuGsZxA9J9XEwX554xNSjWWlVdXVMibaNwI9LF1W5tON5GWy34VJQ13mlAtLUj+7c7inxZjUjAUrgR6cocZdj3rAAgK2EiKbEab+NT7m6pnqEVACxXy41ImyjciHRlWV9hM2vZ5Uqlz8Dhvq5G3N1SCWYJdpvB7qJK9hVr3I1IaynciHRl7vE2SzUF3D+4Z0vZq4oY3jMegG81a0qk1RRuRLoq08S5Tevb+BX3/lJU5DPR/d9Dg4pFWk/hRqSrKtqJveQnakw7eUnjSY4J93VF4u6Woq6KE3pb45807kak9RRuRLoqd5fUKtdgxgzo5eNiBICwGLBbIXNcshO7zWDvgSr2FFX6uDCRwKJwI9JVucPNEtcodUn5C8OA2DQAoqpzGdnLPe5GrTciraJwI9IV1VZjZn0FwJeuUUzoq/Vt/Ea3TOvjgV0Nxt1ovRuR1lC4EemKdn+DUVdFjtkNug8jSeNt/EeDcONpUftmRwGmafquJpEAo3Aj0hVtXwTAl85RnDIkxcfFSCMNws34volEhtrJLqlmw75Sn5YlEkgUbkS6INe2zwBrCviZw9N8XI000iDcRITamTzEmkH18Y/ZvqtJJMAo3Ih0NcV7sBVswWkabI0ey+heCb6uSBpqEG4Aph5jhc+PNuSoa0qkhRRuRLqaHVaX1PfmQE4YPgCbzfBxQdJIYl/rY3ku1FRy2pAUwuw2duZXsD2v3Le1dYTvnoOXz4OqA76uRIKIwo1IF+PyrErsVJeUX4rsBhHWFHCKfyI2IpSTBlrbMny0IceHhXWA2ir47G7YuRg2vufraiSIKNyIdCXOWlw7FgOwJmws4zM1BdwvHdQ1daa7a+rjYAs32z6FmjLr870rfVuLBBWFG5GuZO93hNSWU2jGkjHsBELs+hHglw4KN1OGpWK3GWzMLmV3YRCtVrz+zfrP9yjcSPvRTzaRLsTcas2S+so1gqkjevi4GmnWQeEmMTrMu9Bi0Myaqi6FrZ/Uf12wFSq1WKG0D4UbkS6kapP1y2SF7VhOGKAtF/zWQeEG8I6PCppxN5s/BKcDkgZCYn/r2L7Vvq1JgobCjUhXUZ5PVNGPANgHnkZ4iN3HBUmzmgg3ninhJXs2klNY3OkltbsN7i6pEb+EjPHW53tW+K4eCSoKNyKdIfsH2LnEpyWY7l8cW1y9OHHUMJ/WIkfQMNy417ZJjYvg6tQtfBF+E0Xv/sVnpbWLigJwD2xn+AXQ6zjrc427kXaicCPS0SqLYN7Z8PK5sOMLn5WR/aO1UeY6BnHK4O4+q0NaID4DDBvUVVvr3bhdGGl12yTsXRzYC/ptfBdMJ6SNhOSBkDHBOr5vNbicvq1NgoLCjUhHW/F/4HDvC/TBjVDjm9ku5du/AcDe+ziiwkJ8UoO0kD0U4ntZn3u6pkyT/hVrAOjh2s/arVm+qa09bHjL+jjil9bHlKEQFgs15ZC30Xd1SdBQuBHpSNWlsOJp6/OQCOsX1dKHOr2M3fmlZFRtBmDciVM7/ftLGxw87qZoJ/ay/d6nv/ry004vqV2U7IOfrKDNMb+wPtrs0PNY63N1TUk7ULgRac7mD+HVX8M3T0Dx7rbd47vnoLoEkgfDBc9Zx775D+RsaL86W+CjLz4n0qihwoimz+Axnfq9pY26ubdh8ISbg8Zs1e5exZ4iP1zzproE5p4EjwyFhTfDrmVWV1PJXuv/pVcvAkzoPRESMuqv83RN7f3OJ2VLcFHbtEhT6hzw/p+gIh+2fgSf/hVHyii29bmY9d3PpqDMQUG5g4hQO2nxEaTFRZAWH8GwHnH1s5BqKmH5k9bnk2bD0OnWY9P71r2v/NT6i7WDlVbXkv/jV2ADR+oYom36myYgeFpuitzdT1lLrY8xaVCew0hjJ/O+2cUd5/jZ4PCFt0DOeuvzlc/AymeosscS6SzznuIyQnAcfwORDa/TjClpRwo3Ik3Z8BZU5FMdnsRee2/6Vq4jPG8dw/PWMaemhK9dIxqdPtH2Iw+HzuU9RvNd/z8yccRAppa9Q1RlAST0geHusQXTHoIdS2DfKlj1Aoy/usNfyuvf7WGYuRWAboNP7PDvJ+2kYbeUywW7rAHhTPwDfHYno23bufG73dwwZSCxEaGdW1vWUljzMkz4PfQah2mabMwuZdfSVzh783ycpsE9dTMYbuziDPsqEpxluEyD78zBfOicwEfOCRz4n4vRGd9w4oBkTh7UnWN7jsMAKNppzaaKTu7c1yRBReFGpAHTNFnzUxFpHz9KT+Dx8tOZ6/w5yZRwR9grnGtbxpzo13l68FQSYyOprHFSVFzMX7OeJcVVyIUs4rQdK3lwy6+ZGPImUQYUjJlFst39v1pcD5hyFyy8Cb64D0b+qn6TxA5Q53Tx4te7eMXYBoDh+etY/F/DcJP3I1QWQmg0jP0t5uf30J0SYqvzeH3VXq48qW/n1fX9K/D+9eCqw9z4LmtH38PtO4ZzIOcnPgn/OxjwtPNcvkm6gKz4CFbFhjDC2EmBrTtZNfEUlDsIL6qk7kAVq346wKqfDvCvRdvo1z2at6L60q0yy+qaGjyt816TBB2FGxGsEPDu2v08+9VO4nNX8lr4NqrMMFYmTud3w/pxYv9kjkuZAk+Po7djB3P6b4Qxl1oXL7oPduRhxvWk2hZNUvFWHgp9FoBsM5HJH6cwNft7Zk0ewMDUWBh3Bax8Fgq2wIpn4JSbO+x1ffJjLtXFOWRGuKcT9xzbYd9L2pkn3JTnwNaPrc/7nAAR8RipwyBnPSNtO3jx61789oRM7DajY+sxTVjyAHz5AAAVkT2IrtrPmNW3c17dOYwI/4kEo4Ki+GP49RVPMis+psHFh/6721NUydfbC1i2vYDFm/PYmV/BpyF9uCgki0Wfvk+fbicxICW2Y1+TBC11vkuXVlPnYv7K3Zz2yJf8+Y11bM4p46ow6xdJ+dALeWv22dw+bSgnD+pOZEIKnHyTdeEX90FNBRRsg6//BYAx7SEi//gNnH4vhEYBsCjxYqrNUN5du5+pjy/l3vc3Ul5rwim3WPdZ/oQ1o6oDOF0m/7d0B6Nt260DyYMhsluHfC/pAJHdINzdqvf9K9bHvidbH3tYM4uOD8ti74Eq3l+3v4kbtKPqEljwB2+wmR9+ISMOPMS/684D4PchH3CisR5CIkm87CWSGwWbpmUkRvHr8b154pJjWfHXKdx//ghy463u3uj8NZz+2FKue3UNm3M65v8PCW4KN9K82irI+srq7w8ydU4Xr6/aw+SHl3Db2+vZXVRJUnQYfz8lhinGKgC6/+xPGMZBfw2PvwYSekNZtjXzY+FN4KqFgWfAkLOt9UlO/BPMWgWXvsVvrv8HH/zxJE4florLhBe+zuJnjyzhQ9fxmMmDoLrYGnTZAV5YlsUPe0s4PnSHdSDjuA75PtJBDAO69bE+P+AeVNzvFOujuwVuSoIVau79YCMF5Y72r6HqACyeA4+PgHWv4sTO7bVXclvJ+cRFhRN++p1Un/uctcwBwBn3WYvytVJMeAiXTOjN9TN/A8Cx9h2MYzMf/pDNmY9/xbX/W832vPLGF+Vtsmp79zprerlIA4YZ0Mtctl5paSnx8fGUlJQQFxfn63L82/t/gtXzYNyVcM6jvq7m6FWXYK59lS17C7gzaxgrC60fyN1jw/ndyf24ZEJvor64A759CgZMgd+81fR9NrwFb14Bht1aZdUeDtd9C4n9Dvvtv9yaz53vbuCnQmv67u09f+B3hQ9Yf6HfsB7C268JfnteGWf9exk1dS5W9HiM1KLvYPq/YOxv2+17SCd47TLY9J71eUQC3JIFNps1G2nuSZjhsZwd+Qobc8qZekwqcycUYuRvgYnXHd1MPNOEb/4NSx/2LkC5k17cUXMZX7tG8Ktxvbht2lASo8Os8wt3WAOBB0yxQllbuVzwwhmw9ztMw85byb/j5r0nYpoGNgOuPaaO3yWvJ27HB5C/qf666O7wq5etbrtAt2uZ9bMgfdSRzy3cYf2hlXlSx9flB1rz+1tjbqRpFQWw9v9Zn6963lqDYtRFvq2prcrzMb99CueKZwipLWcI8Ipp54uI8dQcexWnTz6BCLMairfBmv9a1xx/bfP3O+YXsPwpa8YTwKQ/HzHYAJwyqDuf3HAyTy3ZwdwlO3hw33DOiEinb1U25spnMSbNPvrXitUq9ec3fqCmzsVpgxJJyXav+NpLg4kDjmfcDUDfSVawAeg+FEIiMRxl/PvcWM78XwUbN/6AK+tW7K4aK3SfdGPbv++OL+CzOwHYF9aPf5Sfw0eu8QxMjeON80dwXGZi4/OT+luPo2WzwWUL4P0/YWx4k1/mP8UZQ3ez+EASQwoXMXj7XnD3spr2MIz+p1nr5+RugJemw5kPwHFXHV3A8qXti+B/7oUNJ1wLP7sTwqKaPremEl48yxqTdcnrMEiLczakcCNNW/MSOB1Wc3NdtdWKkzYcUo+pP6c0G4p2WP3x1aXWeYPPgthU39V9sFUv4ProNmxOByHAVldPyoxYxhqbmcpyWLMc1hx0TfJg6P+z5u9pGDD1fnhxmhVqTvxTi8uJCLUz+/RBTB+Zzk1vrOM/+8/l0bC5lC9+jKqhM0lJTmrTy2zohS9+YN2eYuIiQnhoUijGKxUQHgfdhxz1vaWTNQo3p9R/bg+x/rLf8y0DarZw/c/GMXTJP61gA/DFP6xWlLTGSxa0iMuFuegeDOANYyq3lF6G3Wbn+p8N4LrJAwgL6eDRDOEx1oKXGRPgk78Qt/MDzgWwQR0hfOkcwYfOCSyzj+cX3Y7h92elkvDZbPjxbaubOGc9nPN4fRDsbKYJOT9YH1OPsbqqW6LOYS166LHiadj2CZz3NPQ+/tDzV8+zgg1Y1/U9GUIjDz2vi1K3VFdgmo3+kqmsqWNPURV52bvptvb/yLal8mX8zzlQWUtpVS3Oulr+nTuTZFcBT8XNZlLNUkZUr6IgPIP5o19mUFgBY/a8RPJPCzHMg8bjxPWEqxZBXHonv8jGTNNk5XfLGbtwOiHU8b1rAM+Y59Fz/Pn8fvJAksu3WjOWfngd6qqsTQpDoyEyAc76Z8umoRZsh6hE69EGdU4Xz3y5lbO+PJdMI4ddZjohKQPp2XuANRtmzGUQGtHyG1aXUDb/KmJ3fcpnzmOxT76N02L3woezod9kmLGgTXWKDzX8S/6676D7oPrnPv4LfPskjL+Gun6nETL/19SadrIihjLIsQFSjoFrFkNIeKu+ZdHK+SQu/B1lZiSnOB4jJa0nD184iuE9O27JgmbtWWkFlphUOOZ8zMHT+Hqvk4c/3cLaPcWANV7n0vEZXBfxEXHL/g6mCyb/rUNnITbJ5YLNH8Cyx2C/+y+mkAhIH01N+rHsyTiP7bY+7D1QRU5JFaVVdZRW11JWXUd1rZNfVvw/fl3+Mgdsify/pFlccuBpEuryMTFYN/x2HMdeRUZiFKlxEdid1fCvUdamqrYQcNXBKbfC5ADfLf4IWvP7W+GmPX3xd6vLItU/Vgytrq6m/N2bSdjyOpvjJ/F26Fl8eCCDgrIqLrV/zk0hbxBnWOM/fldzA5+4rG6LabYVPB32LwrMOE5w/Idoqvgg/K/0NArJMbuRZhzwfo8sMw1HaDxExNOzbjex1Tk4U0div+Ij6y+wTuaoc7JwfTb/t2QHdx+4leNtm1jsGsMXY/7DdacNJC3+oLDgrLWWhg8J91lT9v5l/6PH59cd+sTwX1p/wbakrrzNVP7310SVHbSZYmQ3a1DoybfAaX9tn4Kl85TnWb/EEnrDH75t/G9h/Zvw1pWQOsLacPJAFs86pzO39iwWR91OnKsYTrwBTr+nRd/K6TL579fbOPXzn5NpZPN43S9xnnwLfzxtYMe31rSSaZp8sTmPhz/dyqZsa0xQqN3g/j7fc+H+hwDD3VVzRscXU1EAP75jbZBbaK0n5bSFU2uEEdFgVWaA953H83jdBewwezY63svI4/Owm4kwarm+5jrec51IHBXcEfJfLgxZSp1p41c1d7LGHESY3cYNcV/wh6pnKA1PY9fomxi54iZMezjGH5a3T/egnwq4cPPkk0/yz3/+k5ycHEaNGsV//vMfxo9vfnzAG2+8wR133MGuXbsYOHAgDz74IGeddVaLvleHhZvNH8L8S6wUffwf4NTbICy6/e5/BMWVNfy4v5Qf95ewcX8pP+3bz00l/+BE24+NzvvR1QcDGGb7CYBSWwJxrmKq7bG8O/F1QhN7c8rXM0gqXM32Idey9Zg/UV5dR3ju95yz+nLsZh1ObCwOOYn/VJ3FOmdv770zjFzeCbuLZKOUb0PG8+6QhxjdJ4kxvbsxoHsMtg5ah8M0TdbvK+HN1Xt5b91+iitrucC2lEfC5lJrC6do5lek9hncId+7vTgLdrDo62/4evUPpLpyuMb+ASGGi4KT7iV5yuG7vcyN71L71u8Jc1ayz0zi5fhrmd1rC+Gb3rL+igW49E0YeHonvBJpd0U7ISwGYlIOPf7vBvuExabz4Snvc+OC7ZzqWsEzYY9hYmBc/hH0mdjs7eucLj5cn81Ti3cwpuBdHgh9jhIjnrwrVjAww7ctsEficpks3pLH/325k5W7igD4e8jz/CZkETUhsTguX0RsT/f/+yV7rdmfKUOtLr2j+WOmphI2LqDuhzexZy3BMJ0AlBLNvLrTealuKkXE0tfIYbSxnZ/Z13C23doQ1IWNDUlT2ZZxIRUpxxIXGcYJK68jJXsxBd0n8M2JL1JZ4+RAZS3FlQ7O3HIHY0o+J8foztmO+yl3hfJl+I2kGQf4a+0VvOL8GS+FPsgp9h9YGz6W9475DyPSwjm+cgmpecuwZYy3VkJvafdYe8jbbM30a+dusoAKN6+99hozZsxg7ty5TJgwgccff5w33niDLVu2kJKScsj533zzDSeffDJz5szhnHPO4dVXX+XBBx9kzZo1DB8+/Ijfr6PCjVm8Bz6+DWPzB9aB+AyY9qA1PbglnLVW18gRZjiUO+rYVVDBjvxyNueUscX92Fdc5T0n08jm+dCH6W/LppII5idcw2h7FiMPfEqIyz1dNCLBGqw2+lKYdxbsWw0Zx8O0B+CZU62QdsN6a0Vdj+2fw57vYPQl0K0PLpfJngOVbMouY3NOKZuyS3HtXskTNXcSbtTyQt2Z3Fs3A4DY8BBGZSQwolc8I3paj17dIg+dau3hKLPeE88/z7CoRv+jVDjqWJlVxLLtBSzZkseO/ArvcwNja3iPG4isLYYpdx/dwMpOtr+4irve+5GMLfO4M/S/1Jp2/p3xKJN+9nNG9oonIrT+30eNw0HR+38jbYM1lXy5cxifHvMAt/3yJGt/q8Id8PXjUFcDP/8PhIT56FVJhzBNeKiv1TIHcMHzMOKX/LC3mN/9dzWzKx7nwpClVIUnUz79GboPrx9HZjprKf/iMerWvcZn1UN5uWIC28xeLI2YTSpFuKbOwTbxDz56YW2z+qcDPLN0B19u2s//Qu5jnG0rW81eLEr5Lee4FtOr8BsM3D9PYnvA4DOtMYJ9T2nR/xvVtU42ZpeS/8NnjFt3J0m12d7nfnD1ZYHzJF5znkqNPYqh6XEc0yOeYemxDE6LY3BaLPElW2DxP2DLwvqbJvSGPifBuletn7nXfgPdD/pDrLoU/u9kOJCFOeQcilNPoNuXf6E8PJV7Mv/Hj3kOavK28mHILYQbdXzgnMBJtg0kGPU/E3PCM1k17C/EDj2NfsnR9EiIrF/0sWAbbFwAmxdaNQw83XqkjWrbuKWN78E7v8ccchbGL55t1xbxgAo3EyZM4LjjjuOJJ54AwOVykZGRwR//+Eduu+22Q86/6KKLqKio4IMPPvAeO/744xk9ejRz58494vfrqHCzLbeM6U8s4+dRG/hz7bOkuqwVYXNjhrGz9y8p7v9zYuISiAqzExFqJzLUTrijkIidnxK582Mi93yFCVQlDKYkfhCFUf0pdMVQVGMnv9rGrqpIPi/pSX55TTMVmJwRv58LI1cyqfxjIurKcMb2xHbJfIz0kdYplUWw9hWrGfWEP9bv3VKUBXMnQU2ZNaWyIt/qXrvwxTa9F8XfvUbCh9cA8H70BdxeegHltYeeFxMeQt/kaPomR9MvKYKhrq30L1pGet5Soos3NzrXZYSwKXkqH0ady1flPdmUXUqdy/qna8NFSIidqcek88uxvZi06R5s3//XmlHy+6869y+WdrJ6VxHONy5nfMUS8swEznb8g0KjG/26xzAkLRazZB9X5NzHWGMLAM87zyL8zPu49IT+zQdGCT7/u8D6o6PPSfDbD7y/SPLLHNz036+4PecGhtj24DQNXgy9iB/6XU1EyQ5m5D7EcM+0I7fK0ESiaousP8z+uLrVY3X8RX6Zg09XrOXMr39NklnU6LnN9CXT2E+EWb8mUF1oLAcyTqes/zlUR6Vh27eKyJxVxB3YQLGRwPKQ4/igeiQbSsK5xf7/+E3IIgD2m4nMrzuNVbGnktxnOGN6JzA6I4Gh6XGN/gg5xN7V1tpWmz+wuhM9DteFuP97eO50a00te7g12eOsh71709XUuSheeDcpa/7lvWSf2Z2PnOM4376MJMPqHvvUOZZ9ZjLRtlq6hbkYYvxERm1Wk9+yLiqV6oyTqe49CUfGJFwxaVTXOqmqdVJV46TcUeduXarhQGUNhaXVTNj7POeXvAzAjxHHcszsD5uf7dUGARNuampqiIqK4s033+S8887zHp85cybFxcW8++67h1zTu3dvZs+ezQ033OA9dtddd7FgwQLWrVt3yPkOhwOHo/4fcmlpKRkZGe0ebpZtK+A3z1u72Ubg4I8h73CVfSHhRh0A5WYES1yjCaOWJKOUJErpbeRhM1r+9q9yDeK+2t+wN2oYfZOjGZwWy7i4YiYULyR1z0LsxbvqT+45Dn79astnLnn67z2u+BR6T2hxbYf4+l/eqaRmnxPZOulfrCoMY8PeYsyfvuak4vfox36iqCbKcBBLJZFGc8GtsZWuwXzvGsCw8DyG2HNIqt2PYRgYkd2sMSYF1iaRXP7xYZvj/Z6jHMfcyYQf2Mo2evNF3QhyzURqCOHGkDdJMsooI4pXUm/huLNmMrZP2wY2SwD7abk1q2bKPZDYeH+pWqeL57/YwIDV9zKl+jMANrgyGWjsJdyoo9SM4o2IX3JaQg6ZhUsx6qqtC899qn5rkQBm7v4W139/gcOI5POIKTxx4Hi21qYQTg0TbRv5mW0NZ9hXkWoUt+h+1WYoEYb1V9ra1AsonPgXRvTvRUpsKwb9N1RTCVs/gvVvWaH0F88cfijD8qfgk9utz2PT4fq1jScc1FbB21eDsw7GXY6r38/YXexg++69JH/3MCOz38TGoQuy1pp2vnYN5yP3mMvJtrWcZFtPjFHd6LwdrnQ2mn3Y5urFVrMXu8w0is1oSojGAB4OnctZ7q63F+rO5M2k37Pwxslte2+aETDhZv/+/fTs2ZNvvvmGiRPrfwndcsstfPnll6xYseKQa8LCwnjppZe4+OKLvceeeuop7rnnHnJzcw85/+677+aeew5Nw+0dbmrqXGSXVJFX5iCv1EFuaTVlRTkM2P8e4wrfJbW26RU015v9+MI8jiXGcRASzoiQPQyx7aGfuYd4m4NoWy2RRg3dKrMIcbq7nkZeZM1+WftK/U7BYC35P+hMGP4L62NrWywWXAdr/wfpo+GaJUffnPjjAnh3ltUiFJNq7SC84S1rTYomVNhiWBs2lq+McSx1DqfKHkdYiI2wUDsjbVn8svYDRpYuxm7WHfl7HzvD6ooJdAXb4JnJ1nt4kOrk4YRf/F+MpCOvsSNdm2PVK4R8fBP2OmsCQWnGaYSd928ikjKsE6pLYNMH1szBsVf4bhp1e6uptFqgbHYcdU427CthR34Fuwoq2FVYQfaBSvpXb+AEx1ecXPsNUVSxLXQwOyOHkxtzDP1D8hlRsZyUA6uxuWoxu2Vi/Pw/9dtgdCbThPmXwpYP4ZzHrD3qWiNnPWx4GxcGZU47RQ472c441kUcx+6qCHJLqympqqW8ug5HdRUDHRuYYP7ABNZzDDsP+4e4CwMbJk4jhK3H3Ytz1G9Ii48gOaZ9W/8UbhrorJabwzJNa9XJvSutsS7R3a1HYl+ITWvZPUr3Wxs0rnv1oCcMGPAzaxzMoDOPbhBzbRWsedlaH6O9RtwXbLNWWW24mmhIJIz6tTUeKTzWCmVh0Vb/85ECWel+a6G9ygJIHgRJA6yHYbPGHlQdsNaL6Hty8IwxKdxhDVgvy4Gy/dYMml7Hwam3t26quHRtBdvgq0eh36nWbvTqvjzUQctmeFWXWu9f6jDfriXjrLW2nUgb0bn//aoOWNPy8zdD/harhuLd1vYxLvcfmzFp1irRR9PifwQBs0JxcnIydrv9kFCSm5tLWlrTv/TT0tJadX54eDjh4T7uOzYMa3XRvpPafo+4HnD+01Yf66J7rV/yx5xvNR8n9D7y9S0RGgkTftc+9/JIHghXL4KPb4e938Goi+HYy9q+gWNcDzj11qafi+/Z9PFAl9QfTrze11VIoEseaP0MkeY1Fxgi4qDXoTubdzp7KHjGUHamyG7WCsgHr4Jsmta4oapiq3Xej/6g9Gm4CQsLY+zYsSxatMg75sblcrFo0SJmzZrV5DUTJ05k0aJFjcbcfPbZZ41afoJaz2MDbzG2sGj4+b99XYWIiLQnw7Ba39txX7z24vPtF2bPns3MmTMZN24c48eP5/HHH6eiooLLL78cgBkzZtCzZ0/mzJkDwJ/+9CdOOeUUHnnkEc4++2zmz5/PqlWreOaZjtlZWURERAKLz8PNRRddRH5+PnfeeSc5OTmMHj2ajz/+mNRUa5bP7t27sTUY3HbCCSfw6quv8re//Y2//OUvDBw4kAULFrRojRsREREJfj5f56azdcm9pURERAJca35/B8l8PxERERGLwo2IiIgEFYUbERERCSoKNyIiIhJUFG5EREQkqCjciIiISFBRuBEREZGgonAjIiIiQUXhRkRERIKKwo2IiIgEFYUbERERCSo+3zizs3m20iotLfVxJSIiItJSnt/bLdkSs8uFm7KyMgAyMjJ8XImIiIi0VllZGfHx8Yc9p8vtCu5yudi/fz+xsbEYhuHrcnyutLSUjIwM9uzZo13SO5De586h97lz6H3uPHqv65mmSVlZGT169MBmO/yomi7XcmOz2ejVq5evy/A7cXFxXf5/nM6g97lz6H3uHHqfO4/ea8uRWmw8NKBYREREgorCjYiIiAQVhZsuLjw8nLvuuovw8HBflxLU9D53Dr3PnUPvc+fRe902XW5AsYiIiAQ3tdyIiIhIUFG4ERERkaCicCMiIiJBReGmCyoqKuLSSy8lLi6OhIQErrzySsrLy1t0rWmaTJs2DcMwWLBgQccWGuBa+z4XFRXxxz/+kcGDBxMZGUnv3r25/vrrKSkp6cSq/d+TTz5JZmYmERERTJgwgZUrVx72/DfeeIMhQ4YQERHBiBEjWLhwYSdVGtha8z4/++yzTJo0iW7dutGtWzemTJlyxP8uUq+1/6Y95s+fj2EYnHfeeR1bYABSuOmCLr30Un788Uc+++wzPvjgA5YuXco111zTomsff/xxrezcQq19n/fv38/+/ft5+OGH2bBhA/PmzePjjz/myiuv7MSq/dtrr73G7Nmzueuuu1izZg2jRo1i6tSp5OXlNXn+N998w8UXX8yVV17J999/z3nnncd5553Hhg0bOrnywNLa93nJkiVcfPHFLF68mOXLl5ORkcEZZ5zBvn37OrnywNPa99pj165d3HTTTUyaNKmTKg0wpnQpGzduNAHzu+++8x776KOPTMMwzH379h322u+//97s2bOnmZ2dbQLmO++808HVBq6jeZ8bev31182wsDCztra2I8oMOOPHjzevu+4679dOp9Ps0aOHOWfOnCbP/9WvfmWeffbZjY5NmDDB/N3vftehdQa61r7PB6urqzNjY2PNl156qaNKDBptea/r6urME044wXzuuefMmTNnmueee24nVBpY1HLTxSxfvpyEhATGjRvnPTZlyhRsNhsrVqxo9rrKykouueQSnnzySdLS0jqj1IDW1vf5YCUlJcTFxRES0uV2SjlETU0Nq1evZsqUKd5jNpuNKVOmsHz58iavWb58eaPzAaZOndrs+dK29/lglZWV1NbWkpiY2FFlBoW2vtf33nsvKSkpatU9DP3E7GJycnJISUlpdCwkJITExERycnKave7GG2/khBNO4Nxzz+3oEoNCW9/nhgoKCrjvvvta3GUY7AoKCnA6naSmpjY6npqayubNm5u8Jicnp8nzW/rfoCtqy/t8sFtvvZUePXocEiylsba818uWLeP5559n7dq1nVBh4FLLTZC47bbbMAzjsI+W/mA62HvvvccXX3zB448/3r5FB6COfJ8bKi0t5eyzz2bYsGHcfffdR1+4SCd54IEHmD9/Pu+88w4RERG+LieolJWVcdlll/Hss8+SnJzs63L8mlpugsSf//xnfvvb3x72nH79+pGWlnbIQLW6ujqKioqa7W764osv2LFjBwkJCY2OX3DBBUyaNIklS5YcReWBpSPfZ4+ysjLOPPNMYmNjeeeddwgNDT3asoNCcnIydrud3NzcRsdzc3ObfU/T0tJadb607X32ePjhh3nggQf4/PPPGTlyZEeWGRRa+17v2LGDXbt2MX36dO8xl8sFWC3DW7ZsoX///h1bdKDw9aAf6Vyega6rVq3yHvvkk08OO9A1OzvbXL9+faMHYP7rX/8yd+7c2VmlB5S2vM+maZolJSXm8ccfb55yyilmRUVFZ5QaUMaPH2/OmjXL+7XT6TR79ux52AHF55xzTqNjEydO1IDiI2jt+2yapvnggw+acXFx5vLlyzujxKDRmve6qqrqkJ/F5557rnnaaaeZ69evNx0OR2eW7tcUbrqgM8880xwzZoy5YsUKc9myZebAgQPNiy++2Pv83r17zcGDB5srVqxo9h5ottQRtfZ9LikpMSdMmGCOGDHC3L59u5mdne191NXV+epl+JX58+eb4eHh5rx588yNGzea11xzjZmQkGDm5OSYpmmal112mXnbbbd5z//666/NkJAQ8+GHHzY3bdpk3nXXXWZoaKi5fv16X72EgNDa9/mBBx4ww8LCzDfffLPRv9uysjJfvYSA0dr3+mCaLdU0hZsuqLCw0Lz44ovNmJgYMy4uzrz88ssb/RDKysoyAXPx4sXN3kPh5sha+z4vXrzYBJp8ZGVl+eZF+KH//Oc/Zu/evc2wsDBz/Pjx5rfffut97pRTTjFnzpzZ6PzXX3/dHDRokBkWFmYec8wx5ocfftjJFQem1rzPffr0afLf7V133dX5hQeg1v6bbkjhpmnaFVxERESCimZLiYiISFBRuBEREZGgonAjIiIiQUXhRkRERIKKwo2IiIgEFYUbERERCSoKNyIiIhJUFG5EREQkqCjciIiISFBRuBEREZGgonAjIiIiQUXhRkQCUkVFBTNmzCAmJob09HQeeeQRTj31VG644QaWLFmCYRiHPH7729/6umwR6QQKNyISkG6++Wa+/PJL3n33XT799FOWLFnCmjVrADjhhBPIzs72Pr744gsiIiI4+eSTfVy1iHQG7QouIgGnvLycpKQk/ve//3HhhRcCUFRURK9evbjmmmt4/PHHvecWFhYyfvx4zjzzTJ588kkfVSwinUktNyIScHbs2EFNTQ0TJkzwHktMTGTw4MGNzqutreWCCy6gT58+/Otf/+rsMkXER0J8XYCISEe59tpr2bNnDytXriQkRD/uRLoKtdyISMDp378/oaGhrFixwnvswIEDbN261fv1o48+yuuvv867775LUlKSL8oUER/RnzIiEnBiYmK48sorufnmm0lKSiIlJYW//vWv2GzW32uff/45t9xyC08++STJycnk5OQAEBkZSXx8vC9LF5FOoAHFIhKQysvLufbaa3n77beJjY3lz3/+Mx9++CGjR48mISGBe+6555BrZs6cybx58zq/WBHpVAo3IhI0Tj31VEaPHt1otpSIdD0acyMiIiJBReFGREREgoq6pURERCSoqOVGREREgorCjYiIiAQVhRsREREJKgo3IiIiElQUbkRERCSoKNyIiIhIUFG4ERERkaCicCMiIiJBReFGREREgsr/B/k36ZnBVfFYAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "intensity_noisy = intensity + np.sqrt(intensity) * np.random.normal(0, 300, intensity.shape)\n", "plt.plot(qzs, intensity, label='original')\n", "plt.plot(qzs, intensity_noisy, label='added noise')\n", "plt.ylabel('Intensity')\n", "plt.xlabel('qz')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have our experimental data good to go for the fitting we need to provide a list of initial values and bounds for the parameters. The bounds are necessary to make sure that the fitting algorithm doesn't lose its way looking for impossible values. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "initial_params = {'heights': {'value': height+2, 'variation': 10E-1},\n", " 'langles': {'value': langle-0.2, 'variation': 10E-1},\n", " 'rangles': {'value': rangle-0.50, 'variation': 10E-1},\n", " 'y_start': {'value': y1+1, 'variation': 10E-1},\n", " 'bot_cd': {'value': bot_cd-1, 'variation': 10E-2},\n", " 'dwx': {'value': dwx, 'variation': 10E-5},\n", " 'dwz': {'value': dwz, 'variation': 10E-5},\n", " 'i0': {'value': i0, 'variation': 10E-5},\n", " 'bkg_cste': {'value': bkg, 'variation': 10E-5}\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create an instance of the Simulation class and pass it to the Fitter class along with data to fit" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from cdsaxs.fitter import Fitter\n", "Simulation2 = StackedTrapezoidSimulation(use_gpu=False, qys=qxs, qzs=qzs, initial_guess=initial_params)\n", "\n", "Fitter1 = Fitter(Simulation=Simulation2, exp_data=intensity)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we'll call the `cmaes` method of the Fitter class to fit the data. It's important to know that this method provides a normalized set of values between -sigma and sigma and it is rescaled to the parameter values according to the model used.\n", "\n", "*sigma* - the standard deviation of the values generated by the `cmaes` method.\n", "\n", " Note: Similar to bound denoted by variation the parameter sigma helps to control the exploration of the algorithm. The larger the value of sigma the further the algorithm explores the parameter space and vice versa. Thus they need to be adjusted for the algorithm to converge to the optimal solution. \n", "\n", "*ngen* - the number of generations to run the algorithm\n", "\n", "*popsize* - the population size\n", "\n", "*mu* - the number of parents to select for the next generation\n", "\n", "*n_default* - the number parameters to fit\n", "\n", "*restarts* - the number of times to restart the algorithm\n", "\n", "*tolhistfun* - the tolerance for the history of the function value\n", "\n", "*ftarget* - the target fitness value\n", "\n", "*restart_from_best* - whether to restart from the best solution found so far\n", "\n", "*verbose* - whether to print the progress of the algorithm\n", "\n", "*dir_save* - the directory to save the results of the fitting" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration terminated due to ngen criterion after 100 gens\n", "Doubled popsize\n", "Iteration terminated due to ngen criterion after 100 gens\n", "[]\n", " height1 langle1 langle2 rangle1 rangle2 y_start bot_cd \\\n", "0 19.999831 1.397397 1.39493 1.395469 1.397299 0.15517 39.996793 \n", "\n", " dwx dwz i0 bkg_cste \n", "0 0.099187 0.100135 9.999407 0.099898 \n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/nobackup/nd276333/Workspace/cdsaxs/src/cdsaxs/fitter.py:290: RuntimeWarning: invalid value encountered in scalar divide\n", " if(n_infs / fitness_arr.shape[0] > 0.5):\n" ] } ], "source": [ "bestfit, fitness = Fitter1.cmaes(sigma=1, ngen=100, popsize=500, mu=10, n_default=11, restarts=1, tolhistfun=10E-5, ftarget=0, restart_from_best=True, verbose=False)\n", "print(bestfit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have the bestfit parameters in the form of a pandas dataframe. We will look at the statistical information generated by the `mcmc_bestfit_stats` method." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "11 parameters\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ " 0%| | 0/100 [00:00