diff --git a/docs/source/Examples/Reading_FoF.ipynb b/docs/source/Examples/Reading_FoF.ipynb new file mode 100644 index 0000000..40e427e --- /dev/null +++ b/docs/source/Examples/Reading_FoF.ipynb @@ -0,0 +1 @@ +{"cells":[{"metadata":{},"id":"e5d755f9","cell_type":"markdown","source":"# Reading FoF halo catalogs\n\n[![Binder](https://mybinder.org/badge_logo.svg)](https://binder.flatironinstitute.org/v2/user/fvillaescusa/Quijote?filepath=/Tutorials/Reading_FoF.ipynb)"},{"metadata":{"trusted":true},"id":"8ef62591","cell_type":"code","source":"import numpy as np\nimport readgadget\nimport readfof\nimport redshift_space_library as RSL","execution_count":1,"outputs":[]},{"metadata":{"trusted":true},"id":"88fd93b8","cell_type":"code","source":"snapdir = '/home/jovyan/Data/Halos/FoF/EQ_p/420' #folder hosting the catalogue\nsnapnum = 4 #number of the catalog (4-->z=0, 3-->z=0.5, 2-->z=1, 1-->z=2, 0-->z=3)","execution_count":2,"outputs":[]},{"metadata":{},"id":"7448f84a","cell_type":"markdown","source":"In Quijote, snapnum={4,3,2,1,0} corresponds to redshift {0, 0.5, 1, 2, 3}, but we recommend reading it directly\nfrom the header of the corresponding snapshot. The header of the snapshot also contains information that is needed\nfor some particular operations such as move halos to redshift-space."},{"metadata":{"trusted":true},"id":"1495f131","cell_type":"code","source":"# get the name of the corresponding snapshot\nsnapshot = '/home/jovyan/Data/Snapshots/EQ_p/420/snapdir_%03d/snap_%03d'%(snapnum,snapnum)\n\n# read the redshift, boxsize, cosmology...etc in the header\nheader = readgadget.header(snapshot)\nBoxSize = header.boxsize/1e3 #Mpc/h\nNall = header.nall #Total number of particles\nMasses = header.massarr*1e10 #Masses of the particles in Msun/h\nOmega_m = header.omega_m #value of Omega_m\nOmega_l = header.omega_l #value of Omega_l\nh = header.hubble #value of h\nredshift = header.redshift #redshift of the snapshot\nHubble = 100.0*np.sqrt(Omega_m*(1.0+redshift)**3+Omega_l)#Value of H(z) in km/s/(Mpc/h)\n\nprint('BoxSize = %.3f Mpc/h'%BoxSize)\nprint('Number of particles in the snapshot:',Nall)\nprint('Omega_m = %.3f'%Omega_m)\nprint('Omega_l = %.3f'%Omega_l)\nprint('h = %.3f'%h)\nprint('redshift = %.1f'%redshift)","execution_count":3,"outputs":[{"output_type":"stream","text":"BoxSize = 1000.000 Mpc/h\nNumber of particles in the snapshot: [ 0 134217728 0 0 0 0]\nOmega_m = 0.318\nOmega_l = 0.682\nh = 0.671\nredshift = 0.0\n","name":"stdout"}]},{"metadata":{},"id":"97a1d704","cell_type":"markdown","source":"To read the halo catalog we do"},{"metadata":{"trusted":true},"id":"17e31f7d","cell_type":"code","source":"# read the halo catalogue\nFoF = readfof.FoF_catalog(snapdir, snapnum, long_ids=False,\n swap=False, SFR=False, read_IDs=False)\n\n# get the properties of the halos\npos_h = FoF.GroupPos/1e3 #Halo positions in Mpc/h\nvel_h = FoF.GroupVel*(1.0+redshift) #Halo peculiar velocities in km/s\nmass_h = FoF.GroupMass*1e10 #Halo masses in Msun/h\nNp_h = FoF.GroupLen #Number of CDM particles in the halo. Even in simulations with massive neutrinos, this will be just the number of CDM particles","execution_count":4,"outputs":[]},{"metadata":{},"id":"20069f3c","cell_type":"markdown","source":"Lets print some information"},{"metadata":{"trusted":true},"id":"e94f40bf","cell_type":"code","source":"print('%.3f < X_h < %.3f Mpc/h'%(np.min(pos_h[:,0]), np.max(pos_h[:,0])))\nprint('%.3f < Y_h < %.3f Mpc/h'%(np.min(pos_h[:,1]), np.max(pos_h[:,1])))\nprint('%.3f < Z_h < %.3f Mpc/h'%(np.min(pos_h[:,2]), np.max(pos_h[:,2])))\nprint('%.3f < Vx_h < %.3f km/s'%(np.min(vel_h[:,0]), np.max(vel_h[:,0])))\nprint('%.3f < Vy_h < %.3f km/s'%(np.min(vel_h[:,1]), np.max(vel_h[:,1])))\nprint('%.3f < Vz_h < %.3f km/s'%(np.min(vel_h[:,2]), np.max(vel_h[:,2])))\nprint('%.3e < M_h < %.3e Msun/h'%(np.min(mass_h), np.max(mass_h)))\nprint('%d < Np < %d'%(np.min(Np_h), np.max(Np_h)))","execution_count":5,"outputs":[{"output_type":"stream","text":"0.003 < X_h < 999.997 Mpc/h\n0.004 < Y_h < 999.997 Mpc/h\n0.001 < Z_h < 1000.000 Mpc/h\n-1920.503 < Vx_h < 2176.723 km/s\n-1891.656 < Vy_h < 2023.595 km/s\n-1776.037 < Vz_h < 1855.787 km/s\n1.313e+13 < M_h < 4.270e+15 Msun/h\n20 < Np < 6504\n","name":"stdout"}]},{"metadata":{},"id":"d3c2ffb7","cell_type":"markdown","source":"By construction, we only keep halos that contain at least 20 dark matter particles. We can verify that the minimum mass of a halo corresponds to that"},{"metadata":{"trusted":true},"id":"0d4ddbc8","cell_type":"code","source":"Minimum_mass = 20*Masses[1] #This is 20 times the mass of a single DM particle\nprint('%.3e should be equal to\\n%.3e'%(Minimum_mass, np.min(mass_h)))","execution_count":6,"outputs":[{"output_type":"stream","text":"1.313e+13 should be equal to\n1.313e+13\n","name":"stdout"}]},{"metadata":{},"id":"53e67704","cell_type":"markdown","source":"To get the information about a particular halo, just select its index"},{"metadata":{"trusted":true},"id":"45198628","cell_type":"code","source":"index = 45 #the index of the halo\nprint('Halo position:',pos_h[index],'Mpc/h')\nprint('Halo velocity:',vel_h[index],'km/s')\nprint('Halo mass: %.3e Msun/h'%mass_h[index])\nprint('Number of particles in the halo: %d'%Np_h[index])","execution_count":7,"outputs":[{"output_type":"stream","text":"Halo position: [557.547 854.5584 387.07922] Mpc/h\nHalo velocity: [ 231.89758 -459.18256 655.8717 ] km/s\nHalo mass: 1.910e+15 Msun/h\nNumber of particles in the halo: 2909\n","name":"stdout"}]},{"metadata":{"trusted":false},"cell_type":"markdown","source":"In some cases, we may want to work with halos in redshift-space rather than real-space.\nWe can move halos to redshift-space along a given axis as follows"},{"metadata":{"trusted":true},"cell_type":"code","source":"# move halos to redshift-space. After this call, pos_h will contain the\n# positions of the halos in redshift-space\naxis = 0 #axis along which to displace halos\nRSL.pos_redshift_space(pos_h, vel_h, BoxSize, Hubble, redshift, axis)","execution_count":8,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets finally compute a simple summary statistics from the catalog: the halo mass function"},{"metadata":{"trusted":true},"cell_type":"code","source":"min_mass = 2e13 #minimum mass in Msun/h\nmax_mass = 1e15 #maximum mass in Msun/h\nbins = 30 #number of bins in the HMF\n\n# Correct the masses of the FoF halos\nmass_h = mass_h*(1.0-Np_h**(-0.6))\n\nbins_mass = np.logspace(np.log10(min_mass), np.log10(max_mass), bins+1)\nmass_mean = 10**(0.5*(np.log10(bins_mass[1:])+np.log10(bins_mass[:-1])))\ndM = bins_mass[1:] - bins_mass[:-1]\n\n# compute the halo mass function (number of halos per unit volume per unit mass)\nHMF = np.histogram(mass_h, bins=bins_mass)[0]/(dM*BoxSize**3)","execution_count":9,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nplt.xscale('log')\nplt.yscale('log')\nplt.xlabel(r'$M_{\\rm halo}~[h^{-1}M_\\odot]$')\nplt.ylabel(r'$HMF~[h^4M_\\odot^{-1}{\\rm Mpc}^{-3}]$')\nplt.plot(mass_mean, HMF)\nplt.show()","execution_count":10,"outputs":[{"output_type":"display_data","data":{"text/plain":"
","image/png":"iVBORw0KGgoAAAANSUhEUgAAAlEAAAG7CAYAAADjS8CnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSjUlEQVR4nO3deVzUdeLH8dfMcCmXBwoi4H2EyCEi4pFalFl5pt2JVlZmZbluaVu5ne7WrrmtY1aWZlmpZWZpx2qZpSYKgnhrauIB3lzKNczvjzY2fx7JcHxn4P18PHhsfGfm+33D1sybz/fz/XxNdrvdjoiIiIhUiNnoACIiIiKuSCVKRERExAEqUSIiIiIOUIkSERERcYBKlIiIiIgDVKJEREREHKASJSIiIuIAN6MD1FZlZWUcPnwYX19fTCaT0XFERETkMtjtdvLy8ggODsZsvvRYk0pUNTl8+DChoaFGxxAREREHZGZmEhIScsnnqERVE19fX+DX/xP8/PwMTiMiIiKXIzc3l9DQ0PLP8UtRiaomv53C8/PzU4kSERFxMZczFUcTy0VEREQcoBIlIiIi4gCVKBEREREHqESJiIiIOEAlSkRERMQBKlEiIiIiDlCJEhEREXGASpSIiIiIA1SiRERERBygEiUiIiLiAJUoEREREQeoRImIiIg4QCUKGDp0KA0bNmT48OHnbN+5cyfR0dHlX/Xq1WPJkiXGhPwvW5mdLYdyDM0gIiIiKlEAjB8/nnnz5p23vUOHDqSlpZGWlsaPP/6It7c311xzjQEJf2W323lqSQZDZ67hi82HDcshIiIiKlEA9O3bF19f30s+Z+nSpVx99dV4e3vXUKrzlZbZyT1bSonNzsMfbuLdtfsNyyIiIlLXOX2JWr16NQMHDiQ4OBiTyXTB02lWq5WWLVvi5eVFfHw8ycnJVZ5j4cKF3HLLLVW+34pwt5h57bYY7ureArsdpizdyj++3ondbjc0l4iISF3k9CWqoKCAqKgorFbrBR9fsGABEyZMYMqUKaSmphIVFUX//v05evRo+XOio6OJiIg47+vw4cs7JZabm8vatWu5/vrrL/qcoqIicnNzz/mqDhaziecGd2LCNe0BmPHdHiZ9kkGpraxajiciIiIX5mZ0gD8yYMAABgwYcNHHp02bxpgxYxg9ejQAs2bNYtmyZbzzzjtMmjQJgLS0tEpl+Oyzz7j22mvx8vK66HOmTp3Ks88+W6njXC6TycQjV7ejia8nf/k0gwUbMzl5pph/3xaDl7ulRjKIiIjUdU4/EnUpxcXFpKSkkJiYWL7NbDaTmJjIunXrquw4l3Mqb/LkyeTk5JR/ZWZmVtnxL+a2bmHMvCMWDzcz/9mWzV1vryfnTEm1H1dERERcvEQdP34cm81GYGDgOdsDAwPJysq67P0kJiYyYsQIli9fTkhIyDkFLCcnh+TkZPr373/JfXh6euLn53fOV024LiKI9+7uhq+XGxv2n+LmN9aRnVtYI8cWERGpy1y6RFWVFStWcOzYMc6cOcPBgwdJSEgof8zf35/s7Gw8PDwMTHhp8a0bs/D+BJr4erIzO49hM9fy87F8o2OJiIjUai5dogICArBYLGRnZ5+zPTs7m6CgIEMyWa1WwsPDiYuLq9HjXtHMj8Vje9AqwJtDp88y/PW1pGWertEMIiIidYlLlygPDw9iY2NZuXJl+baysjJWrlx5zmhSTRo3bhzbtm1jw4YNNX7s0Eb1+fiBBCJD/Dl1poTb3/qJ73cdq/EcIiIidYHTl6j8/PzyVcMB9u3bR1paGgcOHABgwoQJvPXWW7z77rts376dsWPHUlBQUH61Xl3T2MeTD8d0p3e7AM4U27hn7gaWbDpkdCwREZFax2R38pUaV61aRb9+/c7bnpSUxNy5cwGYMWMGr7zyCllZWURHR/Paa68RHx9fw0nPlZubi7+/Pzk5OTU2yfz3ikvLmLgonaXpv66F9fSN4dzTq1WN5xAREXElFfn8dvoS5aqMLlEAZWV2nl+2jTlr9gPwQJ82PHFdB0wmkyF5REREnF1FPr+d/nSeqzFqYvmFmM0mnrkxnMev6wDArO9/ZuKizZRodXMREZFK00hUNXGGkajfW7gxk8mLM7CV2enboQmv3RaDn5e70bFEREScikai5Dw3dw3ljTtj8XI3s2rnMYbMWMOeo1pLSkRExFEqUXVIYnggi+7vQbC/F3uPFzDEuoYV27L/+IUiIiJyHpWoOqZziD9LH+5Ft1aNyC8q5d55G/nXit2UlemsroiISEWoRFUxZ5pYfjEBPp7MvzeepIQWALy6YhcPvJ9CflGpwclERERchyaWVxNnm1h+MQs3ZvLUp1sotpXRrqkPb47sSqsAb6NjiYiIGEITy+Wy3dw1lAX3dyfQz5PdR/MZNONHvtt51OhYIiIiTk8lSogJa8jnD/citkVD8gpLuXvuBmau2oMGKUVERC5OJUoAaOrrxYdjunNbtzDsdnj5q5089OEmzhRrnpSIiMiFqERVMVeYWH4xHm5mpg7rzItDI3C3mFi2+QjDZq4l8+QZo6OJiIg4HU0sryauMrH8YjbsP8nY91M5nl9Eg/ruWG/vQs+2AUbHEhERqVaaWC6VFteyEZ8/3JOoEH9OnynhrrfXM/uHvZonJSIi8l8aiaomrj4S9ZvCEht/+XQLn6QeBGBQVDDXhAfi4+WGr6cb3p5u+Hi64ev16z+7W9TLRUTEdVXk81slqprUlhIFYLfbeXftfp5fth3bH6xs7uVuxsfT/b+lyoKPp1v59z6ebiS0acz1nZvVUHIREZGKUYlyArWpRP1m/d4TzF27n1NniskvKiW/sJT8olLyCkspKi277P2Mv7odjya2w2QyVWNaERGRiqvI57dbDWWSWiC+dWPiWze+4GMltjIK/luo8otKy0tWXlEpBf/9573H8/kwOZN/rdzN2RIbkwd0VJESERGXpRJVxaxWK1arFZvNZnSUGuVuMdOgvgcN6ntc8nntmvry3BfbeHP1Xs4W23h2UCfMZhUpERFxPTqdV01q4+m8qvJh8gGe/DQDux2Gx4bw95sisahIiYiIE9ASB+LUbusWxrSbo7CYTXyccpDxH22ixHb5c6pEREScgUqUGGJoTAjW22Nwt5j4YvMRxr6fSmFJ3ToFKiIirk0lSgxzXUQz3ryrK55uZlZsz2bMvI2cLVaREhER16ASJYbq17Epc0bFUd/Dwg+7j5M0J5n8It30WEREnJ9KlBiuR9sA3runG76ebiTvO8kds9eTc6bE6FgiIiKXpBIlTiG2RSM+GNOdBvXdSc88zW1v/cSJ/CKjY4mIiFyUSpQ4jc4h/iy4L4EAH0+2Hcnlljd/Iju30OhYIiIiF6QSVcWsVivh4eHExcUZHcUldQjyZeH93Wnm78Weo/nc/MY6Dp46Y3QsERGR82ixzWqixTYrJ/PkGW6f/ROZJ88S7O/F/DHdaRXgbXQsERGp5bTYpri80Eb1WXR/D9o08eZwTiE3v7GOXdl5RscSEREppxIlTivI34sF9yfQMciXY3lF3PLGOhZuyCTnrK7cExER4+l0XjXR6byqc/pMMUnvJJN+MAcAD4uZvh2aMDi6OVdf0RQvd4vBCUVEpLaoyOe3SlQ1UYmqWvlFpby7dj+fpR1iV3Z++XZvDwv9OwUxMDqYXm0DcLdocFVERBynEuUEVKKqz46sXJamHeaztMMcOn22fHsjbw+u7xzEoKjmdG3RELPZZGBKERFxRSpRTkAlqvrZ7XZSD5xmadohlmUc4Xh+cfljwf5eDIwOZlBUMOHN/DCZVKhEROSPqUQ5AZWomlVqK2PtzydYmn6Yr7dkkfe7+++1aeLN4OjmDI1pTmij+gamFBERZ6cS5QRUooxTWGJj1c6jfJZ2mJU7jlJcWgaAh5uZ126N5rqIZgYnFBERZ6USZSCr1YrVasVms7Fr1y6VKIPlFZbwzdZsPkw+wMZfTmEywV8HdiKpR0ujo4mIiBNSiXICGolyLrYyO1OWbuH9nw4AcH+f1jzRv6Mmn4uIyDm0YrnI/2Mxm3h+cAR/7t8BgDe+38tjC9MoKrUZnExERFyVSpTUGSaTiXH92jLt5ijczCY+SzvM6DkbyC3UCugiIlJxKlFS5wzrEsKc0XF4e1hY+/MJbp61jqycQqNjiYiIi1GJkjqpd7smLHwggSa+nuzIymPYzDW6wbGIiFSISpTUWZ2C/Vk8tgdtmnhzOKeQ4a+v5ae9J4yOJSIiLkIlSuq00Eb1+WRsD7q2aEhuYSkj307mi82HjY4lIiIuQCVK6rwG9T14/954rusURLGtjIc+2MTsH/YaHUtERJycSpQI4OVuwXpHF0b9dxHOF5Zt5/kvtlFWpmXURETkwlSiRP7LYjYxZWA4kwd0BODtH/fx8EebtJaUiIhckEqUyO+YTCbu79OGf90ajbvFxLLNRxj5djI5Z7SWlIiInEslSuQCBkc3593R3fD1dGP9vpOMeGMte45qCQQREfkflSiRi+jRNoCFDyQQ6OfJrux8Eqet5q6317NiWzY2zZUSEanzdAPiaqIbENceh06f5a9Lt7Jieza//dcS0rAed3Vvwc1dQ2no7WFsQBERqTIV+fxWiaomKlG1T+bJM7y//hcWbMjk9H/nSHm6mRkS3Zy7EloQ0dzf4IQiIlJZKlEGslqtWK1WbDYbu3btUomqhQpLbCxNO8zctfvZdiS3fHvXFg0Z2aMl13UKwsNNZ8pFRFyRSpQT0EhU7We320k9cIp31/7C8owjlP53nlRTX09ujw/j9m5hNPXzMjiliIhUhEqUE1CJqluO5hbyQfIB5q8/wLG8IgDczCYGdG5GUkILYls0xGQyGZxSRET+iEqUE1CJqpuKS8v4emsW89btZ8P+U+XbI0P8eWloZ82bEhFxcipRTkAlSrYezuG9db+wJO0QhSVluJlNPHJ1O8b2bYO7RXOmRESckUqUE1CJkt8czy/imc+2sDwjC/h1VOqfI6JoF+hrcDIREfn/KvL5rT+HRapZgI8n1tu78K9bo/Gv587mgznc8O8fmf3DXi3aKSLiwlSiRGqAyWRicHRzvnnsSvp2aEJxaRkvLNvObW/+xIETZ4yOJyIiDlCJEqlBgX5ezBkVx9+Gdcbbw0Ly/pNc96/VzF//CzqzLiLiWlSiRGqYyWTi1m5hfPXolcS3asSZYht/+XQLSXM2kJVTaHQ8ERG5TCpRIgYJbVSfD8d05+kbw/F0M7N61zGuffV7Pt10UKNSIiIuQCVKxEBms4l7erVi2SO9iQptQG5hKY8tSGfs+6kczy8yOp6IiFyCSpSIE2jb1IdPHkhg4rXtcbeY+GprFv1fXc1XW7KMjiYiIhehEiXiJNwsZh66qh1LxvWkY5AvJwqKeeD9FB5bkEbOmRKj44mIyP+jEiXiZDoF+/PZQz15sG8bzCb4dNMhrnn1e1ZsyzY6moiI/I5KlIgT8nSz8Ph1Hfl4bA9aB3hzNK+Ie+dtZPxHmzhVUGx0PBERQSVKxKl1CWvI8vG9uf/K1phN8FnaYa559Xu+zDhidDQRkTpPJUrEyXm5W5h8/RUsfrAn7Zr6cDy/mLHzU3lwfoqu4BMRMZBKlIiLiA5twBeP9OLhq9piMZtYnpHFNdO+57O0Q1pXSkTEACpRIi7E083Cn67twGfjenJFMz9OnSlh/EdpjJmXQnauVjsXEalJdb5EDR06lIYNGzJ8+PDzHnv11Vfp1KkT4eHhPPLII/prX5xGRHN/lj7UkwnX/Lqu1Irt2Vwz7XsWbczUv6ciIjWkzpeo8ePHM2/evPO2Hzt2jBkzZpCSkkJGRgYpKSn89NNPBiQUuTB3i5lHrm7H5w/3IjLEn9zCUv788WZGzdnA4dNnjY4nIlLr1fkS1bdvX3x9fS/4WGlpKYWFhZSUlFBSUkLTpk1rOJ3IH+sY5MfisT144rqOeLiZ+X7XMa59dTUfrD+gUSkRkWrk1CVq9erVDBw4kODgYEwmE0uWLDnvOVarlZYtW+Ll5UV8fDzJyclVcuwmTZowceJEwsLCCA4OJjExkTZt2lTJvkWqmpvFzNi+bVj+SG+6hDUgv6iUJz/N4I7Z68k8ecboeCIitZJTl6iCggKioqKwWq0XfHzBggVMmDCBKVOmkJqaSlRUFP379+fo0aPlz4mOjiYiIuK8r8OHD1/y2KdOneKLL75g//79HDp0iLVr17J69eoq/flEqlrbpj4seqAHT98Yjpe7mbU/n+DaV1dj/W4PRaU2o+OJiNQqbkYHuJQBAwYwYMCAiz4+bdo0xowZw+jRowGYNWsWy5Yt45133mHSpEkApKWlOXTsFStW0LZtWxo1agTADTfcwE8//cSVV155wecXFRVRVPS/NXtyc3MdOq5IZVnMJu7p1YqrOzbliU82s37fSV75eicLN2byzI3hXH1FoNERRURqBaceibqU4uJiUlJSSExMLN9mNptJTExk3bp1ld5/aGgoa9eupbCwEJvNxqpVq+jQocNFnz916lT8/f3Lv0JDQyudQaQyWgZ489F93Xn1liia+nryy4kz3PPuRkbPSWbvsXyj44mIuDyXLVHHjx/HZrMRGHjuX9WBgYFkZWVd9n4SExMZMWIEy5cvJyQkpLyAde/eneuvv56YmBgiIyNp06YNgwYNuuh+Jk+eTE5OTvlXZmamYz+YSBUymUwMjQnh24l9ub9Pa9wtJr7beYz+01cz9cvt5BeVGh1RRMRlOfXpvJqwYsWKiz724osv8uKLL17Wfjw9PfH09KyqWCJVysfTjckDruCWrqE898U2Vu08xhvf7+XT1ENMvr4jQ6KbYzKZjI4pIuJSXHYkKiAgAIvFQnZ29jnbs7OzCQoKMijVr1cLhoeHExcXZ1gGkYtp3cSHOaPieDupKy0a1+doXhGPLUhn+Kx1bDmUY3Q8ERGX4rIlysPDg9jYWFauXFm+raysjJUrV5KQkGBYrnHjxrFt2zY2bNhgWAaRSzGZTFx9RSDfPHYlf+7fgXruFlJ+OcXAGT8yeXEGJ3RTYxGRy+LUJSo/P5+0tLTyK+z27dtHWloaBw4cAGDChAm89dZbvPvuu2zfvp2xY8dSUFBQfrWeiFycp5uFcf3a8u3EPgyKCsZuhw+TD9DvH6uYu2YfpbYyoyOKiDg1k92JlzRetWoV/fr1O297UlISc+fOBWDGjBm88sorZGVlER0dzWuvvUZ8fHwNJz1fbm4u/v7+5OTk4OfnZ3QckT+UvO8kU5ZuZfuRX5fn6Bjky5SBnUho09jgZCIiNacin99OXaJcmUqUuCJbmZ0Pkg/wz292cvpMCQDXdw7iies60qKxt8HpRESqX0U+v536dJ4r0sRycWUWs4m7urfguz/15c7uYZhNsDwji8Rp3/Pc59s4VVBsdEQREaehkahqopEoqQ22H8ll6pc7WL3rGAB+Xm48dFVbRia0xMvdYnA6EZGqp9N5TkAlSmqT1buO8dLy7ezIygMgpGE9Hr+uIwMjm2l9KRGpVVSinIBKlNQ2tjI7n6Qe5J/f7CQ799dlEKJC/Hny+iuIb63J5yJSO6hEOQGVKKmtzhSXMvuHfcz6/mfOFNsAuDY8kEkDOtK6iY/B6UREKkclykBWqxWr1YrNZmPXrl0qUVJrHc0rZPqK3XyUfIAyO7iZTdwRH8YjV7ejsY9ugSQirqnaS9TSpUsrHOqaa66hXr16FX6dq9JIlNQVu7Pz+NuXO1i54ygAvp5ujO3Xhrt7ttLkcxFxOdVeoszmiq2MYDKZ2L17N61bt67ooVyWSpTUNWv3HOfF5dvZevjXxTqD/b3483UddHNjEXEpNbJOVFZWFmVlZZf1Vb9+fUcPIyIuokfbAD5/qBev3hJFsL8Xh3MKeWxBOo9/vFm3kBGRWsmhEpWUlFShU3N33nmnRmNE6gCz2cTQmBC+ndiXP/fvgNkEi1IO8uD8VApLbEbHExGpUppYXk10Ok8Evt6axcMfbqK4tIyE1o15c2Qsvl7uRscSEbmoaj2dd+rUKU6ePAnAsWPHWLx4MVu3bnUsaS2k276I/E//TkHMHR2Hj6cb6/ae4Pa31nMiv8joWCIiVaJCJWr27NnExsbStWtXXn/9dYYOHcrKlSu59dZbmT17dnVldCnjxo1j27ZtbNiwwegoIk6hR5sAPhzTnUbeHmQcymHErHUcOn3W6FgiIpVWodN5kZGRrF+/nrNnzxIWFsa+ffto0qQJOTk59OnTh7S0tGqM6lp0Ok/kXD8fy+eu2es5nFNIM38v3runG22b+hodS0TkHNV2Os/NzY169erRqFEj2rZtS5MmTQDw9/fXJcwickltmvjw8dgetGnizZGcQkbMWkd65mmjY4mIOKxCJcpisVBYWAjA999/X749Pz+/alOJSK0U3KAeix7oQVSIP6fOlHDbWz+xZs9xo2OJiDikQiVqxYoVeHr+ejsHf3//8u1nzpzhzTffrNpkIlIrNfL2YP6Y7vRs25gzxTZGz9nAV1uOGB1LRKTCKlSiLnbarmnTproa7b90dZ7IH/PxdOOdUXFc1ymIYlsZD85P5aPkA0bHEhGpkEqtE5WVlUVQUFBV5qk1NLFc5I/Zyuz85dMMPtqQCcCkAR15oE8bg1OJSF1WI7d9Abj22msr83IRqeMsZhNTh3UuL05/+3IHU5dvR2sAi4grqFSJ0hudiFSWyWRi0oCOTB7QEYA3Vu/liU90vz0RcX6VKlFa1kBEqsr9fdrw8k2RmE2wcONBxn2g++2JiHOrVIkSEalKN8eFMvOOWDwsZr7ems2IWev4bsdRjXqLiFNSiRIRp3JdRBBz7/71fnsZh3IYPXcDN7z2I8s2H8FWpjIlIs6jUiXKYrFUVQ4RkXI92gSw8k99GNO7FfU9LGw7ksu4D1K55tXvWbQxkxLNlxIRJ1CpJQ7k4rTEgUjVOFVQzJy1+5m7Zh+5haUANG9Qj/uubM0tcaF4ueuPORGpOhX5/K50iZo6dSqBgYHcfffd52x/5513OHbsGE888URldu9yrFYrVqsVm83Grl27VKJEqkh+USnzf/qFt37Yx/H8IgACfDy4p1dr7uwehq+Xu8EJRaQ2qNES1bJlSz744AN69Ohxzvb169dz6623sm/fvsrs3mVpJEqkehSW2Fi0MZNZ3+/l0OmzAPh5uZHUoyWje7aikbeHwQlFxJXVaIny8vJi+/bttGrV6pzte/fuJTw8vPyGxXWNSpRI9SqxlbE07TAzV+3h52MFANRzt3B7fBhjercmyN/L4IQi4opqbMVygNDQUNasWXPe9jVr1hAcHFzZ3YuIXJC7xcxNsSH857E+vH5HFyKa+3G2xMbbP+7jype/Y/LizRzNrZt/xIlIzXCr7A7GjBnDo48+SklJCVdddRUAK1eu5PHHH+dPf/pTpQOKiFyK2WxiQOdmXBcRxOrdx7F+u4fk/Sf5MDmT5RlZPDe4E4OigrU4sIhUuUqfzrPb7UyaNInXXnuN4uJi4NdTfE888QTPPPNMlYR0RTqdJ2Kc5H0nefbzrWw9nAvAgIggnh8SQYCPp8HJRMTZ1eicqN/k5+ezfft26tWrR7t27fD0rNtvVipRIsYqsZUx87uf+fe3uykts9PI24MXhkRwfedmRkcTESdmSImC/92QWMPmKlEizmLLoRwmLkpnR1YeAAOjgnluUCca6io+EbmAGp1YDvD2228TERGBl5cXXl5eREREMHv27KrYtYhIpUQ09+ezh3ryUL+2WMwmPk8/zDWvruY/27KNjiYiLq7SJeqZZ55h/PjxDBw4kEWLFrFo0SIGDhzIY489VqfnRImI8/B0szCxfwcWj+1B26Y+HM8vYsy8jUxYmEbOmRKj44mIi6r06bwmTZrw2muvcdttt52z/cMPP+Thhx/m+PHjlQroqnQ6T8Q5FZbYePU/u3jzh73Y7RDo58nfboqkX4emRkcTESdQo6fzSkpK6Nq163nbY2NjKS0trezuXY7VaiU8PJy4uDijo4jIBXi5W5h8/RV8/EACrQK8yc4tYvScDTzx8WbyCjUqJSKXr9IjUQ8//DDu7u5MmzbtnO0TJ07k7NmzWK3WSgV0VRqJEnF+Z4ttvPz1Duas2Q/8emPjl4dH0rNtgLHBRMQwNXp13sMPP8y8efMIDQ2le/fuwK/3zTtw4AAjR47E3f1/NwX9/0WrNlOJEnEdP+09wZ8/Tifz5K/34ruzexiTB1yBt2el1yMWERdToyWqX79+l/U8k8nEt99+W5lDuRSVKBHXUlBUyt++3MF7P/0CQFij+rx2WwzRoQ2MDSYiNcqwdaLkf1SiRFzTj7uP88Qnmzl0+ixuZhMT+3fgvt6tMZu1/p1IXVDj60SJiNQWvdoFsHx8b27o3IzSMjt/+3IHI99J1s2MReQ8Do9E3X333Zf1vHfeeceR3bs8jUSJuDa73c7CjZlMWbqVwpIyGnt78I+bo7QUgkgtVyOn88xmMy1atCAmJoZL7eLTTz91ZPcuTyVKpHbYczSPhz7YVH7bmHt6teLx6zrg6WYxOJmIVIcaKVHjxo3jww8/pEWLFowePZo777yTRo0aORS4NlKJEqk9Ckts/O3LHcxdux+AiOZ+vHZrDK2b+BgbTESqXI3MibJarRw5coTHH3+czz//nNDQUG6++Wa+/vrrS45MiYi4Gi93C38d1InZI7vSsL47Ww7lcuO/f+TjlIN6vxOpw6rs6rxffvmFuXPnMm/ePEpLS9m6dSs+PnX3rzSNRInUTlk5hTy2II11e08AMDg6mBeGRODr5f4HrxQRV2DI1XlmsxmTyYTdbsdms1XVbkVEnEqQvxfv3xvPn/t3wGI28VnaYa5/7Qc2HThldDQRqWGVKlFFRUV8+OGHXHPNNbRv356MjAxmzJjBgQMH6vQolIjUbhaziXH92rLw/gSaN6hH5smzjJi1jtdX/UxZmU7vidQVDp/Oe/DBB/noo48IDQ3l7rvv5o477iAgQPeb+o1O54nUDTlnS3jy0wyWbT4CQK+2AUy7OYqmfl4GJxMRR9TYEgdhYWHExMRgMl18Jd/Fixc7snuXpxIlUnf8tqbUX5du42yJjUbeHvz9pkiuCQ80OpqIVFBFPr8dvrvmyJEjL1me6iqr1YrVatW8MJE6xGQycUtcGLEtGvHwh5vYfiSXMfM2Mjw2hGcGhuOnSecitZLunVdNNBIlUjcVlth49T+7ePOHvdjt0Mzfi5eHR9K7XROjo4nIZdC980REDOLlbmHy9Vew6P4EWjauz5GcQu56O5mnlmRQUFRqdDwRqUJVVqLWr19fVbsSEXF5XVs2Yvn43iQltADg/Z8OMOBfP7D+v+tLiYjrq7LTeWFhYRw4cKAqdlUr6HSeiPxmzZ7jPP7xZg6dPovJBPf0bMXE/h3wctf990ScTbVdnXfzzTdfcLvdbufLL78kPz+/YklrMZUoEfm9vMISXvhiOws2ZgLQpok3/7w5mujQBsYGE5FzVFuJatSoEe+99955C2na7XZuueUWsrOzHUtcC6lEiciFfLsjm0mfZHA0rwizCR7s25ZHrm6Hh5umqIo4g2pb4qBv3774+vpy5ZVXnvdYZGRkxVKKiNRBV3UM5JvHGjJl6VY+SzvMjO/2sGJ7NtNujiY8WH9wibiSal3ioKCggF27dhEaGlrnVjPXSJSI/JEvM47wlyVbOFlQjLvFxPir2/FAnza4WTQqJWKUGlmxHCArK4ugoKALPvbSSy/x/fff061bN3bu3ImPjw9Wq5V69eo5ejiXohIlIpfjeH4RTy7O4Jttv06HiArx5583R9G2qa/ByUTqphpbJ+raa6+94PY5c+Zw8uRJvv76a55//nkWLlzIgAEDePTRRytzOBGRWifAx5M37orl1Vui8PVyI/1gDje89iOfpR0yOpqI/IFKlaiLDWItWLCASZMmAXD33XeTnZ3NiBEj2LBhQ2UOJyJSK5lMJobGhPCfx/rQu10ARaVljP8ojanLt2Mr000lRJxVpUrUxe6d5+7uTnFxMQA9e/akfv36AJSVlVXmcCIitVqQvxdzR3djbN82ALyxei93z91AzpkSg5OJyIVUy+zFe++9l8cffxybzcY999yDr68vr776KgMGDKiOw4mI1BoWs4knruvIv2+LwcvdzPe7jjHY+iO7s/OMjiYi/0+Flji4XIMHD6agoICrr76aFi1acOTIEXr06MELL7xQHYcTEal1BkYF07qJN/fNS2H/iTMMsa7h1VuiubbThS/mEZGaV6mr82JiYti0adMln3Py5EkaNmx40VN/tZWuzhORqnAiv4gH56eyft9JAB5LbM/DV7XFbK5b76kiNaXGrs77owIFv65yXtcKlIhIVWns48n798aX38j41RW7eHB+KvlFpQYnExGt6CYi4uTcLWaeHRzByzdF4mEx89XWLIbNXMMvJwqMjiZSpzlcol544QWWL1+u++WJiNSQm+NC+fC+7jT19WRXdj6DZqzhx93HjY4lUmc5PCfKbDaXn6YLCgqiS5cuxMbGlv9v8+bNqzSoq9GcKBGpLtm5hdz/Xgppmacxm+DJ66/gnl6tNHVCpArUyJyouLg4mjdvzlNPPcWkSZMICAhg8eLFDB8+nLCwMIKCgrj++usd3X2NGTp0KA0bNmT48OHnPfaPf/yDTp06ERERwfvvv29AOhGR8wX6efHRfd0ZHhtCmR1eWLadCQvTKSyxGR1NpE5xuEStX7+e5557jrfeeosVK1bw1FNPsXnzZvLy8li7di1TpkwhJCSkKrNWi/HjxzNv3rzztmdkZPDBBx+QkpLChg0bmDFjBqdPn675gCIiF+DlbuGV4ZFMGRiOxWzi002HuPmNdRzJOWt0NJE6o1ITy0eNGsWuXbvo0KEDXbp0YfLkydhsNuLj4xk7dixvvvlmVeWsNn379sXX9/wbfW7fvp2EhAS8vLyoV68eUVFRfPXVVwYkFBG5MJPJxOierXjv7m40rO/O5oM5DPz3j6z9WfOkRGpCpa/O8/Hx4eWXX2bjxo1s2bKFtm3bXnBkxxGrV69m4MCBBAcHYzKZWLJkyXnPsVqttGzZEi8vL+Lj40lOTq6SY0dERLBq1SpOnz7NqVOnWLVqFYcO6YagIuJ8erQNYOlDvegY5Mvx/GJuf2s9Ty3J0DIIItWsSpY4KC0tpaioiNtuu42QkBBGjx7NyZMnK73fgoICoqKisFqtF3x8wYIFTJgwgSlTppCamkpUVBT9+/fn6NGj5c+Jjo4mIiLivK/Dhw9f8tjh4eE88sgjXHXVVQwbNozu3btjsVgq/TOJiFSH0Eb1WfxgD26PDwPg/Z8O0P/V1azedczgZCK1l8NX5/3tb38jIyODjIwMduzYgZeXF5GRkURHRxMTE8OoUaOqtHSYTCY+/fRThgwZUr4tPj6euLg4ZsyYAfx6g+PQ0FAefvhhJk2adNn7XrVqFTNmzODjjz++6HPuvfdehg4dyg033HDBx4uKiigqKir/Pjc3l9DQUF2dJyI1bu2e4zyxeDOZJ3+dH3Vz1xD+ckM4/vXcDU4m4vxq5Oq8J598knXr1nHTTTexZcsWcnNz+fHHH5kxYwb33HNPtY/aFBcXk5KSQmJiYvk2s9lMYmIi69atq5Jj/DaitXPnTpKTk+nfv/9Fnzt16lT8/f3Lv0JDQ6skg4hIRfVoG8DXj17J6J4tMZlg4caDXDPte/6zTev6iVQlh0tU7969OXHiBM8++yyxsbH07NmThx9+mDlz5pCeno7NVr2X2h4/fhybzUZgYOA52wMDA8nKyrrs/SQmJjJixAiWL19OSEjIOQVs8ODBhIeHc+eddzJnzhzc3C5+v+bJkyeTk5NT/pWZmVnxH0pEpIrU93BjysBOLLo/gdYB3hzNK2LMvI088uEmThYUGx1PpFa4eCv4A99//z0Au3fvJiUlhdTUVFJTU5k/fz6nT5/G09OTzp07V9lE7+qyYsWKiz5WkREtT09PPD09qyKSiEiV6dqyEcvH92b6it28ufpnlqYfZs2e4zw7uBM3dG6mBTpFKsHhEvWbdu3a0a5dO2699dbybfv27WPjxo2XdYNiRwUEBGCxWM677Ux2djZBQUHVdlwREVfj5W5h0oCODIgI4vGPN7MzO4+HPtjE550O8/yQCJr6ehkdUcQlVfrqvBdeeOG8ba1atWLEiBG89NJLld39RXl4eBAbG8vKlSvLt5WVlbFy5UoSEhKq7bh/xGq1Eh4eTlxcnGEZREQuJCq0AZ8/3IvxV7fDzWzi663ZXDNtNZ+kHMTBa4xE6rQKXZ33+OOPn/O93W5n9uzZjBkzBoCXX365SsPl5+ezZ88eAGJiYpg2bRr9+vWjUaNGhIWFsWDBApKSknjjjTfo1q0b06dPZ+HChezYseO8uVI1TffOExFntu1wLo9/ks6WQ7kA9O3QhJeGdia4QT2Dk4kYqyKf3xU6nbdw4UISEhIYMGBA+V8tbm5udOrUyfG0l7Bx40b69etX/v2ECRMASEpKYu7cudxyyy0cO3aMZ555hqysLKKjo/nqq68ML1AiIs4uPNiPJQ/25M0f9jL9P7tZtfMY1766mr/ccAW3xoVqrpTIZajQSNTZs2d58cUX2blzJ1OnTqVt27a0bt2avXv3VmdGl6SRKBFxFXuO5vHnjzez6cBpAIZ1ac5LQzvj5a4FhqXuqcjnt0OLbe7Zs4eJEyfSoUMHPvroI3755ReHw9Y2VqsVq9WKzWZj165dKlEi4hJsZXbe/nEvf/9qJ7YyOzFhDXjjzlia+mnSudQt1V6ifrN06VJ++umnap1A7qo0EiUirujH3ccZ90EqOWdLCPLz4s2RsUSGNDA6lkiNqbESJRenEiUirmr/8QLunbeRPUfz8XQz88qIKAZFBRsdS6RG1MhtX0REpHZqGeDN4gd7cFXHphSVlvHIh5t45esdlJXpb26R36t0idqwYQNXX301kZGRDBs2jOeee46lS5dy4MCBqsgnIiIG8PNy562RXbm/T2sArN/9zH3vpZBfVGpwMhHnUenTeR07diQsLIxBgwaxb98+0tLSSEtL49SpUzRs2JATJ05UVVaXoInlIlLbfLrpIE98kkFxaRntA32YPTKOsMb1jY4lUi1qdE6Ut7c3mzdvpk2bNuds/+WXX0hLS2Pw4MGV2b3L0pwoEalN0jJPc9+8jRzNK6JBfXdm3tGFHm0CjI4lUuVqdE5Uz549OXjw4HnbW7RoUWcLlIhIbRMd2oClD/UiMsSf02dKGPl2Mu/9pOVtpG5zaCRq2LBhREZGEhUVhd1uZ+bMmSxatIiGDRtWR0aXpJEoEamNCktsPPHJZj5LOwzAHfFh/HVQJ9wtuk5Jaodqu+3Lb9q0acOaNWuYOXMmx48fB6B9+/YMHjyY7t27ExMTQ+fOnfHw8HBk9yIi4qS83C1MvyWajkF+vPz1DuavP8Ceo/m8fmcsjbz1ni91S6XnRB06dKh8MvlvX3v37sXNzY0OHTqwefPmqsrqUjQSJSK13crt2Yz/KI38olJCGtZjdlJXOgbp/U5cm+GLbebn55OWlkZ6ejrjxo2r6t07NV2dJyJ1ya7sPO59dyMHTp6hvoeFf4yI4vrOzYyOJeIww0uUaCRKROqOUwXFjPsglbU//7qkzYjYEJ4ZGI6vl7vByUQqrtqvztu8eTNlZWWX/fytW7dSWqoF2kREaqOG3h68e3c3xvZtg8kEi1IOMuBfP7Bh/0mjo4lUK4dKVExMTIUW0UxISNAK5iIitZi7xcwT13VkwX0JhDSsx8FTZ7n5jXX8/asdFJde/h/dIq7Eoavz7HY7Tz/9NPXrX96KtcXFxY4cRkREXEy3Vo34cnxvnv18Gx+nHOT1VT+zetcxpt8STbtAX6PjiVQph+ZE9e3bF5PJVKHXfPDBBzRrVncmG2pOlIjUdV9tOcLkxRmcOlOCp5uZSQM6kpTQErO5Yp8fIjVJE8udgEqUiAgczS3kzx9v5vtdxwDo3S6AV4ZHEeTvZXAykQur0du+yLmsVivh4eHExcUZHUVExHBN/byYOzqO5wd3wsvdzA+7j9N/+mo+Tz9sdDSRStNIVDXRSJSIyLl+PpbPYwvS2HwwB4Ah0cE8OzgC/3paCkGch0aiRETE6bRp4sMnY3vwyFVtMZtgSdphBkxfzdqfjxsdTcQhDpeoM2fOVGUOERGpA9wtZiZc24GPx/agReP6HM4p5I7Z63lx2TaKSm1GxxOpEIdLlL+/f/nNh0VERCqiS1hDlj/Sm9u6hWG3w1s/7GPwjDXsP15gdDSRy+ZwibLZbOesWt67d2+ys7OrJJSIiNR+3p5uTB3WmdkjuxLg48GOrDyGvb6WlF9OGR1N5LJU2ZyotLQ0Cgr0F4SIiFRMYnggy8f3JjLEn5MFxdz+1k98mXHE6Fgif0gTy0VExHBNfb346L7uJF7RlKLSMh78IJW3Vu9FF5CLM6tUifrggw9ITU2lpKSkqvK4PK0TJSLimPoebrxxV1eSElpgt8OLy7czZelWSm269544J4fXierTpw9paWnk5eXh7u5OaWkpt99+O7179yYmJobIyEg8PT2rOq/L0DpRIiKOsdvtvP3jPl5cvh27Ha7u2JR/3x5DfQ+HbvcqUiE1etuX3bt3k5KSQmpqavnX6dOncXNz44orriA9Pb0yu3dZKlEiIpXz1ZYjjP8ojaLSMjo39+ftpK409dPtYqR6GX7vvH379rFx40Y2bdrESy+9VNW7dwkqUSIilZd64BT3vruRkwXFNG9Qjzmj42gf6Gt0LKnFaqREvfDCC3Tp0oXY2FgCAwMdClqbqUSJiFSNX04UMHrOBvYeL8DXy4037oylR9sAo2NJLVUjJcpsNmMymQAICgoqL1S//W/z5s0d2W2toRIlIlJ1ThUUc997G9mw/xTuFhN/GxbJTbEhRseSWqhGSlR8fDxHjhxh9OjRBAQEkJqaSkpKCjt27MBms9GkSRO6dOnC8uXLHfohXJ1KlIhI1SossfHnjzfzefphAB5LbM8jV7ct/4NepCpU5PPb4Usd1q9fz9y5c3nyySeJi4tj2rRptGnThqKiItLS0khNTWXTpk2O7l5EROQcXu4W/nVLNKEN6zFz1c+8umIXmafO8NLQzni4adlDqXmVnlien5/Pc889xxtvvMGDDz7I008/Tf369asqn8vSSJSISPX5YP0Bnv5sC7YyOz3bNmbmHbH413M3OpbUAhX5/K50dffx8eHll19m48aNbNmyhbZt2zJv3rzK7lZEROSibo8P4+2krnh7WFiz5wQjZq3l4KkzRseSOqZKxj9LS0spKiritttuIyQkhNGjR3Py5Mmq2LWIiMgF9e3QlIUPJBDo58mu7HyGWNew7ucTRseSOsThEvW3v/2NO+64g8jISOrXr0+PHj2YOXMm3bp1480338Tf378qc7oM3fZFRKTmdAr2Z8m4noQ38+N4fjF3vr2eN1f/rHvuSY2o1BIHLVu2JCkpidtuu4327dtXdTaXpjlRIiI152yxjb8syWBx6iEABkQE8cqIKHw8dasYqZgaWeLg9/fO8/b2JjIyki5dupR/RUREYLFYHPoBagOVKBGRmmW325m//gDPfr6VEpudNk28eeOuWNo21Qrncvlq/N55v93i5ff3zvP09KRz584kJydXZvcuSyVKRMQYmw6cYuz7qWTlFuLtYeHl4VHcENnM6FjiInTvPCegEiUiYpzj+UU8/MEm1u39daL5mN6teOK6jrhZtJ6UXFqNlChfX19iYmLKb/XSpUsXwsPDtXLsf6lEiYgYq9RWxivf7OSN7/cCEN+qETNu70ITX0+Dk4kzq5ESNXPmTFJSUti4cSPbtm2jrKyMevXqERkZec499CIjIx36IVydSpSIiHP4MuMIExelU1BsI9DPk5l3dCG2RSOjY4mTqvHTeWfPnsXb25snnniCkydPkpqaypYtWyguLsZms1V29y5JJUpExHnsOZrPA++nsOdoPu4WE0/dEM7IhBY6eyLnMWROlNlsJi0trXzkqbS0lG3btmkkSiVKRMQpFBSV8vjHm1mWcQSAoTHNeWloZ+p51N0ryeV8NXrbl4txc3OrswVKREScj7enGzNuj+GpG67AYjbx6aZDDJ25hv3HC4yOJi5KlymIiEidYTKZuLd3a+bfG0+Ajwc7svIYOONHVmzLNjqauCCHS9S9997L66+/zoYNGygqKgLQuWUREXEJ3Vs35ouHe9MlrAF5haXcO28jM1ft0e1ipEKqZMVyNzc3SktLGTZsGH379qVLly5ER0dTv379qs7rMjQnSkTE+RWXlvHCsm3MW/cLAHd2D+PZQRFYzBoUqKtqfMXylJSU8tXKf1ux3GKx0L59e7Zu3VqZ3bsslSgREdfxzo/7eH7ZNux2SLwikH/fFqMJ53WUVix3AipRIiKu5cuMI4xfkEZxaRnRoQ14O6krjX20MGddY3iJqsusVitWqxWbzcauXbtUokREXMiG/Se5992N5JwtoWXj+rx7dzdaNPY2OpbUIJUoJ6CRKBER17TnaD5J7yRz6PRZGnt78PaoOKJDGxgdS2qIU6wTJSIi4oraNvXh03E9iGjux4mCYm59c52WQJALUokSERH5f5r6evHRfQn0ad+EwpIy7ntvI/PX/2J0LHEyKlEiIiIX4OPpxuykrtzcNYQyO/zl0y288vUOrSUl5VSiRERELsLdYubvN0XyaGI7AKzf/cyfFqVTXFpmcDJxBipRIiIil2AymXg0sT0v3xSJxWxiceoh7p67gbzCEqOjicFUokRERC7DzXGhzE7qSn0PCz/uOc6IWevIyik0OpYYSCVKRETkMvXr0JQF9yUQ4OPJjqw8hs1cw67sPKNjiUFUokRERCqgc4g/nz7Yg9ZNvDmcU8hNr69l3c8njI4lBlCJEhERqaDQRvX55IEedG3RkLzCUu56ez1z1uzTlXt1jEqUiIiIAxp6e/D+vfEMjAqmtMzOs59v4+EPN1FQVGp0NKkhKlEiIiIO8nK38Nqt0UwZGI6b2cQXm48w2LqGPUc1T6ouUIkSERGpBJPJxOierVhwf3cC/TzZczSfQTPW8MXmw0ZHk2qmEiUiIlIFYls0YtkjvUlo3ZgzxTYe+mATz36+VQtz1mIqUSIiIlUkwMeT9+7pxoN92wAwZ81+bnvrJ60nVUupRImIiFQhN4uZx6/ryFsju+Lr5UbKL6e48d8/sHbPcaOjSRVTiRIREakG14QH8sXDvbiimR/H84u58+31zFy1h7IyLYNQW6hEiYiIVJMWjb359MEeDI8NocwOL3+1k/veSyHnrO67VxuoRImIiFQjL3cLrwyPZOqwznhYzKzYns2gGT+y9XCO0dGkklSiREREqpnJZOK2bmF8PDaB5g3q8cuJMwybuZaPUw4aHU0qoU6XqMzMTPr27Ut4eDiRkZEsWrTonMe/+OILOnToQLt27Zg9e7ZBKUVEpLaIDGnAskd60bdDE4pKy5i4KJ3JizMoLLEZHU0cYLLX4Rv9HDlyhOzsbKKjo8nKyiI2NpZdu3bh7e1NaWkp4eHhfPfdd/j7+xMbG8vatWtp3LjxZe07NzcXf39/cnJy8PPzq+afREREXElZmZ0Z3+3h1RW7sNuhS1gD3rirK018PY2OVudV5PO7To9ENWvWjOjoaACCgoIICAjg5MmTACQnJ9OpUyeaN2+Oj48PAwYM4JtvvjEwrYiI1BZms4lHrm7Hu6O74eflRuqB0wyxrmFHVq7R0aQCnLpErV69moEDBxIcHIzJZGLJkiXnPcdqtdKyZUu8vLyIj48nOTnZoWOlpKRgs9kIDQ0F4PDhwzRv3rz88ebNm3Po0CGH9i0iInIhV7ZvwpJxPWkV4M2h02e5aeZaVm7PNjqWXCanLlEFBQVERUVhtVov+PiCBQuYMGECU6ZMITU1laioKPr378/Ro0fLnxMdHU1ERMR5X4cP/++eRidPnmTkyJG8+eabDmctKioiNzf3nC8REZE/0rqJD58+2IMebRpTUGzj3nkbmf3DXurwbBuX4WZ0gEsZMGAAAwYMuOjj06ZNY8yYMYwePRqAWbNmsWzZMt555x0mTZoEQFpa2iWPUVRUxJAhQ5g0aRI9evQo3x4cHHzOyNOhQ4fo1q3bRfczdepUnn322cv5sURERM7RoL4H797djWc+28qHyQd4Ydl29hzN57nBEXi4OfV4R53msv/PFBcXk5KSQmJiYvk2s9lMYmIi69atu6x92O12Ro0axVVXXcVdd911zmPdunVjy5YtHDp0iPz8fL788kv69+9/0X1NnjyZnJyc8q/MzEzHfjAREamT3C1mXhoawdM3hmM2wUcbMhn5znpOFRQbHU0uwmVL1PHjx7HZbAQGBp6zPTAwkKysrMvax5o1a1iwYAFLliwhOjqa6OhoMjIyAHBzc+Of//wn/fr1Izo6mj/96U+XvDLP09MTPz+/c75EREQqwmQycU+vVrydFIePpxs/7T3J0Jlr2HM03+hocgFOfTqvuvXq1YuysrKLPj5o0CAGDRpUg4lERESgX8emfDK2B/e8u4H9J84wdOYaXr8jll7tAoyOJr/jsiNRAQEBWCwWsrPPvYohOzuboKAgg1L9erVgeHg4cXFxhmUQERHX1yHIlyXjetK1RUPyCktJmpPMez/9YnQs+R2XLVEeHh7ExsaycuXK8m1lZWWsXLmShIQEw3KNGzeObdu2sWHDBsMyiIhI7RDg48n8MfEM69IcW5mdp5ds4a9Lt1Jqu/hZFKk5Tn06Lz8/nz179pR/v2/fPtLS0mjUqBFhYWFMmDCBpKQkunbtSrdu3Zg+fToFBQXlV+uJiIi4Ok83C/8cEUXbpj68/NVO5q7dz97jBcy4PQY/L3ej49VpTn3bl1WrVtGvX7/zticlJTF37lwAZsyYwSuvvEJWVhbR0dG89tprxMfH13DS8+m2LyIiUtW+2pLFYwvSOFtio21TH95O6kqLxt5Gx6pVKvL57dQlyhVZrVasVis2m41du3apRImISJXaciiHe9/dSFZuIQ3ruzPrzljiW1/efV3lj6lEOQGNRImISHXJzi1kzLyNbD6Yg4fFzKy7unBVx8A/fqH8Id2AWEREpBYL9PNiwX0J9O8USLGtjPvfS+GbrZe3RqJUHZUoERERF1TPw8KM27twY2QzSmx2HpyfypcZR4yOVaeoRImIiLgod4uZ6bdEMyQ6mNIyOw99uInP0w8bHavOUImqYlpsU0REapKbxcw/b45meGwItjI74z/axJJNh4yOVSdoYnk10cRyERGpSWVldp78NIOPNmRiMsErw6MYHhtidCyXo4nlIiIidYzZbOKloZ25Iz4Mux3+/HE6CzYcMDpWraYSJSIiUkuYzSZeGBJBUkIL7HZ44pMM3tf99qqNSpSIiEgtYjKZ+OugTtzTqxUATy3Zwrtr9xsbqpZSiRIREallTCYTT91wBff3aQ3AlKVbmf3DXoNT1T4qUVVMV+eJiIgzMJlMTLquIw/1awvAC8u2M+v7nw1OVbvo6rxqoqvzRETEGdjtdv61cjfTV+wGYOK17XnoqnYGp3JeujpPREREgF9HpB5NbM/Ea9sD8I9vdjF9xS40hlJ5KlEiIiJ1wENXtWPSgI4ATF+xm39+oyJVWSpRIiIidcQDfdrw1A1XADDjuz387asdKlKVoBIlIiJSh9zbuzV/HRgOwBvf71WRqgSVKBERkTpmVM9WPD8kAvi1SM1cpav2HKESVcW0xIGIiLiCu7q34C/X/3pq75Wvd/KeVjavMC1xUE20xIGIiLiCf36zk39/uweTCabfEs3g6OZGRzKUljgQERGRyzLhmvbl99r708J0vt2RbXQkl6ESJSIiUoeZTCamDOzE0JjmlJbZGft+Kj/tPWF0LJegEiUiIlLHmc0mXh4eSeIVTSkqLePedzeScTDH6FhOTyVKREREcLeYmXF7F7q3bkR+USkj31nPnqN5RsdyaipRIiIiAoCXu4XZSXFEhfhz6kwJd85OJvPkGaNjOS2VKBERESnn4+nG3NHdaNfUh6zcQu56ez1H8wqNjuWUVKKqmNaJEhERV9fQ24P37oknpGE99p84w8i3k8k5U2J0LKejdaKqidaJEhERV/fLiQKGz1rHsbwiuoQ14P1746nv4WZ0rGqldaJERESk0lo09ua9e7rhX8+d1AOnuf+9FIpKbUbHchoqUSIiInJRHYP8mDM6jvoeFn7YfZxHP0qj1FZmdCynoBIlIiIil9QlrCFv3tUVD4uZL7dk8eSnGWg2kEqUiIiIXIZe7QJ47bYYzCZYuPEgLy7bXueLlEqUiIiIXJbrIoJ4eXgUALN/3Me/v91jcCJjqUSJiIjIZRseG8IzN4YDMO0/u7B+V3eLlEqUiIiIVMjdvVox8dr2ALzy9U5e/c+uOnlqTyVKREREKuyhq9oxaUBHAP61cjcvf72zzhUplSgRERFxyAN92vD0f0/tvb7qZ16oY5PNVaKqmG77IiIidck9vVrx/JAIAN7+cR9Tlm6lrKxuFCnd9qWa6LYvIiJSlyzYcIBJizOw2+HWuFBeGtoZs9lkdKwK021fREREpEbdEhfGP0dEYTbBRxsy+fPHm7HV8hEplSgRERGpEsO6hPCvW2OwmE18knqQxxbU7lvEqESJiIhIlRkYFYz19hjcLSaWph/m4Q83UVxaO4uUSpSIiIhUqesimjHrztjye+09OD+FolKb0bGqnEqUiIiIVLmrrwjkraSueLqZWbH9KPfNS6GwpHYVKZUoERERqRZ92jdhzqg46rlb+H7XMe55dwNnikuNjlVlVKJERESk2vRoG8C7d3fD28PCmj0nGDVnA/lFtaNIqUSJiIhIterWqhHz7onH19ON5H0nGfn2enILS4yOVWkqUSIiIlLtYls0ZP6YePzruZN64DR3zV5PzhnXLlIqUSIiIlIjIkMa8MGYeBp5e5B+MIeRc5Jdeo6USpSIiIjUmE7B/nw4pjsN67uTnnmase+nUuKiC3KqRImIiEiN6hDkyzu/u2rviY83u+RNi1WiREREpMbFhDVk5p1dsJhNLN50iL9/tcPoSBWmElXFrFYr4eHhxMXFGR1FRETEqfXr0JS/3xQJwBur9zL7h70GJ6oYk91ud73xMxeQm5uLv78/OTk5+Pn5GR1HRETEab2+6ufykah/3RrN4OjmhmWpyOe3RqJERETEUA/0ac3oni0BmLgonR92HzM20GVSiRIRERFDmUwmnr4hnIFRwZTY7DzwXgoZB3OMjvWHVKJERETEcGaziX+MiKRn28YUFNsYNSeZ/ccLjI51SSpRIiIi4hQ83SzMujOWTsF+nCgoZuQ7yRzNKzQ61kWpRImIiIjT8PVyZ+7oboQ1qs+Bk2cYPWcDeU56nz2VKBEREXEqTXw9mXd3NwJ8PNh6OJcH3k+hqNRmdKzzqESJiIiI02kZ4M2cUd3w9rCwZs8J/rQw3elWNVeJEhEREafUOcSfWXfF4m4x8cXmIzz3xTacaXlLlSgRERFxWr3bNeEfI6IAmLt2P69//7PBif5HJUpERESc2uDo5jx9YzgAL3+1k0UbMw1O9CuVKBEREXF69/Rqxf19WgMwaXEG3+7INjiRSpSIiIi4iEnXdWRYl+bYyuw8OD+V1AOnDM2jEiUiIiIuwWQy8febIunboQl2O5wqKDY0j5uhRxcRERGpAHeLmZl3dGHP0XwiQxoYmkUjUSIiIuJS6nu4GV6gQCVKRERExCEqUSIiIiIOUIkSERERcYBKlIiIiIgD6nyJyszMpG/fvoSHhxMZGcmiRYvOeXzo0KE0bNiQ4cOHG5RQREREnFGdL1Fubm5Mnz6dbdu28c033/Doo49SUFBQ/vj48eOZN2+egQlFRETEGdX5EtWsWTOio6MBCAoKIiAggJMnT5Y/3rdvX3x9fQ1KJyIiIs7K6UvU6tWrGThwIMHBwZhMJpYsWXLec6xWKy1btsTLy4v4+HiSk5MdOlZKSgo2m43Q0NBKphYREZHazulXLC8oKCAqKoq7776bYcOGnff4ggULmDBhArNmzSI+Pp7p06fTv39/du7cSdOmTQGIjo6mtLT0vNd+8803BAcHA3Dy5ElGjhzJW2+95VDOoqIiioqKyr/Pzc11aD8iIiLiGpy+RA0YMIABAwZc9PFp06YxZswYRo8eDcCsWbNYtmwZ77zzDpMmTQIgLS3tkscoKipiyJAhTJo0iR49ejiUc+rUqTz77LMOvVZERERcj9OfzruU4uJiUlJSSExMLN9mNptJTExk3bp1l7UPu93OqFGjuOqqq7jrrrsczjJ58mRycnLKvzIzMx3el4iIiDg/px+JupTjx49js9kIDAw8Z3tgYCA7duy4rH2sWbOGBQsWEBkZWT7f6r333qNz584AJCYmkp6eTkFBASEhISxatIiEhITz9uPp6Ymnp2flfiARERFxGS5doqpCr169KCsru+jjK1asqME0IiIi4ipcukQFBARgsVjIzs4+Z3t2djZBQUGGZLJarVit1vKJ7JpgLiIi4jp++9y22+1/+FyXLlEeHh7ExsaycuVKhgwZAkBZWRkrV67koYceMiTTuHHjGDduHAcPHiQ0NFTLJYiIiLigvLw8/P39L/kcpy9R+fn57Nmzp/z7ffv2kZaWRqNGjQgLC2PChAkkJSXRtWtXunXrxvTp0ykoKCi/Ws8owcHBZGZm4uvri8lkqvDr4+Li2LBhQzUkqxnOkr8mc1TXsapqv5Xdj6Ovr8jrcnNzCQ0NJTMzEz8/vwofS87lLP8dOsqZ8tdUFr2PVM3rKvNeYrfbycvLK18C6VKcvkRt3LiRfv36lX8/YcIEAJKSkpg7dy633HILx44d45lnniErK4vo6Gi++uqr8yab1zSz2UxISIjDr7dYLC79IeIs+WsyR3Udq6r2W9n9OPp6R17n5+fnFP/+uDpn+e/QUc6Uv6ay6H2kal/n6HvJH41A/cbpS1Tfvn3/8LzkQw89ZNjpu+oybtw4oyNUirPkr8kc1XWsqtpvZffj6Oud5d+FusjVf/fOlL+msuh9pHqOW11M9suZOSUiUgNyc3Px9/cnJyfHaUYgRMT11NR7iUsvtikitYunpydTpkzRmmsiUik19V6ikSgRERERB2gkSkRERMQBKlEiIiIiDlCJEhEREXGASpSIiIiIA1SiRERERBygEiUiTm3o0KE0bNiQ4cOHn/fYmTNnaNGiBRMnTjQgmYi4iou9j7Rs2ZLIyEiio6PPuTvK5VKJEhGnNn78eObNm3fBx1588UW6d+9ew4lExNVc6n1k7dq1pKWl8d1331V4vypRIuLU+vbti6+v73nbd+/ezY4dOxgwYIABqUTElVzsfaSyVKJEpNqsXr2agQMHEhwcjMlkYsmSJec9x2q10rJlS7y8vIiPjyc5Ofmy9j1x4kSmTp1axYlFxNlU5/uIyWSiT58+xMXFMX/+/ApnU4kSkWpTUFBAVFQUVqv1go8vWLCACRMmMGXKFFJTU4mKiqJ///4cPXr0kvv97LPPaN++Pe3bt6+O2CLiRKrrfQTgxx9/JCUlhaVLl/LSSy+xefPmCmVzq9CzRUQqYMCAAZc83TZt2jTGjBnD6NGjAZg1axbLli3jnXfeYdKkSRd93U8//cRHH33EokWLyM/Pp6SkBD8/P5555pkq/xlExFjV9T4C0Lx5cwCaNWvG9ddfT2pqKpGRkZedTSNRImKI4uJiUlJSSExMLN9mNptJTExk3bp1l3zt1KlTyczMZP/+/fzjH/9gzJgxKlAidVBl3kcKCgrIy8sDID8/n2+//ZZOnTpV6PgaiRIRQxw/fhybzUZgYOA52wMDA9mxY0f594mJiaSnp1NQUEBISAiLFi0iISGhpuOKiBOqzPtIYGAgQ4cOBcBmszFmzBji4uIqdHyVKBFxaitWrLjk46NGjaqZICLisi72PpKenl6p/ep0nogYIiAgAIvFQnZ29jnbs7OzCQoKMiiViLgSo99HVKJExBAeHh7ExsaycuXK8m1lZWWsXLlSp+tE5LIY/T6i03kiUm3y8/PZs2dP+ff79u0jLS2NRo0aERYWxoQJE0hKSqJr165069aN6dOnU1BQUH6VjYiIM7+PmOx2u73ajyIiddKqVasueD+qpKQk5s6dC8CMGTN45ZVXyMrKIjo6mtdee434+PgaTioizsqZ30dUokREREQcoDlRIiIiIg5QiRIRERFxgEqUiIiIiANUokREREQcoBIlIiIi4gCVKBEREREHqESJiIiIOEAlSkRERMQBKlEiIiIiDlCJEhEREXGASpSISA0aOnQoDRs2ZPjw4UZHEZFKUokSEalB48ePZ968eUbHEJEqoBIlIvI7ffv2xWQyYTKZSEtLO+exiRMnMmTIkErv39fXt1L7cMSoUaPKf64lS5bU+PFFaiOVKBExVJ8+fTCZTLz00kvnbLfb7cTHx2MymXjuuedqNNOYMWM4cuQIERER52xPS0sjMjKy2o5bnb+Lf/3rXxw5cqQqYorIf7kZHUBE6i673c6mTZto0aIFGRkZ5zz27rvvcvjwYQC6dOlSo7nq169PUFDQedvT09MZO3bsRV8XHR1NaWnpedu/+eYbgoODL3nM6v5d+Pv74+/v79BrReTCNBIlIobZvXs3eXl5JCUlnVMc8vLymDx5MqNGjQIgNjbWoIT/c/DgQY4fPw7ANddcQ/369enQoQPr168vf05aWhpbtmw57+uPChS41u9CRH6lEiUihklJSaF+/frcdttt7Ny5k+LiYgCef/55unbtSpMmTQgKCqJZs2YGJ6V8fpTVauXJJ58kPT2dsLAwJk2aVCX7d6XfhYj8SqfzRMQwqampREZG0qFDB7y8vNixYwf16tXj9ddfJzU1lRdffLHGT+VdTFpaGo0aNWLhwoUEBAQAMGjQIN54440K7ScxMZH09HQKCgoICQlh0aJFJCQkuNTvQkR+pZEoETFMamoqXbp0wWQyERkZSUZGBo899hhjx46lXbt2pKSkEBsby/79++natWuF9u3Iay4lLS2NwYMHlxcogH379tG2bdsK7WfFihUcO3aMM2fOcPDgQRISEoDL/12IiPNQiRIRw/xWHODXSdnTp09n48aNPP300xQWFrJjxw6nGX1JS0uje/fu522Ljo6ukv1X5nexZMkSbrzxRq699lrefvvtKskjIn9MJUpEDLF3715Onz5dXgxiYmLYuHEjU6dOxdfXl/T0dEpLS8tHX0pKSkhKSuKKK67glltuwW63l+/rxhtvJDY2loiICObPn3/esf7+978TERFB586dL/j4H8nLy2Pv3r3ExMScs72qSlRFfxe/N3/+fBYuXMjrr7/O+++/z7Zt23jxxRcrnUlELoNdRMQACxcutHt4eNiLi4vtdrvdXlxcbD927Ji9rKzMbrfb7TNnzrQ3adLEbrfb7fv27bO7u7vbt27dai8rK7P36dPHvnr16vJ9nThxwm632+35+fn2K664wl5YWGjft2+fPTY21p6cnGzv0qWL/ezZs/YTJ07YW7dubT906NBFc/Xp08c+fvz4c7b98MMPdjc3N/vZs2fLt+3fv98O2Pft21ejv4v/78orryx/3W8SExPtubm5F3w+YP/0008rnVlE7HaNRImIIVJTU4mIiMDd3R0Ad3d3AgICMJlM5Y//fuSnQ4cOhIeHYzKZiImJYf/+/eWPvfrqq0RFRdGjRw8OHDjAgQMHyh9bs2YNN910E15eXjRq1Iirr76aDRs2VChrWlpa+YTv32zatIkGDRrQsmVLB376c1X0d/GbEydOEBYWhru7O7Nnzy5fBiE+Pp5du3ZVOpeIXJpKlIgYYurUqaSkpFz08bfeeouvv/66/HtPT8/yf7ZYLNhsNgC+++471qxZw/r160lPT6djx44UFRVVadaHHnqILVu2nLNtyJAhnDp1qkr2X9HfxW8aNmxYvgr5iBEjmDp1KgA7d+4kNDS0SrKJyMWpRImIS8vNzaVx48Z4eXmRlpZGenr6OY/36tWLxYsXU1RUxKlTp/j222/p1q3bJfc5c+ZMfHx8zls53NmYzWZ69OjBtGnT8Pf3p1mzZixcuBBvb2+aNm16znMfeOABfHx8DEoqUjtpnSgRcWnXXXcdr7/+OuHh4XTq1Om8ydddu3ZlxIgRxMbGYjKZePbZZy+5YOX8+fM5e/YsAGFhYdWavSo888wzPPXUU0RFRWEymYiLi8NqtZ73vOeee46JEycCaMFOkSpistt/d4mLiIiIiFwWnc4TERERcYBKlIiIiIgDVKJEREREHKASJSIiIuIAlSgRERERB6hEiYiIiDhAJUpERETEASpRIiIiIg5QiRIRERFxgEqUiIiIiANUokREREQc8H/wIQ3vSzWiPwAAAABJRU5ErkJggg==\n"},"metadata":{}}]}],"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":5} \ No newline at end of file