From 035830cc0cf57c9e3dd08f8c2fe2b13db01f0a2b Mon Sep 17 00:00:00 2001 From: Hongwan Liu Date: Wed, 25 Sep 2019 16:15:10 -0400 Subject: [PATCH] Successful compilation + HyREC FULL mode * Successful compilation of the python wrapper based on settings in Makefile and setup.py. Do not mess with these anymore! * Turned on HyREC FULL mode so that we're no longer solving using the RECFAST mode, which was dumb. --- Makefile | 6 +- hyrec/history.h | 2 +- hyrec/hydrogen.h | 8 +- .../recfast_fudge-checkpoint.ipynb | 274 +++++++++++++++++- notebooks/recfast_fudge.ipynb | 147 ++++++---- python/cclassy.pxd | 5 +- python/setup.py | 7 +- 7 files changed, 376 insertions(+), 73 deletions(-) diff --git a/Makefile b/Makefile index c246071b3..55f656a90 100755 --- a/Makefile +++ b/Makefile @@ -16,9 +16,13 @@ vpath .base build ###### LINES TO ADAPT TO YOUR PLATEFORM ################ ######################################################## +# Successfully compiled with gcc 9.2.0, in the python3 environment, +# with all the openmp lgomp flags turned on. Note that the GCCPATH +# seems to have to be specified manually in python/setup.py. + # your C compiler: # CC = gcc -CC = /usr/local/Cellar/gcc/9.1.0/bin/gcc-9 +CC = /usr/local/Cellar/gcc/9.2.0/bin/gcc-9 #CC = icc #CC = pgcc diff --git a/hyrec/history.h b/hyrec/history.h index 4b7fff0d4..0ab98d290 100644 --- a/hyrec/history.h +++ b/hyrec/history.h @@ -17,7 +17,7 @@ #define FULL 3 /* All radiative transfer effects included. Additional switches in header file hydrogen.h */ /** here is the switch **/ -#define MODEL RECFAST /* default setting: FULL */ +#define MODEL FULL /* default setting: FULL */ /***** Switches for derivative d(xe)/dt *****/ diff --git a/hyrec/hydrogen.h b/hyrec/hydrogen.h index 815037b12..948c36e75 100644 --- a/hyrec/hydrogen.h +++ b/hyrec/hydrogen.h @@ -39,12 +39,8 @@ double rec_HRecFast_dxedlna(double xe, double nH, double H, double TM, double TR /************* EFFECTIVE MULTI LEVEL ATOM *******************/ -// #define ALPHA_FILE "Alpha_inf.dat" /* Contains the effective recombination coefficients to 2s and 2p */ -// #define RR_FILE "R_inf.dat" /* Contains the effective transfer rate R_{2p,2s} */ - -// HL: Changed to variables in the hope that they can be accessed. -char * ALPHA_FILE = "Alpha_inf.dat" -char * RR_FILE = "R_inf.dat" +#define ALPHA_FILE "Alpha_inf.dat" /* Contains the effective recombination coefficients to 2s and 2p */ +#define RR_FILE "R_inf.dat" /* Contains the effective transfer rate R_{2p,2s} */ /* Boundaries and number of elements of temperature tables */ diff --git a/notebooks/.ipynb_checkpoints/recfast_fudge-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/recfast_fudge-checkpoint.ipynb index 2fd64429b..a20aaf99f 100644 --- a/notebooks/.ipynb_checkpoints/recfast_fudge-checkpoint.ipynb +++ b/notebooks/.ipynb_checkpoints/recfast_fudge-checkpoint.ipynb @@ -1,6 +1,276 @@ { - "cells": [], - "metadata": {}, + "cells": [ + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# import classy module\n", + "%autoreload 2\n", + "from classy import Class\n", + "\n", + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "did it set correctly???????\n", + "b'./hyrec/fake_Alpha.dat'\n" + ] + }, + { + "ename": "CosmoComputationError", + "evalue": "\n\nError in Class: thermodynamics_init(L:388) :error in thermodynamics_recombination(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination(L:2601) :error in thermodynamics_recombination_with_hyrec(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination_with_hyrec(L:2736) :could not open fA with name ./hyrec/fake_Alpha.dat and mode \"r\"", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mCosmoComputationError\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 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 11\u001b[0m \u001b[0;31m# run class\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mLambdaCDM_new\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute\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[0m", + "\u001b[0;32mclassy.pyx\u001b[0m in \u001b[0;36mclassy.Class.compute\u001b[0;34m()\u001b[0m\n", + "\u001b[0;31mCosmoComputationError\u001b[0m: \n\nError in Class: thermodynamics_init(L:388) :error in thermodynamics_recombination(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination(L:2601) :error in thermodynamics_recombination_with_hyrec(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination_with_hyrec(L:2736) :could not open fA with name ./hyrec/fake_Alpha.dat and mode \"r\"" + ] + } + ], + "source": [ + "from classy import Class\n", + "\n", + "LambdaCDM_new = Class()\n", + "# pass input parameters\n", + "LambdaCDM_new.set({\n", + " 'recombination': 'hyrec', \n", + " 'Alpha_inf hyrec file': './hyrec/fake_Alpha.dat'\n", + "})\n", + "LambdaCDM_new.set({'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n", + "\n", + "# run class\n", + "LambdaCDM_new.compute()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "did it set correctly???????\n", + "b'../hyrec/Alpha_inf_orig.dat'\n", + "did it set correctly???????\n", + "b'../hyrec/Alpha_inf.dat'\n" + ] + } + ], + "source": [ + "# create instance of the class \"Class\"\n", + "from classy import Class\n", + "\n", + "LambdaCDM = Class()\n", + "hyrec_mod = Class()\n", + "# pass input parameters\n", + "\n", + "hyrec_mod.set({\n", + " 'recombination': 'hyrec',\n", + " 'Alpha_inf hyrec file': '../hyrec/Alpha_inf.dat',\n", + " 'output':'tCl,pCl,lCl,mPk',\n", + " 'lensing':'yes','P_k_max_1/Mpc':3.0\n", + "})\n", + "LambdaCDM.set({\n", + " 'recombination': 'hyrec',\n", + " 'Alpha_inf hyrec file': './hyrec/Alpha_BB_n_250_fine.dat',\n", + " 'output':'tCl,pCl,lCl,mPk',\n", + " 'lensing':'yes','P_k_max_1/Mpc':3.0\n", + "})\n", + "\n", + "# run class\n", + "LambdaCDM.compute()\n", + "hyrec_mod.compute()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/hongwan/anaconda/envs/python3/lib/python3.6/site-packages/ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in true_divide\n", + " del sys.path[0]\n" + ] + }, + { + "data": { + "text/plain": [ + "Text(0, 0.5, '$(C_\\\\ell - C_\\\\ell^\\\\mathrm{ref})/C_\\\\ell^\\\\mathrm{ref} [\\\\%]$')" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEMCAYAAAAS+xsDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFNZJREFUeJzt3X+MZWWd5/H3x25AmTH8bKQF2m6kJy4yCWKl0YizZkagMdHGVTPN7GDvLG7HDST+mslCDIFRJ+rssuwadTY9oxGJiviDoWfHWeSXa8YsSLWiggjdMs7SQgCBYYZFQOC7f9xTWpRV3aeq79O36vb7ldzce57zVN3vU+emP33Oee45qSokSWrpeaMuQJI0/gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5paPuoDF4vDDD6/Vq1ePugxJWjK2bdv2s6pa0aevYdNZvXo1k5OToy5DkpaMJP/Yt6+H0SRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkppbtGGTZH2SO5PsSHL+LOsPSPLFbv3NSVbPWL8qyWNJ/nhv1SxJmt2iDJsky4BPAGcAxwNnJTl+RrdzgEeq6jjgUuCjM9ZfCvxd61olSbu3KMMGWAfsqKq7q+op4Apgw4w+G4DLutdfBn4vSQCSnAncDdy+l+qVJO3CYg2bo4B7pi3v7Npm7VNVTwOPAocl+Q3gPwF/uhfqlCT1sFjDJrO0Vc8+fwpcWlWP7fZNks1JJpNMPvjggwsoU5LUx/JRFzCHncAx05aPBu6do8/OJMuBg4CHgZOBtyb5c+Bg4NkkT1TVx2e+SVVtAbYATExMzAwzSdKQLNawuQVYm2QN8FNgI/AHM/psBTYB/wd4K3BDVRXw2qkOSS4GHpstaCRJe8+iDJuqejrJecA1wDLg01V1e5IPAJNVtRX4FHB5kh0M9mg2jq5iSdKuZLAzoImJiZqcnBx1GZK0ZCTZVlUTffou1gkCkqQxYthIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktTcog2bJOuT3JlkR5LzZ1l/QJIvdutvTrK6az81ybYkP+ief3dv1y5Jeq7lfTolObRHt2er6p/2sJ6p91sGfAI4FdgJ3JJka1X9cFq3c4BHquq4JBuBjwK/D/wMeGNV3ZvkBOAa4Khh1CVJWpheYQPc2z2yiz7LgFV7XNHAOmBHVd0NkOQKYAMwPWw2ABd3r78MfDxJquq70/rcDjw/yQFV9eSQapMkzVPfsLmjql6xqw5Jvrur9fN0FHDPtOWdwMlz9amqp5M8ChzGYM9myluA7xo0kjRafcPm1UPq09dse1A1nz5JXs7g0Nppc75JshnYDLBq1bB2yiRJM/WaIFBVTwyjzzzsBI6Ztnw0g8N4s/ZJshw4CHi4Wz4auAp4e1X9eBc1b6mqiaqaWLFixRDLlyRNt6DZaEleleSGJN9KcuawiwJuAdYmWZNkf2AjsHVGn63Apu71W4EbqqqSHAz8LXBBVX2rQW2SpHnqFTZJjpzR9F7gTcB64IPDLqqqngbOYzCT7A7gyqq6PckHkryp6/Yp4LAkO7p6pqZHnwccB1yY5NbuccSwa5Qk9ZeqmadCZumU/DWwDfjPVfVEki3AJPAs8EdV9Zq2ZbY3MTFRk5OToy5DkpaMJNuqaqJP377nbM4EbgX+Z5KzgXczCJoDgRaH0SRJY6T3OZuq+hvgdOBg4KvAnVX1sap6sFVxkqTx0PeczZuS/D1wA3AbgxP2b07yhSQvbVmgJGnp6/s9mw8x+B7NC4CvVdU64L1J1gJ/xiB8JEmaVd+weZRBoLwAeGCqsaq2Y9BIknaj7zmbNzOYDPA08AftypEkjaO+ezZfr6qTdtUhyXd210eStG/qGzb/Ksn3d7E+DC4XI0nSr+kbNi/r0eeZPSlEkjS+dhs2Sd4C/LSqbtoL9UiSxlCfPZu3A8uSfG2qoao+2a4kSdK46RM27wH+I4O7XkqSNG+7DZvu1sx/shdqkSSNqb6Xq3lXkr/sXl/YtiRJ0rjp+6XOlwL3dK9f2KgWSdKY6hs2BbwgyQnAixvWI0kaQ7sNmyQBHmLwxc2zgQtaFyVJGi+7DZsa3MpzLfA94JvAbwMkOSXJuUmOneqbZE2rQiVJS1ffKwhcB+wPHM7gkBrACmAdsC7JQ8DngfcBZw27SEnS0tYrbKrqslnarkqyFXglcALwWuCu4ZYnSRoHvcKmm+78eFVdMr29qp4Bvt09JEmaVd/DaGcDJ85sTPIOYEVVfXioVUmSxkrfqc8/r6rHZ2m/HPjDIdYjSRpDvcMmycqZjVX1JIO7d0qSNKe+YXMJcHWSl0xvTHIE8OzQq5IkjZW+s9G+lORAYFuSm4BbGQTV24CL25UnSRoHffdspqY/HwtcCewHPAGcVVWfa1SbJGlM9J36/Grgpqr6Z+CzbUuSJI2bvns2mxgcQrsiyb9LcmTLoiRJ46XvOZt3AiR5GXAG8JkkBwE3Av8L+Fb3BU9Jkn5N73M2AFX1o6q6FHgLcDrw9wwmCdw87MKSrE9yZ5IdSc6fZf0BSb7Yrb85yepp6y7o2u9Mcvqwa5MkzU/fczbPAzYC/xaYAJ4CDgAeBL7GkC++mWQZ8AngVGAncEuSrVX1w2ndzgEeqarjkmwEPgr8fpLju1pfzuDeO9cl+S33vCRpdPru2dzI4G6dFwArq+qYqjqCwcU3bwI+nGSYVxJYB+yoqrur6ingCmDDjD4bgKkLhH4Z+L3u3jsbgCuq6smq+gdgR/f7JEkj0vfaaK+vql/MbKyqh4GvAF9Jst8Q6zqKX92GGgZ7NyfP1aeqnk7yKHBY137TjJ89aoi1/dJTTz/Lu674botfLUl7xZEHPZ+L3vjy5u+z27BJ8hbgpzz3H/BfM1sY7YHM9hY9+/T52cEvSDYDmwFWrVo1n/q6X1r8+MHH5v1zkrRY/OKZWf95HLo+ezZvB5Yl+dpUQ1V9sl1JwGBv5Jhpy0cD987RZ2eS5cBBwMM9fxaAqtoCbAGYmJiY91/8gOXL+Pp7/vV8f0yS9jl9ztm8B7gDuH3ao7VbgLVJ1iTZn8EJ/60z+mxl8P0fgLcCN3S3sN4KbOxmq61hcEtr77cjSSO02z2bqrob+JO9UMv093w6yXnANcAy4NNVdXuSDwCTVbUV+BRweZIdDPZoNnY/e3uSK4EfMrgi9bnORJOk0cpgZ2A3nZJ3ASdU1X9IcmFVfbB9aXvXxMRETU5OjroMSVoykmyrqok+fftOfX4pv5od9sIFVSVJ2mf1DZsCXpDkBAZflJQkqbfdhk33RcmHGEwpPpvBFzslSeptt2HTzfBaC3wP+Cbw2wBJTklybpJjp/p2s78kSXqOvlcQuA7YHzicX31BcgWDy8CsS/IQ8HngfQz5OmmSpKWv7y0GLpul7aokW4FXAicwuE7aXcMtT5I0Dvpe9flC4PGqumR6e/f9lW/jlyYlSbvQ9zDa2cCJMxuTvANYUVUfHmpVkqSx0nfq88+r6vFZ2i8HhnlrAUnSGOodNklWzmysqicZXBJGkqQ59Q2bS4Crk7xkemOSI4Bnh16VJGms9J2N9qUkBwLbktwE3MogqN4GXNyuPEnSOOi7ZzM1/XkNcCWwH/AEcFZVfa5RbZKkMdF3NhoAVfUvwGcb1SJJGlO992wkSVqoBYVNkjcOuxBJ0vha6J7Nnw21CknSWFto2GSoVUiSxtpCw2b395KWJKnjBAFJUnOGjSSpuYWGzf1DrUKSNNYWFDZVdeqwC5EkjS8Po0mSmjNsJEnNGTaSpOYWermaU5Kcm+TYaW1rhleWJGmczOuqz9OsANYB65I8BHweeB9w1rAKkySNjwWFTVVdlWQr8ErgBOC1wF3DLEySND56hU2SC4HHq+qSqbaqegb4dveQJGlOffdszgZOnNmY5B3Aiqr68FCrkiSNlb4TBH5eVY/P0n458IdDrIckhya5Nsn27vmQOfpt6vpsT7Kpazswyd8m+VGS25N8ZJi1SZIWpnfYJFk5s7GqngSeHm5JnA9cX1Vrgeu75edIcihwEXAyg4kKF00Lpf9SVS8DXgG8JskZQ65PkjRPfcPmEuDqJC+Z3pjkCODZIde0Abise30ZcOYsfU4Hrq2qh6vqEeBaYH1VPV5VNwJU1VPAd4Cjh1yfJGmeep2zqaovJTkQ2JbkJuBWBkH1NuDiIdf0oqq6r3vf+7pAm+ko4J5pyzu7tl9KcjDwRuC/D7k+SdI89Z2Nlqq6LMlXgTcDLwf+H3BWVU1O69PrpmpJrgOOnGXV+/uVPeudQn/53kmWA18APlZVd++ijs3AZoBVq1b1fGtJ0nz1nY12Y5KvAFdX1WenGpPsn+R3gU3AjcBn+vyyqnr9XOuS3J9kZbdXsxJ4YJZuO4HXTVs+GvjGtOUtwPaq+m+7qWNL15eJiQnvPipJjfQ9Z7MeeAb4QpJ7k/wwyd3AdgZXDbi0qj4zpJq2MggvuuerZ+lzDXBakkO6iQGndW0k+RBwEPDuIdUjSdpDfc/ZPAF8Evhkkv2AwxlMh/6nBjV9BLgyyTnA/2VwXogkE8A7q+odVfVwkg8Ct3Q/84Gu7WgGh+J+BHwnCcDHq+qvGtQpSeopPU+zjL2JiYmanJwcdRmStGQk2VZVE336eosBSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaW3Rhk+TQJNcm2d49HzJHv01dn+1JNs2yfmuS29pXLEnanUUXNsD5wPVVtRa4vlt+jiSHAhcBJwPrgIumh1KSfwM8tnfKlSTtzmIMmw3AZd3ry4AzZ+lzOnBtVT1cVY8A1wLrAZL8JvBe4EN7oVZJUg+LMWxeVFX3AXTPR8zS5yjgnmnLO7s2gA8ClwCPtyxSktTf8lG8aZLrgCNnWfX+vr9ilrZKciJwXFW9J8nqHnVsBjYDrFq1qudbS5LmayRhU1Wvn2tdkvuTrKyq+5KsBB6YpdtO4HXTlo8GvgG8Gnhlkp8wGNsRSb5RVa9jFlW1BdgCMDExUfMfiSSpj8V4GG0rMDW7bBNw9Sx9rgFOS3JINzHgNOCaqvqLqnpxVa0GTgHumitoJEl7z2IMm48ApybZDpzaLZNkIslfAVTVwwzOzdzSPT7QtUmSFqFUefQIBofRJicnR12GJC0ZSbZV1USfvotxz0aSNGYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRJDVn2EiSmjNsJEnNGTaSpOYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRJDVn2EiSmjNsJEnNGTaSpOYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRJDVn2EiSmjNsJEnNGTaSpOZSVaOuYVFI8iDwjwv40cOBnw25nMXOMe8bHPO+YU/G/JKqWtGno2Gzh5JMVtXEqOvYmxzzvsEx7xv21pg9jCZJas6wkSQ1Z9jsuS2jLmAEHPO+wTHvG/bKmD1nI0lqzj0bSVJzhs0eSLI+yZ1JdiQ5f9T1DFOSnyT5QZJbk0x2bYcmuTbJ9u75kK49ST7W/R2+n+Sk0VbfT5JPJ3kgyW3T2uY9xiSbuv7bk2waxVj6mmPMFyf5abetb03yhmnrLujGfGeS06e1L4nPfpJjktyY5I4ktyd5V9c+ttt5F2Me7XauKh8LeADLgB8DxwL7A98Djh91XUMc30+Aw2e0/Tlwfvf6fOCj3es3AH8HBHgVcPOo6+85xt8BTgJuW+gYgUOBu7vnQ7rXh4x6bPMc88XAH8/S9/juc30AsKb7vC9bSp99YCVwUvf6hcBd3bjGdjvvYswj3c7u2SzcOmBHVd1dVU8BVwAbRlxTaxuAy7rXlwFnTmv/bA3cBBycZOUoCpyPqvom8PCM5vmO8XTg2qp6uKoeAa4F1revfmHmGPNcNgBXVNWTVfUPwA4Gn/sl89mvqvuq6jvd638B7gCOYoy38y7GPJe9sp0Nm4U7Crhn2vJOdr1Bl5oCvp5kW5LNXduLquo+GHyggSO69nH6W8x3jOMy9vO6w0afnjqkxJiNOclq4BXAzewj23nGmGGE29mwWbjM0jZOU/teU1UnAWcA5yb5nV30Hfe/Bcw9xnEY+18ALwVOBO4DLunax2bMSX4T+Arw7qr65111naVtXMY80u1s2CzcTuCYactHA/eOqJahq6p7u+cHgKsY7FLfP3V4rHt+oOs+Tn+L+Y5xyY+9qu6vqmeq6lngLxlsaxiTMSfZj8E/up+rqq92zWO9nWcb86i3s2GzcLcAa5OsSbI/sBHYOuKahiLJbyR54dRr4DTgNgbjm5qFswm4unu9FXh7N5PnVcCjU4colqD5jvEa4LQkh3SHJU7r2paMGefX3sxgW8NgzBuTHJBkDbAW+DZL6LOfJMCngDuq6r9OWzW223muMY98O4965sRSfjCYuXIXgxkb7x91PUMc17EMZp58D7h9amzAYcD1wPbu+dCuPcAnur/DD4CJUY+h5zi/wOBwwi8Y/C/unIWMEfj3DE6q7gD+aNTjWsCYL+/G9P3uH5OV0/q/vxvzncAZ09qXxGcfOIXBoZ/vA7d2jzeM83bexZhHup29goAkqTkPo0mSmjNsJEnNGTaSpOYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRFrHuag4f777NLi1Zho20uL2TwX1GThl1IdKeMGykxW09g8uF3DrqQqQ9YdhIi1SS5zO4W+JJwP8ecTnSHjFspMVrLYOw+VFV/WLUxUh7YvmoC5A0pxXAb7FIb7kszYd7NtLi9WIGN8B63rRb+EpLkmEjLUJJljM4V3Mk8D+AZ0ZbkbRnvJ+NJKk592wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKa+//6Mq33aMve2QAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "cls_new = hyrec_mod.lensed_cl(2500)\n", + "cls_std = LambdaCDM.lensed_cl(2500)\n", + "\n", + "ll_vec_new = cls_new['ell']\n", + "ll_vec_std = cls_std['ell']\n", + "\n", + "plt.figure()\n", + "\n", + "ax = plt.gca()\n", + "\n", + "# ax.set_yscale('log')\n", + "\n", + "plt.plot(ll_vec_new, (cls_new['tt']/cls_std['tt'] - 1.)*100)\n", + "# plt.plot(ll_vec_new, (cls_new['ee']/cls_std['ee'] - 1.)*100)\n", + "\n", + "\n", + "# plt.axis([0, 2500, -1e-10, 1e-10])\n", + "\n", + "plt.xlabel(r'$\\ell$')\n", + "plt.ylabel(r'$(C_\\ell - C_\\ell^\\mathrm{ref})/C_\\ell^\\mathrm{ref} [\\%]$')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.max(\n", + " new.get_thermodynamics()['x_e'] \n", + " - LambdaCDM.get_thermodynamics()['x_e']\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%autoreload\n", + "# LambdaCDM_hyrec = Class()\n", + "# LambdaCDM_hyrec.set({'recombination': 'hyrec', 'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n", + "LambdaCDM_newAlpha = Class()\n", + "LambdaCDM_newAlpha.set({\n", + " 'recombination': 'hyrec', \n", + " 'hyrec_Alpha_inf_file': './hyrec/Alpha_inf.dat',\n", + " 'output':'tCl,pCl,lCl,mPk',\n", + " 'lensing':'yes',\n", + " 'P_k_max_1/Mpc':3.0\n", + "})\n", + "# LambdaCDM_hyrec.compute()\n", + "LambdaCDM_newAlpha.compute()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cls_hyrec=LambdaCDM_hyrec.lensed_cl(2500)\n", + "cls_newAlpha = LambdaCDM_newAlpha.lensed_cl(2500)\n", + "\n", + "ll_vec_hyrec = cls_hyrec['ell']\n", + "\n", + "plt.figure()\n", + "\n", + "ax = plt.gca()\n", + "\n", + "# ax.set_yscale('log')\n", + "\n", + "plt.plot(ll_vec_hyrec, (cls_newAlpha['tt']/cls_hyrec['tt'] - 1.)*100)\n", + "plt.plot(ll_vec_hyrec, (cls_newAlpha['ee']/cls_hyrec['ee'] - 1.)*100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.loadtxt('../hyrec/Alpha_inf.dat')\n", + "a = a/1e5\n", + "np.savetxt('../hyrec/Alpha_inf.dat', a)" + ] + }, + { + "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.6.9" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, "nbformat": 4, "nbformat_minor": 2 } diff --git a/notebooks/recfast_fudge.ipynb b/notebooks/recfast_fudge.ipynb index 5616b90f4..b46f458f9 100644 --- a/notebooks/recfast_fudge.ipynb +++ b/notebooks/recfast_fudge.ipynb @@ -2,21 +2,21 @@ "cells": [ { "cell_type": "code", - "execution_count": 15, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ - "load_ext autoreload" + "%load_ext autoreload" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# import classy module\n", - "%autoreload\n", + "%autoreload 2\n", "from classy import Class\n", "\n", "%matplotlib inline\n", @@ -26,27 +26,54 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "did it set correctly???????\n", + "b'../hyrec/Alpha_inf.dat'\n", + "did it set correctly???????\n", + "b'../hyrec/Alpha_BB_n_250_fine.dat'\n" + ] + } + ], "source": [ "# create instance of the class \"Class\"\n", - "LambdaCDM_std = Class()\n", - "LambdaCDM_new = Class()\n", + "from classy import Class\n", + "\n", + "LambdaCDM = Class()\n", + "\n", + "LambdaCDM.set({\n", + " 'recombination': 'hyrec',\n", + " 'Alpha_inf hyrec file': '../hyrec/Alpha_inf.dat',\n", + " 'output':'tCl,pCl,lCl,mPk',\n", + " 'lensing':'yes','P_k_max_1/Mpc':3.0\n", + "})\n", + "\n", + "LambdaCDM.compute()\n", + "\n", + "hyrec_mod = Class()\n", "# pass input parameters\n", - "LambdaCDM_new.set({'recfast_fudge_H':1.18})\n", - "LambdaCDM_std.set({'recfast_fudge_H':1.14})\n", - "LambdaCDM_new.set({'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n", - "LambdaCDM_std.set({'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n", + "\n", + "hyrec_mod.set({\n", + " 'recombination': 'hyrec',\n", + " 'Alpha_inf hyrec file': '../hyrec/Alpha_BB_n_250_fine.dat',\n", + " 'output':'tCl,pCl,lCl,mPk',\n", + " 'lensing':'yes','P_k_max_1/Mpc':3.0\n", + "})\n", + "\n", "\n", "# run class\n", - "LambdaCDM_new.compute()\n", - "LambdaCDM_std.compute()" + "\n", + "hyrec_mod.compute()" ] }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -65,13 +92,13 @@ "Text(0, 0.5, '$(C_\\\\ell - C_\\\\ell^\\\\mathrm{ref})/C_\\\\ell^\\\\mathrm{ref} [\\\\%]$')" ] }, - "execution_count": 19, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEMCAYAAAAS+xsDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNXd+PHPd7JCVgIBwr4rCIgQENzFfaXuS13rUp/W1ta2T+3T9dc+bdVWbX2qrVStS23drVTcUNG6gYDs+w4hIQkh+56Z8/vj3IEhTJI7WzIh3/frlddM7px759wE5ptzzvecI8YYlFJKqVjydHUFlFJKHfk02CillIo5DTZKKaViToONUkqpmNNgo5RSKuY02CillIo5DTZKKaViToONUkqpmNNgo5RSKuYSu7oC8aJfv35mxIgRXV0NpZTqVpYtW7bPGJPbUTkNNo4RI0awdOnSrq6GUkp1KyKy0025uO1GE5FzRWSjiGwRkXuCvH6TiJSKyArn69aA124Ukc3O142dW3OllFKtxWXLRkQSgEeAs4ACYImIzDPGrGtV9AVjzJ2tzs0Bfg7kAwZY5pxb3glVV0opFUS8tmxmAFuMMduMMU3A88Acl+eeAywwxux3AswC4NwY1VMppZQL8RpsBgO7A74vcI61dpmIrBKRl0VkaIjnKqWU6iTxGmwkyLHWG+/8GxhhjJkMvAc8HcK5tqDI7SKyVESWlpaWhl1ZpZRS7YvXYFMADA34fghQGFjAGFNmjGl0vv0rMM3tuQHXmGuMyTfG5Ofmdpi5p5RSKkzxGmyWAGNFZKSIJANXA/MCC4hIXsC3FwPrnefvAGeLSB8R6QOc7RxTSinVReIyG80Y0yIid2KDRALwpDFmrYj8ElhqjJkHfFtELgZagP3ATc65+0XkV9iABfBLY8z+Tr8JpZSKlDFQVQgl66G6COrKwNsEnkRITofsoZAzGvqOAU+8th0sMSbocEaPk5+fb3RSp1I9kDGwdxXs/AwKl9sP9/pySEiGtH7QfwKMPAVGnAyJybGvT3UxbHnPfm37EOpd/K2cmg3DT4CjzoejL4DeOTGvpp+ILDPG5HdYToONpcFGqU5SthU2vgm7Ftm/2BurIKk39B9vPywnXQHJvWNfj/KdsOwpWPsalG+3xzIGQZ/h0CsHvI1QvRdKN4KvGXr1gem3woyvQ/qhY7xen6GstpGSqkaKqxqoaWyhscVHU4uPpAShd3IiaSkJ5KanMig7lZy0ZEQCcpkaKmHdPFj9Imz/GDCQPgBGz4bB02zAyx4KvftCYip4m+3PrWI3lK6H3Yth64dQuQskATPiROpHnUNB/9PZ6e1HYUU9JdUN1DV5qW/yUt/spdnrQxAQ+O2lk8hMTQrrx6jBJkQabNQRxxjY+SmseRWKVkBLo/3QmnoDjDy5c+vSUAkrn7cf7iXO3OycUTBgov0rvLEa9iyD8h2Q1h/O+bUNOhIsuTRCuxbB54/AhjcAsa2WiZfCmDMhc9Dh5ZsbYNtCWP53zIb5mKTebD3667yffTnrSprYuLeabftqaPa6/yxNSfQwPDuR81PXcGbzhxxd/RmJviYaM4djJl5ByqQ5yMBJbd5/TWMLxVUNFFc2UFjZQGFFPYXldSSWrmZ8+UJmNn7OaNkDwHrfUJb7xrCekVQk9qMxIYPeSYZMaaC/r4QBvhLOuusxstPDC/AabEKkwUYdUbb/B97/JRQsgaQ0GDLNth52f2G7ZaZcBxc+CIkpri/Z0OyluqEFj0BSooeMlMRD/zoPprEGFj0Knz4MTdUwaCpMvgqOOs+2IAIZY7uy3vu5rfekK+Di/4OkXgFFDOV1zezeX8fu8jr2VTdS2+SlprEFY+yHeEqSh35pKQzISmVgZioDMlPI6pWEFCyFD34F2z+y3U75N8P02yAr+DS8msYWNhVXs76oyvmqpmnvBu70Pcc5CUvZ7cvliZTrKBh0HqMHZjIkuxf9M1Ppn5FCZq8kkhM8pCR6aPEZ6ppaqG5oobSqgfrdK8jb/ioT9r1Nuq+KMrKY1zKTf3lPZKUZDQgZKYn0SUsmNclDalICxkB9s5eGZi/ltU3UNnkPq2+/9BQGZ6cyKLsXg7J7cVRSMZNrPmVw2eekla3G01gZ/HeU1Bu+tSx4oHVBg02INNioiHmbYd3rsPld24oYexYcew14EjqvDg2V8M7/wPK/Q+YQOOX7MPlKSE6zrzfXw39+Bx8/AGPOgmv+CQmHd5/UN3lZunM/n28tY9G2MjburT7sAy45wUNuRgpD+vTimEFZTBycyazRfcnLcoJD8Vp4/qu2i+roC+Hk78HgqR3fg89r67fw1zQMnMZHUx/my30JrCqoZE1hJdUNLYedkuAREkRo8voOe22C7OAHyS9zunxJpSeL9/tex9Zhl5PcK4PEBCHBI9Q1eamqb6a8rold++vYVVZHWW3TgWtkpCRydF4G4/MyOXpgJvm+1Yxe/hsSStbY1tnJ37NjJcGCtzFQusF2Ha55FYrX2PGgoy+AY6+F0bOpajY2gO6vp6C8jt3766isb6ah2UdDixcBeiUnkJqYQGavJAZm2SA6IDOVQVm9GJiVSmpSO//OjIGqPVBTYv+NJCTZIJM9zHbNRdCC1GATIg023Yy32X6gFiyxA7fHXh2bLhe3itfBa1+3A83pA+2HSeUuGHW6/UAP+Os8ZnYthpdvtllLJ94Fp94DSanByy79G7zxHTjuOrj4TzS0+PhyVzmLtpbx+bYyVuyuoNlrSPQIk4dkcezQbHIzUshIScRnoNnrY19NEyXVDWzfV8v6oioamu0H/dj+6dzRdwVf2fVbJDUTzxVPwoiTOqx+eW0TK3ZXsHx3BSt2V9Bv19v8xjxMicnm694fkpQ3nomDsxidm87QPr0YmtOb/hkppKUkkpLoQUTw+QwNLV7Kapqo3LWWnCW/Z9Cet6lPyGBBn6t5yXM+O2uEvVUNNLUcDEwiNqBk9U5iSHZvhvftzbC+vRmdm86EvEyG9Ol1eCvO54O1r8LCX8P+bc4g/YmQO85mijVW2fGpwhVQVWDPGZwPU66BYy7t1EH8WNJgEyINNt1Ic739i3nr+5CSBY2VcMp/w+wfd019Ni+AF2+wfyle8HsYP8d+en35NPz7OzDpcrjs8bAu7fMZRGi/u8oYWPI4vP0j2yV02ZO226wdjS1e9r3+Mwav/hN/zrqbh8pm0NTiwyMwaXAWM0f3ZdaovkwfkUNaSsczJLw+w6biaj7duJdBS+/l/JpXWOIbxzea7qLPgKGMHZDBkD69yOmdTHKiB6/PUFHXTFltEzv21bK1tIaSajtH2yMwbkAGU4Zmc3rGbmYvv4tEXyNy5dN2wLwjFbvgw/tg5T/s72TmN2DWN6FX9mF1bvH58PoMKYkJJHjC/GPF57VZY2tehV2fQ8VO8LXYPziyhsKgKTYIHXVe2F1V8UyDTYg02HQTxtgWxKoX4cKHYOqNMO9OWPEPuONjGDipc+uz6iVbnwHHwFdfhowBh77+0f32L9+r/g7jL2rzMmU1jSzZUc7yXeVsLqlha2kNxVUNNDT7SE7wkJOWzKjcNCbkZTJxcBaTh2Qxom8aHuOF+d+FL5+BsefApY/ZrKlWmlp8rN5TwedOy2XZznKamlt4Nvm3TPNs4emJTzJ24nTyR+SEnZVEdTG8cgvs+JiW/NtZMu67LCuoZdnOcrbvq2VPRf0hg+gegT69kxnutCDG9E9n8pBsJg/JOjTAVeyCf1xtu6JmfRNO/SGkpB/+/kWrYPFjsOoFEA/MuA1O+q5NX+5Mxthu1MSUrm1tdxINNiHSYNNN+Lt/Tv8xnPrf9lh9OfzhWBh9Olz5dPvnB2ho9rJydwXri6oorGygsdlLSlICw3J6M2lwFpMGZ+Fp76/djW/ZFtawWbarLDXz8DLeZnjsVGiph29+cWB8xOczLN9dzvxVe/loUwlbS2sBSE70MKpfGqP7pzMoK5VeyYk0tfjYV9PI5uJq1u+tPtD90y/Vx2Opf2JawyI2j/s6+4//AX3SU/H6DLWNLezaX8eWkhqW7ixn5e4KGp3zjh6YwSyn5TIzt4XMp0+D3v3g9oXhd/dteR9eu8NmlV30B9ut2YrPZ6hv9jotKCEjNbH9n2+ghip498c2qKZmwYQ5dqxEPFC2BbYuhH0bbUtmyrVw0t1tDvyr6NJgEyINNt1AVSH8aYbtlrhh3qEzpt/+ke1K+t7GdvvCaxtbeHfdXv61vJDPt5YdGFBOTvTQOzmB+ibvgQ/lfukpzJkyiOtmDmdkv7RDL7TzM3j2Ejs35MZ/Q0pG2/Xe+Bb882p8Fz3M8n4X8+bqIt5cXURRZQPJCR5mje7L8aNyOH5kDhMHZ5GS2PZAb7PXx5aSGtZv38nUT7/BsNrV/Mp7I39rPjto+QSPMHFQJvkjcpg+og8zRvYlJ63VxMQt78PfL4VpN9tAEYqyrbDgZzaNuN9RcMVTMGBCaNcIRcFSm9225T070A2Q2AuGHQ9HXWCTIVp1l6nY0mATIg023cDzX7UfMv/1GfQdfehre1fDX06CCx6E6bccduq20hqe+mwHLy8roK7Jy+DsXpw3cSAzR/Xl2KHZ9Eu3k+yMMRRVNvDF9v28s3YvC9YV0+IznDoul5tPHMEpY3PxlG2Gx8+E9P7wtXcgrW+bVfb5DCt3lzPwhXNpqK/l9Pr7SE5I4JRx/bhgch5njB8QerdVfTk8fbGdEHnpXBqOmsP2fbWU1TRRUd9EokdISUpgaJ/eDMvpTXKii2VMFvwcPv2DDRbHXNJx+cZq+M/v7Qe/JwlOvtt2cXVGIgTYwfm6MjA+SMuN+6VajmQabEKkwSbObX4PnrsMzvyF7YdvzRh4eArkHg3XvuAcMvxn8z7+9ul2PtxYSnKCh4uOHcTVM4YybVgfV104JVUN/OOLXTy3eBel1Y0c28/Hs77/IZ16PLd/YFNHW2lo9vLlznLeXVfMO2v3UlTZwBWJH/O7xD/zycy5TD7t0vDHRRoq4Zmv2PTZq56DccFbNCHzNsOT58K+TXDbQug3Jng5Y+x42YKfQc1em9p9xs8hMy94eXXE02ATIg02ccznta2WZmfco631qeZ/H1Y8R913N/PKqjKe+nQ7W0tryc1I4brjh3Pt8cPIzXA/iTFQU4uPt1buYsTbN3B001qubf4JTXnTGds//UC3VGlNIzv21bKuqIpmryEl0cOp43I555iBnDkum6y/HAd5x8J1L4f3c2ishmcvhcIv4cpn4ejzw7tOW8p3wl9n24HtG/99eOuxcDm8+d9Q8AUMOg7O+x0MnR7dOqhux22wictVn5U6xIp/2CVOrniq3YUQ9+WdQr8lf+Wu+//CgobxTB6SxUNXHcsFkwa560pqR3KihznFj0DzSnaf+gCn+E5l8Xab2VVZ3wzYMZ4hfXpxy0mjmD6iD7NG96V3csB/sem3woe/gX2bod/Y0CrQVAvPXWmXdLnib9EPNGBn9N/wOjx9Icw9zWZ9jTkDaortwPyaV2yX1ZxH7GRE7bpSIdCWjSMuWzb7tkBDBQzp8I+GzuHvQqnbZz84Q1jqxG9fTSM7y+rYW9lAfbOX1CQPI/qmMW5ARvCA0FQL/zcNsobALQsOSyU1xvDZ1jKe+XwHn6/bzvLk23i7340MvPjnTB3Wp+PlVNxa+bxNcZ51p123Kxw1pfCQszbZBQ+4P6+pDv5xpV3n7LLHYeJl4b2/W+U7bTr59v8cPJbUG2b+l50smpoV2/dX3Yq2bLo7nw8en2376L+/2Q5Gd7XPHrZ99WD79i/6Y4enVNY18/6GYj7YUMLyXRXsqagPWi4tOYFTj8rlzPEDOP2o/vTxZ0wtetTOiL/iqUMCza6yOt5aU8QLS3ezrbSWPr2T+Oqpk/BtmcAFmTtgeBRnZ5duhDe+ayfmnfn/wr9Oeq5d72vFP2D2T4LOhzlMcwO88FXY8Qlc8ljsAw3YFs6N/7YJCMVrISXTLl8fbG6LUi5psIlXpesPpnbu+LhzPmTaU1kAC39j17jKGGhX7z3tfw6fxOhYX1TFE59sZ96KQpq8PvpnpDBjZA43nTCCMQPSGZiZSu/kBOqavGwtreGzrWW8t66YN1fvxSMwbXgfZg1K4M5Vf6Rs4Gw+LR1C6batbC6pZsXuCrY581KmDM3mgSuO5YLJeXZtqPkn2A9zbwskROGfd1MdvHij/cv+siciv+bxd8CK5+DLZ+HEb7dftqURXrwetn7gdF1dFdl7h6r/ePulVBRosIlXuxcffL5vS9QuW9vYQkl1I9m9kg62HtxY9Ge7BMe5v7UD9Uset4tOHn/7IcV276/jvrc38MaqInolJXDV9KFcOnUwxw7JbjP7a3xeJhdOHsT/zpnImsJK3ltfwsINJaQunUuyp5qbd57Nhh0rATsuMmVoFtfOGMY5xwxkaE6rZdGHzYQlf7XZWoOmhPSzCeqtH9iZ69e/Gp2Mq7zJdi23L+bawNPWGFRLo10CZ/O7cOEf7BpmSnVjGmziVckGSM6w/eP7t0Z0Ka/P8O+Vhfx90U6W7SrHP0w3Pi+Tm08YwaVTB5OY0M5gb3O9/Wt8/EUHU32zh9ml2p1gU1nXzCMfbuGpT3fg8cC3Zo/hlpNGkt3bfUDzeMRZriSbu0/oi/njAuqHX8yjZ19Pgkfol57S8Tpdw2bax12LIg82G9+2i32e/H13a3K5deJd8NzlNmDP+sbhr7c0wUs3waa37dhO/s3Re2+luogGm3i1fxvkjLD95ZUFYV9mVUEFP3hpFRuLqxmdm8a3Th/D8L5plFQ3Mn91If/9yiqeXbSTh66awpj+bfTJr3vdTiTM/9rBYyNPhfXz8Hq9/HNJAQ+8u5GK+mYunzqE7519FAOz2lht2K1P/4A019Hr7J8wKjeEsYKsIXZp/d2LYeYd4b9/QxXMv9tuNnbqD8O/TjBjzrTB68N77TLzgfu61FfAC9fZrtPzf28TMZQ6AsRtsBGRc4E/AgnA48aYe1u9fjdwK9AClAJfM8bsdF7zAqudoruMMRd3WsWjpXy7XdwRObizYQiMMTz64VYeXLCJ3PQUHrl2KudNHHhIV9Ydp47ijVVF/Oz1NVzw8Mf85ILxXDdz+OEZXGtetS2ZEQG7Ow6ZDsuf5bY/vswHJenMGJnDzy+awDGDopCpVF0MX/zVDqbnHhX6+UNn2E3CIvH+L+3yOFc+E/1950VsIJl7OvzjKruuWp8Rdmxm/t1QuQcumdv5YzRKxVBcBhsRSQAeAc4CCoAlIjLPGBP4qbscyDfG1InIfwH3A/7/nfXGmCh02HcRn9emnx59oU39rS0N6fTGFi8/emU1ry7fw4WT8/j1JZPI6nX4jHUR4aJjB3H8yBx+8PIqfvr6WhZuLOW+yyYfnPzYWG23xJ1+24FssC93lfP6kgT+H5BXv5lHrr2F8ycNjF6a8ScPgbcp/BbF0OPtPiOVBbalE6pdi2wX1/F3xC7tvO9ouOoZeP46ePg424JtrLRB56Y3DnYHKnWEiMtgA8wAthhjtgGIyPPAHOBAsDHGLAwovwg4ckZQq4vA12y7V2pKbReWtznojoqtldc28fVnl/HFjv3cfdY4vjV7TIdBoH9mKk/dPJ2nP9vBb97awOwHPuSGWcO5cPIgxu57j0RvE3vyZvPB5zt4dfkelu+qYGDvPvycBH4x3UvS5CguVVJVCEuftBtMtZ7B7tbQGfZx9xehB5uWRpj3bXve7J+E9/5ujToNvrnYjgvVFNvAdswlnbe+mFKdKF6DzWBgd8D3BcDx7ZS/BXgr4PtUEVmK7WK71xjzr+hXMYaq99rHjEEcGM2vK7Mpx+3YWlrDLU8tobCygYevOY6Lj3W/UZOIcNOJIzlpbD8eXLCJRz/cyiMLt/LHpL9ykieDk/9Zj4+1jBuQzk8uGM81M4bheeJoPKVrw73L4D75Axiv3QwtXAMn2ZWAd38BEy8N7dyPH7RL1X/15c6ZV5I1GE6L8piQUnEoXoNNsD/Fgy51ICLXAfnAqQGHhxljCkVkFPCBiKw2xhyW0iUitwO3AwwbdviCil3mQLAZYLuTwO4d3k6w+XxrGXf8fRmJHuGftx3PtDAnNY7pn8GjX51GcVUDn2ws5Ny3V7K57xn8dvqxTBuew+jctIMtpdxxULQyrPcJqnafXRZl8tWHDpqHKiEJBk87NH3cjZL18PEDMOlKGHtW+O+vlDpMvC5uVAAMDfh+CFDYupCInAn8GLjYGNPoP26MKXQetwEfAscFexNjzFxjTL4xJj83Nzd6tY9UdZF9TB94cG8O/wTPIF5cupvrn1hM/4wU/vXNE8MONIEGZKZyWZ/tpHhrmXjGtVw1fRhj+qcf2iXXd4wdW2ppivj9AJsU0FLf8WRHN4bOgL2r7KRMN3xemPctuy/Nub+N/P2VUoeI12CzBBgrIiNFJBm4GpgXWEBEjgMewwaakoDjfUQkxXneDziRgLGebqGm2O5AmJZrB44haLDx+Qz3vrWB/355FbNG9+WVb5xw+CTHSGx4A5LS7NhCMH3H2i6v8h2Rv1dTrZ3oeNT54WWgtTb0eDsJtfBLd+WXPAEFS+C8+zp/G2GleoC4DDbGmBbgTuAdYD3wojFmrYj8UkT8acy/A9KBl0RkhYj4g9F4YKmIrAQWYsdsulewqd5rA01C4sGthhurDilS3+TlG899yV8+2sq1xw/jyZumh79HSjA+H2x4067629aAdV9nz5OyKKxwsPw5qN9vJzxGwxBn6Xs3XWkVu+H9/2fnv0y6Ijrvr5Q6RLyO2WCMeRN4s9WxnwU8P7ON8z4DJsW2djFWUwzpzppjKc68lYaDwaa4qoHbnlnK6j2V/OSC8dxy0sjopR37FX5pN8c6+sK2y/izxco2R/Ze3hb4/P9g6Mzopfym9bUtr10dBBtj7NwWY+DChw5bVVopFR1xG2x6tOoim4kGh7Vslu7Yz3899yU1DS3MvT6fsyYEXwgzYuv/DZ7E9neC7JVtW2D7Igw26/4FFbvg3Psiu05rI06yWyI0N0BSGysarHzerj927r1Bd91USkVHXHaj9XjVxQdXU05IgsRemIZKnvl8B1fPXURacgL/+uaJsQs0ABvm2w/rjpbB7zvGLq0TLmPg0z9Av3Ew7tzwrxPMUedDc61dnj+Y6r3w9g9ti2rG16P73kqpQ2iwiTfeFrtiQPrBNGeTksni9Tv42etrOWVcLq/feRJHDcyIXR1KN9qusfa60PxyRkNZBAuFblsIe1fDCd+O/s6PI0+xWwNseOPw14yBN+62kzjnPKK7TioVY/o/LN7UlgLmwJyagvI6CuqT2FdWyl1njOXxG/KDLj0TVf4P56NcbD3cd5Qd22msDu+9Pv0jZOTB5CvDO789Sakw/mJY/fLh9Vv8F9g4H874GfQbE/33VkodQoNNvPHPsckYyCeb93HR/31CuTeVWYOT+O5Z49rcEyaqNsyHQcfZ2e0d8WekhdOVVrgctn1otxsOY4tpV6bfCk3Vdg6P39aF8O5P4KgLYGaQJf6VUlGnCQLxpqYYgOfXN/E/XyxmTP90xmQOprevtnPev6oI9iyD2T91Vz7Hn5G2FfKODe29Pn3YziOadlNo54Vi6HTbHfjR/TZ7rr4c3roHco+Grzyq2WdKdRJt2cSZyhK7JNwfF1czZ8pgXvvGifTOyA6/mypUG51sczfjNQA5o+xjqOM2+7fbLLT8m+0GcbF0wYOQPdTufPnvu+ymatf/6+DqDEqpmNOWTRxZuKGEjR8s4Q7g+5edzGXTR9oXktOh2eWyK5HaMN8GELez+JN72zTtUHcT/fwRm1p9/H+FXsdQZQyAr39skxFSMmDYCZoQoFQn02ATB6oamvnN/PU8v2Q3j2ZW4U3oezDQACSnQVNN7CvSUAnb/2N3uAyle6lviBlptfvssvqTr4LMKG5P0J6kVDjqvM55L6XUYTTYdLGFG0r40aurKalu4OunjuKc/ZBQ1eoDODnNrh0Wa5sX2H103Hah+fUdDevmdVzO74u5dsHNE6Kw4KZSqlvQYNNFKuqa+OUb63j1yz2MG5DOX64/kSlDs2Fu8eFbCSSl2a0GXG6gFrYN8+2KAP51xdzKGW3XNasv73gS6IEFNy+wWxQopXoEDTZd4O01e/nJv9ZQXtfEt2aP4c7ZY0hJTLAvVhdD/2MOPSE5zT421cZuULul0bZsJl4CnoTQzj2wRto2GDKt/bLL/26D0knfCa+eSqluSYNNJ9pX08jP561l/qoiJuRl8tTN05k4OCATy+e1qc8ZrZah6Yxgs+NjOx/lqAtCP9ef/rx/a/vBxtsCn/0Jhs06uHWzUqpH0GATqVUvQWIyTJjTZhFjDK8t38Ov3lhHTWML3ztrHHecNpqkhFYZUXVldn+Y9FbdaIHBJlY2zHf2rjm147Kt5YwEpOOtBta+BpW74Pz7w6qiUqr70mATqcV/htTsNoNNQXkdP35tDR9tKmXqsGzuu2wyYwe0sa5Z4HbQgZLT7WOsMtLc7F3TnsQUu2Jy6ca2yxhjl6bJPRrGnhN+XZVS3ZIGm0h5Em1rJIiPNpXyjb8vwwC/uGgC188aQUJ7y83UOBuOHtaycXbfjFXLxs3eNR3JmwxFK9t+fesHULwa5jyqc1yU6oE02ERKEuxYSyuvLCvgh6+sYuyADOZeP83dds01bbVsnG60WE3s3DDf3kd7e9d0JG+K3QOnrYy0Tx6yC27qTphK9Uj6J2akPIcGG2MMj364he+9tJLjR+Xw4tdnugs0cLAb7bCWTYy70dzuXdOeQcfZx2Ctm12LbQLCCd+241tKqR5Hg02kPAngaznw7Z8/2sr9b2/k4mMH8bebZpCRGsK8mJpiu05Y610lY5kgsG8z7NsIR4eRhRbIH2z2LDv8tY8fgF45MO3GyN5DKdVtxW2wEZFzRWSjiGwRkXuCvJ4iIi84ry8WkREBr/3IOb5RRGI7Gi0JB8Zs/rV8D/e/vZE5Uwbxh6umkJwY4o+3ei+kB9l9M5bBZsN8++hm75r29M6BARPt8v2BClfA5nfsUv7++1BK9ThxGWxEJAF4BDgPmABcIyITWhW7BSg3xowBHgLuc86dAFwNHAOcCzzk2kyVAAAe20lEQVTqXC82nG60jzaV8oOXVzJzVA73Xz45vH1naoqDB5skf7CJQTfahvl2a4DsoZFfa+xZsOtzu8Ya2Ay0BT+1rZoZt0V+faVUtxWXwQaYAWwxxmwzxjQBzwOtc4vnAE87z18GzhARcY4/b4xpNMZsB7Y414sJnySwr6qO255Zytj+GTx2ff7B1QBCVb338KVqwI5zeJKgKcoJAtXFULAksiy0QEddYLsUV79sv18/zy7sedo9upy/Uj1cvGajDQZ2B3xfABzfVhljTIuIVAJ9neOLWp0bdMtJEbkduB1g2LBhYVV0x/4GmqrrOWVsLr+7fHL4WzYb03bLBmKzGOemtwAT+XiN35B8GDQVPn7Q3se8b9tWU/7XonN9pVS3Fa8tm2B9UMZlGTfn2oPGzDXG5Btj8nNzc0OsolXVBL09zTx+Yz590iLItGqsgpaG4C0bsBlp0Q42G+ZD9nDo37qHMkwicP7voW4fvPBVW+crnort4qFKqW4hXls2BUDgIMIQoLCNMgUikghkAftdnhs1hfRnIv+xLZNIthiutttBH5b27BftPW0aq2HbhzD9tuhujTxkGnxjERQuh9GztftMKQXEb8tmCTBWREaKSDJ2wL/1hinzAH8u7eXAB8YY4xy/2slWGwmMBb6IVUX3tqSTiDfybZvbmtDpF+1utC3v220LotWFFihnJEy8VAONUuqAuGzZOGMwdwLvAAnAk8aYtSLyS2CpMWYe8ATwrIhswbZornbOXSsiLwLrgBbgm8a0sZ5MFFw08xj4ADtzPjUz/Au5atlEMdhsmG+zxIa2HgpTSqnoi8tgA2CMeRN4s9WxnwU8bwCCrn1ijPk18OuYVtCR298JDvXl0Gd4+Bfyt2zS+wd/PTkNqovCv34gbzNsegfGXwgJcftPQCl1BHH1SSMiOS6K+YwxFRHWp/tJdGb7e5siu071Xnut1Kzgr0ezZbPjE2isjE0XmlJKBeH2z9pC56u9keQEILz84e4swclAa2mM7Dr+tOe2BuujGWw2vgmJvWDU6dG5nlJKdcBtsFlvjDmuvQIisjwK9el+/MEmGi2bttKeIXqpz8bY8ZrRsw9uXaCUUjHmNhttVpTKHHn8c0i8zZFdp6ak7QmdcDD12QSdMuRe0Qqo2qNdaEqpTuUq2DiD8RGXOSJFq2VT01HLJg2Mz078jMSGN0E8MO7cyK6jlFIhCGuejYjMFJEPRORTEflKtCvVrSSm2MdIgk1TnV28st2WjbOnTWOEEzs3zIdhsyCtb2TXUUqpELgKNiLS+k/uu4GLsasq/yralepWDnSjRRBs/CnNWUPaLnNgA7UIJo/u3w4la7ULTSnV6dwmCPxFRJYBv3O6yyqAawEfUBWrynUL0ehGq9pjHzPy2i6TEoWWzUZn2lKke9copVSI3I7ZfAVYAbwhItcD38EGmt5Az+5GOxBsIkgQqHKWbssMuji1FY2toTfMh/7H2OVklFKqE7keszHG/Bs4B8gGXgU2GmMeNsaUxqpy3YK/Gy2SeTb+lk1mey2bDPsYbsumtsxubKZdaEqpLuB2zOZiEfkEuwrYGuw6ZJeIyD9FZHQsKxj3PE6w8bWEf42qQkjNbn/b5EjHbDa9bbPZNNgopbqA2zGb/8XOo+kFvGmMmQHcLSJjsWuQXR2j+sU/j7MrZ0TBpqj9LjSIfMxm09uQMchuZqaUUp3MbbCpxAaUXkCJ/6AxZjM9OdAAeJwfofGFf42qPZA5qP0ykYzZeJvt3jXHXBLdvWuUUsolt2M2l2CTAVqwWWjKT5wfYaTdaG6DTTgtm92L7U6gY88K/VyllIoCty2bd40xU9srICJfdlTmiCRiWzfhBpuWJqgt6TjYJCTaxTPDGbPZ/K4dWxp5anh1VEqpCLkNNuNFZFU7rwt2W+aeKZJg45/Q2VGwATtuE07LZvMCGD4rss3dlFIqAm6DzdEuysRsN8y450kEX5i3f2COjYtgk5we+phNZQGUrIOzevZCD0qprtVhsBGRy4A9xphFnVCf7kkSwm/ZHJhj00E2GoTXstm8wD6OPTu085RSKorctGxuABJE5MAWzcaYR2NVIWdX0BeAEcAO4EpjTHmrMlOAPwOZ2BbVr40xLzivPQWcis2gA7jJGLMiVvUFbPpzp7RsMkJv2Wx5D7KGQe5RoddNKaWixE022neB9cDagK9Yugd43xgzFnjf+b61OuAGY8wx2MVA/yAi2QGv/8AYM8X5im2ggcjHbJLTIcXFeEpKOjSGkCDQ0mhTnseeqSnPSqku1WHLxhizDfhBJ9TFbw5wmvP8aeBD4Iet6rQp4HmhiJQAudgFQjtfJMGmao9dgNNNMEhOh6Zt7q+963PbEtIuNKVUF3O7XM1dIvJX5/lPY1slBhhjigCcx/4d1G0GkAxsDTj8axFZJSIPiUhK7KrqiDRBwE0XGoQ+ZrN5gV0odOQp4dVNKaWixO2kztHAbud5RqRvKiLviciaIF9zQrxOHvAscLMxB6bw/wibPTcdyKFVq6jV+beLyFIRWVpaGsF6oh5PBC2bQnfJARD6mM3mBTD8xPbXXFNKqU7gNvXZAL1EZCLg8s/wdi5mzJltvSYixSKSZ4wpcoJJSRvlMoH5wE8CM+X8rSKgUUT+Bny/nXrMBeYC5Ofnm9DvxBFuN5q32Y7ZZLkMNilOsPF5D67J1pbynbBvI0y7KfR6KaVUlHXYshERAcqwEzevx7YcYmkecKPz/Ebg9SB1SgZeA54xxrzU6rU851Gwe+2siWltwQYbE0Y3WlWhXVMta6i78r2cHIhGF/vVbdGUZ6VU/Ogw2BhjDDAWWAn8B5gEICInicg3RWSUv6yIRGNXrnuBs0RkM3CW8z0iki8ijztlrgROAW4SkRXO1xTntedEZDWwGuiHXbE6tsIds6l0eiaz3QabPvaxvrz9cmC70PqMgL49ewcIpVR8cNuN9h52EL4ftksNbPbXDGCGiJQB/wC+B1wTSYWMMWXAGUGOLwVudZ7/Hfh7G+fPjuT9w+IJc1JnhRNssoa5K+822DQ3wLaPYOr1mvKslIoLroKNMebpIMdeE5F5wDRgInAysKl1uR4h3DGbil32MWuIu/KpTjdaR8Fm56fQUq9daEqpuOEq2DjpznXGmAcCjxtjvMAXzlfPFe5yNZW7IK0/JKW6K3+gZdPBdKLNCyAxFUacFHqdlFIqBtx2o10PTGl9UERuBXKNMb+Naq26m3DHbCp2Q7bLLjRw3422ZQGMOBmSeoVeJ6WUigG382zqjTF1QY4/C1wXxfp0T+GujVa5231yABzMRmsv2JRthbIt2oWmlIorroONP6U4kDGmEbt7Z88WzpiNz2eX/3eb9gyQkGQndrYXbLa8Zx/HtjmVSSmlOp3bYPMA8LqIDA88KCL9AV/wU3qQXYugIMRhq9oS8DaF1o0GtiutvWCzeQH0HQM5o9ouo5RSncxtNtpLItIbWCYii4AV2EB1BfCL2FWvm2ipt48NVe53wzyQ9hxCywZsV1pbwaapDnZ8DNNuDu2aSikVY25bNv7051HAi0AS0ABcY4x5LkZ1635aGtyXrdhpH0MZswFI6we1bazjtuMTW4exZ4V2TaWUijG3qc+zgEXGmCrgmdhWqRu6+E8w787Qgk1lmC2bjDwobWM605YFkNTbLr6plFJxxG3L5kZsF9rzInKTiAyMZaW6nURnnkxLk/tzKnZDapb7bje/9AFQs9cmGAQyBja9AyNPdT9vRymlOonbMZs7AETkaOA84CkRyQIWAm8DnzoTPHumxGT7GGrLJtTkALAtG18L1JVBeu7B46UbbdfcSd8N/ZpKKRVjrsdsAIwxG4wxDwGXAecAn2CTBBbHoG7dh8eJ2aHE24rd7tdEC5ThNCqriw49vvkd+6jza5RSccjtTp0eEblWROaLSDGwAdgC/A4nUSCGdYx/4uwt43aujTGhT+j0y3CmO1XvPfT4pndg4CT3e+MopVQnctuyWYjdrfNHQJ4xZqgxpj928c1FwG9FpOeuJOBv2bQeR2lLfbndBC3U5AA4uGhn+Y5Dr7drEYw7N/TrKaVUJ3C7NtqZxpjm1geNMfuBV4BXRCQpqjXrTjxOzHbbsvGv9hxWy2agXUVgX0BG2qZ3bBeeBhulVJxys1PnZdhtBNoVLBj1GP5uNLdjNuGmPYPdnyZ33KHBZvXLdvxncIe/JqWU6hJuWjY3AAki8qb/gDHm0dhVqRs60I3mtmXj36EzjAQBgH7jYNuH9nntPtj6AZzwLd0oTSkVt9yM2XwXWA+sDfhSgTz+BIEQWjZJvaF33/DeL+9Ym41WvhOWP2tbVJOvCu9aSinVCTps2RhjtgE/6IS6dF8HUp9dJghU7LJdaOG2REadZh/XvQ6LH7MTOQdMCO9aSinVCdymPt8lIn91nv80lhUSkRwRWSAim53HPm2U84rICudrXsDxkSKy2Dn/BRFJjmV97ZuGmCAQbtqzX+7RkDcFFvzUtnBO+1H411JKqU7gNvV5NOAMNJARo7r43QO8b4wZC7zvfB9MvTFmivN1ccDx+4CHnPPLgVtiW10CxmxcdqP5WzbhEoE5f4Jx59l12YbPCv9aSinVCdwGGwP0EpGJwKAY1gdgDvC08/xp4CtuTxQRAWYDL4dzftg8IUzqbKyx82IiadmAncB57fMw9frIrqOUUp3ATeqzAGWAANdjJ3bG0gBjTBGA89i/jXKpIrJURBaJiD+g9AUqjDH+T/0CoM0p9SJyu3ONpaWlbSzb70Yoy9UcmGMzvP1ySil1BHGTIGBEZCzwJlAFTAJ2i8hJwLHAW04SASIy0hizvaNrish7QLCVo38cQt2HGWMKRWQU8IGIrHbqd9gttHUBY8xcYC5Afn5+m+U6JCFko2mwUUr1QG5XEHgPSAb6cfDDOxeYAcwQkTLgH8D3cLFOmjHmzLZeE5FiEckzxhSJSB5Q0sY1Cp3HbSLyIXAcdjWDbBFJdFo3Q4BCd7cYgQMrCLgINpURzrFRSqluyO0WA08HOfaakwU2DZiIXSetjV29QjIPu3/Ovc7j660LOBlqdcaYRhHpB5wI3O+0whYClwPPt3V+1IUyqbNip93/Jr2t3kGllDryuN2p86fYD/cHAo87e9h84XxFy73AiyJyC7ALu4UBIpIP3GGMuRUYDzwmIj7suNO9xph1zvk/BJ4Xkf8FlgNPRLFuwYWyXE2kc2yUUqobctuNdj0wpfVBEbkVyDXG/DZaFTLGlAFnBDm+FLjVef4Zduwo2PnbsN17nSeU1OeKXdqFppTqcdymPtcbY+qCHH8W6LlbC/iFslyNBhulVA/kOtg4g/WHMMY0Ai6nzR/BPC670Rpr7HbOGmyUUj2M22DzAPC6iBySrysi/QGXC4Idwdzu1KmZaEqpHsptNtpLItIbWCYii4AV2EB1BfCL2FWvm3A7ZqNzbJRSPZTblo0//Xkk8CKQBDQA1xhjnotR3boPt8vVHAg22rJRSvUsbrPRADDGVAPPxKgu3deB1OcOehR1jo1Sqody3bJR7QilZaNzbJRSPVBYwUZELop2Rbo1EbunjZsxG+1CU0r1QOG2bH4d1VocCTyJHac+a7BRSvVQ4QYb7QdqTRLa70bTOTZKqR4s3GAT/nL8RypPIvjaSRDQOTZKqR5MEwSixeNpv2VTWWAfNdgopXogDTbR0tGYjX+OTdaQzqmPUkrFkXCDTXFUa3Ek6GjMprLABqT0AZ1XJ6WUihNhBRtjzFnRrki350loP/W5sgAyBx2ck6OUUj2IdqNFiyex42CTNbTz6qOUUnFEg020iKf9MRsNNkqpHkyDTbS017LxtkDVHk0OUEr1WOEuV3OSiHxTREYFHBsZjQqJSI6ILBCRzc5jnyBlTheRFQFfDSLyFee1p0Rke8Brh21nHROedhIEavbaVo8GG6VUDxVuyyYXmAH8XEQeFJF84DdRqtM9wPvGmLHA+873hzDGLDTGTDHGTAFmA3XAuwFFfuB/3RizIkr1al97qc/+OTbajaaU6qFC2mLAzxjzmojMA6YBE4GTgU1RqtMc4DTn+dPAh8AP2yl/OfCWMaYuSu8fHmknG+1AsNGWjVKqZ3IVbETkp0CdMeYB/zFjjBf4wvmKpgHGmCLnPYqcrafbczXwYKtjvxaRn+G0jIwxjVGu4+HaS33WCZ1KqR7ObcvmeuCwsQ8RuRXINcb8NpQ3FZH3gIFBXvpxiNfJAyYB7wQc/hGwF0gG5mJbRb9s4/zbgdsBhg2LcBmZ9sZsKgugVx9ISY/sPZRSqptyG2zq2+imehb4Eggp2BhjzmzrNREpFpE8p1WTB5S0c6krgdeMMc0B1y5ynjaKyN+A77dTj7nYgER+fn5ki4tKQvtjNtqqUUr1YG4TBOqdD/5DON1THWxPGbJ5wI3O8xuB19spew3wz8AD/nqKiABfAdZEuX7BtZf6rHNslFI9nNtg8wDwuogMDzzojKe0s65+WO4FzhKRzcBZzveISL6IPB7w3iOAocBHrc5/TkRWA6uBfsD/Rrl+wbU3ZqMtG6VUD+eqG80Y85KI9AaWicgiYAU2UF0B/CKaFTLGlAFnBDm+FLg14PsdwOAg5WZHsz6ueRLA23T48YZKaKzUlo1Sqkdz1bIRETHGPA2MBF4EkoAG4BpjzHP+MjGrZXfQ1qrPmvaslFKuEwQWisgrwOvGmGf8B0UkWURmY8dWFgJPRb+K3URbYzY6oVMppVwHm3OBrwH/dJalqQBSgQTszP2HOm2mfrxqa8zGvx20tmyUUj2Y2zGbBuBR4FERScIOvNcbYypiWbluxdNG6nNlAXiSdNM0pVSPFvJyNc6clqIOC/Y0bY3ZVOyGrMHg0QW2lVI9l34CRkub3Wg6x0YppTTYREtbqz7rHBullNJgEzXBVn32tkB1oQYbpVSPp8EmWoJ1o1UXgvFpN5pSqsfTYBMtwVZ91gmdSikFaLCJnmCrPuuETqWUAjTYRE+wbrQDEzoPW8JNKaV6FA020SIJYFptiVNZAL1yIDmta+qklFJxQoNNtIgneDeajtcopZQGm6jxeIJkoxVB5qCuqY9SSsURDTbREixBoHovZAzsmvoopVQc0WATLZ4EO6fGz9sMtaWQrsFGKaU02ESLtOpGqymxj9qyUUopDTZRIwmAOZiRVr3XPmbkdVmVlFIqXsRdsBGRK0RkrYj4RCS/nXLnishGEdkiIvcEHB8pIotFZLOIvCAiyZ1ScU+CffS3bmr8wUb3sVFKqbgLNsAa4FLgP20VEJEE4BHgPGACcI2ITHBevg+7c+hYoBy4JbbV9VfK+VH6x22qnS1/tGWjlFLxF2yMMeuNMRs7KDYD2GKM2WaMaQKeB+aIiACzgZedck8DX4ldbQMcCDZOy6Z6rz2Wltspb6+UUvEs7oKNS4OB3QHfFzjH+gIVxpiWVseDEpHbRWSpiCwtLS2NrEaHdaMV20DjP66UUj1YyNtCR4OIvAcES9P6sTHmdTeXCHLMtHM8KGPMXGAuQH5+fpvlXBEnqPhbNjUlkN4/oksqpdSRokuCjTHmzAgvUQAELqU8BCgE9gHZIpLotG78x2PP34Lxj9nUlEC6JgcopRR03260JcBYJ/MsGbgamGeMMcBC4HKn3I2Am5ZS5PxjNr6AYJOmLRullII4DDYicomIFACzgPki8o5zfJCIvAngtFruBN4B1gMvGmPWOpf4IXC3iGzBjuE80TkVD0gQMAZqtRtNKaX8uqQbrT3GmNeA14IcLwTOD/j+TeDNIOW2YbPVOldggkBDBXibNNgopZQj7lo23ZYEjNn4l6rRMRullAI02ERPYDeaP9joHBullAI02ERPYDdaTbF9ri0bpZQCNNhET2A3Wq0zQVTHbJRSCtBgEz2B82xqisGTBKnZXVsnpZSKExpsokWcxQt8zphNWq7dKloppZQGm6gJXK5Gl6pRSqlDaLCJlsAEgdoSzURTSqkAGmyiJXA/m7py6N23a+ujlFJxRINNtAR2o9WXQ68+XVsfpZSKIxpsosXfjdbSCE3V0Duna+ujlFJxRINNtPi70Wr32Udt2Sil1AEabKLFH2zqNNgopVRrGmyixd+NVltmHzXYKKXUARpsosWfIHCgZaOrByillJ8Gm2jxt2zq9ttHXapGKaUO0GATLf4xm/py+6jdaEopdYAGm2jxd6PVOy2blMyuq4tSSsWZuAs2InKFiKwVEZ+I5LdRZqiILBSR9U7ZuwJe+4WI7BGRFc7X+cGuEXX+RTfr9kNyOiTE3Y7bSinVZeLxE3ENcCnwWDtlWoDvGWO+FJEMYJmILDDGrHNef8gY8/tYV/QQgS0bbdUopdQh4i7YGGPWA4h/yf7gZYqAIud5tYisBwYD69o8Kdb8YzYNldBvXJdVQyml4lHcdaOFSkRGAMcBiwMO3ykiq0TkSRHpnJF6fzYawOgzOuUtlVKqu+iSYCMi74nImiBfc0K8TjrwCvAdY0yVc/jPwGhgCrb180A7598uIktFZGlpaWmYd+O/WECwyR4W2bWUUuoI0yXdaMaYMyO9hogkYQPNc8aYVwOuXRxQ5q/AG+3UYy4wFyA/P99EVKHAlo3uZaOUUofolt1oYgd0ngDWG2MebPVaXsC3l2ATDjqjUgefp/XrlLdUSqnuIu6CjYhcIiIFwCxgvoi84xwfJCJvOsVOBK4HZgdJcb5fRFaLyCrgdOC7nVNxbdkopVRb4jEb7TXgtSDHC4HzneefAEHT1Ywx18e0gm3RbjSllGpT3LVsui1PQNxOSu26eiilVBzSYBMtiRpglFKqLRpsoiWpV1fXQCml4pYGm2hJSO7qGiilVNyKuwSBbksEzv899B3T1TVRSqm4o8Emmmbc1tU1UEqpuKTdaEoppWJOg41SSqmY02CjlFIq5jTYKKWUijkNNkoppWJOg41SSqmY02CjlFIq5jTYKKWUijkxJrINKo8UIlIK7Azz9H7AvihWpzvQez7y9bT7Bb3ncAw3xnS4r4oGmygQkaXGmPyurkdn0ns+8vW0+wW951jSbjSllFIxp8FGKaVUzGmwiY65XV2BLqD3fOTrafcLes8xo2M2SimlYk5bNkoppWJOg00ERORcEdkoIltE5J6urk80icgOEVktIitEZKlzLEdEFojIZuexj3NcRORh5+ewSkSmdm3t3RGRJ0WkRETWBBwL+R5F5Ean/GYRubEr7sWtNu75FyKyx/ldrxCR8wNe+5FzzxtF5JyA493m376IDBWRhSKyXkTWishdzvEj8nfdzv127e/ZGKNfYXwBCcBWYBSQDKwEJnR1vaJ4fzuAfq2O3Q/c4zy/B7jPeX4+8BYgwExgcVfX3+U9ngJMBdaEe49ADrDNeezjPO/T1fcW4j3/Avh+kLITnH/XKcBI5997Qnf7tw/kAVOd5xnAJufejsjfdTv326W/Z23ZhG8GsMUYs80Y0wQ8D8zp4jrF2hzgaef508BXAo4/Y6xFQLaI5HVFBUNhjPkPsL/V4VDv8RxggTFmvzGmHFgAnBv72oenjXtuyxzgeWNMozFmO7AF++++W/3bN8YUGWO+dJ5XA+uBwRyhv+t27rctnfJ71mATvsHA7oDvC2j/F9rdGOBdEVkmIrc7xwYYY4rA/oMG+jvHj6SfRaj3eKTc+51Ol9GT/u4kjsB7FpERwHHAYnrA77rV/UIX/p412IRPghw7klL7TjTGTAXOA74pIqe0U/ZI/1lA2/d4JNz7n4HRwBSgCHjAOX5E3bOIpAOvAN8xxlS1VzTIsW5330Hut0t/zxpswlcADA34fghQ2EV1iTpjTKHzWAK8hm1SF/u7x5zHEqf4kfSzCPUeu/29G2OKjTFeY4wP+Cv2dw1H0D2LSBL2g/c5Y8yrzuEj9ncd7H67+veswSZ8S4CxIjJSRJKBq4F5XVynqBCRNBHJ8D8HzgbWYO/Pn4FzI/C683wecIOTxTMTqPR3T3RDod7jO8DZItLH6ZY42znWbbQaX7sE+7sGe89Xi0iKiIwExgJf0M3+7YuIAE8A640xDwa8dET+rtu63y7/PXd15kR3/sJmrWzCZmz8uKvrE8X7GoXNPFkJrPXfG9AXeB/Y7DzmOMcFeMT5OawG8rv6Hlze5z+x3QnN2L/ibgnnHoGvYQdVtwA3d/V9hXHPzzr3tMr5MMkLKP9j5543AucFHO82//aBk7DdP6uAFc7X+Ufq77qd++3S37OuIKCUUirmtBtNKaVUzGmwUUopFXMabJRSSsWcBhullFIxp8FGKaVUzGmwUUopFXMabJRSSsWcBhul4pizmsOfnJnsSnVbGmyUim93YPcZOamrK6JUJDTYKBXfzsUuF7KiqyuiVCQ02CgVp0QkFbtb4lTgoy6ujlIR0WCjVPwaiw02G4wxzV1dGaUikdjVFVBKtSkXGEccb7mslFvaslEqfg3CboDlCdjCV6luSYONUnFIRBKxYzUDgb8A3q6tkVKR0f1slFJKxZy2bJRSSsWcBhullFIxp8FGKaVUzGmwUUopFXMabJRSSsWcBhullFIxp8FGKaVUzGmwUUopFXP/H5nYb78fetXFAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEMCAYAAABtKgnyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNX5+PHPk8lKSAKBBMK+L7FYFMSl7lYFN6xL1bpgq7W22tZqF/1aqz+7aTf7bWvbL1brXqVaFRVF3JeKAoossu+BAGEPCVnn+f1x7oRhmCSTydxMEp736zWvmbn33DvnJjBPzjnPPUdUFWOMMSYZUpJdAWOMMYcuC0LGGGOSxoKQMcaYpLEgZIwxJmksCBljjEkaC0LGGGOSxoKQMcaYpLEgZIwxJmksCBljjEma1GRXoL3r2bOnDho0KNnVMMaYDmXevHnbVLWguXIWhJoxaNAg5s6dm+xqGGNMhyIi62IpZ91xxhhjksaCkDHGmKSxIGSMMSZpLAgZY4xJGgtCxhhjksaCkDHGmKSxIGSMMSZpLAi1hU2fQsm8ZNfCGGPaHbtZtS1MPdk937U7qdUwxpj2xlpCxhhjksaCkDHGmKSxIGSMMSZpLAgZY4xJGgtCxhhjksaCkDHGmKSxINSWypYnuwbGGNOutNsgJCITRWSZiKwUkVuj7D9RRD4RkToRuShs+1gR+VBEFovIAhG5JGzfwyKyRkTme4+xbXU9ANx/VJt+nDHGtHft8mZVEQkA9wOnAyXAHBGZrqqfhxVbD1wN/DDi8ErgKlVdISJ9gHkiMlNVd3n7f6Sqz/h7BU1QBZGkfbwxxrQn7TIIAROAlaq6GkBEngImAw1BSFXXevuC4Qeq6vKw15tEZCtQAOyiPdAgSCDZtTDGmHahvXbH9QU2hL0v8ba1iIhMANKBVWGbf+l1090nIhmtq2YcgnVt/pHGGNNetdcgFK2/Slt0ApEi4DHg66oaai3dBowCjgLygZ80cux1IjJXROaWlZW15GObF6xP7PmMMaYDa69BqAToH/a+H7Ap1oNFJBd4Gfipqs4ObVfVUnWqgX/iuv0OoqpTVXW8qo4vKCiI6wIaZS0hY4xp0F6D0BxguIgMFpF04FJgeiwHeuWfAx5V1X9H7CvyngU4H1iU0FrHwoKQMcY0aJdBSFXrgBuBmcASYJqqLhaRu0XkPAAROUpESoCLgf8TkcXe4V8FTgSujpKK/YSILAQWAj2BX/h+MXU1B7637jhjjGnQXrPjUNUZwIyIbT8Lez0H100XedzjwOONnPPUBFezeXP+EVEJC0LGGBPSLltCnUp56YHvrTvOGGMaWBDyW33tge8tCBljTAMLQm3NxoSMMaaBBaG2ZkHIGGMaWBDyW0bOge+tO84YYxpYEPJb/pAD31sQMsaYBhaEfOfNNtTnSO+tdccZY0yIBSG/qReEjrjcPduYkDHGNLAg5LfQ3KmBdPccmbJtjDGHMAtCvvNaQilp3ltrCRljTIgFIb+FuuNSvBmSrDvOGGMaWBDyW0N3nNcSsuw4Y4xpYEHId15LqCEIWUvIGGNCLAj5TSPGhKwlZIwxDSwI+S4UhALeW2sJGWNMiAUhv4VaQqEUbWsJGWNMAwtCflMbEzLGmMa02yAkIhNFZJmIrBSRW6PsP1FEPhGROhG5KGLfFBFZ4T2mhG0fJyILvXP+SUTE/yuJSNF+9hqoqfT/Y40xpgNol0FIRALA/cAkoBi4TESKI4qtB64Gnow4Nh+4EzgamADcKSLdvd1/A64DhnuPiT5dwn6hFO2UsJXUd2/w/WONMaYjaJdBCBc8VqrqalWtAZ4CJocXUNW1qroACEYceyYwS1V3qOpOYBYwUUSKgFxV/VBVFXgUON/3K4kcEwIbFzLGGE97DUJ9gfDmQom3rTXH9vVex3POVogYEwKo3uv/xxpjTAfQXoNQtLEabeWxMZ9TRK4TkbkiMresrCzGj22ERqRoA9SUt+6cxhjTSbTXIFQC9A973w/Y1MpjS7zXzZ5TVaeq6nhVHV9QUBBzpaNqGBMKbwlZEDLGGGi/QWgOMFxEBotIOnApMD3GY2cCZ4hIdy8h4QxgpqqWAuUicoyXFXcV8IIflT9QlDEh644zxhignQYhVa0DbsQFlCXANFVdLCJ3i8h5ACJylIiUABcD/ycii71jdwA/xwWyOcDd3jaAbwP/AFYCq4BX2uBi3HP4mFCNBSFjjAFIbb5IcqjqDGBGxLafhb2ew4Hda+HlHgIeirJ9LvCFxNa0GaHuuNTM/dusJWSMMUA7bQl1LmGJCQWj3WtLTDDGGMCCkP9C3XEI3DAbsvKhpiKpVTLGmPbCgpDvvCAUmiEoo6t1xxljjMeCkN/CW0IAEoAFT8HSl5NWJWOMaS8sCPktFITE+1H3GOqen70W6muTUydjjGknLAj5LqI77ry/wDHfgdpKKJmTvGoZY0w7YEHIb6EU7VAQyi2C477nXm9emJw6GWNMO2FByG8aZXq6nN6Q1R22ft729THGmHbEgpDvdP94UIgI5A+BneuSUyVjjGknLAj5TRuZwDuvH+wuOXi7McYcQiwI+U2D+8eDwuX1hz0bo3fXGWPMIcKCkO8aaQnl9nUZcvt2tnmNjDGmvbAg5DeNMiYErjsOYNf6tq2PMca0IxaE/NZYd1w3b9293RsO3meMMYeImJZyEJH8GIoFVXVXK+vTCTXSHdd9kHveubYN62KMMe1LrOsJbfIeUb5NGwSAAa2uUWfTWHdcVnfIzLM0bWPMIS3WILREVY9oqoCIfJqA+nRO0brjALoNtJaQMeaQFuuY0LEJKhMzEZkoIstEZKWI3Bplf4aIPO3t/0hEBnnbLxeR+WGPoIiM9fa97Z0ztK8wkXWOSoM02oDsPgh2rPa9CsYY017FFIRUtSoRZWIlIgHgfmASUAxcJiLFEcWuAXaq6jDgPuBerx5PqOpYVR0LXAmsVdX5YcddHtqvqlsTVedGqTbeEiosdkHIFrkzxhyi4sqOE5FjRORNEflARM5PdKWACcBKVV2tqjXAU8DkiDKTgUe8188Ap4kc9G1/GfAvH+rXAk0Eod5j3P6tS9q0RsYY017EFIREpHfEppuB84CJwM8TXSmgLxCeu1zibYtaRlXrgN1Aj4gyl3BwEPqn1xV3R5SgBYCIXCcic0VkbllZWbzX4DTVHdd7jHvevKB1n2GMMR1UrC2hv3tf2pne+13A13Bf8nt8qFe0b+3I+W2aLCMiRwOVqroobP/lqjoGOMF7XBntw1V1qqqOV9XxBQUFLav5wSdrIjFhAGTkweZF0fcbY0wnF+uY0PnAfOAlEbkSuAkIAl0AP7rjSoD+Ye/74VLEo5YRkVQgD9gRtv9SIlpBqrrRey4HnsR1+/mskfuEwAWngpGwbbn/1TDGmHYo5jEhVX0ROBPoBvwHWKaqf1LVVvZXRTUHGC4ig0UkHRdQpkeUmQ5M8V5fBLyp6mYDFZEU4GLcWBLetlQR6em9TgPOAfxvgjTVEgLI6QUVfvwIjTGm/Yt1TOg8EXkfeBP3xX0p8BUR+ZeIDE10pbwxnhuBmcASYJqqLhaRu0XkPK/Yg0APEVmJG6MKT+M+EShR1fD85wxgpogswLXqNgIPJLruB9F6kEDj+7v2gr1bfK+GMca0R7HerPoL3H1AWcAMVZ0A3Cwiw4Ff4oJSQqnqDGBGxLafhb2uwrV2oh37NnBMxLYKYFyi69msYD2kNPFjzi50M2nX1UBqetvVyxhj2oFYg9BuXKDJAhrurVHVFfgQgDqVYD2kNNUS8u6XrSiDvMgEQGOM6dxiHRP6Ci4JoQ6XFWdipfXR544LaQhC/t83a4wx7U2sLaHXVPXIpgqIyCfNlTkkBeua744D2GvJCcaYQ0+sQWi0N6DfGMGlSJtIzXbHefchWUvIGHMIijUIjYqhTH1rKtJpaXOJCV4Q2mtByBhz6Gk2CInIhcBGVZ3dBvXpfILNpGinZ0NaNlRsa7s6GWNMOxFLS+gqICAiDenSqvpX/6rUyQTrIaWZ/I+uBdYdZ4w5JMUShH4AfBtY7HNdOqf6agg0c/9PdqF1xxljDknNBiFv1oEftUFdOqfafZDWpekyXQttcTtjzCEp1ml7vi8iD3iv7/C3Sp1MTaUb92lKdoHNH2eMOSTFerPqUPav75PjU106p/pqCKQ1XSa7ACq3u/EjY4w5hMQahBTIEpEvAH18rE/n09zcceC64zToApExxhxCmg1C3uqj23E3pF4J3OZ3pTqV5mbRhv1T95SX+l8fY4xpR5oNQt4aPcOBz4B3gTEAInK8iNwgIkNCZUVksF8V7bCCwaZnTADo7v3Ydqzxvz7GGNOOxDpjwutAOtCT/UtoF+BWJp0gIttxK5XeAlyW6Ep2aNrMtD0A+aEgtMr/+hhjTDsSUxBS1UeibHtORKbj1uj5AnACYOtURwrWNd8dl5HjFrezNG1jzCEm1hTtO0Tklsjtqlqvqh+r6kOqep+q3pmoionIRBFZJiIrReTWKPszRORpb/9HIjLI2z5IRPaJyHzv8fewY8aJyELvmD95413+am4C05D8obDdgpAx5tASa3bclcDfIjeKyLUikvBEBREJAPcDk4Bi4DIRKY4odg2wU1WHAfcB94btW6WqY73H9WHb/wZchxvjGg5MTHTdDxJLYgJA/hDrjjPGHHJiDUL7VLUyyvbHgCsSWJ+QCcBKVV2tqjXAU8DkiDKTgVA34TPAaU21bESkCMhV1Q+9ZItHgfMTX/UIsSQmAPQYAnu3QPVe36tkjDHtRcxByPsSP4CqVuNWW020vuy/ORagxNsWtYyq1uGWIO/h7RssIp+KyDsickJY+ZJmzpl4zS1qF5LvJRnauJAx5hASaxD6PfCCiAwM3ygihUAw4bVy9yRF0hjLlAIDVPUI4GbgSRHJjfGc7sQi14nIXBGZW1bWyul0mlveOyR/qHu2LjljzCEk1uy4f4tIF2CeiMwG5uMC2MXAXT7UqwToH/a+H7CpkTIlIpKKW9l1h9fVVu3Ve56IrAJGeOX7NXNOvOOmAlMBxo8fHzVQxSzmxIRQmra1hIwxh45YW0KhNO0hwDQgDagCLlPVJ3yo1xxguIgMFpF04FJgekSZ6cAU7/VFwJuqqiJS4CU24N1IOxxYraqlQLmIHOONHV0FvOBD3Q8Ua2JCKE3bMuSMMYeQmFpCInIsMFtV9+AG9H2lqnUiciMwEwgAD6nqYhG5G5irqtOBB4HHRGQlsAMXqABOBO4WkTrckuPXq+oOb9+3gYeBLOAV7+Hnhbg54WIZEwLXJdfa7rjNi+CF70DBaDjvz5DazFpGxhiTRLHOmDAFuF9ElgOvAq+q6mb/qgWqOgOYEbHtZ2Gvq3DdgZHHPQs828g55+JurG0boVmxY+mOA5ecsHJW6z7zlR/D5oVQ+hl0GwCn3t668xljjI9i6o5T1etV9Ujc+E934GER+VBEfiUiJ4a6v0wE9YJQLIkJAN0HujTt2n3xfd6eTbDuAzjlf+CwC+DD+6FyR/PHGWNMksQ8JgSgqktV9T7gQuBM4H1ca+QjH+rW8bW0JdTNSz7ctaHpco1Z8557HjERTvox1FbAJwfNuGSMMe1GrNP2pIjI10TkZRHZAiwFVgK/xUtQ8LGOHVeoJRTrmFC3Ae5517r4Pm/zAghkQMEoKBwNA46F+U+6sSljjGmHYm0JvYVbXfU2oEhV+6tqIW7S0tnAr0XEj5kTOragdx9vrL2V3UMtoTiD0JZFLviEVnId+zXYthw2zovvfMYY47NYg9CXVfXnqrpAVRtuTlXVHar6rKpeBDztTxU7sJZ2x3Xt7VoyO+MNQp9Dr7C8i+LJEEiHxc/Fdz5jjPFZLCurXohbrqFJqlqbkBp1JqEEg7Ss2MqnpEC3/vG1hGr3QcVW6D5o/7bMPBh8Eix50brkjDHtUiyDFVcBARFpSJdW1b/6V6VOpNab8zWtS+zHdBsAu9a3/LP2eJM/5PU7cPvoc+DF77uuut5jWn5eY4zxUSzdcT8AlgCLwx4mFm//2j23KAgNjK87breXURcZhEaeBQgseanl5zTGGJ81G4S85RR+pKrvhB5tUbFOITQWk96CINR9IOzbAdXlLfus3d4E4ZFBqGuhy5JbakHIGNP+xJqi/X0RecB7fYe/VeqEUjNjLxtK025paygUhHKjrE4x6mzXHbdzbcvOaYwxPos1O24o+9f3yfGpLp1XRgt+ZN0GueeWjgtVlEFW9+hzxY06yz0ve7Vl5zTGGJ/FGoQUyBKRLwB9fKxP5zLS+/LvdVjsx8R7r1DldujSM/q+/CHuBtZlL7fsnMYY47NYUrQF2I5bFO5K3A2rJhYZOQemTMeiSw+XyNDS7riKbe7Yxow8C9Z+APt2tuy8xhjjo1gSExS3Js9nwLvAGAAROV5EbvDW7MHbNtivinZI9bWxT9kTIuIy5FraHVe5o/kgpPWw4vWWndcYY3wUa3fc60A60NN7ABQAE4A7ReQPIjIe+FXiq9iBBetaHoTAdcm1uDtuG2Q3EYT6joPsQlg2o/EyxhjTxmJd3vugqZhV9TkRmY6bTeELuHnklie2eh1csB5S0lp+XLeBrutM1bWMmqPqjQk1EYRSUmDkRFj8PNTV2GJ3xph2IdaVVe8AKlX19+HbVbUe+Nh7mEjButjnjQvXbQDUlLvxmy75zZev2u0+q6kgBK5L7pNHYd37MPTUltfLGGMSLNbuuCuBv0VuFJFrRcSXRAURmSgiy0RkpYjcGmV/hog87e3/SEQGedtPF5F5IrLQez417Ji3vXPO9x6FftS9QWu64yD2LrnK7e65sey4kCEnQ2oWLPN3VXNjjIlVrEFon6pWRtn+GJDwJRy8lVrvByYBxcBlIlIcUewaYKeqDgPuA+71tm8DzlXVMbhlyR+LOO5yVR3rPbYmuu4HCMaRmAD7F7eLNUMutHpqc62mtCzXAlo6wyY0Nca0CzEHIREpityoqtVAXWKrBLiEh5XelEE1wFPA5Igyk4HQWNUzwGkiIqr6qap6s3myGMgUkQwf6ti8YH2cQaiFi9tV7XbPmd2aLztyEuwpgc0LW14vY4xJsFiD0O+BF0RkYPhGrzsrGP2QVunL/hkaAEq8bVHLqGodsBuIHBS5EPjUC5Yh//S64u7w7oHyT7AOAnEEoaxubhmGWNO0q/e451hmZhgxERDLkjPGtAuxZsf9W0S6APNEZDYwHxfALgbu8qFe0YJDZP9Rk2VE5DBcF90ZYfsvV9WNIpIDPIsb63r0oA8XuQ64DmDAgAEtq3m4eMeEoGWzadfsdc8ZXZsv27UA+k9wQejkg4bajDGmTcXaEgqlaQ8GpgFpQBVwmao+4UO9SoD+Ye/7AZsaKyMiqUAesMN73w94DrhKVVeFXcNG77kceBLX7XcQVZ2qquNVdXxBQUH8V9GaINSSe4WqQ0EoxjnqRp4FpZ/tn/TUGGOSJOYgBO7LW1UfVdWfqOrdqjrXp3rNAYaLyGARSQcuBaZHlJmOSzwAuAh4U1VVRLoBLwO3qeoHocIikioiPb3XacA5wCKf6u/Ut7IltGt9bAkEoWUf0lsQhMCy5IwxSdeiINRWvDGeG4GZuAX1pqnqYhG5W0TO84o9CPQQkZXAzUCob+lGYBhwR0QqdgYwU0QW4LoTNwIP+Hoh8d4nBC4I1VXB3hgS+GrKXep1rONPBSOgxzALQsaYpIvrz3QROVdVX0x0ZcKp6gxgRsS2n4W9rsKNSUUe9wvgF42cdlwi69isYF18MybAgfcK5fRqumx1eWzjQeFGToLZf4eqPZCZG18djTGmleJtCf0yobXorFqVmNCCxe2q90J6S4PQWe4+plVvtLxuxhiTIPEGIX9TmzuLeO8TgpbdK1Rd3rKF8wD6Hw1Z+e7GVWOMSZJ4g5Ddbh+L1owJpWdDdkFsQahmb8uDUErA3TO0YqZbcsIYY5KgXSYmdBrxTtsT0m1AjN1xe1oehMAt+121G9Z/2PJjjTEmASwI+am+DgJxJiaAl6bt05gQwJBTIJBhWXLGmKSJNwhtSWgtOqF563ZQX1cFgVas29N9oLuhNFjfdLl4uuPAZdQNORmWvmwTmhpjkiKuIKSqpye6Ip3NP99fg9TXQGor5k7N6+fGlZq7VyieFO2QkZNca2vrkviON8aYVrDuOJ9kpNSTgrYuCOV6c7buiZyxKEx9rbupNSPOe31GTHTPNqGpMSYJLAj5JAOvCy3QmiDUxz3v2dh4mYYpe+JsCeUWQd9xFoSMMUlhQcgnGeKlPfvdEmrJDNqNGTkJNs6D8s3xn8MYY+IQVxASkeNF5AYRGRK2bXDiqtXxZaR4a/21JjGhSw93fCwtoXgSE0JGnu2eLUvOGNPG4m0JFeCWQbhTRP4gIuOBXyWuWh1fBgloCYm4LrmmWkKhZRxinUE7msLRLh3cuuSMMW0s3uy454BvAPfjlkM4AViewHp1eOkkoCUErkuuySCUgJaQCIw6B1a/vX+pcGOMaQMxBSFvKexbwrepar2qfqyqD6nqfap6pz9V7JgyG8aEMlt3otw+UN7UmFAoCLViTAjgsPOhvgaWz2zdeYwxpgVibQldCfwtcqOIXCsityW2Sp1DeiK642B/d1xjN5O2dFXVxvQdDzl9YPHzrTuPMca0QKxBaJ+qVkbZ/hhwRQLr02kkrDsup49roVRuj76/tSnaISkpUHwerHx9/zmNMcZnMQchESmK3Kiq1RD6tjXh0ryWUDCltWNCzdwrVJOglhBA8flQX21dcsaYNhNrEPo98IKIDAzf6C2bHUx4rdy5J4rIMhFZKSK3RtmfISJPe/s/EpFBYftu87YvE5EzYz1nIoW644KJSEyAxpMTqve4cafWTJQa0v9o6NobPrcuOWNM24hpnQFV/beIdAHmichsYD4ugF0M3JXoSolIAJd5dzpQAswRkemq+nlYsWuAnao6TEQuBe4FLhGRYuBS4DCgD/C6iIzwjmnunAnTq2IZAPUp6fGtoR7SXEso3hm0owl1yX3yqDtva5MdjDGmGbFmx4mqPgIMBqYBaUAVcJmqPhEqk8B6TQBWqupqVa0BngImR5SZDDzivX4GOM2rw2TgKVWtVtU1wErvfLGcM2GOWvt/ANRJK1soXQtBAk20hOJYVbUpxZPdXHQrXkvcOY0xphGxdse9JSLfBbqr6qOq+hNVvRtYICKnisgjwJQE1qsvsCHsfYm3LWoZVa0DdgM9mjg2lnMCICLXichcEZlbVlbWisuAX89YQk1dK3osUwKQU9R4EKpJcItlwLGQXQifv5C4cxpjTCNiDUITgXrgXyKySUQ+F5HVwArgMuA+VX04gfWK1qqKzFFurExLtx+8UXWqqo5X1fEFBQVNVrQ5b6ws5+H/rmnVOVyadmPdceXxz6AdTUoARp/rWkI10RIiW2HNu/DaHbB+dmLPa4zpsGIKQqpapap/VdUvAQOB04AjVXWgqn5TVecnuF4lQP+w9/2AyKZAQxkRSQXygB1NHBvLOROmftw1AGwL9ORXM5Yyc3ErJgdtauqe6vLEjQmFHHY+1FbCylmJO+fKN+CR8+C/f4KHJsKCaYk7tzGmw2rxtD2qWquqpaq6y48KeeYAw0VksIik4xINpkeUmc7+LsCLgDdVVb3tl3rZc4OB4cDHMZ4zYQLpWZCWzSvfP4Eu6QG+9dg8fv5SnDkQuX1gT2n0G1bjXVW1KQOOgy49E3fjajAIr/wEegyFH62CQcfDCzdAmc30ZMyhrl0u5eCN8dwIzASWANNUdbGI3C0i53nFHgR6iMhK4GbgVu/Yxbjkic+BV4EbvCmGop7T72sZVpjDf75zHIN6dOGfH6zhzaVxrIyeUwS1FdHndWvNqqqNCaS6LrnlM6F2X+vPt/5D2L4CTvwRZPeEix6CtC7w4vdsWXFjDnHtMggBqOoMVR2hqkNV9Zfetp+p6nTvdZWqXqyqw1R1gqquDjv2l95xI1X1labO6SsvYXBU71xeuPF4ivvk8s1H53HX9MVUVLfgHt9QmnZ56cH7qn1oCYHLkqutcDMotNbnz7ugM/pc975rIZx+twtOS19q/fmNMR1Wuw1CHV7EX/h5WWk8dd2xfHV8fx75cC3XPz4v9qy5hnuFIsaF6uugbl/rlnFozKATICs/MVlyq9+BgcdBevb+bWMvh54j4Y2fu+swxhySLAj56sCEvK4Zqfz6gjHce8HhvLdiGxP/911unjafDTuayULL8WZMimwJJWoG7WgCqTD6HFj2KtRWxX+e8s2wbRkMPvHg8592h9u34OnW1dUY02FZEPJN42MdXz2qP7+96HD6dsti5qLNnPnHd3n0w7UEg40cEwpCkS2hRKwl1JTiyS7QrXoz/nNsnOeeBxx78L5R50DvMfD+fRCsj/8zjDEdlgUhPzUxicTF4/vz2DVH89rNJzFuYHd+9sJiLntgNuu2VxxcOC3TLfV9UBAKrarq0/Q6g0+CzG6t65Irc9MXUTDy4H0icPzNLmlhyYvxf4YxpsOyIOSXGLO++nbL4tFvTOA3Fx7O55v2MPGP7/HwB2sObhXl9Dm4O66hJZTAm1XDBdJca2XZDKirju8c21a4llxmXvT9xZMhfyi8/wfLlDPmEGRByFexTacnInz1qP68dvOJHD0kn7te/JxLH5jN2m1hraLcooNnTQgFoUyfghC4G1er97ilv+OxbRn0HN74/pQAHP8DKP0MVr0R32cYYzosC0K+aflf9UV5Wfzz6qP4zUWHs2TTHk6/7x3umr6Y7Xur99+wGq7au2/IrzEh8Lrk8uK7cVXV3ZDaM0pXXLjDL3FLVrz3h/jqaIzpsCwI+SmOecVFhK+O788bt5zkxo1mr+PE37zF/F1ZULntwG4xv7vjAFLTYeTZsOxlqKtp2bHlm11iQ7TxoMjPOPZGWPcBrP8o/roaYzocC0J+aeX4RmFuJr/6yhhm3nQi4wbl8+RSdy/NnrKwicD9zo4LKZ7sZmtY807LjtvmJSU01R0XMm6Kuy/p/ftaXj9jTIdlQcg3jU3c3TLDCrvy8NVHccpRYwH48UOvMnftDrezao/7DL+y40KGnuJaWy1dcTU0N1xz3XHgbmQ9+nq0iBdPAAAfRklEQVRY/gps8WWdQWNMO2RByE8JWucvJUWYdOyRAPSWnVwydTZ/eXMFweo9rhWU4vOvMTUDRk6CpS9DfW3sx21b5oJXTu/Yyk/4JqRlwwd/jK+expgOx4KQXxKdbpzrbli99fhczh5TxO9eW867C1dTn9ZGS3AXnw/7dro1gWK1bTn0HBF7MO6SD+O/DgufgZ1r46qmMaZjsSDkqwSueJ7ZDVKzyKzcwv9eOpbfXHg4NRW7WLM3wDvLW7f6a0yGnurmqGvJjatlXhBqiWO+A5IC//1Ly46LV8lcN3/dnAe97k1jTFuyIOSbBLeERFyadvmmhvuKju+fQU0gmykPfcw9ryyltr4Vy4g3Jy0TRk50s17HMuFo1W7YuxkKWhiE8vrCFy+FTx+DvVvjq2usPn0C/vFleO/38PLNcP/RtuqrMW3MgpCfEjQm1CDiXqEuWsnIgX25bMIA/v7OKi75vw8p2ZngJbnDFU+Gyu2w7v3my25b4Z5jSUqI9KWbXCr67L+1/NhYbV8FL93kJla9bQNcMwvSsuDR8xOzfIUxJiYWhPzixxQ0OUUHzh9XtYdAVi6/vmAMf77sCFZs2ctZ//sery5qxVLiTRn2ZZc4EMuNq6E541raHQfQcxgUnwdz/gH7fFrA953fgATggqkuuaP/BLjmNffZT10BJfP8+VxjzAEsCPnKh5ZQealbLhu8VVXdPULnfrEPL3/vBAb1zOb6x+dxx/OLqKpN8MzUaVkw4kw32Whzs15vWwaBdOg+KL7POuGHbrqg//45vuObUrkDFj0L464+MHMvuydc8R+36N6TX4Udqxs9hTEmMdpdEBKRfBGZJSIrvOfujZSb4pVZISJTvG1dRORlEVkqIotF5J6w8leLSJmIzPce1/p7JT60hHL7QLDWdYmpQtUul7DgGdCjC89cfxzXHj+Yx2avY/JfPmD5lvLE1qF4spu5Yd0HTZcrW+4mJg2kxvc5RYfDYRe4LrlEjw0tetb9HMd+7eB9XQtdINIgPPHV6EuqG2MSpt0FIeBW4A1VHQ684b0/gIjkA3cCRwMTgDvDgtXvVHUUcATwJRGZFHbo06o61nv8w9ercBVN7PkaFrfbBLWVUFfllngIk56awk/PKebhrx/F9opqzv3z+zw+ex2aqO7B4We4pbqby5LbtrzlSQmRTv2pu8Z3f9e680Ra+AwUFru1jKLpOQwueQx2roFnv2lrHRnjo/YYhCYDj3ivHwHOj1LmTGCWqu5Q1Z3ALGCiqlaq6lsAqloDfAL0a4M6H8yPMaHwZb4rt7vXEUEo5OSRhbzy/RM5ekgPfvr8Ir712Dx2VbZw7rdo0rvA8NOb7pKrq3Zf4PGMB4XrMRSOuALmPgQ717XuXCH7dkLJx26Jiqb+SBh0PEy6F1bMhDd/kZjPNsYcpD0GoV6qWgrgPRdGKdMXCJtEjRJvWwMR6Qaci2tNhVwoIgtE5BkR6d9YBUTkOhGZKyJzy8pacw+OD2NC4AUhb+qeLvmNFi/IyeDhq4/i9rNG89ayrUz63/eYvXp76+tRfD7s3QLrZ7O1vIpXF5Xyh1nLuf25hfzt7VVsXrPYdWfFkxkX6aSfuPuG3v51688FsOY9V7ehpzRfdvw1btzo/T+4LjxjTMLF2WHfOiLyOhBtLpfbYz1FlG0NTQ8RSQX+BfxJVUOjyy8C/1LVahG5HtfKOjXayVV1KjAVYPz48XE2aXxoCWUXui/k8tJmW0IhKSnCN08cwjFDevC9pz7lsgdmc/7YvlxxzEDG9u9GIKX5QFkfVKbN3cCMhaVs3LUPqU7nZdJ57uE/c1vVle5zBPKy0thZWcvy1Nnclwq1PUeS1tprzusLR1/nbl495jturKg1Vr/l5trrd1TzZUVg0m9dpt/zN7gxrj5jW/f5xpgDJCUIqeqXG9snIltEpEhVS0WkCIg2Kl0CnBz2vh/wdtj7qcAKVW2YhExVw5sADwD3xlH12KkmfkwokApde0W0hJoOQiFj+uXx0neP54+vL+eJj9bz3Kcb6dYljaMG5TNhUD4TBuczpm8eKRFBaWHJbv7nuYUs3Lib4YVdGV2US1Zad1ZuOJazK+dQPunnjB/ck+KiXDLTApTsrGTJE69QV5bCFc/v4L6v7aNPt6zWXfcJP3Q3lr52O1w1vXU/11Vvua62QIzhMTUdvvooTD0FnrocrnsbuhbE//nGmAMkJQg1YzowBbjHe442Aj4T+FVYMsIZwG0AIvILIA84IPstFNi8t+cBSxJf9UgJDkLg3bC6CSq8bsIYgxBAdkYqt59dzI2nDued5WW8t7yMOWt3MOvzLQD0zs3knMOLOPvwItICKTz58Xqe+ng9Pbpm8JevHcHZY4qQUABYeBU8ew3XDd4OA/Yv1dCvexf69dhOefVgFm2u4qw/vcdvL/oipxf3iv+as7rBKf8DM34Iy16BUWfFd56da91Y1dHXt+y4roVw6RPw0Jnwn2td9lxKIL46xCoYdMkdGnRja/FmGRrTzrXHf9n3ANNE5BpgPXAxgIiMB65X1WtVdYeI/ByY4x1zt7etH65LbynwifeF+RcvE+57InIeUAfsAK729zJ86I4Dd99NyRy31HdaF8iKmsHepLysNM77Yh/O+6IbY9q6p4r/rtrOSwtKeeTDtfzj/TUApKYIU44bxA9OH0FuZkTLYcSZEMhwN64OOObAfVs/J2fAEbx0ygl891+f8M1H53L2mCJuPmMEQwvinHB13NXw8VR47afuptnU9JafY9Vb7jmW8aBIfcbCpN/Ai99z2Xon/6Tl54jV8tfg5Vtg93r3PivfdUUee4NLDDGmE2l3QcjrNjstyva5hLVuVPUh4KGIMiU00vxQ1dvwWkttJtHdcQCFo90gedlStyR2Aj6jMDeT84/oy/lH9GV3ZS3vrigjqMqxQ3pQmJsZ/aCMHBcMlkyHM3+1fzmJqj2uxTH2cgb3zObZbx/HX99axQPvrWbGolK+PLoX3/jSYI4Zkr+/VRWLQBqc8Ut48mKY+yB1R32LNdsqKNtbDeqSMIYUdG16jGv1W5DTJ/6svSOvgnX/dUkSA46GISfHd56mLHwG/vNNKBgFk/8KKaluHae3fgELnnIzPPQdl/jPNSZJ2l0Q6jT8SNEGd38LuPnNhkbNq2iVvC5pnOu1kJp12Plu2e+Nc920N+BaaQD9xgOQkRrgB6eP4IpjBvLoh2t5fPY6Zn2+heKiXK45fjDnfLGIjNTmu7aq6+pZmDaOwu5Hk//aLzljRgGbag4ca8pOD3DiiAIuHt+Pk0cUHji+Fax3y1CMPCv+wC0C5/wBNn0Kz14L178f+1pJsdi8CF64AQYcC5f/2y30B/DFS2D1O/D8d+Afp8Pp/88th+7HHznGtDELQr7yoyVUvP91ry8k/vwtMeJMNzXP4uf3B6ENH7sMvr7jDyhakJPBLWeM5IZThvHcpxt56P013PLvz7jn1aVcecxALjiyL/267+9q2llRw6JNu5mzdicfr9nOp+t3UV0XZKR8hRkZc/hd4atsPu4ueudlIgilu/cxb91OXl20mVcWbWZU7xxuOWMkXx5d6FpcpfPdPUJD4uiKC5ee7RIVHjgF/v11mDI99iSHpgSDMP27bhHAix/eH4BChpwE3/4Apt/ouiTXz4bJ97vxMmM6MEnYnfSd1Pjx43Xu3LktP/CFG2Dlm3BLgvMfVOH3I919Opc8DqPPTez5W+rpK2HNO3DTIsjMdX+p19fAt95p8jBV5b0V23jw/TUN6yH1ys0gLyuNXZW1bC2vBlzq92F98pgw2GXwHTUon/w3fwSfPg7f/vCgWRlq6oK8vHATf3x9Beu2V1JclMu3ThrCpJ2Pk/7Or+BHq9wcca21YJrrNjvuu3DGgTez7q2uo6y8ml2VNaSIkJ6aQq/cTLp3SWu8C/KTR10Q+spU1/JpjKqbymjWHZDXDy5+pG3TxutqXEAvL3WBsmCUq4cxEURknqqOb7acBaGmtSoIrXoLbv488ZVa+z4snwmn3Zn8rKlNn8LUk+GUn7p1gP44xmWynfTjmE+xdlsFry/ZwpLScvZW15Kbmcawwq4c1iePL/bPIycyKWJvGfxlHBSNhateiNotVVsf5LlPN/L3d1axuqyCael3kxuo4cacP5KZlkJmaoC+3bO48Mh+nDC8Z8vGp0JevgXm/IPqCx7mdY7htc838+n6XazfEX05jcy0FAb37EpxUS7FfXI5rE8uX+ibR9fgXvjzOOgxDL7xamzdbBs+hn9fDRXbYNI9MO7r/nbP1VTAe3+Ajx+A6oj59HqOgDFfhSMu339DtTnkWRBKkLiD0PM3uIFwP4JQe/PU5bBilvuLfOM8+N586NbohBSJ8fEDLmX7wgdhzEWNFgsGlY+XreOoaUfydo9Lebb7NVTVBqmqrWdJ6R52VtZy3NAe3HFOMaOLcmP++GBQ+XhlKUXPXUjPytVcXvM/bOhSzNFD8jmsTx59umXSLSudoCpVtUG27Kli0659rNi6l89L91AW1tL7fc5TTK55kZnHPc3ALxzLiF5dSQ3EMJlJxXbXGlv1hgsC59wHGT4s977lc5h2FWxf6ZbY+MJFkD/YzeK+6VOXNr/2Pbc0xmFfgeNvanxePnPIsCCUIK0LQm/DzYsTXqd2Z28ZPDoZti52rbMTbvb/M4P18MCpUL4ZbpzjugIb8/kL7kt0yksw+ISGzTV1QZ6as577Zi1n975avnb0AL536vADMgIra+rYtKuK0t37KN1dxZbdVazZVsEHq7axZU81QzP2MC3jbvK0HLl8GoFBx8VU/a3lVSzetIcNyz7h8k8u4zlO44dVXwcgKy3AmL6uFThuYD7jB3WnZ9eMRn4OQbcy7Nu/gh7D3XhV4ShUlW17a9iws5Ide2uorgtSXVdPIEXISgvQJT2VnMxUeudl0rNrRuNZhRvnwWMXQGqmy8wbclL0cttXuTn+5j0MNXtd5uTxP4CBX2rbBIqaCpedWbsPMvPcpL9+BGbTLAtCCRJ/EPqOy2g6FIIQuKBQtbvJuewSrmQe/OM0d/PppHsaL/fvr7vMuFuWRe2+3FVZwx9fX8Fjs9cRVGVAfhdSU4TtFTXsqqw9qHzPrhkcPSSfL48uZOJhRWRVbnRBeOc6OOU2OO57kNpI0AinCo9fACXz0O/OY311F+Zv2NXwWLxpDzV1bu2oIT2zGV2US7/uWRTlZdIlI5XMNJdVWFldR07pfzlp4a2k1e/j3qyb+Ff5EeyLcT2pQIpQmJPBkIJshhfmMLxXV0b0ymF07WK6/vsy9zudMj22taH27YQ5D7pxq8ptLp382Bth9Hn+dR1X7oD5T8LCabB5obvBN1yP4a4eQ06GYae5m4+N7ywIJUirgtCad+EHixJfKbPfSz+Auf+EKS8e0MppUFMBvx0Gh18C5/7x4P1hVpft5eUFpazYupe6YJAe2Rn0zsukT7dMivKy6JOXRWFuRsOX/wGqdsP077l7evKHwGk/g9GT998/Fc38f8Hz18PEe+CYbx+0u7qunkUbXYbg3LU7WVW2l40791FTH4xyMhiUvpu/pv0vxfVLeaPP9Wwo/hb9e2RTmJNJRloK6YEU6oJKVW09+2rr2VVZy+Y9rnW3adc+VpXtZcXWvVTW1HO0LOGh9N+wLaUHf+r7ewr7DWFU7xxG9c5lSEE2ac11F9bug/lPwIf3u8UBuw2AY26AI64gmJbN5j1VrC6rYPW2vd5zBWXl1Q1LjnTNSKVbl3Tys9Po260LA3uEHtn7EzxqKmD2X+GDP7kFEPuOg2GnQ8FIlzSxbxfsWue6DEvm7J9lpGisS9U/7HxX1k+VO2DVm1D6mZuDsGoX1Ne6VlrXQug5HApGu67sTpbgYUEoQeIOQs992/WTWxDyV/VemHoS1FS6FObIltjch1yg+sbMg2d28MPK12Hm7e5m4sLD3MwKo849OBhtW+G6E3t9Aa5+KeZpgIJBZWdlDftq6xtWzu2Snkp2Riq5malIXbXLsls4zQXes37rvvBiFAwq2xe9Tv4LV7I7ozf39fkdc7als3LrXuqC7rsiLSAMLXDzCI7sncOo3jkMK+xKfnY6WWkBRIS6+qDrDtxeTt2SGQxY+iB9yz+jXLJ5rv5Enqg9iWU6AHD3dw0p6Eqv3ExCvYJ7q+vYUVHDjoqahkzJkJyMAFNyPubafY/QrX4bJb1OYdtRP6Rw2Dh652YeNP+hd2GweQGsnOXGLzd8DChaMJqakeeyY+BZlGUNZofX+t1ZWUNFdR1pgRQy0wJ0zUglLyuNbl3SyMtyj9ystIP/IFGFLYth+auw4jUX/DTobmXoMdxlZqakuj9a9m6B3WGLAeQNgIHHwsDjYNAJ7o+ZturKDNa7QK7qWvGpma2emsqCUILEH4Suh7UfwA8WJr5S5kCbPnWp4QOOgSue3d8VVlcDfzvWTW/0rXfb9j/04ufg7Xtg+woXaE76iVvDKCXF/VX81OWutfDNN6H7wMR+viq8+1t461fur+0J18Gos93SGk21zFRdt9ZLN7kZw6e82DBZa01dkNXb9rK0tJylm8tZunkPyzaXU7q76oBTpAdSUJTa+gO/VwIpwpm567k6ZQbj9n1AQOvY22MMwTGXkHP4OUj+4EarVVVbz4YdlazbVkHdyjcoXv5XBlQsYmnKMO6suYKP6van6aenptC/exa98zJdcE4PkBpIoa4+SG29UlMfZPe+WlL2ljKu4n1OqvuA8bKMFFFWBPvySvAoZtWPZ6EOJpb7/DJSU+iRnc6orJ1ckPpfjq98k26Vbtorisa6e+mGnwFFX4x+P1l1uWshlcyF9f+FdR9ChTdnc04fGHQ89QO/xK5ex7AlpYjtlS4wV9e6Mb7quiCqkJkeIDt9f7DMCwuWmakB6oJKUJX6oFJXX09g+3KyN89BNsyG9R/Crg0cONWYuCnBrvwP9Dmi2Z9DNBaEEqRVQWjdB3CTBaE28dnT8Nx1MPJsuPAB1x3z9j1uip2vTXNfBm0tWO+mWHrnXpdZlj/U/cfeOM/NtPC1p92Xk182zoNZd7oWOUBatvvruscQV5ceQ91zXl+3HPtHf3MtucEnuRtmYxjf21VZw7LN5azeVsGuylr3BS+QmRage3Y6/btn0T+/C/27dyE91QuAFdtdS+3Tx2GL11NQWOxmiij6oqtjdoH70q6tdIkGJXNh6Uvu55hTBKfcDmMvp06hdHcV67ZXsm5HBeu3V7J2ewXb9tZQWVNPZU0dtXVB0lJTSAu4R25mKt27pNM9O41uXdLpE9jFYbvfZdCWWXTfNhfRIMGuvdERk6gdegYVhUeyR3LYva/2gEfN7q302TSTEWUzGVq5AICPgqN4of5LvB48kq49+zGiMIeeOen0yM6ga0YqmWkpZKQFyEwLkJUWIDMtBUHYW11HeVUtOyuqqd6ynO5bP2Lgnk8orllAD3a5X6f2YHZwNLODxcwLjmCt9iYYw5JwWVQxRtZwZMoKxqcsY3zKcrpJBQDb6M6y9MPY0WUwKV26EQikQd0+AnWVZNTsZOgFd9F34NAW/bMLsSCUIHEHof98y/1lY0Go7Xw0FV79iZtTr8dQl514+CUuqyuZ6utg0TMuINVVQf9jXDJFduwzoLfKznXuD6LSz1wW247VbqwkWHdgucw8OPHHrm5tdf/Z9lUuxXvFa65FW70nejkJuCU4Dv8qjLk4tsSPeFRsd3VZ9rK72bzWfVmTP9QlZnTJd7/D7auhbMn+Wc69em1N7c2ijbtZtHEPCzfuZnXZXtfFt6825pm8cjJc1mJRtyyKcjIYnb6Z0VXzGbDnE3pun0NalVuVRlMzCfYYTrDnKGqyi6hO706VplG9r4L6iu2k7VlH9p7V5FesJgXXdbuzyyC2dDuCzXljWZJ+GCtqelK2t4ate6rZWl5FdV2QzLQAmamuG/IvXzuS4j6x37oQzoJQgrQuCH0INy1IfKVM49a851KWy0vdbBIn/ji+Gbc7u/pa2LXeBaTdJZDX333JpzUyYW1bCAbdUhu7S1wSQbDeBZtu/V1XYlunWtdWwYaP3NyImz6F3RvdYpJpWe7n1W+862LtdVizXb31QWVfbT37atxYXnVdPVW1QfbV1qPqEjFyMlPp1iXt4Juzw6nC1iWuPls/d6/LlrrxpfA/KlJSXTJI/lDXndbvKJe40VZ/+BB7ELK543xjwT0pBp8QPUvOHCiQ5lqLPeLravFFSkr7qlNaprsvqrF7o1ogkCJ0zUila0Yrv3JFoFexe4RTdckO9TUuqSA92/81rxLEgpCfbJZjY0xbEOmwk9nGMDeIiYt1cxpjTLPaXRASkXwRmSUiK7znqEuHisgUr8wKEZkStv1tEVkmIvO9R6G3PUNEnhaRlSLykYgMaoOr8f8jjDGmA2t3QQi4FXhDVYcDb3jvDyAi+cCdwNHABODOiGB1uaqO9R5e0j3XADtVdRhwH3CvnxdhY0LGGNO89hiEJgOPeK8fAc6PUuZMYJaq7lDVncAsYGILzvsMcJrENX9/C9iYkDHGNKk9BqFeqloK4D1Hm22wLxA23wUl3raQf3pdcXeEBZqGY1S1DtgN+JevaGNCxhjTrKRkx4nI60DvKLtuj/UUUbaFvvUvV9WNIpIDPAtcCTzazDGR9bsOuA5gwIABMVYp1moaY4wJSUoQUtUvN7ZPRLaISJGqlopIEbA1SrES4OSw9/2At71zb/Sey0XkSdyY0aPeMf2BEhFJBfKAHY3UbyowFdzNqi26uP1nie8wY4w5hLTH7rjpQCjbbQrwQpQyM4EzRKS7l5BwBjBTRFJFpCeAiKQB5wChaazDz3sR8Kb6OV2Eqo0JGWNMM9rjzar3ANNE5BpgPXAxgIiMB65X1WtVdYeI/ByY4x1zt7ctGxeM0oAA8DrwgFfmQeAxEVmJawFd6v+lWBAyxpimtLsgpKrbgdOibJ8LXBv2/iHgoYgyFcC4Rs5bhRfQ2oZ1xxljTHPaY3dc52HdccYY06R21xLqNPodBV2jJQAaY4wJsSDkl2NvSHYNjDGm3bPuOGOMMUljQcgYY0zSWBAyxhiTNBaEjDHGJI0FIWOMMUljQcgYY0zSWBAyxhiTNBaEjDHGJI34OZF0ZyAiZcC6OA/vCWxLYHU6ArvmQ4Ndc+fX2usdqKoFzRWyIOQjEZmrquOTXY+2ZNd8aLBr7vza6nqtO84YY0zSWBAyxhiTNBaE/DU12RVIArvmQ4Ndc+fXJtdrY0LGGGOSxlpCxhhjksaCkA9EZKKILBORlSJya7Lrk0gislZEForIfBGZ623LF5FZIrLCe+7ubRcR+ZP3c1ggIkcmt/axEZGHRGSriCwK29biaxSRKV75FSIyJRnXEqtGrvkuEdno/a7ni8hZYftu8655mYicGba9w/zbF5H+IvKWiCwRkcUi8n1ve6f9XTdxzcn7XauqPRL4AALAKmAIkA58BhQnu14JvL61QM+Ibb8BbvVe3wrc670+C3gFEOAY4KNk1z/GazwROBJYFO81AvnAau+5u/e6e7KvrYXXfBfwwyhli71/1xnAYO/fe6Cj/dsHioAjvdc5wHLv2jrt77qJa07a79paQok3AVipqqtVtQZ4Cpic5Dr5bTLwiPf6EeD8sO2PqjMb6CYiRcmoYEuo6rvAjojNLb3GM4FZqrpDVXcCs4CJ/tc+Po1cc2MmA0+parWqrgFW4v7dd6h/+6paqqqfeK/LgSVAXzrx77qJa26M779rC0KJ1xfYEPa+hKZ/yR2NAq+JyDwRuc7b1ktVS8H9IwcKve2d6WfR0mvsLNd+o9f19FCoW4pOeM0iMgg4AviIQ+R3HXHNkKTftQWhxJMo2zpTCuKXVPVIYBJwg4ic2ETZzv6zgMavsTNc+9+AocBYoBT4vbe9U12ziHQFngVuUtU9TRWNsq1DXneUa07a79qCUOKVAP3D3vcDNiWpLgmnqpu8563Ac7hm+ZZQN5v3vNUr3pl+Fi29xg5/7aq6RVXrVTUIPID7XUMnumYRScN9GT+hqv/xNnfq33W0a07m79qCUOLNAYaLyGARSQcuBaYnuU4JISLZIpITeg2cASzCXV8oI2gK8IL3ejpwlZdVdAywO9TN0QG19BpnAmeISHeva+MMb1uHETF+9xXc7xrcNV8qIhkiMhgYDnxMB/u3LyICPAgsUdU/hO3qtL/rxq45qb/rZGdrdMYHLotmOS575PZk1yeB1zUElwXzGbA4dG1AD+ANYIX3nO9tF+B+7+ewEBif7GuI8Tr/heuSqMX9xXdNPNcIfAM3kLsS+HqyryuOa37Mu6YF3hdMUVj5271rXgZMCtveYf7tA8fjupAWAPO9x1md+XfdxDUn7XdtMyYYY4xJGuuOM8YYkzQWhIwxxiSNBSFjjDFJY0HIGGNM0lgQMsYYkzQWhIwxxiSNBSFjjDFJY0HImA7Im73iL96d+8Z0WBaEjOmYrset8XJ8sitiTGtYEDKmY5qImzJlfrIrYkxrWBAypoMRkUzcypZHAu8kuTrGtIoFIWM6nuG4ILRUVWuTXRljWiM12RUwxrRYATCCdrx0tjGxspaQMR1PH9yiZClhyzAb0yFZEDKmAxGRVNxYUG/g70B9cmtkTOvYekLGGGOSxlpCxhhjksaCkDHGmKSxIGSMMSZpLAgZY4xJGgtCxhhjksaCkDHGmKSxIGSMMSZpLAgZY4xJmv8PVo76Y48gFyEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] @@ -83,8 +110,8 @@ } ], "source": [ - "cls_new = LambdaCDM_new.lensed_cl(2500)\n", - "cls_std = LambdaCDM_std.lensed_cl(2500)\n", + "cls_new = hyrec_mod.lensed_cl(2500)\n", + "cls_std = LambdaCDM.lensed_cl(2500)\n", "\n", "ll_vec_new = cls_new['ell']\n", "ll_vec_std = cls_std['ell']\n", @@ -105,56 +132,45 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.max(np.abs(\n", + " (hyrec_mod.get_thermodynamics()['x_e'] - LambdaCDM.get_thermodynamics()['x_e'])\n", + " /LambdaCDM.get_thermodynamics()['x_e']\n", + "))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%autoreload\n", + "# LambdaCDM_hyrec = Class()\n", + "# LambdaCDM_hyrec.set({'recombination': 'hyrec', 'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n", "LambdaCDM_newAlpha = Class()\n", - "LambdaCDM_newAlpha.set({'recombination': 'hyrec', 'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n", + "LambdaCDM_newAlpha.set({\n", + " 'recombination': 'hyrec', \n", + " 'hyrec_Alpha_inf_file': './hyrec/Alpha_inf.dat',\n", + " 'output':'tCl,pCl,lCl,mPk',\n", + " 'lensing':'yes',\n", + " 'P_k_max_1/Mpc':3.0\n", + "})\n", + "# LambdaCDM_hyrec.compute()\n", "LambdaCDM_newAlpha.compute()" ] }, { "cell_type": "code", - "execution_count": 23, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/hongwan/anaconda/envs/python3/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in true_divide\n", - " # This is added back by InteractiveShellApp.init_path()\n", - "/Users/hongwan/anaconda/envs/python3/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: invalid value encountered in true_divide\n", - " if sys.path[0] == '':\n" - ] - }, - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 23, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VFX6+PHPSQ9JSAiBJARCgASQXiIdRenoCnZQEbCgrnUtX3Vdd112/a266qqLq6KCNEFARKQoFqSIIAktdEIoCUlISAjpbeb8/riDRkifmUwy87xfr7xm5s7JPc9lwn3mnnuK0lojhBDC9bg5OgAhhBCOIQlACCFclCQAIYRwUZIAhBDCRUkCEEIIFyUJQAghXJQkACGEcFGSAIQQwkVJAhBCCBfl4egAqhMSEqKjoqIcHYYQQjQZ8fHx57TWrWpTtlEngKioKOLi4hwdhhBCNBlKqVO1LStNQEII4aIkAQghhIuSBCCEEC7K6gSglGqnlNqolDqklDqglHq8kjIjlFIXlFJ7LD9/tbZeIYQQ1rHFTeBy4Cmt9S6lVAAQr5T6Vmt98JJyW7TW19ugPiGEEDZg9RWA1jpNa73L8jwPOAREWLtfIYQQ9mXTewBKqSigL7CjkrcHK6X2KqXWK6W627JeIYQQdWezBKCU8gc+B57QWude8vYuoL3WujfwX2BVNfuZqZSKU0rFZWZm2io8IYRoGk79DD+/Cw2wXK9NEoBSyhPj5L9Ya73y0ve11rla63zL83WAp1IqpLJ9aa3naK1jtdaxrVrVajCbEEI4h8JsWHEP7PwISgvsXp0tegEp4GPgkNb6zSrKhFnKoZQaYKk3y9q6hRDCqax/Fgoy4Ja54O1v9+ps0QtoKDAVSFBK7bFs+zMQCaC1fh+4BXhIKVUOFAGTtW6A6xshhGgqkndCwjK46hlo07dBqrQ6AWittwKqhjKzgdnW1iWEEE5Ja9jwF/BrDUOfaLBqZSSwEEI4WtJGSN4OI55rkKafiyQBCCGEo237L/iHQt+7GrRaSQBCCOFI6fvh+A8wYCZ4eDdo1ZIAhBDCkX75ADybQew9DV51o14QRgghmprScjMJZ3I4djaf/JJylFJEBPkQ3TqATq38sPSItxQuhP1fQPcboVlwg8cqCUAIIWwgMSOfOZuPsz4hnbyS8krLhPh7c3XnVkzs04YhnVricXgtlOZB78kNHK1BEoAQQlghr7iMf60/zJJfTuPt4cb1vdowulso3cKbE9TMk3KT5kxOEfvPXODnpCw2HEzn810phPh7scT3QyL9IvCMHOqQ9nhJAEIIUU97knP446J40nOLmT4kioeviSbE//IbuS38vOgREcjkAZGUlJv48Ugmm+L20THpF941TWTpaz8ytkcY/du3IKZ1AMF+XrQKsP8NYUkAQghRD6t2n+H/Pt9HaHNvPn9oCH0jW9Tq97w93BnbPYyxOZ/BCc0VY2fS9bgXi3ecZt5PJwEI9vNi14uj7Ri9QRKAEELU0dytJ5i15iCDOgbzvzv7E+znVbcdaA17lkDbAYwePpTRw6HMZOZwWh4nswooN5vtE/glJAEIIUQdfLz1BP9Yc5Bx3cN4Z0pfvDzq0XqftgcyD8F1v82f6enuRs+2gfRsG2jDaKsnCUAIIWrpoy1J/HPtIcb3ME7+nu71vHW7dym4e0GPm2wbYB3JQDAhhKiFiyf/CT2tPPmXl0LCcugyHnxrd9/AXuQKQAghavDh5iReXneI63qG89bkPvU/+QMkfgeFWdD7DtsFWE+SAIQQruvMLtj+Pzj5E5TkQkA4dLoWrrwPWnUG4H8/JvLa10e4rlc4b91u5ckfYO+n4NcKokfa4ACsIwlACNHklZvMZBeWYjJrAn09aeZVw6lNa/jhn7DlDfBpDjFjwS8Eso5D/Cfwyxx0/xm86zGV1zelMbFPG964tTce1p78C7LgyNcw4H5w97RuXzZgdQJQSrUDFgBhgBmYo7V++5IyCngbmAAUAtO11rusrVsI4bqSswtZHpfMD0cyOJqeT6npt66TLf286NMuiIEdgxnTLYyoEL/fflFrWPsUxH1sTL889l9GErio4Bzlm/6N2y9zuM68luLur/Kn2/rg7lbtule1s38FmMugz53W78sGbHEFUA48pbXepZQKAOKVUt9qrQ9WKDMeiLH8DATeszwKIUSdnM0t5pX1h/lyzxkAYtsHM2NYFBFBvni6u5FdUMrJcwXEnz7P94cz+H/rDtMjojnX9WzD+B5hRB1fbJz8hzwGo2eB+v2JPf6cG88cGk+LklAWBMzmqdMPo062gI4jrA9+9yII7w1hPazflw3YYknINCDN8jxPKXUIiAAqJoCJwALLOsDblVJBSqlwy+8KIUStrIhP4e+rD1BiMnPvsA7cO6wjYYE+VZZPzSliXUIaa/al8erXh/nsm41s8P4zZ4KHkdbxUaLzSvD1cie7oJQ9yTms3HWGTUczaRPowxP3TMUvdAosugkW3wa3L4TOY+sffHoCpO+D8f+u/z5szKb3AJRSUUBfYMclb0UAyRVep1i2SQIQQtSozGRm1lcHWbj9FAM7BPPqzb1+36xThTZBvtw3vCP3De9ISnYBauEkynK8uD1tChkf7bysfKsAb54Z24VpQ6Lw97acHqevhYWT4LO74M7l9b8S2LXQ6Pvf85b6/b4d2CwBKKX8gc+BJ7TWuZe+Xcmv6Cr2MxOYCRAZGWmr8IQQTVRJuYmHF+/mu0NneeCqjjwztku9bsa2zdoG53+BCa+zvtutHEjN5WRWAcVlJoJ8vegcFkCviEDcLm3rbxYMd38J8ybA0jth2lcQ0a9ulRdfgD2LoftNDpn3vyo2SQBKKU+Mk/9irfXKSoqkAO0qvG4LpFa2L631HGAOQGxsbKVJQgjhGorLTDywMJ5NRzP5x8TuTB0cVb8daQ0bX4ag9tB/Oi3dPbmqcyuuolXtft+3Bdy1EuaOgcW3wD3fQEhM7evftRBK82HQQ/WL306sHgls6eHzMXBIa/1mFcVWA3crwyDggrT/CyGqYzJrHl+6m83HMnn15p71P/kDHPsWUnfDVc/Uv/tl83CYugpQsGAS5CTX+CuAMfJ3xwfQfii06VO/uu3EFlNBDAWmAtcqpfZYfiYopR5USj1oKbMOSAISgQ+BP9qgXiGEE3t57SG+OXCWF6/rxu1XWtkcvON9Y5CXtStvtewEU1cag8YWToL8zJp/J34eXDgNw5+0rm47sEUvoK1U3sZfsYwGHra2LiGEa5j30wnm/nSCGUOjuGdYB+t2di4Rjn8PI/5sm8FX4b3hjmWw8EZYdCNMWwO+QZWXLcyGTa9B1HDo5PiRv5eSyeCEEI3KhgPpzFpzkLHdQ/nLdd2s32Hcx+DmCf2nW7+vi9oPhsmLIOOw0U204NzlZbSGtU8aN4DH/euy8QaNgSQAIUSjsSc5h8eW7qZ32yDeur2v9aNvy0th7xK44noICLVNkBdFj4LbFsDZA/DRSEiJ++09sxm+nwUHvoBrnoewnrat20ZkLiAhRKNwKquAez/ZSasAbz6aFouvl7v1Oz3+PRSdh15Wtv1XpesEY5zAZ1Pho1EQMxpadYVT2+BMHPSbBsMaX9v/RZIAhBAOl11QyvR5OzFpzSczBlS6sHq9JCwH32D7zrzZNhYe3gHb3oF9y+D4DxDcCW74L/Sd2iibfi6SBCCEcKiCknLum7+T1JwiPr1/IJ1a+dtmxyX5cHgd9Jli/5k3fZrDtX8xfrRu1Cf9iuQegBDCYfJLypk+7xf2plzg7cl96N/ehqNkD6+F8iLoeZvt9lkbTeTkD3IFIISoQX5JOZuOZLLjRBYnzhWQlV+Ku5siqJknHUL86NEmkKExIUQE+dZpv+kXipm5MI4Dqbm8M7kv43qE2zbwhOUQ2A7aycTDVZEEIISoVPqFYmZvPMYXu85QUGrCz8ud6Nb+tAnywWTWZBeU8sWuMyz4+RQAHUP8GNGlNdd2bc2VHVrg7VH5TVyttTHA68v9FJaU8/5d/RndzcY9dPIzjbb4oY+BmzR0VEUSgBDid8pNZt7fdJzZGxMxm+GGPm24LbYdfSODLlsOUWvN0bP5bDmWyeZj51i04xRzfzqBn5c7Q6NDuKZra3pGBBLo68mFojL2puTweXwKu07n0DUsgMX3DaRzaIDtD+LgKtAm6Hmr7fftRCQBCCF+dSaniMeX7Cbu1Hkm9Azj+fFX0C64WZXllVJ0CQugS1gA9w3vSGFpOdsSs9h4JIONhzPYcPDsZb/TIcSPWRO7c8eASOuXWKxKwnJo3R1Cu9tn/05CEoAQAoAT5wq448Pt5BWX89btfZjUN6LO+2jm5cGobqGM6haK1ppjGfkkZRaQW1xGgLcHXcObE9WyGcqeN0qzT0DyDhj5N/vV4SQkAQghSMzIY8qHOzCZNcseGEy3Ns1r/qUaKKXoHBpgnyae6iSsMB4b0cIrjZUkACFc3KG0XO76aAdKKZbOHNTwJ2xb0hoSlkHkEAiSBaVqIrfHhXBh+89cYMqH2/F0d2PZA0385A+QthfOHYVecvO3NiQBCOGi9iTncMeH2/Hz8uCzBwbR0VYjcB0pYbkx82e3SY6OpEmQJiAhXFDcyWymz9tJsJ8Xn94/kLYtqu7pYxPpCXD0G8hKBG2GwLbQdgB0ugY8bDTvT3mpMRdPzJhGte5uY2arNYHnAtcDGVrrHpW8PwL4Ejhh2bRSaz3LFnULIerm5+NZ3Dt/J2HNfVh8/0DCA+s2grdOso7D2qcgaaPxunkEuLnDhTNGP33vQIidAYMfAf9ars9blcNroCDDtvP+OzlbXQF8AswGFlRTZovW+nob1SeEqIfvDp7l4U93ERncjMX3D6R1gI/9Kju+EZZNM9YLHPNP6H0H+LU03isrhpNbYc8i+Olt2PmRMZHagJlGgqiPuLkQGGnfmT+djE3uAWitNwPZttiXEML2tNYs+eU0DyyKp2tYAJ89MNi+J//TO2DJFKOp58GtMOTR307+AJ4+EDMKbv0EHtkJkYPh6+eMOfWzjte9vrMH4eQW6D+t/gnEBTXkTeDBSqm9Sqn1SikZnidEA7lQVMbTy/fx/MoEhkaHsPj+QQT7edmvwvxM+OwuaB4Od39Zc3fMkBi4czncMheyk+CDq2H/yrrVueV18PKH2HvqH7cLaqibwLuA9lrrfKXUBGAVEFNZQaXUTGAmQGSk9OMVQmvN+cIyCkvLKS0308zLg+a+Hvh6ulc7ovZCYRnL45P534/HOV9YyuMjY3hsZIz1yyxWHyx89TgU58DUL2rfrq8U9LgZ2l4JK+6BFTOMb/Rj/2VcLVTn7EEjYQx7Qm7+1lGDJACtdW6F5+uUUv9TSoVorS9bSVlrPQeYAxAbG6sbIj4hGpPScjNbEzPZeiyLnSezScrMp6DUdFk5L3c3Wvh5EuznTbDlMdDXg9yick5lFbA/NReTWTOwQzAvXt+NHhGB9g/+2AY4shZGz4Kwy/qD1CwoEmasN9bT3fYOJO+E2+ZDy06VlzebjYXXfYNg8KPWxe6CGiQBKKXCgLNaa62UGoDR9JTVEHUL0VScySnioy1JrNp9hvOFZXh7uNEvsgW3xrajXXAzArw98PRQFJWauVBURk5RKecLSskuKCO7oIT9ORfIKSwlwMeT8EAfHrq6E+N6hDXMiR/AVAYb/gIto2HQH+u/H3dPGPMPiBoGXzwAH1wFf3i78qkdtrwOp3+Gie/+/h6DqBVbdQNdAowAQpRSKcDfAE8ArfX7wC3AQ0qpcqAImKy1lm/3QgDn8kt4Y8NRVsQnAzCmexg394tgaHRIlXPqN0oJK4xRuLcvts0SjJ3HGjeQV9wDn99r9PG/+lmI6Ael+bD530YPol6Toc+d1tfnglRjPg/HxsbquLg4R4chhF2YzZolO0/z2tdHKCwtZ/KVkTw4olOdV9ZqFLSG94YYzx/aZttlEU1l8PNs2PIfKLkAnn5QXmyMI+h3N0x4AzzseFO7iVFKxWutY2tTVkYCC+EA5/JLeGrZXjYdzWRQx2D+OakH0a2b8Dw8id9BxkG48QPbr4nr7gnD/mQM8Dq8Ds7uB89m0GUCtO1v27pcjCQAIRrY9qQsHluym5yiMv4xsTt3DWpv3/nxG8LP7xqjfHvcbL86fFtAX2nqsSVJAEI0oJW7Unj28320a9GMT2YMsMm8+w53/qQx1cOIP9um7V80GEkAQjQArTWzf0jkjW+PMrhjS96f2p9AXyc5We5eDCj5dt4ESQIQws7KTGZeXLWfpTuTubFvBK/e3AsvDyeZid1sgj2LIXqUMe2DaFIkAQhhR/kl5fxx8S42H83kkWuieWpM56bf3l/R8R8g9wyM+5ejIxH1IAlACDtJv1DMPZ/s5MjZPF65qSeTBzjh1Ca75kOzEOg83tGRiHqQBCCEHRxKy+WeT3aSW1TGR9NiuaZLa0eHZHuF2XDka2MKZ+mH3yRJAhDCxjYcSOfJZXvx9/Zg+YNDnKOnT2USVoC5DPrc4ehIRD1JAhDCRorLTLz69WHm/XSSnhGBzLm7v31X23K0PYshrFf9Jn0TjYIkACGspLXmh8MZ/P2rg5zOLmTG0CieG9+1ac3jU1dnD0LaHhj3iqMjEVaQBCBEPWXkFrPh4FkWbT/F4fQ8OrXyY/F9AxkaHeLo0Oxv76fg5gE9b3V0JMIKkgCEqEFRqYnjmfkcz8zn2Nl8jmXkkZiRz/HMAgA6h/rz2i29mNQnwnn691fHVG7MzBkzFvxcINk5MUkAQlRQVGoi7lQ2cSfPczg9lyPpeZzKLuTipLnubor2LZsR09qfG/tGMLpbGJ1D/Z2rb39NjqyF/LPQ9y5HRyKsJAlAuLzMvBLW709j7b40dp0+T5lJ46YgKsSPbm2aM6lvBJ1DA4hu7U9USz/X+JZfnR1zIDDSmK9fNGmSAIRLMpmNG7cLt59i67FMzNpoyrl3WEcGdQzmyqhg/Lzlv8dlzh6AU1th1N/BzYlvcrsIW60INhe4HsjQWl/WJ0wZ18dvAxOAQmC61nqXLeoWoi7yisv4dMdpFm4/Rcr5IsKa+/DHEdH8oXcbuoQ14fn4G8rW/xhz8fe729GRCBuw1VecT4DZwIIq3h8PxFh+BgLvWR6FaBB5xWUs+PkUH25JIqewjIEdgnlhwhWM7haKh7uLN+nUVuYRY/DX0MegWbCjoxE2YJMEoLXerJSKqqbIRGCBZR3g7UqpIKVUuNY6zRb1C1GVS0/8I7u25vFRMfRqG+To0OxPayi+AGVFgAafIPD0rd+KXVrDN38GLz8Y8rjNQxWO0VCNnBFAcoXXKZZtkgCEXRSWlvPJtpPM2exCJ/6yIjj2rTFD55k4yDoOZYW/L+PZDEJiILQHRA2HTtdAQFjN+967xFj2cdyr4NfSPvGLBtdQCaCyrxyVrkavlJoJzASIjHTC2ROFXRWXmVi84zTv/ZjIufxSRnRpxZOjOzv3ib8w22ib37UAinPAKwDaDTBO8M0jwKuZUa74AuRnQOZhOPq1MZUDQNsrodftxnKOlTXtJG2Cr56A9sNgwP0Nd1zC7hoqAaQA7Sq8bgukVlZQaz0HmAMQGxtbaZIQ4lI5haUs3nGa+dtOkpFXwtDolnwwugv927dwdGj2tXcprH8WSnKh20Rj4fT2Q2temtFshrMJxhXD/s9h3dPw9fPQ8WroMh7C+4CpDA6vgR0fQMtOcNt86fnjZBoqAawGHlFKLcW4+XtB2v+FtUxmzS8nslm99wyrdqdSVGZieEwIb03uw5BOTj5C1WwyTtpxcyFyMFz3BoR2r/3vu7lBeG/j56qnIT3BGN17eA2sfapCQQW9JxsLvvg6eTJ1QbbqBroEGAGEKKVSgL8BngBa6/eBdRhdQBMxuoHOsEW9wvVordl1Ooev9qayLiGNjLwSfD3dub5XOPcO70DXMCederkiswm+eAASlsOQx2Dk38Ddyv/KYT2Nn9GzjHsHWceMuX7Ce4O/E65lIADb9QKaUsP7GnjYFnUJ16O15kBqLl/tTWXNvjTO5BTh5eHGNV1acX2vNoy8ojXNvFxo0NZ3Lxkn/5F/heFP1Vi8TpSCkGjjRzg9F/pfI5qarPwSlu5MZkV8CifOFeDhphgeE8JTYzozulsoAT41tHM7o4QVsO0duPI+25/8hcuRBCAanSPpeXyw6Thr9qVRajIzqGMwD1zVkXE9wghq5sJLD+alG+3zbQfIPPzCJiQBiEbjSHoe73x/jLUJafh5uTN5QDumDmpPTKhM0QDAmiehvBgmvVdzLx8hakESgHC4MzlFvLL+MF/tTcXf24NHr43m3mEdXPvb/qUSvzemYR49S9rnhc1IAhAOU1xm4oNNSby3KRGt4ZFrorlvuJz4L2M2wYYXoUUUDHzQ0dEIJyIJQDQ4rTXr96fz8tpDnMkp4rqe4Tw/oSttWzRzdGiN0+5FkHEAbp0PHt6OjkY4EUkAokEdTs/l76sP8nNSFl3DAlhy/yAGd5K5ZapUkg8bX4Z2g4yRvkLYkCQA0SByCkt589ujLNp+iua+nvxjUg+mXNlOpmKuyc+zjeUXb19cv1k8haiGJABhVyazZskvp3l9wxFyi8q4a1B7nhzdWdr5ayM/A356x/jm3+5KR0cjnJAkAGE3O09m87cvD3AwLZeBHYJ56YbuXBHuAlM12Mqm14xun9f+1dGRCCclCUDYXPqFYl5Zf4hVe1IJD/Rh9h19ua5nOEqaMGov6zjEz4P+06Tbp7AbSQDCZs7ll/D+j8dZuP0UGnj02mgeGtHJtebpsQWtjSme3b3h6uccHY1wYvI/U1jt2Nk8Fvx8is93pVBcZuKmfm15fGQM7YKlW2e9HPwSEr+Fsf8PAkIdHY1wYpIARJ2Vm8wcPZvP5mOZrN+fzt7kHLw83Li+Vzh/HBFNdGt/R4fYdOWmGfP9hPWEAQ84Ohrh5CQBiN/JKSzlZFYhZ3OLycgr4XxBKecLSzlfUEp2YRnZBSUczyigqMwEQM+IQJ4b35Vb+7elpb8MUrJKaSEsn26s43vTR9bP8S9EDeQvzMWVm8xsO57FuoQ0fk7K4lRW4WVlArw9aOHnRQs/L0L8vYltH0zfyCBio4KJCPJ1QNROqOi8cfJP3gG3zIXWXR0dkXABtloRbBzwNuAOfKS1fuWS96cD/wbOWDbN1lp/ZIu6Rf0UlpazbGcyH209Qcr5Ivy9PRga3ZLJV0YS09qfsEAfWgV406KZF14eMljLbkxlRpv/dy8Z0z1P+h/0uMnRUQkXYXUCUEq5A+8CozEWf9+plFqttT54SdHPtNaPWFufsI7ZrFkRn8Jr3xzhXH4Jse1b8MKEK7ima2t8PGXBb7sqLTSWWsw8CueOwrkjcGKz8e0/tAfcMk8GfIkGZYsrgAFAotY6CcCy8PtE4NIEIBws7mQ2f//qIAlnLtC/fQvev6sfsVHBjg7LOeWmQfJ2OBMPmUeMn5zTgDbeV27G7J7Ro41v/DFjjYXahWhAtkgAEUByhdcpwMBKyt2slLoKOAr8SWudXEkZYQeplvn2V+81Bma9PbkPN/RuIwOzbKm0EJJ+hKPrIWkT5Jwytrt7QUhnaBsLfe6EVp0hpAsEdwRPH4eGLIQtEkBlZxF9yeuvgCVa6xKl1IPAfODaSnem1ExgJkBkZKQNwnNdRaUm5mz+bb79x0bG8ODVHWVglq2UFRuLtOxbDkkbjWkbvAKg49Uw8AFjBs/wXrJ6l2i0bHEmSAHaVXjdFkitWEBrnVXh5YfAq1XtTGs9B5gDEBsbe2kiEbWgtWZtQhr/Wnf41/n2nxvfVQZm2Urqboj/BPZ/ASUXoHkE9J8OncdB+6HgIRPdiabBFglgJxCjlOqA0ctnMnBHxQJKqXCtdZrl5Q3AIRvU63TMZk3qhSJKys2ENffBz7tuH4/Wmk1HM3nru2PsSc6hW3hz3rytNwM7ynz7VjOb4PAa2P4enP4ZPJvBFTdAnykQNRzc5Aa6aHqsTgBa63Kl1CPANxjdQOdqrQ8opWYBcVrr1cBjSqkbgHIgG5hubb3O5OS5AuZsSWLN3lRyi8t/3d41LIAx3UIZ3S2MHhHNq2yzzy8p56u9qXy64zQJZy4QEeTLKzf15NbYdri7STu/VbSGo1/D97Mg4yAEtYex/4K+d4JPoKOjE8IqSuvG28oSGxur4+LiHB2G3ZSbzMzemMj/Nh5HKbiuZzixUcH4ermRnF3E1sRzxJ3MxqwhPNCHYdEhdA4NIKiZJ+VmTXJ2IbtOnyf+1HnKTJrOof5MH9KBW/q3lb77tpCbCl89Ace+geBOcO0L0G2SfNsXjZpSKl5rHVurspIAHONCURmPLtnN5qOZ/KF3G1687gpaN7+8V0h2QSk/HM7g24PpxJ08T1ZB6a/vebgpYkIDGB4TwtjuYfSLDJKePbZydAN8fh+YSmHkizBgptzMFU1CXRKAdAdxgOOZ+dw/P47k84W8clNPJg+ourdTsJ8Xt/Rvyy392wLGXD35JeW4KUXrAG9ZUtEedsyBr581Bmfd+gm07OToiISwC0kADWzz0Uwe/nQXnu5uLL5vEAM61G0gVlAzL1lO0Z52fgTrn4Eu18HNH4KXn6MjEsJuJAE0EK01c386yctrD9I5NIAP746VbpmNzcHVsPZp6DwebpsvTT7C6UkCaAAl5SZeXLWfZXEpjOkWyn9u71PnLp7Czs4lwqo/QkR/o9lHTv7CBchZyM7O5hbz4KJ4dp/O4dFro/nTqM64SdfMxqW0EJZPM076t82XKRqEy5AEYEfxp7J5cNEuCkrKee/OfozvGe7okERl1j0DZw/AnSsgsK2joxGiwUgCsINyk5l3Nx7nnR+O0baFL4vuHUiXsABHhyUqs3sR7FkEV/0fxIxydDRCNChJADZ2OquQPy3bQ/yp80zq04a/T+xBoK+0JzdKafuM9Xc7XA0jnnN0NEI0OEkANlJuMvPJtpO8seEoHu6Ktyf3YWKfCEeHJapSlAPL7gbfFnDzxzK6V7gkSQA2sC8lh+dXJnAgNZeRXVsza1IPWSu3MTOVw8qZcCEZpq8F/1aOjkgIh5AEYIWzucW8ueEoy+OTCfH35r07+zGuR5hMx9CYaQ31y5DAAAAT80lEQVRrLPP7THgdIgc5OiIhHEYSQD3kFZcxZ3MSH25JwmTWzBjagcdHxdDcR9r6GzWzyejxs3shDH8aBtzv6IiEcChJAHWQV1zGwu2n+HjLCbIKSvlD7zY8M6YLkS1lRG+jV1oIqx6Cg6tgyGNw7V8cHZEQDicJoBYy8opZvP008346QW5xOVd1bsWTozvTp12Qo0MTtZF51BjolXEIxvwThjzq6IiEaBScMwGc2AwBbSAkut67MJs1WxPPseSX03x78CzlZs2YbqE8cm00vdrKib9JMJsh7mP47iXw8Ia7PofokY6OSohGwyYJQCk1DngbY0Wwj7TWr1zyvjewAOgPZAG3a61P2qLuSi2+DQbcZ3zbqwOzWbM7+Txf709nXUI6Z3KKCPbz4p5hHZh8ZTs6tvK3U8DC5jIOw1ePQfIO6DgCJv4PAqVbrhAVWZ0AlFLuwLvAaIwF4ncqpVZrrQ9WKHYvcF5rHa2UmoyxKPzt1tZdJQ9vKC+tsVi5ycyxjHx2n85h2/FzbE/K4lx+KZ7uimHRITw7vitju4fi7SF9xJuM0kL46S3Y+h9jKudJ70PvySA9s4S4jC2uAAYAiVrrJACl1FJgIlAxAUwEXrI8XwHMVkopba/lyDx8oLz4d5suFJZxMC2Xw+m5HE7L4/DZPI6k51JcZgYgtLk3V8W04uourbima2vp0dPUaG3c4N3wotG/v8ctMO4V6eMvRDVskQAigOQKr1OAgVWVsSwifwFoCZyzQf2XKcGDQ0mpLFmxj7TcYk6cyyc5u+jX91s086RrWHOmDIikd9sgercLIqplM+m/31Sl74evn4OTWyC0J9z4PkQNc3RUQjR6tkgAlZ01L/1mX5syRkGlZgIzASIjq14qsTopxT54lh7lh4IMQpt706ttEHcMaE/3Ns3pGhZAqwBvOdk7g8Js2PgyxM0Fn0C47k3oP12mdRCilmyRAFKAdhVetwVSqyiTopTyAAKB7Mp2prWeA8wBY1H4ugZjMmu2lcUw2f17dj5/jZwMnJGpHOLnGSf/4ly48j4Y8Tw0q9vymkK4OlusKL4TiFFKdVBKeQGTgdWXlFkNTLM8vwX4wW7t/8CQ2Fg8dSkUX7BXFcIRtIbDa+G9wbDuaQjrCQ9uhQn/lpO/EPVg9RWApU3/EeAbjG6gc7XWB5RSs4A4rfVq4GNgoVIqEeOb/2Rr662Ku5uiU2RbiAeKzsuJwVkk74RvX4TTP0PLGLh9MXS9Tnr3CGEFm4wD0FqvA9Zdsu2vFZ4XA7faoq5a8W1hPBblNFiVwk7OHoQf/wWHVoNfa7j+Leg7FdydcwyjEA3JOf8X+QQaj8XnHRuHqL8z8bDlTTi8Brz8jTb+wY+AtwzGE8JWnDMBeFrm4i8rrr6caFxykuHIeti9ANITwDsQrn4WBj4oTXlC2IFzJgAPSwIolwTQaJjNkLITTm+D3DQoyYWSPONGfUkeFGVDzmmjbFhPY67+Xrf9djUnhLA550wAnj7GY1lR9eVEw0jfD6sfgdTdxmvvQPANBO/m4B0A/q2hZScYMBOiR0Prro6NVwgX4ZwJQK4AGo+jG2D5dKPt/g/vwBV/kOYcIRoJ50wAcgXQOJzeAcumQqsuMOUzaB7u6IiEEBU4ZwKQKwDHyzoOSyZD8zZw10rwC3F0REKIS9hiJHDj4+4Bbh5yBeAoBedg0c3GIK07V8jJX4hGyjmvAMC4CpAE0PBKC41v/nlpMG2NcXNXCNEoOW8C8PSBckkADcpsgpX3Q0oc3L4Q2l3p6IiEENVw4gTgCyX5jo7CtXzzgjFyd9yrRm8fIUSj5pz3AMCYMCzrmKOjcB3b/gs73oNBD8OgBx0djRCiFpw3AQSEGTcjhf3tmAMb/gLdJsKYfzo6GiFELTlvAqhkXWBhB798COufga7Xw00fgZvz/kkJ4Wyc9x6Ahw+Ulzg6CuelNWz+t7EqV5cJcMs88PBydFRCiDpw4gTgLVcA9mIqh7VPwq750Gsy3PBfOfkL0QRZlQCUUsHAZ0AUcBK4TWt92ST8SikTkGB5eVprfYM19daKhw+Yy42uibIusO0UnYfP74fEb2H4U3Dti7IqlxBNlLUNts8B32utY4DvLa8rU6S17mP5sf/JH4wrAJBmIFtK2wsfXA1JP8L1/4GRf5WTvxBNmLUJYCIw3/J8PjDJyv3ZjodlQjhpBrKe1hA3Dz4eY1xVzVgPsfc4OiohhJWsvQcQqrVOA9BapymlWldRzkcpFQeUA69orVdZWW/N5ArANvLSYfVjcOwb6DgCbv5Y5vYRwknUmACUUt8BYZW89UId6onUWqcqpToCPyilErTWx6uobyYwEyAyMrIOVVxCrgCsozXsWwZfP2vMqTTuVWPBFunmKYTTqDEBaK1HVfWeUuqsUirc8u0/HMioYh+plsckpdSPQF+g0gSgtZ4DzAGIjY3VNR5BVdw9jUdTWb134bLS98O6Z4zlGyP6w6T3oVVnR0clhLAxa7/OrQamWZ5PA768tIBSqoVSytvyPAQYChy0st6aXWwCMrlgE5DZDMd/gLi5cCbe+DZfG+eOwcqZ8MFwyDwMf3gb7v1OTv5COClr7wG8AixTSt0LnAZuBVBKxQIPaq3vA64APlBKmTESzitaa/snAPeL9wBK7V5Vo1KSB5/dZfTUuahNPxjyqDFB28Uro4tM5ZC0EeI/gSPrjKazwQ/DsCdl6UYhnJxVCUBrnQWMrGR7HHCf5fk2oKc19dTLr01ALnQFoDWs+iOc2AITXoeY0ZD4vTFR24oZ4B8Gna41VukqK4SsRDj1M5TmQbMQGPq4MZmbfytHH4kQogE490hgAJMLXQEcXAWHVsOov8OA+41tV94L/adD4newa4HxWHjO+KYf1B563gLRIyFmrIzmFcLFOG8CcLUmIFMZfD8LWnczmnsqcnOHzmONn4u0lkFcQrg4J04ALtYEFP8JZCfBHctqN/WFnPyFcHnO26nblZqASvJg06vQfhjEjHF0NEKIJsL5rwBcoQlo22woyIQpS+WbvRCi1pz3CsDdRcYB5KUbvXyuuAHaxjo6GiFEE+K8CeDXJiAnHwm88WWjmWvUS46ORAjRxDhvAvi1CciJrwDOHoDdi4wuny07OToaIUQT48QJwMmbgMwm+Opx8AmEq55xdDRCiCbIiW8CWwY1OetN4O3vQcpOuOlDmbJBCFEvznsF4OYGbh7O2Q00PQF++Ad0Hg89b3V0NEKIJsp5EwAYzUDOlgCKzhuTvfm2gBvekW6fQoh6c94mIDDmtnGmBFBeCivugQtnYMY68K9qATYhhKiZcycAdy/n6QVkNsEXM415/m+YDe0GODoiIUQT5wJNQE4wDsBUZkzzfOALGPNP6DfV0REJIZyAk18BeDb9bqClBbB8OhzbANf+5fKZPoUQop6sugJQSt2qlDqglDJbVgGrqtw4pdQRpVSiUuo5a+qsE48mfhM4PwMWTDTm8L/+LenvL4SwKWubgPYDNwGbqyqglHIH3gXGA92AKUqpblbWWzvuXk13HEBKPHxwtbFA+63zIXaGoyMSQjgZqxKA1vqQ1vpIDcUGAIla6yStdSmwFJhoTb215u5l+yYgswkOrTFO0PagNcTPh3njwN0D7t0A3W6wT11CCJfWEPcAIoDkCq9TgIENUK8lAZTbbn9mMyyfBoe+Ml7ftgC62TCXFWbDmifg4JfQcQTcMk9G+Qoh7KbGBKCU+g4Iq+StF7TWX9aijspGKulq6psJzASIjIysxe6r4e5p3ES1ld0LjZP/1c/BsW9g3f9B9Cjw8rN+30mb4IsHoSADRv7NWKC9Nit7CSFEPdWYALTWo6ysIwVoV+F1WyC1mvrmAHMAYmNjq0wUteLuBeYcq3bxq8Js+O4liBwCI54zvqHPGwdx82DII/Xfb2mhMaXzz7OhZQxM+Q7a9LVNzEIIUY2GGAewE4hRSnVQSnkBk4HVDVCvpRuojcYBbH4dinNgwmvG9AvtB0PUcGMxlvoONkvaBO8NMU7+sffAA5vk5C+EaDDWdgO9USmVAgwG1iqlvrFsb6OUWgegtS4HHgG+AQ4By7TWB6wLu5bcPW3TDfT8Kdj5IfS5E8J6/rZ9+FOQnw57Pq3b/orOw5cPwwLLzd27V8P1/7FNU5IQQtSSVTeBtdZfAF9Usj0VmFDh9TpgnTV11Yubja4AfvgnKDcY8fzvt3ccAW36wU9vQd+pRq+d6mhtjOZd/ywUZsHQJ4zmJE9f62MUQog6cvKpILysTwCpuyFhGQx6CAIjfv+eUnDV03D+JBxYWf1+zh6A+X+AFTOgeTjM3Aij/y4nfyGEw7jAVBBWNAGZTbDmT+DX2vi2XpnO46HVFbDlDeh+429LUV6Um2rcP4ifZ6zeNeF16D+j5qsFIYSwM+c+C7l7gtmKK4BfPjSuAG7+GHyDKi/j5mbM0fPZnca3+wlvGCf61N2wd4nxo80Qey9c82fp1y+EaDScPAFY0QR07hh8Pws6jYQeN1df9orrYczL8O1ffxskBuDhA72nwPAnoUVU/eIQQgg7ce4EUN8lIcuKjRk4Pbxh4uzarbo15BHoPBaOrIfyYgjpDNEjwTug7vULIUQDcO4EcPEKQOvaL52oNax7Gs7uhzuWQ/M2ta8vJMb4EUKIJsD5ewGhjZu5tbX1TWPKh6uegc5j7BaaEEI4mpMnAMsFTm2bgfYtN9r9e94K17xgv7iEEKIRcPIE4GU81qYn0OG1sOpBaD8UJr5b+yYjIYRoopw7AbhZ+uTX1BPoyNewbBqE94EpS42bv0II4eScOwFcHJRVXRPQ8Y2wbCqE9YCpK8GnecPEJoQQDubkCcDSBFTVFUDqblh6J7SMhrtWGgO4hBDCRTh5AqimCejQGlgwCZq1NE7+MkJXCOFiXCQBXNIEtG+Z0ewT3AHuXmVMziaEEC7G+QeCwW+9gEzlxnQN2981FnOZshS8/R0XnxBCOJBzXwFc7AVUbrkC+GGWcfK/8n64c4Wc/IUQLs25E4BvC+Nx8c2wawH89LYxFfN1r4Onj2NjE0IIB7N2SchblVIHlFJmpVRsNeVOKqUSlFJ7lFJx1tRZJxfb9osvwOpHIaI/jHulwaoXQojGzNp7APuBm4APalH2Gq31OSvrqxv/UPANhqJs4/Udy+WbvxBCWFi7JvAhANVYp01w94RnTxiLtrt5gl9LR0ckhBCNRkP1AtLABqWUBj7QWs+pqqBSaiYwEyAyMtI2tfe5wzb7EUIIJ1JjAlBKfQeEVfLWC1rrL2tZz1CtdapSqjXwrVLqsNZ6c2UFLclhDkBsbKyu5f6FEELUUY0JQGs9ytpKtNaplscMpdQXwACg0gQghBCiYdi9G6hSyk8pFXDxOTAG4+axEEIIB7K2G+iNSqkUYDCwVin1jWV7G6XUOkuxUGCrUmov8AuwVmv9tTX1CiGEsJ61vYC+AL6oZHsqMMHyPAnobU09QgghbM+5RwILIYSokiQAIYRwUZIAhBDCRSmtG29Xe6VUJnCqnr8eAjTs1BOO5WrHC3LMrkKOuW7aa61b1aZgo04A1lBKxWmtq5ygztm42vGCHLOrkGO2H2kCEkIIFyUJQAghXJQzJ4AqJ5xzUq52vCDH7CrkmO3Eae8BCCGEqJ4zXwEIIYSohtMlAKXUOKXUEaVUolLqOUfHY0uVLa2plApWSn2rlDpmeWxh2a6UUu9Y/h32KaX6OTb62lFKzVVKZSil9lfYVudjVEpNs5Q/ppSa5ohjqa0qjvklpdQZy2e9Ryk1ocJ7z1uO+YhSamyF7U3ib18p1U4ptVEpdciypOzjlu1O+zlXc8yO/Zy11k7zA7gDx4GOgBewF+jm6LhseHwngZBLtr0GPGd5/hzwquX5BGA9oIBBwA5Hx1/LY7wK6Afsr+8xAsFAkuWxheV5C0cfWx2P+SXg6UrKdrP8XXsDHSx/7+5N6W8fCAf6WZ4HAEctx+W0n3M1x+zQz9nZrgAGAIla6yStdSmwFJjo4JjsbSIw3/J8PjCpwvYF2rAdCFJKhTsiwLrQxkJB2ZdsrusxjgW+1Vpna63PA98C4+wfff1UccxVmQgs1VqXaK1PAIkYf/dN5m9fa52mtd5leZ4HHAIicOLPuZpjrkqDfM7OlgAigOQKr1Oo/h+5qbm4tGa8ZelMgFCtdRoYf2RAa8t2Z/q3qOsxOsuxP2Jp8ph7sTkEJztmpVQU0BfYgYt8zpccMzjwc3a2BFDZ6vTO1M1pqNa6HzAeeFgpdVU1ZZ393wKqPkZnOPb3gE5AHyANeMOy3WmOWSnlD3wOPKG1zq2uaCXbnOWYHfo5O1sCSAHaVXjdFkh1UCw2pyssrYmxDsMA4OzFph3LY4aluDP9W9T1GJv8sWutz2qtTVprM/AhxmcNTnLMSilPjBPhYq31Sstmp/6cKztmR3/OzpYAdgIxSqkOSikvYDKw2sEx2YSqemnN1cDF3g/TgC8tz1cDd1t6UAwCLly8vG6C6nqM3wBjlFItLJfUYyzbmoxL7tfcyG/LqK4GJiulvJVSHYAYjJX2mszfvlJKAR8Dh7TWb1Z4y2k/56qO2eGfs6Pvjtv6B6PHwFGMO+UvODoeGx5XR4w7/nuBAxePDWgJfA8cszwGW7Yr4F3Lv0MCEOvoY6jlcS7BuBQuw/i2c299jhG4B+PGWSIww9HHVY9jXmg5pn2W/+DhFcq/YDnmI8D4CtubxN8+MAyj2WIfsMfyM8GZP+dqjtmhn7OMBBZCCBflbE1AQgghakkSgBBCuChJAEII4aIkAQghhIuSBCCEEC5KEoAQQrgoSQBCCOGiJAEIIYSL+v+doGI0hZBIfgAAAABJRU5ErkJggg==\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ - "cls_hyrec=LambdaCDM_newAlpha.lensed_cl(2500)\n", + "cls_hyrec=LambdaCDM_hyrec.lensed_cl(2500)\n", + "cls_newAlpha = LambdaCDM_newAlpha.lensed_cl(2500)\n", "\n", "ll_vec_hyrec = cls_hyrec['ell']\n", "\n", @@ -164,8 +180,19 @@ "\n", "# ax.set_yscale('log')\n", "\n", - "plt.plot(ll_vec_hyrec, (cls_hyrec['tt']/cls_std['tt'] - 1.)*100)\n", - "plt.plot(ll_vec_hyrec, (cls_hyrec['ee']/cls_std['ee'] - 1.)*100)" + "plt.plot(ll_vec_hyrec, (cls_newAlpha['tt']/cls_hyrec['tt'] - 1.)*100)\n", + "plt.plot(ll_vec_hyrec, (cls_newAlpha['ee']/cls_hyrec['ee'] - 1.)*100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.loadtxt('../hyrec/Alpha_inf.dat')\n", + "a = a/1e5\n", + "np.savetxt('../hyrec/Alpha_inf.dat', a)" ] }, { @@ -192,7 +219,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.6.9" }, "toc": { "base_numbering": 1, diff --git a/python/cclassy.pxd b/python/cclassy.pxd index d99da503b..15deea684 100644 --- a/python/cclassy.pxd +++ b/python/cclassy.pxd @@ -29,6 +29,8 @@ cdef extern from "class.h": cdef struct precision: double recfast_fudge_H + FileName hyrec_Alpha_inf_file + # FileName hyrec_R_inf_file ErrorMsg error_message cdef struct background: @@ -69,6 +71,7 @@ cdef extern from "class.h": cdef struct thermo: ErrorMsg error_message + int thermodynamics_verbose int th_size int index_th_xe int index_th_Tb @@ -86,8 +89,6 @@ cdef extern from "class.h": double rs_d double YHe double n_e - char * ALPHA_FILE - char * RR_FILE int tt_size diff --git a/python/setup.py b/python/setup.py index 74a9b7c94..60614c49c 100644 --- a/python/setup.py +++ b/python/setup.py @@ -11,8 +11,12 @@ GCCPATH_STRING = sbp.Popen( ['gcc', '-print-libgcc-file-name'], stdout=sbp.PIPE).communicate()[0] + GCCPATH = osp.normpath(osp.dirname(GCCPATH_STRING)).decode() +# For some reason the path must be corrected to the below: +GCCPATH = '/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/lib/clang/11.0.0/lib/darwin' + liblist = ["class"] MVEC_STRING = sbp.Popen( ['gcc', '-lmvec'], @@ -33,6 +37,7 @@ VERSION = line.split()[-1][2:-1] break +# also need the extra_link_args below. setup( name='classy', version=VERSION, @@ -43,7 +48,7 @@ include_dirs=[nm.get_include(), include_folder], libraries=liblist, library_dirs=[root_folder, GCCPATH], - extra_link_args=['-lgomp'], + extra_link_args=['-lgomp', '-Wl,-rpath,/usr/local/opt/gcc/lib/gcc/9/'], )], #data_files=[('bbn', ['../bbn/sBBN.dat'])] )