{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Batchfitting of multiple datasets\n", "\n", "This notebook is a demonstration of how to batch fit multiple datasets using `refnx`. Batch fitting is essential for bulk analysis of hundreds or thousands of datasets. Such situations are becoming more common as faster instruments (such as those being built at the ESS) come online.\n", "This example is based on a deuterated polymer film that is gradually being swollen by an hydrogenous solvent. 314 datasets were acquired over a period of a couple of hours." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# some initial imports\n", "%matplotlib inline\n", "import time\n", "from multiprocessing import Pool\n", "from copy import deepcopy\n", "import glob\n", "from tqdm import tqdm\n", "\n", "import matplotlib.pyplot as plt\n", "from natsort import natsorted\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.1.18.dev0+c02b356\n" ] } ], "source": [ "# the refnx imports\n", "import refnx\n", "from refnx.dataset import Data1D\n", "from refnx.analysis import Objective, CurveFitter\n", "from refnx.reflect import SLD, Slab, ReflectModel, Structure\n", "from refnx.dataset import ReflectDataset\n", "\n", "# we print out the version so that others can repeat our analysis\n", "print(refnx.version.version)\n", "\n", "# natsort is used for alphanumeric sort, but it's not essential\n", "# An alphanumeric sort typically ensures that the datasets are loaded in\n", "# the order in which they were run.\n", "files = natsorted(list(glob.iglob('./*.dat')))\n", "\n", "\"\"\"\n", "initial model setup\n", "\"\"\"\n", "# set up the SLD objects for each layer\n", "air = SLD(0.0 + 0.0j, name='air_sld')\n", "polymer = SLD(6.14127029941648 + 0.0j, name='d-polymer')\n", "sio2 = SLD(3.47 + 0.0j, name='native SiO2')\n", "si = SLD(2.07 + 0.0j, name='Si')\n", "\n", "# set up Slabs corresponding to each layer\n", "polymer_l = Slab(520.5277261491321, polymer, 8.762992087388763, name='polymer layer')\n", "sio2_l = Slab(15.305814908968332, sio2, 5.020239736927396, name='sio2 layer')\n", "si_l = Slab(0.0, si, 3.0, name='Si')\n", "\n", "# SLD limits for polymer film\n", "polymer.real.setp(vary=True, bounds=(0, 7.0))\n", "\n", "# limits for thickness and roughness of the layers\n", "polymer_l.thick.setp(vary=True, bounds=(200, 1020.0))\n", "polymer_l.rough.setp(vary=True, bounds=(2.0, 20.0))\n", "sio2_l.thick.setp(vary=True, bounds=(5.0, 10.0))\n", "sio2_l.rough.setp(vary=True, bounds=(2.0, 10.0))\n", "\n", "# set up the Structure object from the Slabs\n", "structure = air | polymer_l | sio2_l | si_l\n", "\n", "# make the reflectometry model\n", "model = ReflectModel(structure, scale=1.0, bkg=1e-10, dq=8.6, dq_type='constant')\n", "model.scale.setp(vary=True, bounds=(0.5, 1.5))\n", "model.bkg.setp(vary=True, bounds=(1e-10, 1e-5))\n", "\n", "\"\"\"\n", "Create lists of all the objectives to fit\n", "\"\"\"\n", "filenames = []\n", "models = []\n", "objectives = []\n", "\n", "for idx, file in enumerate(files):\n", " data = Data1D(data=file)\n", "\n", " # make the objective for each dataset. Deepcopy is used so that all the\n", " # models are independent of each other\n", " objective = Objective(deepcopy(model), data)\n", "\n", " filenames.append(data.filename)\n", " models.append(objective.model)\n", " objectives.append(objective)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since all the objectives and curvefits are independent we can fit the datasets in parallel. To do this in Python one can use a [`multiprocessing.Pool`](https://docs.python.org/3/library/multiprocessing.html#multiprocessing.pool.Pool) object to map over all the datasets. Use of a `Pool` object requires us to create a `fit_an_objective` function in a separate Python file to fit each objective (this is because the function needs to be pickleable). The `tqdm` package can be used to display a progress bar." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 314/314 [02:31<00:00, 2.08it/s]\n" ] } ], "source": [ "from parallel_curvefitter import fit_an_objective\n", "\"\"\"\n", "# parallel_curvefitter.py\n", "\n", "def fit_an_objective(objective):\n", " # make the curvefitter and do the fit\n", " fitter = CurveFitter(objective)\n", " fitter.fit('differential_evolution', verbose=False, tol=0.05)\n", " return objective\n", "\"\"\"\n", "\n", "with Pool() as p:\n", " obj_out = list(p.map(fit_an_objective, tqdm(objectives)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We fitted all 314 datasets in 2\"31\"\" (using a 4 CPU processor), which is 0.48 seconds per fit. If we had fitted in a serial manner this would have taken 4 times longer." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAEGCAYAAAANNmA4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOydd3hUZfbHv4cEkhBSICQIREINqJTQYUEUUVgVZe1dUcHFVWy/ta264qqra2+rFAuLXVlFVxBwXREQkGbovSe0FEJCaCnn98eZl/e9d2aSmZBMCu/neea5d26bNzB3vvec9xRiZlgsFovFUtupV90DsFgsFoulMrCCZrFYLJY6gRU0i8VisdQJrKBZLBaLpU5gBc1isVgsdYLw6h5AZVOvXj2Oioqq7mFYLBZLreLw4cPMzLXayKlzghYVFYXCwsLqHobFYrHUKojoSHWP4WSp1WpssVgsFovCCprFYrFY6gRW0CwWi8VSJ7CCZrFYLJY6gRU0i8VisdQJrKBZLBaLpU5gBc1isVgsdQIraDWBkhLgvfeA48ereyQWi6WW8tRTwM8/V/coqhcraDWBDz4ARo0C3nijukdisViqmaIi4NCh4M7Zvh0YNw6YN68qRlR7sIJWE9iyRZa2wonFcsrz1FNAz57BnTNlCkAE3Hxz1YyptmAFrSZw4IAs4+KqdxwWi6XaWbAA2LgRKCgI/JwPPwTOOw9o1arqxlUbsIIWaj75BMjMdG7Ly5Pl0aOhH4/FYqlySkqAb78FmMs/ds0aWW7dGti1MzKAzZuBSy+t+PjqClbQQsmBA8ANNwDJyc7Hr/37ZZmfXz3jslgsVcr//geMGAEsWVL2cdnZ+ucgUEFbuFCW/ftXfHx1BStooSQnR6+/9ppe37FDllbQLJY6SVaWLHNzyz5u3Tq9rqbWy2PhQiAyEujWrWJjq0tYQQsl5rf5o4/E/1BaCuzaJdusoFksdRJ1a5uOmZ9/Bs4/35mts3atLMPDA7fQFiwAevUCGjSonLHWZqyghRIlaLffLrO+S5YAe/dKnC4AHDxYfWOzWCyVzvr1Yj39+qu8N8PxR4wAfvxRO2gAYMMGoGFDsbYCsdDWrAEWLxZhtFhBCy1K0EaPBsLCZJZ482a931poFkud4p13gGPHJKwecFpo6vlVxYQBwJ49QPPmQPv2vgVtxgyZj1M88wwQHQ3cfXflj702YgUtlChBa9sW6NNHHs9WrZJt3bqJoBUVAZ995gyH+uYb4PTT5c6wWCy1BiVK4eGyVIJWXKyPMafW9+0DmjUDUlJkJqK0VH4Kpk8H5s4FLr4YGDJEji0qAqZNA265BUhIqPq/pTZgBS2UKEFr3Fh8BIsXA/PnA/HxwBlniKDNmgVcdx2wfLk+b8UKic1V+WoWSx1g8WIpklOXUYKm5sl+/FHEx7SyfAna6afLOfv3A3/+MzB8ODB4sPPaq1dLps/AgVX7N5gQUTwRTSWi9US0joj6u/bfQEQriWgVES0gopCGqlhBCyW5uUBsrDyuDRkij1+ffQZ06SJJ1QcPatEzxUut20oiljrEhAnAQw9V9yiqjqNHnTMKgARw5OYCn36qt2Vnyzzaxx9rQVMJ0v/+N/DKK7JeWqrP+eQT4NlnZb1376r7G3zwOoCZzNwJQDcA61z7twE4h5m7AHgawMRQDs4KWijJzQWaNJH13/0OaN1a1pWg5ef7DodSTvbDh/W24mIrcJZazaFDofsKHz8u1lFVU1gIDB0qDpZly5yuRUDPGixerLfl5Ei5qxtvlHVloQHAww8DSUnAf//rvM4NN4jYATKDEQqIKA7AIADvAQAzH2fmPPMYZl7AzOppfBGA5NCMTrCCFkpycrSg1a8P/OlPst6ypVhux47phBVT0JSFZgratdcCjRpV/Zgtliri0CHgyBGn5VFVfPKJePm3b3duZ/ZfvWPVKuA//5Hb7qefgKVLvY85ftw5/gULgB9+kLmtH34A6tUDWrTwPm/9eiAxUX4OsrJkjkxhWmiFhcAFFwCDBgFRUd7XadJEajhWEuFEtNR43eHa3wZAFoAPiOg3InqXiKLLuN7tAL6vtNEFgBW0qqakRN8xpoUGAPfdB7z6qghbbKxsU2WxfFlo5uOsejyzVpqllqJC2M3ntKpCzWXt2ePc/qc/yfOkm82bRUQuvRS45hrJtLnvPucxzEBEBDBmjN62aJEsV60CZs8Wd2D79t7XLy2Vz01IkONUdRBABK1JEy1gvXrJ8+8FFziv0b8/MHNm+X97EBQzcy/j5XYXhgPoAeAdZu4OoBDAI74uRESDIYL2cKWOsBysoFU1nToB118v625Bq19f7pL4eF2YWCVZl2ehKXburPwxWywhQAlasM9kzMDLLzuravhCpXfm5elcL+UAAeR2HD9eRK6kRLa9/z5w//3AAw+I5TN0qEQXbtsmbkR1TUDfopMm6W2qDNXcueJWvOAC5y1v0qKFCNrWrfJZKhKyWTN5r9yOvXrJ8osv5KUYOzbk82cZADKY2ZNVh6kQgXNARF0BvAtgBDPnuPdXJSEXNCLa7omASSeipcb2sZ7ImTVE9IKx/VEi2kxEG4hoWKjHe1IUFcmj3mefAb/84i1oJspCy8iQpS9BM+989fhmBc1SS6mohfbLLxL599xzzu3Tp+v63unpkp81fry49ubMke3Z2fp4M8JS3Vq33y5V6X76Cbj6auCSS/S09pEjunAw4BTHw4dFaBctEjdjbq6I5KWX+g+pVxYaAJx1FpCaKuvNmsmyVSu5Vvfu8j4iQnLUFKed5vefqEpg5r0AdhFRR8+mIQDWmscQUSsAXwG4iZk3hnaEYkJWB4OZ+cRXy2OejgDQjZmPEVGSZ/uZAK4FcBaAFgD+S0SpzFxSHYMOGvPu+fxzEabyBM200EpLRRB9BYU0bix3mFlmwGKpRSgRCdZCe+stWU6fLkEX4eFiPQ0fLu7Bzz4TASsqAiZPlmPUbWWK0KZNev3QIX0LqveDBmkrSbF4MZCWJuumm3DhQpkCP3AAuPJKYOpU2d67N/Dll77/jpYttQD37y8iuHatFrTzzhNRjjZmqUxxNMUthIwF8DERNQCwFcCtRDQGAJh5PIC/AkgA8DbJ5F4xM/cK1eCqS9Dc3AngeWY+BgDMrL4qIwB85tm+jYg2A+gDYGFIR7d1q5SpuuaawM+ZPt1Z+WPtWnlkS0z0fXxSkiyVZVZQIDPZN92kjzEFTQWEWAvNUkvx5XL85BO5TcyvvUlpKfD110CbNuIGnDFDrBsVNfj552KV/fabvF+2zHm++Yy5e7deV8+PJmefrUUuKkpKUr30klhMvXs7xXHhQpnWbtsWePFFEbTx42WfeoatV8/5GS1a6HH37y//HsuX61v70Ue9/35T0EJtoQEAM6cDcAvUeGP/KACjQjoog+qYQ2MAs4lomRFFkwrgbCL6lYh+JiLlGW4JYJdxboZnW2gZOlSiCoPpiz58uJ47A3SitK8ZaMD721lQoAvAKcw7X4mbtdAstRBmfTvt2SMlTQEJR/fXdfngQQkUPn5cqmMkJko9xPbtgcce08d9+qm+3dxh85s3SwQi4GxLeOiQHgMAtGsn1llcnGTXnHkmMHGiWFF//ascY1poa9aIm/NPf5LjjxwB/vhH2acETd3iKiqxZUs9J9erl5Sv2rKl7KhFda3ISNsP2BfVYaENZOZMj1vxByJa7xlHEwD9APQG8AURBZxd4RHGOwCgQVWUnFYuv3XrKjYL26WLLnHlK4YXkG9qeLi+A/PznY+AgNNCU5actdAstZBjx3QgxmWXybK81iq9egEDBsh6q1ZSKeOLL4C335Zbc8gQeb77/HP/ASPTpslr9255tWsnInLokFh8APDee+JuVPzjH2KhXXKJlF9VkYVK0E47TWqNA/p5NTJSn69EqGVL/ZmbN8v7CRPEsuvcObDw+/BwEbJKDtevM4TcQmPmTM9yP4CvIS7EDABfsbAYQCmApgAyAZhe7GTPNvc1J6pQ0/DwKtDoZE9uoDkjXBbuxJauXfW6PwutXj3tPAdEsFaudB6jBI3ZCpqlVuNr3myin5oSt9wi2S2bN+uw+CZNxEt/993AXXfJtt69gT/8QdqylJSIixDQ3vx27fQ1t2+XqhwdPeENBQU6R+13v3OG2qvgEEDcjfv2iVWZlSXuwZYt9XycrwAQVT/hjDNkOWoUMG6ciFibNhLgEow4JSRUj7uxNhBSQSOiaCKKUesAhgJYDWAagMGe7akAGgDIBvAtgGuJKIKI2gDoAGCxr2tXKU2bytKXoE2erGepFWqmF5Cq+uqbDJQ9k2t+S9evd5bhBvSvgJmNum9fmUO3WGoKr70mP9ylpb69948YGU3KFccs81GTJ8t7lU/WuLE+9vrrgX79RMyuukq2DR0K3HqrrI8dK1PYffvqc5Yvl2uryMJDh7SgpaT4/xt69NDn79+vk6PV86UvQevVS+e1AZLJ8+ST8tNQEQYPloARizehdjk2A/C1J/olHMAnzDzTEzHzPhGtBnAcwC3MzADWENEXkNDQYgB3VWmE46ZNwF/+Ir0eGjQQH0hiohYSX4Km7hqzf4N5tyYkaKssIUFib/1hip1bzABtoalgk5QU8bEUFjpDoSyWGsj998uysLD86ejcXHFY5OXJ117desojbwYLN26s878AiWhs0QJ44w15n5Iiz5RmpY0lS2SpBC0vTyp7NG/uuyKHQnWFXr5cLLSkJOdY/IXot2unK42cbGX8d989ufPrMiEVNGbeCilo6d5+HMCNfs55FsCzVTw0YcYMeRz8859lhnfMGKlloxoXrV7t/1xm7Tcw/SlJSdry8uduVKjjzHCoc84RH0rDhlrQ1ONg27YiaFlZVtAstYa8PB1Gr7j+enH1rVol80pK0NRxJa7HWNNCc6NmCJTbULkazaAPFV2oXI7PPCNBIhMmlD322Fjp/PTmm/K+b1+nQJUlVgMGSF6aOQNhqVxspRATdffs2KGDOEaO1NbSrl3+e5Ll5srrwgvFQa5ITNRC5S8gRKGOU8uUFCkmN22a+CmUUCpBU3eqO3jEYqnBXHWV3CYmqakyF3b55fJeBYioOgNu/KVzmlx4IfDddxISD+h4LiIdNKIstMxMEcA73NULfTB5styKWVlATIweS4MGZT9XJidLa0Mz381SuVhBMzEFTVlIGzfKDLD61rofLRUrV4pze+ZM4F//0ttNCy1QQVOWXNeucseMGCF3ii8LDbCCZqkR7N0rYlVetKI7GwXQlo26zVSPMF+3W0RE2W5BRViYNMRUjpPHHpPbuU8feR8T43QxludAUZxxhiRu9+ol2TlqzAkJNvKwurGCZqIiBnfs8L4rlZ/AX97X738vgRxuevUSKy0+XqysslCCpnwmpm+iYUN5LHzrLeDcc2WbFTRLDeKLL8RjP2eOPAOOGRN4yxZTFICyLbRArDNfhIcDHTpo4brgAhE9lcgcTOWN3r1lHu76673Hbqk+rKCZmBaaW9C6dJHl1KnSiQ9wpv0fPy6ZlFdfrbctWybzcWFhEmZ1771lf/6gQVKNRAmWKWjKQlNZoYAWNLP8gcVSiXzzjbNEVFko8dq4UaajJ0yQli3/+x/w4Ye+z1ExUkoUArHQypo/CwQ1g6DcnkrQynOg+MMKWs3BCpqiqEj3lti5UwTNDLc/6yzxJ4wfL534AO9kmvPPd95tpkO9eXNxspdFYqIUouvbV87tb3Q3b9hQglNUIg4gj5r161sLzVIlHDokdQmff778Y4uLdQHgTZucX8khQ/xX/1Cpl0oUYmPl+e+zzyQvbeZM7dhQt1NFLTSFCgQZOlSWKm20orURraDVHKygKXbvFosrOlpbaGeeqfcnJup8NIU79njQIHEtKiragLNvX5knMyujNmwoj6umiMbEyLisoFmqgPnzRahUFQxm4PHHxdngJj1dsknCwuT4rCyZmxo+vOzPUKVNlSgQSZWNpUvF4bFvH9Czp3z9O3eWY07WQnvpJXGeqCaaakq6ooKmhMwKWvVTU4oTVz/Kt9G/v/Q7P3hQygPExcl6XJxTOI4fd7Z4AUTM/FloweKeXfZ1rehoK2iWKkNZXErQMjKAZ58VZ8Y//uE8Vrkl+/aV49u2la/mvfdKSSl/RXZUxJ9pdalntn/+UwJ9+/SRCvdt2wJXXHHyFlrDhjpBGtDPpdblWPuxFppCCZqavwLkm6oe29yVQPPztaD985969rqyBM2NErhBg4DRo2W9Xr2yBW3tWlu82BIUjz6qc7F++kmW+/fLM93mzfJ+/Xr5+nfsKFGNBw7oKhsXXCDHb94sX83zz/efvtm6tdwiRE7HhppXu/NOiVJMTJSp6MsuEzFSpawqC5WJczIux+7dnZVILNWDFTTF8OHi57joIr3NLWjTp+v3pqB16qRDp5SgRUTI/FZlocppjRsnvzgqIKUsQbvxRqkn1KqVFkFLnWXWLPFClxc2P2eOFOBdsQKYN09vLy2V+bIxY6S62vLlehp50yYtaOvWiRW2caPESD3yiAhaYqJuRrl0qdNDbxb7BYB77hHLLTpabpl6xi/Rxo3yfOl2UhDJ3/jAA4H+iwRHRQVN9WP7wx8qdzyW4LGCpoiJEWe9OW/mFrSLLpLS3oAImvJVxMToc9SjZmVX7njySWkENXiw3Nnqbj/tNEkAOn7cu/VvVpY8Pu/aZevlnALcead8JU333r/+Jc9bZhuVRx+VNifXXy9dkVRQhKo2D8gPdHGxDubYuFEL2tatzkK+774rQtO6tU5UPn7c2frvP/+RaEdVO1xNL599tneSdatWOnPFzcCBlV+YV01V24Tn2o8VNDdmrUVfLkf1rV++XE8ymIKmLLSKBoT4IynJ9yPgaafJpMNNN3n/MuTny6O25ZRACZJpoX3/PbBhg96XnS2JzcePi0d6924pXTpvngR2KFQw7fXXy7PTuHHARx/JtpISYO5cWX/jDRHEHTtE0Nq21daWKWixsfIspp7z1C1z1136utXFokUSAGOToms/VtDKokkTuaOfeEILnRK0sWOBl1+WdVO8lKCFqraiEtxZs6RC6/Hj8r60VFyibkHr3l1a6lpqDcXF0kLF3U0IEKupfXtnNofZgEGJlCr1NGuWCJBZ6X3kSHEJmoI2f74IVKtWwOuvyzPT7t366/a//4kAdO2qAyxat5bbRFWrdwcFA/pWMZ8Bq5sWLXSfNUvtxgqaL9SjZXy83K1/+5vepwTNdO+FwkLzh/K/HDwo4WeqWklhobNvGiAil54OPPRQaMZmqRTS08WKuu02732LF0tLldde09tUEd7CQh2hqARt6lTJ/brtNql2YcY6mYL2ww86yGHsWOCFF2S9dWsRsjVr5DapX1/Pm6lrKbejaaEp1HNeqG4Py6mFFTRfPP64LH0563052s27U93VobbQFCtWyFJ1CFAlFwBn3pxqOGWpdvbu9d0tSKFciL6+ert3y1JV6WjQQFtoq1bp+bF168Qt+O23ImbvvCONJH7/e32tZct07n9hoS6OA8hc2x//CLzyiq6Jrb56114rS9VaRQmaLwvN7XK0WCoTK2i+GDtWJgrcofqA718V038TFibHhNpCU6xYIVVY//xneW9OqGQazb6XL6/6sVnKhVmEYeBA/8fs3y9LX/lX6r80O1ssr7ZttaApi6t9exG099+X92PGyNc0PFxKUk2aJNv37HFGI5qFcsLCpEhOv366IpsStCFDRFhVEnUgFpoVNEtVYAXNF0TOOGKTqKjyW80mJZ18OYNAadJEh46Fh4ugTZsGfPmlbDMbSZlhbP/7nyyZgaee8t+nw3LSPP20RAP6QrkC3YnHZplQXymOCmWhASJmKugVkApu4eFS4mndOrHIunXTFTIAcRmqvmGAf0EzcQuae33YMAkAMS08hRU0S1ViBS1YiHxbbiaffioiEQrq1dNWWs+e8iu2f7/zF1GhYq0BaerELL9048Y5iypbKpUff5R4HVUq9KWXgL//Xda/+06WZpWzggJ5ZlJNJJWg+YrCcwtas2baQtu/X6ykXr3kmnPmaLegiSlwStDCw51CZ+JL0Ew6dJDnJV+3iZ1Ds1QlVtAqQnkJK716AW3ahGYsgAgakQhaZqZz3sxECdrIkRIt8NNPWvhMd6SlUlG1D3/5RVyDDz4oXmFAiu8CTotFBXioih3qv8adZmjuA7SFtmOHJDbv2yfOAmUdFhf7FrSWLbVYdu8u1Tg6dPBfFyAtTZYqmjEYrIVmqUqsoFUEJWgLF5ZfliEUNG8ur1atdHSjL5Sg3X67PD7/6186UMRdaNlSKWRl6UIuv/wi81CKo0e1IOXny5JZcrsAbQEpC+3wYfkvZJbnkC+/lPcqkKNtW/FAFxVJBOPKlWKxpabq+TdfghYRIULYvLl8tVu39n2cok0bsfZuuim4fwvACpqlarHFiSuCErT27UM3V1YWDz8s/ix3Oxs3StASE6XX/dSpUowP8C60bKkUlHUWGQl89ZWzdd2WLToYVf3z5+ToY5TYKUHbskW+cp98IgKkvMS9esn8WNu2Tlfezp3iQiQSK+2775wt9kxSU3Wq5TfflC8455xT/t/uC+tytFQl1kKrCLGxMslwsmW/K4sBA6RxVXnF6FRQSEyMxFoXFOjgERvGf1Iw+zZylaC98YZUIYuIkNB5QLy+yjLLz5drbNki79u0ESE7flzPiW3ZIjE+CxboyEdAkq7fflu+Bldc4azPqAr5jh4NjBrlvyL8xx/LtCogoqn6lFU2554rRYbLaw1osVQEK2gVoXFjeUT2FwlZXZRX5E49+sfE6K4CCxdW6ZDqEqoIiy+eflr+WQ8cEFfizp2yfdUq2T5qlBT4/e03bd2sWSOFXGJjRcwKC7WgDRokzxs//KA9yMoAX7FCi1z9+hI2f+edEkiiplIVStAuvVSH5/uiZcuKF+cNhuHDxVK1WKqCGvaLXEt47DGZf6pplCVoZshZdLQ8IkdHO8P1T+G6jytW+I+lAaQgb0yM1DCcMkU3PwDEanrySVlfuRK4/34JmDh4UKyl/v1FaJo1k+2xsbK+dKmcoyIc8/O1oJ19tizfe0+ESnVXVp+xb588Tx05ohOdFVFR+qtQ2a1WLJaajBW0inDGGcB551X3KLxp2lR+/XzlyalfvUaNtGUZH+8UMRVXXodZsUJXB1Mwi8HqblppsmyZWGh//rO4+FQVegCYPVuvr14N/PyzrH/0kbw3W+wpOnTQue2qsrwStBYtZD8gDRZ+9zun9ZSbK+eq/25fqCBbK2iWyoSI4oloKhGtJ6J1RNTftb8TES0komNE9OdQj88KWl2iXj159E9KcnZMBHRSka9WN4o6LmhHjkjIudmtGBBXXl6euAnvvluiEd2ojsxLlsjyyy91seC5c2VKNSpK3Igqr+vuu2XpK4DitNN0T1klaM88I4Z/mza6vR4g7sKGDZ3n//RT2fNcrVvL0gqapZJ5HcBMZu4EoBuAda79uQDuAfBSqAcGWEGrezRvLr9izZo5Q8mUhRYVpbedYoI2ZYos3Z5VFVa/ebM0H1dBGyZK0ADtzpsxQ6Ylf/sNOOssmbtavdoZsBETI1GIbkyhUYL28ceyTElxCtrdd3sL2uHDgQlaVQV3WE49iCgOwCAA7wEAMx9nZkcVUmbez8xLAFRLlJkN269r3HKLZNA2aCAmxwsviPmgzAYzO9edcmDGlNdBPv1UlurHXqEEbfVqWZqphczi9lM1nwHg4ovFUnv0UaljXVIiueoREcAXX8g//bBhkhkxYIDviD5T0MwqIRdeKF2jIyPlc3v2lHUV7h4RIc8s27eXLVaXXCIRli1alPEPYrE4CSeipcb7icw80XjfBkAWgA+IqBuAZQDuZeZy8oVChxW0usbYsXp90yYRtKgo/chv5pspCy05WYJDyir5XgdQZaIOHHBuV4J27Jj3/o8+0vNljRvLvl69JLpw5UpdKrNHD/H4quoeffoAd9zhfyy+LDRAohWVwJn9XJWFFhcnSc/lCVr//lLS02IJgmJm9uFPOEE4gB4AxjLzr0T0OoBHADwRktEFgHU51mWUezEqSj+qm8nXStCaNpVHf/cvvZvly4GLLtK//LUMFep+8KCzZrMSNIWy0I4edbaOu/JKsc4uuUSqZNSvr4MyunVzzs2VZxn5EzR/FdOUoMXG6ioe1p1oCTEZADKY+VfP+6kQgasxWEGryyhBa9jQ9y+sErSYGDE/yrPQZs0Cvv9eJptqCe+9BzzwgIhTfr5uaaIqdADegqZ0fc0aqVx///3y/oorpNpGy5YSeXjsmPxT3HOPs60KUH5Ol9laxZwvc7tDFW4LDbCCZgktzLwXwC4i6ujZNATA2mockhfW5ViXMS00X+FuZnft+HjfFtrixZIEtWaNrsW0a5dEQdQCRo2SperR1amTCFhuri704s9C27BBlrffLonT7p6tRCJAr78u7825skAttPBwZ4qgv5JQpqANGiT90wYMKPszLJYqYCyAj4moAYCtAG4lojEAwMzjieg0AEsBxAIoJaL7AJzJzPmhGJy10OoykZGyVD3cHn5YrCyFstCUoOXlyUSTWdx47lwxZ5Yt05NQ1dw7bdEiSQU0rSx/KItHTS126iTLefNkjmryZP8W2oYNIlrt2wffgLw8C00JWlyc77YwbkyXY9OmMn5/7V0slqqCmdOZuRczd2XmPzDzAWYez8zjPfv3MnMyM8cyc7xnPSRiBlhBq9vUqyeipiy15593lpwwBa1xY4k/P/10SXJSqCzkjRu9q+VWEz/9JMNyJ0i7YdbzZsXFsuzocZbcdpsU4Z04UQRNaT8gcTNFRSJorVvror2B8PbbIj7lVSGLjxfrzOxEVNY5poVmsVh8YwWtrhMV5Z3EpDDn0OLjJWy/tNSZdKX8bps2OV2O1cjWrbJUw9m/35k7Vlws1tuBA5JzduaZep+y0BTbt4ugubfn5cmfrgQwUO68U+Juwstx5terJ5aWEqjdu8sWaGUhWkGzWPxjBa2uExXlTKY2cVtoCtMHpwRt3TqJkACqzULLzASuvx5IT9fvAQnGGDFCajEeOybzZWedpYdpGqWmcHXtKrnkW7bI9ogIHTKfnS1Gad6U+H4AACAASURBVLCCFgxJSdpCa968bLEyXY4Wi8U3VtDqOnFx/n8p3UEhCiVoubmyTgTMn69j3avJQnvrLUmOVkV9MzNlak+5Fb/6SgIltm2TfcriMQXNjChUeV45OTIf9eOPwLPPyra1ayUHXdVUrAoeegi4997AjrUuR4ulfKyg1XU+/BB46inf+5o1E1FLTfUtaOs8Zdr69dOBIikpTgvt+PGQWWzu9nOZmc4Yl7vvFqvqssvk/dy5sjTD6c25shEj9PoVV4gYqmhIZQX6C6OvDG64QaqJBEJCguS9mTlrFovFScgFjYi2E9EqIkpXZVaIaBwRZXq2pRPRRcbxjxLRZiLaQETDQj3eWk/PnmVn6+7fL7+qvlyOkyaJApjVR/r2lagJlbN21lnO2k1uVq2SkP9KIN8VK7V7txiRkZEyH3X8uBTyffRR2f/zzxLc6S/YwhQ6ldulRFMJWkpKpQz9pElIEItTNRi3WCzeVFce2mBmdhcOfJWZHRWaiehMANcCOAtACwD/JaJUZi6BpXJQ0QtuC23bNrHu7rsPuO46icQYP16KDX7xhWQUq8xiQIJJfDQ8LbpzLFC/Aer/NNtrX3kwA1Onikd02DBnqcnERLHQUlLkxz41VaIfhw7VdZhXrxatDQsT4Sst9f7TH3kE6N5dh84rXVe1G1UJzJpA27bVPQKLpWZT012OIwB8xszHmHkbgM0A+lTzmOombgvt119FAUaOlG2PPSZzZ308//wbNkgPFYWf5qA7VuVjy4oCn/t8UVwMvPyyGI4LFgBXXy3VtubO1YI2cKCUn8rMlEjGJk20hXXBBfKnqATlgQP1n5eQIOvp6eKaBIDnnpPPMP8ZGjaUPzU+3gZhWCy1ieoQNAYwm4iWEZFZvvVuIlpJRO8Tkfp1bQnAjEDI8GxzQER3ENFSIlparBKOLMFhRjwqCw3wdle2ayeW2MaNOL5ui95e6LvgNh05gtLDR33u88U770gTzTfe0P3GAIlGzM4WgZo3TxKrDx2SsPsmTYD/+z/g3/+WaEEi2Qfo+TSTbt38B3uEhel2LzXJOrNYLOVTHYI2kJl7ALgQwF1ENAjAOwDaAUgDsAfAy8FckJknerLXe4WXlwBk8U1qqnSivOoqySpesUL8eu5aTBEREimxYQPyftuqt5ttaQzqFx9B2PEjjuIjbrKy5CNLSnQZqeJi59Rbbq4IWtOm8l5FK65bJ1ZVcrIzwEJFBV54Yfl/upt+/WRZU+bPLBZLYIRc0Jg507PcD+BrAH2YeR8zlzBzKYBJ0G7FTABmxEGyZ5ulsomJAebM0e2VFy/2H+KXmgqsX4/43K3Yqf57fFhox44BDfgoIvioo8dYcTFw441iURUXS67XhAkyTbfFY/RlZoqg9ewp73Ny5OUWtGPHvKMfAWDhQukn5q82YlkoQbMWmsVSdRDRBUQ0iYjSPO/LaLgUGCEVNCKKJqIYtQ5gKIDVRGRWvrsMgKfVIr4FcC0RRRBRGwAdACwO5ZhPOVQZ+G3b/EdHduwIrFiBBiVHsQpdZFthIX574issGKZTBLKygCgcQSSOOiL7X3pJujM//bRYXgcOyMdt366PyciQoI60NJnHysnxbaEBvgWta1dnP7Fg6NdPvKpVmYNmsVhwG4AHAdxIROdBPHQnRagttGYA5hPRCogwTWfmmQBe8ITyrwQwGMD9AMDMawB8AWlRMBPAXTbCsYox/Wz+BK1v3xOrJwTt8GEcefdjtJ39zon86/37RdCicOSEoDGLoAESpJGTI+t5eVrQ+veX1mvZ2ZIV0KQJsGOHuCWVoJnV7H0J2snQvLkYqKNHV+51LRaLgwJmzmPmP0OMm94ne8GQTjgx81YA3Xxsv6mMc54F8GxVjsticOaZYv5kZvr3uRmTVaaFFpm3FzEowJYt4pXcv6cEDVAEBmHXLuDFFyVnTInYnj16/cABEa2wMNHLhQtle1qaCJaqwKUELSpKtpttYCoT5eq0WCxVxnS1wsyPENHYsg4OhEqx0IgovvyjLLUCImkzA/iPioiIAP74RwDAWkjl35KCw4g7uhfROIxV6WKi5WZKKH8EjiNzZwlefFEiGAGpnbhnj+49piy00093Tt317SuWnAqzV4IGaLejmXFgsVhqB8z8DQAQUVPP+zdP9prlChoR9SSiJ4moMRHFEFE/IrqdiF4hollElAlg+8kOxFKDuPtu8blddJHXrm+/lZD5Y6+9gx7t85EHeZbJ2VmIZpCiihuXH8LcucDnk3Vu2rb1x5CVJdU8AGDIEBExVWBYCVpKii7vFBMj0YpNmuikaLNPqRK0qrDQLBZLyHi/si4UiIU2AcB3AHYC2ADgacjk3WYAXQB0Z2ZrodUliIDevQEi7NsH3HOPVLtavx649lpZ7thJ2FsYg8OQ+Pjc9fvRCBLpuHVFAc45B1i+UOefLZ2v12NitEtvtSf8Jy9PXI6tW+sEaNXA0hQss1q+FTSLpU4QQIvbwAhkDm0BJBJlOYCGACYx8xcAQEQPesLvLXWU8eOBN9+UEPs9e3SN4r17VW1FadRVtEHnpO1aW4BXcD/OxrwT2wqytLV25pk6qGPVKllmZUnkf0qKaOlFFzmDRwDZFxOjx2YFzWKpE5SRpRoc5VpozHwPgNuY+RwAwwD0I6KFRHRhZQ7EUj0wSyPr3bvl/ccfi8vv3nslR+zYMdm+ezccydG7dokAHYGn19oWXTXkwM4CdMMK9MKyE9siIRZacrJUtVeCpiy0Q4fk+q1bSzPL6dPFtQlowXLHqJxzjtRhbN4cFoul9hJSCw3MfNizzAXwABG1hrgemxHRYGb+qbIGZKka8vIkMjAiwrl9wwapTh8VBQwfLgnPrVtL6anBg3UFLECqY6ki+ypIoxRhOIJINNynLbRGKEAjHHJ8TiSOIixMPi8iQqIaAb1U+MrlVsnRLV1Fz847T8L7LRZLrebRyrpQwFGORNSfSGqSM/N2T6j9AAAPE9HPlTUgS9UwYIDvtmjKsMrK0usTJ0pi8W+/SfWO/v2BTz4BXnhBC6JqngkAhYhGO2hBi0EBouGsHBKFI0hOliCPsDBxI3apvx6xOOg4zldg5UHPIW5Bs1gstR9mXl3+UYERTNj+zQCWEdFnRDSSiE5j5nRm/j2AcZU1IEvVsH2709pSKBHbv193fm7TRoIv0tNF0Lp0kQ4yo0aJq7FdO6egqcAQRUJ9b0FrHn/UIVbEpVjA/fA4nkFYmGyrV893A8uRI6XE5KOV9hxnsVhqA8GmhAWcWM3Md3o+oBOksPBkIooD8BOAmUQUZqt41ExKSqR2sHIXZmcDcXHSAdm00JSgNWsmCc0zZsg5qg8XkVTGatFCovoVhZ7AENSvDxQVIbVFAWJ2HnLMsN523REcG2QMas8eNCo+iDSkn6gs0rIl0KCB9/ibNpUWbBaLpe5ARD0BDAfwBoBiSN9L89UZEnUWsKgFnVjNzOuZ+VWPZXYegPkArgLwa7DXsoQGVTf44EERt06dZI4McFpoe/fKXFqjRiJoSgDdjSWbN9fBIv/4BxB7msdC69wZADC41yHEhjkttMsuPIprrzU2bBUX5RlYd6Iyvs887gMHpCijUluLxVJXqPSUsJOqFMLMR5h5BjOPZeZeJ3MtS9VR4OmvefCgFP3NyQHWrpXWKtM9xWeUhXbaaWKJ9e+vz+/a1Xk9s47ihRcCLTt4LLS0NKBePfRJzUP9YlcPtKOu9x5BS0YmZn6RD8BPcf+lS4FvvtG1sCwWS13BTAnLhKSEjWXmtwEcq0hKmG0edgqgBC0vT8+jLV8uc2SKrCyx0Jo1k/cDB0pEYlyc3qYwaxbHxAAnfIYdO4p558uacne03qqDSNIi1wPo49tCU7Wx3OGQFoulVsPM9xBRQ2Y+TERNADxORPcD+BsqmBJmBe0UwLTQlKCtNuKKevYEli0T661jR709NdX39a67TvLUAGntcqJ+VWqqKNyePd4n+bLQPHNuMRnr8MorfTBihI8PU0JmNlSznJIUFRUhIyMDR93fJUtQREZGIjk5GfXr16/uofhKCUsB8AwqmBIWtKAR0VUAZjJzARE9DqAHgGeY2WYE1VCUoBUWAps2yXpxsSzT04FFi0TQ1q/X/T3LIjFRigb/+qvHQnML2t693if5ErR+/SS6ZMUK3P+Kse/AAcmYfuEFLWRW0E55MjIyEBMTg9atW8OTQWQJEmZGTk4OMjIy0MZfe6hqhJl3ALiJiF4G8DwRjfMU9QiIisyhPeERs4EAzgfwHoB3KnAdS4g4ZOQ4m25GQCILVU9PwNu96I+5c8UlWb8+tDq2a+df0EyXI7Moa4cOImpz5zqPXbdOCjtec402Ka2gnfIcPXoUCQkJVsxOAiJCQkJCjbdyK5oSVhFBU6H5FwOYyMzTAfgItrbUFJSFBjgFLSJCEpzNCvaBClqDBoZLctYs4KGHpNmZ2+UY7QkYMW+gLVtk0q53bylH8ttvOqQScJ7/2WeytHNoFsCKWSVwMv+GRBRPRFOJaD0RrSOi/q79RERvENFmIlpJRD3KuNZNRDSYiL4kok+J6E73McG6HCsiaJlENAHANQBmEFFEBa9jCRGmoO3ZgxNh8i1aSESjKWK/+10FPmDoUInfB3SdKkV0tKifaaHNmSPLc8+VV2mp00pThSUBbV5aC81SzeTl5eHtt98GAMyZMwfDhw/3edyoUaOwdu1av9cZN24cXlKVt2sfr0OmnDpBmjWvc+2/EEAHz+sOlO296w3gYma+ipmvA9CpjGMDoiJCdDWAWQCGMXMegCaQ0EtLDcUUNEC8fIAuJdW+PTB5srgQ09JO8sPcpe+jouRlWmhz5oiKduyoB/Pbb3r/7t3iyzTF0QqapZoxBa0s3n33XZx55pkhGFFo8RTSGASZZgIzH/dogMkIAFNYWAQgnoj8lQ/PB9CUiEYT0ZVQrTtOgookVh9m5q+YeZPn/R5mnn2yA7EEB7OeuioPt6BdcIEslaARAbfc4j+qMShOP12vh4WJmEVGOgVt/nxg0CD54IgIOcYc5O7dkr1tlte3gmapZh555BFs2bIFaWlpePDBB3Ho0CFceeWV6NSpE2644Qawpx3Fueeei6VLlwIAZs6ciR49eqBbt24YMmSI1zUnTZqECy+8EEeOHMG5556Lhx9+GH369EFqairmzZP2SyUlJXjwwQfRu3dvdO3aFRMmTAAA7NmzB4MGDUJaWho6d+6MefPmoaSkBCNHjkTnzp3RpUsXvPrqq8H8ieFEtNR43eHa3wZAFoAPiOg3InqXiNwi1BLALuN9hmebL54AMA1iFDUAMDaYwfr8A4I9wUY51gzuvBOYMEG8deW5xN2Cdv75UhexSor9miJ02mkiVkeOaJdjcbEUhLzpJn1co0YSgrlnj2R6r1kj/tDYWMkAB6ygWRzcd593gNPJkpYGvPaa//3PP/88Vq9ejfT0dMyZMwcjRozAmjVr0KJFCwwYMAC//PILBg4ceOL4rKwsjB49GnPnzkWbNm2Q6/oOv/XWW/jhhx8wbdo0RHiqfhcXF2Px4sWYMWMGnnrqKfz3v//Fe++9h7i4OCxZsgTHjh3DgAEDMHToUHz11VcYNmwYHnvsMZSUlODw4cNIT09HZmYmVnvycvLy3AZUmRSXUyAjHPJ7P5aZfyWi1wE8AhGmoGF5AphGRE2ZObsi13BjoxxrAStXAnfcofOXAREzwDnd5KZTJ+C222QaKipKb+/eHbj5ZvjO+zpZTEHr0EEKMUZFiUjt2yeiVVrqtOSioyVIJDUVGD1acghatnQeo+p2WSw1hD59+iA5ORn16tVDWloatm/f7ti/aNEiDBo06ER4fBPDHT9lyhR8//33mDp16gkxA4DLL78cANCzZ88T15s9ezamTJmCtLQ09O3bFzk5Odi0aRN69+6NDz74AOPGjcOqVasQExODtm3bYuvWrRg7dixmzpyJ2NjYyvyTMwBkMLMqczgVInAmmQCMGxfJnm1l8X7lDK9iidVeUY5E9ExlDcjizfffA5MmAePG6bJTDRtKweHffpPf/hdflO7PU6bI/t27ZU5swwaJfm/RQtdtDAsD/vWvKhqsKUIffSTLHj0kya17d+C552SbWVa/USMJ0zfzC1q00PNxSUlSbDIvT8Iyi4uBcFsT4FSmLEsqVJhCFBYWhuJA5wAAdOnSBenp6V75YOqa5vWYGW+++SaGDRvmdZ25c+di+vTpGDlyJB544AHcfPPNWLFiBWbNmoXx48fjiy++wPvvV45eMPNeItpFRB2ZeQOAIQDc0S/fAribiD4D0BfAQWb2UWnBQaWFrp5MlOO1sFGOISFfSh0iN1c8dzt3anehiqV46CHgww9FwADgu+/0+QUF0pwT8FT2qEpMQWvZUl5PPAH86U8yj3bnnd7HNWqkQ/V79NCDVqLXrp0sc3OBn3+WgJH//lefv2ABcOutzpbaFkslExMTgwK3/74M+vXrh7lz52KbJ5fSdDl2794dEyZMwKWXXordZblZAAwbNgzvvPMOioqKAAAbN25EYWEhduzYgWbNmmH06NEYNWoUli9fjuzsbJSWluKKK67AM888g+WV3wF3LICPiWglpJDw34loDBGN8eyfAWArpMDwJAB/CuCalXbjVuQx92oAvwfwEjPneSJYbJRjFaLuodxc0YZp07Sl9txz0rxTdZPu1Enem+76zEygcWPgl1+cdRirBHfYPgDcfbcsDxwAPv1U1k0LLTpa13+87TYpNHnRRZ4yJJAwzIULpaqyZ7Idw4fL9aKigG+/lTDNl1/2jrK0WCqJhIQEDBgwAJ07d0ZUVBSalZO0mZiYiIkTJ+Lyyy9HaWkpkpKS8MMPP5zYP3DgQLz00ku4+OKLHdvdjBo1Ctu3b0ePHj3AzEhMTMS0adMwZ84cvPjii6hfvz4aNWqEKVOmIDMzE7feeitKS0sBAM8pj0glwczpANzzbOON/QzgriAvW3nJhcwc1Mvz4TcB+KvnfSsAfYK9TlW9GjZsyHWNW25hBpi/+oo5PFzWzVf9+szx8bLeuTNz69bMzZszX3utPuaSS0I4YPWhbl58UbZHRzOXlurtI0boc+bPZy4qku379zP36sX8r3/Jvtmzmf/6V33s4sVy3MiR8n7duqr/2yzVxtq1a6t7CHUGX/+WAAq5Gn6zAXSurGtVxFX4NoB+AK7zvC8A8M+T1FVLGZgux8aN9fYxY2RaqqhIrLOnn5Z5tG3bxCp70LCblbFTrfTsKcvkZGdoZrQR+RsTo+fHEhOBJUt0/5pDh2QuTbF5syzVtv1Bd5uwWCzVDDOvJqJORPSwp8rIG571M4K9VkUErS8z3wXgqGcwB2BLX1UppsvRFLRevXTVe0Ci5BVEzsacvXtX7RgdZGbqKsgm3bvL0pw/A5xuSl+TfGp/YaGIVrt28gcqQVPuSitoFkutg4geBvAZxPu32PMiAJ8S0SPBXKsic2hFRBQGz0QeESUCKK3AdSwBYlpo8Ub/1qZNnb//bpe+eeztt1fd+LwwO4CaxMeLCrvLkbgtNDdK0JSF1qoVcPy4t6DZrtYWS23kdgBnMXORuZGIXgGwBsDzgV6oIoL2BoCvIf1qngVwJYDHK3AdS4AoCy0nR0LuFU2bOq0yX3PUTz4JpKTUEJcjIJEp5h8BOC00XwNVgqcETVl6mzfLbJqyzPbtA378UfLcVDkUNx9/LNf7wx/0tmPHpN6kLXxrsVQHpQBaANjh2t4cQRpLQQsaM39MRMsgOQgA8AdmdheotFQipoVmRg03beoMFjTFTTFuXJUOLXga+PBOK8GKjPS9X1VTLiwU0UpKEtP0m28k4fr4cdn/9NPyAvyH8N94o/f+tDSJllywQMZgsVhCyX0AfiSiTdBls1oBaI8gy2FVpPRVBCQ7PM5z/lVEBGb+W7DXsgSGOYdmClpiopRCbN5c0rjMNjC1CmWh+TMjw8JE1HJyRMCSkkT49u/Xbkc3R444y6MAvotfFhVJ0jcgORBPPVWxv8FisVQIZp5JRKkA+kDXfcwEsISZgyoPVJGgkG8gFZWLARQaL0sVUFrqW9Dq1dNzZK1aAXFxtdi4UBZaWVnf0dG62WdSkqg4AHhq1p1wF9bzfKVzcryvsXOn9zZz3q3yk1AtdRyzELHJ5MmTcbfKv7SUCzOXMvMiZv6357WImUuI6NZgrlORObRklk6ilhBQWKi9Y0rQzjlHWpCp3+6zzqrlZQ7Ls9DUMVu3ynpSkv7jN26UpfpHuvhi4D//AbKznf5YwLc1p6o0REToa1kslprCUwA+CPTgigjaAiLqwsyrKnCuJUiURaaqQxUXi5j95S/6mNdek7iGWksgFlqjRlqQmjXT82ZuEbrwQi1obnwJmiq5de65ElBSVCSltSwWF4WFhbj66quRkZGBkpISPPGEs8j8Bx98gOeeew7x8fHo1q2bo9ajxT+eMlo+dwEouxyLi4oI2kAAI4loG4Bjng9lZu5agWtZymDdOuA6T/p669bau+Y2ZGJialAUY0VQFlp5LkfVgiYpSSu9Eqnx44HHHgNU+46cHImmeeIJKWj81786BU0Jl7LQzj0XmDVLrMCOHSvtT7NUEdXQP2bmzJlo0aIFpk+fDgA4ePAg3nlHGo3s2bMHTz75JJYtW4a4uDgMHjwY3VU0rqU8mgEYBuCAazsBWBDMhSoyh6ZabA8FcAmA4Z6lpZJ57DFgxQpZ72Q0J6/V4uULZaGV53JUJCXJpCEg82qxscAf/yhWmcpdyM4GHn4YeOMNyV0AnIKmBHHPHnFfnn22vFfVnS0WF126dMEPP/yAhx9+GPPmzUOc+g4C+PXXX3HuueciMTERDRo0wDXXXFONI611fAegETPvcL22A5gTzIUqYqE1ZeZl5gYiGg7vHAKfENF2SLmsErgayhHR/wF4CUAiM2cTEQF4HcBFAA4DGMmnUCPRxES93qULMHWqrNc5QQvEQlPHREbKuqfyOPLznRWXVXHi7Gz9NKDOzcjQx+Xny7G7d4sInnmmbLfzaLWDaugfk5qaiuXLl2PGjBl4/PHHfXagtgQPM/st+8DM1wdzrYpYaJOIqLN6Q0TXIfiOpYOZOc0lZqdDrD4zFE1Zgx0A3IFTrJGoOS/WubNer7OCFoiFlpQkEY2m+JkV9sPDJfwzJ0dHRR46JHNu+/bp81Ry3+7dUtmkcWN5grAWmsUPu3fvRsOGDXHjjTfiwQcfdLRm6du3L37++Wfk5OSgqKgIX375ZTWOtHZBRP09xstJUxFBuxLAFE8xydGQfjdDK2EsrwJ4CM7eOCMATPEUg14EIN7TruaUwIwoVwYE4LtDS60mEJejOkYl24WH620JCc5jExIkRH/vXh3pmJ0teWsdOsh7JWh79ugUgNRUK2gWv6xatQp9+vRBWloannrqKTz+uC6Q1Lx5c4wbNw79+/fHgAEDcMYZQdfVPZW5GcAyIvqMiEYSkY8SEYFRkUohW4noWgDTINbUUGY+EswlAMwmIgYwgZknEtEIAJnMvMIl1C2hM8cBaQHeEkB5HVBrNUeOAF9+qQPwAGd5xDpnocXGAs8+C1x5pf9jTAtNERcneQ1uQWvaFFjm8Yr36iWuxi1bJES0fXvZl58vof47dgD9+smxHTs6O6NaLAbDhg3z6ho9Z86cE+u33norbr01qLQpCwBmvhMAiKgTxCs3mYjiAPwEYCaAXwJNsA7YQiOiVUS00hNiORVAEwBtAPxaRtilLwYycw/PwO8iokEA/gLgr0Fcwz22O4hoKREtDaYNek3lq6+AW26RKaDbb5fkatPDVucEjUjyEFJT/R/jT9AA76aeTZvq+bJeHq/2Wk+neGWhvfWWNAU9cED7czt2FCsuKwvo0wcoo+miF0ePAosWBX68xWJxwMzrmflVT57zeQDmA7gKwK+BXiMYC214kOPzCTNnepb7iehrAOdAhFFZZ8kAlhNRH0j5E7PXSLJnm/uaEwFMBIDo6OhKa+ddXaj8YUDqM7q9y3VO0ALB7XIEdKkUXy5HhRK0dZ5yo0rQvv9eC5bqt6bC9efPlz5sS5b4L3LsZvRo4KOPZE6u+SnjFbdYqgSP12+G5xUwAVtoKpQSwN8AHDTe5wN4MpBrEFE0EcWodcjc2xJmTmLm1szcGuJW7MHMewF8C+BmEvp5PrfOuht37ZJCF4sX621mwWHV3+yUFLRgLDSzkr6ay1AWWvv2ep+y5k0LDQB+9TwQFhSIebxvn/9ix4rZs2V5+HDZx1ksliqjImH7XZk5T71h5gNEFGgGYTMAX3sssXAAnzDzzDKOnwEJ2d8MCduv0w7q994DZrieR0xBmzsXWLjQd0H6Oo8SNLNHjhI0t4V22WUS1r1li7gfAW2hmV1PAWk2qrqmtm0r5rCaf9u7V84/cACYPl3ckdddJ2Wy3OR5bon0dPEZm+3CLZUGM6OSAuJOWbi8h7NaTEUErR4RNfZ0qgYRNQn0Osy8FUC3co5pbawzgLsqMMZaxZw5MmXjK9JX/R4DQMuWZcdN1Gl8uRz9WWiAs5V3ZKTMqdWr50zuCw+XBD9FgwYibioBe/16ETMA+Pe/gfffl/0jRnh/nirFNWmSVBwZM+YUNaWrjsjISOTk5CAhIcGKWgVhZuTk5CCyhlYy9zSMBjNnVeT8igjaywAWEpH6+b0KwLMV+XALsHQpMHiwVGEqMvq1xsVJp5RWrapvbDWKnj2B/v11c0/Av4XmRiVQN23qbC76+efe/8AJCWLZAc5E7B2eugEHD5b9WaqUVkGBFbRKJjk5GRkZGcjKqtBvncVDZGQkkt2Fu6sRTw7akwDuhkyDEREVA3gz2LZkFQnbn+Jp8DnYs+lyZl4b7HUswpo1siwtld6TJSXAp58CL74onjPTQjulSUmRBpwmgQpaQoKuCGJy+eW+j920SdbNvIlMTyyS2ZBu5UrgmmuAefP0NnWOeZylUqhfZJgHlQAAIABJREFUvz7amFVhLHWF+wEMANCbmbcBABG1BfAOEd3PzK8GeqGKWGhg5jUA1lTkXIuTjRvF83X4sFhp77wjgtamjRWzcklNlXwGX626TZTI/N7T9ejLL3WEpBvzH93sybPLkw5ZUCBRO717i3m9fr1MbipUlX8raBZLoNwE4AJmPtEiw5PvfCOA2ZCiGwERsKAR0XxmHkhEBXBW81DV9ssoxGfxx4YNEougOpZcdplEi6tcX0sZXHmlhIU2bFj2cc88I1GOzz2nz/OHL2uPSBK4AQn4ePRR4NVXxawGdBCJiRU0iyVQ6ptipmDmLCIKqpdTwILGzAM9SzsxUIls3OjsVnLaaRJ7YAkAovLFDJBK/IHiS9CSk7WFpubVXn4ZGDlS1q2gWSwnw/EK7vOiIrUcLZVEaalM15RVIMMSYtx+XiIJL1WoApsZGbru45Il3texgmaxBEo3Isr38SoA0KXcsw2CnkMjoggAVwBobZ4fbDSKRR76jx61glajcFtoMTHO+TblZgR0MEhurliKzLoJ6aFDEvGzfDlw0036nI8/lny2+++vmvFbLLUMZg7zt4+IgrpRKmKhfQOpgl8MoNB4WYJEdTdp1656x2ExcAtaXJz/8Pu9e/V6u3bO4woKpLno6NHOc268EXjggfIrj1gsFkAiIAOmIlGOyZ7ikZaTZPt2WaakVOswLCZul2NsrHfjUSJvQWrXTgJH9u+X9wUFkipw7JiY4e5E1owMqVJisdQiymrQ7NnfGMD7ANoBOArgNmZefTIfGczBFbHQFhBRUH5Ni2927JDfRvu7VoNQFpoqbxUX5y1oKSne9cfatXM2qiso0Llr//d/wMCBzuOXnzKN1y11D68GzQZ/AZDOzF0hfc5eP8nPCsqVEXT7GAADIdXwN3jayajtliDZsUMKs/sqDWipJpSgqUoKvgStcWNnCS5ABE2V5wK0hQZIbbMVK3QxZMC/oP30k/ZFWyy1jzMB/A+QdjAAWhNRs7JOIKICdzCIegFoUda5boJtH0MAXoD0MjsxHs82S5Bs327djTWOpk2lHFb//lICy5egxcdLXciMDMmKV41DTQstN1e7H7duFbejOee2fLmU0SJyXv+882Rp59gsoSeciJYa7yd6WnOZeDVodu1fAeByAPM8LcBSIG2/9vn70MpMBQsmD20HABBRe7Wu8HQatQTJjh1A377VPQqLg/r15T9m7lzpbxYb6x0UEh8PREXJer9+UjGkc2enoG3apEXp6FFZbtyo9+/aJddp0EDm2SyW6sdrTswHA5k5k4iSAPxAROuZ2SiVg+cBvE5E6QBWAfgNMt/mFyLqDWCXp2UYiOhmSCT9dgBPMXNuoH9AMC7HO4loFYCOqnO157UNgHU5BklJifymWQuthqJEzLTQ6nlul/h47XK85hqxxJo3dwqaylEzWb9elikpQH6+rB838kbNlID583XOW1ns3SvVUnIDvuctlgpjNmgG8DWAPq79+cx8KzOnQebQEgFs9bqQkwnwJFAT0SCIKE6B9Np0W4BlEkxQyCcALoE03bzEePVk5huD+VCL/A4VFVlBq7H4EjRVMzI+Xhc6btxYH6vm0BITnUKlUCLXoYMWNMDbkgOAs8/23S27sFAamKqOAEuWSBO99PTg/j6LJUj8NGhe7TomnohUxNQoAHOZOR9lE2ZYYddAXJ3/ZuYnALQv4zwvgulYfZCZtzPzdapbtedlHw0rgKqX2yKoKU9LyGjcWCyyZs20oKmKIY0ba0Ezk66VhebvP1UJWvv2TkHLyZGlu9v1Gh/1v9euBb75Bpjp6Yurakzml/ebYbGcNM0AzCeiFQAWA5jOzDOJaAwRjfEccwaA1US0ARJrca+fa5mEEZGa/hoCT1CJh6BSyypUbd9y8qh4AXewnKWGkJAg7Wq6dtX1G9u3F4soPl5HQ6pu14AWtJYtJaoxMlJ8y6rR3YYNMmfWqpWz+d3OnRKM4hY0X41LlWtRRUK6BY1ZAk0slkrGX4NmZh5vrC8EEGzto08B/ExE2QCOAJgHSLwGgHIaEDqxtRyrCTU94m7RZalB9O0rwR8dOgALFwLXXivb4+PFHThqFJCWpo9PSJAGoqqO45gxTsHbvl2eYFQfN8XOneKiPHTIud08V6Gsua2eaQklaAUFwIcfilWpnpYslloAMz8L4P8ATIYEnagQ33oAxgZzLStoIYZZqumvWiXvrYVWS+jXT7sX4+NlnmzSJGe1/1tvlfqOt90m7594wrvvWnKydxrAzp2SjKhC9hW+BM1toSkRzM/XZbbMTtsWSy2AmRcx89fMXGhs28jMQVUgsC7HEDNzJnD77bIeHe3MxbXUcPr3B559Fjj/fN/7o6PlmP79gVtuEdefW9A6dvQWNNUh2x3VqFIDTJSgbd0qT0fKQsvO1uH/ZnCJxXIKYS20EPPPf+p1626sZdSvD/zlL76Fxo2ax3JbWarLtsns2b6v4Z5TA7TLMT8fOHBAC9qMGfoYt+uyLF56CfjrXwM/3mKpwVhBCyGHDgHTp+v31t14ClCehdaunTPh2sSXoJn5Ztu2aUFTOW6A3hYI//63vCyWOoAVtBCi5upVXVtroZ0CJCXJfJsiNdUZFDJihP9zfQlTbq7+Au3Y4fuYYARt716x9CyWOoAVtBCivEX9+8vSWminAI8+6jTL27d3Wmh/+IP3Od26SbCJP0Fr78k1zc727V4M1OXILIJmq4xY6ghW0EJIdrYslaBZC+0UoHlzoHdvCf0HZP7NFLQePSQvzeTHH4E779Qux2XLdOWRnBx9rZwcp+ipL1SgFlp+vgSQHDumO21bLLUYK2ghRFlov/udLK2gnUIsXqzD6SMixG0YFiaW2IwZktOmaNhQXocPA3/7G9CrF/D447IvN1cSt6Oj5QnJFK/WrWVpWmivvy6fZdaJVJjV/62VZqkDWEELIaagvfsucP311TseSwiJj9elswCZR4uLk2jIs84Cuhg9cyMjRbCYgSeflG0ffCBW2oEDUkEkIcFb0Jo3Fwvw4EEdaHLffVKVxJfVZgqanUez1AGsoIWQnBwp5NC4seSi+apsZDlFiI11BocoN2TDhiJyZsL25ZeLeH3+uVhaTZpIqSy3yzEpScpvvf++tLMxRcpXrUe3hfbEE8BrrwX+N8yaZVvfWGoUVtBCSHa2rnlrOcWJjXXOpZmCBjgz7q+5Rpbz5skyIUEEzR0UkpQk5x08KFaZcgkAwLp1wP/Mmq/wttAmTgQmTw5s/GvXAr//PXDXXYEdb7GEAPvTGkJycuR3yGJBSoqe8wLKFrTWrWUeTLWMcbscwz0Ff5SFpjDF7oILgCFDnBadKWgbN0peyYYNvufb3CjLbMGC8o/dvBl47DHbhdtS5VhBCyE5ObpIu+UUZ8oUKSascAua6XJs1kxM+82b5b1yOWZlSeCIii5SFppC9SgyWbxYr+/dqz93/nxZHj0q+W3loQTNtAL9MW0a8Pe/OwXUYqkCrKCFECtolhPExOjGoIAWFlVWyxSmxEQRsV275L1yOap5sebNZem20Nat8/7cX37R65mZkugdFqbdmf7Oc2PWkCwP1X1ALS2WKsIKWgixgmbxS1kWWsOGImjKZadcjoohQ8SlmJbmFMK1a70/xxS0HTvEndm4scyhRUTIdlPQtmyROTk3KkcuEPek2RHAYqlCrKCFkOxsO4dm8UNZc2iAs8hx48bOL1KnTlLguHFjp4XmS9AWLhRhZJa2NSkp2lLs0UOsQSVopaVSleSii7yvY87FlSdqStDcFtrcucBbb5V9rsUSBFbQQkR+vhRjsOWuLD6JjnaG67sFTeV4xMZKEIgpaOax5rq76DGRWFvbtsn829GjImiqt9p990m/NtXGZudOWfoK/DALJ2dllf23+bPQ3nsPGDeu7HMtliCwghYiVDH01GCbk1tODYhErHy5HAFtoSlXY7duep95rGmh5eQ4r6kaiP72mw78SEkBXn1VohCvvlrOVwK0erUs3e1uAKeFpub2/OFP0PLyfHcUsFgqiBW0EKG8OGecUb3jsNRg2rUTgQH8W2hq2bQpMGeOuAk7d9bHuc9LSNCCdM45EgDiFrT77gOeeUbeN2qkxUoJ2mmneY/VFDQVfekP5Wr0JWhHjjhdlkVFQHFx2dezWPwQckEjou1EtIqI0oloqWfb00S00rNtNhG18GwnInqDiDZ79vcI9Xgri3XrJJWobdvqHomlxjJ3rtRuBLRV1aaNLN2CBohALVumRRBwWmjqeCVoLVrIfJtb0Eyio7VFtWaNLEtKpIrIr7/q45Rl1bRp+f3U/M2h5eXJUhVGPnxYxnvFFWVfz2LxQ3VZaIOZOY2Ze3nev8jMXZk5DcB3AFQL3QsBdPC87gDwTuiHWjmsWydF0lUOrMXiRXS0dMUGxJL6z390VKLb5VjWNUwSEnSJrSZNgJ49JTBk9WoROncDUl8ux8xMseD69dPHFRbKZ914I/Dtt2Xno5XlclTXAsTteeiQXM9iqQA1wuXIzOY3PRqAKikwAsAUFhYBiCei5iEf4Elw5Ig8SH/7rW5jZbEExPDhOsfMl4Xmi06dRKhatNDHKwutSRNg7FgJ0Z88Geje3ft80+W4dassjx7V+5UoHT4sVuQ110jR5GnTpEDp9u3e1yxP0JS1t3Sp/hsslgpQHYLGAGYT0TIiukNtJKJniWgXgBugLbSWAMwZ5wzPNgdEdAcRLSWipcU1zP/+j3+IJwnQbWMslqBRFlp5gnb++SIUp58u7805tCZNpBXNmDFAx47ARx95n69cjkePigC5exx9950slYWmXJbffCNFkdu08e6t5ha0wkLJbzPFEdDRkjZQxFJBqkPQBjJzD4g78S4iGgQAzPwYM58O4GMAdwdzQWaeyMy9mLlXeA3z6b3xBnDZZVL15/77q3s0llpLoBYaIBGTai7NFDQlim+/LTlqycne5zZqJBbX7t3y3h3FtGyZLAsLxUJT6QNmMvbChbIsKpIKJGruLDtbLLinn9ZNStW11H4gNIJWVGRLcdVBQi5ozJzpWe4H8DWAPq5DPgagZoUzAZxu7Ev2bKsVFBWJdyctTR501fSIxRI0KSniglSh9+WhkqXdLkdABM9fywclhMrdaLr/2raV4sVXXAH8/LOe84uL04WTAS1M99wDDBqkoxa//x4480xpO2MWKj58WAJPVJPRUAjamDHizrXtb+oUIRU0Ioomohi1DmAogNVEZDyuYQQAT9YWvgVwsyfasR+Ag8zso+JqzURNEZhFHiyWChERIUEiZv5ZWZgWWrduIiSqTmRZqKASJWjKQmvaVK7x44/AV19JEIiKxExMdAqUChDxFf145AiQnu7cVlgoYsYsQSqHDwdXmX/qVN/lucri669laetL1ilCbaE1AzCfiFYAWAxgOjPPBPA8Ea0mopUQkbvXc/wMAFsBbAYwCcCfQjzek0L1V7SNPC0hR1loCQkSrLFmjVhm5aGEUFUPURZacrJYaGaAiBI/5XZMTJRlTo4IUnkVRBSHD2urrlUrWZqf44YIePRRWV+7FrjqKuCPfwzssxSqbqWtL1mnCOmEEzNvBeD1iMnMPhNPmJkB1NoOgkrQrIVmCTlKmIJ9mlLnKReistBattQ5cQolaErITj9dxCknB9i0yXlsZKR/kSos1IKWkgKsXCnX8WVRKvfl888Dzz2ny3RlBjkTERkpy4wMmTO0kZV1ghoRtl9XUVMCVtAsIce00ILBdDk2bChCRqQtNBPT5QhIodKEBBG05cudxyox69DBO1fu8GFtzSkLzd88mjuCUt1kZiueQFCCds451Ve+JysLmDSpej67gvgqjOHaH0dE/yGiFUS0hohuDeX4alZIYB3DWmiWaqN5cylNE2w1bNPlmJQkwSPPPgsMHuwtRL5cjkrQ/PVJ+/BDOd5MyjQttPIEzb1dNTENVtCUy1HBHJhLtjK56ioJrjn/fG/rt2YzmJn9NcK7C8BaZr6EiBIBbCCij5n5eCgGZgWtCrGCZqk2br5ZIgx9FRYuCyVoubk6tF7NVxUWinAkJuqgEcC3haYsp/r1Jdz3ySclmTstzVtMNm/Wofu+BC03V9yPUVHegqbSC5TFFSju448eDSxopjJR85R1CwYQQ0QEoBGAXAAhSw62LscqxAqapdpo0KBirR1MK8xt3UVHy4/wE0/Ie/UFVxaaW9BiY3Vi9rXXSg6aErOzz5ZlZCQwYYJO8lZzfqZwJSQAZ53lvX3XLpkDA5zFkgPBl6CFGvWZx0NivARCuCpQ4Xnd4eMYn4UxDN4CcAaA3QBWAbiXmQPoAls5WAutCjlwQH4DGjSo7pFYLAFiFjdW5bNMEhK0gClB82ehNW4M3HAD8Pe/ez/VzZolIfNnnOEUEzUv57bEtm2TwA9zDk1Zc0Dw4fduK7E6BE39Le55weqj2Kiv64+BzJxJREkAfiCi9cw819g/DEA6gPMAtPMcM89V3rDKsBZaFXLggLXOLLUM00Lzl/OmrC71pKaCRdq3F0HLy5M5sSZNpCpIZqZ3Ca2oKBFAs5fbo496C5ppvUyc6H9uzRS0Q4ekLI9KBPWFO8/NrPivhLqqUSJaHWJaQQIojHErgK889Xc3A9gGIGQhpFbQqhAraJZah2m5+CpeDEg9yFdeEYEBxMrasQMYOFAEjVnC/ps0kaASX5aeQpXPGTdOLDm3oJlV/N9803+5KlPQZs8GXntN6kv6w10hRIlKWlroEkeLipyfXcPxVxjDddhOAEM8xzQD0BGSSxwSrKBVIcrrYrHUGsxIv65d/R9z//3/396ZR1lRX3n8c9maVRBQD6MtyiKEOEaQEDUmuBwRUVwSTdQk6oyaQTNq3EZckmg0xgWX48TlmGiIc1TcxoxRY7txNJKjGJVmU6EVThDQRgRXELXv/HF/P6ve473u10332/p+znmn6lX9qt6trtfvW/f+fr97M72uGP6L0wQaGgoThuhFxWTK2YIWRz9On25PiDfckPs8UdBUYcECW3/ppfyfmy1oGzbYHLc4f641mUq2lAoRNPIkxhCRaSIyLbS5DNhbRBYATwPnNzMist3xPrQOZN06L+jpVDDpcGChxOrWTU2FCVoM70VBjCHPbA9t0iQb8h/nt02aBBdeCH/8o4nQ66+bZ3bCCck5WiNoGzdmtl+3rnieWoUIWjOJMW5Nra/CPLeS4B5aB+IhR6di2WOPth03alSy3povf0se2uDBlgQ5CtEdd9ik6JkzbXrC++/D1KkWkoxZTurrkz441STLCGw+snDDBhPESGszj+Rj2TKb6hCrf+eiQgStEnBB6yA+/9zmfG6/WfU2xylzNm5MSsC0lnRJmtZ4OFHQ4nD6bA8tXXkbMueMxUnVmzYlnt6YMSZ+M2bAvvta32BNjfXzbdyY20Orr0/ex+kALXHFFUlh0lwsXGiDVBZmdTXF/rP42U674ILWQSxbZhUx0mWfHKciqKlpe60jkWQYf2sELXpmIrae7aFlC1o6HJrOEhLnyF1yiY3SvOgi89hOP93mws2ZYyVwcvWhrViR5HQsxENTtfN/85uZ2zdtSvrgYq7J9OAWyByB6YLWbrigdRCxb9kFzel0xMEihQjaiy/C3XdnbuvdO5kovXatzY2rqbHSMmCilx6NmRa0n/zEPK2jjoKnnoKrr7ZyNddeC+ecY23eeiu3h7ZiBUyYYOcvxEPLJUTr1pn4Pv64vY+Clp0KzAWtQ3BB6yBc0JxOS/Y8teaYMAGOPTZzW9pDW7s2mcgdPbTevTNHY6YFrabGRmeK2HHnnZeMvIwjtHIJ2gcfmPgMG2b2F+Kh5ZoTt2KFhRgXL7b32YKmahlT0nPdymdidcXjgtZBLF1q/3/xf9FxOg0//KEta2ubb5eP7JBjFKQoaNk5F6Og5avCHRkwwAaqREE76yx47TXb19CQ2Lz99oV5aLkELYYWo2Blhxyvv94SEaf7KN1Dazdc0DqIpUvNOyt2Am/HKTmnnGIeTr55bC0RBW3BAnjssdweWprWlMoZNsz61D77zEKZcRBJDKnU1tq2BQtazrHYGkGLHtrMmbaMAgouaO2IC1oH8dprbcsN6zhVQXPZQQo59tlnLSMJJFMB8glaHBlZqKC98Yat19Qkxy5ZYsvaWqvwvXIl3H578+dqi6DFKQWNjckxzQnae+9ZGNMpCJ9Y3QEsX24Riz33LLUljlOB3HwzfP/7lq3/lltg+HDbHgeFZAta9OB+/OOWzz1sGNx/v63X1FiYskePzJDjqFE2vP/KK03oevSwPJPdu2fmukwLWlOTnaulkGM8JvbR9e3bvKDtvru1LWbmkgrGBa1A1q2z79Wuu26+74sv7LscQ/izZ9tyv/2KZ5/jVA21tTb6MTten68PrbbW6qLFLCXNkc7QH0dK9upl4cWtt04E68wzrQDnDTdA1642WvKgg+Dhh82uoUNh3LjkXOvX26jOtKBt2pQIW2Njpse3cqWdZ9CgTEHbsMHONWRI0g7MU/MO+RbxkGOBXHWVlXDK9aA0diwceGDyfvZsm4oTSzg5jtNKcnU+5ws5gglAIR3W6bBkFLQYdkwL4uGH22jH3r2Tf/q6OhOx7bazrAmPPpq0j0KWFrQ1a2x9xx1NqE4+OWm/apX1/fXqlSloJ55oIdfs/rtcabx++1t44IGWr7kT4YKW4p//tPmX6Un8kdWr7cHpgw8ytzc0WBKAZ55Jtv3tb5aVxweEOE470pygFUp6bly2oKULmnbvbgM47r4bLr3Uwo5du1oKq1zlZWKF7rhct87yTEJm1YLf/c6WmzZZSLVnz0xBi/PX4ijI6I3Onbv5Z954o+W3BJtrFz+7E+OCluKYY+z7liuTTZwHuXp15vZ0FOG99+x7vHx50p/tOE470R6Cls4vmQ45wuY12yZPNk/t4ovhD3+w1D/5yPbQ1qyxrCU/+AEcdphtmzgRTjstedLNJWixBl1dnfXLxflyc+Zkfp6qCVhjo62PHWseZifva3NBSxG9+hgpSBMfylatytz+178m64sWwfz5tp6vNqLjOG0kDgrJ7kNrDYV6aNmkky6vWAEHH5y5f+1aeOSRZBRj5Be/SKYvnH++iVmsCh4FLT2xOoaH6ursR6epydo9+aRlPol88ol5eY2N8OGHybHpMGgnxAUt8OGHSTLuXHMqc3loGzeaiJ1wgr1fuNA8f7DBSY7jtCP9+iW5HttKLkHr2tWWhQjawIE28TqdVxLg73+3bP/vvpuZIWWXXSxc8+mniQjGeXO5PLQYNly8OHmyvu46G+l5ySWbt2tszPxRev31/NfQCXBBCzz7bLKea9pH9NDS3535800Ep061h8dFiyyN3LbbFjbgynGcVtClC0yZAnvt1fZzbLVVsh4FLRYHbU7Q+ve3f+pvfMNENX2evn3hzjuT9+mKA1Hc0l5ltoeWFrQYsty4MfH2amth//2TuXLpdh9/nLSrq4Nzz81/DZ0AF7TAyJE21WS77RIP7S9/sTD6l18mHlo65Pjyy7bcYw8LYT/9tIW63TtznA7ikUfguOPafnw6PVZrBA3gppvgsstsPd2fN22aeWDR04vL7D65SPTQskc5NjXZk/POO9v7WHJmm21s25o19gP00UeZA0BiWCg9JaGT4oIWGD3aShuNGJEI2mGH2UPPqlX2IASZHtrcuRaBGDrUkgssWWKvQuZ3Oo5TYqKgxT6olgTte9+Db3/b1qOH1rs3nH22eV1XXmm5Gi+/3PblK4aY7aEtXWrnfucdE7U432fBAlsOHgw77WTrEyZYH0e6HE2s4xbnrnVifGJ1FrW1NsoxfschMywdBe3tt2HWLBvEJGJzMM87zx7ctuQB0nGcIhEFLZaqyedR5SLtoQ0ZYnN++vc3D3D9ejvXtdfmPjbdhxbnmz30UNI3+PWvmycaPbTBgxOvbeVKePXVzImv9fV2bDoM2klxQctihx3gz3/OHCwUE3J36WIPTSedBM8/bw9Tl15q+3r0sJBjz55JxMFxnDImXVMNWvbQ0qQ9NMicDjBggHlb+Uh7aDH0A3DXXbZMe2h9+lhYMgoamHimQ0VLltigEZ/46iHHbHbYwULaV1+dbIuljaZMsWVdnYW1Z81KIgEAX/ta5vfOcZwyJg7YiMt0XbWW2JI5cWkPLfbppacFjBhhE7ubmhKR3XbbZGBJU5OFkdJVxT3cCLigbcYhh9j3e948mzYCiYd2/vnWF/v22+ahHXlk6ex0HKeNxBBK9NDq6+Hee1vn4WR7aK0hemj9+lltNoBTT032DxqUhD8nTLClSObT8wsvmIjF3JNbUt2ginBBy2LECEsMUFNj37Gtt04ELR1VcBynQokiED2z0aOtM7w1tJeHFnPp7btvsn/QoGT7AQck248/Hs45x9bXrbN2kydDt27JYJVOjgtaDi6+2OaiDR1qD0FxfmNMVOA4TgUTh95HT6kttIeHttVWcNtt8JvfZBZDHTAgmUqQFrTp02HGjOTJeuBAS078+edwxhmtt6MKcUHLgYj1kUHmBGn30BynCjjjDMt5GFNetYX28tCGD4cLL7QfnUWLTOC6drW5bSJWvy2biRNtOXJk22yvYnyUYwvEvtY42MhxHGeLPLRRo+y4dEYRgDFj7AVW2PTmm3Mf/+CDNg8tncbLAVzQWiQm2P7Vr3xUrOM4gZ49rQ+uLU+5++1nc9XSoxRzke8Hp0uXJITkZFB0QROR5cBHwJfAF6o6XkSuAaYCm4A3gX9T1fWh/QXASaH9GapaV0x7p083z/7ss4v5qY7jlD0zZsDee7ft2JbEzGkTokWunxMEbbyqvpfaNgl4RlW/EJGrAFT1fBEZA9wDTAD+BXgK2EVV8xYm6tOnj34SZ/47juM4BSEin6pqn1LbsSWUxaAQVX1CVUPxFl4AYnD5cGCWqn6mqsuABkzcHMdxnCIjIstFZIGIzBORzUohi8h5Yd/Z4/CUAAAID0lEQVQ8EVkoIl+KSNE6+0ohaAo8ISIvi8hPc+z/dyCWzdweSBdzeTtscxzHcUrDfqq6u6qOz96hqteEfbsDFwDPqur7m5+iYyjFoJB9VHWliGwLPCkir6vqcwAichHwBXBXa04YhPGnAD3SxfUcx3GcUnEs1mVUNIruoanqyrBsBB4ihBBF5ETgUOBHmnTsrQRqU4fvELZln/M2VR2vquO7dfOBm47jOG2gm4j8I/XKFUFrKcIGgIj0BiYDD3aUsbko6q+/iPQBuqjqR2F9EvBrEZkM/BcwUVU/TR3yMHC3iFyHDQoZCcwtps2O4zidhC9yhRGzyBthy2IqMKeY4UYofshxO+AhsfkV3YC7VfVxEWkAarA/EMALqjpNVReJyH3AYiwU+bPmRjg6juM4HUc6wiYiMcKWS9COocjhRijBsP2OxoftO47jtJ6Whu3niLA9CfxaVR/PatcfWAbUqmpRf4y9w8lxHMcphHwRtmkAqnpraHck8ESxxQyq0EMTkSZgQxsP74aFNisZv4bSU+n2g19DOVBs+3upalnMTW4rVSdoW4KI/KOATtGyxq+h9FS6/eDXUA5Uuv2loKLV2HEcx3EiLmiO4zhOVeCClsltpTagHfBrKD2Vbj/4NZQDlW5/0fE+NMdxHKcqcA/NcRzHqQpc0BzHcZyqwAUtICKTReQNEWkQkemltqdQctUnEpGBIvKkiCwNy61LbWdERO4QkUYRWZjaltNeMW4M92S+iIwrneUJea7hEhFZmaoFNSW174JwDW+IyEGlsTpBRGpFZLaILBaRRSJyZtheMfehmWuoiPsgIj1FZK6I1Af7Lw3bdxaRF4Od94pIj7C9JrxvCPt3KqX9ZYuqdvoX0BV4ExgG9ADqgTGltqtA25cDg7O2XQ1MD+vTgatKbWfKtu8C44CFLdkLTMFq4wmwJ/Biqe1v5houAc7N0XZM+D7VADuH71nXEts/BBgX1vsBS4KdFXMfmrmGirgP4W/ZN6x3B14Mf9v7gGPC9luBU8P6acCtYf0Y4N5S34NyfLmHZkwAGlT1LVXdBMzCqmVXKocDfwrrfwKOKKEtGahl5s7OwJ3P3sOBO9V4ARggIkOKY2l+8lxDPsqu6rqqrlbVV8L6R8BrWOHcirkPzVxDPsrqPoS/5cfhbffwUmB/4IGwPfsexHvzAHCAhBxUToILmlHJlbFz1SfaTlVXh/V3sBxs5Uw+eyvtvvxnCMndkQrzlvU1hNDVWMxDqMj7kHUNUCH3QUS6isg8oBFL9PsmsF5VY7qrtI1f2R/2fwAMKq7F5Y8LWuWzj6qOAw4GfiYi303vVItRVMzcjEqzN8UtwHBgd2A1cG1pzWkZEemLFWD8uap+mN5XKfchxzVUzH1Q1S9VdXescPEEYHSJTap4XNCMgipjlyOauwL4uzEkFJaNpbOwIPLZWzH3RVXfDT9QTcDvScJZZXkNItIdE4K7VPV/w+aKug+5rqHS7gOAqq4HZgN7YeHcWAUlbeNX9of9/YG1RTa17HFBM14CRoYRRj2wTteHS2xTi4hIHxHpF9exCuALMdtPCM1OAP6vNBYWTD57HwaOD6Ps9gQ+SIXEyoqsPqUjsfsAdg3HhFFqO1MGVddD38vtwGuqel1qV8Xch3zXUCn3QUS2EZEBYb0XcCDWDzgbOCo0y74H8d4cBTwTvGgnTalHpZTLCxvJtQSLY19UansKtHkYNnKrHlgU7cZi608DS4GngIGltjVl8z1YKOhzrI/gpHz2YiPBbgr3ZAEwvtT2N3MN/xNsnI/9+AxJtb8oXMMbwMFlYP8+WDhxPjAvvKZU0n1o5hoq4j4AuwGvBjsXAr8M24dhQtsA3A/UhO09w/uGsH9Yqe9BOb489ZXjOI5TFXjI0XEcx6kKXNAcx3GcqsAFzXEcx6kKXNAcx3GcqsAFzXEcx6kKXNAcJw8hc/u5LbQ5QkTGtPPn7iQixxXQ7p4wd/LnInJse9rgOJWIC5rjbBlHYJnc25OdgBYFDdhJLdHuROC5drbBcSoOFzTHSSEiF4nIEhF5HhiV2n6KiLwU6lc9KCK9RWRv4DDgmlB7a3iuduH4o0VkYdj+XNjWVUSuCe3ni8h/hI+7EvhOOOdZOWy8S0QWA6NDcttJwKMicnLH/nUcp7zxidWOExCRPYCZwLeAbsArWA2qGSIySFXXhnaXA++q6n+LyEzgEVV9IOzL124BMFlVV4rIAFVdH6ojbKuql4tIDTAHOBoYitX0OrQZW48GdsRKicxQ1aM74E/iOBVFt5abOE6n4TvAQ6r6KYCIpPN57hoEagDQF6jLc4587eYAM0XkPiAmA54E7CYiMXdffyzH4KYCbB2HpanaDUt95jidHhc0xymMmcARqlovIicC+7amnapOE5FvAYcALwdvUIDTVTVDHEUk37kRkSnAFVjV5UOBbYBPROQAVd2vjdfmOFWB96E5TsJzwBEi0itUMZia2tcPWB1Klvwotf2jsK/ZdiIyXFVfVNVfAmuwUiB1wKmhLSKyS6iakH3Or1DVx4A9gIWq+q9YUuqxLmaO4x6a43yFqr4iIvdiIbxGrKxQ5BdYReQ1YRkFZxbwexE5Ayvrka/dNSIyEvPKng6fMR8b0fhKKIeyBhs1OR/4UkTqgZmqen2WqWOB+lDqqLtmFed0nM6KDwpxHMdxqgIPOTqO4zhVgQua4ziOUxW4oDmO4zhVgQua4ziOUxW4oDmO4zhVgQua4ziOUxW4oDmO4zhVwf8D564nql3B6ckAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# now process the output\n", "thickness = []\n", "sld = []\n", "for objective in obj_out:\n", " model = objective.model\n", " thickness.append(model.structure[1].thick.value)\n", " sld.append(model.structure[1].sld.real.value)\n", "\n", "fig, ax1 = plt.subplots()\n", "ax1.set_xlabel('dataset #')\n", "ax1.set_ylabel('thickness /$\\AA$')\n", "l_0 = ax1.plot(thickness, 'b', label='thickness')\n", "ax2 = ax1.twinx()\n", "l_1 = ax2.plot(sld, 'r', label='sld')\n", "ax2.set_ylabel('SLD /$10^{-6}\\AA^{-2}$');\n", "\n", "# setup the legend\n", "lines = l_0 + l_1\n", "labels = [line.get_label() for line in lines]\n", "ax1.legend(lines, labels, loc='right');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we see an induction period during which the film thickness and SLD remains constant (even a possible loss of solvent?). After this period the the thickness grows rapidly before with the rate than slowing down. The SLD has a continuous decrease, but the decrease is not significantly faster during the period in which the thickness of the film grows rapidly." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }