{ "cells": [ { "cell_type": "raw", "metadata": {}, "source": [ "---\n", "title: \"Join\"\n", "teaching: 3000\n", "exercises: 0\n", "questions:\n", "\n", "- \"How do we use `JOIN` to combine information from multiple tables?\"\n", "\n", "objectives:\n", "\n", "- \"Upload a table to the Gaia server.\"\n", "\n", "- \"Write ADQL queries involving `JOIN` operations.\"\n", "\n", "keypoints:\n", "\n", "- \"Use `JOIN` operations to combine data from multiple tables in a databased, using some kind of identifier to match up records from one table with records from another.\"\n", "\n", "- \"This is another example of a practice we saw in the previous notebook, moving the computation to the data.\"\n", "\n", "---\n", "\n", "{% include links.md %}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Joining Tables\n", "\n", "This is the fifth in a series of notebooks related to astronomy data.\n", "\n", "As a continuing example, we will replicate part of the analysis in a recent paper, \"[Off the beaten path: Gaia reveals GD-1 stars outside of the main stream](https://arxiv.org/abs/1805.00425)\" by Adrian M. Price-Whelan and Ana Bonaca.\n", "\n", "Picking up where we left off, the next step in the analysis is to select candidate stars based on photometry data.\n", "The following figure from the paper is a color-magnitude diagram for the stars selected based on proper motion:\n", "\n", "\n", "\n", "In red is a [stellar isochrone](https://en.wikipedia.org/wiki/Stellar_isochrone), showing where we expect the stars in GD-1 to fall based on the metallicity and age of their original globular cluster. \n", "\n", "By selecting stars in the shaded area, we can further distinguish the main sequence of GD-1 from younger background stars." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Outline\n", "\n", "Here are the steps in this notebook:\n", "\n", "1. We'll reload the candidate stars we identified in the previous notebook.\n", "\n", "2. Then we'll run a query on the Gaia server that uploads the table of candidates and uses a `JOIN` operation to select photometry data for the candidate stars.\n", "\n", "3. We'll write the results to a file for use in the next notebook.\n", "\n", "After completing this lesson, you should be able to\n", "\n", "* Upload a table to the Gaia server.\n", "\n", "* Write ADQL queries involving `JOIN` operations." ] }, { "cell_type": "markdown", "metadata": { "tags": [ "remove-cell" ] }, "source": [ "## Installing libraries\n", "\n", "If you are running this notebook on Colab, you can run the following cell to install the libraries we'll use.\n", "\n", "If you are running this notebook on your own computer, you might have to install these libraries yourself. See the instructions in the preface." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "# If we're running on Colab, install libraries\n", "\n", "import sys\n", "IN_COLAB = 'google.colab' in sys.modules\n", "\n", "if IN_COLAB:\n", " !pip install astroquery wget" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reloading the data\n", "\n", "The following cell downloads the data from the previous notebook." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import os\n", "from wget import download\n", "\n", "filename = 'gd1_candidates.hdf5'\n", "path = 'https://github.com/AllenDowney/AstronomicalData/raw/main/data/'\n", "\n", "if not os.path.exists(filename):\n", " print(download(path+filename))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we can read it back." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "\n", "candidate_df = pd.read_hdf(filename, 'candidate_df')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`candidate_df` is the Pandas DataFrame that contains results from the query in the previous notebook, which selects stars likely to be in GD-1 based on proper motion. It also includes position and proper motion transformed to the ICRS frame." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEJCAYAAACdePCvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACAiklEQVR4nO29f5Tl11Ef+Llm9GxZDwKztiPAyEKa2HInmhmsHqIR2ATkEIfjCKJMHIVA7JDEwCaEs4ZdBzkCspCTEC9hc2bXA4yDdxdDFEdRZK092N32TpvMaCS6R55pm9YIjYxHlmnFcnfA/bB4zxPu/vFevalXXVW37nuv1aPpb53zPe99v99769b9VVW3qr73ppwzGmiggQYa2Hnwou0moIEGGmigge2BRgA00EADDexQaARAAw000MAOhUYANNBAAw3sUGgEQAMNNNDADoVGADTQQAMN7FDYVgGQUvralNJ9KaVzKaXHUkoHt5OeBhpooIGdBLu2ufx/C+AjOedDKaUWgJduMz0NNNBAAzsG0nZ9CJZS+hoAZwHckINEvOxlL8vXX3/9ltLVQAMNNHClwenTp7+Yc365fL6dK4AbADwL4H0ppX0ATgP48ZzzH/NEKaW3A3g7AFx33XVYWlp63gltoIEGGnghQ0rpgvZ8O30AuwC8DsCRnPO3APhjAP9UJso5/2rOeTbnPPvyl28SYA000EADDYwJ2ykAngbwdM75kcH9fegLhAYaaKCBBp4H2DYBkHN+BsDnUkqvGTy6HcDKdtHTQAMNNLDTYLujgH4MwG8MIoA+A+DvbTM9DTTQQAM7BrZVAOSczwCY3U4aGmiggQZ2KjRfAjfQQAMN7FBoBEADDTTQwA6FRgA00MDzAL1eb7tJaKCBTdAIgAZ2DHQ6nang8Zi59q7X6+HMmTNu+ZMIiF6vN5J/GvVsBNbOgEYANLDtQMxGMrJSei2/lX59fR333XffJuZYy+gkM5e0nzlzZlM9Wq0WZmZmsLKyopbX6XSG+bTySgJncXERS0tL6HQ66HQ6I/Uch5HzekwCjRC5/KERADscLKbjvZ8Uv3xPDHVpaQmLi4tFhseZE88vn9PvqVOn8Nhjj+G7vuu70Gq1TFyRenBmTmXSLwDs378fADbhbbVaw3ccZ6/Xw9mzZzEzMzNCG71/6KGHsLS0ZNLYarVw4MAB7N27FysrK2i1Wjh06BBarZZav0h/aLTWgtW2lgBsYHugEQA7GEh7lAyCa7FSo60RCBqz5u+AS4yx3W5jdnYWBw4cUBkhAaWnNDw/PZfl7tq1CzfeeCPm5ubw8Y9/fFPZsjwJXLj0ej20223MzMwMy2y329izZw9ardbw4nh5ftm2nU4HKSWThlarhb17924SXACwvr4+ZNhEE5XPBRLlLQk8SV9JGFvAhQj1h1Z+r9fbtGLxcEbKbaAOGgFwGYGmvW41XLx4ceReMq+ZmRkAfY12fX29yEAsXBbzo3T0qzF/yTQ4Q6N89JyecWa4b98+XHvttXjzm9+Mz3/+85sEigVE68rKykg7dDodnD17dljm+vo6HnzwwREmxmlcWVnBnj17hgKD8s3MzOD8+fObGDyv5+zsLNrt9iYm+swzz+AXf/EX8du//dsjdErmy+koCTyu+Z89exYXL14MafDyPfWX7HeOn4QirVja7fYIfkspobxeuSX6GkHBIOf8grluueWW/EKGbrfrvnvkkUfyxsbGyO/a2trUy+JpTp48qabtdrt5bm4uLyws5G63O0KXhevEiRO52+2aZfPna2tr+ZFHHgnRSWVSG1EZ/P/JkyfzxsZGPnHiRJ6fnx/e83QcVwkkfnpGOKndKN3q6qpZF952Gxsb+X3ve9+wX6225zRr99SGvD4WLpnfqq/MQ3TTvdYmVnke7pMnT+bjx48P21JrJ63tKQ1vv1LZ8v3JkyeH41TSdiUDgKWs8NRmBfA8AddQuAajmUK4WUFqljVllfKRhqlphJ1OB48//viQPk6XBSklt2yujZ4/f35oNinVxdJsSSsH+ooMAOzbtw9XX301brjhhqGdnmu9Hv2WqanVaqHT6WBxcRHLy8sAMGw3ouPxxx8fsdVzXGSeIRv9HXfcgXPnzqmrGVm2vOd52u32iGPZWj0BcDV/aeKiPNTnAIbbsFs+Do5Hgix37969SCkN8VK5y8vLI6Y1Xg6v85ve9CacP39e1eS9VQCNdzIzaibOElxpq4dGAGwBaIOEJnGv1xvaPDVmyQf77t27hw69GihFnci0ku5er4dz587hNa95DQ4ePDjCiDw8ZK7gTlKPPmsSy7TcIcnt6qdPn8bi4iKAPuMnBrtv3z489dRTuO666/Dggw+aE9wzM/CySAjt27dvhPFz2LVrF2666SbVB0F+Ay5QSGCVyqb3Gl5qx5JJR9IrfTpk4pL+Cz4Wc84j95p9n0xdFnDhsmtXfxear3zlK6p5T5bDcTz11FNqvSM+Hf6Oj62o+SiSjua2h+eyAW1ZcLlel4MJqLRk9MwqBHzZTqaQjY2NTctTwhc1lXi0lpbGcnlfWsqXIGLiieLTzDgnTpzIx44dM00gPK1mQtHatcasQcBNddxExsvi44Gb+iJtwU1rVJeSacjDSfgkTVo6jpvXTzPDURptDGs0UV0WFhbyiRMnhuM/Ujetntr7SL08nBb9pbl99OjRoelUy1/iD1sBMExA287Ua67tFgAlZkyDI2IL5xOIM43IgJ8m3ZJeizFajMuqWyQtpY/WgddldXU1Hz16NGzTX1tbG8uO7TEMXk9uM5dCSRPq0XoTY1xbWxvavjXGyBmzRbcUgCV6LKFlCVreFhZYyoaGtyQ06R0JEdknXr954PVZZGzTWCvhjoyBaQmKRgBMCJ5WQc9LjJTS8QFrObOiWn9J0yrdEz2WwOED3ppUGt0R2kv4LOBMx3rHcUqtrMSguBZrMUvLOSqFgJU/4ogkvDRGNAbM00nBo/WjRotXX69dtfLle06T5tzVQEvnMXRaQWhCUaMzIvz5nOh2u/n48eNDzd0bt1YZ3hiyoIYPlKARAAKiS3B6F+ksubTT0pNGxyNCtGiIEgP36CoNUDlIreWqhs/S2CT+ElB7TTPKiTNMWX9q7yhtlFYyaovRybI0LZTnl1FEFj4yd5XoLjF7Lw2/SBh4iok2d3hd6D0f31YUmRxLso2l6c/qd2onXnZp5aYJDCpb9jvNWY7DAkl7tD80mAbzz7kRACPABybdR+z2EY3NA42ByMEiaeP5NEZvaRaWdmxN+iiD4RqXNaEtWmW6+fn5oiYcEXwW85V5SvZpiV8bF7UM2Zr03uSX/eT1z6RCWMtH5iZtDFmhoRbj5eNb4pP+E00YSEHIBSfHQ+G/3tjU/EWcNr7akv2jhSJbbVpaAXh5tgoaASCAd+iJEydcLZgzj0nNFpYGwvFZ2pekn09SKVy8eH2ePsoseH5tKWzVxaJV06o0E4W3kiHcJWeuRluEgfNLq38JIkqDV76kRUsjFQqrHSNl8Hp5Nm7Z95oDV2PmGu3SeS7fSyEjBYr8b5WnKVayHayxJPmD1ralZyXYSuafcyMAXIgwQamNRHDy/5pJgr/jmkdUI5eCgmsxkuGWaKx5rzEnz/FGtFF6S/B4jN4qu+SQsya1ZChauVadosz/yJEjww+ePBxevWX7yvQUeXP8+HFVaFlCVFNCZLSRBdqY1UxlpUiw0nySgo3THPUn8LJq5y2nQzPvlsaINo+jvGba0AiAPFnDliaUTOs55eTg4AybJlNJU7FWLTyNZkeld1oeq5588pSYWHQy1qSxaCrh8gTExsZGPn78+JBJl4RWLXS73SHz1+qj2YelJsp/Nb8CZ7T8S2jN92DRKO83NvwwRp7OK8Myw2i4tLmysbExEiyh0Rztl+ic9XB6/aflk0Ix5zysU6ltt8IctOMFQKlhPcmsaThaGg5S+5ETVxsc/L+ltRJemvSeALAczJrpSDIjSkuM7OjRoyN10jTISLt4YDF6jaaoGaY06UsrgBqGY/UFgWaDpufdbt8ZLxkn1z618UfP6dJs5JKmUj28FYDG/CyI9j03A9HcmJ+fzwsLC25IZW0ZVl4uhCKCM1qu7KeTJ0/mY8eO5fn5+Yn6ZxzY8QIg583Lav58YWHB3RfHese1M46fa/GWlhOhVZYlTQGSkUjmoeGRk0FjcnKic5xaWq39tHp79bUEijZRvVjqSHklkIynVBdLwSgJZC7Y6bsGqQxoTmern0pMa2FhwVzxRKFkIqrFRdFg3LxDkUjyozUOUSGt2f+t8cbva3wiJaB6HT9+3P1IbKucwY0AyKPao9TordA3Amtg0XJZDi4vJHSSwaQxD64pr66uuuFnNU5sjelYTM4yW8lVQ7Q8i06p8Wr1s5x9Xt0kzffcc88mIaAJYIs2urciYrSxQZu7aSYeWY7GvCJmCY8BRYCYsvZNhcXULKB5wv0XMl9EqBFdHgPVVhIct9a3/FudSFSURR/nB5H6NAJgCwSA1MDkYKthjJSH8GoTzmP+mjZCv7UaAM/LncAWHRZDiYKnAWvPZLRGhEnwfFYdPM0s4vwuafMWQ9Fs8TxfdJsF7RnlLZkDOQ7J8HndNHPGJAyGGOGFCxc2ObctwVYaz9HVhEU39ZUmsHk+HsXHy+ZtxZUnnlebM7x+pVWhtYrTYNyggxLseAGQs20HrHUs1gx2i4lY78cxX3AcWsQCB3qv+QBqyqmhTftfwunRZwlMbeJ67RDpL41BlZh7tC01fFGBQc+5Bm3RN06/aWU98sgj+cKFC/nkyZNTWwF4SgEfA9Y+WTTWV1dXXfwbGxvD1YsUlDn7349Y/aIpkzKfVEa8NtJomEbf5dwIAJfRS+YcsTlqjNqbqCXaqOxJzDN0zc3Nmf6GbnezGaCkoZXqMi4jtMqnia35Hawy+ISe1EnJcZa+EYng1+pnMaJaKAl8ScskIJ2144InrOj//Pz8yF5PWn9SO2ofksn/3LFs9Ys21q3VHq3ktVBXGr/S36L1tVyF1PCVGrAEwI7ZDprvH8/PoJV7si8vL7v7idP2sdpWy3x7XMJHz+UzufWv3JbXA2tbWqrXk08+ObLPutxCeHZ2drgdr6SPAx3X98wzzwzLLZ3YRc9OnTqFU6dOjdTX2n7Zap/l5WV0Oh21rfnWwlQHfp6CV69SvfmWxPv27RtpK06jfKbVzzrwnfbDt7bujm4ZTEdpam2o1alUD+0dnWFx/vz54X79fM6UttqWz+V41MbQ1VdfjTvvvHN4foO2BTfQ3wac6CJci4uLI3Ob+vGlL30pbrnllk1bTlMdtHEmjxslaLfbOHToEHbv3r3pHc0xXhY9533N57xWBs+3ZaBJhcv1mtYKgIebWZo0v7fwlcqwzEQlB2e0PhJIe9AiWCL5NVhdXc3ve9/7NjmXPVzd7mi0CbV1aa8kuUrzlt9eHcZdNkfK0XBr9SutNLVxpq0EI2OR57dWFBouz/zCy9fs51a9S7itemi4vag5GQnmta1XrqxrtA4WeH2u0TaNMiOAy9UEBOCrAHwSwIdKaacRBqoxGS/9ODbvyCS0BqN0AtaU6wkXmVb7b9FFzF9z5nqTXC6ZI+1oTYhpMXPrvVcvr4040+Z5S+c6yDFCebgjX9bfMxFwJq05gOWOn/TcqluEoXrvaueMF83V7W7enVNGIlkmoGj5tc89YSIFuDUuSsJhmnA5C4B3APjN50MAENQ4QCNCIvI88ll8t2uH6UXKKQkX+rUiN6yJpDmm6Fc6lC0aI0LDyjttzciix2sji/FaETdaJJCkXzIKuceSFJqecO92u27ILWdGWv9b7Vxi7pG54dnxPTzcB0QfIsr6yPzEZL2PJD0GXAISRp6fkNMnN7CTu8BGeNG4/iEOl6UAAPBKAB8H8F3PlwCIaAtyYHmTLuKAI21Gi1SQZciQMY8BRh2e2kSn55aTVw5WCy+Zekq0jBt1xNtgmkB1LwkDor0Uwuoxm4g5RNbRKtMSRCVBZd1rgkB7ztNrwtGinx9kEzUDam0VOWFO4tDqYoXxWv2kCRrr/AGrba33UmBb8ybyXUsJLlcBcB+AWwD8JUsAAHg7gCUAS9ddd91EjUDgMXhtIHhMMhrfG90YS+K37JNy8pWiM7SJru1yWKqbxpD4Urzk+9BwaPjlhK0VAqX00sSg0WnRpeEu1akkPCzmI+9LDNnCXWIgFuP36M55dF8qL8RSawdt7pW+/NXotepDq9NI9J91b/l2rPFpra4lXm3+WnW5IlcAAN4M4D2D/6YA4NdW7gZq3UtGKDurpnPGYRJWaJpMTxPR+ghJeyYdaaWTobSQQ8qnTXwOEQ2u9PFS1IzU7Y6e4OSVZ9Fgxf8Ts/POr5X0RRi3x8i8e3ommawco5ExWxIm8p5++XGbHp1e/WRfl/JI5i5Bju0oHRpdVj9rv57jWuIotVe0DSNwOQqAfwngaQCfBfAMgC8DeL+XZzuOhOQMjj/LuS5u38JNv5rt32IcHrOQg9XTMLRBaGkvpOlr+8hEBnJkhVLSVD2GK+tgOdMjfUWCzosSK7VByZbutTVPEwkK4DHplN4qr2RS0PJx8yTvJ27ikWXWCACNsZba12PupbHk1ZHuS4qUhjP6zUhUIEWDVSJw2QmAESK2cQUwjnZAz6z9b6LlcqYot4nVzCWW9mDRywWUt5LwcNB/om+S+nqCSJZp0SLxldqnxJit8ixhFR0vpf6JRK2QYhD1M1E/WYLLGguR+nDBIU/NkkpGZG7I8SCZr7fq5mV7+K3nGi4tiqsm6o3aIXLWc4RWPr49s1gUGgGggJT0tXm9Xf34rwWSQfMJ5Z1e5Gk4FnOJ2Ji9ukrHF8etpfdwaXXxhIL13toegEAzKZVos/BqtGjl1ggZ+pWMQ76vGZ9e+mi/WGlpszpabUjgNn9LEBFuzy6v1UFLp9EaFdD8P61s5e6fGo2lvpVCJDIOPPMg3U/C/HO+zAVA9JqWAIh2UqmjNS2H25Cjn+hrg0DaEiUeL0rEm0AafosuDa/U8mq2gLYGubU6sZyHvI5e31l7HkWYUonByPRRLVGjUduaY5Jzh2shEhHD01jbH8goG48xe3m9Pi05Xmu37eh2Lx3cI2nSaK5xpk9yZsW0+pagEQADsCarFrFQMrdYE4cGU2Tf9VJHc0ZrTSpeB21/Eo2Za1ocZ35yKauFhFptY0UCaQesW2fCapN8Y6P/kZM8scqaQJpW7WlxJcFgpY8IXG8Mae0a2VGUbPA14ZFa+VqElsaguenHY5C15kZZbimdNQ75mQoRoHyRL5kpnWfamkQRsP5HowxLsOMFgKeNEEOkbQ+siAptwFsT3WKaMm0pLpro06J8KD0fJJKhWJNPExJc0Fga6TiDmuPgbbuwsJCPHj26aRteDYdlLqtl6N7zSF0izE1jBFb9IsKlFFnlnWPB8WnnVkjm70XVcIGjCXkSzvwoTFkfWTeNTlnfGua8uroaGqel+WnN9cgqpVSuvLdWMxp/GBd2tACIaCA04Kwto7VJ7XUO70Try17JvEvMxBI4NSGJMq9Fc0l79pi1BdqKydI8rUmpvZ/GBNFAE5qab8ZimDJNjZDi7zVnPj2zhKdGz5EjRzZt4yz/WxFEvExvtUDKiifYCIe1jTV9NKaZHa2xR7TTLpylfb6kQuKBnFM1Kwx5X1pJy3mmBYOMAztaAOQcYxiepqbhiDi5SCvyVhXWYCjRG00fCa30nkttLGIztvB5/UB4LBMWZx6TCCJLyGiCV2rHMgLGGxuRlUqpzTjDkCY+y3xmAdFjmdfIH+EpFFafeWNYG498Xkg6aOtmeYqW5py3HP0anYSXC9GSqY2eR1ZJkXbh+z15Y0Prm0mEwI4XAASlDowMcJnWe5ez71fgnSyFw7h1kWWPe5ydll8yRq3MSDkW7TTYrQmnaUo1k9PSnLVytIkq20Azh3S7l8JmS5t/lcaixvAkLTWg9Z/FPKNleP3q3WsavkYHx+2ZsEqwtra2aRUkv17WhJcWlqr1idUumjC02lrmpzpP6gtoBAADryFLzFx7F3H8eO+lfb+moz3mROCdhBaxH9cwe6scSwBp6b3QTg2H1f5yktJkJl+PdTSmxQystBqNfPJ6IPOTdmwxQVn3qKLg3Vt5rH6rEQg1915/cJq8scHbzPKbWMxbCgINp6RDMwtqabU5RIJF+65CKheReepBIwAC4E18776kgUYGrjWZrTxSe+B0a0t8S/DUOlBL77z0nqY4btkWk+QTj5tOSLh47V1ivNH6e3XU0pGGevjw4RETiAbc4eqNzRIN3nMtwiiyarTwlcqic315OZ6A1/BwpspNTPQ+6pwvAV+9yF1KvTaTfcM3yyNcFOUm26K2zTk0AmAAkckoO1E6v7hNThtUFk5+Xxslwhlat2tvUmU5maLtwQe25oCKCDMNf0SzKtFn9YvV5prztESP9RGYV38tjeyXaD3lCsBrD8thG+kfOXa1dqCoIdlmUW2U57O2SeA0cE1YmkU9/Fq9ut1L8f08KqlGgfDqo5mtCDS6NV7B03McmiN9EuafcyMAcs6xpT1Py/Nwh5s8iNvTICJahpxk9CsHWGn72G53c8RRaYBLOmg5a2lgWv20ENGIcPT6w5owHJf2AQ/VQzu0XLaHZd4oTboSU5KHlZSEShSi47dUDu9bKSh5Gs1hGWFMWt/RPkVWX2r5PWEbcYBT3Y4cOWJ+k2ONMd4+9MxyGEtFoxSJZLWX9nxa0AiAAUjHkwXeQC0xgBoNmS8B+cCTH3/Je2vASGeVNhB5OdrEl/9LjFDbdpebWmQemZ/TQ+V6jI5HkWgRQwsLC/nIkSOqfZfao7SsljRp76w+1qJLPAFqgTUGS4KxhEcqNJq27QlC/txz7st0JVOJ9lwbg6SoXLhwIWQaofFigTXGpF/O2g5FhqyWVi48L5/bkwj1EjQCINdpUNbSOFpO5DmVo523KzUzftXsmV7aFVOb5CXti9OuMQpLKMg8lpnNY7xylWK1qxUtIj9G8yDiIykx5FK7efi1dPzS6ucxZI12Xgf+pes4vjAvTWklmbP9EaNmgpUrAK1MTVB5c1+jRypEmulV3msreo22CE2S/nGhEQAD0BidlS46WUvl0L1nCooOyhLdkfxaGl5PKxopMvkJ+MT0nGDy16obp7Fkgy4xbCutls/THEsrFQuv1Qaa4OJlcSZdM46i40UyLjlmvbFhgRxbXtlkquGMnvJwJ6l1kA+/93xklslS4tB2YrX8LvSMVqbkP6lxCJfacRJoBIAAbZBY6bT/JdzjRhdo6T0aanBG8JCWExVgVjma89pKaznUiR7+36ODT3KPcXt0lpiCfFcqJ8IAaSzOz8+re9lIHJ6A055H+o5oKJnsrBVCRBnwaCaGazF6uQLwcPF6y/e8Db2VkTaOvHkt/R01Y6bb1X13Vt3GgUYAZHuQRKRvDfPXyqoFPli4XdHSyCSUmIjG2LQBWjsQS0JL6wPLLMPNPSVGymnnmluJZmuy0r1sm3EVgojPgdJZ7+Qz/t+KLqkJepDHOlrpOE3ctOb1i5xrJQZM9Gv1rQWPcZfo1d7z/5E5VCpTbmInaajd4VSDHS8AotJbyxc9kCGiadXQyyeCtWzUHLyScfJ38r8lDMapT2milRiVNmmk49ArWzKnEv1kdvDCEycRhh6uSB6eV2O8Hj0120TUtpvMaykiHK9lP98KKDHu2rRWGivSpyTsJR7pIOdl1u5wasGOFwA524MgsoyXXxDWlhGBWm1Taow8n7ViiNLMhU8ts5P31mqDg7fHfEmD9RhWSZu1TjmjNrSW7DUaJccXAa/NIvZjUgA04akJSq3dI3Wz0ku8XACOK0Q1sBQGTWBadEdMY1Zf87FT6hMPSvRN+hVwzo0AyDnrnRSR1uTYqdnffxzatPA3jRZKXzoRi99rzENjHDwthTOWYtgjDNqra7dr78zItUdPkHmTqPbQH17/qNPOYqRcIHunvEmcUXqt8etFQZFtO+IMpWeWX0ULebVMPqX2skBLo7UnF3zdbnf47YE1tiLhmrIe/BlvP8tKUDNWvfInhR0vAGQn8XtPwnKGXJqUNYOa5/Pw8DRy0Je0G+0ZZwZS4Ejc9OuFTZImVPOhE8dJ99rkoTanD6tKTmCrvnJzOQmamYcYGwkBTWh5eCQ+rxwrhLO2bqX0PJ/V5lo9ut3uplPL6Ln8KFKWVTKv8l/tv4WDGLxGO41J+k5Eti0XyNoHarwM2TfW3NTmohX6WrPXVUmBicKOFwA5168ArM620uRct0FXxCmoDdwo7pIjSr63tCSPQdEAjZjRePrSQTjSj1Hj0CR8VIZXBy3ahMogpie/Ki4diRhpA1ra0+Z09LxUNynEa5lDifFYgk5q+pKemjI4yH6Vu3Nq8ymyMpUCQrY9CXZ+9oCFp7ZOvBwppLhiQelKez5N41CYRgA4UJLAVifJwUHpSl57T4vUyqb/NZM9EnUStX+WlsolmrQtLCLtUwpJ1N5x/JrZSGPalqOZmIg8cc07FD2iTFDaEydO5GPHjo1sWhbts+gY4nmjqwzJLGvKsupgtQlnvpLZ8Xs5D0r+BK8tqR3IwRo5kU6j3QJeH9l+VK4l9LQ6TOoAzrkRABOB1knaZMx5NI7eGnxRjUMOnkg4WI2mEhVA42jfBNIuO67TzUvDJ0rtttqlcvh2vfQsouVH3hEuzsxq+zciTKkeJXORRl+UkZdotBgs7yst7l5+d6BpxCWhY9VrY8PfMyoqzDnI7aGpfy2aeB9GgiXGhUYAGFDDKLRO0gYDDS7tMA/S+qxdO72QzIgzNsI8IyC1laj9XUujTQAPamza8j7KFCO4Seh6H6pFQKOJ44to5hpj0E7wknk0hyzHo6XV7q00UUZp1Y3jKn3hzNvLwuPRwO/5KsBrEwu8+vJdTckXYZlXtW9gJJ4SLRFoBIAC0SgXDhZz1vDKg9WpU7XzTimNxiSscjz6rHc1S3puO6+lIcpktf/jMFmJu7T891ZgpT6ppUUyao0xeng5zdrJVOMyKlm+N/YkvghzL9WrNm2kLKvPtXLkymLSccdxc3o0PwYJh/n5eZX/WL6pcWFiAQDgGgBfFU2/Fdc09wLimkDtIC1pInzXTqts+UyzF06yjXBE44ponR4TsMop4ebvrUgJrYxJmAnvN82J6+Wtfa+lrxUk1jiNbAJYQxPXhrkmOu7YmAS8frdo4mm04ALLfKnhthSdaD0tRUajvdsd3c68JLgmbetqAQDgRQC+H8CHAXwBwOcGv78L4N0A/pyVd6uuae4GSvc1eSXj0r5qpIFYw1wsgTLuJItqDKUJHrHbesv/UtmcIVtlSyY1aXtYfV8zyb2oDG08WGlK9GqhqzUCN1InTwGxhPIk88ijQ/pvtBBMmccKLpA4tP9Wu1rjz6JbzgWtXbRxwRVFb8U6DRhHAHwCwD0A9gJ4EXu+G8DfAPCfAPyAlX8rrmmFgUYZqzbgNUavDUQLX4lhjjuZJmFoNTikgPMYasRWazEeb/KUaNRgHGGngebbkUyU/0YdlpLGiOPdY2Ae44r0jRzTcjXKNddahYc/94SMJgA1PJpA4ooagRZVZNGp9YWmoGnndESdxzIaSMLGxsamYy3HhXEEwFXWu5o007ymsRdQ1O5vMTqZhn5Lnc/TlWiMCidtwFu0lsrThJysp3evRSxph7Jo9eMT3fuoTBMM07JDe8xUlqfVSTOjWMxRu9fSetqh9lz2o2Ve8xyTGp1SKPM55MWpe+Y9b1VnzVFrTMhIm5yzGnXm9Yen5Vt1oeeab0ajU/6WVhYUpTQN/8RUncAA2uPkEzi+CcBxAI8NzEo/Xsozrb2ALIkqB1tk46oo840w2CgjI9osTbxmsEhB4mlTJTwceIidl062i/fZvnaurMXoNPwcl1ePku3WqntEaFnlaSZFT7hFhaHli+KRKh5Y44Hj5SGPWr0sGqVgIRyab87qC1I05DkCpHxo8yMqrLS24PXU+k4bC54Qt/7LsN1JmH/O0xcAT42TT+D4egCvG/z/agC/B2DGyzMtJ7AWb+xNdktzkOkiZcsJJLUhq8O1iail0RiDR2+JyZUYWKm+tc88Js4/ytIYUoQO+i7BW1ZPq/6cvqhg5xCNVJJ4+L08AIWPOW3jtAh9msIg498tfFq/1djbtWekhROetbW1kdUnn7+eYPXqz5m/tWWKVScPN28LGfapzedxYRwT0DuM6ycArFv5xr0AfBDAX/bSTEsAzM3NmSF5lnYwzufYVnrSZjkNNJC9pXFpEGmrAjkQrWW5hXOSpWdJiI2Di7eTppFF6PBCcKdBc02/We81YVeDk6eja35+friCojHNteaS+VK2l9R6axyWFnOsbTueRgok/l7OtWi/avMqomzUlJHz5g/ICEd0K/oSjCMA/gTAzwH4GeX6QyvfOBeA6wE8BeBrlHdvB7AEYOm6666bqBFytlcA9OsdoC7Ty/+ynNIE1gajpTFENAEtv6aZRGy/Gp4asAb/JMyM4+VMX54/zMETeJKZyY+qxmkDrS+9umrtRFsUzM/PjzDqWm2Q9zNpxceOHRvi0fbdKdXNmytePnkfyRthoJG5w2n3VjyluRqpm4QaE47lmJ+U+ec8ngB4CMAtxrvPWflqLwBtAKcB3FlKO60PwbzOLG3pyieLN0AtxqThKnWypoWU6mc538bReGo0GYlrHFo9kIyQ+xpkfS3HrCWYJBOxfCFRQaz1QcQW/Mgjj+TV1dUhHtqQzhJQVjvJD4lWV1dHdu7UhGKpTrV2aUsAe3ik8lKaG6V+9e75c0/b1uiICHvNJBSlaZowjgB4DYCXG+/+rJWv5gJwFYCPAnhHJP1WCwDvnQwnkxtSedqex0ijk1CLdPDol7hqzFjEPDW8/H5aAzc6IWiZzLen0BzD3MQmVwpa3WroKgksTejK91IgSdCiSrSrpBBofq7oyXbyntp7fn6+uJOrlk9rA8v3xRWnjY0N8xwOWR+t7aPmTp7HKseKILIEgdbXUoiXAhimNc8sAfAiGJBzfjzn/Kzx7r9a+aKQUkoA/h2Ax3LO/2ZSfFHo9Xo4c+YMOp2O+r7Vam161ul0cN999w3ztFotPPfcc1heXh7iXFpaQq/Xc8tutVrYv3//SBn0n+jScPR6PSwvL+Ps2bNuGbxuvIxer4ezZ8/iK1/5ippHg5WVFfR6veF7iW9paQmLi4vFOmtlyHuJW2uHXq+Hc+fO4Y477sDu3buxf/9+tNttHDhwAAcPHhxpx5WVFezbtw979+5Fu90e4so5D99H6OY46V7rQ492OdZ6vR7uu+8+PPPMM8M27HQ6w7amtP3psbl9CGev1xumkWOZym+325vqMzs7q9Lu1YHX+aqrrsLevXtH2obGAeXh7UX9JPHzdzMzMyN9knNGq9XCzMwMAOD8+fNmf/V522g7dTodLC4uAsCQbqJPG7OyfzWYnZ3FgQMHhrh4mb1eD/v37x9pOyrr9OnT6PV6Q7wppWH70VicmZnZVC7lP3XqVIi3jA2aVMiXNPS3AngUwB8PriUAf9fLE70AfDuADGAZwJnB9T1enmmtAMZx6kptX27i5EWnWBK+pGXLdxGt3DI5Wdqot8e7p+3Kulp0SRwlDVrDR8+idmNOu7b8run3Wnq1lZdm4yfa6EyA48eP5/n5+U0hjbzuWqihFpVS8/GZRbOVTvORLSwsDLdULmm0WnnymVytlOLzNRrlPlz8mwUvNFWjyzMFbmxsDA+eIeczpV1dXR1pr7W1NTWM2YLnYwXgMei/C+CTAL4TwJ8B8LUAvgt9e/1UhEDtNU0T0DgHaVDHa+YYjWF6zFdO3CjNMgbZGphaGk6T1g6Scco6yLpadZbLck9AlOpck89rf62+UbzR7wa8NFr/U15iGjS2PIGs0aMJi5JQtsaIZbO3BGi32x2eaSC/LeDpvMi2cQW+Nb41GuXhQFoez1en0czNYoSbTEXcXJnzJcf+3NxcUQBuBYwjAB4GcL3y/HoAD1v5tvKa5m6g4zQ4l/jW3vzaMy9CJUIHDTj5lWyEkcmBbEU0SQHnffxDGqvlfNVirnk9vHpKWq26WnUu4R9nL/kSHeMIJqutS3kjikOJ1tJHhJpAobGnjY3jx4+PaL4Sl4aX0vDttq16eOAJMy2N1y4cT4QmyqMpUTlv/mCMHPtaO5RWTZPCOAJgZZx3W3lt14EwEjwtTYI1+Pm7qJZTYqBR2iODTRNWnI7V1dXwJMk5jxzQ7bWTltfSBrWw1lqBQc8ss4DEGY1YKgmmEkPiIMeOZLZW+bXmO28M8zbieSkKSxtXkegX7cvvmrFuzY9on2lCWGPqHi5Zrta3VIYW78/TWXWbFMYRAKfHebeV16R7AY3TkNHIm1IIWS1+i0HUakXRMrXBbqUhLbTEiKz0NcyZnmv2b9KqOAMo2aG1cjh9VnuUaJT4StEn3r0XtkptEfWHlN7LtrPyaEySzrwgs4bG/Czm6aW3DrrxaNPqwkOGrTHP9zTS2s5SBrR25NFnVt8RvV6deLtETcQlGEcAfHngoJXXpwD8sZVvK69xBQBn0KXG589r7PSWVlXKU/PxUkQjoLqSXTJSpmScJQZKNMpJ55Vj+S6ioNWf02tpgDy/Vz53zHpt5dHH88iwYc8Rzu+5MLKc11Z7R6FWudDqRPXijlU5HiQDJwGrhe/KvJpw8OaEXJXzkFVpj9fCY722iry3BJsnfKw6kVksMg+jMI4AeJV3Wfm28pp0BSCXpBEmK5e9kXI8k48Eyx4otVEvUkECOZwiqwu6l4wzymA8wVTSbkp0yecW4+MCPiL0rDpqB8WUJr+GXxMi3pJf1ksztdTUzyqnlFarl6RJwyPnFaWXB74TEyYnqRY559Ek6SmVzx3TPEonotRFxjUJObmCsAQ60eXNZ+rr6GZ9UagWAJfjNc0TwbT/NFFrnDQaWBq8tc2EfOaZEGpoiKSPMI8Ig48yf68sTZPjz0uTMcLU+L3cJM0S3FpdPaHimS4kndGxxem06q/liUR0WXn4fKjpf8LBjz/N+RLzk/vy1PjBtHrx9vGUN/48Wl6k37UVHs+rCR9PuEZMiLUwzgrg7wP4n9n95wF8CcAGgB+18m3ltdVOYGnD41DbAdrgkSYomU7TZOT/GkFUw1w0mrxBHWHm2juPwfIwOm3waxN63OgJcmJrbe7VR/7XGITn6I7uc6+9sxihRUfOo2YyrT+t/Dlf2pBOO+LUo5mERrd7aRsLyfAkA/Xqze+1eVHKZ5k8owJN5vHmsBSgcpxIBdMyM1v9NC6MIwAWAfwP7P6Tg9+XAPhtK99WXtNeAVjv5WAbF4/sNM3EICeo1GqskD1JhzUZagWZxxTkcre2fayBXIoW4Xlr0vI0sh/uueeevLq6GmoTmcZiwBEznUe3Vw+NDrqXbeKFeWq0WhonXwFwwebRRaYXvuso5eMrAPrYTdsVl5dtMW76L48Ttdpb0inneXRcybpbbaG1NbWPDDjwlKIa068HE0cBAbib/V+08m3lNakA8JwqEc0h8s7C2e12h+FyEo/HUKX2oNHhTeDSsxLd8rmmSdbisOpQAjl5vbpQGnI2Wg5Jj1aPDo1RSKbi5a8Zh14a2X7cIW6l5e8t5zf9csXE23fHOi6S4yFNlzuOpYChMi5cuDCicGjtXbLna8xXChAuhLiA8j6m0557Y1OmLx0F6ZU1LowjAM4bz18E4DNWvq28pnUkpPWuZvJFGSCfaNIxa2knGkSXilH6S1BiYB4TrmFwtXSW2suanB6+qACi9FaYYiSKyhs3nBaL0VqCncqX20J4KwhLG5VfzZZWfBsbm09/08rmYcM8EogUI6ofX6FRuvn5eXVTOM+5zu9Lgs4S5F595LuIOZG35TjzZ1wYRwC8B8DPK89/HsAvW/m28prWkZDas0hja5Ld6kC5o6Fcrnq4rEHuMalpDpoSLo9WnkbeRwWtZg6jdKWT2Twhb0HEyanVRWp+1hbHpfEiaeEaaa2JSXMwemNcMiRiykRHKT8XiF77EW5u7tnY2BjujyO/KpdCaW5urvj1vdVGXNBFBb0GpfFeUjY0gWvVpySQamEcAXANgH8P4DyA/zS4zgO4F1M4E3ica9pO4ChD4s9KGgeliYS4eUySM/2abxFq3tWm996X2kxzbvN3Umu9cOGC6rSU9mIJ0lTBcdbUN5qH18v7WlkbB95/mZ4YppfPosmLKdeEC29rrrxYqwlOr+cs5oyfm+Zo0zQ+bySNCwsLqr1f9pOlaJCpKbqNhvXOC24orf5y1v2ApfpM41SwscNAAdwA4K8NrhtL6bfymsYKoMS4tE7mTDkSPVBrUpAgGeY4zF8OfjmoSgM1Qic9j7SZtvS1mAvZdvmhKJJBSbu+nDDac69uFoMuCU6NMfPfiNnRaj8Ca1XgrZQ4/VYbaO2q1Sv6kZx0AHvjgsxBvA5cwMvxopnctLbynhGNVrqo0NeeTeLLi9A4CfPPeQIBcDldk/oAShE1PK31zItMoPsapurZFa332mT1zElycFqfv9fQqTFmLZ/8tVYCMj+FIHobzvGPZbQ+KLUvx2Wd3lbSLqUjUstjgWeP1vpPS7O6uqqaRkrCRyuXl8nHSI1ZjTTtUhrpnOXpuamE1z0yryLtztNa49Gj3yvXeifbUyuP6KjdqTgCO14A5LxZ09OYlzfQJLPgzyXOWnq8NPJXE2TewNVw1A4wS+CUwvA8PF4aYgTWZJB+FV43yTg0/JbQttrRYhbSt2PVkepj0WCVIfHx/2tra/n48eP5yJEjKoOuXflZAQrR+UHvvbOwZT2sOcNNoBG/CcdVs7r1FAWaa54CEQWtPS1BP41D4CU0AkABTcPyPtSh9LSzpTUwIlqIJywsbU5OZG/wWuVN4gST+Gr8E5Sn9E4KWS+Ppn0vLCzko0ePbvIheLTUapfeM0sYE128rSymbZWtOYXJkerRQ7/exmceE+T/tT18ZHraf6eGCde8521MZifLDBfFa/WtdITT82kGXMj/05qjHCYSAIPTu/7e4P/LAXxzJN+0r634EphPLOpoTVPjaWTUgzY4vKUeT6s9i0xELqgiTC7n/sCyzle16LGgdrCWBJ70eUQjc2T5GxsbI6dTWWlJkGtMTSsnUj9LgaCyIysAi0FrH+FJxq0pChKPZGZeTLqmjcoAB0t7jyob0TbXFASioeZsDv4uot1vbGxsOoejhHtcsNp0GmVN4gT+GQD/L4DfG9x/A4CTpXxbcU3qA9Ce0UCVXnzLxEO//L9lTrDswiW6Spovn/xeGlmmt0lcjVbjCanIJLbeaZO7hkZPy+WrlW63m+fm5vLdd989DA2UJrVaUx43iUUZolY/fui6ZAIlP4M1Lvl45DvFUpvQF9Ecl6bo8PoQWCZAjX6rP6lPPNu452ez5pCnQNBY02iUZdGKZhIojSWtv+TzSWASAXAGQKKtIAbPlkv5tuKaZDtoT8vRJoqmHXjLc60cy9GllU35acBpKwyatOMsrYme6ISN4pYTZVIHVkm4WRM/why0FYCWR+aL0q1t9kXtTl+2euV0u93hFgmcOdI7nl8zS8jvECgNP7BldXV1RAmQ/hKpBGlCi0chWSZATq83NrT5R8y55Ez2+kKLypFM32K0Ho2yHA2PRU+pDlY7Tsr8c7YFwItQht4AQQaAlNI1gTyXFbRaLezfvx8AcObMGfR6veG7drs9fNdqtYa//Spfgk6ng/vuuw/PPPOMip+X02q1Rv73ej2srKyMlEvPZmZmhmnOnDkDANi7d+8QJz3vdDpYWlrC2bNnh3kkEH5eDof9+/dj9+7daj4qOwqc5l6vh1arhdnZWRw4cMCkK4oXGO0nnl+2S6/X20QHT3Pq1Klh/na7PYJn9+7dI3knhfPnz+OOO+5Aq9UaoWfPnj2Ym5vDF7/4xWFaSTc9I/pbrRZmZmZGxg3/PXv2LJaWlkbo/uIXv4jl5WV0Oh089NBDOH36NDqdDp544gkAwMzMDHbv3o277roL7XYbvV4Pu3fvHpYD9McIzYlWq7Vp3PJx3W63cejQoZF25fXjc2h2dhb79u3blI6Xxcf0xYsXcfbs2ap+obREI9WRt3FKaaSttXrJucVp5GWdOXMG6+vr6jjkwHF7tGs8onZeVoMmFfgF4CcB/AqAzwD4hwBOAfixUr6tuLZqM7jSEpNgdXW1uPdIbbkcrCgMafe18HMfhbTterb/KH2RfFpbjrOE9fJr8daWGYC0XU9701Zi45y4xfuPm1XonqJ1eF9Js44VUUSavDVG19bW8jvf+c7htxP861oyS9EJXladtXp5KxYvyorCdHm95AdN1q+snzamtDJLIdX8vxe+bNVXQsRcXAPjzsEIYNwVQM75fwNwH/pfAr8GwE/nnA9vlUDaatCksNTWLal77bXXDjUeKeU7nY6rSVrlEmgaAGn+9NzTIGQdzp07N7K6OH/+vLtCoHScnohmLGmS2k5E++HAtTgtP63EOp3OJg2u1+uh0+kMn7Xbbdx1110AMKItU1/RL4BNqyquKVp0Su2d+ml9fR3vec97sL6+Pky/e/du3HjjjSPa/enTp7G8vDzUHmdmZkb6SdZr165dI885tNtt3H777di9ezdarRYOHjyIW265BUB/ZQL0V5a7du0ariB4neXKitMu69zpdNDpdPD+978fDz300KY26vV6uHjxIh577LFh254/fx6vfe1rMTs7OzLPZF/w/uc0yfkpx7Ic/xLkClGzCGj18OYArQxoBURjr2ZVqa1uOUTnzdigSQV5oX8K2BsH/18K4Ksj+aZ9bcVeQOOkkz4CaXOM4IhIe01DikRXkIan3VuakqdpR+pTosfCYdlhPbBspfPz8yOnQNFz6fjkbUiOUW37iFJ7e9qx3HdmY2Nj5JhObuemlQp3IEvN/+jRo5uctV60CKddq5e1s2fOl1a6WoQQb7/5+flNNGnlUrt74dUaHVSetomi1wbaXPRWzvQ77jkgVL+jR48OV9qR1XZ0vE+yqiDAuCuAlNI/RH8F8CuDR98I4IGtEEZbCb1eD4uLi0U7OeBL3V6vN9TaeHquCZQ0AGu1IFcEUjvT7MI8PQdKQ7R4GjmVxVcMXKOVbRcB3s5LS0tqnWU7RW2l1vtdu3bhzW9+84gW3Wq18NrXvhZXX331sE579uwZ9lW73ca+fftGNFPKR20ibdFylWLRRSuVXq+H06dPD+mi1cuBAweG5bziFa/AuXPnRtqMcB48eBB33nknnnzyybDGK8eKpFlq4pSm1+vhqaeewh133IHdu3fjuuuuA3BJU5b28AceeADr6+sjK2dtLqyvrw9XVZxG7Vem2bNnz6bVHq+jfEf1ppWKNRdln2Xh96M0HvD5ddddd+G2225T/Qh8PvDVCK1ErPksV0rT8FVxiDiB/xGAb0P/NDDknJ8A8IqpUvE8AV/We4NCAz6BZmdnhx2tgcfIiAlxZivL8OjzHFIac9Bo0fKSqYmbRXg+OXmt9pHOsF6v7wy06PDaSQI5N6Xzk+PjTn0q//z580MHJJkjyJlP/SHbgu5brdZI3TWGyR2BlIYcwu12e8jE3/KWtwC4JBiAS2a+L3zhC7jpppuGbcGZEdVL0kHjZ3FxEZ1OZ/hL0G63hwKM2owzat4HRAu1Ubvdxvr6Og4fPoxOpzNMQ3larRZuvfVW3HjjjUPBJU11S0tLAPqm0/e+97244YYbqk0arVYLt91220i9uSlPU4ba7Tb27NmDs2fP4vTp09izZ4/pYOVzhgQyQafTwb333jsU4pSP3hFwpzif0/yXFChNaC0tLW1STrX0NabUMGjLAn4BeGTw+8nB7y68wMJACSwHGn/Pl7z8Kn1g45lPPDoicdClOkWcWSUc2i+HUn2oHtY5u9GlNDkMKWSRl8/3jtfopS+0pVmn2x0N0VxdXc333HOPuo2E/B7kwoULI+YiXjcyLUnzjdaGvH+leYS3b8RcxtuV2sX6YInSSzOg1f5UL2rHubm5TXNAXtIMRB/Wzc/Pj4Sf8jbkUPo+gveJ3LeJ4+MnvHHHM/W/VV8Lut1+8AQPlSW8ZCKTpjTeN7y9+FjUzFa8HtzEzPtjK0xAEQHwrwHcDeAcgL8M4D8D+BelfFtxTSsKqGQ75YNY29aZT2QNhxUfrXWg9SVmRBBY9utxwbJ/WsJJ1lkOdEuAauXxwU+Ml5h0t9u378/NzW3yZfBJubCwMLTBSmYtBQoxBM03Qs+eeOKJ/IM/+IP52LFjw4/o+OZ0lJaEimdDlvRY35lEvgKlck+cOJE/+MEPjkT2WOPAOj2Ljzv5DQEvXwq+hYWFPDc3t0mYSvo4/rW1tXzkyJGRaCY+VyQDlf0saeNj7sKFC/mHfuiH8oULF4bv6L0sg9NZssHz8cXbg3BpX0XLNuNty3HwNJpg0/poXJhEACT0wz//I/q+gH8IIJXyRS4AbwLwOPrnDPzTUvppCQDZSZ7W5E1Aj8lrHWg5rOTk0bQIiZ+niQyOUhqNyVtMX9NwpHZohfzJ8qzVT7d7SbMkxsGdqDwdab+kcWo0klZ14cKFYT75ZTQ9p4+xTp48mR944IEhM+FMik/W1dXV/M53vnOYltNH4ZC0SqB7S7GQtGjtQs8uXLiQ77777mGdvK+z+dessh4nTpwYCXGWDluqB9G2uro6rBP9t/qYYGOjvw3J/Pz8Jm2cQld5X1lnT/M5I1edFy5cGOabm5vb9MEbF7zauPZot8JGPWHtzWP+XuNJFs5xwRIArg8gpfQiAJ/KOR/NOf/NnPOhwf/N3pJKSCl9FYD/E8BfBTAD4G+nlGYmxesBt/tpH4dI0Jw59Lzdbg8/+tHw89AwzaZOeHjon+ZI43mkc5DS8Pf8l/5HwtJK4Ztk2+Tlan4Gug4cODB0NFJ+bjcFsOkDJN42i4uL+MxnPoP19XWcO3cO3/d934c3vOENm/wmDz300PB/znkTPfwjv2effRZ33303PvKRj+DUqVP41Kc+heuvv37EDkz2e6DvKH3961+P9fX1oS2Y+vQrX/nKSD/ddtttuOaaazY5jC9evIhWq4VXvepVaLVaw/BI8knQh2q8XeU4ILrIpj4z058mzzzzDH74h3946M/Ys2fPSPvyfs85D+369957Lz7ykY8M67hv376h47fVamFpaQnPPffcCJ6LFy9iaWkJTz31FN7znvcAAN785jcP8XQ6HZw6dWqkj3lAwfLyMlJK+At/4S+MOOmfeeYZfOhDH0Kv1xvauqmdlpaWNoWjklP4scceG6bft28fWq0Wrr322uFcvuaaa3DnnXcO25PPMemIlfND3nO/kpxHnm9NmyuyLnLecDyReTspuAIg5/ynAM6mlK7bgrK/Ff1zhz+Tc+6hf9LY925BOUOwmBdnVECs4Xu9Hh577LERh6RkmtLhS88kSEZKuPgvAZe9Gs1PPfXU0LEoGUmkfax7rU08h5QUnpyJ8sgJ7etLmtR79+7Fo48+ij/6oz8aiZDhTlia2K95zWtwzTWXPlLnTlESErfffju+//u/H9/+7d+OgwcP4uDBg3jd616HBx98EOvr6+h0Omi327jzzjuHjOE3f/M3cfjwYbzkJS/BAw88gI9//OP4wAc+gC9/+cvDaLDl5WXs2rULt95664gjsdfrDePuU0pDZ/C+ffuGjOlTn/rUiGOT6q+1S855xGm9f/Bl9/79+9Hr9XD//ffjE5/4xJBp8vF+4MCBYXvddttteOihh4bMmsezE7Pn/Xb27FncfPPNuHjxIp5++mn8g3/wD9But7F7927ccccdeOyxx/Dwww9jcXERDz/88CZGCvS/Qbjlllvw1FNPDedDp9PBhz70IXzDN3wDAAyDEFqt1rC8D3zgAyO09nr9iKovfelLAPqKAo/Ko7k8Ozs7VMK09iwpO1p0kTa/tegdTcBYc0XOMfnOi/SaBkSigL4ewO+mlD6eUnqQrimU/Y0APsfunx48G4GU0ttTSksppaVnn3124kI1xupJZsD+cEpqubzDJR7S4LiWtL6+jvvuu28Tfkv4aJEK/N11112H9773vbj22mtx9uxZANi01cQ4A4lrMwA2ReFoGpTMf+7cOVx//fWmoOPAoztIG+erNNKAW60W7rzzTjz66KN44IEHhlEmnU4H999/P9bW1nD69GmcO3duyNzf8IY3DD9QarVaQyb26KOP4t577x2uOD72sY+h1WrhbW97G37iJ34CN9xwA97ylrfg9ttvxw/8wA/g9a9//bDvb7jhBvzBH/zBSH2IUb/2ta/F8vLyiAbKI1JIQCwuLo6MDd4uPGyUKy+9Xm8YUUSCK6WE+++/X8Vz9uxZ7NmzB88++yy+4zu+A61Wa6TfTp06NVxlfOUrXxm243PPPYdWq4Wrr74aN9xwA55++mmcOnUKi4uLQ7pe97rX4SUveclQQSGc1M8rKytDxspDcO+880688Y1vBIBhtE6n08H58+dx66234i1veQtmZ2eHz3u9Hh5//HE8/vjjADCcg/KDMi0qqtVqYe/evWpEEJ+3so80pWd9fR2Li4vDEGdKR+3J+8iaG/y/Fxq+pSsBzS6UR+3036FdpXwBvH8TwHvZ/Q8COOzleb4+BJN5oh9ryI20ZH4eqSGdSxxPlE7NZshtnTKiJeKULtnqtQ96eD7PqS3vte2spb2ZHx1IURdkSya79sLCQv7gBz+4KeqH2oBHa/AIC27XP378+NDnQDZ96TTW7PXkENXs4Dw9p0N+IEa/8phLKlceJM5BHtdIdddooLJlJBW3+xOdZD+X0Sk0ZqkcOTa4c9SKkJE0UXncPyDTy/7kznjyx9E24ORjoTHC+5zXg5ehzVvZ17zP+CHzPCKLyuLtJsuVfi/u1+AwTV8ALrcDYQAcBPBRdv9TAH7KyzPpdtDeToClvNEyci4fBC6dPvxdKVqGgA80LY82kEtlyTbiuIjRWYdzUzoeEuk5sDlD4qGG0lFHTtZjx44NcdMEp8PEiSHIicXbgDt2aeIeOXIkHz58eJifT9xjx46NCIojR46M7KckGbfmzOZ1pbbhkUz0Xjp+eXtS/S28PJ1kmpqTUfYx0SWFGDEzrY/5XlPEmDXnthb5Yo231dXVPD8/P7yXbXLy5Mk8Nzc30ic8vxe1pI1r7b929rFsX/6M/6evtemUNu3Mau68tvpapp+U8ROMLQAAbKD/ERi/Pod+OOgNpfwO3l3obzD3zQBaAM4C+PNenmmcCWwxGS3iYBLQJg+nxercKKPmIY9RQRIpyzskgwa6dSJUt9vdxMz5O41JrK2t5bvvvnuEyfG8tInaL/3SL41Eu9BZuMeOHVOjViTzI+ZOk5M0SYrzJlp4Hj4mKI+MMqJ3Wr2pznLFQfckeLStDrxxyfHK9zIqTNZf63eiS4t111aNfBXL68W33LAEES+fM3Ja2VG/yGMueVrOTCMh0NpY1N5LBi3TlMqQQkurrzZfPbqi5ZdgEgHwzwH8MICvBvA1AN4O4KcB/C0AC6X8BdzfA+D3ADwJ4F2l9NP8EExqQV64VhQvgRf2ZeWR9HnPif6SkIkImFLZEoiBaxovZzSa5kmhmvIEJ/6RkDRHLCws5NXV1U3a38mTJ/Ov//qv53e/+915dXVVbQteNmdWpGkSs5EHkFNeKeh4Ws7IrZ0+LUFI5iaOi7+3Qg45cKbN66gd+MOZdgknbzdrfPF6SQHnKQdUBrXP8ePH1e9GFhYW8uHDhze1oVxdRMyammnHM9FqOCYBTxB6dFs4xgVLAEScwG/KOf9Kznkj5/ylnPOvAvienPN/APB1cW/DZsg5H8s5vzrnfGPO+V9MgisC3AFMzlQK2aQoiNJe+9pz6eSRkT9aOvmO3pPzSjqeuJOIomS0fdg5eM5az4lr1ZOcsjfeeCOWl5cBjIYp8ugLctSSo6zX6+GWW27BS1/6Utx8881DvJ1OB08++SQ6nQ4efvhhvP/978fDDz+Mixcvotfr4YknnhhGz5DTtd1u45WvfCUeeOABPPbYY7j33ntHok9kG9De+IuLiwCAq666CjfeeOPQMXno0KFhNA0P333iiSdGcJ4/fx5vetObho7pnPtROU8++eQwBJMcnbT9guxroO/QfvLJJ4f77ciABB6xwvtKjjFKy6Ntrr/++k1tQO1nBRBwnFQGRWlpDn05j6hfWq2WuWMpjWE+N1JKw7Mv+HXw4EEcOnRoGDhB8yGltCmCjM8BbUsO2vuJ76jL90jiAQmRXXd5P0RABpRo/EWb87wciWOqoEmFPKqlnwLwFvQjhl40+P/w4N2ZUv5pXtM8E5hrEVxjKX2wwZ9p/7V777nUQORynafRtBivLO8DN8umbOWRH1TJE8ukxs/v+bJasxmTeYWe8XrJNPxjpne/+935gx/84MhXuFL7lFokaeDcgaytsOSXwxKHTC9NB1p6jsvSSGU6opWvijRauTlK29k0ArJu2grPM5NYZXFcmllNAvlcjhw5Yu6QKsFyllt5rX6kvvTmlefbsoD6T/oFuH9N4rT8iOMAJjAB3YD+mcBfBPDs4P8eAFcD+PZS/mle0xIAfOmq2UytPPy/NWnHpYfbsTUTkra0LR0MX1pWyv/crFFKbzE1ab7gzmrebjyqQ5oeZN3JVn/s2LHhOb6EW2NY0sxgCUFuMpHlao5VzoSt8SCfWweKS9ONHE/0jreR5X/x+kZrT4vpakxHG3c8EqgENNc0urgJSOaxnNDW2IxE1Ml2l2bfbrfvw3rnO9+5aetu2QZa2Z5Q5JFCvF00HiTpp3EwiRAYWwBcTtc0TwTjjetFbmh5c45/Hm49k3RwLSDK0IlGTyvxgMrhZ8aWhFpJM5KTQ6sXf25pzgSc4fAVgEWT1K40mzgHKUQ8vNRWMmKIv5d1sXw1kvFrAkeOUXJ8W+PWAilMpXN3Y2Njkz/Cwy/t9155HvPSxrpXbumdvJeOclmmLJ/+ayG/1j21hzb2tbqWBKyFd9KztidZAbwawMcBfHpwvxfAPyvl24prGk7gknOWNFa5l4iclPQsMkCsNLyDPYFEeeREsoQQDX5LkOR8KeSN71NjTQqtPKt9tWdWOxE+7gTWGJ88TtBqW62+niNT5rdAtgHRzpkg78NSXl4uhbJq9ZL1kKaCWobA21Pi4KY5TyBSXTXGSnRGVhz8XY1JpVRnSZNm5qV6WkK3VGZE6E3q7OV9Fal3CSYRAJ9Af9uGT7Jnny7l24prK88Epl/OkLUVQMk8YzEi6z46YYipR2zHNEF5yCUHaS/lwoAidawJoi1XPUEj6ZKgReFowpYzWEkHPZcan2Xe04QR/5X14rRwZs3NM3ynUKvuMg0xfyufRk/ULGKBbENLI9XKpjBPuaqR+LR3Wj94q06LiVr00erQ26CxVM9JmaxWZi1uOeanZWqeRAAsDn4/yZ6dKeXbimuaTmACyUhLEtoLv7QmVgSvxKGl1yasRcPa2lqem5vb5Gij97RtLsdNzJ/isTXaOBPkbSK/tJT5vNWD9AHwZ1qbSP8Hx8PPEpBH9HlMX3uvrfgsDVzTqrV6cuWCBJ+lMGjO3MjqUgNLuFmMUOsrHgJs4bfGpGzXCHOT7WLNLy5crY8UvXlVA5H8WlvXmOssX9QkMIkA+C0ANwJ4dHB/CMBvlfJtxbUVAiDneCNrk1JLw3/pv2Wq4TRozy38kn6Zhkw7Fy5c2MQgtIgJLtwsYcg1X00wEGOT++bTR2taHS3mEDk/Vr7n0UWao7I0MTWtUqbX6JL9qwkOSSO1lZfH23LDahPLn6U9t/aytwRSrY/IE64yrYZHaxdrbHq+HI02K12Jpoi5sGRmLoGs46TMP+fJBMANAD4G4MsAPg/gBIDrS/m24toqAaANeC+td2+ltxgFMVxt6VoCKy/hpq9rNUajRS1QOn6Sk2RIXEDIenLcPJ8WvaLZ9nk7WdEhHhPSmG4JPGZRmsga48x5NKLGc3Jr+WtMMiWmqjFQ/k5btXnnOETGuja+I3k5DssXZK2sOa01ypzW1pEVdqQO1n2Urmk4fjlMHAUE4BoAXx1NvxXXpFtBWM9rnIDyeSkM0/P688Gq2dY94PZOzYQiy9QOIJGRICQ0aGsGiYM0+VJUh2Wj50DlWeGEVltoDshSnhqomaw8jcakotEhXvmaENE0RM+8xhlpVBCV8HrtwOnUBJAnwL0tLrT24en59yIR8OovBeS4oAnFiBCg+bZtAgDAO7zLyreV17gCoNTwJTuk59H3NmSjNFSGNjhLk80TLNKEo9nF+TsaWMREpX2caOP78mjamMbM5QS3aJH5au2b3a4eP17DrDzck+DQmFWJxlJZUojIFZ8sU9Jz8uTJ4QqMBHgp9Fark6RJpuGrWQ03Xz1qxzNyXJ5Q9dLTuBun/3iZRGvJDOnhsZ5F25ro2E4B8DOD6zcBPAHgFwfX74Ft4/x8Xlu5ArCcR17eiN2R0h0/fjwfPnzY/ZBHm9QlwVIzuEiroAgh78hGbsbwBKQVAePRotWjJp2FexqTZRx65HOPaUlBGVl9an3D8/NfTehwwS/9MTVCj3DwfY84TTy8UpbPv1DWxlaJKUb7e9IxIJUl/pzotb60LrUlV7a0ceEpB5PWaxIfwBw3/aC/KdxHSvm24toqH4AMH4s0NjG+iFOMa0eWduP9lpiNx2Tkc207AUujk9qYho/yT0tTseiOLp35bxR3bbrIRC+1CWmY8uMricdTAPj44NsQa0yYlzvOVhFcUaItGrQxYgk4fsn3tStBbaXh0a09q9XieT4ZBReNwrL8MbLcmnEQhUkEwDkAL2b3LwZwrpRvK66tEgA5j2pV0S8ctW0GrLTau8hviQloE9BbBmvapBc7Tf9L4YdbDR59nJ4aZhKtg+Uw92i1lAPORObn5/P8/PzwQzxPmGvPueOWQng5c6T31iE0NUDl8fzW+OTPpSCqNavI5/zwGMsvwOupRRFJR3ctlMai9W6SMbmdK4B3ob9X/88OTEJnUDi4ZauurRQAOV9a/ln7x0S0YpnWexcdHBYe/rEW//W+5tTwRAa01h6RdBEaSqBpWTXMchLarJDZ0gpAc+DJvqcrYgaR+bi5hXBYW5zLbZdr/WE01iymKZl8NPJIvpO/Vr/TlhjclCTTcwEklYPoirU0f0vpx1WSps38c55AAPTz4nUAfnxwfUskz1ZcWyEANObnTcRop3qdSBNWDtAaIcAnkzZ5Ingt7S1CS40dsxSuGIFJJ4VFW6Se1v5DtULbex55J/vccrhKfJ6ywf97DlpP4Jcc2yVFh7ezZNRaxFO3290k1KL1jI6jGkbvpa8Zp1x4S6E3qRAYxwnctt7VpJnmNe0zgUsRCVb+SZmPHPBRRhQdXLJetUwnwiwsWrQ2kvUeJ/rHo9dKVzJredsw8HS1mqBFZ6n9IpFkWjgvgbcjpQZaP0iNucS8rbpYDExrT+orHiXH9yayBJ2kWzM5eTBpH9a8K6Uh+snHMsmeTxqMIwA+Poj6eQOAa9jzGwD8fQAfBXDIyr8V16RRQFpna8vImvz0rmSPlIPWmxBW+RKnVo5GT0TIWGVqy2sCqwyvrWoc7Ro9EdvtxsbG0LZestWXzGUkKGqYgVb/UlsSPVb51oqRCwfv8HiLRkv54VE9VE70+EW5cR9vXxkpxJUhuTW3FmXH0/My+TnF4zqJrXawhFktWPk5/d1u193SYlwYywSE/pGNvwHgs+ifBbwG4KGBX+BaL+9WXNNeAXDQTCbWr8TJD4Smwanhsr56HcdpKfF6zJc/r3UAanioveRuqSWmMu4E0urqpeVbTdfgtpiyFhYYxSufeWPLYw6yfK2vLQHi0WiNefrPd6Ms7TAr8Wq4STnRytS0d40uLWqPb68R/SDMEtAENMbl4UeTgNa/3e7owT+TrJAtaM4DcMBirBFzBQ06Ho4XsaNqEQoefRIHz6udw6uBVpeSBiTp4GXwQzN4GZ5ZZZznERMCf6cx2ghoeGv8NBH8cvuFiH/E69MSgy8JfW3Mc0YtGa0nqDz6PEETob8UhSWjnLxQWI0eDSe1Qa3wt8rU3sv2Hne1XoJGABTA0tisd/Scf4ouHVcWSHtlZBJ7E87S3ry6aIIk6o/gdZdp+VI+CoTLmrCepiajn+QKaxpMW/s/DnDGZDFnjeZx6yGZivZes5vz8iKrCinYLBg3/DQyJjVmLiOf5IrQil6ie22zvBq6I0EHnrCfphBoBEABPAFQymdpdBoeOekig6tEmzbINKezpemV8lp0WFqVtM+WoNvtDveZ5zSUBILWlt6EqgEv3zj4uJOvJGi1/ZFqmUJUoEfGrEwv30tfSoke630kn1UfbRxzmqTZqOQ0lkItKryscUr0RIIOuEIzLSHQCAAHarWuqIZsaQ+RyRadtHKiW5NeY/aROmjx4x690TaUeeUKQKujNgm1MEFZ51L7a/R459V6DNNqm5KgIqAVFN9Ou1S2V48IRJkbMVG+VxSvW6m8km/IYsbS1l/60NEapxwPdxrLvtN+o+NYBkPwlYjH1LVxO43N6Agm+RDsVmUriL9YyrcV1+WwAvCYEU/Dt0COdKI20S2GZy3bNZprB7HMIyM2piksqT4ynabtWf4LbiPmaaQAsdrBopkYBH9maYOcDu/bjFJ70XN+SIzmm/LGlMb0SsI+0nfUJnS+hIyO8SKqiF6LeWv9ZUVMlWzypfnGaaFyeTvUrogs/PRBKe295I2HScuMwCQC4JMAErt/ER0O83xfWxkFVJM+4himzo9oDdpELzE8iaNkT9Ymhifg6J0WWx7FVQLpZCQ8pVO8tHIlg/TMbF4bc9p420lfj1Z/LjClYIoyZq0PZT96jFQySxlTznFabchx8TakMa2FotKuoxa9miCWNGs0aWm9bTNK81JzukolxFPuakxB1ryLrm4lvklgEgFwRnm2XMq3Fdek3wHU2NQig6mUP5JGMugow5PAY6A1OjQnWUnziwgWja6aNtaiOWR4aQ1OiwaLWVgTTzIrjXlxmmU5nJloh7lrTJwYaanfPa1e0qwdCWqVXRoj8oxpjyaeplSWNles+nFfild/C6RSwM1ZVjvQcy0qKlKmxBPxb2l5JhECkwiA+wH8EwBXDa4fB/BAKd9WXJMKgNoDFsZp8GkzqiiOkmbkaXc19GnPpSCL2NxLE41DRFiV6PPSeJpkSShrqxgJxDgtkxJ/Zm2WptGm4ZJapVZ3vgOpRY+srxwzJeEhcXvPeL2jbT4OQ5TtK1cAVlQep8cyGWnt6dFhPY/M33FgEgHwCgD3AvgCgP86OB/gFaV8W3FNKgC8/e8nAU97KJUzLXr4oJbPS4OqhrFa4Gk00UFttZ/8qEtO4BoGrtXTMqfI/xaukj/okUceyRcuXHC/D7H+l2zD5HegPf6lD0QLfZROZo3hen02zjj3gPpAmkaibRDBz8vR6C4Jck0x4XmnFX5sCd5tEwBbcQF492Cb6WUA/xnA10byTcMH4MVgl/Jaz63lrGREOeuxx5NquB4tpbSeryE66CJMVuKP7gpJe6Nop5eVTDjafy1d6fB1DpFoF4mHM2M59iJt5+EmvFrAgbaTKa8DLzuiwcrx7CkRkTayhKz1bhyomRMl27+11YY2hsalmdq0ZkffCEyyAnj1YF+gTw/u9wL4Z6V8BZzfDWDX4P8vAPiFSL5JBQCX8jXLSE0z4eAxTz7RZCwyzxth6KXnRGspTbTsUr4aXLJNvHA4yeRlXPw4TCJi6y+1HQmumi0RNFOMHHvjMg7OhPkHZhyk3V6WwwWUR4ccz7wtNGdqSWvXlCP+LrrvUARK88eqt9aWtYpaDfDxMs0Q0JwnEwCfAPCtAD7Jnn26lC96AfjrAH4jknYaYaCWqcIC3ilRoeENkJryPYGjDV4ZdWJNrknA06iiE8uiQ2MKJS0owsCs5b00mXANzytPqyvhk+Xz8MgSo+eKQgn4uOTlSMZBzMRaecoyZd24gJHvuMDxhBrdazuZWu1IZxdbEInE0/7zZ1pdS+cmTIspW3ion56vzeBehDK8NOf8O+LZxUC+KPwQgN+yXqaU3p5SWkopLT377LMTFdTr9dBut4f3rVarmKfVamH//v1otVpYWVnBzMxMMZ/2np5R+b1eD2fOnEGv19tEI/0uLS1twtPpdHDmzJlNOHq9HlZWVrBnzx4sLy9jcXFxUxmlMjWQ73h7ED5Ok3xH+WWbWG1E+em9hovT9olPfALLy8sj/SLTt9ttHDp0aKTve70ezp49O8zXbrdxxx134Ny5c+j1emYfE22cVqDfL/fddx86nc5I+a1WCymlYnvQsz179uDBBx9Ep9Mx+4DSzszMYGVlZTiuDxw4gNnZ2ZF2WFlZwU033TRMJ/MCwJ49e3D+/PlhP/Z6Pezfvx8AhvcppZE+AjDMT23M+0CrX0pppAytHalfnnvuOdx///0j7cDbw5uLcl7IsUP/tTHb55XYRBdvd1mWd6/R5o3nlZUVrK+vbxoDWwaaVMijGvpvAbgRg9h/AIcA/FYg38cAfFq5vpeleRf6PoBUwpcnXAFYEr0WR83zEh5vaco1Jvm+FEJmrRC09BpOjZ7auljPavKX0tFXs97WCR4+bWfWqHlHA96GXDuuXf3Vnj7mgTYm+DttJRGNjedptNWbVq+I34ZwlJzrHmjjnP5rq0xvHlhlaiveUl+V2shr70kAE5iAbhgw8y8D+DyAEwCuL+UL4H0rgFPorzBCeab1IVjUvBDFWRMRYaXXaKlhqrV1kUt5i6aSAJmEFqIhcpCJxmSkY3USQSxNHuOY+jg+7ivQyvB8OONEkpXGChd4JWGpmXOsQAHed15fkInKcqKXmOMkc1TWxzOHWWAJCElvqd+mzdijMLYAyJcY9jVgW0JMcgF4E4AVAC+vyTcNHwDvyFIETg1Oid9jjlI7lMxCw1lTpxohYP0vfSxk1TNKC5+MEc1LYw41bVSCElPUyvOYJDFELb/FhKndI5FN/Fnkwy6OO3K+QUTwywNYtFPoeL3oozTLd6KNK3pmnZBVMze0dqHyPf8FT2etyq15UMMbIvUZF6oFAIB3eJeVL3IBOA/gc+gfMH8GwC9H8k1rLyDZ0ZM6XDxNUD6nPYKkdugx2sgXzFKwRGn10kWYoFeG9z4y4axya99PY1JZjEMrh/pCczxrwowzKU9L90x11nN53+12hw7WcSJaNAWKgDuSJbPnNMo29CLjaJUhTwvjwqKkVUvc2riLRM91u/43I/J5abUQxa/hrgVLAHhO4K8eXLMAfhTANw6uHwEwE/UxaJBz3pNz/qac8/7B9SOT4KsF7qjqdDoTOVykw4mcW5ajJ/cFGQ4cOGA6zKQD6uLFSz53y8nU6/Vw3333FR26Gk0acEesfB7J66WV9Svh7HQ6bhqq1/r6+sgzzeEWqXuEZnLYcXy831ut1ibHM4F0sBLudruN2dnZkXccNzlutbEqaaE82n1KCY899pjqRC2NEU4r/VK+8+fPD9v8iSeeGKkrp5GXqbUTHxf79+/H7t27cdddd21yGvM2ob6X/cz7ieom68nTcCe3VnfuZOfPtXHGy+Vlc/pk2/d59ebn0XlbDZpUyKPa+hw27wb6kVK+rbi2YjdQqT2VJLQlnbWlOP1y7aJmi2GZp7SELGmBJQ3EAw2PZRMdF6/2rvSFJoE8iY1v3iY1Stkv49Dt5ZnGysfSMumw9Ii2aZXhjQNLQy7RK8vUxnx03GvlaGYauvgHb9L8KOtjrWA1v0DNaldLI+esNQ6pfBmuW2qTGsAETuBzAF7M7l8M4Fwp31Zc0xQAGkMtMeLIXuRyQEgGZjHpKDOoGQjjLJejeKhuJTt4DV4OMkZf5tNATmLPoaoJDI8RleoSbVfOlKJ9yftObi/M3/N2l6fURcqy+pnjLDFCOc663e5wmwovn0WDNna73e6wHfhe+5IOrz4azZ6gsHxQpagerY4anbSdhzXmtlMAvAvAWQA/C+Bn0LfZ/1Qp31Zc0/QBWB1a6siaSasNzpJjqEYQRNJKLSMyYDU6ogMzIkxleksQR0JTS4zXmnz84yfty8sS8+DgrVK0MUYO0Zq9qWgcaZE0VhmUxwsyqOln7+tt2Se87bjPK0oDCQ5t9UC01Pq8tP7Vwl29iDhNc48KAo02ruSMO/4iMLYA6OfF69DfBfTHAXxLJM9WXNNcAWiOIWky4O9LoOWJ7jETYZpyQGiTmtfNY6DeZOb5J9U6+ISICNaoxkNpo+YhC0fOfmx6lB6KgtHyW8xXtk0kEoozPeu8Wk0wam1rRWp5TKzUBlZaah8pGEo08EN5PEVjXLBWsDySyWpfnj6S1mobi8HLfpg0fHQiAXC5XFu1AqBOpwMvrPhtiUPDZ6XR7i3aSu80hilpkbZiSa8XbSDNBjUMIcpoch7VnniIXUlwSM2JPy+1sZxUpe0fSvhKQihqfil9C8H7TPafHKvetuCcLvmMBEvJ1ONFNsnn3W43z83NDVc8PK9HA5XjjR+rXA8ic1QKbm9cyHfaaW4Rmj1+4s3XKDQCYAAWQ6RfrcMjnWj9J/AGcY02Exn0hJMP4JrdCvkArjl9yWsXjfnLuHFiEKXTtzR83W53k1lFQikGfdy4bW35zvNHVhfcJOWVpf3nIFc0nBbtxC6tDP5ltWyr48ePj4RkWv3L68S/1pZto+WzhBtP4723wOsPDa/FByQeer6x0d9qm3wenhCN0Mr/155lIqERAHk8jcECbxBZjESblBFmp5XBf+VzLd048cWcmUQFV4RJ8Qlj0Rud2LItPU08YmaJPOPvtKgTyx4u80nmWloBlGjiQlSWw2PpPRwbGxvD7be5b4QzdBlp5dHDzWNa22i0lNpci7azVpjyv6eIefWR+fjc4HnW1tZCX7Z79bPq3AiAKa0AJpHMlM5z+njP5GCVzM6yJ3s0aPde+R4O+U5zcsl7jcFZDHdc/0JkUnp1LzH/aJkaDTKtx9AiGqZGh9YXGk5tpUerKu1QHa0srsFqJ2gR7SXzlmYek+NVMktvfnIcnkIiy5Z0lASW9lwqcLw/5EpbO4ZT+2+VU6rbONAIgAFYzKCWSVCnjHPKmJVubW0tHz58OB85ciREh8ZcS+llHeWE9AauFUUkJx5pUtoZsuMyYmki8XBFhTJ/p2mPJeFSErpe+khe2a5W28ttEuR77UxizwTEGaHFjEraP6dbqw+vE/fByHDIcdpWK7tWqYoIGOu/tTrx8PNySvUaBxoBwMBzSkZtzvTMisawwOv4hYWFPDc3N/Z+8B7DsjQuDnLweuWVaKI2kbbLcbWZ6GE6tf2gMTNL46sZFxZwnJ7w0sq0+owrDNHIHq+P+YdV3nj1cEYZmaRXMzvVMMtIOg+sciXz9mj0hAfPb63IagMbItAIAAE1g4Uzeis976ga/wA9J205qqlQOdwOKR2q1gTxnpc2pYswVg+fJ2SsehJYx/GNY0Io5am9J/B8DZzBalEmGl0l4EJWCuBSPu85jcfSecl0Lxmk1HS9fpHvSkKPnkcFvRd1U8sH6Nc7srEk/LrdfmSUtr+RDI6Yls+yEQBBsJhBydmldSB/r/2ne4qskBqXxkQ0TZ/MLnxJXXIaRZiONtAjZgup0fH8J0+ezB/+8IerJq6laWl1qTl9a1zwmAppdt4qS66SSiu3Et2akC6ZeCJftFvOf8/EodERcXCP4+QszSvvnUVXRPh476yVuEaDxSM8s9m40AgABt5AtD6wsjQYzVYqO7DENMlRFREkmp1emoKsPUWiGoxsD22pGommIVp4WTIaxQPJzEqTnOiNmhHG1f68NNq3DDK95n+ymG0posiCkvCxVpvUhvJ7GK2OVllav5TabZI4d43plgQOT6u1r4ajNF74XCQcnmCrfT4JNAJgADWDQz6Xk0Yu17Q8/DdCT4lRlMrR6JTldLujZgMPNNt7JKpCTgYOGm0WWIxR1keuFkoCzhIsWn94GqOGV8NPgomYL0+3sTEaPy7xyPaMtJvWXlwx8LRfSyBFyuXlyLaz5l5JQJTqpo0zj3YtVFb7ItwbF1r5vC7cNFsaf9Hnk0AjABiM05jENLVDty1twXOa1g76kkbC31uMipfHTUal8i1TlEWXVseaevEyvFWQN+G5MLDSeH0naRuHQfE8x48fz/Pz88OvYqXgIuGgOeq1sqOMgrdXyb8TiVCyxjMfD/Pz85s2apMC2io34nOQUKJZpuW0Ee2R6DRtzFnzraR8ePWyhP8kQqARABOCNnjlewneRPb28tHKzTmrWoYcdJ5WLmngOL08JbAGdM2ktGiI+EEsnNyJ6YXneRM/YqoqAZVBhwF5ba59yRvxfXhpeV+X2ozj0vpCM9/xVaIMSNBotaKUqM9Kgq5mDmrPtbpFxigfT3IOW8Kydk5Y7TQJ88+5EQAh8Dol2mkyjzWRPWHCcdAg41FClk3cmqBWPaJan8QRicbhz6NQo+lE03DGIs1eEbNETRhiqU6lPudjRjLjCH4vbUTp4AxYjl96Pz8/r45zfh5wyT+kKS3WmJZ1qhkjsk4e1PikrHuJ7+jRo+oWDprg1sbbNJ3BjQAowDhMnueLLANLg9sql5g/aY8ejaU9QzwGHaFlfn4+33333WqMeMkXEoEIDbXAJ5jWX56WVqMUWHWvoZlr0dHImBpmWBqfmmnEeq+VHxVWdFE9JbPzyqhZlVEZ3rzRTI2lukRACnKJm5vmvK3ja4WeBo0AcKBkNiGwGEOJwXtL+GhYJbcR1zD4Ur1qBxUJJE3zjBxOHsHvvRvHpGQ59zzGVct4LQVARmh5uLvdS/ZpTcBKjbGEy3tmCcKShhth9t5zydC0Mr29frRxJ8vVFBNp95d0yfYeZ6xZ7yMmPGL+pbqNC40AYCAnLt8SuJRHG8Aefu2+Nq02ECJOUPqdRAB59Gtl1gqgWnpK7aXhiDjFI05Arx20/BsbG8OPvuSxhVY7Wd8REGPkGqNXL0sYeWMp6nugPDVRLMSES4LL0pplf/KypPlJq6cXGs0P6JF+NosGr65aGSXgfb8V0AiAAWjMm0+4KKPUtBErbYme2rzSvqhpVKUoElkPi67a+kSghsGX8JQYu4VX5i2trmrKInxHjx7Nx44dG/YTZy6RaBqNiXFtUSu/1NcctxwrlgZqtZ0XMinTHz9+PB8+fNhkwlS+d7aCNrb5dyUevd64orw8GsgSwCWhoNEcgdIqZVJoBAADi0l6zM5jDBIi4WSUVzKhmgHjCaSIsKJ3JbNNLVP2no/Txh6+qMM7QmvtZC6Vsba2lhcWFtQvvLXvR2T/RcJT+b0W6qnVQX7oZZVRGlPRdiDmqu2TT8LBOtvXahuZ31NkIoK95Hjvdm3/Wq3SIIELHmqDaUSfcWgEgALjMKHI+9LSWX5EFdF+ojR52ppV35ovn0t0aFEMVEbEKRbpE16nUltrOMapV+SdxTQ0zS7yDUJtNA0x94gpy+tf6i/OlEtKQQkX4dHs+ydPnsyrq6sjeSwzkjU+xlEs+DiNnG2gzZPIyrFmfm5s9Df34xv8jTMXJTQCwIFo40YlemTAaTi1KIRJBA4vL8ooebmaXTtad6nlWoywhNtjBlqeUhSFVa/a8q10VlvW1tmLAvIYXsTfUQLO9GXb8bpGcEqhpDFE/h2BNkcoHfd9WGVFQatD7byOzBFrHMr0kg/w9i+ddBeBRgAYEJ3gPP00yilppCUmUKIpOkmlZiMHomZqiLaBFto3TvtFmZx06Ftt7EXsjKNBEtSeBFXSDEv9GmXynFF5woH+89BMjw6LqdUIJT7W6N5ygmtj0oNSGk8oafWMjA2rPax5XxqjFAI+qXO4EQAO0MCaBniaZS0jHZdhynK9ZzThvYFYyj8NemuYpqXBU9TNuKs0j9FZdMg+LZlt+H/vLAmNaVi/Gj1W2VGGVTKL8HEToSMqRDc2NtSvpXnbRsK2PUEpcUp/CM938uTJ4fYd436XYbWxt/EkF4SR3X1L0AgAByaxv+esmz34u5oJUXo3LliMxiuvxNS0yaJ9yl8qLyIQPQ0z51Ht3xJgUabP35Xix7VfC5fUdq3T5HhbysgUzTQTrUvtu5LN3dpbaJLxS9FTWkitbIfImNEEtVYPbZxQX5U+JItASbmw5gifW1ecAADwkwAygJdF0m/ll8CTMH/N8SnTjIuvlDaKU9PWIvisyBJNC5MDmiZX5ADvKPP0QO6jI/0QtSG71G4lrT4qvIihaBFAnrmMj63S/jNWubWhvF4/WdE4445bb2xpbWQxc6ssrx2sduQCprbdat5bY7MWbwkuOwEA4JsAfBTAhctBANRqh+Oki0LEjOANGm2QRcIDrckhw98sLY2/lw4tqw5RiKblGmLOo6eIjdPHkXK10Eb+X7arpt1qp0NZ9ETaz9N8o+DVqZTegqhfiSsPsl08oVR6LhUUL1+t1l/qO1l3TRmYNi8huBwFwH0A9gH47PMlALyBG9mFMIq3Vmvg4C1tI5NaS+NFp1j14Pk0HwG301r5SIPW9uYfRyu1aI2YC2SeGvyldFIrrzHN0HvrK1mrryPlRN+VoEbZiOCJRpaRAsEZf8SRrT33xqiW3yprXIXBEiqWuZhf04DLSgAAuAPAvx38dwUAgLcDWAKwdN11143dANbgo3fT2odcai61+CImh9LAKAkkTmMEjzc5vIgN7b8nwKIDnpdp+V24luVtm2DVtzaaJ8Ikahlyyd4dVVgmjeay6LMUJys90VJbXkmZidCprcK08EyvrEj5Xh01oWLNAe9UtnHgeRcAAD4G4NPK9b0AHgHwZ3JAAPBr0hVASbvWoCTRNaaldXSUAUW/5Cwdy2fRyp/ValESl3WQvWybkiOWnkeEkubs9ZzuvE1r2kgLJdXyREG2ScTMo0W+eGApN7VadxSk4sTL0sqMhnB67eMJG5mPl1/yzWn9EhWm8r0VUCL7otvt5mPHjpn1jCpEEbhsVgAAbgbwhQHj/yyAiwCeAnBtKe+09gKKvitJ/FLn8I6Oai6RDpeaSkSjKwkLC4cnPOiAE28ALywsVJ8BHGnXKJPQ7kvl0yReXV0Nm+NkWbIe1r3X3rXKg3XYTM2YmnRlwAVWyRRZEoClcrXVERfgXgCAVQetn6z0Fl4vQIK/W11dzW9961tHvoLWYNwAFQ6XjQDYRMDzuAKQYGm2/HlUa9FwR51w3gAtCY1IdM+kZVmTgGv/ss14upoBXDtZIyBtvVZ78bK1M2I9Wnh+crp7B9NT33kbjPG8JTo2NvpbCPB9cbwVqryvHUcaDd6Y1xirxTyj4Z2esH3kkUc2bS+hgSxPbgwZpUV7zvHQjqPyUKcSjaWAiyg0AkCANSmt51HQok6iGpiE0nKSBtY4wkn+H0fQaZqQR3MUrOiLcTRUmkCkHWtb/2p1kv8t3FaeEpOynvF3kTBE2fal8azhjYyjaCisZ/rQGLWWziuLfDoyCMGz8csyrPJorPA5XBrH1pjU6srNllbdteelQ54icNkKgJprq1YANRPTg7W1tXzPPfdsEgKTRLvUTCj+PgI08bVoGU94WcKqRHMNTZZGWmKG2i9njBSWWrLllsrkvgjPTFQTHRQZh15/lH457ZpgsCAi1In5ex/NTaJR0zvS7OkDvZIQ1vJbfdDtbo7GivosSuVE27DmeQ00AsCAbre8HXINWOfl8t9pQERzjGgvktHSM7lU1fDzSxvwWp5Ifbi2FJ2AXJujgzU0R7kmXDzG4NFQCi+VY8sbYxYN0szmRbJJB2zJpxOxx8u6WkC0emfgRsvx+pzXj0d4ae81GsZh5rXKgZXP4w0lGiaFRgAw0CaZpiV59+OUU1qSWvm09xazpfdRU0wJv/deMlNKK8uOTCJ5Fq085MPLKwXS8ePHRxzUsj0jbc73rp9kAvJ2LPW91o/SXBWJZCsxOxoftbbuEmhKh6VEWPRHHMV8nPCtOig/d97L+llt47WVpdxYvIPelb6AjzyLKkElaATAADxtj0AL4ay1a1uDudTZ9MzSgCXj9WirdWCW0mjvrYHq2Wc1PHwibWxsVNk9La0tytwkkEN1bm7O1OBrQBOWWpmaWW9tbc00V43LsCNjUz6PhCdbzEuWU9rqWtJZEhaSOfMVgiYINJNnNPCAm/4kk5d5ZPix1m/a/6hwqoFGADConRw5x3YglGCFwmnaEE/Lz4/1JqWkbRK6JjGDlRhSBAePuy85AmtxjyNIpPYfdZJ7k9sKSeRCn5fDhbtsY8veHhV6GlOywk95WSXhZZleNBoj4KWz5iSfT9Lswueb9oV6BHiEmFZfjtMSil77WzCJIGgEgAE1A9TS0iw82nPtXvtymA4Tr+n0SFrNRqppKxHa6dmkS1S5iqlh/pOUrTHhkn1d5rPo0fJ4jntisvLDOs2XwqOZLGYT1Sy9FSjXsC3BwmmymJvV5uOAZJjanJR1sOiNzi+NWVMfaKeblT7S1NqzJBQnHeuNAFDAatTSQOfaQ63GpT23BEpNiGeUBo/JeThK77wycy6boyLLcW8i1YJWjuVAlGXVmAssba9G8dDwW21mjWXPKSz/83aoZeoSh9ceNaC1PS/Le2fRG21vKRA9ph3BKXFYfhKL7nGgEQAGRCYyf0aMeX5+3l0RyHwRQeNNzknqUpNmHCFmAV+Glz6Pl4M/0lYeY4vUQSvHc1p6fR0dRzR+ItqnV06EMWt5eH041AqrEt3SPu8pHzxfSfnwBFzN/LLMcRpoPgONRqtOFr3avUb3NKARAFMEvgLwthu2JLr2f9xOn+Yg0XCOi19OzFqHdElb5fdRJmUJEaktcnwlxk9pSl/SRhyuHq0W/SUcFq3S/i23kbBoiIIUmiV6NcFYEtYElllT/vIVvOVjizJy+S5qSdBwaeGh04hA49AIAAGTNiwxNWu/EYu5eMu9Gtt3jQbDy4ikk3bscduqlkZJh/bfej9uO/C256GnkhFFmJe1GolEkUW0SklHhCl55ixOn7XdwDTmSRSfJXj5vdWGcjWu9QP3CXj+gUkEXumZRre2m+3x48fzsWPHXri7gW7FNc0ooJpIEw80rUPec43DWu5ZS3CJiwZHTZy8Z9u28lh1ikJEc67FU3o2aTn8uwH5HUDJL1ESFl7bR8xYJQWjFq+ke5pRV1q6cRQKTZBaAlam18aGnAOW2Wgc2kgJ8ASWhVuLVJqfn89HjhyZ+DB4gkYACIja7zWo1TSJCfPn3gqgNNlJg4nQQMt9qfFYgmZaUKuteukimnDNhLNwcoa/sLCw6SMjDcckBwnx9BG6PWZl0Wc9r/GflISOdm+NX0tQRsejNq80DVqrq3yuRV1ZdfLqx/HxZ6VViScU+XicBjQCQAGtA0oMKDLB5YCParBRh7CV36tjzpsFzDhRGZG0mtPMGuRygkY0Vou5WBMuwqBOnjy5ieFLhqXRpNmsLSZngeyj2mV/dIzJPFa9PNxeW3r1lgyaH3bCx0t0pSr73au3FwGlrbwjbWeNZdlG3n8tmqiW10ShEQAOWMzFe1arqXjvo3m8tCXNzIvyoN9aZqABj/zxGAenR27BWyq3JBz5/0jbaAw/Wm9r7PA48RJD40KQ4vs1h2yprhImNXFaY0X7L5m4Nz673Us+ExKino/Mgoh/Q0trjQVPabDqouWrUc54/mg46DjQCIAgRJj1JJI5qhFOqplrg9IKY6MJGTki0ju8olSOlk7midavBmontFZetN9ytrcj8HBxZqAJaxKq3hkDkXoRRIWSxRzl1h2y3yPtqeGy2sejj1/aWCIhY5Xn0anVp5TGE35efSzc04BGAEwZxmX+JXPTJPgpn6exeQzPm7Tdbt8xdffdd7v+hyjdFp01OKJQ0uxkPScV8BozLLW/R1POlzY+k1//evll2RyXdY6yx+h4uVIQSedqrZm0VH+vjnJvLE3g8oNyar6w73a7obN5vX7dKqZeA40AuEyglhHU4N3Y2BjuHhkVLt5/bSKNG5WgMRNLEE2zXTjeGq15kvJlH0c+JIrglOYkybQ1H4+1spN4+IlVUear/d/Y2LxFc2S8Sdw1/SBj5rUxxQVnjXOV2ibql7hcwRIAqf/uhQGzs7N5aWlpu8mYGvR6PbRarangWVxcxMWLF5FzxuzsLNrtdrGcXq+HpaUlzM7OAgDOnDmD/fv3D9PW0ueVI3Fb6Xu9HgBMpV1qaIy8r6WNp+90OlhZWcH+/ftHcNS0McfX6/Vw6tQpAH0lbu/evXjwwQdx6NAhtNvtIV6PZl42peP4SzTK5zSevvzlL+MNb3jDEA/1PTA6xrwyiJ7FxUUcOHDA7ZMzZ85gZmZm2L7amCI8ALC0tIScs4vXq+cLEVJKp3POs/L5i7aDmAb6MK1B1Wq1cODAARw8eBC33XbbJuZ/5syZkQnOgRSAVqu1afLICWnhKJWj4Zb4AaDT6eDMmTN2RZ2yI/DMM88U02hMiP4vLi5iaWnJLI8/p7rQs3a7jZmZGQAYPi/1DeFcX18f3vO8u3btws0334xWq4VWqzVk/rweFuOm+kjGT+9LNFrPb7rpJlx99dUjOKnv+f9SGURLSslsG46/3W6rY4zSELNvtVqYnZ0NM39OyxUJ2rLgcr2uBBPQdsCkvgYyP5ScZ967SBnk6BwnXynt6upqftvb3paPHTtmbhFcwluqo3S6WweMR00imt9Fs7uXPvjSQnI1h7I8SMWiS7YHt8VrEUwWRMtoYHJA4wNoYBKIMr+adxIin/NzvCW6JKyurm7amI6EmzxknP/KMi2QtnD+W6KTGKgmcLTIKy3SxcIvBYQWb8/PoIgC0cwdsLVbmDfw/IAlABoTUAMhoOWz9c5bflvvJPBlvJePTAa1JqNrr70Wu3fv3mQq2bt3L3LOWF5eHjHdSHOQZQqhi5tQVlZWhnb/KJ3c3MHLe+qppzaZo1ZWVjbRYpmTpHmE7rmpkNpl9+7dRToJyLQyOzuLVquFdrs9/N/ACwMaJ3ADL0jgTs5pOdKBzQ5QrUx+T873q666aoT5SfoidGr4rbzRZw00ANhO4F3bQUwDDUwKnpNzEnweTs2JTZElFo4aOi0nubWyijxroAEPGgHQQAMTQMN0G3ghQ+MDaKCBBhrYodAIgAYaaKCBHQqNAGiggQYa2KGwbQIgpfRjKaXHU0q/m1L619tFRwMNNNDAToVtcQKnlL4TwPcC2Jtz7qaUXrEddDTQQAMN7GTYrhXAjwL4VznnLgDknL+wTXQ00EADDexY2C4B8GoAr08pPZJS+kRK6YCVMKX09pTSUkpp6dlnn30eSWyggQYauLJhy0xAKaWPAbhWefWuQblfB+BWAAcAfCCldENWPkvOOf8qgF8d4Hw2pXRhTJJeBuCLY+Z9oUJT550BTZ13BkxS51dpD7dlK4iU0kfQNwEtDO6fBHBrznnLVPyU0pL2KfSVDE2ddwY0dd4ZsBV13i4T0AMAvgsAUkqvBtDCzpPmDTTQQAPbCtu1FcSvAfi1lNKnAfQAvFUz/zTQQAMNNLB1sC0CIOfcA/ADz3Oxv/o8l3c5QFPnnQFNnXcGTL3OL6jtoBtooIEGGpgeNFtBNNBAAw3sUGgEQAMNNNDADoUrXgCklP5DSunM4PpsSukMe/dTKaXzgz2J/so2kjl1sPZauhLrnFL62ZTS51k/fw97d8XVl0NK6SdTSjml9DL27Iqsc0rp51JKy4M+nkspfQN7d6XW+d0ppXODev/nlNLXsneT11k7KPhKvQD8IoCfHvyfAXAWwIsBfDOAJwF81XbTOKV6fieAjwF48eD+FVdynQH8LICfVJ5fkfVl9fsmAB8FcAHAy670OgP4Gvb/nwD45R1Q5+8GsGvw/xcA/MI063zFrwAIUv/E7bcA+PeDR98L4N6cczfn/PsAzgP41u2ib8pg7bV0JddZgyu9vr8E4H8BwCM5rtg655y/xG6vwaV6X8l1nss5XxzcPgzglYP/U6nzjhEAAF4P4L/mnJ8Y3H8jgM+x908Pnl0JYO21dCXX+R8Plsm/llL6usGzK7a+KaU7AHw+53xWvLpi6wwAKaV/kVL6HIC/A+CnB4+v6Doz+CEAvzX4P5U6XxFnAnv7DuWcPzj4/7dxSfsHgKSkf8HExI6z1xJewHUu1PcIgJ9Dvy4/h76p74fwAq4vUKzz3eibBzZlU55dEXXOOX8w5/wuAO9KKf0UgH8M4Gdwhdd5kOZdAC4C+A3KpqSvrvMVIQByzm/03qeUdgG4E8At7PHT6NtQCV4J4A+mT93WgFfnlNKPArg/942Fv5NS+lP0N5J6wda51McEKaWjAD40uH3B1hew65xSuhl9u+/ZvmUTrwTwaErpW3GF1lmB3wTwYfQFwBVd55TSWwG8GcDtgzkNTKnOO8UE9EYA53LOT7NnDwK4K6X04pTSNwP4cwB+Z1uomz48AH2vpSuyzimlr2e3fx3Apwf/r8j65pw/lXN+Rc75+pzz9egzg9flnJ/BFVpnAEgp/Tl2eweAc4P/V3Kd3wTgnQDuyDl/mb2aSp2viBVAAO7CqPkHOeffTSl9AMAK+kurf5Rz/u/bQdwWgLXX0pVa53+dUtqP/hL4swB+GLji+1iFK7zO/yql9BoAf4p+5NOPAFd8nf8P9CN95gervYdzzj8yrTo3W0E00EADDexQ2CkmoAYaaKCBBgQ0AqCBBhpoYIdCIwAaaKCBBnYoNAKggQYaaGCHQiMAGmiggQZ2KDQCoIErHlJK/3tK6Q3K87+UUvqQlmc7IKX0jsHOj59KKZ1NKf2blNJVg3efHTz/VEppJaX08ymlF7O8H0kp/aGsT0rpXhE/30ADQ2gEQAMvOEh9CI3dlNJuALfmnH97i2n6qgnz/wj6WzvcmnO+Gf0tPL4A4GqW7DsH774VwA0YPSLw3QB+UEF9BP0N4xpoYBM0AqCBFwSklK5PKT2WUnoPgEcBfFNK6UhKaSn1zzz450bWQwA+wvC8aaBln0B/exB6fs1gI7nFlNInU0rfO3j+0pTSBwYbzf2HwQZ7s4N3nZTS/5pSegTAwZTSD6SUfmewX/2vkFBIKX13SulUSunRlNJ/TCm1FTrfBeBHc85/CPTPzc45/yuxAyYG7zrofwT1fQMBh5zzxwFsKHj/C4A3DrZDaaCBEWgEQAMvJHgNgP8n5/wtOecL6G+WNQtgL4DvSCntVfJ8G4DTAJBSegmAowD+Gvq7w/INuN4F4P/LOR9A/zyFd6eUrgHwPwL4bznnvehvNMf3k7oGwKdzzn8RwBqAvwXg23LO+wH8dwB/J/UPavlnAN6Yc34dgCUA7+AEppS+GkB7sK1vCAaC4ffR3wLAS/en6G8VvC+Ku4GdA40AaOCFBBdyzg+z+7eklB4F8EkAfx79QzIkfD2AZwf/bwLw+znnJwZbY7yfpftuAP809U+MWwDwEgDXAfh2APcCQM750wCWWZ7/DuA/Df7fjr5wWBzguB19M82tA7pODp6/FcCrBI0JbCfHlNJfSZdOsLvNaQ9tR0gNvgDgG4qpGthx0CwLG3ghwR/Tn8EGWD8J4EDO+b+llP4v9Jm2hOfEc2vvkwTgb+ScHx95ONiAxYA/YfuvJAD/d875p0T+vwZgPuf8ty0kOecvpZT+OKX0zTnn3885fxTARwcO3ZZKbH/VcD2A33PoI3gJ+u3QQAMj0KwAGnihwtegLxD+KKX0ZwH8VSPdYwD2DP6fA/DNKaUbB/ecKX8UwI8Rw08pfcvg+Qn0T5JDSmkGwM1GOR8HcCil9IpB2t0ppVehf4rTt6WU9gyev3SwQ6uEfwngSBqc+TqgQxNoGPgQ3gPggZzzfzPo4fBqAL8bSNfADoNGADTwgoTBSVifRJ+x/RqAk0bSDwP4S4M8fwLg7QA+PHACX2Dpfg7AVQCWB7uo/tzg+XsAvDyltIz+trzLAP5IoWcFfVv/3CDtPICvzzk/C+BtAP794PnD6JuiJBxB/xznRwbpTg7q90mW5viAtt8B8BQGu54CQErpvwD4jwBuTyk9nQaHhA+E43M551WjfRrYwdDsBtrAFQ8DZv9mirCpzPtVAK7KOf/JYOXwcQCvzjn3pkzmlkBK6X8C8KWc87/bbloauPyg8QE0sBPgJ9B36P7hGHlfir7mfRX6dv4ffaEw/wH8IYBf324iGrg8oVkBNNBAAw3sUGh8AA000EADOxQaAdBAAw00sEOhEQANNNBAAzsUGgHQQAMNNLBDoREADTTQQAM7FP5/BOHOfsXpNw0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "x = candidate_df['phi1']\n", "y = candidate_df['phi2']\n", "\n", "plt.plot(x, y, 'ko', markersize=0.3, alpha=0.3)\n", "\n", "plt.xlabel('ra (degree GD1)')\n", "plt.ylabel('dec (degree GD1)');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the same figure we saw at the end of the previous notebook. GD-1 is visible against the background stars, but we will be able to see it more clearly after selecting based on photometry data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting photometry data\n", "\n", "The Gaia dataset contains some photometry data, including the variable `bp_rp`, which we used in the original query to select stars with BP - RP color between -0.75 and 2.\n", "\n", "Selecting stars with `bp-rp` less than 2 excludes many class M dwarf stars, which are low temperature, low luminosity. A star like that at GD-1's distance would be hard to detect, so if it is detected, it it more likely to be in the foreground.\n", "\n", "Now, to select stars with the age and metal richness we expect in GD-1, we will use `g - i` color and apparent `g`-band magnitude, which are available from the Pan-STARRS survey." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Conveniently, the Gaia server provides data from Pan-STARRS as a table in the same database we have been using, so we can access it by making ADQL queries.\n", "\n", "In general, choosing a star from the Gaia catalog and finding the corresponding star in the Pan-STARRS catalog is not easy. This kind of cross matching is not always possible, because a star might appear in one catalog and not the other. And even when both stars are present, there might not be a clear one-to-one relationship between stars in the two catalogs.\n", "\n", "Fortunately, smart people have worked on this problem, and the Gaia database includes cross-matching tables that suggest a best neighbor in the Pan-STARRS catalog for many stars in the Gaia catalog.\n", "\n", "[This document describes the cross matching process](https://gea.esac.esa.int/archive/documentation/GDR2/Catalogue_consolidation/chap_cu9val_cu9val/ssec_cu9xma/sssec_cu9xma_extcat.html). Briefly, it uses a cone search to find possible matches in approximately the right position, then uses attributes like color and magnitude to choose pairs of observations most likely to be the same star." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Joining tables\n", "\n", "So the hard part of cross-matching has been done for us. Using the results is a little tricky, but it gives us a chance to learn about one of the most important tools for working with databases: \"joining\" tables.\n", "\n", "In general, a \"join\" is an operation where you match up records from one table with records from another table using as a \"key\" a piece of information that is common to both tables, usually some kind of ID code.\n", "\n", "In this example:\n", "\n", "* Stars in the Gaia dataset are identified by `source_id`.\n", "\n", "* Stars in the Pan-STARRS dataset are identified by `obj_id`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For each candidate star we have selected so far, we have the `source_id`; the goal is to find the `obj_id` for the same star (we hope) in the Pan-STARRS catalog.\n", "\n", "To do that we will:\n", "\n", "1. Make a table that contains the `source_id` for each candidate star and upload the table to the Gaia server;\n", "\n", "2. Use the `JOIN` operator to look up each `source_id` in the `gaiadr2.panstarrs1_best_neighbour` table, which contains the `obj_id` of the best match for each star in the Gaia catalog; then\n", "\n", "3. Use the `JOIN` operator again to look up each `obj_id` in the `panstarrs1_original_valid` table, which contains the Pan-STARRS photometry data we want.\n", "\n", "Let's start with the first step, uploading a table." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preparing a table for uploading\n", "\n", "For each candidate star, we want to find the corresponding row in the `gaiadr2.panstarrs1_best_neighbour` table.\n", "\n", "In order to do that, we have to:\n", "\n", "1. Write the table in a local file as an XML VOTable, which is a format suitable for transmitting a table over a network.\n", "\n", "2. Write an ADQL query that refers to the uploaded table.\n", "\n", "3. Change the way we submit the job so it uploads the table before running the query." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step is not too difficult because Astropy provides a function called `writeto` that can write a `Table` in `XML`.\n", "\n", "[The documentation of this process is here](https://docs.astropy.org/en/stable/io/votable/).\n", "\n", "First we have to convert our Pandas `DataFrame` to an Astropy `Table`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "astropy.table.table.Table" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from astropy.table import Table\n", "\n", "candidate_table = Table.from_pandas(candidate_df)\n", "type(candidate_table)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To write the file, we can use `Table.write` with `format='votable'`, [as described here](https://docs.astropy.org/en/stable/io/unified.html#vo-tables)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "table_id = candidate_table[['source_id']]\n", "table_id.write('candidate_df.xml', format='votable', overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that we select a single column from the table, `source_id`.\n", "We could write the entire table to a file, but that would take longer to transmit over the network, and we really only need one column.\n", "\n", "This process, taking a structure like a `Table` and translating it into a form that can be transmitted over a network, is called [serialization](https://en.wikipedia.org/wiki/Serialization).\n", "\n", "XML is one of the most common serialization formats. One nice feature is that XML data is plain text, as opposed to binary digits, so you can read the file we just wrote:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\r\n", "\r\n", "\r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n" ] } ], "source": [ "!head candidate_df.xml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "XML is a general format, so different XML files contain different kinds of data. In order to read an XML file, it's not enough to know that it's XML; you also have to know the data format, which is called a [schema](https://en.wikipedia.org/wiki/XML_schema).\n", "\n", "In this example, the schema is VOTable; notice that one of the first tags in the file specifies the schema, and even includes the URL where you can get its definition.\n", "\n", "So this is an example of a self-documenting format.\n", "\n", "A drawback of XML is that it tends to be big, which is why we wrote just the `source_id` column rather than the whole table.\n", "The size of the file is about 750 KB, so that's not too bad." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-rw-rw-r-- 1 downey downey 396K Dec 29 11:50 candidate_df.xml\r\n" ] } ], "source": [ "!ls -lh candidate_df.xml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you are using Windows, `ls` might not work; in that case, try:\n", "\n", "```\n", "!dir candidate_df.xml\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "There's a gotcha here we want to warn you about. Why do you think we used double brackets to specify the column we wanted? What happens if you use single brackets?\n", "\n", "Run these code snippets to find out.\n", "\n", "```\n", "table_id = candidate_table[['source_id']]\n", "print(type(table_id))\n", "```\n", "\n", "```\n", "column = candidate_table['source_id']\n", "print(type(column))\n", "```\n", "\n", "```\n", "# This one should cause an error\n", "column.write('candidate_df.xml', \n", " format='votable', \n", " overwrite=True)\n", "```" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "# Solution\n", "\n", "# table_id is a Table\n", "\n", "# column is a Column\n", "\n", "# Column does not provide `write`, so you get an AttributeError" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Uploading a table\n", "\n", "The next step is to upload this table to the Gaia server and use it as part of a query.\n", "\n", "[Here's the documentation that explains how to run a query with an uploaded table](https://astroquery.readthedocs.io/en/latest/gaia/gaia.html#synchronous-query-on-an-on-the-fly-uploaded-table).\n", "\n", "In the spirit of incremental development and testing, let's start with the simplest possible query." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "query = \"\"\"SELECT *\n", "FROM tap_upload.candidate_df\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This query downloads all rows and all columns from the uploaded table. The name of the table has two parts: `tap_upload` specifies a table that was uploaded using TAP+ (remember that's the name of the protocol we're using to talk to the Gaia server).\n", "\n", "And `candidate_df` is the name of the table, which we get to choose (unlike `tap_upload`, which we didn't get to choose).\n", "\n", "Here's how we run the query:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created TAP+ (v1.2.1) - Connection:\n", "\tHost: gea.esac.esa.int\n", "\tUse HTTPS: True\n", "\tPort: 443\n", "\tSSL Port: 443\n", "Created TAP+ (v1.2.1) - Connection:\n", "\tHost: geadata.esac.esa.int\n", "\tUse HTTPS: True\n", "\tPort: 443\n", "\tSSL Port: 443\n", "INFO: Query finished. [astroquery.utils.tap.core]\n" ] } ], "source": [ "from astroquery.gaia import Gaia\n", "\n", "job = Gaia.launch_job_async(query=query, \n", " upload_resource='candidate_df.xml', \n", " upload_table_name='candidate_df')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`upload_resource` specifies the name of the file we want to upload, which is the file we just wrote.\n", "\n", "`upload_table_name` is the name we assign to this table, which is the name we used in the query.\n", "\n", "And here are the results:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "Table length=7346\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
source_id
int64
635559124339440000
635860218726658176
635674126383965568
635535454774983040
635497276810313600
635614168640132864
635821843194387840
635551706931167104
635518889086133376
635580294233854464
...
612282738058264960
612485911486166656
612386332668697600
612296172717818624
612250375480101760
612394926899159168
612288854091187712
612428870024913152
612256418500423168
612429144902815104
" ], "text/plain": [ "\n", " source_id \n", " int64 \n", "------------------\n", "635559124339440000\n", "635860218726658176\n", "635674126383965568\n", "635535454774983040\n", "635497276810313600\n", "635614168640132864\n", "635821843194387840\n", "635551706931167104\n", "635518889086133376\n", "635580294233854464\n", " ...\n", "612282738058264960\n", "612485911486166656\n", "612386332668697600\n", "612296172717818624\n", "612250375480101760\n", "612394926899159168\n", "612288854091187712\n", "612428870024913152\n", "612256418500423168\n", "612429144902815104" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results = job.get_results()\n", "results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If things go according to plan, the result should contain the same rows and columns as the uploaded table." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(7346, 7346)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(table_id), len(results)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['source_id']" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "table_id.colnames" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['source_id']" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results.colnames" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we uploaded a table and then downloaded it again, so that's not too useful.\n", "\n", "But now that we can upload a table, we can join it with other tables on the Gaia server." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Joining with an uploaded table\n", "\n", "Here's the first example of a query that contains a `JOIN` clause." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "query1 = \"\"\"SELECT *\n", "FROM gaiadr2.panstarrs1_best_neighbour as best\n", "JOIN tap_upload.candidate_df as candidate_df\n", " ON best.source_id = candidate_df.source_id\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's break that down one clause at a time:\n", "\n", "* `SELECT *` means we will download all columns from both tables.\n", "\n", "* `FROM gaiadr2.panstarrs1_best_neighbour as best` means that we'll get the columns from the Pan-STARRS best neighbor table, which we'll refer to using the short name `best`.\n", "\n", "* `JOIN tap_upload.candidate_df as candidate_df` means that we'll also get columns from the uploaded table, which we'll refer to using the short name `candidate_df`.\n", "\n", "* `ON best.source_id = candidate_df.source_id` specifies that we will use `source_id ` to match up the rows from the two tables.\n", "\n", "Here's the [documentation of the best neighbor table](https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_crossmatches/ssec_dm_panstarrs1_best_neighbour.html).\n", "\n", "Let's run the query:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INFO: Query finished. [astroquery.utils.tap.core]\n" ] } ], "source": [ "job1 = Gaia.launch_job_async(query=query1, \n", " upload_resource='candidate_df.xml', \n", " upload_table_name='candidate_df')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And get the results." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=3724\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
source_idoriginal_ext_source_idangular_distancenumber_of_neighboursnumber_of_matesbest_neighbour_multiplicitygaia_astrometric_paramssource_id_2
arcsec
int64int64float64int32int16int16int16int64
6358602187266581761309113851876713490.0536670358954670841015635860218726658176
6356741263839655681308313884284887200.0388102681415775161015635674126383965568
6355354547749830401306313783776573690.0343230288289910761015635535454774983040
6354972768103136001308113804456319300.047202554132500061015635497276810313600
6356141686401328641305713959221401350.0203041897099641431015635614168640132864
6355986079743697921303413920912795130.0365246268534030541015635598607974369792
6357376618354965761310013993335021360.0366268278207166061015635737661835496576
6358509458927486721320113986549341470.0211787423933783961015635850945892748672
6356005321197136641304213922858936230.045188209150430151015635600532119713664
........................
6122417812491246081297513437559955610.042357158300018151015612241781249124608
6123321473614430721301413414585387770.022652498590129771015612332147361443072
6124267440168024321305213468524656560.032476530099618431015612426744016802432
6123317393403417601301113412177938390.0360642408180257351015612331739340341760
6122827380582649601297413404459335190.0252932373534968981015612282738058264960
6123863326686976001303513545702197740.020103160014030861015612386332668697600
6122961727178186241296913380061687800.0512642120258362051015612296172717818624
6122503754801017601297413464758974640.0317837403475309051015612250375480101760
6123949268991591681305813551997517950.040191748305466981015612394926899159168
6122564185004231681299313490752973100.0092427896695131561015612256418500423168
" ], "text/plain": [ "\n", " source_id original_ext_source_id ... source_id_2 \n", " ... \n", " int64 int64 ... int64 \n", "------------------ ---------------------- ... ------------------\n", "635860218726658176 130911385187671349 ... 635860218726658176\n", "635674126383965568 130831388428488720 ... 635674126383965568\n", "635535454774983040 130631378377657369 ... 635535454774983040\n", "635497276810313600 130811380445631930 ... 635497276810313600\n", "635614168640132864 130571395922140135 ... 635614168640132864\n", "635598607974369792 130341392091279513 ... 635598607974369792\n", "635737661835496576 131001399333502136 ... 635737661835496576\n", "635850945892748672 132011398654934147 ... 635850945892748672\n", "635600532119713664 130421392285893623 ... 635600532119713664\n", " ... ... ... ...\n", "612241781249124608 129751343755995561 ... 612241781249124608\n", "612332147361443072 130141341458538777 ... 612332147361443072\n", "612426744016802432 130521346852465656 ... 612426744016802432\n", "612331739340341760 130111341217793839 ... 612331739340341760\n", "612282738058264960 129741340445933519 ... 612282738058264960\n", "612386332668697600 130351354570219774 ... 612386332668697600\n", "612296172717818624 129691338006168780 ... 612296172717818624\n", "612250375480101760 129741346475897464 ... 612250375480101760\n", "612394926899159168 130581355199751795 ... 612394926899159168\n", "612256418500423168 129931349075297310 ... 612256418500423168" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results1 = job1.get_results()\n", "results1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This table contains all of the columns from the best neighbor table, plus the single column from the uploaded table." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "['source_id',\n", " 'original_ext_source_id',\n", " 'angular_distance',\n", " 'number_of_neighbours',\n", " 'number_of_mates',\n", " 'best_neighbour_multiplicity',\n", " 'gaia_astrometric_params',\n", " 'source_id_2']" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results1.colnames" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because one of the column names appears in both tables, the second instance of `source_id` has been appended with the suffix `_2`.\n", "\n", "The length of `results1` is about 3000, which means we were not able to find matches for all stars in the list of candidates." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3724" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(results1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get more information about the matching process, we can inspect `best_neighbour_multiplicity`, which indicates for each star in Gaia how many stars in Pan-STARRS are equally likely matches." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "<MaskedColumn name='best_neighbour_multiplicity' dtype='int16' description='Number of neighbours with same probability as best neighbour' length=3724>\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
1
1
1
1
1
1
1
1
1
1
1
1
...
1
1
1
1
1
1
1
1
1
1
1
1
" ], "text/plain": [ "\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", "...\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results1['best_neighbour_multiplicity']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like most of the values are `1`, which is good; that means that for each candidate star we have identified exactly one source in Pan-STARRS that is likely to be the same star.\n", "\n", "To check whether there are any values other than `1`, we can convert this column to a Pandas `Series` and use `describe`, which we saw in in Lesson 3." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "count 3724.0\n", "mean 1.0\n", "std 0.0\n", "min 1.0\n", "25% 1.0\n", "50% 1.0\n", "75% 1.0\n", "max 1.0\n", "dtype: float64" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "\n", "multiplicity = pd.Series(results1['best_neighbour_multiplicity'])\n", "multiplicity.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In fact, `1` is the only value in the `Series`, so every candidate star has a single best match.\n", "\n", "Similarly, `number_of_mates` indicates the number of *other* stars in Gaia that match with the same star in Pan-STARRS." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "count 3724.0\n", "mean 0.0\n", "std 0.0\n", "min 0.0\n", "25% 0.0\n", "50% 0.0\n", "75% 0.0\n", "max 0.0\n", "dtype: float64" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mates = pd.Series(results1['number_of_mates'])\n", "mates.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All values in this column are `0`, which means that for each match we found in Pan-STARRS, there are no other stars in Gaia that also match. \n", "\n", "**Detail:** The table also contains `number_of_neighbors` which is the number of stars in Pan-STARRS that match in terms of position, before using other criteria to choose the most likely match." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting the photometry data\n", "\n", "The most important column in `results1` is `original_ext_source_id` which is the `obj_id` we will use to look up the likely matches in Pan-STARRS to get photometry data.\n", "\n", "The process is similar to what we just did to look up the matches. We will:\n", "\n", "1. Make a table that contains `source_id` and `original_ext_source_id`.\n", "\n", "2. Write the table to an XML VOTable file.\n", "\n", "3. Write a query that joins the uploaded table with `gaiadr2.panstarrs1_original_valid` and selects the photometry data we want.\n", "\n", "4. Run the query using the uploaded table.\n", "\n", "Since we've done everything here before, we'll do these steps as an exercise." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "Select `source_id` and `original_ext_source_id` from `results1` and write the resulting table as a file named `external.xml`." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "# Solution\n", "\n", "table_ext = results1[['source_id', 'original_ext_source_id']]\n", "table_ext.write('external.xml', format='votable', overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use `!head` to confirm that the file exists and contains an XML VOTable." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\r\n", "\r\n", "\r\n", " \r\n", " \r\n", " \r\n", " \r\n", " Unique Gaia source identifier\r\n", " \r\n" ] } ], "source": [ "!head external.xml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "Read [the documentation of the Pan-STARRS table](https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_external_catalogues/ssec_dm_panstarrs1_original_valid.html) and make note of `obj_id`, which contains the object IDs we'll use to find the rows we want.\n", "\n", "Write a query that uses each value of `original_ext_source_id` from the uploaded table to find a row in `gaiadr2.panstarrs1_original_valid` with the same value in `obj_id`, and select all columns from both tables.\n", "\n", "Suggestion: Develop and test your query incrementally. For example:\n", "\n", "1. Write a query that downloads all columns from the uploaded table. Test to make sure we can read the uploaded table.\n", "\n", "2. Write a query that downloads the first 10 rows from `gaiadr2.panstarrs1_original_valid`. Test to make sure we can access Pan-STARRS data.\n", "\n", "3. Write a query that joins the two tables and selects all columns. Test that the join works as expected.\n", "\n", "\n", "As a bonus exercise, write a query that joins the two tables and selects just the columns we need:\n", "\n", "* `source_id` from the uploaded table\n", "\n", "* `g_mean_psf_mag` from `gaiadr2.panstarrs1_original_valid`\n", "\n", "* `i_mean_psf_mag` from `gaiadr2.panstarrs1_original_valid`\n", "\n", "Hint: When you select a column from a join, you have to specify which table the column is in." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "# Solution\n", "\n", "# First test\n", "\n", "query2 = \"\"\"SELECT *\n", "FROM tap_upload.external as external\n", "\"\"\"\n", "\n", "# Second test\n", "\n", "query2 = \"\"\"SELECT TOP 10\n", "FROM gaiadr2.panstarrs1_original_valid\n", "\"\"\"\n", "\n", "# Third test\n", "\n", "query2 = \"\"\"SELECT *\n", "FROM gaiadr2.panstarrs1_original_valid as ps\n", "JOIN tap_upload.external as external\n", " ON ps.obj_id = external.original_ext_source_id\n", "\"\"\"\n", "\n", "# Complete query\n", "\n", "query2 = \"\"\"SELECT\n", "external.source_id, ps.g_mean_psf_mag, ps.i_mean_psf_mag\n", "FROM gaiadr2.panstarrs1_original_valid as ps\n", "JOIN tap_upload.external as external\n", " ON ps.obj_id = external.original_ext_source_id\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's how we launch the job and get the results." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INFO: Query finished. [astroquery.utils.tap.core]\n" ] } ], "source": [ "job2 = Gaia.launch_job_async(query=query2, \n", " upload_resource='external.xml', \n", " upload_table_name='external')" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=3724\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
source_idg_mean_psf_magi_mean_psf_mag
mag
int64float64float64
63586021872665817617.897800445556617.5174007415771
63567412638396556819.287300109863317.6781005859375
63553545477498304016.923799514770516.478099822998
63549727681031360019.924200057983418.3339996337891
63561416864013286416.151599884033214.6662998199463
63559860797436979216.522399902343816.1375007629395
63573766183549657614.503299713134813.9849004745483
63585094589274867216.517499923706116.0450000762939
63560053211971366420.450599670410219.5177001953125
.........
61224178124912460820.234399795532218.6518001556396
61233214736144307221.384899139404320.3076000213623
61242674401680243217.828100204467817.4281005859375
61233173934034176021.865699768066419.5223007202148
61228273805826496022.515199661254919.9743995666504
61238633266869760019.379299163818417.9923000335693
61229617271781862417.494400024414116.926700592041
61225037548010176015.333000183105514.6280002593994
61239492689915916816.441400527954115.8212003707886
61225641850042316820.871599197387719.9612007141113
" ], "text/plain": [ "\n", " source_id g_mean_psf_mag i_mean_psf_mag \n", " mag \n", " int64 float64 float64 \n", "------------------ ---------------- ----------------\n", "635860218726658176 17.8978004455566 17.5174007415771\n", "635674126383965568 19.2873001098633 17.6781005859375\n", "635535454774983040 16.9237995147705 16.478099822998\n", "635497276810313600 19.9242000579834 18.3339996337891\n", "635614168640132864 16.1515998840332 14.6662998199463\n", "635598607974369792 16.5223999023438 16.1375007629395\n", "635737661835496576 14.5032997131348 13.9849004745483\n", "635850945892748672 16.5174999237061 16.0450000762939\n", "635600532119713664 20.4505996704102 19.5177001953125\n", " ... ... ...\n", "612241781249124608 20.2343997955322 18.6518001556396\n", "612332147361443072 21.3848991394043 20.3076000213623\n", "612426744016802432 17.8281002044678 17.4281005859375\n", "612331739340341760 21.8656997680664 19.5223007202148\n", "612282738058264960 22.5151996612549 19.9743995666504\n", "612386332668697600 19.3792991638184 17.9923000335693\n", "612296172717818624 17.4944000244141 16.926700592041\n", "612250375480101760 15.3330001831055 14.6280002593994\n", "612394926899159168 16.4414005279541 15.8212003707886\n", "612256418500423168 20.8715991973877 19.9612007141113" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results2 = job2.get_results()\n", "results2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "Optional Challenge: Do both joins in one query.\n", "\n", "There's an [example here](https://github.com/smoh/Getting-started-with-Gaia/blob/master/gaia-adql-snippets.md) you could start with." ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "# Solution\n", "\n", "query3 = \"\"\"SELECT\n", "candidate_df.source_id, ps.g_mean_psf_mag, ps.i_mean_psf_mag\n", "FROM tap_upload.candidate_df as candidate_df\n", "JOIN gaiadr2.panstarrs1_best_neighbour as best\n", " ON best.source_id = candidate_df.source_id\n", "JOIN gaiadr2.panstarrs1_original_valid as ps\n", " ON ps.obj_id = best.original_ext_source_id\n", "\"\"\"\n", "\n", "# job3 = Gaia.launch_job_async(query=query3, \n", "# upload_resource='candidate_df.xml', \n", "# upload_table_name='candidate_df')\n", "\n", "# results3 = job3.get_results()\n", "# results3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Write the data\n", "\n", "Since we have the data in an Astropy `Table`, let's store it in a FITS file." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "filename = 'gd1_photo.fits'\n", "results2.write(filename, overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check that the file exists, and see how big it is." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-rw-rw-r-- 1 downey downey 96K Dec 29 11:51 gd1_photo.fits\r\n" ] } ], "source": [ "!ls -lh gd1_photo.fits" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At around 175 KB, it is smaller than some of the other files we've been working with.\n", "\n", "If you are using Windows, `ls` might not work; in that case, try:\n", "\n", "```\n", "!dir gd1_photo.fits\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "\n", "In this notebook, we used database `JOIN` operations to select photometry data for the stars we've identified as candidates to be in GD-1.\n", "\n", "In the next notebook, we'll use this data for a second round of selection, identifying stars that have photometry data consistent with GD-1." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Best practice\n", "\n", "* Use `JOIN` operations to combine data from multiple tables in a databased, using some kind of identifier to match up records from one table with records from another.\n", "\n", "* This is another example of a practice we saw in the previous notebook, moving the computation to the data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 4 }