From 14e415dfc51dd80d6c9f6ad64dd6492380b8603b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 23 Jun 2020 08:13:52 +0300 Subject: [PATCH] Add the size vs. lMax script from triton. Former-commit-id: b2670f5fb1c54242dd670a4e71eb018a85e2ea6c --- .../size_vs_lMax/AuSphere/_mode_r_dep.ipynb | 334 ++++++++++++++++++ .../AuSphere/parallel_generate.sh | 19 + examples/size_vs_lMax/AuSphere/serial.sh | 13 + examples/size_vs_lMax/AuSphere/serial1.sh | 11 + examples/size_vs_lMax/_mode_r_dep.ipynb | 300 ++++++++++++++++ 5 files changed, 677 insertions(+) create mode 100644 examples/size_vs_lMax/AuSphere/_mode_r_dep.ipynb create mode 100644 examples/size_vs_lMax/AuSphere/parallel_generate.sh create mode 100644 examples/size_vs_lMax/AuSphere/serial.sh create mode 100644 examples/size_vs_lMax/AuSphere/serial1.sh create mode 100644 examples/size_vs_lMax/_mode_r_dep.ipynb diff --git a/examples/size_vs_lMax/AuSphere/_mode_r_dep.ipynb b/examples/size_vs_lMax/AuSphere/_mode_r_dep.ipynb new file mode 100644 index 0000000..279b8b2 --- /dev/null +++ b/examples/size_vs_lMax/AuSphere/_mode_r_dep.ipynb @@ -0,0 +1,334 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import re\n", + "import numpy as np\n", + "import matplotlib\n", + "from matplotlib import pyplot as plt\n", + "from scipy.constants import hbar, e as eV, c\n", + "eh = eV/hbar\n", + "import glob\n", + "def ri(z): return (z.real, z.imag)\n", + "#m = re.compile(r\"([^_]+)_r([0-9.]+)nm_\")\n", + "#removek = re.compile(r\"(k\\([^)]+\\)um-1_)\")\n", + "remover = re.compile(r\"r[0-9.]+nm_\")\n", + "\n", + "\n", + "markerdict = {\n", + " 4: \"3\",\n", + " -4: \"4\",\n", + " 3: \"^\",\n", + " -3: \"v\",\n", + " -2: 'x',\n", + " 2: '+',\n", + " 1: 's',\n", + " -1: 'd',\n", + "}\n", + "\n", + "prop_cycle = plt.rcParams['axes.prop_cycle']\n", + "colors = prop_cycle.by_key()['color']\n", + "colordict = {i: colors[(i+1)] for i in range(-4,8)}\n", + "\n", + "def markerfun(b):\n", + " if b in markerdict.keys():\n", + " return markerdict[b]\n", + " else: return 'X'\n", + "\n", + "def colorfun(b):\n", + " if (b+1) in colordict.keys():\n", + " return colordict[b+1]\n", + " else: return colordict[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L1_cn410.npz 203\n", + "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L2_cn410.npz 251\n", + "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L3_cn410.npz 251\n", + "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L4_cn410.npz 251\n", + "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L5_cn410.npz 251\n" + ] + } + ], + "source": [ + "allfiles=glob.glob('*sph*k(0_0)*.npz')\n", + "allgraphs=dict()\n", + "for f in allfiles:\n", + " base = remover.sub('', f)\n", + " if base in allgraphs.keys():\n", + " allgraphs[base] += 1\n", + " else:\n", + " allgraphs[base] = 1\n", + "for k in sorted(allgraphs.keys()):\n", + " print(k, allgraphs[k])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['sph_r103nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L4_cn410.npz', 'sph_r238nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L5_cn410.npz', 'sph_r257nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L5_cn410.npz', 'sph_r147nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L5_cn410.npz'] collected_XXXsph_rXXX_p580nmx580nm_mAu_n1.52_bX2_k(0_0)um-1_LX_cnXXX.pdf\n" + ] + } + ], + "source": [ + "projectors = dict()\n", + "projectors_list = dict()\n", + "lMaxes = [lMax for lMax in range(1,6)]\n", + "for lMax in lMaxes:\n", + " proj = np.load('projectors_D4h_lMax%d.npz' % lMax)\n", + " irlabels = sorted(proj.keys())\n", + " proj = {f: proj[f] for f in irlabels}\n", + " proj_list = [proj[irlabels[i]] for i in range(len(proj))]\n", + " projectors[lMax] = proj\n", + " projectors_list[lMax] = proj_list\n", + "globpattern = '*sph_r*_p580nmx580nm_mAu_n1.52_b?2_k(0_0)um-1_L?_cn???.npz'\n", + "filenames=glob.glob(globpattern)\n", + "plotfilename = 'collected_' + globpattern.replace('*', 'XXX').replace('?', 'X').replace('npz','pdf')\n", + "print(filenames[:4], plotfilename)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "#projectors\n", + "#glob.glob('cyl_r100nm*L3*3100.npz')\n", + "#glob.glob('sph_r100*m5*.npz')\n", + "#dat['meta'][()],list(dat.keys())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "inpure result detected [1. 0.99999999 1. 0.97991334 0.99999996 0.9999989\n", + " 0.99999983 0.99999966 0.99999322 0.99999721 0.99999653] [3.28735741e-04 2.66532534e-05 2.47011478e-05 1.45012420e-01\n", + " 2.44785416e-04 7.05405359e-04 1.60203586e-03 1.71245137e-03\n", + " 1.03244480e-02 9.18732728e-03 1.18651583e-02]\n", + "inpure result detected [1. 1. 0.99999998 0.99999999 0.99999996 0.96608887\n", + " 0.99999852 0.99999397 0.99998951 0.99999912 0.99982435] [2.66223026e-04 2.12357147e-05 3.54211968e-05 1.06651057e-04\n", + " 2.79595790e-04 2.41939163e-01 2.17645058e-03 3.41541473e-03\n", + " 1.14507609e-02 1.49639498e-02 2.33483138e-02]\n", + "inpure result detected [1. 1. 0.92521572 1. 0.99999627 0.99990293\n", + " 0.99946049] [1.59712906e-05 3.60193407e-05 2.48341492e-01 1.21848930e-03\n", + " 3.81805601e-03 2.42649228e-02 2.99534246e-02]\n", + "inpure result detected [1. 1. 0.99999998 0.99999961 0.93267685 0.99999964\n", + " 0.99999822 0.99921774 0.99995547 0.99997301] [5.22490396e-04 3.01556792e-05 4.88795563e-05 6.29703960e-04\n", + " 2.34414238e-01 3.72766210e-03 4.72444059e-03 7.62106094e-02\n", + " 6.32796684e-02 5.63231562e-02]\n" + ] + } + ], + "source": [ + "plotdata = {}\n", + "for file in filenames:\n", + " dat = np.load(file, allow_pickle=True)\n", + " kx = dat['meta'][()]['k'][0]\n", + " radius = dat['meta'][()]['radius']\n", + " b = dat['meta'][()]['band_index']\n", + " eigvals = dat['eigval']\n", + " lMax = dat['meta'][()]['lMax']\n", + " residuals = dat['residuals']\n", + " ef =dat['empty_freqs']\n", + " eigvecs = dat['eigvec']\n", + " irweights = []\n", + " #for proj in projectors_list[lMax]:\n", + " # try:\n", + " # irweights.append(np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)))\n", + " # except ValueError as err:\n", + " # print(proj, len(proj))\n", + " # raise err\n", + " irweights = np.array(irweights)\n", + " #print(irweights)\n", + " irweights = np.array([np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)) for proj in projectors_list[lMax]]).T\n", + " irclass = np.argmax(irweights, axis=-1)\n", + " purities = np.amax(irweights, axis=-1)\n", + " if (np.any(purities < 0.98)):\n", + " print(\"inpure result detected\", purities, residuals)\n", + " #print(purities)\n", + " \n", + " #for i in range(len(residuals)): \n", + " # if residuals[i] < 0.01:\n", + " # vec = eigvecs[i]\n", + " # for irlabel, proj in projectors.items():\n", + " # print(irlabel, np.linalg.norm(np.dot(proj, vec))) #maybe some conj() here?\n", + " # print('--->', irlabels[irclass[i]])\n", + "\n", + " \n", + " plotdata[(lMax,radius)] = (eigvals, residuals, b, ef, irclass,)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "collected_XXXsph_rXXX_p580nmx580nm_mAu_n1.52_bX2_k(0_0)um-1_LX_cnXXX.pdf\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+8AAAGFCAYAAACMisIwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeXxU5b3/309WAkQCBlxQQLFQQ7UqLiVc6zKKFpdC29va1rYuXbzW0nv1F6ytiTWxrSZdY229tq63i61toSJoNQO2CmpFxQ2XIiQUVEAhECAh2/f3xzMnORlmOTOZs8zM8369zuvMnDkz88zkkzPP9/luSkQwGAwGg8FgMBgMBoPBEFwK/B6AwWAwGAwGg8FgMBgMhsQY491gMBgMBoPBYDAYDIaAY4x3g8FgMBgMBoPBYDAYAo4x3g0Gg8FgMBgMBoPBYAg4xng3GAwGg8FgMBgMBoMh4Bjj3WAwGAwGg8FgMBgMhoBjjHeDwWAwGAwGg8FgMBgCjjHesxyl1L+VUsf7PQ6DIRWMbg3ZhtGsIdswmjVkG0azBkNyjPGexSilxgIHA2v9HguAUuoqpdRqpdQ+pdQ9fo/HEEyCpFulVKlS6k6lVJtSqkMptUYp9TG/x2UIFkHSLIBS6jdKqXeUUruUUm8qpb7s95gMwSJomrVQSn1AKdWllPqN32MxBIugaVYp9XhEq7sj2xt+j8lgAGO8ZzvHAG+KyD6/BxLhbeAm4C6/B2IINEHSbRHwb+A0YAxwPfBHpdQUH8dkCB5B0izAD4ApInIAcCFwk1Jqps9jMgSLoGnW4jbgWb8HYQgkQdTsVSIyOrJN93swBgMY4z3bORZ42e9BWIjIX0RkMfC+32MxBJrA6FZE9ojId0WkVUT6ReQhYANgDCGDncBoFkBEXrVNcCWyTfVxSIbgESjNAiilLgLagbDfYzEEksBp1mAIIsZ4z26OAV5y68WVUg8ppdrjbA+59b6GnCewulVKHQRMA151a3yGrCRwmlVK/UIptRd4HXgHWObW+AxZSaA0q5Q6AKgHrnZrTIasJ1CajfADpdR7SqmVSqnT3RqbwZAKRX4PwDAsjgVcM6JF5Hy3XtuQ1wRSt0qpYuC3wL0i8npmR2XIcgKnWRG5Uin1DWAWcDoQpFBTg/8ETbMNwJ0iskkp5caQDNlP0DR7LTr/vhu4CFiilDpORN7K+OAMhhQwnvcsRelfvw/h4iqlwZBpgqpbpVQB8H/oH+mrfB6OIUAEVbMAItInIk8ChwH/5fd4DMEgaJpVSh0HnAX8xO+xGIJJ0DQLICLPiEiHiOwTkXuBlcBcv8dlMBjjPXs5AugTkTal1OlKqb8ppRYppV5USn0mcv+fSqkDlVIHKaVWKKWeUEr9SSlVqJSao5S6C0Apda9S6szoN1BKPWyrshm9Pez5JzbkAoHTbWTScCdwEPBJEelx+0swZBWB02wMijA574ZBgqbZ04EpwEal1LvA/wM+qZR63t2vwZBFBE2zsRDAhI0YfMcY79nLscArtvsFIjIfXcn1IhE5Bx0C/HFgB3C2iJwKbAbOFJFHgS6l1G3AThFZHv0GIvIxW5XN6G2/dlpKqSKl1AigEChUSo1QSpnUDIOdwOkW+CVwNHCBiHRm9NMacoFAaVYpNUEpdZFSanRk0noO8FlMETDDIIHSLHAHenHpuMh2O7AUOCejn9qQzQRKs0qpCqXUOdY8Vin1eeCjwCNufHiDIRWMYZW9RBf2sG6/bbu9GZgMHAj8UukemocC1mr3rejCXJMzNKbrgRts9y8GbgS+m6HXN2Q/gdKtUmoy8DV0vvC7ajAX82si8tvhvr4hJwiUZtHen/9CG0AFQBvw3yLyYAZe25AbBEqzIrIX2GvdV0rtBrpEZNtwX9uQMwRKs0AxuvXxB4E+dGHQeSLyZgZe22AYFsZ4z1JEpCH6UJzbCvgc8JCI/FopdSs6UlgBTcAVwC2Rc4Y7pu9iDHVDAoKmWxFpw4TBGRIQQM1uA04bzmsYcpugaTbG+L6bydczZD9B02zkOnvScF7DYHALEzafH4SBbyql/gqMjxy7AnhCRO4ANindf9VgCBJGt4Zsw2jWkG0YzRqyDaNZQ16jRCT5WQaDwWAwGAwGg8FgMBh8w3jeDQaDwWAwGAwGg8FgCDjGeDcYDAaDwWAwGAwGgyHgGOPdYDAYDAaDwWAwGAyGgGOMd4PBYDAYDAaDwWAwGAJO3rWKq6yslClTpvg9DEOW8txzz70nIuOTn5k5jGYNw8Fo1pCNeK1bo1nDcDHXWkO24YdmDcMn74z3KVOmsHr1ar+HYchSlFJtXr+n0axhOBjNGrIRr3VrNGsYLuZaa8g2/NCsYfh4FjavlLpLKbVVKfVKkvNOUkr1KqU+ZTv2JaXUvyLbl2zHZyqlXlZKrVNKNSullJufwWAwGAwGg8GQPmHeZjYPEeZtv4diMDjCaNYQJLzMeb8HODfRCUqpQuAW4FHbsXHADcApwMnADUqpsZGHfwl8BfhAZEv4+gaDwWAwGAwG/6jjeVaxlQU8bQwiQ1Zg1+wM/sIM/mJ0a/ANz4x3EfkHsD3Jad8A/gxstR07B3hMRLaLyA7gMeBcpdQhwAEi8rSICHAfMM+FoRsMBoPBYDAYMkA9J1DNBABWsZU6nh94rImXOYD7aOJlv4ZnMOyHXbNraWct7VzAY0zmDxzE7xjJvUazBs8ITM67UmoiMB84AzjJ9tBE4N+2+5sixyZGbkcfNxgMBoPBYDAEkBCHEuJQwrxNHc9TzwkDjzXwAh30ci3Pcj3PMZJCOunjIMq4i1MJcaiPIzfkK3bNLuBpNtBBJ31sZM/AOdfyLN9mNf0IBSgOZaTRrMEVAmO8Az8FrhWR/kynriulvgp8FWDSpEkZfW2DwQ2MZg3ZhtGsIdswmvUXyyCyU8vxfItn6Qe66aebfgA2sodzeQRQCKCAcZRSyQia+UheGUhGt/4R4lBe5RMDRvxueuiij210IUAvAkA/wkb2cBaPUIhCIYymmL30cbBZiDIMkyAZ7ycC90cM90pgrlKqF9gMnG477zDg8cjxw6KOb471wiJyB3AHwIknnigZHrfBkHGMZg3ZhtGsIdswmg0eNRzDCRzI5TzJO+xlJIXspCdiGAEM/pm20sVWujiLRyhGMYpi9tHLGRzKE7xLLcdTwzE+fRL3MLr1H8uItwjzNpfzJJvZM+B5twz5vsi+nR5AL0SdzSOMoZjdEVWPpsgY9gbHBMZ4F5EjrNtKqXuAh0RkcaRg3fdtRermANeJyHal1C6l1EeAZ4AvArd6PW6DwZCjNLVBQyvUToGayX6PxmAwGPKCEIfSyqcH7ltezvfoYjv7BjzvMOjp7EFopxuAZZGMSiuMGWA0xeymG4XiUEZxJ/9hDCRDxojWLOj6Ddezmj4Y8LxbC1HCoDEPQw37s3iECorpoJc+hEIUB+ZplIkhNp4Z70qp36M96JVKqU3oCvLFACJye7znRYz0BuDZyKF6EbEK312JrmJfBjwc2QwGg2H4NLRCR5/eG+PdYDDkOuEw1NVBfT2EQn6PZoBoL6fFoLdzNyrK8/4Im+hn0Li3DHsQ2tjNWTwy8DrGODJkhPB2qFsP9UdCaBw1HLNf5EeYt7mMJ9hCJ2UUDvG8W4Y9DDXs+5AhUSYAIyigHziEkWYhKg/xzHgXkc+mcO4lUffvAu6Kcd5q4EPDHpzBYPCH8BNQ1wj1CyF0qt+jGaSpDbr7oURpz7vBkA4BNYYMhpjU1cGqVTB/PixaFHjNxvJ2WtjDmEF73gcN+KFEG0cjKKA34vFsYGZOht4bXKBuPazapfehcTFPCXEobXwm5mOWZq10EcvzHouuSC0I+0JUAZhCeXmCl33eDQaDYSgLrodVq/U+GeEwzJ6t927T0Ar7BEoLjNfdkD4LFmhjaMECv0cyQPg9mL1K7w2GIdTXQ3k5dHRoQz6LsQz7Hi6lh0vZwcU0cpIjj1UX/fQi7KOfhTxLOfdRzN0UczeT+YPp722ITf2RUH2A3qeBpdl9XMIOvkAvl9LISZRRyCgKkz7fijSxQu8LuMtoNkcxxrvBYPCftW/CeRcnPsfyCnkxqTx1zNC9wZAj1L0Jq9r13mAYQigEtbXagJ83z+/RDCETi041HEMPlyGRrYVzmcQoSihIaBztppdeZMAwmsMjFHM3RdxNiTGODBahcbDyxAGve6Y0u5cvsZsv7afbZAacVf1+I3tYwNPpD8IQOIzxbjAY/KP5psHby8LQ9Iv459bXQ1UVtLe7631vaoNlkbIaT+x07328IBzmaPig38PIScJPwOwL9D4ezc1as+BNxEg8wtth9moIb6d+GlSNgvbeAHvfjW7dwYlmFy/WnveGBn81G4Ubi05WCPM+LhkwjizDqIj4LYstD2cfQo/dODr6EKPZTNP0CzjgqMRzg4BiaXb+c5m91lq67Ysy5otRFKFiKnct7RRyF6XcQxMvZ24wBl8wxrvBYPCP0Kkw15ZXef3NCc6NnLd2rbthyA2tg7ezPd99wQJGwii/h5GT1DXqlI85F8WfWHql2WRYuZjzXyb04nZQsHY3LFjr35DiEg7D/PlGt25gaXb+pfENeHvo/Pz5gTHg66dBdYXeu4llGPVwKS2cSxUVTGBEQmNeoY0jRpYYzWaahh9Dxx5Y2JB1Bnz9NCgv1HVv3Yx0sjTbHUkRecy2AGU35vuBbvpp4AX3BmPwBGO8GwwGf1n6Gygt0be7e5z9QG/Y4M6ksqkN9vVDqYLGqdmd7x4O6+/J4A71C/W+vz/xopOFW5p1Qv2Rtlnk+oFW2Rv2Bsz7HjHc6eign0hFJkPmqF8IJSXaGLrs6tjnhEK6WF3Act9DldoYqnvTO81aVe638LmExrxptO4itTadZpkBH6qE2qP0pXfeQR6+r20BqodLuYWTKKGAAqCEAmo53rvBGFzBGO8Gg8F/Gq4dvJ3oB7q5WU8qOzvdmVQ2tEK3QEmWF6oLh+GCC6CzEzFGkDuETnW26OS2Zp0QGqejSMoLYd54mmfom539mQ/pTBtLsx0dUFbGeljn95ByjtCpUBiZ9m3clDhiJIC5726FITsl2pi3iolVUEwJCnr6YpezN6RPzZXQWDt4v/YW/8aSBou3RDrOrvPvOlvDMezjEvq4jH1cYron5ADGeDcYDP4T/QMdz5Pp9qQyVwrVLVigjUWgmzj9kQzDx77oFG9SaXkyvajXkIjF2/QscvE2QpWwaKbuhNjRB5e95M+QhmDTLEccwU7o8HdAOcqNNYO3Ey2UBjD3vX4alBVozQYh5cMqJraDL7CPS+GlTSaZ2A1qroSSYn17X3dWed+HaPZVv0djyBWM8W4wGIJBzZXOPJluTiqtAnXZXKjOHi5fVsa/YaO/A8phnE4qg5D7Xn8kTC6FF3ZDUxuhysGHNnb57H2P0izNzT4OJsdx6skMYNu4UCUcUaZvBy7lw+AuN31r8HYWed9DlTA+Mq1ZZzRryBDGeDdkF7aqyYYcxIkn061JZVMbdPdrd2S2FqqzcoY7O/V3tGSJ8WC6jX1S6ST33S9C42B7r46VjxRlPKR08GFfve+W172sDJYsGVzsMLiDk0WngIbON8/QnszOfuPJzCuy2Ps+ukjvuwUuWG0MeMPwMca7IXsIb4f5L+uqyXXr/R6NwQ2cTirdCENuaIV9AqVZnO++YMFAzjCLFgXCCFpBmDOYzQqCEXqbcewRIyp+RepAtI2rnaIXp7r7oamNO49loOzWv/3yvtu97kccEQjNxiWXFo/ti04NP459TgBD50OVcMRIfXtDpzGE8gonmg0gzVV6wQn0olMg0pQMWY0x3g3ZQ916nThUXqhDQA1DaWqDA/6u99mME0+mG2HI2Z7vHlAjqJ46nmYV9QQj9NYVGq7VBrxIsEPnaybrxal9Ag2thCrh8BH6IcHddkYxiY4UCXq4vNVyb86a7L/O1lw52Kbz1I/EPieAofOgjaHAFVw0uI8TzQaQUCUsOTEAC6WGnCGvjfec9wjlEuHt0N4LVSNh0TE6BDQPSajZhtZIWdNWr4eVWZzmvltkogVXUxssi3jTsjXf3R56HCAj6ELmMZpyLiQ4obcZx4oY6e5x1s7Iz7ZxUYtUdx0LE4r1xHLGaI/HEsBIEYuY19r6I/WsqZ/sv84CPPG03i8LZ1XovFVw0Yse2tlEXsxpk2k2oPi+UGrIKfLaeE/oEcoVL2YuYIXLr90LFUWODfdc/CFLqNls9xzbsee+xwuPy2QLLvtEPBvz3QPqdQd4kMXspoMHWbz/gz5cZ8N/h9nn6H1GsfcjjlevIQht46KKMoYqtQEkwK83eegRCrBmIc61NjQOzo38/uTCddaJZgMYOg/+9dAOMkGa03pynQ163/eoNBtfF0oNOUVeG+911PMRqqmjfv8Hc8WLme00tekQxTTC5XMxXDehZnOhUrqFVRG5tMSbgkrWRHzuuOzMd6+rC2zocdCus3Xfh1X/1PuMki1FwKx+77ZFKqsasuBhQaUAaxYS6DbXrrPJNBvQ0HnIcA/t8BMw+wK9z1KCdK119TqbLX3frTSbSI2mWAul7xFmJbN5Lx0nUw5o1pAeeW28n0GIOuqpp25/72wueTGzlaY2WPiWDlEsIKVw+RWE2Uk7R1MV+4csS0mo2dopUKpgX39uRIzUXKn33T3ue4UCMiFPO1pk3jw9wa6tDZwHM6lmbQXUvKD+21A1Hdp3uuAVclKvwW9PZs1k/b03tA5853cdO7Sgkic9tAOs2YT4oFlw0ZNp12wsT6ZbBUIzQKy+72kbQ3WNsGq13uciHs9pXb3OZkvl+XnjI6Eh4wcOTbAtlF72ErxBHTtYxbPMN5o1OEdEXN+Au4CtwCtxHv848BKwBlgN/Efk+BmRY9bWBcyLPHYPsMH22HFOxjJz5kyxc7pUywhBTpfqIcel/HERwnpv8IfS5fpvQFiksTWlp8b9uw4TYLV48D8j6WhWJPd0W3K4CAfrrfG2/R9vaREpLxcBkeo0/9aNrSIly7XeUtRZpklLt0m+A6PZ/ak6RYQKvc84pZPc1+xwifGdt2wTKVkmwlK9b9nm4vs7+A681m1ea7bxtkHNlhweZwBV+u9VVeXCANKn6nGt2bKHtWaflGpZIsiTqf72t/xDZPKJImVTYv/fOsDva23QNFs9R2u2/DCRlky/rV2zpZMy/OIZovpZ/Z1XPztwqGWbiFqqNauWivytvUWWSXleadZsw9+88rzfA5yb4PEw8GEROQ64DPg1gIisEJHjIsfPBPYCj9qeV2M9LiJr0hlY3GJKPq2uGyI0temKyKD/DimGMudykayEn82HiBHXPEKQ3JOZiTDkhlbdgLXE3xZxaUWLWNW6Ozr0d1AfzCiToGnWVZLVawiCJzPGdx6qHOz73i0uF1SyF6oLqGbrqOdoqthJeyAiRlwl1SKhAeKSw3QOsdX3/WDmUUg5B6f62x86Fd7ZAp1d8aNmAk5SzXocmVf/bSgfDR27fUpT8pv6I2FSKbywe+A7D1XClw/TDwvwqadDFO5cxGiq6KE9Ne976FTYvkNrNuj5/4aM4onxLiL/AOI2RhWR3ZEVIIBRaE1H8yngYRHZm8mxxS2mZG+ps/Ct3PiBzibseVk3pd4WLmGRrCwn4Wezwr6XbfdMs1Zu2/yLXQqPSzapHE4YclObnoCXKN8L1dVTx2usZQwVnIHDMOIAV+u2EzTNNt+iQzrBJc06qdcA/rWNi5MmMqpw8LZrBZUCXqjO4gxCjKGC11i7f92UmsnaYtwnULvek/G4qllIvujU3KwXnCBQofOLt9gmjAreZTF9dPBuOr/9KtLMK8sWMCySalbQK3PXe6PZ0GlQW6MN+HnnufAGQe/7HhoHO3r1ypJtTvvq7sFTOvrgf9aEKKaC3azljVRrNDkpOmnIOQKT866Umq+Ueh1Yiva+R3MR8PuoY99TSr2klPqJUqo0wWt/VSm1Wim1etu2bUMeS7pSaWEK13mH3aBqnJqWNzRh4ZYsIJs06+rqOgydVMZaXR5OQaWGVj0BL81Sr3uAjKBs0mzoNKgYA2vfcEmzTuo1+Ekcz3HzjMFJwZ1uVJ4PWF/3RJqFJBEjErV3GU80m6iHtt8LTnGonwZVo6BqtO7/Pp369LyY4KzLSQBIpNuEmlVRew9YvFTPDRYvdeHFs6Hve4wop/ppOhXeYk/fMCJG7BEI4tHFyOA/XsXnA1OIk/Medd5HgZaoY4cA24DiqGMKKAXuBeqcjCM6r00kSZ7Q3BdECOu9wRsCnLdNAPKHRYKn2caf6by2xp+59Qa2/Lbyqfs/3tKiczGrqvRtpwTk/zutXHcrB7WsLOFnzgrNNrbqmgMl3tUdaHlc52RmPBfTwkm9hnQ0myniXGe/8qIIkZzMqkx/Nw41K+K9btPSrMe1Mloe1znvVae4pNvyqYk1G9C8d4uWbSLVK0Ue6k4z711Ef+6Sw3UedYp5xEG41hrNBow419mWbbpOg1WvIZ80a7bhb4HxvFuIDrE/UilVaTv8aWCRiPTYznsnor19wN3Ayem+Z96038oWYrQySoVc7O8eTdA06+rqOrjjFWpq06Ha4Pv/d8o1GgLmdXdCQs3WTNY1B7rFU+/7yr/pvSvYQzpjed/99mTGyYG1h3S+tTeD3vdc1KyVWpcrESPJQnADGjpvseBVWNUOf1z/hfS8mKB/a5TSKS9BjJpJQl5rNoh533Gus6FKuPEDg/Ua/nftMCJGslyzhtQJhPGulDpKKZ1spJQ6Ae1Nf992ymeJCplXSh0S2StgHvBKuu9/BiFWsDJ2nmmuFaYJMk1tcEAkmW/XaWmHMedif/eU8KkwjWttYSyeeFrvl4UT/0Bv2OBsYmmfvPic755SjYaAhR47JWHLOMi9wnWpFFRyqtlMEmfBpH6a/skDPcef/1yGDPgFC7Rmy8qMZofBvPNczCFOptlQCCoq9ILT/PnBM+Ajuj12wv+ln/cOg+HH1j6LyEvN2vu+B63YYIKFaXu9hkffDvHie6H08t4hqzVrSB1PjHel1O+Bp4DpSqlNSqnLlVJXKKWuiJzySeAVpdQa4DbgM5FwDpRSU4DDgWiT4LdKqZeBl4FK4CZXBu/DSmXecv16Xb1jGMVUcrW/ezQJFyh88mK6uroOQ1fYY/1ANzdrY7az01nuuzV5mTsuu/Lds6RIXSwS6jYXo5ySFVRKVbOZJk7V+aNG6dsKfUkeduX5LPS6WwRNs65HOSXr+15fr689HR2Byn0HuGSi1ux9b9bz7p40vZigv4OSYu3NDJon1wF5p9mAdUvYHIYHZ+s9EHfBJHqh9Pcbzqez5+D0IkZu+tbgd5CFmjWkhlfV5j8rIoeISLGIHCYid4rI7SJye+TxW0RkhuiWb7NE5Enbc1tFZKKI9Ee95pkicoyIfEhELhaR3dHvmwoJQ61zzSMURJratLEJwyqmklbF7iwkYQEw8MX77urqOiT/gU61bVxAjMWUNJvFRhAkSQ/wQbPgcqtDJ+kew211OBzi/A80V+kCYOOLoawA5h00jPfI0kgRi6Sa9Tgyz/Uop2hPZvSiUyikrz0esTkMf5qht81J7HDLk/nS+yF2dqdZvRsGf2v2dQe6cF08gtgyzvXIvAAVG3xqAWxdpfdA3Ots9ELpJ6feSFnxu+lFjFhRM909vn9+g/sEImw+CARtpTLvsHvbG1JvDWeRy/3d7SRsCQO+eN9dX12HoT/QsXK7nLaNC1CLuJQ0m4Whx3YSpgf4oFkYbHXoWsRIsnSP4bQ6HC4J8jErimBrj87HvOFfwwidz+JIEXCg2VzLIYbki04e5r4/Vwfta/X2VBJHf/00mDxi0Pu+aVc109ONwLM+d1CrmCcgacu4XIzMC2Dl+Y4NkQWnBHWcmqt0t4TxxfDAunr2duanZg2pYYz3CAmLfAyzgJohCXave6kaVghzLvd3jyZpOzyPI0Y8WV1P1hbFaUhnQFrEQQqazXKvOwRPs+BBxEiygkr2Vode5xEnmMjXT9OXY9AG/IK1aby+0awruK5ZSLzo5GGxxZn1UFimbw8YQ3EIVcKoQu19f/n9EDt33Mob1KUXOu+0xkpASTqnzUXve0D+ZrOaoagc+jrhsfmw+YTJ+jtvaN3v+w5VQkWxXih99r0QX3nibzzTk5+aNTjHGO9OqIn/j2fIABnyukP293dPhaSFaTyOGPFkdR10blf5qKG5mQODcBjSGaBUGMeazXKvOwRPs+BRPmayMORFiwYNeA9y34fkZMb5XwhVwtSRg/c3pFN53mjWFTyJcrIvOvkYhjsxBHOWRBlDieyayIKTAN99fTqPvzcivdD5gHx+V8hV73tA/mYTQ3D2Ir3o1NsRiRhpaNUFRGJ83/bc9z19xazZaTRrSIwx3iMkrVCe4B/PMAwy6HVfQZh66qijPqfz3e0k1K0PESOeed93rdP7WDgJ6QxIKoxjzeaAB9MibzXbWKsXnewTLAuPc9+fq9M5mY/Nh80rIrPGGP8LzTNg0ojBdkYpFa4zmnVvPEHQrIeh85YxVGA3huLQXKXrNAB09Y3iztd+xNadX0j9TQMYhp0KSee0PujWk7o4iTTrIRNDUB7xI3RsgM0XfTBufQx77nuflPLr136cvmYbawdz3433PXfxu9G819vMmTMlFsulRU6XalkuLTEfl8ZWkdLlIiXL9W1DZih/XISw3ob5vZ4u1TJCkNOlOkOD2x9gtQREsyIOdVv+uKearZ4jQoXe+0Z1tQiIlJeLtER9NwH6X3ak2ZYW/TnifZ4kGM0mJ/CazTCbWkTuLhf5FSJ3l/bJprI1cb/vqr+LsFRELRVpXJfCm1RV6c9TVpbW5/Fat0azaWD9jauqPHm7B6q0Zu8q0xqOR+M6kbKHRYof3issFfngky+l94blU0U4WO8dEKRrrdGs/2xq0Vr9FVq7A/Pd8sf3O7dlm8iksIha2js8zYqkpFs/NGu24W/G8x4haWicT8WUch5r9bdx6rC97vnQIi5lfIgY8SQfMxmJct8DlO/uqFhdlhf8ShmjWdfziC1PZlE59O4roH18Rfy0MBncNaxzGDqfQ153cDA/yFfNWmzY4EmthlnNOhS5rxMevSB++K9HTigAACAASURBVPziLTpSpL+/hAkj1jOudGN6OcS1V2etFzOImvUkYiRA2L3v7a/B24dG2nbESNcLVcLoQhAKKVadnD/lJ+lpFkzhujzAGO82koYZBShPNutpaoMDIlfvXacN24jKlxZx0QRRs57kYybDnvtun1g2tekiPaX+V5kHB8XqcswIAqPZuMTTrEtYBnxFFUz7d2vcibw9dL6jD+Y/l8SAz/LWcPFIGjrvccu4QGi2uVn/jTs7PanVYDeG+jrjh8/XT9M+gT4K2dY1meljV/BMT/61jAuaZj2rixMgZjWjLS2BZW9MYzNj46frRTKYDhzxDqdNvDs9zcJg4Tprb8g5jPFuI2nv7IDkyeYEtev1TLB2ffJzHZAvLeKiCaJmA7O6HmtiWbt+sMaCz153R9EiOVDwKxqj2QT4YAyVVMDbMkY72BN6hAYN+ISV53M0UiThb4wPLeMCoVmr2GJVFbS3e+p9h/jV50OVsGimNuCFQn73rxv50Yu3pveGtVcPGvBZ5n0PmmYhYBEjHjAxBCfdTMTaKuAxjmFz15iYCybNVVA1Gvb2TuCXr9zKJf9Yml6LziyOGDE4wxjvNpL2zjYt4zJDU5v+wYCBkMzhkk8t4uw40qzHLWECs7oeqwiYpbcM6W44JI0WyUGvOxjNJh6I98bQzHo4tGCndvok8QgN/NvE+//JUc2Cg98YjyNGAqVZ8KRtHDivPh+qhNqjYHRRO6HD7mZN+5HpGUJWe9LunqzzvgdNsxCQiBGP+XANfOzRSPV5iniuZ0rMBZNQpd7v6hnN5j1Hs/D48/jZ22n8BmR5xIghOcZ4jyJh26aayRkJ8c577Betm4bXGs4in1rERZNUsz7UagjM6vrixdoL2NAAXw1rI6REZUx3wyGpZnPQ625hNJsAH4yh4pun0F+k6OvoZ+dXY3uEygv17bICQMUJnc9XzUJ+R4xYeJT77rT6/OItUHfieazfdTw7uiuSp3zEI0tziI1mg4OV8jGGPcxkA51Hx1kwiSyMfv4DdRw9bhWnHV6XV5o1OMMY7wZvaWrTOVYlathF6izysUVcSuTz6np9/WD/7F9dBfueC0ShuqSazWEPJjgoppSvOcTReGQMUTOZvv4C3qWCx39VsZ8n0wpDrhoFfQJrd8OCV6New2g2f3OIPU73AG0MHWBvxRXjT1I/Df7yVj2nHvIHRhe1J0/5iEeW5hA70my+Rjn5wKxm+DirmcgOip6NvWDSPEMvlP7mX/W8tr2a+96sT61Fp0WWatbgDGO8R5G0mJIhfZraYOFbGa/0ne9/s6SfP59X160w5LIy4HXg7kAUnEz4N8vRgl/RJPwO8jWH2KK5WWu2s9MT7zvA3sunsIORfIw1tF+2/0Q+VAmowZIRGzqjvJg57HW3CJpmISARI7FSlDwgWfX5UCV8YlyIl9/7OA0nn8exB4bTS5nK4hzipJrN5ygnj5kYgtcnTKGNA/lz4SkJ6zXs6zqFHzz/IG/uOJ15B6XxZlmsWUNyjPEehaPwa6tSukcrlTnD9bbidBmsG5Cvxeoskmo231fX7VW82Q0r3vd1OJBEszla8CuapLrN1xxi8LzyPMCYOyZTVfYuJfRx9LbW2CdFDB+FbsU1EIac4153i6BpFgIUMWJPUfIiWgRn1ecXb4FPTK3jqIpVXD69Ln7KRyKyOIfY0fzARDl5xoG/m8yWogo+0fsM738u9vcdqgRVtIEd3QfS1V/ovEWnnSzWrCE5xniPImmYEfjSHzPraWobdNmUqoyGLedrsTqLpJo1q+swuxm4CSiG8jv8Hk18zeaJEeSIfI4YAV+87wU3ToHyQjounsKDs/f3ZDbP0KHzh4/Que8dffDX3+VHpIgj8lmz9fVarx0dnukVhlafjzmsSOj8a9urue2VH8ZO+XBC7dVQPkrvc4l8j3LymIkhOL6vlRL6+ODW1pjed4BDSns4uOxfFKp9Ot3DaNZgwxjvMUgahmyqzqeOvSVcQ2aLheVzsToLR5rN59X1+4uA3wFvQeUrfo8mvmbzIPTYIojX2cB63z3m1b+OYOuq/T2ZVjXkjV36T1NeCNfeXpMXkSJgNJt4IN5HiwBMfL6NObzIhMlduqd29LAq4aYjQtz+6u3s6D5EDy865SOHcZRW6LFuA6NZn+g9SbfmfIcxcbslfGfKCXznhEuYOOpfQH5p1pAcz4x3pdRdSqmtSqmYM2el1MeVUi8ppdYopVYrpf7D9lhf5PgapdSDtuNHKKWeUUqtU0r9QSlVkomxJjUGTdX51LC3hivJrNfdFKvTONKsD6vr1Sfrve+cOgb4KYy6St/3aGIZi7iazTOvuyPN1k7RevUwRSlQum1u1m3jwBvNRqLK3ts2AohTCCzSNm5bj/a8X/GfN9I+8gCjWTCa9aFwHQ2tTOx8nwu3PwMQM2IkVAlXT93KJ478McWqj87+NDyZDT+Gjj1ZF4IcVOdGYDTrA2WvRVpzouJ2SwhVwl/+dStdvSMpRGs25cJ1WapZQ3K89LzfA5yb4PEw8GEROQ64DPi17bFOETkusl1oO34L8BMROQrYAVye4TEbMoELreEs8r1YnYWjdA8fVtfrv61X1n0NjWtqg2XbgULY+zHPWnDFI6Zm86RInR2TouSAUAgqKrRm589334CPROjMLNrA6JKumH20rbZxAhT3dnPWC49R+8WbjGYt8lmzVoHQqipob/dmwal2is7hGFfEcwt62LoKHp4DLzYNPW3lDugX6LcOqFTf5+rBHOIsKgAWVM0GZn7gB5E6RM+j58PxuiWMLtnBx4/4KWVFnUwaAe29KXrfTdG6nMUz411E/gFsT/D4bhGx6oCOIklNUKWUAs4E/hQ5dC9kpmKZI4PQFK1zhgut4ezke7E6O0FcyKj7Pqz6p8+hcfYJyfiH9d7DsM5oYmo2T4rURZNUsz4UAAuEZu14mUscidCZ2LOdTxc+Q2Gkj7a9krdVDfnY/nb+fNMn+OaSW7l4xW8Jf9hoFvAlRSlQmrWuXV4tktZMhuNHQ9s+ZrJez2r74dlvDTWIPj+tjkc2fpU+KaREAZKiIVRz5aAhlGWeTEea9bioLQRMt14SqUN0MusopC9ut4SvHV3HE+98ht29o9nWHadFZ8L3MUXrcpVA5bwrpeYrpV4HlqK97xYjIqH0TyulrBnvgUC7iPRG7m8CJsZ53a9Gnr9627ZtScfhKMwon1fXU6GhNeOt4ezkarG6VDULDnTrg2Z9L0zT1KYnJKWRxaPfzRsM6/TCkxmD/TSbI+Hyrmg2nwuAWXidSxxZKCk4Y0zcSt6hSrj/v0/jgmeXsmXMeL79xe+l14vYZ1zRrCkANohXi6TzxkN5IRMvGcFJNzNgwD9ns1VPKa6nUI0EoFDB2j22bglOCUgBsIzPaU1RW++pncLE8l2UT+gDiBnldEpxPWce9AKjC3spL9LHUs59N9733EREPNuAKcArDs77KNBiuz8xsj8SaAWmApXAOts5hzt57ZkzZ0pGaGwVKVshMulJkZb3M/OauUhjq0j543rvAsulRU6XalkuLa68fjTAavHwf0YyrdmS5SKly137e8Sieo4IFXrvOeWPixDWe4uWFpGyMhEQqaryfEhDNGsfS1mZvp9hsl6zpcu1bvNFs7HwUrO2/5lNLSJ3lon8Cr3f1BIZS1WVbB53iLx8eJWc/92HpPwRkcZ1mR2G17rNmGZFROa+oL/DuS9k7jWTECjNen2NrX5Wf9/Vz4qI1ulfqyN6tQ9rm0jV4yKTwyIly0RYqu87pvE2kfKpeh+DrL7W5rtm/aCxVTaVrZEHJuyTO0v0dfaBqH+Xk1buEJaKHPBoR3qaFdGa5WC9j8IPzZpt+FugPO8WokPsj1RKVUbub47s1wOPA8cD7wMVSqnIehSHAZszNYYVhDmD2fFzhKxQrY37oG597HPynaY2vYpbO8UVr7spVjcUR5r12CMEPq+uxwq7tnsyX3sNmpr2f55L7KdZq7o8ZLXXPV0cadZ4hLzVrC2EduLzbZyzBArKoN/yvi9YAGvXcuj2d+gpHUH4pHN1QE86vYizkKSaBRMx4nW0SP2RUH2A3qPbcV24Uu/tfLgyzF7W09Y1eCwlT2YWFwBLqtt816wfRIotfqrzKQ44Sh+Kzn///LQ6Plixkv845B56rGTidOo1BCBixJA5khrvSqnblFKz3R6IUuqoSB47SqkTgFLgfaXUWKVUaeR4JTAbWBtZMVoBfCryEl8C/pqp8TjKH54xSn+DM0Zl6m1zh6Y2WPiWq2HaQczx9hNH34cPOcS+toyLNyFpboaCAhCBb33Ls/D5IX8je7h8HrSGi4XRbAp4pdmoBZOJITggYoft+lcff33tLjZzJpSVceuVjXRSiEJf6lMOQ47DCsKoo/ng8F8p8zhuveVxDnHg2m81N+vrWmen+7nvoXHacK9bD+GhpZU2hwcr0L9BHf1qLwCHlOo6dylVnk9QtC7ImoVg1moInGa9xvbbNqtZL5JGpyh9vvLj/PCkhfxzy8UIWrOQHwulhvg48by/CfxQKdWqlGpUSh2fzhsppX4PPAVMV0ptUkpdrpS6Qil1ReSUTwKvKKXWALcBn4kY6EcDq5VSL6KN9ZtFZG3kOdcCVyul1qFz4O9MZ2yxcJT3fv9WXbr0/q2ZetvcwW6wu1DdfAVhdtLO0VQFrgWKXzjSrE+r6760hLHnu0drMBSCm28GpaC/37Pq81axuq+HZwytLr9kSd553cFoNiUszRYUaM26WbMhqjPFrGYoKofCng62ySmsVt+HJUt45gSt2WKl5/0dfbBgbfyXdcIKwnya+TCSQK6KO9KsTxEjgdKt1973uvWwatd+kZBPLYCtq/R+OvVcXXUH00fvZlQhTLCaCzv1ZMYpWhd0zUIwazVAwDTrNbbfNvsiqd37XkmI84pXsnBKBeWF+rKccuG6LI4YMcQmqfEuIj8TkVnAaehQ9buUUq8rpW5QSk1z+kYi8lkROUREikXkMBG5U0RuF5HbI4/fIiIzRLeDmyUiT0aOrxKRY0Tkw5H9nbbXXC8iJ4vIUSLynyKyL+VvYDh43Hora3C5wjzoVeTXWMsYKkzIfATH7eJ8WF33pSXM9ev15FmIrcGaGjj6aH3bo8JKD7KYE8MdnH/Br/Oyunw0QWxxCAFuY1RTA48+6l31+QgTCfPZcR9mFt9gAis5Sn7DZkI0z9B/mm7RRcAANuxN3ytkGUG76YD+wa5eQcKRZsHMD2Bo33e3C4TWHwmTS+GF3TF/23a9Bf+4DI4peo4TDnyGtXtgazdMGqEfd6zZGCHI17BAa1aCqVnHGM16S1SEzqxmKIzyvr9HmJXM5oEt7XT0wbYefTyldA9TtC73SCdRHp1z/gLQ53fSfqqb0+Iep0u1jBDkdKl2dL4hQmOrCGHZr0hYhvG6UJ0FAS9I40i3sYq4uYwvhWlKl+vPWbo8/jktLSLl5SIgUlAg0tjo6pCWS4u8fHKkkJNHxZxyQrMuF76MRaCLKVVVaf0o5Y5mo68R1vuBbOMEWUrLQAG7xnUi5Y+IfOVFvWepSPXK9N72eKmSEYKMlTJhDG+I0WxKBFKzPhVbtNjUInJ3uS4G9ucnq2WJIA91V0vZw1qrA0XA/p7C+9gK1y2XFhkrZTJCEPUhOiWg11qj2YASpdkHqoYWCH1StGZ/uu0bA9fXgqV6X/6ILsLo7H1iF63zY35gtuFvjgvWKaWKlFIXKKV+CzwMvAF8IvPLCcHAUWicYX9cDpc3JMaRbn1YXfe8AJjleSlR0HBk/PNCIe35tkKRXcwlXkGYN5ou4+jnItWS8jTPPZqgtuYMdDElt/Pf7R6h85p0kTyAsjIen/x33iZEf6S10R/e1H+a+9+Biw7Rl5Z5B6X+lisI04quATGFI2AnHZn7QJnFaDYFvAyfj1EfY2IIzl6kPZq76+rpfqqaqS/Xc8TIoU9NKWIkEoa84h838Wnm00knoylH/s3GzHyQzGM0G1DipChZ19cDn69nLNV8vvLj1B6lT738MJ373tGXYr0GU7QuZ3BSsO5spdRd6D7qX0H3YJ8qIheJSMYKxAUNx6FxTW1wwN89C0EONB6Ey1uYYnXZhecFwBpade5eaUFyHXqU/17fvoDPNWyksE+04ZWnee5pYQqADcVtzdpztpfdpRcJIpqddedoKqqgoAR6O+C82/WE0jLg0608fw0L6KSTMsr4EcFe1ApyilIgNetV8bo49TEmhmDOEuh9KsSuK2/l6Usn8IP256kaDQfbCtc5Lrh46kcAuObne9hNByMo448sCvSCk9FsdmBfbOrtgKf+I8SBz9fzBnW80BUeuM6mXK/BkFM48bxfB6wCjhaRC0XkdyKyx+VxBQJHBqIPK5WBpXa9NpgUrhruYCIj4hFUzXpalCZRobp42PPf3WjF1dRE3X++xu++AL3FShtexnAHHGrWFADbH7tm162DGTMy69W8qBdoBS4GKgc0OzEEn3qVgdZGhy+BX74IVaNgXHF6hetWEGaDzeueDXVMkurWFAAbxCvve4KoMqsgWH97Bb0vHUP7dRUgsLELxpcMLkDVvengfZ54mhVn7qN1vLbVjzCaHd64gqhZr4gxH5sYgvLIv0tfJ7yp6tjBKj4/rW5Ap6MK9TUXHC44maJ1OYWTgnVnisivgXal1MVKqToApdQkpdTJro/QRxyHIJcVwLii/VqU5BVNbfqCD7pAmIuY/u7xcaxZH7yYnhUAa2jVRl6JA6+7HXso8rXXZs4Yamqi97qF7B4pfPb3iqKHH9OGlwFIYSHOh5ZxgccqCNbdDWvXwgUXZEyz3DkHGA/MAnXvfpqd1YyeQQh0Xge9exnaPzuFMORrWEBXlnjdLRzp1mh2EC+K19VM1r9vDa0xf9tmNYPaOQGA3s0T6H5FH9/WnWLKR+3VXHPrbjpHCiNyTbM+zA/ymjjRDlbxOoD3z1lKSfgbnFJcPxA6f8lhgIpUnneyUGrC5nMKxznvwC+AWcBnI/c70C3dchZHYUY1k+H40dC2b78WJXnF9bbPflOCHOMMYELm4+NYsz54Meu+D6v+6UFoXLoTZnsrLpHMGENNTbBwIUV9cNoKmPfyJONxj8JxipIPLeM802y6WDUbyiKzvM7OjGmW/n7gN6C64Mul+502MQQn3YyOtOqH8Y8OPmaFITvxYkbnuufUgqzR7CB2rbrZKSFBZNnEEMx5YLTOKd45mm8shDHtWquppHysqJnO+UtHsuWACdQ1HZxbmvVhfhBYzXpBnGgHK9WjoAz6tlWwa0EzlYRYvGVQp7t79bmOFkprrtSGe8OPTcX5HCAV4/0UEfk60AUgIjuAksRPyX4cGYr1R0L1AXqfr1h5N6XK9ZB5q1f2hcxz9X2yFUea9SG3zZPQuKY2WBaJgElnwmxvxQXpG0NNTTBypPbgo4NRflRbwoUTv576mPKAIGs28MWUQiFdP8Gu2TlzUk/9CIdh8mRtuFvMPRFGj4IPHBXzKR+ugYpI5P5Fv4TrroYjX4AjyqBqNLT3JJ5UWq3hrIJf2eLBhBQ067EXM9CatYfPu5GeBEkXb0tDYSZtnM2Ij4UZuRcu/SGM3aM9706LgF3DAv5fQykHdBTwjQaXQw0ziGPHh8cRI4HWrBfEuU7Ye793TQ6zvH02N8wID4TOb+tOsV6DCZ3PGVIx3nuUUoVEgqKVUuMhy3taOsBRmFFoHKw8Ue/zkaY2rYrSJJW9M8SDLGY3HTzIYtffKxtxpFmfcttcJxPdDmIZQ2edBaWlySeb4TBMmaINoM7OgUJf/9s4icaabqPZOARVs1lTTMnSbFXVYBG7hQuda3byZDj7bNhoK5bd2AhPnJy0PsasZqiogpECR74B37gRvhx5mbV7Eod0Wv2xyyIFv7LJg+lYsx57MQOvWbc7JSSJdniDOvZWrOLQ39VBAcx4AW78Ijy3bLAIWKIe2lZ9hh/U7mZPORTVZk8KlOMUJY8jRgKvWbdJcJ2wqs+P/I7WbZHUsWimDp3v7IcjRqZQr8GEzucMqRjvzcAiYIJS6nvAk0DO/5udQYgVrMyqSYWnNLXBwrfSyzFOE1OsLjGONevx6ronoXHWZ5k7bnhajDaGQOcVL1wIY8dqr7plFDU1aSOpsFAb+W02D9vkyfDoo0yvuctoNgFB1SxkUTGlUAhefRVuuWXwWHe3jv4oLh5qyFuRIWPHas1u3KiNKdCLVo2NOgrFQbRDdAG7dw6D9hug8p/6fryQzmwsUmfHsWZ98r4HVrNud0pI0gp1Orr11rEV9Zx0M/QXQlkXfOF7MOmtQaMolifTihTpopMSSillRGbH7jJGswEmjm6t6vO7/9+tbDvmJbZfeiuhSgZz3ycO3k6nRachO3HSKm6WUkqJyG+BhcAPgHeAeSLygNsD9JsVhDmD2YlzMUEXq5u9Ov+K1tXact096BtuitUlx7FmPV5d9yQ0LpOfyW4MldgyhNrbtVd94UJtFC1cqI2k/kggklKDBlBrKytCGM06wJFufcghzjpqarT2LM2KQG/v0MUnKzKkvX3weUrp5+3dO1icLoVoh1nNoEZAYR8cvgHOvHMwpDNWGHI2FqmLxpFmfaoxEmjsnRIyXX0+SdG6SkJMR7femlgTZtx0fbynBD7zDbhpRXxPpj1SpLZhDEUdXVkXghxUzXpa1DbLmBiC0UfvYMztV7C3YwcvNjEk9/2eTQ7rNZiw+dxBRBJuwC+B54H7gUuAg5M9J8jbzJkzJRVOl2oZIcjpUp34xOpnRQjrfb7Q2Ko/M2GRkuWevKXjv4dLAKslVzTb2CpSulz/7RpbU3qPdKmeI0KF3mecxlb9WUpd+jyNjSJlZSIVFSLaJBrclBIpKBApKdHn2TCadYaj76mxVaT8cc/0KuKyZt2mpUVk8mSRoqL9NQtayyUlIpMm6XNjMfcFfY2f+0LSt9vUInLdMSI3TBa57kMiR/1chKUiZQ+LtGwbPG+5tEiFlMkIQY6Xqpiv5bVuXdOsiOe6zQrNtrSIlJdrHRYU7HfdHBblj2vNlj8e8+EnpVqWCPKkVMumFpEHqkR+OEbkV+h9440i5Y+INK4bfM5+mm28TaR8qt63tIhUVcmHoFMCfq0NqmZFskS3bpFEs+EdWrN/frJaflUgcv+jWqMsFan6u77GWrfjEhDNmm34m5NWcf8lIicA3wXGAvcopZ5SSn1fKfXRSB58zuI4RDsfi9bVeldh3sIUq0uOY836sLruamjc9ev1ZxHcSd+oqdFeyR07Bj2bRUUwaRI89hj09cG+ffu11DKadYbjHOIEXjU3yOpwzlAIWluhp0drtqwMKioGI0N27NCabWuL3wUhhWiHiSF488vQVQbPVMOc++CI16HyzaFezFzwukNwU7iyQrNW9fmCAh21lMn89yTpNVbo/HTqB9I+Hq2BdR+Ev34Byuuh4o2hnkxLswOt4WquhIsugGsbYM58WLuWUoIfRx9UzUKW6NYtkqQqHFtRT8nmanbfUA/9UPjfQ0Pnj4iU6ElYed6qOH99PZx1TtZo1rA/jnPeReR1EfmJiJwLnInOef9P4Bm3BhcEHLcxyseidVaR1RL3K8xbmGJ1yXGsWUiaH5hpXA2NU1F7N6mp0UZPT09iwwej2YyToBWUG+RMOKd98ckeGp+MFK8R//U5+NVNMP4d2HYwfL0BjloLc7fqx3O6NVw8jGZjE53/Ptw2hxZJFpwqCTGblVTatLfwa3D3d+DYZ+DlE6G3OFJ5fu1QzR5habapCX71c5B90K8tJ8mlIs4eazYrCYdh9uzMaDaJM6WSEHMmrmTkO1qzHRvgD28OhstfcthgvYa4heuammDhVdC9DegDckyzeYRj410p1ayUelQp9TfgRuBtEfmGiJzo3vCCgeP2GvmS997UBgf8Hc4aq68WHnndIdirxkHCsWZ9wJXCddZKdYk3HQ9SwWjWGY4160MxpbzvQ5xCtEOoEv54Gjx0BRz5pu6jPf9e+FsYfvBE9raGi4XRbAaw57+n25ozGgff93uEWcls3osscIcq4f/Ngnu+q9cSZj4BtVfBIQ/AV9/7LZ1WpMjzl9haKvYCHVDQCVVVvAXrhjdw9zGazSB1dbBqFcyfnxkDPknEyHuEOaB5AYXlvfR1wmVfhBNeGjTgLzokTuG6IW1ABW2vF8KkSVmhWcP+pFJtfi3QBPwM2Ar8Ril1lZMnKqXuUkptVUq9EufxjyulXlJKrVFKrVZK/Ufk+HGREP1XI49/xvace5RSGyLPWaOUOi6Fz5ISjiffdeth1S69z2WuX6+vFi07YNdpnnndTbE65zjWrA+r666ExjW06qJapd50PHCK0axzgpzuMe88KB+t93lJiteJUCXcMB0e+xL0jtDVvD9xN3R+cxyHhk9iRBa2hotFkDWbVb2zm5uHtuacM2d4PeAdfN9vUMcOVvGGzYhdvAX2AL+7Gs7/I0zaABfct5uLzvkm08Jn8ul7ezjjxGttLRWV3i6/BF59lZ3Qkf6gvcFoNoPU10N5OXR0DF+z4KjNYXfoViYuupqCMih5D668cdCAv/+dqMJ14TDMmAFz5w5tA1oxFsomwFXXZoVmDfuTStj87SLymIgsE5EfAicCX3P49HuAcxM8HgY+LCLHAZcBv44c3wt8UURmRJ7/U6VUhe15NSJyXGRb4/SzpIrj9hr5kPfe1KYv5uBNeLKNIHuTg0ZKLWGStIIKPE1tevwlyrPwf6cYzTonpbacHreMW7wUOnbrfV6Sxve9eAs8fyz86gboPQr+dWaYP916FR/6y2eZ+O4Hst5wh2BrNqt6Z0e35uzv117C4RhDDlvGTbcZsfXT9FN2jIKll8CuynZ2H/RvJj//YT5xzS3cehm61KPVUaT0MOBAuH9Z+uP0GKPZDBJdt+Haa7WxnK4X3qFmTwh9nAOO0MekCy6/Rz/N8rx39MH8v3cQvub7lGew8wAAIABJREFUsHat7jQCgx1F+sZBp5iq81lMKp53AJRSVyilfoj2wO9y8hwR+QcQN5ZcRHaLiJVBPYpINrWIvCki/4rcfhvt8R+f6piHi+PWW/mQ925fgfU4PNmEH6eG45YwDltBZYqMh8bVrtfjVwTK6w5Gs6kS5DaHeVtICdL6vi1D6Plj4f/uhT/+4jtsmL2KZdf9kk+edzd/mgGbM9ghzC+CqlnIMt3aW3NaDMeAd9AyLjrvPVQJi2Zq3U6TJs48dDwfqriKbR9aycRxP6O4v5cuxrLzy7fouhEN34LyUboIWBZhNJtBrLoNBQV6YWft2vRTP1LQ7KxmKIgEq+zapQ32dx56lDe/cDg/uuN/ePj6c/nrKRfqE0pK9MLYY4/pNJXaq/nJt7sZ9/5G1LEck/ZnN/hGysY7sAx4DTgM3fM9Iyil5iulXgeWor3v0Y+fDJQAb9kOfy8STv8TpVRppsYSTUreMysfPFu9mMmwVgYbp3pqKJnw49RJKbfNQ+97xkPjJGpvyFqCqtmsKQDmFmkUtgxVDlZDPmbi83SVK9R7J3LSz29h/Ksn0L4WHj4bXhxmpKnfBDmHOCupqdHeQYuFC2HkyPSM+ATpHtE57wCEw4ROm8FT/3MS3/nLtRz3Ui8ztj7NnK4rGNl2MluoZjkP8Off1OiFJ6t6d8OPoekXqY/PJ4xmM0xNDTz66NDUj7PPdlWzE0NwzhIoKocP7F7CH+Y9x29vuJiDt21iwZJbmf36Kup+36Bz3Zct0wtjoRArCHML3+PyW4v42k+LoZiSYX12gz847SkHPAAcbbtfALyYwvOnAK84OO+jQEvUsUOAN4CPRB1TQClwL1CX4DW/CqwGVk+aNElSZbm0yOlSLcslTh9cO0l6NWY1PvT9tPC7V7YFHvUeHq5mRYKt24z1c3W7t/swCYJujWYzQ173IBZJ6/pfvVKEpSIFTz4lIwQp/PtaYanI7P/WPbWt7XeTdI94O17o1mg2wDQ26v7v2p+pt0mTdH9qp8x9QX/fc1/Y76GBXu87qkSqqkQmTBBRauC9vvC/U2XVyWXy8y9NlW/VfUU++D2RW0ttulUiv58ssqnsYhEOFimfmjXXWqNZl2hpESkrG9SrUvp+Y6Pz13CiWanWr1lSIn2FxfIEt8qvEPkHv5Q9aqx0VU6QNydXyaU/bpGWbfq5y6VFpslkKRMlS86vEOFg2Vk+QdSx7BMPNGu2zG6peN7/D/iDUuplpdQ9wB9xocWA6BD7I5VSlQBKqQPQ3vjviMjTtvPeieh5H3A3cHKC17xDRE4UkRPHj0896j7l1lu5ulJpFaq73vuCfPkWfjxczaaMx7ltGQuNc7u3+zDJJ90azbpAUxuMWAFq+dDtPNdKvMQnjeKW9dOgoHwNhdO/A/vGc5gcDsBLH4XCwwbP27NRe+Hvn+JtKL3RrEuct2aoXguWQ+mK1OZEljdz0qTBYxs3wllnQXGx9igmC02OF/YdDjP9c+sYuwo+dMZaHeq8das2uYA3phfwl0+M5sxnOrnuZwfy4Jl38PpxcP1v4N8Xot1GArvb4OHOe3mBBey8yLvQec/ntB5H5vkWNm/ptmh5ep/Vqt0waZLOLxfRXviFC7Vmp0xJX7NNTcwq+yezqmHWiGf0a3Z3U9DXw0xuZCwv8h4nc9+Yl3nytS3Mu+9V7p4eGmh1+EkuYCNtCMI3b99FR3k/a2pPQ17i5dQ/qMFvUilY96CIHAt8HlgOPASck4lBKKWOUkqpyO0T0N7095VSJcAi4D4R+VPUcw6J7BUwD4hZyT5TOA4z8qFCpyf4WKjOkB4ppXv4kNuWEbzs7Z4iJtUjdYKsWU9D55va9ARy4Vu6nkM0y7YPb5KZDmlM4gsqw5Sc8ClUybuof32PO2eM1gWVSuGbv4axN0BhJNJ0wBg6C34/GcoYU+7aZ8kgRrMRLM2q5VqfdgQ9f1j4FhSmoNlQCNradBi9FZIM0NurDflzz9VGUWEhjBixf5jyqf8ELoXd58PYsfrcoiI4+2wqf7+V2bNhjLUOVlQEpaV0TprA6U8LnZuvon97NZWt36d5BpQVwM6R8Kdr4KRbsP3mFPAc1/HAr6/kMI7NmvzhlOa0HtbF8TxFKZZu+9BaTecaa2n2scf0ApNFb68+fvbZWq/FxVqT0SkhlmZ7Pjmo7eJiuPZaCrp6OfApKNjXN3h+URHdE6azdfSHePLDJ3BL3URuXMRAGuFbXbu5cOf1dNIJgEJx1cSbKa+9gY82PM+hFEQ3ljNkASnnvIvISyJyn4jcIyJbnTxHKfV74ClgulJqk1Lq8kjhuysip3wSeEUptQa4DfiMiAjwaXQY/SUxWsL9Vin1MvAyUAnclOpnSYWUPGi5UMHbTlObvpBZ+NBH21TtTp2UNethxEhGitY1tekfqNLg9XYHo9l0CPp11pM+xNb1ti/5qQOTTC888WlM4q9hAYx6iwIpZcL7nwV0QbCyAu3ErzsTLt0LJzUyZAFuz0YYz5FHZf5DZJ4gX2chgJrtJ3UjvqZGF4hrbNTFtyx6e/XW3w/79ulq33aDZ9lCYCPIbmhv1+f29Q1UjJeiQvqKoG/yBHjkEVZ0LWVSWyftFULxxD8w9pW/UtauF16PiKwdbNgL710KH3sMRk2CglKAXhAopDhr8odT0q3HESOeadZaIE2kW+sam4peQRvxra1as6WlenHI8sb392sttrcPeuatxSVLs117BrXd2zug2X0TFL0jYffsCdDSAj09jNnyJBcuLmTNLPjMHVDyKjTPgNKD/oqceA7dxe8CUEYZS3mM/6FG12no2MNBFBw6jG/R4BfJ4urRBv6RQGHU8ZnA+X7H/ae6zZw5Uzwhl3LfS5brz0LYt7zilHK0XASPctokxzXb8rjOaWsZzluVRnRZujxj48oUy6VFjpcqOV6qjGbdxOPrbEZ0mwgr39G+FUZddxtb9bHo80o8qPuQQt77cmmRCimTEYKMfuUOYalI+SMiLdtEqh4XYaneW2xq0TnEd5bobTIzRYxmh43rmm1s3V+LBTE0WxBDs+nOKVpaBvPUi4p0brwtX33opkSo0FtRkUhhoUhpqUhj49AcYhE5XqpkhCBjpUyWS8sQnbZs0/plqa7lYOffJRfLA7wmRzJDJBevtfmgWWub+0Ls63Am5sCRPHUpKNBarKhIoNmvijBh8NxIvYdozdq5s0rXY/jxGJHm6e1yzMMXyghBSrtH7z8XabxNpHyqHErBv8VjzZpt+JsTz/vdwAvAc0qpC5RS9yml1gP/BD7jxoJCUHHcXgPSqs4bSOzh8qXKt7zilHqTGgZISbMerq5nJDQuwFXm66njNdYyhgqj2RQJqmZd57w1Q8ONS5Tu6tF75tDrbs1kfaxx6tDYOSssOSDRXtewgC46GUEZ3z1o5kD/4bo3tVeoahSgIPyePn9iCC5qhcv2wTnLoJu9e/wcfyqkPDfIlbo4560ZGpUHWpd9MTTbdybImTA3qpVuOpq12spt2QI9PdqbboUpFxXptl1FRTr3+ORfAn+GuSv0ub290NUFNTVDer2vIEwrGwCYwhH6um1Ly7J3T5gXFWh82E2z+FT5+XTx2qbUPoh/5K1moyNJQWtSIvpcepzeJIN6taip0REifX1aizt2DKaEVFRozZaUwNxvQclFUHo/3Lxen9vWBqHQEM1Gc04zdJeBUrspe2MM591wI9JfxNjXf0PTe6+auUgukcy6B/4XmBC5fTGwETgdOMTvlYd0tuGsUgahcrTnBMDrHiTIMi9mSprNpqqyAa8yH5RIERGj2UziWiXkaE9QjErDjp/r5rXa4fe9XFpkbMTrfrxU6WGu057LxnX6nFjedzte69bTuYGHuvVMs6lciwOo2Wivu4j2tlevlIGK3Vb3BCuCZPDz+OfFTFe3RrMp6C7W81K5RqdDGt/3cmmRae8fK9Mev1Cuq35Sjl30FTnlz0tia7Z8qggHy/EU9YnHmjXb8DcnnvefAxdGDP3fAL8GnhSRdzK9kBB0Uq4cne093wPidU9phdgwhCDnY6ZdVdZaOe8WXRwygFXmDekT5Lx3VyohR3sv547Tnh+n1EzW3k47bnngHXzfP6GJ85lDJ52Mppwf0QzA4i2RYvXrIt72iFdzQ+eg9z1bSXlu4GHEiCuajfZezh0HXWc4vxZbmi21FTpwU7NxoiCtvtktNLEh2uuO9ravrNZ70N0T7BEkA2Rh/nDea7Y0EtnkRLOWXktsel223d1aI3F0a+/1bmegovy4l2ib3cJPHvo6H1z2Gb588fmM2as1u2Ct/fWvhvJRbKH/bfc+hMEtnBjv3UC9UqpDKdUN1AF7lVJvK6W+4O7wspw0WusEilpbSzgfC4KZwl8e4XGnhLRD5+1jC2haitGsR2R7JeSmtqGh8o1TUzPcLbwyhhx839+ngX76KaCAP7JowBCqnzZYqG7BWmiu0vc7+2HBq5kdptek1HoLPK0670r1bvvcoFSlr9muM4YaRG5ptnaK1mvUa79BHTtYxct8i66oxaZYhCp1wcXqCq3nAbLQEMp7zaay2AT63H1nDA2jX7bdc+ecpdk3ouYW17BgoKJ8gRIOevhHnHZ3CDph/n36nA17bQulNVfCrnW8Tf8WD4dvyBBOjPfvA9cAY0WkBBgd2X8U+LqbgwsaKU/Iszkfs6ltsEVRiX9ed8ivXtmZJuiaTauqrDW2ueMC63U3mk2flDWbjZ0SLK63TSiHq2fLGLIb8PYJa6ZI4MlcQZhxjKWMMm7i5iE5lv+fvXuPj6o+8P//+szkQkJCLgYRogYpkhC5CUhdqd9qI22trSn1++tqd3vBS+uXWn+7tbFdlpvAWiv7rbvsVlup2EpLra61amtXJMXWlWoBEUEgCAhqAAEJ4Z6QzOf7x5lJJkOuM2dmzmTeTx95TOZ2zpnw9sz53KtK2mfsxgbv5zp3+0Pre59ym8qZDb82gNgr9hdFvH92HDLbRUOKM9Y9h4UEGEBOh8qmrkS2xgNOQWjOt1Oq5R2U2aj8fkLHAnw8zrHQbWYjx7w/wGK2sRVwZpRfePw5SoqqMMHz64U7YSBORWmHHiOLH+RSMi6NzweQeOpN4X0LsNda22KMuQr4UvDx08Cr8TowL+rzBXmqrp0NHU8YkV+uCaS1smMTdWYTVKP8+esgP8+57ZXwlkqP/n+lzMamz5lNcI+RPme2K9e90XFYUjStl50JvzBtsgltGbqLO3mXdxnORc5yRBGWXAKVec7vtYec1vd8fycXlSmoT7lN5cyGdz3ubbfj7kQO+2iOQ2a7KHhuwpmR+RRwUVh3+ags/CG+KJZfTiYvZ3bBLKgshyONMba+R3aXdyOzv5/Q3mOkycKA1fHJbCdDlEqoYiqvUEIVq6llIpfwz3wXi8WHj6d4jrsKqlh3LVz/G8jIhwu2wPf/Fq7fGNFjJAUzK47e/KPNBc4zxnwLyLXWLgs+Pha4O25H5kF9nvE8Vdd7X7zHOebQbMdJbN1U9+PYRJXZkAR8Qf/293DsuHPbK+oy3+9FtbJEAnuM9DmznYnsLu/msKTIwpDbXZG7aBEKn627K1UlUJgBW07A9PXOY09Pcgr0R86kdut7n3Ob6pl189qgs8y6OZ64i4JnqKtxDjnddpfvlTnfJuCsYp8yvHxNW/VxKCyALXUxtL53Ns7drcwuiqgkdbvHSC+GKN3FnWxlS1vBPbK3U2kVTHsa8EHOcfjcbKjYGLaBFMysOHosvAfnJfw1cBD4pjHm18aYOUAWUBrvA0xpCR6P6YrQya7JOsee5G7J6n4cuz5N+Bc+djYB3eP6PDGNusynhT5PUpnAHiMxT6bU2WRfbme5pix+49+7aMUMFYQG9FAQipz0q6oEsE6BPpXHvvfrzELH7sHxymx4Ad7tv0sXwz3GAItwWqNiUjOTDbRsiHUzidTnzIb+zeNRWO1EzLmNHOcez0rSePQY6aKC7xC1/AeX8A7O90gOOfyOlZ32diqtgsvuw5kgNAB/uTP8M6ReZsXRY+HdGBPqaPUD4AFgCXAA+BTwWPwOrZ9ItXHv4SfkJLdsqvuxO/rcEpzA7nF9mpgmBbrMizuiGvce4qXMdib8gjLaCep6Y+GI+Ix/7+T8EN7q3lP3407Xy+4HM8/368xGzoETr8zWlHUcTxzHAmKowDqDHC7k1FkTgKWDqHqJmYhbr4rMbF8nqOuNyAK823ntYujtX7mTJWzhNE3kkMNTPNftOXd8DRSOdn4/shU2Lnb3MCXxetNtPvTPvMNau8pa+4q19ifW2pnW2ivjeXD9QqqNew+dkJO4NFyIuh+7I6qW4ARWOvV6YhoPVSx1R7mNXVTj3kMX/V7KbKRETgQaOZu3m+Pfw1oxV1PLF5l+1tJw3YlcNq4/zDzfrzMb3lMk3nPg/H5Ce6WTm62ZEcM9Ql2O/8jgsyYASxdRXRuEKgUtce8xEtOkdeGVlfHMbHgvJ7db3zvpLbKaWr7OOxwP3h/ey7ka/mYJTonPwtrvQb1WXk5pvSm8lxpj/gnYaoy515jQ/IXSKwmeoTMmi/c4J2S3uxdF6Xo+Tx75XM/nk30o6SeBlU697hrnoYqlrqymlkaOMJpKdZuPQZ+XMQJvZjZSoi4ow4Xvx83W9+DyW/9UfyfHOUZOL2frhohl497qHzPPp0VmEzUHzsI4ZDZsvPbqX69oW9d9F3ltE4Clm6gym8CeeVFnNtGrJYXn1c0hSp0sc3gXd/IhpzhJFqOp7PVcDW3d531AAP7wSbXAp7KMXrzmDeA4kAl8HLjVGPM+sBF401r7QByPL/XVlLXX+IZqrz1a8GD2LueEnOSl4UKe5bcc5xjP8ttOx/JI74S3BPdpYprZu9onpvFAHlg4wvl/yeOt7lvZwuVcoaEeMepzbucMdy70mzyU2XDJWn6zpqz93B5qfXdj38HvtTm3f4n//dzsXrcAQXthfctx2irlllQ6k9iFxsKnIi9nturjzm2oFTN0v1v9KbPBa7HVl63li5+Z27au+/18jVeYSjkL0rIAH9X1wZUFzhC2OPcYiTqziewpAu3ZDO139i73/l8J6zGyumZHW6VTGefxHxRS3odNja+BkolOwZ2A0wKfQ0G+OweaXOvXrz83IyPjpzjTWPSXGfQDwOaWlpZbJ02adCD8CWOt7eI9XTPGXIDzBxpjrU2pupvJkyfbdevWJXan4SeTfD8c7c0ZKAkGrA5OVBccH5RkXhzzboxZb62dnMh9xprZqP+OCcrD1E85XeOumAKvvBC33SSEMutw4zwb1d9y0J+cC504n2ejymz26val4RK9ikf4d1CWgSYX/n9evIc//PmXfPPHP6Sh9CS/6WHcZaTaQ8Eu8sYpuFeVwOKdTlf6OSPh7pGJzW1/zyxEkdvQdwAkN7Nu7X/xHiZ+9uNsHb2nbaxwFnNpYA1FXMFUXolp86l4rlVmXRSPfV/3Bjx/mNX/9DY33PuPbasj3MNFXMSWqHK7cTGs/S5gYT5jT79vN+XEfqDJtXHjxmfPO++80YMHDz7q8/n6XrD1oEAgYA4ePFiwf//+LePHj78+/Lmoaieste9Za/+QagX3pEnw2LY+W7zHORlXFTknYw90mfdiIShVRbX0FjhDKMJv48SVmZA9QJl1l5eXjOtTZhfv6VhwT8awjziMy9zBDj72p3H87YpPRLVGdlUJFGY6re+hlvbwsfAMTL0Woagym8Dlt/qU2+veSE6re0gcJgPbwQ4C1lkZK9RTpJwFaTvmHfpZZpPVUyQkHsM9gkNq7vry9znFKcDJ7k0siTq34RPY+cnKcudAk25Mfyq4A/h8Pjt48OBGnMbyjs8l4XhSXp+X1wBvT1w3e5dzxbSqwalF9UB3U0365a6oMrtohPMFaIjrF3TMMyF7hDLrrpjOs3FefqtPmZ2zq73gDsmrHHV5XObQhU3kH8tl1sKvRL1G9oJRcEWhcxu6H1pGjvNGjIzpAJMgqswmePmtXolc0z1R8zNEcrnSaejCJh741p1c9mpFW2ZLqErbMe8hUS0ZZ3Ay61YhNVbJ6C4fqabM/QlC5wynJd9yeLAzRV1oKc5Yc/s3SyAjHwymv5QDff2p4B4S/Exn/Rsl7B/NGLPMGHPAGLO5i+erjTFvGmPeMMasM8Z8LOy5rxpj3g7+fDXs8UnGmE3GmB3GmCXGmIQsXhHVRXoCayr7ZPGe9gtLDy39obWy3RVVZmvKINvnfAnFeWKamGaV9Qhl1l1Rn2dDvJDZ8JYgSHw3znCRLZkxXHSvppaJW2/h+ev+Qv35B7l6cXTl7KoSeOWK4FrvwftPT3IK8Bhfyl1URl2Bl6Dlt3qd2fBCUDzWdO+L8EqnGCo3VlPLfXN+yeX/cwkvXvVvUWe2P4oqtwnqmdfra4PINd2TldlF7uQ1ZHXNDsZu/QrfeOBzDK8f1ufhSV0prYJpT0MzJ0/EvLF+Ijc399JkH0NvJfLL8WfAp7t5vhYYb62dANwM/BTAGFMMzAM+CkwB5hljioLveQi4Dbg4+NPd9l0T1UW6F2sqoeOxeKC7fEjUXb2lU1EXLL3YDVnSQtTn2fs/kpAVPnqV2fCLt2QW3ENcaBkKLQ33bul+7pv9Cyq2lrlaURIqwNOUeheVUZ9nF45o/3fxWmbjtaZ7b7nU+n4Xd/KvNStoyW4luymzQ2YPUcsrTOVQX3pM9CNR5TZBPfN6ldnIStJkXsu63FvkLu5kd+lenrxpNZtGP9Z2TexGZkurYB9bt8V0gP3cmTNnOtwPBAK0trYm6WjaJazwbq39M3C4m+eP2/bZ8wbSXp/3KeBFa+1ha20D8CLwaWPMUGCQtfbV4Pseg8SsKeb1McS9luzxQZIwUS0JAwnrhtwfqNu8u6LObE2Zc45tjm835B67zof3avLS8oYxLh23gLkc5xh55HPPnjlnrUPshqoS4L3Uu6iMKbMJqNxP2cy6OOTjp3NeOCuzdcFJ6+p07u698J55bi6PFqFXmXV7YsNYuZTX1dS2zS7/5Z99ioxmX9u2lNn4+d3vfpc/adKk8k984hMjL7744jF1dXVZw4cPHzN9+vTho0aNumTnzp1Zv/nNbwZNmDChorKycvS11147orGx0QdQWlo69vbbbz9/1KhRlWPHjh29efPm7Hgco6e6pRljphtjtgG/x2l9BygF3gt72fvBx0qDv0c+3tl2vx7sir/u4MGD7h94byWoprJXvDA+SLoUj8zG3A05jl/Q6jaf+jyTWfBGN+TZ3uzVFGvr+/V8njzymcUcrn53ahwOMHE8lVkvdEP2aE+8WCevC1WkjKaSz3XSxpNqk9a5nduoM5ugYUq9Ps96oeAOrk22eBd3cjo4u/y3Hv5ihyGMqZbZVLNly5bcBx988N3du3dvBnj33Xez77jjjoM7dux4Kz8/P3DvvfcO/fOf/7x9y5YtWydOnHhy4cKFQ0LvLSgoaNm+ffuWb3zjGwe+9a1vXRCP4/NU4d1a+7S1tgKnBX2hi9t92Fo72Vo7efDgwW5ttu+80nXeizWVYaKa9KefiUdmY+qGHBKnL+hU7zavmeY9lFlwCh7ZxikIJatLZ6jiwEstmCGLomsZWk0t97KQ4xzjWX7bYQ3iVOSpzHqhG3Ko4sCLPfGi7I4cGuaxlS0UUMjIhRlnZTbVJq1zO7dRZzZBKyml5Hk2PK9RVCKHt7oP5yIy5ozs0GMk1TKbasaNG3eioqKiOXR/6NChzVVVVScAXnrppYE7d+4cMGXKlIqKiorKxx9//Jx33323bdb+r371q4cBbrvttsMbNmzIi8fxearwHhLsYj/CGFMC1APhNRfnBx+rD/4e+bi3hb4c3ZqJMhrhF1oeK7iDuh/HS9TDPRLwBZ3qM84rs/ERU2YT0HW+WwtHeGbpzbNE2TJ0F3dynGMMIMe50PfqRKxJFFNmE9ANuVuLgpn1ak+8KJbiCmU2JzyzCZgTI5XENMdQAobWdXt94OXzbAxzWYRa3UOzy0ti5ebmBrq6b63lYx/72NFt27Zt2bZt25adO3e+9cQTT7T9A/vC5lo1xsSlL5VnCu/GmJGh2eKNMROBbOBD4AXgk8aYouBEdZ8EXrDW7gOOGmMuD77vK8AziTreqFuHw78Uk9FasXiPc6GVZTxZcAd1P46XmHo0JOALOpW7zoe6El+fmGk30kZMmQ21dri0rnlnus1sTZlnlt7sVB9bMldTy+5gS1Dbuu4JXJEilUSd2/BuyHGqdEr5zPZhyEdk62VbZrN8TuY1aV0bL2cWusmtlzMbZW/bTs+1nfRySvfMJstVV111Yt26dXmh8exHjx71vfnmm21j2x977LFigEceeaTo0ksvjcvEq4lcKu5XwF+AcmPM+8aYW4wxtxtjbg++5AZgszHmDeBHwN9ax2GcLvRrgz8Lgo8BzMSZlX4HsBP4Q6I+T9QtbaEWj2S1VszZ5ZxIDN482aGZ5uMlptbhBHxBp3LX+Wf5bXtXYnFNTJldGP+K0lTOLNDrlsxQ1+NTnCKP/I4tQXOGx2XSulQW0/VBDF1teyPlM9uHpbi6bL3sJLPpPgGYlzMLKZzbPva27fJc20mPkXTPbLIMGzas5Sc/+cnuG2+8ccSoUaMqJ0+eXLFp06YBoecbGhr8o0aNqnzwwQeHLFmy5L3uthWtRM42f5O1dqi1NtNae7619hFr7Y+ttT8OPv8Da+0l1toJ1tq/sdb+T9h7l1lrRwZ/Hg17fJ21doy19iPW2jvCZquPu5hah5PVWhE+u7xXZryXhIk9s/H9gk7lrvPqLRIfMWc2zsvGpXJmgV63ZIbPMP8ET6titQcx5TbO8zX0i8z2cex7W+tl+DbmDHeuv4LvT/cJwGLObAKWOkxJfext2+W5tpMeI+meWbedPHlyA8BnP/vZY6tXr94Rery8vLz57bfffiv8tddff/2xzZs3b92+ffuW7du3b/m7v/u7xtBzc+cuzV8aAAAgAElEQVTO/WD79u1bNm/evHXMmDFN8ThWz3SbTzVRLwkTkuixgppdPu3FnNkETAKWil3nNVld/MSc2QSMfU/FzHbQi4vLDjPMR2Y8xSet8xxltme96FUTPsN8p2OGI3KrCcBiEOrF2aTMnqWPcwZ1e66N6DGizKYvFd5jEFOXzkTPPO/xSeokMWLOLDgXlXHKbCp2jdNkdfEV8983zmPfUzGzHfRwcXnWDPORNAHYWZTZOOuhV03kDPOdVqqGsh7HWdJTiTIbR72cM6jHc614Xn19/aahQ4e2xHs/KrzHIOausgla1xVor7FTwT2tuZbZZK6W4DHqMh9fMf99o5ihOu10c3F51gzzkbqYACydeT2zKd91HrrtoXDWDPOdCWU+eJvuk3+5mtk4rJaQ0pkNnw+km3Nkj7mN6C2S7plNZyq8xyDmidUWjWgfuxXvgpCXZ+QM0vru8edKZkPi8AWdsl3jJG5izmwfZ6juq36R2fCLyy7+vz5r3HDk+zVpXRuvZxb6SW57aO0d3lNmw1ru033yL1cyG778ZBwqnVI2s72YqLrTlREiRQy3TffMpjMV3mMQc2HTC2sRe4i6H8efK5mN4xd0KnaNU27jy5VKvUXxa8lMxcyepYv/r3scNxz+/ogJwNJdzLmNc0Vpv8htJz0U+pTZsB4jmvzLpeuDOFY6pXRmexgqG1oZIae7dd0jJrtWZtOXCu8xcOWiPQFrEacKrZUdf65kNvIL+ro33Dk4Uq9r3GpqaeQIo6lUt/k4iUtmXTzXplpmuxT+NwoOj7mLO9nKFoCeW+Q0cV0HMec2AS2ZKS/WzIaNe9fkXy6da+NYUZryejFUttveIqDMCqDCe0xcGesaOU7IxYJQqtFa2fHn2vjs8C/oHiZh6atU6hq3gLndT4okMYtLZl1uyUylzHZrUXDJJwM7FteyO9iNs1c0AVgHruQ2jpVO/TGzq3+9om+ZjRj3nu6U2TgLy2r436XXvUVAmfUwv98/qaKiojL0M2vWrPMApkyZUl5XV5fl5r5UeI9BzGOE4Oza9ecPx35gKUoTf8VfzEtvhYTPUA3wvZ1Q6052U6lrnDIbf65mNvxc6+JQpVTKbLfClny6YHYLpzhFHvk9X1CCLirjJbzS6bvpeZ7tVlhm615+tW+ZTfSSvekiTq3vKZ/ZTpbU69XKCOEi5hfRpHXekZ2dHdi2bduW0M+99967P177UuHdC8ILQuEFojSitbITx7Ux2r+f4BSGfEAAmJteXeSU2cRxLbM1Ze2ThJruX5q2gn+XlowWRh8ayRM83bt8a8m4DlzNbHjXcJfOs/1muAe0ZfYbP6rm7sVf7nyN7M5EjCFO94JQ3DLrkn6RWdPxtlcrI4SLmF9Ek9Z5X2FhYYvf73d1XTEV3r3i9xPAfsK5TUOa9CtxXG0trimDlRPgikGwYETPr++FVOkap8wmjquZXTjCablY6E5eIXUy2ysLR3Ai9zSL5v+cQTuyel8xpSXjOnA1s4tGQI4PyrJdO89CP8rtwhG0+FsxGGb/09+Ts/hg798bNtwj3QtCrmc239+xFd4FKZ/ZhSPaKzbCKjl7HOveYRu72+YX0aR1nbh56wVMWVve40/5XyrJe2kC5X+p7PG1N2+9oKfdNjU1+cK7zS9durQIYOXKlTtHjhx5xs2PqMJ7jLS8mTvU/TiFVRXDK5OdWxekStc4ZTZxXBmiFBKHZTNTJbO9UlNG7ec3MHv+1/j2f3ypb+/VknFtXM/syatg91TXzrPQj3JbU8ZPv/8iR/NP8Mg3nuPWhZ/q/XtDwzyeP0z5r7+Z1gUhr59noR9kNmzW+dV//C+gl2Pdw4UN9yhZPFKT1kWrvimLEwE/9U2ujEeP7DZ/2223Nbix3c5kxGvD6SK89U1dZyUVeD2zVR93bkM166H7kt68PEyhP2V2NbX88FsryN+fzedWXA4T9rh+AZ4uvJzZ/mQ1tTxc83se/uJT/N+b/w8ZV57T+zfPGe5MYAmU3FZKyd++Ep+DTBHKbAIEO1DfM3cZW9nCaCr79reuKeu4uofOzx0tG/1er173zMF85r8zjPkX7aV68LE4H5Wr1PIeI7W+uUNdkBMnFTKbCl3jlNnE8vrfOxUy2xt3cSevXb6FexY+isH0bWZ+LRfXgTKbGKHl4QbtG8DVf5zYtxVQQhNZauI6QJlNiEUjINswce2o6LehXk6xqx58jA1T6lKt4A4qvMfM1W5GaSwVCpT9RSpkNhW6ximzieX1v3cqZLYvjo5qbr/T28K4Zu/uQJlNrIl7xrTf6UsFUtjEdS0L30rbCetAmU2ImjJasgJ87pmpfPTVPnaZ70S6T7ToFZFj3mfOnFkar32p8O4CjXuPjbppJZ4yGxtlViL1h5mQO6w3XPLj9tVPert2e9h4TjeXiJL46G+Z/dzffr3vmQ25sgCL5dCVG9J2wjpwcWnOOOkvmR279Ss0ZZ5h5VUPcPXikX3fSFgvp3SfaNErWltb14ePeX/wwQfr47WvhBTejTHLjDEHjDGbu3j+74wxbxpjNhlj1hhjxgcfLzfGvBH2c9QY8w/B5+YbY+rDnvtMIj5LZ7zezcjr9PdLPK//zb3eNc7rf7/+KBX+5l7PbU9C3Y/BuYgPn8yr1y3pNuI2jSmz8edKZgFebsRgKHnpUsp//c04HGnq8Hpu+0Nmd5fu5Xs/fIispozohhmFdZvXjPPpJ1Et7z8DPt3N8+8AH7fWjgUWAg8DWGvrrLUTrLUTgEnASeDpsPc9EHreWvt8fA69Z17vZuR11/N58sjnej6f7ENJG17PrNe7ximzief1zIL3c9tn4eMpe3uBGRzPCaR913llNgnCMzu7D70/gu/LODmAkq8MdfWQUo3Xc9tfMvvVJz7r/NLXXiLQYTb/Eqo043yaSchs89baPxtjhnfz/Jqwu68C53fysipgp7XWc1cDV1OlrrMxeJbfcpxjPMtv+Udqkn04aSGU11DNutfy6/XZu5XZxPN6ZlNdhy7zoTGYNWXwUoPTitmXrvOaCRlQZuOty8zO2eUM3TB92Fi075OE8/r1QXfCM/utH97gPBjqJZLG50rpGy+Oeb8F+EMnj98I/CrisTuC3e2XGWOKutqgMebrxph1xph1Bw8edPNYxQVer+VNhkRkVl3joqfMnk2Z9XZmu7OaWr7IdLayhQIKOxYyo+mGPGe40/re5O2J65TZ1M0sdNJlPmRhsPeHpW/5i/Z9CRbv3Ho9s5C6uQ3PbMacsLHuWqFD+sBThXdjzNU4hffvRjyeBVwPPBn28EPAR4AJwD7g/3a1XWvtw9baydbayYMHD3b9uCV6mvirc4nIrNcLoF7tGqfMdk6Z9W5me7KAuRznGHnkn/23jaYbcqgFqdnbE9cps6mb2W7VlEGWz8lfX5Y6rCkjkBWAZou9e4dnC/Dxzq3XMwv9JLc1ZdFPsBhGs82nH88U3o0x44CfAtXW2g8jnr4WeN1a+0HoAWvtB9baVmttAFgKTEnc0YpbUqGGt7/y+pJxXp1VVplNHmU2PkJzOMxiztl/25qy9jHsfelOrInrAGU2XjrtMh8umvkagD1zXsBiMRhPVzzFk9czm6o6zWyoZ1PoNgqabT79eKLwboy5EPgN8GVr7fZOXnITEV3mjTHhM4pMBzqdyT5RtPRWdFKhhre/SoXMerFrnDKbPMqs+1ZTy70sbJvDoVPRdCdeNMKZDXnRCNeONVV5PbepmNkuh3mERNmqmV/zKXZ/8xkC/ta0rnhSZt23gLlnZ3bOcMgy0Bz9ECPNNu8dy5cvLzTGTNqwYcMAgDVr1uRMmDChYuTIkZeMGjWqcunSpW1DvKdMmVJeV1eXBVBaWjq2L/tJ1FJxvwL+ApQbY943xtxijLndGHN78CVzgXOAB4PLvq0Le+9AYBpO4T7c/cGl5d4Ergb+Mf6fpGtqjYuOaniTJxUy2y+6xolrlFn3ddtlPiSabshhsyGnO6/ntl9mFqKar6GEKi4quwtfblZaVzwps+7rdJWamjLI9jkTJUY57l2zzXvH448/Xjxx4sTjjz32WDFAXl5eYPny5e/s2LHjrZUrV749a9asCw4dOuSPdT8JKbxba2+y1g611mZaa8+31j5irf2xtfbHwedvtdYWhS37NjnsvSestedYaxsjtvlla+1Ya+04a+311tp9ifgsXVFrnKQaZTY6Xr+o6c+UWfd122U+XJTdkMX7uU21rvNxz6wqnpRZl3XbwylFJviU7jU2NvrWrl2b9+ijj+5++umniwHGjRvXNHbs2CaA4cOHnykuLm7Zt29fBkBhYWGL3++3AEVFRS192VdClopLB1ouTlJNKixjFN41zgvLwaymlkaOMJpKz17U9GfKrPt6vexhqCCzcHfHQpH0C6mUW2VWILUy221vEReW1zxELXXMpZwFad8C/w1uvuAtNuf29LqjNPr3Up81jNLmQRS0dvfaSxhz8icse6+716xYsaLwqquuahw3blxTUVFRy8svv5x75ZVXngw9v3r16twzZ86YysrKJoCVK1fuDD23efPmrT1/snaeGPMuIsnh9VZkr3WN63TMmiSUMuuePldGqUUyKl7PLKRObhORWc3ercy6KTyzT/B059cOc4Y7c4REWcmkSev6bi/1WSc44d9LfZYb23viiSeKb7rppgaAG2644fDy5cuLQ8/t2bMnc8aMGSOWLl262++Pude8Wt5F0tlcFrQte+ZFodr00KQ0ya5d9/rfKx14/d/Aa5ntTqgy6nKuUGVUHHk9s6kkEZkNLwilayumMuueXmW2piymStFyFrS1vKe7nlrIQ57jmfxFzB82m/l7P0f1sVj2+cEHH/hfffXV/Lq6upw77riD1tZWY4yxgUDg/SNHjviuvfbakfPmzauvqqo6Ect+QtTyLpLGUmHCQC/NKpsKf6/+LhX+DbyU2e54fVxrf6HMuicRmdXs3cqsmxKRWU1a13efo/rYa2yoi7XgDrB8+fKi6dOnH967d++m+vr6Tfv373/z/PPPb37hhRfyrrvuupE33njjhzNmzGhw47hBhXcR8bhU6RonEpIqmU2FC3RJDGW2nQpCqUGZFa948skni7/whS90KJxXV1c33HrrrRetXbs2b8WKFSUVFRWVFRUVlWvWrMmJdX/qNi8inlb1cW93PRaJpMxKqlFmJdUos+IVr7322vbIx2bPnn1g9uzZB+KxP7W8i4iIiIiIiHicCu8iIiIiIiIiHqfCu4iIiIiIiIjHqfAuIiIiIiIi4nEqvIuIiIiIiIh4nArvIiIiIiIiIh6nwruIiIiIiIhIlJYvX15ojJm0YcOGAaHHrrzyyovz8/MnXH311SPDXztlypTyurq6rGj2o8K7iIiIiIiISJQef/zx4okTJx5/7LHHikOPfec739n/k5/85B0395OQwrsxZpkx5oAxZnMXz/+dMeZNY8wmY8waY8z4sOd2Bx9/wxizLuzxYmPMi8aYt4O3RYn4LCIiIiIiIiIAjY2NvrVr1+Y9+uiju59++um2wnt1dfWxQYMGBSJfX1hY2OL3+200+8qI5UD74GfAfwKPdfH8O8DHrbUNxphrgYeBj4Y9f7W19lDEe74H1Fpr7zPGfC94/7vuHraIiIiIiIh43s03X8Dmzbk9vq6x0U99fRalpc0UFLR2+9oxY06ybNl73b1kxYoVhVdddVXjuHHjmoqKilpefvnl3CuvvPJkV69fuXLlzh6PsQsJaXm31v4ZONzN82ustQ3Bu68C5/dis9XAz4O//xz4fEwHKSIiIiIiIv1bfX0WJ044BXgXPPHEE8U33XRTA8ANN9xwePny5cU9vSdaiWp574tbgD+E3bfASmOMBX5irX04+PgQa+2+4O/7gSEJPEYRERERERHxih5ayNs880w+8+cPY/78vVRXH4tllx988IH/1Vdfza+rq8u54447aG1tNcYYGwgE3vf53G8n91Th3RhzNU7h/WNhD3/MWltvjDkXeNEYsy3Ykt/GWmuDhfuutvt14OvBu8eNMXVuH3sPSoDIbv/9UTp8zvJE7ESZTZh0+JzKbP+SLp8z7rlVZhMmXT6nzrX9Rzp8RkhQZj2puvoY1dWu/L+zfPnyounTpx9esWLFntBjl112WfkLL7yQd+211x53Yx/hPFN4N8aMA34KXGut/TD0uLW2Pnh7wBjzNDAF+DPwgTFmqLV2nzFmKHCgq20HW+sf7ur5eDPGrLPWTk7W/hMlHT5n+KSJ8aTMJkY6fE5ltn9Jp88Z730os4mRTp8zEftRbuMvHT4jJC6z/d2TTz5ZXFNTsz/8serq6oZf/OIXxbNnzy7dtWvXgFOnTvmHDBky7sEHH9x9ww03HI1lf54ovBtjLgR+A3zZWrs97PGBgM9aeyz4+yeBBcGnnwW+CtwXvH0msUctIiIiIiIi6eq1117bHvnY7Nmzu2xUjlVCCu/GmF8BVwElxpj3gXlAJoC19sfAXOAc4EFjDEBLsMZrCPB08LEMYIW19r+Dm70PeMIYcwuwB/hiIj6LiIiIiIiISKIlpPBurb2ph+dvBW7t5PFdwPiz3wHBrvVVrhxg/CWte1OCpcPnTIfPCPqc/Uk6fEbQ5+xv0uFzpsNnBH3O/iYdPmc6fEZIn8/Zrxhro1ofXkRERERERCRpNm7cuHv8+PH9coLBjRs3lowfP354+GMJWeddRERERERERKKnwruIiIiIiIiIx6nwLiIiIiIiIuJxKryLiIiIiIiIRMHv90+qqKioLC8vr6ysrBz94osvDgw9d+WVV16cn58/4eqrrx4Z/p4pU6aU19XVZQGUlpaO7e2+PLHOu4iIiIiIiEiqyc7ODmzbtm0LwFNPPTVo1qxZ50+bNq0O4Dvf+c7+EydO+JYuXTrYjX2p5V1EREREREQkRo2Njf6CgoKW0P3q6upjgwYNCkS+rrCwsMXv91uAoqKilsjnu6KWdxEREREREUlpN7/JBZuPkdvT6xrP4K9vIqs0m+aCTFq7e+2YfE4uG8d73b2mqanJV1FRUdnU1GQOHTqU+fzzz2/v6RhWrly5M/T75s2bt/b0+hAV3kVERERERCQt1DeRdaLVKcAXZHIq1u2Fd5tftWrVwBkzZly0ffv2t3w+9zu5q/AuIiIiIiIiKa2nFvKQZz4gf/52hs0fxd7qIRxz8xiuueaaEw0NDRn79u3LKC0t7XV3+N5S4V1ERERERETSQvUQjlUPoS4e296wYcOAQCDAkCFDXC+4gwrvIiIiIiIiIlEJjXkHsNby0EMP7c7IcIrZkyZNKt+1a9eAU6dO+YcMGTLuwQcf3H3DDTccjXZfKryLiIiIiIiIRKG1tXV9V8+tX7/e1RZ+LRUnIiIiIiIi4nEqvIuIiIiIiIh4nArvIiIiIiIiIh6nwruIiIiIiIiIx6nwLiIiIiIiIuJxKryLiIiIiIiIeJwK7yIiIiIiIiJR8Pv9kyoqKirLy8srKysrR7/44osDQ89deeWVF+fn50+4+uqrR4a/Z8qUKeV1dXVZfd2X1nkXERERERERiUJ2dnZg27ZtWwCeeuqpQbNmzTp/2rRpdQDf+c539p84ccK3dOnSwW7sSy3vIiIiIiIiIjFqbGz0FxQUtITuV1dXHxs0aFAg8nWFhYUtfr/f9nX7ankXERERERGRlHYzL1+wmYbcnl7XSLO/npNZpeQ2F5DV2t1rx1B0chlXvtfda5qamnwVFRWVTU1N5tChQ5nPP//89p6OYeXKlTt7ek1nVHgXERERERGRtFDPyawTtPjrOZlVQNapWLcX3m1+1apVA2fMmHHR9u3b3/L53O/krsK7iIiIiIiIpLSeWshDnmFP/nw2DJvPpXurKTvm5jFcc801JxoaGjL27duXUVpa2tLzO/om7QrvJSUldvjw4ck+DElR69evP2StdWXCid5SZiUWyqykokTnVpn1nmOcYTfHaaa9R+sAMriQgeSTmZD9v8sJTtN+7W0wXMygTvevc610ltnuMhMP+zlFPSd6tf9kZNYrqik7Vk1ZXTy2vWHDhgGBQIAhQ4a4XnCHNCy8Dx8+nHXr1iX7MCRFGWP2JHqfyqzEQpmVVJTo3Cqz3rKYTXyPtYTP8HQ/l1HD2KQcy3dZiwUs8B4+nuOTVDGsw+t0rk1vtezlc6ykOZhaA4ymkCVcflZWEnUspwhggRLO5RU+e9brkpHZ/io05h3AWstDDz20OyPDKWZPmjSpfNeuXQNOnTrlHzJkyLgHH3xw9w033HA02n2lXeFdRERERLyplr0dCu45+LmHiUkpuAPUMJaJnMPneJFTtHKKAJ/kv7kvSZUJ4j1OYflFToUV3H+QxHxUMYzn+CR38ioAC5iYlONIJ62treu7em79+vWutvCr8C4iIiIiSVfLXqazigDOWsZeKSA7haFp3MmrbOUIAeBu1gJ44vgkeSJ7ieTg5zmmJby1PVIVw3iLLyT1GCQ+tM67iIiIiCTdnbzKMVrIwc9KPu2pgnGoMPQDLmt77G7WsphNSTwqSabwXiIGqKTQEwV36d/U8i4iIiIinnER+Z4tAIUqFEIt799jLRM5J5mHJEng1V4i0v+p5V1EREREkm4Jl3MF57KEy5N9KN2qYSz3cxk+IIDTY4DRQyuSfVySGKGC+zFayCfDc71EpH9Ty7uIiIiIJF0Vwzzb4h4pNJHdXF7nCM2QmzUw2cck8RdZcH+aa1Ims9I/qOVdRERERBKqlr1M5XfUsjfZhxK1KobxCp91egqcbD7R8zsklangLl6gwruIiIiIJNRcXmcNB5jL68k+lJhVMQy27tuW7OOQ+JrL6yq4S6f8fv+kioqKyvLy8srKysrRL7744kCANWvW5EyYMKFi5MiRl4waNapy6dKlRaH3TJkypbyuri4LoLS0tNfjLtRtXkREREQSopa9zOV1Pk8ZoDWoxfs6y6wK7hIuOzs7sG3bti0ATz311KBZs2adP23atLq8vLzA8uXL3xk7dmzT7t27My+77LLR06dPP1pSUtIa7b5UeBcRERGRhAi1uAO8wmeTfDQiPVNmpS8aGxv9BQUFLQDjxo1rCj0+fPjwM8XFxS379u3LKCkpaS0sLGzx+/0WoKioqKW321fhXUREREQSYgETmcvranGXlKHMpo4/38wFhzeT29Prmhvxn6wnK7eU5qwCum0FLx7Dyf+1jPe6e01TU5OvoqKisqmpyRw6dCjz+eef3x75mtWrV+eeOXPGVFZWNgGsXLlyZ+i5zZs3b+3pmENSovBujLkAeAwYAljgYWvtvxtjioFfA8OB3cAXrbUNyTpOEREREelaKs0oLwLKbH90sp6slhNOAT6rgFOxbi+82/yqVasGzpgx46Lt27e/5fM508vt2bMnc8aMGSMeeeSRd/x+f0z7SonCO9AC3GWtfd0Ykw+sN8a8CHwNqLXW3meM+R7wPeC7STxOERERERERSbCeWshD9jxD/vr5DJs0n71l1Rxz8xiuueaaEw0NDRn79u3LKC0tbTl8+LDv2muvHTlv3rz6qqqqmFelSInCu7V2H7Av+PsxY8xWoBSoBq4KvuznwEuo8C4iIiIiIiKdKKvmWFk1dfHY9oYNGwYEAgGGDBnScvr0aXPdddeNvPHGGz+cMWOGK73DU6LwHs4YMxy4FHgNGBIs2APsx+lW39l7vg58HeDCCy+M/0GKxEiZlVSjzEqqUWYlFSm3It4TGvMOYK3loYce2p2RkcHDDz9ctHbt2ryGhoaMFStWlAAsW7bsnSuuuCLqrvopVXg3xuQBTwH/YK09aoxpe85aa40xtrP3WWsfBh4GmDx5cqevEfESZVZSjTIrqUaZlVSk3Ip4T2tr6/rOHp85c+bhmTNnHnZzXz43NxZPxphMnIL7L621vwk+/IExZmjw+aEQXMdBREREREREpB9JicK7cZrYHwG2Wmt/GPbUs8BXg79/FXgm0ccmIiIiIiIiEm+p0m1+KvBlYJMx5o3gY7OA+4AnjDG3AHuALybp+ERERERERETiJiUK79ba/wFMF09XJfJYRERERERERBItJbrNi4iIiIiIiKQzFd5FREREREREPE6FdxEREREREZEoLV++vNAYM2nDhg0DANasWZMzYcKEipEjR14yatSoyqVLlxaFXjtlypTyurq6rGj2o8K7iIiIiIiISJQef/zx4okTJx5/7LHHigHy8vICy5cvf2fHjh1vrVy58u1Zs2ZdcOjQIX+s+1HhXURERERERCQKjY2NvrVr1+Y9+uiju59++uligHHjxjWNHTu2CWD48OFniouLW/bt25cBUFhY2OL3+200+0qJ2eZFREREREREurKRmy84yubcnl53hkZ/E/VZ2ZQ2Z1LQ2t1rBzHm5HiWvdfda1asWFF41VVXNY4bN66pqKio5eWXX8698sorT4aeX716de6ZM2dMZWVlE8DKlSt39vYzRVLLu4iIiIiIiKSFJuqzWjnhb6I+qnHnkZ544onim266qQHghhtuOLx8+fLi0HN79uzJnDFjxoilS5fu9vtj7jWvlncRERERERFJbT21kIfs55n87cwfNor5e8+j+lgs+/zggw/8r776an5dXV3OHXfcQWtrqzHG2EAg8P6RI0d811577ch58+bVV1VVnYhlPyFqeRcREREREZG0cB7Vx/4XG+piLbgDLF++vGj69OmH9+7du6m+vn7T/v373zz//PObX3jhhbzrrrtu5I033vjhjBkzGtw4blDhXURERERERKTPnnzyyeIvfOELHQrn1dXVDbfeeutFa9euzVuxYkVJRUVFZUVFReWaNWtyYt2fus2LiIiIiIiI9NFrr722PfKx2bNnH5g9e/aBeOxPLe8iIiIiIiIiHqfCu4iIiIiIiIjHqfAuIiIiIiIi4nEqvIuIiIiIiIh4nArvIiIiIiIiIh6nwruIiIiIiIiIx6nwLiIiIiIiIhIFv98/qaKioheen8kAACAASURBVLK8vLyysrJy9IsvvjgQYM2aNTkTJkyoGDly5CWjRo2qXLp0aVHoPVOmTCmvq6vL6uu+tM67iIiIiIiISBSys7MD27Zt2wLw1FNPDZo1a9b506ZNq8vLywssX778nbFjxzbt3r0787LLLhs9ffr0oyUlJa3R7kuFdxEREREREZEYNTY2+gsKCloAxo0b1xR6fPjw4WeKi4tb9u3bl1FSUtJaWFjY4vf7bV+3r8K7iIiIiIiIpLSb7+CCzVvJ7el1jUfx1+8lq3QYzQWD6LYVfMxoTi77T97r7jVNTU2+ioqKyqamJnPo0KHM559/fnvka1avXp175swZU1lZ2QSwcuXKnT0dZ2dUeBcREREREZG0UL+XrBMnnQJ8wSBOxbq98G7zq1atGjhjxoyLtm/f/pbP50wvt2fPnswZM2aMeOSRR97x+/0x7UuFdxEREREREUlpPbWQhzzzPPnz72PY/O+xt/ozHHPzGK655poTDQ0NGfv27csoLS1tOXz4sO/aa68dOW/evPqqqqoTsW5fhXcRERERERFJC9Wf4Vj1Z6iLx7Y3bNgwIBAIMGTIkJbTp0+b6667buSNN9744YwZMxrc2H7KF96NMZ8G/h3wAz+11t6X5EMSERERERGRNBAa8w5greWhhx7anZGRwcMPP1y0du3avIaGhowVK1aUACxbtuydK664Iuqu+ildeDfG+IEfAdOA94G1xphnrbVbkntkIiIiIiIi0t+1trau7+zxmTNnHp45c+ZhN/flc3NjSTAF2GGt3WWtbQYeB6qTfEwiIiIiIiIirkr1wnspdJiY4P3gYx0YY75ujFlnjFl38ODBhB2cSLSUWUk1yqykGmVWUpFyK5LeUr3w3ivW2oettZOttZMHDx7c9vhqaimnjAIGUE4Zq6lN4lGKtOs+s8MpINszma39EwwfB5klkD0EFi9J9hFJMnSV2QdYzCCyySfTE5kN5TX7XOd28RKY+inncUkvXWUWnNwWkE0B2TzA4iQdYbvFSyB3KBSVtWc3HTJbXwuPD4dl2fCrMud+uuvuXOuVzNb+CS65HIaMgszBzk/Z2PTI7MbF8Ggu/GII/Nclyqy4z1hrY9uAMT8CVlhrX3HnkPq0778B5ltrPxW8/08A1trvd/WeyZMn23Xr1gFwNVN5lTXt28NQQAFNNDGHe/hHauJ6/JJ6jDHrrbWTE7nP7jKbEZy2woeP+SxKSmanfgrW/LX9vjHg94EFDJCXB01NcM8/Qc2dCT+8tJfszA5mEMcjVmHJIINiijmHEv4vS7iaqoQdW2RecwbAqdPB48qAvIFw/Dj4/bDwn/tnZutr4S93QssJsBYy8+Dir8HbP3OeHzIVdvwCfNnQcpy2/5kz8oL3af/dBsD4ur5/1nsNZAwMvtY654uMPDhzHEzwtdnF8L0DY0+/bzflJOpvEp5ZODu3hRRygpMMZSg/5pGEZhZg0AVw7PjZjxcWOHnFQHERlBTDkh9A1ccTenhx8+xUOLCm42PGT4dM2QD4MuHir8Cel/eTff+tZH/8T2RnFlP47DLemlFFRjY0B/NoTDd55ezsZgbzeVaWLeAL5vkkzG8e3/ye3ZiduL9O1+dag8GPs5b0MIbxY5Yl9TwbYozz4/NBXq5z7h1yLiz7z/6T2Z8NgpbIhccynBxFZq7Vd4ScO39I5qVvYE0A/7B9nLjnHgIvf7bDa3OHwf9aBqUu/xMm4/ogHjZu3Lh7/Pjxh5J9HPGwcePGkvHjxw8Pf8yNlvftwL8aY3YbY+43xlzqwjZ7ay1wsTHmImNMFnAj8Gxv3zyXBVzIhWThnGstliMc4RSnmMXd5JPJUIooIjfptZgiEMpsWdv9luB/zTS3ZXYg/oTWvC+YBWUXQIbf+VK2FlpaobXVuT3S6HxB3z3PaZ0vGh6siS/pY+tRbS1MnercSudqa+GSS2DIEMjMhIwMJsLEZB7SLOaQSRYZZGAwgJPbAxxgK1v4DNeQTyZlDGEil8S9ZT6U16xM53ZwSftzLS1OXltaoak5mNnB4D/HuS0a7rR+9rp3yeLFMGiQc+sh6+fCkS1wfA+ceNf5fcNC5/bIFqh7BFpPwZkjYFvAtjq3bffDfifQ/f2z3nsm7LVhzw1r+W8+2zqVspYnufDAUvI4NiCZf6NZzCGLrLb7RzjCGZp5lz1cxzSGUhTsATU8Ib1J5tQ4FU2Fgzo+HsprSwscOAhb6uCazwfzWubcZp/by8zW1kJZGeTmeiazkxZAXhn42v8pzsoUAQg0Obk9ve08MobvwTfoOGdy3uWD668jd8G3yfvdVEzukQ6Z6zSvnWS3uasst4bluRn8ZGZ19TkSIZRZg8Fi264P3uVdrufT5JNFHpkJyeyCWVBZDucOdipFMzLarw8CgeC59qhznn33fZg2vf3awH9Oavfiu3QO+HPAhE8J3kXmAocLaXp6Ov6L3iFj+B6yJr3OwHtm46t6kqKXpmLyjmBbnPP0H6bBT/3wSCb8vAgeyQr7PRN+mhFx39/J/az2+49kwvmMG5usv5NEL+aW97YNGVOGU3i+EcgBfgX8ylq73ZUddL3fzwD/hrNU3DJr7b909/rI2vWQB1jMQuaRTTZHONLpewsp5DRNfIm/ZwtvMZcFCa99l56tppbbuYX97OM8znO1xjnZrZghq6nlLu7kQw5xmMO00HLW+0I17xaLwccwhsa99r32T3DLt6B+b8eW9yONXb+nQ038QDh+or22edh5YTXyU6fCmjXBpn1/WJP+PVCT2r1kQueffPJ73xq9eDHMnu1cCeXlOc1vLWfnYDKwzloTnyPvXHfn2fnMJkAAoNPcAmSQSR4DE9ILKjyzmPaW99aAc6HZlYwMp+Xo+EnnnyCU3w69TAYNgmPBJpjCQmfDPh8sWpTUzCar5d139DBZgcOM4FcUshWAIrbxFncwibnkUc9pijC0cgXHWW9bE5bbWK4NQr32jnM8Ib2gQpndt98pZ4da3qHTU0CbwgLn/IoNnjKC//5+HyyaDTW/Dp5jof08Gwg4FYELFyb9PLtxMaybDbTSq5Z3X/5xMGBb/JiMVk7/8eM0v/ApTv3oW2SN3EfzWxefndfWFibbf2Eof6GI7TQwKnh7MYXsYIO5i20FMz3Z8h4SugbaS32351jn2oCE9YJavARm/4uT0VDLe+Oxrs+zoQhC+7VB6FxbXAgl53i7l0l9Lfz5ZjgZvB7qVct7aT0AGWV78J3TwOk/fhx/8RE+/Ju/wOmuOiIFAB/jeYDx/IgtzKCSR9tu/4uXOcnQTt+5iMnstusSen0QD+nW8u5a4b3DRp3W92XAOGut3/UdxKCrL+hw4ReYeeSd9YXtw0eAQFutfCGFDCCHgQxMeBfQdBf6t2qhBR8+8sijkUYs7bm+nCtYjTujOrxSeI/kfFnfzF72EiCADf7XGT9+DD7yGMhJTrpewdGZ0Jd2oLX9ohHbcwEJnC/wC8+H33ytlok1n3S+vSNlZLQXYK11vvGTXEDqygMsZh6zscHzy3GOn3WBdfmRSlZPOAH19e19i/Py4ETw6sUYp2tDV3+8zEznOWuZ1Npq11ub0PlN+prZkK4vNDPa/laQmG6gtX+Cm++Avfs7Fs67q4gKV3YBPPKxX1P10I2dvyBUmLfW2XhREZSUwJIlUJWC3yG1tXDLLU5moeP/j6H8Hum88NvGH+y+U1zM2AMHTm+yNmnd5jsTyux+9pNLbpeFeXAyGyDQ9r10mibmJroiKli51F2BPsQYGJDZyu8DN3D1h892fm7JCDYlhv5tPVAR1Z2dLGZb8FzrlPjDP1MGPjIp5x4+Qg3Uvgz/+7be/Q9+/xyomXnWw168PuhY0d9AC2e63V6o4jSRFVE33wH7D0BuTnsBPfj11SNj4Afz+8/wpp0spo55+MimJez8cuq568i44H1OPvo1Tj90BxkDM2g92kp24CDjWUI5v8JHC35aCODHR2vb7XtcxYu+5QTIwPjbhy0BzG9JfIVTPHih8O73+yddfPHFbWu3f+ELXzh877337p8yZUr58uXL3ykvL28uLS0dW19fv6kv241r4d0YkwFci9PyXgW8hNPy/owrO3BJb76gIz3AYhYwjwFkt7W8P8HjZ43jBHcLitLxi+cIR8gll+McD9YYG1pp7bSQajBkktlvW957El7zHmp57+5LO1m5PasmPqLlPfyi84op8Mo1Ea3N3RUGOlTbBy8220phXdx3oeC/mlqW/+5/c8OPj3DeB/Cjb0GrHwacgqYB8NZYODgYPhgGRR9Cwznw0VfgRB5MfRnemATzZsPVf+zFzsJ7IRwPfhsPGwbLlrUVAFMrs6HCvCGPgd0WjpKa2UUQsJ23vIdf+7dldt48yM7usnfEWUIF+0DA+fcdOPDsgrDbBf/aWrj5Zti7t/P9hN8HyMuj9fhR1l8aILsZdo6AUW/D3mHw/bnwz/fAxHVQ3ABvjIc9w+HfauDvfg7vXwDlTsM7FdvgR/+/4a9XGHZdbDAYiinmvbEHTgc2eavwHqm9MP8BueR0WhEXLo98DnI01kPts8VLYN73gxHspOW9tbX9tfl5cPRnYVnoTUkqvCKqs4rGiNwko+AfKsg30EwhcBgoBv4K/Bfw3CXn8JEtmRw6p4XCIz52XdTER3Zl8/aIJi7Ym83bFzdx8dvZvHd+M++ff4blN59gewXc/p/wxkQ4XAz/tQh7Yqv3KkrDhV9TfcjhUIeNLnObrMyGV57C2S3vBIfnQTCz73W5qZS1k8X8gu9xDgGKgALgKE521wINK2HePwOBgUx6PZ/3hp1m2P4s3ilr4qI92W23O0c4Wa4bdZzdI07yzkXwq6/A7f8B85/F7mhMbGbjwQuF99zc3EtPnjy5IfLxeBTesdbG9ANMw2ll348z3vxLwMBYtxuvn0mTJlk3/NGusqNsmR1ks+yF9lw7ypbZS22l/aNd5cr200Xo75hnM+xA67eDbJadaW+zl9pKe6E91+ZYYwdYuv3JscbmWp/Nsxn2PFtoi2yO/aG9Py7HC6yzKZrZH9r77SCbZQdav82zmfY8W2gH2Sw7yl7o2dze/+/W5pxnbdlYa1e91NkL7rc2K8vajAxrCwudW78/dKkZ209GRvt2MzOt9fns8UF+++I07IYJ2GUzsK9fiv3ll5zb313nPLb2MuzvP4PdNMbZzporsJdtaP8ZYLE/vRV71f9gf1vt3O4/t5P9G9P+ecKPI3S/rMzaVT3/u/WHzIb+386zGTbPZqR+ZnNyOuY1M9P5PdbMGtO+XZ/v7P8vwu9nZNiTJXl2x0ewy76GfWM89t3zsUcGOds6Mgj7YdHZ90O/h37eHOs89sZ453bULifjo3Zhb3zCyf+nV2GHHXQeP7eh+/N56MdMwlqvZXbVh9Zesda57YLznXahzbMZHb6XCuP4vRSrVS9Ze+EYJ7f3/3tnL1hl7YUXdjwnupFX6FNebWGhfbvcZz8sxL79Eeyym51z8U9uc865r090zsO/vCl4Tv5s2P2Jzvn5H17C/uFN7J0vO7e3/LeT0W3lWdZSbOvKs2z9sFz72hTn+F79qJPH0Ln7qv8x9njuuWH3nddctiHxmbUunWtD12EDbUbbtUGezbCDbFbqZjYhB7HK2spKa88917kO6SLLrX6f3VruZPH1S50cvjEe+/DXnfz+8qb2rIauJULXFX+/wsnnP7zi5PXbr2A/vAL71H9j35joZHTjpVjLue3n4nEdb1t851rLebbFnGs3jcG+NqU9u5UDEp/ZePy88cYbu62165L5k5OT09rZ49OmTWt4++23N1pr111yySUn+rrd4Gfr8HndmG3+j8AK4ClrbUNMG0uAaGrXJX4iZ0+H9mEJIQbDYAZ32vI+jNKEzvzr6VbMxXtg4W6YMxxqynp6df/WVRfe3rS8hzdDdWLzGBizGTaOh/EbYeM4GP+m83irP/jYeKgfBkMOnt3yXjca/v2bMGV9RtfH4WKrlKczK+1Cmd23r30gc29b3q3tXYt+F0JZ/uaPIeCDr/wMHvsaHCmEwiMd79/5b87vISUHoPq3sG4KTNgAt/3csPtif9swh9C5un3Yg2nrlhvA4gs+dwxnELGnW96nroM1R52pfu/7iM6ziyN6lYTns6eW9yjzeiIHBgY7pradg8eDv9U5L28eA62+9nMytN/fOB4ueNfpDfJhEZzTAAdLYXA9/NejMHIpvH07nMqHM8/DuDfgl1+Gv14B5W/B2+Uwe/5AztubwV+mnuSxm89w+49Sq+VdXDQ1bI6IHrQa2HpJxLVD2DWEP3jJG5nhD86FIQfa8xq6PTQcXv5nuOhR2PxdGPMDeKMaxv8atlbA6K3B221gAgMZu2kge89r5nh+M3+85mT/bnm/+R8vYPO23B7f1HjUT/3+LErPa6ZgUPcXfmMqTrLsgW77d0R2m7/rrrv23XbbbTGXiztrec/o4rW9Zq39BIBx/D0wwlq7wBhzIXCetbaTxSJSgApCCTGXBR26d/vx8/d8lb/wCh9yiGMca5+wqvYw3Bmc/3DJKKgqTu7Be83sXdBsndt0z2xVFezeHd17w7sPQ/uFaGsrJ/IM+4e20pIBGy51CjpbR0PA73QXPnAuWJ+hOdsy6urbGFnzMMs628erUX6u/kbn2XZuZPaDDyAnp+uKqYiC1KnCAewtOM7rwTViGoYN4LnPtbD8to4F7vD7v/v/nPfmkdc2Z8Yn/nkZtwcrUPvUH7AbZrN5y6VNuWfBCJj2hjM/1BydZ6mpib6CMbLg34u8kpfHviFHKf4gwOFz4PVJYA389TI4vx7OZMHeodA4KHhOLgVs8H6GM3xj8AFoyfRxjh0EGccZvC+AzTBM+muAdx+yDGmBoxOBz0MzV/BvoaE5k4LH/QdgQBkT3szh/zwE3D+Hr37LGf/+2N+b16P/Y8aZzrXuW7AA7rwTDh1yhvCFV7qGZTlgA7w9MsCeMjiTCW9McN7+1486+d0yGgqOAdZ5fu8waCxwMhwAzmloz+s5jQHI8DGgIJtL/3qCt5/yU3ie4f3rWygBTt59BV8KZrbD0l8Dyrhgrw8YwOhbFvPNKTNhOXzZeDiz8Va/P4sTJ50CfMGgUz2/oXvZ2dmBbdu2bXHj0Hri5pj3h3By9glr7WhjTBGw0lp7mSs7cEmvaykHrIYmC9kGTl8d/wOT7tUehumb4FiwcuyKQfBK4pem9HQrZiizAPerVUgcns7soD+1/z+tzEqYROe215nNXu1UkoIy2w8dopbNOLOfjWEJJZ316lv8INy90Pk9OwtO7wE8fq4N5TbLQJOuafuTPmc2KxOa3gXSfJ33Z/47n/n/Ooz539lL9afPnsSsj7oa8x6reK3zHvJRa+03gdMAwS70SV3zMiahOo0m69RYSnLN3eVc5Of4oDLXaQGRjhaG/U1m70recYj01pzh7b8rs5IKFoWdZ+/eqeuDfqaEKq7iLa7irc4LQeDMNJ+V6fze1OwUjLwuNDNds65p+5teZzY7WCRrPpMamY236k8fY8OqOjcK7onmZuH9jDHGT7DYa4wZDGEDl1ONvqC9YfEeyH0JdpxyCu3PjYO3LleX+c7UlDk9RUBf0JIalFlJNTVlTot7yBxVOqWlRd9r/33OD5J3HL2lyn1Z+N323+9eqAK8y5qamnwVFRWVoZ+ZM2eWxmtfbhbelwBPA+caY/4F+B/gXhe3n1j6gk6+xXucipNTAThwBgozVGjvib6gJdUsVEWppJiaMqf7Mah3XrpKtdZ3VZRKzUy4f077/dn3Je9Y+qHW1tb127Zt2xL6efDBB+vjtS/XCu/W2l8CdwPfB/YBn7fWPunW9pNCX9DJEyq4h1yYra7yvaEvaEk1kRWlqnSSVBDeO0+V++kpvPV94Q+Tdxy9pcp9Ce8+b0z3rxXPcnV5AGvtNmvtj6y1/2mt3ermtpNmkU52SbFwd/vvWQb2TFWre2+pJVNSjSqdJNWocl9CLZnZWdDUzDB8Q5J9SN3SeVbA6T6fnQXWej+z0inXCu/GmMnGmKeNMa8bY940xmwyxrzp1vaTRie7xFu8B5oDTjqzTccKFOmZWjIlFalVSFKNKvcl1H2++QxD8A1L9uH0SJX7kmqZlbO42fL+S+BR4Abgc8Bng7epTye7xFq422nJGOh3lunTUjx9F17ppJ5RkgpUUSqpRpkVgCsvB+AE1vuzVqtyXyC1Mts7gUAg0O+udoOf6azJ390svB+01j5rrX3HWrsn9OPi9pMn8mQX3qVb3LN4j7Pu85UFkO/vuIyU9N3CEc6FpUUXlZIawitKdZ6VVKDMysuvAjAQk5/kI+kdVe5LqmW2Z5sPHjxY0J8K8IFAwBw8eLAA2Bz5XIaL+5lnjPkpUAs0hR601v7GxX0kT00ZvNQAzx92CpfirvAJ6lY1QNPVyT2e/qCmzKlVb7bOrXowiNfpPCupJnRenb3LGe61eI/Otelmzrdh9n2Y5hQqCi8c0d7qrsymnznfhjk/wDSlUGa70dLScuv+/ft/un///jG4PJ9bEgWAzS0tLbdGPuFm4X0GUAFk0t7Eb4H+UXgHeLmx4624J7zrVr84lXhE6G8Z6tKpL2jxutD59fnDyqxbag/D3F3Oih2a+NN9NWXOjPNN1rlVZtNLzUxY+MPUKryHKvebVLmflkKZ7SeF90mTJh0Ark/2cSSKm7UTl/2/9u49zq6yvvf455eEGS4GMYKIgBMSboKXaCNFKKUkiki0lFZfek7P8doXtdhyWttJ9ZjJsSS0kmm1pR7aco7K0fZYtYqlGlAIeDkgIEIwEASTmAGNchGBCMkMZH7nj7VWZs3Onj37svZaz1rr+3695rX3rL1n77Wyf1l7/Z7n9zyPuy9193e4+7vin3dn+PrFG1mocu5+GB2LksvEGk1QlxlNAiZlkz6/ao6RbKzeBjc/CWdv1L9nvyRfYZp5vp5G3s9kk7GpQWts3JfebHgMTr89ui2DMsasANkm7zeb2UkZvl54hofgyTPVQpm1dFK5brH+fbOksW1SNppjJHsXL4q+7SdRg0i/rNXY91obvpA7efbOonejI5qvIVtJI+nqknSUlDFmBcg2eT8V2Ghm91VqqbhmkonVdAHUu3Sv+6Apce8HTVwnZTM8BOfG5d0a+9675QvgI6kGkZVbYcXG4vanihSzUjZJQ+mATc3XIN27eBGcdnB0K9JHPSfvZvYaMzPgHOA44GyqtlRcozXbYecetVRmIf1vqHL5/hgegoE5USOJet2kLBrHvktvGisa1j8GdkPzHyX23VHM9teKjc3jde4N+vfuVtJhkox9l+4tXwA3LZ2aV6RsZfRSGln0vL8d+B5wKXAWMF65peIajSxUS2VWknkEVC7fX+lxxGp0kjJQzBZn/WNKiLqhmM3e6BjMi5P09TMkQclwEDU6dccabiUb530/KqN/bUOj0zydW6U3PSfv7v4H7v4q4MPA84Arzew7ZvaXZvbrZja31/cIzvAQDM6JWir1Bd2dZOgBaB6BPKikU8pGMZud0bEoGU+W42yHxsd3TmXI2UqWkN3T5vPXPwaDN+rfvVMaWpetpELkqRnmgtuDGpukJ5mNeXf3H7j7x9z9HGAZ8P+AtwC3ZvUeQRlZCAfMgefNU0lMp5Iv5J17omV1JB8q6ZSyUcz2LjnftppTeMDAl0U/5zYsJacEvjMqQ+5d0tve2Ng0h6hxJIlVXxb9nr6S1fCwzg0PRYn7hGK2J0nczlQh0mj9Y7C/Gpukc31ZyN7dd7n7enf/I3df2o/3KNzwEBy2HzwwDu++t+i9KZd0tYLP+CzJWrqkU40mUgbpmNVFZeeSxD3t3AVTSU8ivVTUV5dEjw2mamiVDHVGS3D1ZmTb9N72QYtics+yfav0hoei7el4Bp0vOqWY7U2zKpG5TG8MHbB9G5vG1dgkncsseTezy8zs62b2NTNbZ2ZLMnrdUTP7QTyD/VVmdkjqsQ+a2ZZ4hvvXZ/F+HXnkmem3MrvRsaiccA7RF/JaTVKXm+Gh6MsDwmk02fBtOP1N0a1IIy112L1mifu6xVFyDvtOYJdu0Bsegt1nTZ0vGh+X1rQEV/dWbIwSmsS5C6JYnG1oXWM8KwntjGK2eys2Nm8kfXZZdL5NvsOSlZXU2CQ9yrLnfTMwCvwd8DDwz2b2hxm87nXAS9395cD9wAcB4jXl3wacTDTT/eW5j6//i2Oiydb+4phc37a0kovJcYeD5rb3hSzZWrsouiA3wriwWb0Obr49uhVpRuMxO9eYuCc9Ps16LZMEfbxJspNuXG32uDSn+Rq6Mzo2veQ43djUjsYEXglR+xSz3WmM2eRcm47bdMNI0suuxibpQZZj3v/R3a+Ly+X/GlgK/H4Gr/t1d382/vUW4Kj4/nnAv7r7uLv/CNgCnNLr+3VkeCgq61yzXf/p2pH+Ik2Xw0p+hoeixH3cw+hJ+603wPyDoluRZrTUYefS59oBg/EWDaXpBL3xnKBkqHuar6Ezjb2X5y7ornE/Xa2jhKgzitnONDaSnrug+bl2piqnxu36fpM2ZT7m3czea2Z/TdQD/2TGL/9u4Jr4/pHAg6nHfhxvy5fWfG/P6NhUydCgqce9SN5wW6QvXwM7n4puRWaiJbjalz7XwuxDk2YbTqNkqDuK2fb12uPeqFlPp8xOMdu+ZtVNrWJ2piqnxgRe/+7Shn5MWLceuJeoh/yv2vkDM7vezO5u8nNe6jkfAp4F/qXTHTKzC8zsdjO7/ZFHHun0z1tLZp1foFnnW0qfkNZonPts+hqzaxdNXYwXfVFz8UoYOgruvBtGLy92X6QnfY1ZlXS2L9173qxUvpnZzgkVTYYUs4FIV3R02+Oe1mo+hwroW9wm/26DBuNa6rCldMy2O39T+jmrGuYYWbc4GoarqlRpQ5YT1n3BzF7i7g+4+yeANwGXtPO37v5ad39pk59/j1/7ncAbgd9196RvqYvUbwAAIABJREFU4CfA0amXOSre1uz1r3D3pe6+9LDDDuvyCGcwPASvfA6MjcPqan1BZCq5eMnii7kG+h6zoZQhLz8DHvsF7NoNaz5a3H5Iz/oas6CSztmMjkXLDiWTfQ10UOE02zmhor1DitkArNg4vSqvlx73tNnmcyixvl8faNm41horSdudv6lVFdPwEDx5pq6PpS1Z9rx/BvicmW0ysyuBz9N6Zdm2mNk5wErgN9396dRDVwNvM7NBMzsGOA64rdf368rFi+C0g6Nb2Ve6JC65mJFihVQed8ap029FmgkpZkO0atv0Wbo7XcljtqUk1TvUOcVsa43l8t1U5Y2OwcHfbJ6cr9UM6l2xhluZLn1+7DRm17SYY0SkTVlOWHd1PCP87wI3AF8Bsli+7ePAfOA6M9toZv8Yv989RA0Em4Frgfe5+56ZX6aPli+IEvfV21Q630z6S1MXfWEIqaTz27dMvxVpJqSYDU3jOPd2y+XT2umpVO9QZ1SGPLNmk311E1fJvEMf2Lrv9ZcanLoT0gofoS0nm17KsJPqpkSFK0IkPz0n72Z2oJm9IL5/KHA6cAjRTPAP9/r67n6sux/t7kvin/emHrvE3Re7+wnuXuyMV6u3wc1PqnS+UbKu+0xLFUlxQinpHHl/NOP8yPuL2wcph1BiNjTdjHNvZqYxmdI9lSE31xiz3ZbLjyyMeogngYvu3/dxNTh1LqSYDWk52cZKkU6rm5r9XdFDF6WUek7e41L2vzOzO4AdwB8AhwEfN7Mjen390lDpfHNrtketi4Nz9OUZmlBKOocvjBL3NR/VpHXSWigxC+H0CI2O9dYTlKaZ5fsjlDLkqsbsSw6M7v9ot6ofsxJKzIaynGxWlSKgJTilZ1n0vB8NvAx4NXCcu7/c3UeAi4CP9Pr6paHS+X2NjkWlgoOmkrUQhVSGvOaj0ZJxmrROWgmpDDmEHqHGC8pue4LS1qj3PXOhlCGHELMwvdc9i5i97PioNH7XJJy9UY1OWQglZkNZTjbdWNzrUoaghlLpSRZj3l8JbHX3Pe6ejr6DyWbMe3modH5KclE54dEsxup1D1NShlz0RIKatE7aFUpJ58Ur4bSl0W1RGi8oszjPpi8qi+51q4qQYvak4+HxJ4vtfU+mZ+i11z2xfAFc9bKp8nk1OvUu+VwmvNiJ1S5eCS8OYDnZfqyYpIZS6VIWyfvBQLNvgQuAz2bw+uWh0vkpmqSuHEYWRhdQEwX3YmrSOulECCWdy8+ISjnPf1dxF5X9WoIzlF63KgklZg85GDbfX1zv++hY9G8w0Oba2O1avmBqIjD1ZGbDG26LsPwM+EXBy8n2a8UkNZRKl7JI3tcDDzTZ/ijwxQxevzyWL4DfOgzO31TvLw5NUlcew0PRfATjBa/5PvJ+GByA8QmNe5fZrVk0daFe5Lm2yOEe/VyCc7Z136VzoTSIFF0xMhIvaWhkf22gnsxsrV00lVzWeVLbdOVB1p1RoXyXSalkMWHdYzRP0v8B2NLr65dOsmxJXS94knJ5TVJXHiFMAjZ8IQzsBxPPaNy7zC7d6FTkxHVFDvdY1ccLysbX1HrEvUs3iBQZs3dsgk33Rrd5S09U14/eXPVkZiuUmC1yUtssJ1dsJpQOFCmVTNZ5d/c9ZvZWM/uKmX3OzEaA3wQOyOL1SyWERKhIKpcvn1Amriu6dV3KJYnVImO2qOEe6XXdB/twQQlaj7gfQojZoqpF+jG5YjOhVDhURShD64qK26wnV2z6Hgun7qtiRNqQxWzzK+K764C/BS4DHiaarO7Tvb5+6SSzIYdwsivCyMJo1leVy5dLCOtna8k46UQIMVvUcI/0BeWaPs6xkr5YVe9770KJ2aTKKc+YXdWwrnu/rg8KmBxwAzs4mS/BS488OZc3zFPyOY0XPNliEVVO/e51T2jmeelQFj3vo/HtFne/3t1vcvd/cvcL3f2MDF6/fEIp6czb6Fh0vCMLlbiXTSgVI1oyTtoVQswWNdwj69m6Z5LufS9y0qqqCCVmkwanPGM2KWPvV6VIs/fqJhEaHYN5N4DdACs2tvUnF3ELm3kcBuft39mblUQIky0WUeWU/j/ar173ve+l+RqkfVkk70ea2QeBe83sL83swAxes/ySHug6lY6PbIvG+6uHpnxCKZ3XknHSrlBiNu/hHv2arXsma+MJlQz1CPUqqcwbNBgvsDKviJh1ouPuZ6VIIv0enYwjXrExev6e+PdkQsgZJD3uW3gy2uA+2fnOlkAIQxGKiNk8J15W77t0IIvkfSPwS2A/4Exgu5ndYWafMrM/yeD1y2l4KErc12yvx3/Cfk9EI/0XQkmnloyTToQQs3kqYkLQUMpmqyKENd/zHqI0sm1qjoa8Ynbd4unvP5sVG/dN1pPGwSY2sIM38XU28zgTTDKfebD1kWpO0hzCmu9FxGy/VkWYSbrRqQ+VORvYwel8hQ3syPy1JV9ZzDZ/prv/vbv/vruf7u4vAM4DPg/M63kPy6xOPdGrcpjUQ/orhJJOTVonnQghZvMc6tHPJYta6aUMWfYVQhlyXnFbVMN+44SLM5XAj47B/jdOT9zPXQC+DL66pOmfJIn7LqKO9gOYy1W8Fp7YtTPLQwhKCGu+1yFm+1RNFsXsddzMw1yEOkfKLpPZ5hu5+4Pufo27j87+7AoL4WSXhzxmPpb+C6EMWZPWSSdCKEPOs8Epr7HujbopQx4dg4O/qWS/mRDKkPMaopTnuOFG6febqTpn1bapRA2i78AZknaYSoLSift/8DqW86Ks9jpMIaz5XoeY7UM12VTM7pn9yVIKfUneJVaXsYLpE10e49mkf5IvjuS2CJq0TjpRdBlyXg1OeY91T+u0DDkp79+5R6X2zYRQhpzHEKXRsahRbTCnccON0g3SEMVkugd+xcapjgeI9rFF4j7KJs7m2r1JUG0Sd5i+5ntR65HnNawu6bw4d0H+MZuuqMrg33mqSiSJ2TlchuYUKjsl7/00PBRdbI0X+AXdb3lP6iH9NbKwfpMpSfkVXYacR4NTEWMw09pd971xPe8iS8NDVnRlXh7n2VXxWHenuGuDry6Z3vC0/jE48TvRjPLpUvlZrl9G2cRKvsskUUifxCH1SdwTRQ9TyiNmR8em4qKIToxu5muYwVRjU7pK5Ox6xWxFKXnvt6K/oPut6AtKyVYIresinSq6DDmPi8oQvkvWtrGcUXr7QE6zi5dR0ZV5eVSMFN2olmjsgb9vF9MqiFv0sG5gBwv5PCv57t5tl/Jq7uG365cEFT20Lo+YTTdKFLVaVLsNpS2osanalLz3WwjjhPpFM8xXU9Gt6yqbl06lG52Kmriun5LvjsECSubTGpczGrxx+vdaugx50GD8LDXqzmR4KFoxYLzAmO3nuTaJi1AacL66pPns8S3GuCc9l2P8cu+2dbyaYV7Wr70MX9FD6/ods0UO80hrp6G0iQ3sYIjPqbGp4pS891uVezJHNMN8JRXduq6yeelGEqtFxGy/G5xCKD9OpBOx9DwD6XLTxudJc0XGLPT3XJueCK7omE2kE/gWM8qnE6Ck5/IA5ihxh+KH1vU7ZkM5z3a47nsSs6/jWh7gqb3bFbPVpOQ9D0X3ZPZLUTMfS/8VuX62ZpyXbhQZsyPvh8EBGJ/oT8yGUn4M+47JnHCwG6aPcy9ioqcyKjJm+y2kmE376pKWSfvJfIlz+dq0BOhSXs3TvFNJEBQ/QWg/hRaza9rrfR9l096k3Ukam+Yqca8wJe95KLons1/WLoL5c9XrXkVFNzipdF46VWTMDl8IA/vBxDPZx+zoWHSxPBhI+THsm8CnDVrLGbslparn2RBjtoV0r+VmHmci7pkw1HPZVJFJbp1itkXvexKzg1zJSr67ty/NSBqb3qG4rTAl73mpUgt7sn4vwJNnqoeliopucMprPddWtE51uRQds/0o50xmbp/waPhVSOfaxknAjLAufMsgaQQZsGjVliIqRvpRgjyybWrug5BiNiWZiG4/PjWt1xJgAOMkDuE6zlEC1MyaRVMTqilm+yt1PvWVW7lk9Gr250rOiatDJuKZ5A14MQcpZmtCyXteim5hz9KqbVq/tw6KbHDKaz3XVtZsj+K87P9f66Rqwz1CmPm4laQE2ZfB5DLYrQnqOpb8e40XUIbcryFKIayM0MQomxjkU+zHpziHrzHGL3kW37ubg8zhxRzEel6vCb5aKXIZ5BrF7AZ2MDR8CwdMjPGJd+7kNTfv4IX3OuNM8myqOiRJ2sd4q2K2JpS856XoXqGsjI5NtU6GMi5I+qPIBqeiJ60bHYt6wgYszKRJmiu6kTTrks7ku0JjyKutSmXIoayM0MQa7mQC59n4JzEQT0a3m3cqAWpXkcluP2LWiL7vA4rZ1dzBAzzF7v2cK35/J7e+ZoIrLtiJTcIL2H9vdYhitn6UvOepCqXz6QtilUdWW5ENTkVPWrdme9SrMBhYqbK0lpQhFzUbcpaNTunZ24talknyUZUy5JCHeQAjvJIBjHnY3uTnes5hXJPRdW5tHLNGuWMWouqBcY+OJaCYvZhX8WIOYpA5nDg2yK/eMsB7L5/PM/OGeGj0dFWH1Fhpkncz+1MzczM7NP7dzOwyM9tiZt83s1cVvY+zSvcK5V1qlIWQ1sCUfBTZ4FTkpHVFL+Ek3avKbMihl8xLdoosQ85S4DE7zMsY5108w7t4iP+s5KcXw0NR4/Z4AcsgZ924H2DJPMByXsQYb2U37+T/vPUt3PIbR/Kuz8xnrlu5v9ukZ6VI3s3saOBs4IHU5jcAx8U/FwD/UMCudWZ4aKp1PbCTRFtCWgNT8lHH0nn1eJZfUWXIWTY4qWS+XopKIBSz0q0irw+yittAS+abSle7trH2u1RXKZJ34GPASqZ/rZ0HfNojtwCHmNkRhexdJ4osNepVaGtgSv8VWYY8fCE8uSW6zVPgvUfShjWLoph18o3ZrBqc1IBUP0VdGyhmpVtFDq3LIm6TYR5lGSLXuETnyq2wYmNx+yOFCT55N7PzgJ+4+10NDx0JPJj6/cfxtrAVWWrUi2Q/B7QUUO0UWYY8ejkcfGx+4941NKQaiorZrBqc1IBUP0VdG2RVgqyYracyz+VUxphtTOCTBjOplSCSdzO73szubvJzHvDfgdU9vv4FZna7md3+yCOPZLPTvSh6RuROla11sgKCi9kqlCG3Q0NDuqaYjWXR4KTy41wEF7NFXRtkcZ5VzOYmqLhVzOYvXfGQ3EqtBJG8u/tr3f2ljT/ANuAY4C4z2w4cBdxhZi8EfgIcnXqZo+JtzV7/Cndf6u5LDzvssP4eTDuSlrP5c8vR2pfuuSrD/lZAcDFb9jLkdmloSNeCjNkiZvDu9aJS5ce5CS5miypD7vU8q5jNVVBxq5gtxleXgC+LbqV2gkjeZ+Lum9z9Be6+0N0XEpXGv8rdfwZcDbw9nnX+VOAJd/9pkfvbkeGhKBFesz3sUqP0uu6DVq7WSclOVWbwbkVDQ6qlqDLkXi8qy1jKKdkpogy519J5xWy9KWZFchV08j6L9UQ981uA/wXkPKtVBka2wc49YS8Nk07UlNDUWxG90nmWza+K13oFNVJVRfqirAyNTqNjMDEZNSBpzoV6KipmeznXlrX8WLKhmBXJVamS97gH/tH4vrv7+9x9sbu/zN1vL3r/Ohbo2pJ7qddd0ooonc+rbD4d6yqZr47hoShmIb/PtZcLyjXbNbdI3RURs9DbuTbpeS1j+bH0TjErkqtSJe+Vs3bR1AkvxNL5EfW6S0qSTEx4ftUiWc2EPJt0CZ1ivVryHvveywVl0huU97JLEpYi5mvodqWE9AodKj+uryIa93uJ2aTCSTErJaTkvUihjyNOKgIG1OsusSKqRfIonVcJXXUND0W9QeM5NTp12+BU9gmUJDt5xyx0t0pCshLNhMOAqkVqbXgoioGJHOcY6XZlj5F4iJyhmJVSUvJetKTEaMLD6n1P9mXQogoBESimWqTfpfNKmqov70anbhqcNIGSpClmpWzyXjau24b90IesisxCyXvR0iW6ec6IPJs129WaLvsqonW9n5Keo4QuQKtpbVyGbIRbOq/qD0lTzErZ5L1sXDcxq44pqQAl70VL1nxP5NFaORuNB5JWytK63tZrb5+6r9m9qyv00nlVf0ijvJc6VMxKFvJcNq6bmNUwD6kAJe8hyLu1spXk5KYZj2UmZWhdb5d6juojKZEcz2mIUieNTio/lmbSsZBHo5NiVnqlmBXpOyXvoQhl2Yr0xHk6uclM8m5d72ZG2dmo56he0iWSeVSMtNvoVODa7hsehdNvjm4lQMNDU7PO5zE+t5OGUjV8SjOKWZG+U/IeipGFMH9usQmz1nWXduXdut7trLKtjKihqlaSIUoDFiXLoTQ6FTjz8UX3wM2PR7cSqDzHvrcbs2r4lFYUsyJ9peQ9FMND8OSZxSbMq7Suu7Qp79b1rMe9j45FCRNoKcQ6ST7n8RyW52ynwSkdh0XMfBz/F/7RLvW+ByvP+RrabSRVw6e0opidRhVOkjUl7zIlWbZOve7Sjjxb17Me957+IteMs/ViDbf90k6DU7p8v4A4vOwkOGAO7JqEizbn/vbSrvR8DRse69/7tNtImuyPGj5lJnktx1aCmF19f1ThdP73lMBLNpS8S5R4HfxNWP68qHRfve7SjvRsyB/Y2t+LyqzHvevis77WLMrnPNdOg1PBYzCXHwrHHBDd/9HTurAM1tpFU1drq/vYk9lOzI6ORQ1fA1pqS1pYG59n+x0j7cYsFLY83MXHR/8UO/fA2bfB6NbZ/0akFSXvEpWP7twD1/+i+NJ9KZeRhdFZZJL+XlRCduPedfFZb3kNUZqtwSmQMZiXnTzV+/6m25XAB2l4CL6+BE47GC7u4zlrtqW3tBqNtCvP8+xsy8Wt2hbN5+QU1kh61a9MXSp94D6dZ6U3St7rLj1JXb/LSKV68rqohGzGveviU/LUqsEpkGWL0r3vuyZV2hms5QvgpqXRbT+1Os8GErMi08x2bZDXUKkWlh8KHzkh2oVJNEmo9EbJe50liUxC5fLSjbwuKrMY966JliRPrS4qA1q26LKT4aTnRMUoO/eoB77WWp1nA4pZkb1axezoWNTjPmiFX+MOL4aXPCe6r0lCpRdK3uuqMXHX2F+pOs0wL3mb6aIykJL5xPJD4Z5fh2MPjH7fNameodqaabhHYDGb0Eze0jJmV26NqksHwqi0S08SqkZS6ZaS97pKl79p7K+UQa9l85phXkLQ2HAaUAVIMv4d4N6nNLFSLc001CPQqqVkJu/V9xe9J1KoZnEbYMxqmJJkQcl7HY2OwfhkVEa0bjGMnxVEi6RIS72UzavXXYrQrMEpfUG5bnFQsbj8UPiPpdGFgQN/fh+c/C1dXNZKs5gN+Px58fFw2iHRrdRYY9wGHLMapiS9UvJeNwGWEYm0pZfl4lap110K0NjgFPAFZSKZWClJ4Df/UheXtdIsZtOVIoGdP5cfCjedFt1ueDRqbDr5m4rX2mmM24Ar7fYOUzoo+l3DlKRTSt7rJOByTZG2dLNcXHpFhcEwEyapqMYGp/RwpcAuKNOGF8PXT5kqod81Ca/T+sT10CpmA6sUabT6/qixafNTSoZqpzFu46/8UBtJYWr8O8DWp9XoJO1T8l4XmqBOqqDTce9aUUGKlm5wKtFs3UkJfXJx6cDK++DAa5XEV15JY/bi4zVnQ22lY3Z0LFqTLfD5nJJz7Py5UUHW5qdU5STtUfJeF+nW88GwT2giM+p03PuqcMcXS00kDU6rNgU5W3crycXli/efWiJ512SUxA/doIvMykpiduRvShWzTedsUG9mPaTPsyu3RtnwYPhDQ5cfClf9yvQqJyXwMhsl73UxsjBq3lu3GHZrgjopqU7GvatcXkKQNDjt+e3UtoWF7U6nlh8KY8vgulOiJD7xwO6olH6/9UrkKyeJ2WfPS21bWNjudGKfORvUm1kPe8+zv5PatrCw3elEY5WTEniZjZL3uhgegifPVAIj5dbJmPd0r7vK5aUowxfCyAbYMzf6vaQNSUkSv+6E6aX0zxIl8ufcppL6yhi+EM74N9jzxuj3kg2zazZng5Khitsbs/GHXrKYbZbAn615RmQGSt5FpDzaGfM+OgaDN6rXXcJw4odg5d1Tv5e8IWl4MTx9TpTED9jU9meZKqnf/xpddJba6Bis/zTwHuArpRxmp2SoZioUsycdFA1TmiQ6ny5UZZM0KEXybmZ/ZGY/MLN7zGxdavsHzWyLmd1nZq8vch9FJAftjHlftW0qcYfSJ0tSYqNjcN//Ba4CdlVq3oXhxTD+Brj+lGjN4hfsN/XYuEfjjQevUUl96eyd5PNpYDfM/XJpY3amZEgJfAWNbCM6z5Y/Zu85Ey49YWrb2G5Vjsh0wSfvZnYWcB7wCnc/GfjrePtJwNuAk4FzgMvNbG5hOyoi/TfbmPf0OHeoVLIkJbRqG3A+sAFOuLmSsZisWfzQ66Z6442opH7Cp8bGK2Eqib3DjX4HBg+Ev/qzQnenV82SISXwFTM6FrUYcj6wf+ljFqLG0XUnTJ8oVJUjkgg+eQf+APiIu48DuPvD8fbzgH9193F3/xGwBTiloH0Ukby0Gveu2eUlFCs2xg1Jb4TBT8IPLil6j/ou6Y2/7hQYSs1Qv3fm72+p9yhoe2MWGJwHu7e2NzloCSTJUEIz0VfEio2p5WDfCOu+UamYvS6ubFLliKSVIXk/HjjDzG41s2+a2avj7UcCD6ae9+N4m4hU2Uzj3qddeGqcuxQsWWILYPnN7U+0WAHLD4Xty6Yn8Q5s/qXKP4M2LWZ/VrmYTRJ4zURfIemYPfchWLO8UjGbVDapckTSzN1nf1a/d8LseuCFTR76EHAJcCNwEfBq4HPAIuDvgVvc/Z/j1/gEcI27/1uT178AuCD+9QTgvqyPYRaHAnX4eqjDcZ7g7vP7/SaK2Zm9iDmHH86cFz3E5I4dTD4UbXv+4Ufw/KOS5/yUn/94Bz9/qI2XC/Y4M6SYLcDxHHXsfA587k6efuIgfjZ/DsyZhMk7efbOHl86qONsy0HPnc8LFx+LWdRh4D7Jz7Zu4akndrb4q77HrWJ2utrEbHfx2A6da3NWm5gFWPCiw1lwxN7rHJ4Z382jDz7QY9zmErOSrSCS91bM7FrgUne/Mf59K3Aq8HsA7v5X8favAR929+8Uta8zMbPb3X1p0fvRb3U4zjocI+g4q6QOxwg6zqqpw3HW4RhBx1k1dTjOOhwj1Oc4q6YMZfNfBs4CMLPjgQGi1rCrgbeZ2aCZHQMcB9xW2F6KiIiIiIiI9Mm8onegDZ8EPmlmdwMTwDs8Khe4x8w+D2wmWmL2fe6+p8D9FBEREREREemL4JN3d58A/ssMj11CNCY+dFcUvQM5qcNx1uEYQcdZJXU4RtBxVk0djrMOxwg6zqqpw3HW4RihPsdZKcGPeRcRERERERGpuzKMeRcRERERERGpNSXvPTKzT5rZw/GY/GTbAjO7zsx+GN8+L95uZnaZmW0xs++b2auK2/POmNnRZnajmW02s3vM7L/F2yt1rGa2v5ndZmZ3xcf5F/H2Y8zs1vh4PmdmA/H2wfj3LfHjC4vc/3bVIW4Vs4rZEn6WilnFbNk+S8WsYrZUnyXUI27rErN1pOS9d1cC5zRs+wCwwd2PAzbEvwO8gWhW/OOI1uj8h5z2MQvPAn/q7icRLdX3PjM7ieod6ziwzN1fASwBzjGzU4FLgY+5+7HAL4D3xM9/D/CLePvH4ueVwZVUP24Vs4rZsn2WilnFbNk+S8WsYrZsnyXUI27rErP14+766fEHWAjcnfr9PuCI+P4RwH3x/X8C/lOz55XtB/h34HVVPlbgQOAO4FeJliecF29/DfC1+P7XgNfE9+fFz7Oi973N46tV3CpmFbNF738Xx6uYdcVsmX4Us4rZove/y2OudNxWPWbr9qOe9/443N1/Gt//GXB4fP9I4MHU834cbyuVuJTmlcCtVPBYzWyumW0EHgauA7YCj7v7s/FT0sey9zjjx58Anp/vHmemcp9lQjGrmKUkn2VCMauYpSSfZUIxq5ilJJ9lWpXjtsYxW2lK3vvMoyasykzpb2bPAb4I/LG7P5l+rCrH6u573H0JcBRwCnBiwbuUu6p8lqCYrYuqfJagmK2LqnyWoJiti6p8lomqx61itpqUvPfHQ2Z2BEB8+3C8/SfA0annHRVvKwUz24/oJPcv7v6leHMljxXA3R8HbiQqKzrEzObFD6WPZe9xxo8/F/h5zrualcp9lopZxWysNJ+lYlYxGyvNZ6mYVczGSvVZ1iluaxizlabkvT+uBt4R338H0ViaZPvb41krTwWeSJXnBM3MDPgEcK+7fzT1UKWO1cwOM7ND4vsHEI2BupfopPfm+GmNx5kc/5uBG+LW2jKq2mepmFXMlu2zVMwqZsv2WSpmFbOl+iyhHnFb85ittqIH3Zf9B/gs8FPgGaKxI+8hGiOyAfghcD2wIH6uAf+TaMzJJmBp0fvfwXH+GlH50PeBjfHPuVU7VuDlwJ3xcd4NrI63LwJuA7YAXwAG4+37x79viR9fVPQxKG4Vs4rZ0n6WilnFbNk+S8WsYrZUn2Vd4rYuMVvHH4s/MBEREREREREJlMrmRURERERERAKn5F1EREREREQkcEreRURERERERAKn5F1EREREREQkcEreRURERERERAKn5F1ERErHzPaY2UYzu9vMvmBmB3b493+c/hszW5+siTvD8z9sZn/Wyz63uV9Xmtmb4/v/28xOyuh1rzGzo7J4LRERESmGkncRESmjXe6+xN1fCkwA7233D81sLvDHwN7k3d3PdffHs99NMLN53fydu/+eu2/O4P0PAJ7v7j/u9bVERESkOEreRUSk7L4NHAtgZl82s++Z2T1mdkHyBDP7pZn9jZndBXwIeBFwo5ndGD++3cxg06FnAAADaklEQVQOje+/3cy+b2Z3mdlnGt/MzBab2bXx+3zbzE5s8pwPm9lnzOwm4DNmtjB+7h3xz2nx88zMPm5m95nZ9cALUq/xDTNbmux/avubzezK+P5b4uqDu8zsWzP8+/wG8I0m+/gNM7vUzG4zs/vN7Ix4+zvjf8fr4n+XPzSz95vZnWZ2i5ktmPGTEBERkb7pqjdAREQkBHGv9huAa+NN73b3x+Le5u+a2Rfd/efAQcCt7v6n8d+9GzjL3R9teL2TgVXAae7+6AyJ6hXAe939h2b2q8DlwLImzzsJ+DV33xWX6L/O3Xeb2XHAZ4GlwPnACfFzDwc2A5/s4J9gNfB6d/9Ji7L/NwBfnuGxee5+ipmdC/wP4LXx9pcCrwT2B7YAf+7urzSzjwFvB/62g30UERGRDCh5FxGRMjrAzDbG978NfCK+f5GZnR/fPxo4Dvg5sAf4Yhuvuwz4QpLUu/tj6QfN7DnAacAXzCzZPDjDa13t7rvi+/sBHzezJfG+HB9v/3Xgs+6+B9hhZje0sY9pNwFXmtnngS/N8JzTgZnG6yd/8z1gYWr7je6+E9hpZk8A/xFv3wS8vMN9FBERkQwoeRcRkTLa5e5L0hvM7DeIeo5f4+5Pm9k3iHqOAXbHCXKv5gCPN773DJ5K3f8T4CHgFfFr7O7wfT11f/+9G93fG/f+rwC+Z2a/ElcaAGBmi4AH3X1ihtcdj2/3MP2aYDx1fzL1+yS6dhARESmExryLiEhVPBf4RZy4nwic2uK5O4H5TbbfALzFzJ4P0Fg27+5PAj8ys7fEj5uZvaLNffupu08C/xWYG2//FvBWM5trZkcAZ83w9w+Z2UvMbA5RqT3x+y9291vdfTXwCFG1QVp6SIGIiIiUmJJ3ERGpimuBeWZ2L/AR4JYWz70CuDaZsC7h7vcAlwDfjCe3+2iTv/1d4D3x4/cA57Wxb5cD74j/5kSmeuWvAn5INNb908B3Zvj7DwBfAW4GfpraPmpmm8zs7vixuxr+7hyUvIuIiFSCufvszxIREZFSMbNB4CZ3X1r0voiIiEjvlLyLiIiIiIiIBE5l8yIiIiIiIiKBU/IuIiIiIiIiEjgl7yIiIiIiIiKBU/IuIiIiIiIiEjgl7yIiIiIiIiKBU/IuIiIiIiIiEjgl7yIiIiIiIiKB+/9g2+lwIbCO1wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure(figsize=(15,6))\n", + "axesR = {}\n", + "axesI = {}\n", + "for i, lMax in enumerate(lMaxes):\n", + " axesR[lMax] = fig.add_subplot(2,len(lMaxes),i+1)\n", + " axesR[lMax].set_xlim([50,300])\n", + " axesR[lMax].set_ylim([1.25, ef[1]/eh])\n", + " axesI[lMax] = fig.add_subplot(2,len(lMaxes),len(lMaxes)+i+1)\n", + " axesI[lMax].set_xlim([50,300])\n", + " axesI[lMax].set_ylim([-60, 30])\n", + " axesR[lMax].set_title('$l_\\max = %d $' % lMax) \n", + " axesR[lMax].tick_params(labelbottom=False) \n", + " if i == len(lMaxes)//2:\n", + " axesI[lMax].set_xlabel(\"Particle radius / nm\")\n", + " if i == 0:\n", + " axesR[lMax].set_ylabel('$\\hbar \\Re \\omega / \\mathrm{eV}$')\n", + " axesI[lMax].set_ylabel('$\\hbar \\Im \\omega / \\mathrm{meV}$')\n", + " else:\n", + " axesR[lMax].tick_params(labelleft=False) \n", + " axesI[lMax].tick_params(labelleft=False) \n", + "\n", + "res_thr = 0.00015\n", + "\n", + "ir_labeled=set()\n", + "if True:\n", + " for (lMax, radius), (eigvals, residuals, b, ef, irclass) in plotdata.items():\n", + " for i, (e, res, iri) in enumerate(zip(eigvals, residuals, irclass)):\n", + " #if i == 0:\n", + " if res < res_thr:# and e.real < 2.14e15:\n", + " if iri in ir_labeled: \n", + " label=None\n", + " else:\n", + " ir_labeled.add(iri)\n", + " label=irlabels[iri]\n", + " axesR[lMax].plot(radius*1e9, e.real/eh, \n", + " marker='.',\n", + " #marker=markerfun(b),\n", + " ms=3, #c=colorfun(b)\n", + " c=matplotlib.cm.hsv(iri/9),\n", + " #c = colorfun(iri),\n", + " label=label,\n", + " )\n", + " axesI[lMax].plot(radius*1e9, e.imag/eh*1000, \n", + " #marker='x', \n", + " #c=colorfun(b), \n", + " c=matplotlib.cm.hsv(iri/9),#colorfun(iri),\n", + " marker='.', #markerfun(b),\n", + " ms=3,\n", + " #label=label\n", + " )\n", + "fig.legend(title=\"Irrep\", loc=\"center right\")\n", + "#fig.suptitle('$l_\\mathrm{max}=%d$, residual threshold = %g' % (lMax, res_thr) )\n", + "fig.savefig(plotfilename)\n", + "fig.savefig(plotfilename.replace('pdf', 'png'))\n", + "print(plotfilename)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0. , 1.40635433, 1.98888536, 2.81270865, 3.14470387,\n", + " 3.97777072, 4.21906298, 4.44728287])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ef / eh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3/anaconda", + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/size_vs_lMax/AuSphere/parallel_generate.sh b/examples/size_vs_lMax/AuSphere/parallel_generate.sh new file mode 100644 index 0000000..de3a9ca --- /dev/null +++ b/examples/size_vs_lMax/AuSphere/parallel_generate.sh @@ -0,0 +1,19 @@ +#!/bin/bash +#SBATCH --mem=200 +#SBATCH -t 30:00 +#SBATCH -c 4 +#SBATCH -p short-ivb +#SBATCH --array=0-250 + +cat $0 + +contour_points=410 + +#radii_nm=(`seq 80 1 150`) +radii_nm=(`seq 50 1 300`) +radius_nm=${radii_nm[$SLURM_ARRAY_TASK_ID]} + +for lMax in $(seq 1 5); do + srun rectlat_simple_modes.py -p 580e-9 -m '4+0.7j' -r ${radius_nm}e-9 -k 0 0 --kpi -n 1.52 -L lMax -t 1e11 -b -2 -f 0.1 -i 1. -T .3 -N ${contour_points} +done + diff --git a/examples/size_vs_lMax/AuSphere/serial.sh b/examples/size_vs_lMax/AuSphere/serial.sh new file mode 100644 index 0000000..c756339 --- /dev/null +++ b/examples/size_vs_lMax/AuSphere/serial.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +kx=0.0 +contour_points=410 + +radii_nm=(`seq 50 1 300`) +radius_nm=${radii_nm[$SLURM_ARRAY_TASK_ID]} + +for lMax in $(seq 1 5) ; do + for radius_nm in $(seq 50 1 300) ; do + rectlat_simple_modes.py -p 580e-9 -m 'Au' -r ${radius_nm}e-9 -k $kx 0 --kpi -n 1.52 -L $lMax -t 1e11 -b -2 -f 0.1 -i 1. -T .3 -N ${contour_points} --lMax-extend 10 + done +done diff --git a/examples/size_vs_lMax/AuSphere/serial1.sh b/examples/size_vs_lMax/AuSphere/serial1.sh new file mode 100644 index 0000000..c7086da --- /dev/null +++ b/examples/size_vs_lMax/AuSphere/serial1.sh @@ -0,0 +1,11 @@ +#!/bin/bash + +kx=0.0 +contour_points=410 + +radii_nm=(`seq 50 1 300`) +radius_nm=${radii_nm[$SLURM_ARRAY_TASK_ID]} + + for radius_nm in $(seq 50 1 300) ; do + rectlat_simple_modes.py -p 580e-9 -m 'Au' -r ${radius_nm}e-9 -k $kx 0 --kpi -n 1.52 -L 1 -t 1e11 -b -2 -f 0.1 -i 1. -T .3 -N ${contour_points} + done diff --git a/examples/size_vs_lMax/_mode_r_dep.ipynb b/examples/size_vs_lMax/_mode_r_dep.ipynb new file mode 100644 index 0000000..08adb48 --- /dev/null +++ b/examples/size_vs_lMax/_mode_r_dep.ipynb @@ -0,0 +1,300 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import re\n", + "import numpy as np\n", + "import matplotlib\n", + "from matplotlib import pyplot as plt\n", + "from scipy.constants import hbar, e as eV, c\n", + "eh = eV/hbar\n", + "import glob\n", + "def ri(z): return (z.real, z.imag)\n", + "#m = re.compile(r\"([^_]+)_r([0-9.]+)nm_\")\n", + "#removek = re.compile(r\"(k\\([^)]+\\)um-1_)\")\n", + "remover = re.compile(r\"r[0-9.]+nm_\")\n", + "\n", + "\n", + "markerdict = {\n", + " 4: \"3\",\n", + " -4: \"4\",\n", + " 3: \"^\",\n", + " -3: \"v\",\n", + " -2: 'x',\n", + " 2: '+',\n", + " 1: 's',\n", + " -1: 'd',\n", + "}\n", + "\n", + "prop_cycle = plt.rcParams['axes.prop_cycle']\n", + "colors = prop_cycle.by_key()['color']\n", + "colordict = {i: colors[(i+1)] for i in range(-4,8)}\n", + "\n", + "def markerfun(b):\n", + " if b in markerdict.keys():\n", + " return markerdict[b]\n", + " else: return 'X'\n", + "\n", + "def colorfun(b):\n", + " if (b+1) in colordict.keys():\n", + " return colordict[b+1]\n", + " else: return colordict[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "allfiles=glob.glob('*sph*k(0_0)*.npz')\n", + "allgraphs=dict()\n", + "for f in allfiles:\n", + " base = remover.sub('', f)\n", + " if base in allgraphs.keys():\n", + " allgraphs[base] += 1\n", + " else:\n", + " allgraphs[base] = 1\n", + "for k in sorted(allgraphs.keys()):\n", + " print(k, allgraphs[k])" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "ename": "FileNotFoundError", + "evalue": "[Errno 2] No such file or directory: 'projectors_D4h_lMax1.npz'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mlMaxes\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mlMax\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mlMax\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mlMax\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlMaxes\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mproj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'projectors_D4h_lMax%d.npz'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mlMax\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m \u001b[0mirlabels\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mproj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mproj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mproj\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mf\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mirlabels\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/.local/lib/python3.7/site-packages/numpy-1.17.3-py3.7-linux-x86_64.egg/numpy/lib/npyio.py\u001b[0m in \u001b[0;36mload\u001b[0;34m(file, mmap_mode, allow_pickle, fix_imports, encoding)\u001b[0m\n\u001b[1;32m 426\u001b[0m \u001b[0mown_fid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 427\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 428\u001b[0;31m \u001b[0mfid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mos_fspath\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"rb\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 429\u001b[0m \u001b[0mown_fid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 430\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'projectors_D4h_lMax1.npz'" + ] + } + ], + "source": [ + "projectors = dict()\n", + "projectors_list = dict()\n", + "lMaxes = [lMax for lMax in range(1,6)]\n", + "for lMax in lMaxes:\n", + " proj = np.load('projectors_D4h_lMax%d.npz' % lMax)\n", + " irlabels = sorted(proj.keys())\n", + " proj = {f: proj[f] for f in irlabels}\n", + " proj_list = [proj[irlabels[i]] for i in range(len(proj))]\n", + " projectors[lMax] = proj\n", + " projectors_list[lMax] = proj_list\n", + "globpattern = '*sph_r*_p580nmx580nm_mAu_n1.52_b?2_k(0_0)um-1_L?_cn???.npz'\n", + "filenames=glob.glob(globpattern)\n", + "plotfilename = 'collected_' + globpattern.replace('*', 'XXX').replace('?', 'X').replace('npz','pdf')\n", + "print(filenames[:4], plotfilename)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [], + "source": [ + "#projectors\n", + "#glob.glob('cyl_r100nm*L3*3100.npz')\n", + "#glob.glob('sph_r100*m5*.npz')\n", + "#dat['meta'][()],list(dat.keys())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "inpure result detected [1. 0.99999999 1. 0.97991334 0.99999996 0.9999989\n", + " 0.99999983 0.99999966 0.99999322 0.99999721 0.99999653] [3.28735741e-04 2.66532534e-05 2.47011478e-05 1.45012420e-01\n", + " 2.44785416e-04 7.05405359e-04 1.60203586e-03 1.71245137e-03\n", + " 1.03244480e-02 9.18732728e-03 1.18651583e-02]\n", + "inpure result detected [1. 1. 0.99999998 0.99999999 0.99999996 0.96608887\n", + " 0.99999852 0.99999397 0.99998951 0.99999912 0.99982435] [2.66223026e-04 2.12357147e-05 3.54211968e-05 1.06651057e-04\n", + " 2.79595790e-04 2.41939163e-01 2.17645058e-03 3.41541473e-03\n", + " 1.14507609e-02 1.49639498e-02 2.33483138e-02]\n", + "inpure result detected [1. 1. 0.92521572 1. 0.99999627 0.99990293\n", + " 0.99946049] [1.59712906e-05 3.60193407e-05 2.48341492e-01 1.21848930e-03\n", + " 3.81805601e-03 2.42649228e-02 2.99534246e-02]\n", + "inpure result detected [1. 1. 0.99999998 0.99999961 0.93267685 0.99999964\n", + " 0.99999822 0.99921774 0.99995547 0.99997301] [5.22490396e-04 3.01556792e-05 4.88795563e-05 6.29703960e-04\n", + " 2.34414238e-01 3.72766210e-03 4.72444059e-03 7.62106094e-02\n", + " 6.32796684e-02 5.63231562e-02]\n" + ] + } + ], + "source": [ + "plotdata = {}\n", + "for file in filenames:\n", + " dat = np.load(file, allow_pickle=True)\n", + " kx = dat['meta'][()]['k'][0]\n", + " radius = dat['meta'][()]['radius']\n", + " b = dat['meta'][()]['band_index']\n", + " eigvals = dat['eigval']\n", + " lMax = dat['meta'][()]['lMax']\n", + " residuals = dat['residuals']\n", + " ef =dat['empty_freqs']\n", + " eigvecs = dat['eigvec']\n", + " irweights = []\n", + " #for proj in projectors_list[lMax]:\n", + " # try:\n", + " # irweights.append(np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)))\n", + " # except ValueError as err:\n", + " # print(proj, len(proj))\n", + " # raise err\n", + " irweights = np.array(irweights)\n", + " #print(irweights)\n", + " irweights = np.array([np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)) for proj in projectors_list[lMax]]).T\n", + " irclass = np.argmax(irweights, axis=-1)\n", + " purities = np.amax(irweights, axis=-1)\n", + " if (np.any(purities < 0.98)):\n", + " print(\"inpure result detected\", purities, residuals)\n", + " #print(purities)\n", + " \n", + " #for i in range(len(residuals)): \n", + " # if residuals[i] < 0.01:\n", + " # vec = eigvecs[i]\n", + " # for irlabel, proj in projectors.items():\n", + " # print(irlabel, np.linalg.norm(np.dot(proj, vec))) #maybe some conj() here?\n", + " # print('--->', irlabels[irclass[i]])\n", + "\n", + " \n", + " plotdata[(lMax,radius)] = (eigvals, residuals, b, ef, irclass,)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(15,6))\n", + "axesR = {}\n", + "axesI = {}\n", + "for i, lMax in enumerate(lMaxes):\n", + " axesR[lMax] = fig.add_subplot(2,len(lMaxes),i+1)\n", + " axesR[lMax].set_xlim([50,300])\n", + " axesR[lMax].set_ylim([1.25, ef[1]/eh])\n", + " axesI[lMax] = fig.add_subplot(2,len(lMaxes),len(lMaxes)+i+1)\n", + " axesI[lMax].set_xlim([50,300])\n", + " axesI[lMax].set_ylim([-60, 30])\n", + " axesR[lMax].set_title('$l_\\max = %d $' % lMax) \n", + " axesR[lMax].tick_params(labelbottom=False) \n", + " if i == len(lMaxes)//2:\n", + " axesI[lMax].set_xlabel(\"Particle base radius / nm\")\n", + " if i == 0:\n", + " axesR[lMax].set_ylabel('$\\hbar \\Re \\omega / \\mathrm{eV}$')\n", + " axesI[lMax].set_ylabel('$\\hbar \\Im \\omega / \\mathrm{meV}$')\n", + " else:\n", + " axesR[lMax].tick_params(labelleft=False) \n", + " axesI[lMax].tick_params(labelleft=False) \n", + "\n", + "res_thr = 0.005\n", + "\n", + "ir_labeled=set()\n", + "if True:\n", + " for (lMax, radius), (eigvals, residuals, b, ef, irclass) in plotdata.items():\n", + " for i, (e, res, iri) in enumerate(zip(eigvals, residuals, irclass)):\n", + " #if i == 0:\n", + " if res < res_thr:# and e.real < 2.14e15:\n", + " if iri in ir_labeled: \n", + " label=None\n", + " else:\n", + " ir_labeled.add(iri)\n", + " label=irlabels[iri]\n", + " axesR[lMax].plot(radius*1e9, e.real/eh, \n", + " marker='.',\n", + " #marker=markerfun(b),\n", + " ms=4, #c=colorfun(b)\n", + " c=matplotlib.cm.hsv(iri/9),\n", + " #c = colorfun(iri),\n", + " label=label,\n", + " )\n", + " axesI[lMax].plot(radius*1e9, e.imag/eh*1000, \n", + " #marker='x', \n", + " #c=colorfun(b), \n", + " c=matplotlib.cm.hsv(iri/9),#colorfun(iri),\n", + " marker='.', #markerfun(b),\n", + " ms=4,\n", + " #label=label\n", + " )\n", + "fig.legend(title=\"Irrep\", loc=\"center right\")\n", + "#fig.suptitle('$l_\\mathrm{max}=%d$, residual threshold = %g' % (lMax, res_thr) )\n", + "fig.savefig(plotfilename)\n", + "fig.savefig(plotfilename.replace('pdf', 'png'))\n", + "print(plotfilename)" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0. , 1.40635433, 1.98888536, 2.81270865, 3.14470387,\n", + " 3.97777072, 4.21906298, 4.44728287])" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ef / eh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}