Skip to content

Instantly share code, notes, and snippets.

@smsharma
Last active October 12, 2020 16:18
Show Gist options
  • Save smsharma/fac2c92484bf9d8cbc9b508e2c300be8 to your computer and use it in GitHub Desktop.
Save smsharma/fac2c92484bf9d8cbc9b508e2c300be8 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The autoreload extension is already loaded. To reload it, use:\n",
" %reload_ext autoreload\n"
]
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import arviz as az\n",
"from tqdm import *\n",
"\n",
"%matplotlib inline\n",
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"samples = np.loadtxt(\"sample.dat\")"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [],
"source": [
"m_min, theta_min = np.min(samples, axis=0)\n",
"m_max, theta_max = np.max(samples, axis=0)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [],
"source": [
"# Cuts in theta\n",
"theta_cuts = np.linspace(theta_min, theta_max, 20)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/smsharma/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0\n",
"Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`\n",
" \"\"\"\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "777fa1bcebd346aab090420bf1d82d94",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(IntProgress(value=0, max=19), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"# Calculate median and 68/95% HPD intervals for different cuts\n",
"\n",
"median_ary = []\n",
"onesigma_ary = []\n",
"twosigma_ary = []\n",
"\n",
"for theta_cut in tqdm_notebook(theta_cuts[:-1]):\n",
" samples_cut = samples[samples[:, 1] > theta_cut][:, 0]\n",
" \n",
" median_ary.append(np.median(samples_cut))\n",
" onesigma_ary.append(az.hpd(samples_cut, 0.68))\n",
" twosigma_ary.append(az.hpd(samples_cut, 0.95))"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, '$\\\\log10 M$')"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd5zcdZ3/n5/ps3W295Ye0kmB9NCCCCocKooFBEE90dPzUE/u51nwPM/T81RUsB3YxYYIgkLoJZVAes8m2d1k+07Z6d/P74/PTHYJm2TLzM53dj/Px2Oys9/vzHfemZ35vr6fdxVSSjQajUajORuWTBug0Wg0GnOjhUKj0Wg050QLhUaj0WjOiRYKjUaj0ZwTLRQajUajOSe2TBuQDkpLS2VjY2OmzdBoNJqsYuvWrZ1SyrIzt09IoWhsbGTLli2ZNkOj0WiyCiFE81DbtetJo9FoNOdEC4VGo9FozokWCo1Go9GcEy0UGo1GozknWig0Go1Gc060UGg0Go3mnGih0Gg0Gs050UKh0Wg0mnMyIQvuNBrNREMCrwGvjOkoMSPOH/Zs5Whv55gtKsup5Ia5XyDHXjzmY5kdLRQajcbkdAN/BvYBHkbrCAnFoty3dSNbWk9Q7HaP2ao9Hfto7ns/n7z4XgpdNWM+npnRQqHRaExKHNgEPArYgaZRH8kfCfHtjS9yoKuDJk8VQogxW+dxFdHibeE/nnsf/7LiPspyp435mGZFxyg0Go0JaQN+ADwMVACVoz5STzDA155/nMPdHdQXFqdEJJLUFNTgDfu5+7mbOOHdnrLjmg0tFBqNxkREgCeA7wB9QCPgGPXRTvm9fOW5v9Ie8FFbWJRSkUhSkVeBYUi+8uxtHOx+LuXHNwNaKDQajUk4Avwv8DRQC5SO6WjNvV18+dlHCEYjVOUXjt28c1CSU4LL5uKrz32M107+Ja2vlQm0UGg0mgzTD/wJuC/xez1jDZ/u6zzJV557FKsQlOXmj9G+4VHo8lDs9vCNl/+Vl47/fFxec7zQwWyNRpMhJLAH+CMQBhpIxbXrtrZmvrPpKYpcORQ4x57dNBJyHflUWWx8b/NX8UW6uGLKx9Pi7hpvtFBoNJoM0As8AuxABavfMFRtVDzTvJ+fbHuBirwCcuyjj22MBZfNTV1hDQ+8eh++cDfXzf48FmHNiC2pQguFRqMZRwxgG/AXQKBSXsd+xS2l5JEDO/jNzi3UFnhw2uxjPuZYcFidNHrq+NPe3+EL9/LeBf+FzeLMqE1jQccoNBrNONEO/Bj4HVAMVJEKkTCkwW92beE3O7dQX1iccZFIYrPYaSpq4KmjG7h3y0cIx/yZNmnUaKHQaDTjwHOojKZTwBTAlZKjRuNxfrLtBR45sINGTwl2q7lcPBZhpdFTz5bWLXzr5Q8QiHRl2qRRoYVCo9GkmX2oeEQ1UJ6yo4ZiUb676SmeP3aIJk8pVos5T2dCWKgvrOdA90H+64X30Rs6nmmTRow531mNRjNBCAJ/QAlE6lxCvnCI/37xb7zW3kKDpxiLyTOLhBDUFtTS5j/FV557P+2BfZk2aURoodBoNGlkA6pOIi9lR+wOBvjP5/9Kc2839QXpqbZOF9X51QTCAe5+9gMc79uWaXOGjRYKjUaTJo4DL6BcTqmhzdfH3c8+SlcwQE2BJ2Mi0X7Ei787PKrnludVAHD3s7dzoOuZVJqVNrJCKIQQVwghPieE+LwQoiLT9mg0mvMRBX4PFAKpCTCf8PZw93OPEIlFqcxLb0uOcyGl5Gd3buTe259FSjmqYxS7S8h15PDV5z/O/q6nU2tgGhi3OgohhAf4R+Ay4JCU8nYhRCNwA/BO4A9Syq8M8bxpwNeBxcAlwAPAleNktkajGRUvAR2oauuxE4iE+fbGDQgEpbmpc2ONBiEEb//8hThctjGtaAqchYRiYZ44/BNmlKxLnYFpYNxWFFLKXinlf6DWo47EtqNSyq8BPs4e6boeaJVSxlFdw9YLIVKXOqHRaFJMB/A3UuVyMqTB/21/ia5+P6U5mRWJtgN9RIIxGuaXUDVDrWqkMbpVBUCRq5hXT+0kZvSkysS0kA2up2mo1AmAnkHbNBqN6TBQDf6cpCrL6YnDe3n5xGFqC4pScrzREovE+dXnNvGHu9U4Viklj/zPDv7836+N+ph2q41oPM6xvmdTZWZayIYWHoMdnJEhtgEghLgduB2gvr5+HMzSaDRv5BXgMGOZRjeYg93t/HLHJurSNEtiJNgcVt7x74uxu9TpRwiBu9BBPBJHSjkG+xzsbv8rU4reljpjU0w2CEUHA+OtnIO2vQ4p5X0k+hQvWbJk9GtBjUYzSryoHk5VKTlaXyjIdzY+hcflxmHN7KnKiEssVkHd3OLXbb/0lpljPrbHVcjG1p1cM7OTsc7gSBemdD0JIRYKIZKB7SeBaqHkuh5oAfZnzDiNRnMW/opqHT729hwxI859W58lGIvgceWM+XhjIR4z+MkdL7Dpj0fP+piTB/vY98LJUR0/3+HiRJ+PvtCmUVqYfsZNKIQQFiHEF4CFwEIhxBeEEOWJbY3AOiHEvyQePhVYDyCl/BtKLO4GPgO8T0ppjJfdGo1mOOxHuZ1GP9t6MA/ve40d7S1UZTANNkk0FMdTlUN+6dm7v/7t+3t44r69owpsq2tgJwe7H0cJrfkQo80DNjNLliyRW7ZsybQZGs0kIYhq+OcgFRXYO06d4Osv/o36wmJsFnM1+TsbvSf7cebacOePbgbGKb+XC6tyuW3xfagxsJlBCLFVSrnkzO2mdD1pNJpsYgMQIBUi0RHwcc/mpynPzc+4SEgpee4XB/F2BM/7WE9lDu58B1JKgr7IeR9/JkXuHLa1tRM3Xh2NqWlHC4VGoxkDqWvTEYnH+N7mZxAI8hypaUM+FrpPBHj2gf3sfqZt2M/5w92v8MvPbhqxC8phtRGO2zjh3QDERmhp+smGrCeNRmNKUtemQ0rJb3du4UhPB41F5sj8KanL46P3r6OgbPiiNXNlJUFvBClHM5LJwt7OFho8zagwrXnQKwqNRjNKXkJNrRt7IdymliM8fmg3dYXF53/wONDdEgCUS8liHf5pcu6l1Sy9thGLdeQyUejMYWNLO2pUrLnQQqHRaEZBsk1HzZiP1OLt4Yfbnqc6v9AUw4da9vTw3fc/xc4NLaM+xr4XT/H8rw6O6DkFThdHeoL4wluA0KhfOx1k/q+i0WiyDAN4iFS06eiPRvjupqdwWm247aPLGEo15U0FrL1pBtMvHn2j6oMb29n5RAvx6PAz+VWarOBQTxdwaNSvnQ50jEKj0YyQV1AnsrG16ZBS8rNXX+aU30e9xxwuJwC7y8ra988Y0zGu+PBsrHYLVtvIrsUdVhvbT/axsHIjMGdMNqQSvaLQaDQjIHVtOp46so/njx2ktjCzzf6S9J7s5/8+8SKdx/xjPpbDbcNqsxCPGnQc9Q37ecXuHLa0dmDIg6j32hxoodBoNCMgNW06Dvd08LPXXqa2oMg08657T/bj6wydbvqXCv741e387M6XiYbjw3q802anPxql1ecDzDNXW7ueNBrNMEm26Wgc01G84SDf3vgUBU43Tpt5TkGNC0v56P2XjCpj6WysuGEKC9bXYHeORHwk+zoj1BZsBJamzJaxoFcUGo1mGASBPwDljOW0ETcMfrTtefyREEXuzDb7SxLyR9nxRAtSypSKBED1TM+Ig+L5DjebWtqBVlT6cebRQqHRaIZBatp0PHpgB9tPnqAm35MSq1LB1oeb+eNXX0lJbOKsr/GXZh78wtZhzdgudLk50N1Of9QA9qTNppFgnnWfRqMxKceB51Fd/kfP7o5WHty9jbqCzA8hGsyKG6ZSN7eYsob8tL1GLGIQ7o8RCcZx5pz7tKtiNpJD3ZJ5FS8Dq8n0Nb0WCo1Gcw6SbTo8jKVNR1e/n3s2PU1ZTh52qzk6wsZjBvGogcNto35eetNzl13byLLrGoctkFZh5bVTncyrsKFG8NSl1b7zoV1PGo3mLMSAvzPWNh3ReJzvb3mGmGGQ78x8s78kLz94mHtufppATzjtryUsAiEE/X0R9r146ryPL3bnsLm1GSltwOhncqcKLRQajWYI+oD/A55jLFezUkp+u2sLB7raqcwrSJFtqaF+XjFzLqkmt+jsA4lSzZM/2svvv7ztvK3I3XYH3nCQk34HsBW1sssc2vWk0WjO4BDwy8T9sVVfb2o5wmMHd9HoKTFVXAKgbm7xG2Zgp5t1N89g2XWNwxpwJIH9Xd1U5buAo8D0NFt3dvSKQqPRJIgDTwE/BHKB0fc6AjhhsmZ/SXY/08Yz9+8nHhv9RGVDSrzh8w80OpP8EhcVU9TKKhY5dxFevsPJppZmwE2mO8qa56+n0WgyiA94ANURtp6xpsEGImG+vXEDLpvdNM3+khx9tYt9L54a9QonHItxtKeTaDzO0Z5O4sbIBefF3x7iBx989pxiUeh0s7ezjVCsANiJqmXJDNr1pNFMeo4Cv0AFr8fmagIwpMFPXnmBrn6/aeZLRENxOpp9VM/08OaPzyUSjI2quK431E9vKMhNC5ezqn4av9u9jccP7aIyr5CcEQhi5dRCGheWEI8a2BxDZ4FZLRakhMM93VxQZgAHgXkjtjkV6BWFRjNpMVDB6ntRLcPH3ugP4LGDu9jUepTaAnM0+wP489df5Ref2XS655LDPbJrZCklrb5eJPD/1lzNZVNm47TZuXHeMj627BJ6Q/20+4ffxG/K4lKu+ef5OHPP3abdYhHsaG8BCoBNI7I5lWih0GgmJX7g58AjQC3qRDR29nS08ZudW6gvKM5o8Drkj/Lszw7Q36eyi1a+exrv/OLiEfZcUsSMOEd7u5heXM6X1r2VqcVlp/cJIVha08SXLnkrpbn5HO3twpDDd0W1H/Gx9S/NZ91f7M5lc8tRpCwADgO9I7Y/FWih0GgmHceA76JOPE2MdfhQkq5+P9/d/BQlObkZL6rzdoR4+v/2cXCT6pVUOa2AhgUlIz6OPxLmWF83b5u1gE8uv5xCl3vIx1XmFXLX6qu4tGkmR3u76I+eO/01yZaHm3nyh3sJ98eG3O+22enqD9DRH0BN4c5MR1kxnN4jaTVAiCtQLRJtwL1SyjdUowghnEA+yonql1IO/a4mWLJkidyyZUs6zNVoshgDeBk1T6IIKEzZkSPxGF97/jFOeHupyk/dcUfC3+9VfZGu+NBsQLUN91SOvvHgKb8XA8lHFq9lYdXwakmklGxqOcKPtr2Aw2qlLPfcbUGSK56cwrPHN5p7u/jghatYWV+Jqo7/p+H+F0aMEGKrlHLJmdvHbUUhhPAIIT4nhHhSCHFfYts04OvA14AXUWkXQ+FHDek9DKwYD3s1molFP/Br4GGgmlSKhJSS3+7cwqGejnEvqgt6B67cw4Eo4UD0dOO90YqEIQ2ae7sozcnjS+veOmyRAOWKuqh2Cl+65K0UuXNpPo8rKqfQcVokzjY2NdfhZHPrUVQmWjtw/sruVDNuQiGl7JVS/geqw1hSPq8HWqWUceAIsF4IUT7E03+G6ozVIKV8dlwM1mgmDC3APSi3RSMDX7/U8NKJwzx+aDd14xyX2LmhlW+8/Qm6WwIAXP3JeVzzz/PHZEMoFuVITxerG6Zz15o3UzFK4avKL+Tf1ryZNQ3TOdrbRfAcrigpJb/9/Bb+/PVXh9zvcbnZ2d5KJB5DnbJ3jcqmsZDpGMU0BpKDewZtO5NqlHvqDiFE9VAHEkLcLoTYIoTY0tHRkXpLNZqsQwIbge8l7teg/Nyp41hfNz/e9sK4FdWdOuw93Q68cWExS97WcDqDaawi1R0M0B7wctviVXxg4QpctrHFblw2OzcvXMGHl6ylKxigIzD0SFQhBFUzPVRMHVqUbBYrhpQc7e0CSlF/09EXC46GTAvF4IhXZIhtCCEsqNFa3wLagMfEEJ8IKeV9UsolUsolZWVlZ+7WaCYZAeC3wB9Raa+pn//gj4T4zsYN5NjHp6guFonzf594iWfu3w9AXrGLN310DnnFY+vVJKXkhLcHm8XK59dew5qGGSlbGQkhWFE3lS+uewsFLjfNfd0YQ8SFV79nGitumHr24wC72ltRVdp+4ERK7BsumRaKDtT/HFQid3LbYG4B6qRyPLahKk5Gnr6g0UwKDFS30f9BuSgaSbWrCQaK6npCAUpyxlbFPVxsDivv/OJi3vxPc1N2zGg8zpHeTuaUVfPFdW+h0VOasmMPpqagiM+vuZqVtVM40tNJKPbGJn9SSg5t7hhygFKRO5eNLUcS8Rc7MLSbKl1kWiieBKoTK4R6lDN1vxBioRDiK4nHBIBkXKISFcnpHndLNRrT0w78BPgVKvBZS7q+4o8e2MnmlmZq8tNfVPfED/ew+U9HAWhaVIq7IDXC5wuHOO7t5h0XLObjF12a9hbobruDWy9cxe2LV9MR8NHZ/3pBCAdiPPjFrbz04OE3PDfX7uCU30t3MACUoXo/DS8FNxWMWwuPhAvp88DCxO9fAO5GicXdwFTgfVJKQwgxFVgP3AU8CHxDCPElYD7wdilHUNGi0Ux4IqgK6w2oBfrY23Cci90drTy4ayv1hekPXhtxg44jPsKBc2bEj5iT/j6EEHxm5ZXMKa9J6bHPhRCC1Q3TafSU8N1NT3Gst5vawiIsQuDKs/O+/76YiilvTKkVQiCBA13tlORMAcKo1iszxsfuTNdRpANdR6GZPBxExSH6ULGI1BTPnY3Ofj//76mHyLE5xm0IkRGXgMRiTc3q6HhfN/WFJXx02TpKx8ltNhT90QgPvPoSL584/AaXl5TyDSLcHvAxp7yKjy69BOWhbwJuTKlNGa+j0Gg0qaQXVRfxI1Sos550i0QkHuOeTU8jpUy7SHQe8/PgF7YS8kexWEXKRKI94KWmoIjPrroyoyIBkGN3cMuildTke17nhupo9vGDDz5L677Xt+sodufw2skWYkYcKAZ2o+pj0o8WCo0mq4ijqqu/BexFBavTX+QmpeRXOzdzpKeDyrz0V16fOuzl+K4egr7UTXYLxaJE4nE+vGQNzjGmvqYKh9XGbYtX44+EicZVw8KCMjfOHBuR4OvdbTaLlagRo7m3G5UcKoED42KnbjOu0WQNJ4A/oXI+qhlIFEw/Lx4/xBOH9tDoGZ+EwznrqplxcQV2V2p6RkkpafX2cuviVVTnpz5VeCw0ekp528wF/Gnfdpo8pThzbNzynZVDPlZgYU9Ha6IxYbKj7IK026hXFBqN6elH9We6B5VD38R4ikRzbxc/fiX9RXVSSh67ZxdHt3cCpEwkAFp8vSytaWRNfebGiZ6Lq2fMo/YMF1Q8ZrD/pde36/C43GxsOZr8DRXQ7iHdaKHQaEyLBHagaiI2oeIQ4zsIyBcO8e1NG8i1O9JeVBfyRzm0uYOjr6Y2+90bDuK2O7hp4XLTze1OMpQL6pVHj/Orz22mZe9ArCLP4eSEt4feUD8qNiVQLsj0ol1PGo0p6UA18NuPml09vt0G4obBK23H+O2urfSFguMyhMid7+C2769K6UoiZsTp7A/wr6veRIFz6BbhZuFMF9SC9bUUlLmonjkQE0oK3cHudpZUN6IuHF4GLibV7VkGo1cUGo2p6EOVFn0LFYuYAuSO26sb0uDVk8f5fxse4jubniJqxNMuEq37evn7D3ZjxCUOty2lV/0nvD1cM2Mes8tSM70v3Qx2QdldVmYsr3jD++Gy2dnWdizxWx7QSbo7ymqh0GgyjoHyNf8K+C/gKVSwumLcLJBSsqu9lS88/Re++dIT9MciNBWVnnVQTyo5sLGdXU+3EfKnLsMJVCpsfWEJ185amNLjppOhXFDbHzvObz6/5XT79CJ3Dq+0HSduJOuOrcDOtNqlXU8aTcYIonLhn0FdFeYAdYzn9ZuUkv1dp/jd7q3s72qn0OWm0VMyrr78te+fwbJrG1PWmgNenwrrsGbXae5MF1Q0HCfsjxIOxHDl2XFYbYTjMY57uxOFeqWoGNYlnNFTNWVk1zuo0UwITgJbUV/uGKrHZeO4WiCl5HBPJ7/bvZXdHW3kO13jKhBGXPL3H+xm2XWNFFXnplQkkqmwt1y40nSpsMPl6hnz2NJ6lM5+P0ve0sDStzW+4TF7O04lhMKF6vPVjxoEmnq060mjGRdiwB7gXuDbKJGoABpQfubxo7m3i2+9/CRfeuYvNPd10egpoTQnb1xXET2tAbY/foJDWzpTfuxkKuzahvHpg5QOBrugYonWdiF/lI5mNdOi0OlmY+uRQc9I799Oryg0mrTSi2oJ/Rzqiq8QJQ7jn6bZ4u3hT3u3s6nlKG67fdxdTIMpqcvjjvvXkVuU2nqQbEiFHS5nuqAe+NTLANz2g1UUOF0c7e3EFw6NS88tLRQaTcoxgGPAS6ggowWV3jrUlN/0c9Lfx5/3vcqLxw/jsFpp8JRgydBJ9Mi2TnpPBVl0VV3KRSKZCvvZVVeaPhV2uCRdUF39fq740GzcBfZBAig42N3Ooqr6tNuhhUKjSRlB1LCgZ1AjU9yMd3B6MB0BH48c2MEzR/djs1ipLyzCIjLrbd76cDPtR/3Mu6wamyO1gddkKuwFZUNOS85Kki6oLzz9F2oXFGG3DrxnDouV7SePa6HQaLKDLmAzagWRDE43ZMSSmBFnf9cpnj66ny2tzVgQ1BYUjcs86+Fw3V2LCHojKReJjoAv61Jhh8tgF1Stq4gNP9lHw/ximi4uZWvbMW5auBxLmheIWig0mlEhgWbgRdQqwooKTmemK2mbr49NLUd44vBe/JEQLpudmnyPKQQiHjV49mcHWPGuqThzbOQVp9anHopFCcViWZkKO1ySLqje/n4Ob+3AmWtj1qpKTgUitHh7qUtzQ9+J+a5qNGkjguqtswGVkphLptxL/dEIr508wd8P7+FQTztWYaE0J4+SnPGr5B4Ox3Z08/wvD1I9s5CZKytTeuyJkAo7HAa7oG7+7nLcbpVOLKVkX9dJ6grTkxabRAuFRjMsvMArDGQvFTPetQ+gWmwc6u7guWMHeen4IWKGQYHTRUNh5jKYzkbvyX48lTk0XVjKHQ+so6g69QLW4utlSU0DaxrM2RU2lSRdUA/t206ju5Tek/3k5zvZeOIol0+Zl9bX1kKh0ZyTVmAjqkBOkKnspa5+P5tajvL3w3voCQWwW6xU5BVgs6SnEnesvPLX4zzyPzu47furqJhakBaR8IaDuG0Obl64IuNB+vEi6YLau6ONBz+5jWvunI/vojCBSITcNDb31UKh0byBOHAIeBrVg8kJ1JCu9ghnIxyLsqujjScO72FPRxsCQUlOLg2F4zM8aDRIQyIsgpkrKuht66ekLj1usLhh0Nnv57NZ0BU2lSRdUP/ufZjVN01n2uIyuglwuKeXeWlsDaaFQqM5TT9q/sPTKFdTAZkojjve181zzQd49thBQrEoeQ4ndYXFGat9GC7PPLCf1r29vOsrS8kpdHDJLTPT9lrHvT1cPX1ipcIOl0ZPKdfOWshDYjv5Hhc+X4jtJ09qodBo0oeBci/tQPX1j6OarI1/emswGuFP+7bz+MHd2C0WynLzsyqLx51vJ7fISTxqpDz9dTAdAR91BUVcN3tR2l7D7CRdUAf3t7Pr4TYctwV573xJuq4lMv4pFEJcASxN2HKvlPINjdWH8xiNZvjEgOOozq3bUSsJGyr2kJn01t0drfxo2wv0BAPUFnhMG3sYjBGXvPTbQ1ROL2TqkjKWXtuY9oD6ZEiFHQ5JF9Rtz/yM3RtaqV3jpG99Hx5XQVpeb9zeaSGEB/hH4DLgkJTydiHENODrwGJUj9wHgCvPeN55H6PRnJ8wqu5hR+IWRYlCKeM9PW4wvnCIB3dv5Zmj+yh259HgMW/84UziMYPtj5+gqT3I1CVlaRcJKSUtvl5uXbSSmnGYuGd2Gj2lfORdaymdkYc7P3p6XkU6GDehkFL2Av8hhJgBJOPz1wOtUsq4EOIIsF4IUS6lbB/01OE8RqMZggAqGL0d2IdyK7lQwpCZlUMSKSWvtB3jp9tfoj8aTvRfMn/mjpSS3U+3MWt1JXanlVu+vQJXfvrfSyklx/q6WVo9OVJhh8s1M+ezta2Zw73H0/o6mf5kTkM1yAHoGbRtpI9BCHG7EGKLEGJLR0dHyg3VZAt9KGH4CfBV4JeolUQVKu6QuerpJD3BAN/b/AzfevlJnDZrIlCd6a/i8Di2o5vffWkbr/29BQB3gWNcVhLNfV3Mr6jl9sVrsua9Gg+UC2oRZTnp7SB73hWFEOIy4DUpZTrOvoMdsZEhtg33MUgp7wPuA1iyZEn61mAaE9KJSmfdipozLVADXGrJ/LXQAIY0ePnEYR549WXihkFTUanpiuTORl97kMJyNw3zS3jP15Yxdcn4uOsMKWnu7eLi2iZuvXDVpI5LvJ4Q6nMfo9FTzZ0rvozHVZe2VxvOu/53oF0IYaB6Ju8AXkv83CWlDI/h9TuAZE2/c9C2kT5GM2mIof787ShxOISa+SAAD1BPJmY9nI+OgI/7X32J106doDKvkBx7GqujUsxzvzjIC786yD/+dC0FZW6mLRufgsO4YdDc18XahhnctHB5VgT400sMJQ5hVOuYtcA8oNwUvZ4+BtwK/BbVAW0mKrB8MzCbgZP4aHgSuFKoy6p61OXgfiHEQuAdUsq7zvaYMbymJmuQKBFoR2UpHUT9+WXi5katHDLTqXU4xA2Dp4/u49c7N2MRFpo82bGKCPfHiIbi5BU7mXtpNVJK8opTOz/iXMSMOM193Vw55QLePW+ZKZobZobkd8CLOl0vABahToXjJ5xiOJFyIYQb+DhwA2qO4/1yhCF2IYQF+DxwbWLTn4C7ga+h1lFTUamvTwkhrgc+K6VcmnjuN858zLlea8mSJXLLli0jMU9jCoIoUWhDrRSOoP7soL4U+aixodlx0mjx9vDT7S+yv+sU1fkeXLbMxkaGSzxq8D83PMnMlRW85VPzx/31o/E4x73dvHXmAv5h9qJJGpMIoNrXS2AKcDHq9JfeWIQQYquUcskbto/kfC+EKATuRKWn3iGl3Jg6E1OHFopswGDAhXQEJQydDLiN8lDCkB0n18FE43EeP7SLP+x5BZfVRmyhaB4AACAASURBVFluvulXEQc2tnN8ZzeX3joLgG2PHKNiSj41s8c3DTUSj3G8r4cb5izm6hnzTf++pZYI6jsRQzWdXA5cgHKpjg9nE4rhBLPXALMSt9moqiQfajqLRjMCvMAJVKHbHpSvFVToKR+zxhdGwpGeTn78yguc8HZTne8xdfA10Bsmp1BlLZ3Y3cOOJ1pYdeM0HG4bF16d/qlpZxKORTnh6+WmBcu5bMqsSSIScdQ0xH7U9+BiYD5QjZm+C+ddUSSC2NuBXwO/lVIeHQe7xoReUZiFGKo9xmHgVdTqAVRsoYhsXC0MhZSSo71dbDiyl+eOHSTP4aQ0Jy/TZp2TI9s6+flnNnLTN5dTP6+YaCiO1W7BYs3MySkYjXDS7+XWC1eypmFGRmwYH+KoCyY/yq1kQV2DL0W1rc/sd2LUKwrgI8Bc4GrgX4QQnQyUt+6UUv4ppZZqspxu4BgqQe4gSiwsmDkjabSEY1G2nzzOowd20tzXhcNiM237jXjMYPtjx/FU5DB1aRm1FxRx8fVNFJQpn7fdlTmbA5EwHf1+PrpsHctqmjJmR3pICoMX9dm3oGIOa1Hp25mv6xkO5xUKKeW9g38XQtSicrLmo6qmtVBMasIod9JBVNZ0MlU1F/UlMN9Jc6y0B7w813yQvx/eQzgWpdDlNuXgIIBIMIbDbcNiEbzwq0M0zC9h6tIy7C4rV3z4gkybhy8coicU4JMXX8aCyvTVAYwfMQZWDKA+/1OAdQwIg3ndkWdjWBYLIWYBb0M15QeVo/hnKeXX0mWYxozEUJlJfpQ47EAFog3UVVERZk5VHQtxw2BPZxt/O7SbHadasAhBeW4BTpt5v/RP3LeH3c+0cccDl2CxCm75zkpyi8xTv9EXCuKLhLhzxZXMLqvKtDmjJCkMPtQFkhXVOGIW6nRZTjYKw5kMJ5j9GeDdqBjFpsTmWuBXQohfSyn/M432adJKHHXiDyV+Ju/7Ua0wvINuflRWhkjcJCoAXUO2pKuOBm84yKaWIzyyfyc9oQC5dvPOhvB1htjy52ZWvGsqzhwb9fOKsTutxGMGFqt1XOsgzkdPMEAoFuOzq97EtOLxnxg4dryo9FUXKm11sDBMvFX0cKTuVmCOlDI6eKMQ4pvALkALRUoIM3DCHvzTizpp96NOzmdyrmSEofaFUFc/vsT9M0/yRmKbfdDNwUT9AgzF4OD0iycOYxgGpTl5NHpKM23aG4hF4sSjBs5cO70n+3nuFweonVPE9IvKmbG8ghnL0zjNZpR09vuRUvK51VdlVbfcAdpR360PAk1M5AulJMMRCgOVq9V8xvaqxL5JikRdkRuJW/yMn8n7cV4vAj5ef7XuT9ziDFytJ48vUX8iG6MLeA111WtBnfhLmAhL4lQSikXZ3nacRw/u5FgiOF1l4rnUkWCMb7/nKS68pp5Lb5lJ7ZwiPvHryygoM+9o0PaAF7vFxqdXrc/CVuES1SGgErgRVeswORjOmeITwJNCiAOodwlU+so0VHsPE9IFfCcFx5EoH2TyZqDmGCTFYPCJfTBDbRt84h98te5Gjdyc+FclZuWkv48Xjh3KiuD0rqda6WnrP13vcNH1TdTNUSdcIYSpReKkv49cu5NPr7ySirz0DNhJHzFUNt9CVHOJ9FZIm43hZD09lpghsYzXB7M3Synj6TRu9IRQf9ixIhg4sVsYEIbkfU22EonH2NneyuOHdrGv8xRWEwenu477KalTdRlHt3dxfFcPK26YisUqWP2eN3TcNx1SStp8fRTn5HLnivWUmLzG5I0EUW1lrkSltU6+i7phfSuklAZqoPDrEEJ8QEr505RblRJyM22AxoS0+fp48fghnjyyl2A0Qp7DSUNhsSlXD6BaaTz836/x0fvXUVqfxxUfno3dZTWtvWcSjkVp8fUxo6ScO5ZdgseVk2mTRkgPqu/S+4A5GbYlc4z18umLgEmFQqNRhGPR06uHA13tWISF8tx8ynPzM23aGwj0hnni3r3MX19D06JSZiwv500fm3M6Y8nhNt+KZyiklJz0ezGkwU0LLmZd48ws7ADbCuSgao6zNX03NQwnPfa1s+1CVY9oNKak1dfLC8cO8uSRfYRjUfKdLupNuHoIeiP4ukKUNxXgzLFx5JVOai/w0LSolLxiFxf9Q3ZVK/dHI5z097Ggopb3L1hOmQkF+dwYqHBsE/AuVIPKyc1wLk8qUM65njO2C9R8Co3GNIRjUV471cLjh3ZxsLsDa2L1YMbYQ5Kf3bkRi0Xwwe+vwuaw8vFfXJqxnktjwZCSVl8vdouVjyxZy0W1TVnYIjyCKiZdAVxFNrTXGA+G8+35C5Anpdx+5g4hxNMpt0ijGSEycYJ64fghNhzZRygWpcDpMm3sYc+zbWz641He998XYbFaWP+R2bjyBiqms1EkfOEQHf0+ltdO5d3zlmZhLAJUmnoncB0qdyf7/g7pYjhZT7eeY9+NqTVHoxk+wWiEHada+Nvh3aZePUhDcnhrJzWzPbjy7AiLwIhL/N1hCsrcNC40XyHfcIkbBi2+XvIcTj61/ArmV9SaUpzPTwcq7f02lMtJMxhzfaM0mvMQNwwOdLfzwrGDvHziCDEjTr5JVw9SSoQQnDrs4+ef3sib/2kuS69tZObKCmatGssEYXPQE+ynJ9TP5VNmcf3sC8l1mKdFyPCRqGz/UuA96DE7QzNsoRBC/PMQm/uArUO5pTSaVCGlpM3fx+aWo2w4shdvOITDaqMiL9+UVdPSkDz4ha0U1+Zy+e2zqZxWwI1fXUrThWrlYDZBGykxI84Jbw/luQX825o3M6MkW3Na4qiGE/OBf2CyFdGNhJGsKJYkbg8nfr8G1Vf6w0KIB6WU/5Vq4zSTG284yPaTx3ny8F6a+7qxCkFpTh5FbvPVyHS3BGjd18vcS2sQFkFukRN3/kAgdPrF2XoyfT2d/X4CkTBvnbmAq6fPw5klc8DfSLKI7nLgUiZjEd1IGIlQ1AIXSin9AEKIfwceAdYAWwEtFJoxE4nH2Nt5kmeO7ueVk8eRUuJx5ZjStdTfFyGnUAWhNz/UzJaHjjL9onKcuXau/uS8DFuXWiLxGC3eXho8Jdy5Yn2WNvMDlfqaHD36XtRMNs35GIlQlDMw5BhU06MKKWVQCBE+y3M0mvOS7Nb68onDPNN8gHAsittupybfY9oirf0vneLX/7aFD923moqpBax45xSWv6MJZ262XmEPjZSSU34vUcPg3fOWcvmU2aZ0952bOAMV1gI1cvTNqF6nmuEwEqH4BbBRCPEQ6t2+BvilECIX2J0O4zQTm85+P9taj/HE4T10BH3YhIWy3HwcVvPlWAS9ER6/ZzcXrKtixvIKai8oYsUNU3Al3Ev5pRPDvy2lJBCN4A0HicZVK7cZJRV8YNEKKvMKM2zdSIiiVg4hVHv8GaiGfo3oArqRM+xvpJTyy0KIvwIrE5s+LKXckrj/nrEYIYRYg3IU5gLflFK2nbHfiZqSEwP8UspUdPzTjDOReIyjvV3s6zzJptajnOhTNZwlObk0FJrPldF2oI9wIErjwlKcuTZO7OmhNtGpNafQweW3z86whWNDSkk4HsMXDtEfjSCEQCKpzitkTcMMZpSUU5PvoSq/MEsK50KoztExVCv9eSjXUj06UD02RnrpFkU5+WTi/pgRQkwDHke5ttYBP0atCwfjR9nag+rx+2wqXluTXgxp0Orr42B3O1tam9nbeRJDSgTgceWYsp1GOBA97T565H92IA3JbT9YjcVq4aP3rzOdvSMhZsTxhcP4wkEQAikl+U4Xc8qrmV1aSX1hCdX5hbjt5hmXen4CqJWDgVopLEdNm6tFZ/+njpGkx/4Tqhrl9yjX08+FEPdJKcc6+OEqwCel9AkhjgBXCSHKpJQdgx7zM+AnwKtSSt8YX0+TRrqDAQ73dPDKyeNsbztBMBY5fUKqyis0bcwB4JkH9rPx90f41O+uwGq38NY755NfMnAlmk0ioVxIYbzhENG4gRBgs1iYUlTG5VNm0VRUSk2+B48rJ6v+X+oaNTn1UaLqHi4HpqMGCpn385XNjERybwUuklIGAIQQXwNeYuwTgpIDmAfbMx1VKpmkGlgKrBZC3C+lbH3DQYS4HbgdoL4+24aiZC/BaITDPZ3s7mhlS+sx2gNeQOCy2Shy51BmNa8/uOu4n2ceOMD6j8wmr9hF48IShEUQjxlY7RbKm7Lnc2RIiT8SojcUJPl1qsr3sK6qnunFFdQUeKjILTC1UJ+fUyj3Ui0q2XIqqlBOk25GIhQClT6QJDm7c6z8DfhPIYQLWJTY1n/6RYWwAPuBbwE3AY8JIRZIKV83EFpKeR9wH8D8RRXyUHcHEomU8rQKJe+rn5LkEZL3Bz9eSonNYiXH7njdzW7NtoyP1JLsDHqwu53NLc0c6ukAJBZhocik7qTB+LpCSENSUOZGSjjwcjsLrqwlr9hFw/wSGuabL1YyFIY08IZDeMOh09vqC4tZ2zCD6SXlNBSWZGml9FAku7nWAe9kMo0gNQsjEYqforKe/ogSiGtR7qAxIaXcK4S4Dvhn1JoS4OCgh9wC1EkppRCiDRWhKkF17xqSNr+P/3juUZKLFYEAIUGqYN1phECphXrcYOURp/9VWyXqqs1usVLgdCVubgpdbjyuHIpcbnIdTnLsDty27BeXYDRCR7+f9oCXE95eDvd00NzbTV84iCXh3y50uakt8GRJoBOi4Tjffd9TzL+ilqs/OY/S+jw+9fvLsTnM//eJGwPCIIT6ZDYVlXL5lNlMLy6nvrA4y2ILwyWCEonlqNDlRPw/mp+RZD19M9EtNpn1dFMqWncIIRqAT6JiFTcCvwGmCSHeIaW8CxWtSgavK1Hrz+5z2wp1ham/6pBSYkhJ1IjjDYfoCgaIxuNEjThxwwABltPjUuXrxCXf6SLf4Uz8dCUExk2ew4XbZsdtt+O2OXDb7afFxmaxpP3qPBSL0hHw0R7w0err5WB3B8e83fQG+7EIC4aUWIQg1+Eg1+7A43KbesVwJi89eJj2Iz7e9ukF2J1W3vIv86maMZDmaVaRiBsGfeEgvnAIIQQWIZhWVMZV0+YytbiU+sLiLK6KHi5e1Ff9H9DdXDPLcAYX+WCIi221T0opx+rIPYlaHXwRNfviYygH5HrgLuBB4BtCiC+hmrK8PTGaddwRQmAVAqvFgmuYX9LB4hKIhukNBYkZSlxihoEhDSzCgkAkFjivF5g8h5M8p5M8uxO71YrdYsVhteGwqp92q1Xdt9hw2KzYLTZsFgs2ixWrRaifwoLNYsFqsRCIhGnx9XK4p5Pm3i66QwG14kqQa3eQ63Ca3oV0NmKROAc3dTBzZQVCCMKBKEFvBCMusVgFcy+tOf9BMoSUkp5QP95wCLvVyoziChZW1jK1uIzagiJT1pekj3bUaed2VO2DJpOIM1z9E4LSaTnyjh9dlmkzxkzcMBKion4aUp4WnuRN8vrfkwgGsnTEIBeckdD83IRbzGWzZ6UgnEmyU+srjx7jz19/jVvvWUntBUWnt5sZKSW9oX56Q0Gaikp5+wUXMqu0MgsroFNBMh5RjXIwFGXWnEmGEGKrlHLJmdsn0yVK1mFNrAImSkgyHfi7Q/zqrs0sf8cU5l5awwXrqimsyKFmlgcwf0prb6ifnmA/dYVF3HrhKuaWV2dNzCf1RFEisRTV+EF/8s2CFgpN1nHg5VPEogazV1eR63GSW+jEalcnV2eOjSmLzZ8y2RcK0h0MUJPv4aaLl7OgsnYSCwSomtoO4K2owLW5BX6yoYVCY3qklHjbQxRWuAF44deHiMcks1dXISyCG/9zWYYtHD6+cIjOfj8VeQV8bNklLKqqz/LahlTQgWq78UFUbYTGbGih0JieJ3+0l81/PMq//GE9dpeV6z63iLyi7HJL+CMhOgJ+ynLz+Mela1lc3TBJYxCDkShXUzl6upy50UKhMR0nD3p5/J5dvO0zC/BU5nDBmiqKKnNOeyMKy92ZNXAEBCJh2vt9FLty+dCS1SytbsrKuprUk4xHLESVZOmmfWZGC4Um4xhxg0ObOygod1MxpQBXng1fVwhvRwhPZQ7VMz1Uz/Rk2swR0R+N0B7wUuh0c+uiVVxc2zTJ0lvPRQBVDnUVsBrdn8n86E+uJiNIKQkHYrjy7MSjkt99aRvz19dy9Sfm4anMydpOrcFohFMBH/kOJzcvXMHy2imToDBuJHSiqq0/AMzMsC2a4aKFQpMRfvP/thDyR7n5Wyuwu6zc/L8rKG/MP70/m0QiFIvSnajSz7E7eO+8i1jVMG3YRZmTAwm0AB5U0Lo8s+ZoRoQWCs24sPuZNl559Bg3fnUZwiKYvbqKWDR+uiCuanr2TE8zpKQvFMQXCYKEPKeLVfXTWFhZx4ySCi0QbyAONANzgOuBnMyaoxkxWig0acHfHebVv51g8TX1CfeSQTgQI9AbIa/YyYIrazNt4ogIx6Kqt5dhIIBpxWW8pWYes0qrqM6aCXCZIAycQM2MuBQ1llSTbWih0KQMb0cQi9VCXrGT3pP9PHHvHkpqc5m1qpK5l1Uz73Lz9lk6E0NKvOEgfaEgQghy7U5W1U9jQUUt04rLJ1AL73TiR8Uk3gkszrAtmrGghUIzJqQhERZByB/lf9+9gZU3TuPSW2ZSM9vDx395KUVVys2QDTGHcCxKd7CfiBFHIJlaVM41M+Yxs6SSmixqp24OOlEpsLcBTRm2RTNWtFBoRs3vvrgVi1XwD/92Ia48O2+9cwF1c1UTNyHEaZEwM+FYjI5+H4aU5NgdLK+bwqLKOqYWl5Hn0Ln9o6MVyEeJRFmGbdGkAi0UmmGzc0MLhzZ38rbPLACgfGoBFuvASiFb4g5xw6Ar6Kc/GiXH7mD91Au4qKaJusIivWoYEwZwDJgCvAsw7xhczcjQQqE5K8d2dLP5T0d522cWYHNY8XaEOHmoj2g4jt1pZc17p2faxGEjpcQbDtEb6kcIwaKqOtY1zGRWaaWulE4JyUrrZcBbAJ35NZHQQqE5TXdLgBd+fYjV75mGpzKH/r4Ix3Z003cqSEldHsvfOYUVN2RX0zY1wc+PIQ1qC4q4dtZCFlXVUeDMnjYg5icItKEqrdegO79OPLRQTGKCvggv/vowM5aXUze3GCMu2bmhlVmrKvFU5jBzRcXpSXGQHQFpUK6lzn4/oXiUHJuTN0+fy0W1TdTke7Lm/5A99CVu7wXmZtgWTbrQQjGJkIbkpd8dprhapazaHFY2/fEIuUUO6uYWU1KXy6cfWo/Vpvz0wpI9J1UpJX3hIL2hIFYhWFzdwNqGGcwsrdBdWtNGB6ri+kNAXYZt0aQTLRQTDCklIX8Ud74DgMfv2YXdbePSW2YiLIJtfzlGw4ISZq2qxO60cuef1mNzqBOpEAKrLXvEAZKuJZW11OAp4e0XXMiCijrynTpjKb2cAEpRK4niDNuiSTdaKLKclr299J7sZ866agB++dlNhAIxbv3uSgBC/ihGfGCW9m0/WI0zZ+DPnhSJbMMfCdHZHyDX7uAtM+ezrKaJqrxC7VpKOwaqHcds4B2AjvVMBrRQmBxpSIK+KDmFaoWw/bHjHHi5nXd8QVW6bn24mX0vnjotFIveXEc0ZJx+/ts+s/B1xxssEtmGlJLeUJC+cD8l7jxuXbSSZTWNujvruBFBZTatBt6EPn1MHvRfOgPEYwYWi0BYBD2tAZpf62bOJdXYnVb2v3SKjX84wru+vBS7y8rzvzzIUz/dx78+ehV2p5VwIIa3M0Q8amC1W1h38wwuvXXW6eZ6F6ytzvR/L+VIKens9+OPhmkoLOH9Cy5mQWXtBI09GEA/4EjczEJyhsR1wEXozKbJhSmEQgixBtUxLBf4ppSy7Yz9VwBLUfbeK6U8NZ72SSmTdhCLxPF2hMgrduJw2wh6I5zY3UPN7CJyCh30nQqy57k25qyrJr/UReu+Xl789SEu/9BsPJU57NzQwu+//Ap3PLCOkro8ml/r5qGvvUr9vGKKa3KJRQzC/hiRYAy7y8rUpWU43LbT7qOLrm/iousHWiIUlE3cpb8hDU75fYTjUWaXVvHWmQuYVVo5wdxLQcCX+CkStwrAizo5W1Di4UIVsLkZ/0E/PQlb9AyJyUrGy1CFENOAx4FvAM8CPx5i/9eBrwEvAg+c75i9x0O8+vgJQKWAfu/mp9nxRAsAgZ4w333/U+zc0AqoRnbfee8G9jyrtKmnNcD/vvtJ9r1wElBjOb98+aPse0Fp06lDXr7z3qdofrULgI5mP7/8182cPNgHQHdrgMfv2U3XiQAA0XCck4e8hPxRACqmFrDu5hk4c5VGz1pVycd+fgmFFeqEf8HaKj74/VXkJmZCV8/0cNH1TVntMhopMSPOCW8Px/p6WFBZwxfWvYVPr7yS2WVVWS4SUaAb5b45jqpijqLab78T+Cjw78AdwOcSt9uAtwMLUNdJrYnnNaMCyj0ol9BYMVCdXv0JG08lbDyCymz6CFokJi9mOPtcBfiklD4hxBHgKiFEmZSyI7H/eqBVShlP7F8vhCiXUraf7YAWmzh9IrZYBKUNebjy7af3VU4rJKdQ/W61W6iZ5cFdoJb5NqeVhvkl5HjUiTqv2MmqG6dSXJMLQFF1Ltd+diEVUwsAdeL/4PdWUlqv2hXUzyvmMw9ficOtXr9hfgl3PHDJadvKGvJZe9PAgB5Xnh1XnvaxA0TiMU76+xAI1jbOZP3U2VTmZc+citdjoK7CfSgxEChXUhMwFahGDe85Vz+s3MStHliU2BZDiUM3cBIlGscZevXhTLx2FCUmkcTzYcB1JBM3a+I5hUAlasBQIapnUyNQMPK3QDNhEEm3SsYMEOLjwF1SygohxELgFWCllPLFxP4fAsVSyuuFEMVA1+D9Q1E6LUfe8aPLxsN8TQoIRiO0B3w4rDbeNG0O6xpnUOTOzbRZI8QAelHCkHQhVQPTUDUG5UAR6fPtB1BfjS6UcDQn7MlDneyTJ/5ClPi4z7jptukaEEJslVIuOXO7GVYUfwP+UwjhYuCyqX/Q/sERy8gQ2wAQQtwO3A6QW6av0LMBXzhEV9BPvsPFu+cuY2X91Cyc89CPaqktUSuFK1FX5KWMb7+joVYfGk1qyLhQSCn3CiGuA/4ZFcEDODjoIR2obx4MXPZ0cAZSyvuA+0CtKNJjrWasGNKgI+AnFItSmpPHbReuZmlNIw5rxj+KIyCG+giGUauEq4ALEvc1molHxr+dQogG4JOob9uNwG+AaUKId0gp7wKeBK4UKopZj5rQvj9T9mpGRygWpd3vw0CyqLKOK6bOZmZJJVZLxvMpholEuXK8qK/Nhagr91pMkBOi0aSVjAsFKiLXCXwRlRf4MVQLyvWo2MXfhBBXAnej1vbvk1IaZzuYxjxIKekOBvBFQuQ5XFw7eyEr6qZSmpNNcwqCqI9nHBWIvgb1MdQtQjSTh4wLhZQyjGoYM5jfJ27Jx3xqXI3SjIlIPEa730dMGswsreBNU+cwt7wmi+Y+xFDiEEZl+6xHpbDqnkaayUnGhUIzMUgOBuoJ9eO0Wrl86mzWNEynOt+TadOGiWSgZbYV5Va6EJWxpF1LmsmNFgrNmIgZcdr9PiJGjLqCYt4x50IWVdbjtpup/cTZiKOEwZf4vQEVKpuOdi1pNANoodCMCn8kTFe/H6vFwoq6qVzSOJNGT0kWVE6HUMVqUdRKYRrKtdSIdi1pNEOjhUIzbOKGQUe/73Rq6/sWXMTS6iaTz34wUJlKfahitzxUU7sZqIylbKvb0GjGHy0UmnMipcQXCdET7MdisbCsupG1jTOYUVKORZjVdx9BrRrCKHFoAC5BrRpK0Z1PNZqRoYVCMyTReJxTAS8xI051voebFi5ncVWDSVcPEhVn6E3czwEWArNQweiJ22FXoxkPtFBoTiOlpCfUjzccxGG1sbZhBqsbptFQaMbYg4FyJ3lRK4QaYAUwBdVXyayrHY0m+9BCoXnd3OmpxWW8d95FzK+sxWXKyXHJNtjJ3kpvRolDtjUR1GiyBy0UkxRDSrr6/QSiEXLtDq6eMY/ltVOpyjdjW+8QA9XRZcBbUG4lM9qq0Uw8tFBMMgKRMJ39fgDmV9Rw2ZTZzC6tMmHVdJSB6uh8YC0wF9XlRaPRjCdaKCY4oViUvlCQUCyKEAKPK4d3zV3K0ppGik0388FAuZUCqBbdi1BBaV0drdFkEi0UEwgpJYFohL5QkLhhAJJ8p5sLq+qZU1ZNg6eYqvxCk6W1SlRAuhcVlJ4FLEHFHbKhulujmfhoochi4oaBLxLCFw6pgZZSUpVfyGVNs5hZWkF9YTHF7lwTZiyBijt0oOIOdcBlqCK4bOosq9FMDrRQZBHReJy+cJD+qCokEwimFJVyedMsphaXU19YnAUT4iKozvJ24ApUV9aSjFqk0WjOjRYKkxEz4oRiUUKxGKFYlLhhYBEWJBKn1cas0krmV9TQ6CmlpsCTRZPhYiiBEKjVw0WowjiNRmN2suUsMyGQUhIzDEKxKOH4YCFQriFDSpw2G+W5+TR5yqjOL6QirwCPy02xO5fy3HyTxReGQxwlEHFgFbASlcWk0WiyhQkrFMFo5LTfXiIxpERKkMjENrXPGHRf7SOxTd1Xz1BItYHT/0rUBbIUINQvAkAIhJSQEIBkhMCQEpfdTnluPtNyy6jKL6QitwCPKydxc+Oy2U0aUxgpBnAK5WpahkpvzZbZFBqNZjATUiicNisSsAoLVqvAZrGq+xYLNosFq1A/bZbENmFRj7Gox6rt6qcFEMKCRYjTN6uwYLGoGIEQiWiBUCKR/GkRltP78hzOCSgEZ0MC7agRogtRzfjKMmqRRqMZGxNSKKry8vivK67PtBmTkE5Uc77ZqEB1VWbN0Wg0KWFCCoVmvOlB1UFMQY0/r8+sORqNJqVoodCMAS/Qherc+naUUExkt5pGMznRQqEZBQFUsVwJcBOqRXKyjAAACmlJREFUUC7bsrE0Gs1w0UKhGSZRlDhEUdlLN6Ca9JmtmaBGo0k1phAKIcRiYDVqFNlfpZTbz9jvRCXfxwC/lDI2/lZORuKoJn39qNnSS4H5KFeTXkFoNJOFjAuFEMIOfFxKeVNCEF5AdYUbjB9law9wLfDs+Fo5mZCowLQXJQZzgMWoedNmHGSk0WjSTcaFArVSeKcQ4jngftRZ6kx+BvwEeFVK6RtP4yYPgyfHNQJXAdPQbTY0Gk3GhUJK2S2EeAj4IfBF4JYhHlaN8nusFkLcL6VsPfMBQojbgdsB6usL0mjxRGLw5LhS4GpUDYSuoNZoNANk3NEshMhBpdF8HzW+7LsJF1RyvwXYD3wLaAMeE0OUNksp75NSLpFSLikr01fBZyfZnK8ZtYpYA3wM+ASwHC0SGo3mTDK+ogCuA9qllP8qhHgAeAyYB2xJ7L8FqJNSSiFEW2JfCepSWDMskpPj/KhhQAsZmByns5Y0Gs25MYNQ5KBSapBSviyEeBpwCCG+IqW8C7XaSAavK1Gd5rozYWh2IYG+xG3w5LgmEm+3RqPRDAszCMXPge8LIb6HCmQ/hGoStB64C3gQ+IYQ4kuo3My3SymNTBlrfgYHpRtQb+N0wGzzsTUaTbaQcaGQUgaBm4fY9fvE/hjwT+NpU/YxOChdhgpKzwKKMmmURqOZIGRcKDSjZXCldD6wDlXzUJFBmzQazUREC0VWEUc14QsyUCk9D6jFBAlsGo1mgqKFIivwoorSrShhWISKP+hKaY1Gk360UJgWiVo9+FHJXu9CBaVdmTRKo9FMQrRQmI44apRomP/f3t3GylGWYRz/X5a+UIrokYJSagOpIaUYNNQCUUSjVCQaIBgiia/R9INRI2iimIgflH4xMakR1IqkFAlIohaMwRdMpEpEqZVI+SA0VJTWlJM0IH2zeHr74Z6l2+Oe2T09Z2dmd69fsuHs7HT36s2c3pl5dp4np+++lJxSw+s8mFk93Cga4zDZII6Qk/BdTJ5JmJnVy42idgfIby/NJafTeAtwSq2JzMzauVHU5gVygPqVwJXkIPWJtSYyM+vEjaJSQd4Yt5+8+fxD5DiE/zeYWXP5X6hKtA9Qn0MOUC/DA9RmNgjcKPrmCHn/w7+L56vIAWrfOW1mg8WNYlbtJ8cdjpBnC8vIAeoVeIDazAaVG8WMHCYbw6Hi+WlkYzgbWIKn8zazYeBGMS0T5LeV9pED0ycB55EztS4hv8FkZjZc3ChKBfAi2RyCnHhvOTlL61JySm8PSJvZcBviRnGQHCuYKP7b+nli0nY49h/79p+DPFNYTa4M91o8EZ+ZjZohbRTzyEYwj7w8NJ+cTG9e8XPrcSJZgtZj7qTni8iVWs3MRteQNorFwA11hzAzGwpe7cbMzEq5UZiZWSk3CjMzK+VGYWZmpRoxmC3pAuAS8mtID0TEY5Nev4xcqOEE4HsRsaf6lGZmo6n2RiFpLvDZiPiopPnAw+QMeq3XlwPfIJd9eyewCXhPHVnNzEZREy49nQxcK+mT5B1wz096/Rpgd0RMADuBNZJOqzijmdnIqr1RRMRe4D7g+8DfybOHdsvJ26whZ+BrbTuGpLWStkraOj4+3qe0ZmajpwmXnhaS83N/B1gLfFvSeRHxn2KXOW27H+6wDYCI2ABsKN5zXNIz04xyKrn83CBw1v4YlKyDkhOctR/6mXNZp421NwrgauC5iLhR0ibgF+QC0luL18fJSZbg6LzdpacMEbF4uiEkbY2IVd33rJ+z9segZB2UnOCs/VBHztovPZGTKc0HiIhHgN8C8yTdXLz+G+AMSQJeD+wCnqwhp5nZSGpCo/ghMCbpVknryPGK1wFrACLiV2Sz+DrwReDDEXFkqjczM7PZVfulp4g4CHysw0s/btvn8xVE2VDBZ8wWZ+2PQck6KDnBWfuh8pyKiKo/08zMBkgTLj2ZmVmDuVGYmVmp2scoBoGkVcAV5Dqpd0XEjpojTamYEuVkch3XfRHxUs2ROmpKTXvJUVdNu2VrSg17yTIoxyU0q67dVFbXiBjZB3kgfJz8RtUm4K0d9hkDniK/xrsc2A7MrSHrq4Avk98A21Cy31PFQXMAuLaJORtU055y1FHTbtmaUsNes9R9XPZ6fDasrr38LlVS11G/9LQGWAfcDGwBHpR0+qR9LgNeiogD5BQjK4E3VxkSICKej4h1wD/Jxb+n8mvgUuDMiLi3knBteszZiJpOI0cdNe2WrSk17DVLrcdlSw/HZ2Pq2uPvUiV1HfVG8Sh5NnEIeBxYwNG7v1tenmsqIv4LvEiHuaYaZAy4EFgr6Zy6w0yhKTXtNUcdNe2WrSk17DXLIByX0Ky69qKSuo50o4iIvRFxS+Q53EeA9RHxj0m7TZ5X6nCHbY1Q3L2+B/gm8AjwkKST6k3VUVNq2jVHjTXtlq0pNaTD5x6TZYCOS2hWXUtVWdehHsyWdD7wmZJdvhoRuyRdR17r2yzpjIjY3bbPOLmgUst8usw11c+sXd7mXcDbI2JC0r+A04EVHJ03a8ZmKWdTarqzhxx9r+kUutWokhr2qFuWump4PJpU126qq2sdgzRNegAfAO4EPgHcC5xFzil1K3nGtYLs2nOBxeR6GWM15t0IbGx73p71YvKsCPK65UHgNQ3M2YiaTpWjCTWdIlvjathj1sYcl52Oz6bWtYesldV1pO/MlrSaXFGvdWY1QY5TrALuAt4YEQckXU82kEXALyPiRzVkfQVwE3BVsWkzOb7yclbyQPkauQDUSuDOiNjctJxNqWmR9/9ySLqIBtR0cjbgGRpYw25ZacBx2Zaz0/H5IHAHDatrt6xUWNeRbhRmZtbdSA9mm5lZd24UZmZWyo3CzMxKuVGYmVkpNwozMyvlRmFmZqXcKMzMrJQbhdk0SNo3gz97u6TnJG2ftP1ySX+TtEPSl9q2z5G0XtITkh6XdPZMspsdLzcKs+psBC5v3yBpDnAL8F7gXOA6SecWL98IPB0RK4FvAZ+qLqrZUW4UZsdB0g2SthePz7Vt/0pxdvB7SXdL+kLrtYjYAuyd9FargR0R8XREHAbuAa4sZgG9OiLWF/vtpNnTXdsQG+rZY836QdIF5MqIF5KrJP5R0kPk79M1wPnkpHLbgD93ebsl5MI0Lc8W7/tuYKmkx4rtY+Q8P2aVc6Mwm763AT+NiP0Akn4CXEKeod8XEYeAQ5J+NoPPeBNwU0R8t/iM24C/ziy22fHxpSezeu0ClrY9P7PY9mpyHWQknUAu2zuTxmN23NwozKbvd8BVkha2xhKKbQ8D75e0QNIi4H09vNejwBsknSVpHvBB4H7gSeCiYp/rgZ9HxM7Z/ouY9cKXnsymKSK2SdoI/KnYdFtE/AVA0v3kJaI95DrsL7T+nKS7gXcAp0p6llwN8AeSPk2u4TAHuD0inpC0G3hA0g7gD8DaSv5yZh14PQqzWSRpUUTsk7QQ2AKsjYhtdecymwmfUZjNrg3FfRALgDvcJGwY+IzCzMxKeTDbzMxKuVGYmVkpNwozMyvlRmFmZqXcKMzMrJQbhZmZlXKjMDOzUv8DKQ37WKDHoHEAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(theta_cuts[:-1], median_ary, color='k', ls='dotted')\n",
"plt.fill_between(theta_cuts[:-1], *np.transpose(twosigma_ary), alpha=0.4, color='yellow')\n",
"plt.fill_between(theta_cuts[:-1], *np.transpose(onesigma_ary), alpha=0.5, color='green')\n",
"\n",
"plt.xlim(theta_cuts[0], theta_cuts[-1])\n",
"\n",
"plt.xlabel(r\"$\\log10 \\theta$\")\n",
"plt.ylabel(r\"$\\log10 M$\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment