{
"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",
"- \"Write ADQL queries involving `JOIN` operations.\"\n",
"\n",
"keypoints:\n",
"\n",
"- \"Use `JOIN` operations to combine data from multiple tables in a database, 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": [
"# 5. 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": 28,
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Getting photometry data\n",
"\n",
"The Gaia dataset contains some photometry data, including the variable `bp_rp`, which contains BP-RP color (the difference in mean flux between the BP and RP bands).\n",
"We use this variable to select stars with `bp_rp` between -0.75 and 2, which excludes many class M dwarf stars.\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": [
"## The best neighbor table\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. Use the `JOIN` operator to look up each `source_id` in the `panstarrs1_best_neighbour` table, which contains the `obj_id` of the best match for each star in the Gaia catalog; then\n",
"\n",
"2. 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",
"Before we get to the `JOIN` operation, let's explore these tables.\n",
"Here's the metadata for `panstarrs1_best_neighbour`."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Retrieving table 'gaiadr2.panstarrs1_best_neighbour'\n",
"Parsing table 'gaiadr2.panstarrs1_best_neighbour'...\n",
"Done.\n"
]
}
],
"source": [
"from astroquery.gaia import Gaia\n",
"\n",
"meta = Gaia.load_table('gaiadr2.panstarrs1_best_neighbour')"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"TAP Table name: gaiadr2.gaiadr2.panstarrs1_best_neighbour\n",
"Description: Pan-STARRS1 BestNeighbour table lists each matched Gaia object with its\n",
"best neighbour in the external catalogue.\n",
"There are 1 327 157 objects in the filtered version of Pan-STARRS1 used\n",
"to compute this cross-match that have too early epochMean.\n",
"Num. columns: 7\n"
]
}
],
"source": [
"print(meta)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here are the columns."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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": [
"for column in meta.columns:\n",
" print(column.name)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the [documentation for these variables](https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_crossmatches/ssec_dm_panstarrs1_best_neighbour.html) .\n",
"\n",
"The ones we'll use are:\n",
"\n",
"* `source_id`, which we will match up with `source_id` in the Gaia table.\n",
"\n",
"* `number_of_neighbours`, which indicates how many sources in Pan-STARRS are matched with this source in Gaia.\n",
"\n",
"* `number_of_mates`, which indicates the number of *other* sources in Gaia that are matched with the same source in Pan-STARRS.\n",
"\n",
"* `original_ext_source_id`, which we will match up with `obj_id` in the Pan-STARRS table.\n",
"\n",
"Ideally, `number_of_neighbours` should be 1 and `number_of_mates` should be 0; in that case, there is a one-to-one match between the source in Gaia and the corresponding source in Pan-STARRS.\n",
"\n",
"Here's a query that selects these columns and returns the first 5 rows."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"query = \"\"\"SELECT \n",
"TOP 5\n",
"source_id, number_of_neighbours, number_of_mates, original_ext_source_id\n",
"FROM gaiadr2.panstarrs1_best_neighbour\n",
"\"\"\""
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
}
],
"source": [
"job = Gaia.launch_job_async(query=query)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"Table length=5\n",
"
\n",
"
source_id
number_of_neighbours
number_of_mates
original_ext_source_id
\n",
"
int64
int32
int16
int64
\n",
"
6745938972433480704
1
0
69742925668851205
\n",
"
6030466788955954048
1
0
69742509325691172
\n",
"
6756488099308169600
1
0
69742879438541228
\n",
"
6700154994715046016
1
0
69743055581721207
\n",
"
6757061941303252736
1
0
69742856540241198
\n",
"
"
],
"text/plain": [
"
\n",
" source_id number_of_neighbours number_of_mates original_ext_source_id\n",
" int64 int32 int16 int64 \n",
"------------------- -------------------- --------------- ----------------------\n",
"6745938972433480704 1 0 69742925668851205\n",
"6030466788955954048 1 0 69742509325691172\n",
"6756488099308169600 1 0 69742879438541228\n",
"6700154994715046016 1 0 69743055581721207\n",
"6757061941303252736 1 0 69742856540241198"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results = job.get_results()\n",
"results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Pan-STARRS table\n",
"\n",
"Here's the metadata for the table that contains the Pan-STARRS data."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Retrieving table 'gaiadr2.panstarrs1_original_valid'\n",
"Parsing table 'gaiadr2.panstarrs1_original_valid'...\n",
"Done.\n"
]
}
],
"source": [
"meta = Gaia.load_table('gaiadr2.panstarrs1_original_valid')"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"TAP Table name: gaiadr2.gaiadr2.panstarrs1_original_valid\n",
"Description: The Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) is\n",
"a system for wide-field astronomical imaging developed and operated by\n",
"the Institute for Astronomy at the University of Hawaii. Pan-STARRS1\n",
"(PS1) is the first part of Pan-STARRS to be completed and is the basis\n",
"for Data Release 1 (DR1). The PS1 survey used a 1.8 meter telescope and\n",
"its 1.4 Gigapixel camera to image the sky in five broadband filters (g,\n",
"r, i, z, y).\n",
"\n",
"The current table contains a filtered subsample of the 10 723 304 629\n",
"entries listed in the original ObjectThin table.\n",
"We used only ObjectThin and MeanObject tables to extract\n",
"panstarrs1OriginalValid table, this means that objects detected only in\n",
"stack images are not included here. The main reason for us to avoid the\n",
"use of objects detected in stack images is that their astrometry is not\n",
"as good as the mean objects astrometry: “The stack positions (raStack,\n",
"decStack) have considerably larger systematic astrometric errors than\n",
"the mean epoch positions (raMean, decMean).” The astrometry for the\n",
"MeanObject positions uses Gaia DR1 as a reference catalog, while the\n",
"stack positions use 2MASS as a reference catalog.\n",
"\n",
"In details, we filtered out all objects where:\n",
"\n",
"- nDetections = 1\n",
"\n",
"- no good quality data in Pan-STARRS, objInfoFlag 33554432 not set\n",
"\n",
"- mean astrometry could not be measured, objInfoFlag 524288 set\n",
"\n",
"- stack position used for mean astrometry, objInfoFlag 1048576 set\n",
"\n",
"- error on all magnitudes equal to 0 or to -999;\n",
"\n",
"- all magnitudes set to -999;\n",
"\n",
"- error on RA or DEC greater than 1 arcsec.\n",
"\n",
"The number of objects in panstarrs1OriginalValid is 2 264 263 282.\n",
"\n",
"The panstarrs1OriginalValid table contains only a subset of the columns\n",
"available in the combined ObjectThin and MeanObject tables. A\n",
"description of the original ObjectThin and MeanObjects tables can be\n",
"found at:\n",
"https://outerspace.stsci.edu/display/PANSTARRS/PS1+Database+object+and+detection+tables\n",
"\n",
"Download:\n",
"http://mastweb.stsci.edu/ps1casjobs/home.aspx\n",
"Documentation:\n",
"https://outerspace.stsci.edu/display/PANSTARRS\n",
"http://pswww.ifa.hawaii.edu/pswww/\n",
"References:\n",
"The Pan-STARRS1 Surveys, Chambers, K.C., et al. 2016, arXiv:1612.05560\n",
"Pan-STARRS Data Processing System, Magnier, E. A., et al. 2016,\n",
"arXiv:1612.05240\n",
"Pan-STARRS Pixel Processing: Detrending, Warping, Stacking, Waters, C.\n",
"Z., et al. 2016, arXiv:1612.05245\n",
"Pan-STARRS Pixel Analysis: Source Detection and Characterization,\n",
"Magnier, E. A., et al. 2016, arXiv:1612.05244\n",
"Pan-STARRS Photometric and Astrometric Calibration, Magnier, E. A., et\n",
"al. 2016, arXiv:1612.05242\n",
"The Pan-STARRS1 Database and Data Products, Flewelling, H. A., et al.\n",
"2016, arXiv:1612.05243\n",
"\n",
"Catalogue curator:\n",
"SSDC - ASI Space Science Data Center\n",
"https://www.ssdc.asi.it/\n",
"Num. columns: 26\n"
]
}
],
"source": [
"print(meta)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here are the columns."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"obj_name\n",
"obj_id\n",
"ra\n",
"dec\n",
"ra_error\n",
"dec_error\n",
"epoch_mean\n",
"g_mean_psf_mag\n",
"g_mean_psf_mag_error\n",
"g_flags\n",
"r_mean_psf_mag\n",
"r_mean_psf_mag_error\n",
"r_flags\n",
"i_mean_psf_mag\n",
"i_mean_psf_mag_error\n",
"i_flags\n",
"z_mean_psf_mag\n",
"z_mean_psf_mag_error\n",
"z_flags\n",
"y_mean_psf_mag\n",
"y_mean_psf_mag_error\n",
"y_flags\n",
"n_detections\n",
"zone_id\n",
"obj_info_flag\n",
"quality_flag\n"
]
}
],
"source": [
"for column in meta.columns:\n",
" print(column.name)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the [documentation for these variables]() .\n",
"\n",
"The ones we'll use are:\n",
"\n",
"* `obj_id`, which we will match up with `original_ext_source_id` in the best neighbor table.\n",
"\n",
"* `g_mean_psf_mag`, which contains mean magnitude from the `i` filter.\n",
"\n",
"* `i_mean_psf_mag`, which contains mean magnitude from the `i` filter.\n",
"\n",
"Here's a query that selects these variables and returns the first 5 rows."
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"query = \"\"\"SELECT \n",
"TOP 5\n",
"obj_id, g_mean_psf_mag, i_mean_psf_mag \n",
"FROM gaiadr2.panstarrs1_original_valid\n",
"\"\"\""
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
}
],
"source": [
"job = Gaia.launch_job_async(query=query)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"Table length=5\n",
"
\n",
"
obj_id
g_mean_psf_mag
i_mean_psf_mag
\n",
"
mag
\n",
"
int64
float64
float64
\n",
"
67130655389101425
--
20.3516006469727
\n",
"
67553305590067819
--
19.779899597168
\n",
"
67551423248967849
--
19.8889007568359
\n",
"
67132026238911331
--
20.9062995910645
\n",
"
67553513677687787
--
21.2831001281738
\n",
"
"
],
"text/plain": [
"
\n",
" obj_id g_mean_psf_mag i_mean_psf_mag \n",
" mag \n",
" int64 float64 float64 \n",
"----------------- -------------- ----------------\n",
"67130655389101425 -- 20.3516006469727\n",
"67553305590067819 -- 19.779899597168\n",
"67551423248967849 -- 19.8889007568359\n",
"67132026238911331 -- 20.9062995910645\n",
"67553513677687787 -- 21.2831001281738"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results = job.get_results()\n",
"results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following figure shows how these tables are related.\n",
"\n",
"* The orange circles and arrows represent the first `JOIN` operation, which takes each `source_id` in the Gaia table and finds the same value of `source_id` in the best neighbor table.\n",
"\n",
"* The blue circles and arrows represent the second `JOIN` operation, which takes each `original_ext_source_id` in the Gaia table and finds the same value of `obj_id` in the best neighbor table.\n",
"\n",
"There's no guarantee that the corresponding rows of these tables are in the same order, so the `JOIN` operation involves some searching.\n",
"However, ADQL/SQL databases are implemented in a way that makes this kind of source efficient.\n",
"If you are curious, you can [read more about it](https://chartio.com/learn/databases/how-does-indexing-work/)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Joining tables\n",
"\n",
"Now let's get to the details of performing a `JOIN` operation.\n",
"As a starting place, let's go all the way back to the cone search from Lesson 2."
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"query_cone = \"\"\"SELECT \n",
"TOP 10 \n",
"source_id\n",
"FROM gaiadr2.gaia_source\n",
"WHERE 1=CONTAINS(\n",
" POINT(ra, dec),\n",
" CIRCLE(88.8, 7.4, 0.08333333))\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And let's run it, to make sure we have a working query to build on."
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
}
],
"source": [
"from astroquery.gaia import Gaia\n",
"\n",
"job = Gaia.launch_job_async(query=query_cone)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"Table length=10\n",
"
\n",
" source_id ra ... pmdec \n",
" deg ... mas / yr \n",
" int64 float64 ... float64 \n",
"------------------- ----------------- ... -------------------\n",
"3322773965056065536 88.78178020183375 ... -2.5057036964736907\n",
"3322773758899157120 88.83227057144585 ... --\n",
"3322774068134271104 88.8206092188033 ... -1.5260889445858044\n",
"3322773930696320512 88.80843339290348 ... -0.9292104395445717\n",
"3322774377374425728 88.86806108182265 ... -3.8676624830902435\n",
"3322773724537891456 88.81308602813434 ... -33.078133430952086\n",
"3322773724537891328 88.81570329208743 ... 0.3110526931576576\n",
"3322773930696321792 88.8050736770331 ... 1.0772755505138008\n",
"3322773724537890944 88.81241651540533 ... -6.393939291541333\n",
" ... ... ... ...\n",
"3322962118983356032 88.76109637722949 ... --\n",
"3322963527732585984 88.78813701704823 ... -2.46251296961979\n",
"3322961775385969024 88.79723215862369 ... -6.605711792572964\n",
"3322962084625312512 88.78286756313868 ... 1.3495903680571226\n",
"3322962939322692608 88.73289357818679 ... 1.002126813991455\n",
"3322963768250760576 88.7592444035961 ... --\n",
"3322963459013111808 88.80348931842845 ... -3.3780655466364626\n",
"3322963355935626368 88.75528507586058 ... --\n",
"3322963287216149888 88.7658164932195 ... -0.5046963243400879\n",
"3322962015904143872 88.74740822271643 ... 0.5565841272341593"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results = job.get_results()\n",
"results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adding the best neighbor table\n",
"\n",
"Now we're ready for the first join.\n",
"The join operation requires two clauses:\n",
"\n",
"* `JOIN` specifies the name of the table we want to join with, and\n",
"\n",
"* `ON` specifies how we'll match up rows between the tables.\n",
"\n",
"In this example, we join with `gaiadr2.panstarrs1_best_neighbour AS best`, which means we can refer to the best neighbor table with the abbreviated name `best`.\n",
"\n",
"And the `ON` clause indicates that we'll match up the `source_id` column from the Gaia table with the `source_id` column from the best neighbor table. "
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"query_base = \"\"\"SELECT \n",
"{columns}\n",
"FROM gaiadr2.gaia_source AS gaia\n",
"JOIN gaiadr2.panstarrs1_best_neighbour AS best\n",
" ON gaia.source_id = best.source_id\n",
"WHERE 1=CONTAINS(\n",
" POINT(gaia.ra, gaia.dec),\n",
" CIRCLE(88.8, 7.4, 0.08333333))\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**SQL detail:** In this example, the `ON` column has the same name in both tables, so we could replace the `ON` clause with a simpler [`USING` clause](https://docs.oracle.com/javadb/10.8.3.0/ref/rrefsqljusing.html):\n",
"\n",
"```\n",
"USING(source_id)\n",
"```\n",
"\n",
"Now that there's more than one table involved, we can't use simple column names any more; we have to use **qualified column names**.\n",
"In other words, we have to specify which table each column is in.\n",
"Here's the complete query, including the columns we want from the Gaia and best neighbor tables."
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SELECT \n",
"gaia.source_id, gaia.ra, gaia.dec, gaia.pmra, gaia.pmdec, best.best_neighbour_multiplicity, best.number_of_mates\n",
"FROM gaiadr2.gaia_source AS gaia\n",
"JOIN gaiadr2.panstarrs1_best_neighbour AS best\n",
" ON gaia.source_id = best.source_id\n",
"WHERE 1=CONTAINS(\n",
" POINT(gaia.ra, gaia.dec),\n",
" CIRCLE(88.8, 7.4, 0.08333333))\n",
"\n"
]
}
],
"source": [
"column_list = ['gaia.source_id',\n",
" 'gaia.ra',\n",
" 'gaia.dec',\n",
" 'gaia.pmra',\n",
" 'gaia.pmdec',\n",
" 'best.best_neighbour_multiplicity',\n",
" 'best.number_of_mates',\n",
" ]\n",
"columns = ', '.join(column_list)\n",
"\n",
"query = query_base.format(columns=columns)\n",
"print(query)"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
}
],
"source": [
"job = Gaia.launch_job_async(query=query)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"Table length=490\n",
"
\n",
"
source_id
ra
dec
pmra
pmdec
best_neighbour_multiplicity
number_of_mates
\n",
"
deg
deg
mas / yr
mas / yr
\n",
"
int64
float64
float64
float64
float64
int16
int16
\n",
"
3322773965056065536
88.78178020183375
7.334936530583141
0.2980633722108194
-2.5057036964736907
1
0
\n",
"
3322774068134271104
88.8206092188033
7.353158142762173
-1.1065462654445488
-1.5260889445858044
1
0
\n",
"
3322773930696320512
88.80843339290348
7.334853162299928
2.6074384482375215
-0.9292104395445717
1
0
\n",
"
3322774377374425728
88.86806108182265
7.371287731275939
3.9555477866915383
-3.8676624830902435
1
0
\n",
"
3322773724537891456
88.81308602813434
7.32488574492059
51.34995462741039
-33.078133430952086
1
0
\n",
"
3322773724537891328
88.81570329208743
7.3223019772324855
1.9389988498951845
0.3110526931576576
1
0
\n",
"
3322773930696321792
88.8050736770331
7.332371472206583
2.264014834476311
1.0772755505138008
1
0
\n",
"
3322773724537890944
88.81241651540533
7.327864052479726
-0.36003627434304625
-6.393939291541333
1
0
\n",
"
3322773930696322176
88.80128682574824
7.334292036448643
--
--
1
0
\n",
"
...
...
...
...
...
...
...
\n",
"
3322962359501481088
88.85037722908271
7.402162717053584
2.058216493648542
-2.249255322558584
1
0
\n",
"
3322962393861228544
88.82108234976155
7.4044425496203
-0.916760881643629
-1.1113319053861441
1
0
\n",
"
3322955831151254912
88.74620347799508
7.342728619145855
0.1559833902071379
-1.750598455959734
1
0
\n",
"
3322962118983356032
88.76109637722949
7.380564308268047
--
--
1
0
\n",
"
3322963527732585984
88.78813701704823
7.456696889759524
1.1363354614104264
-2.46251296961979
1
0
\n",
"
3322961775385969024
88.79723215862369
7.359756552906535
2.121021366548921
-6.605711792572964
1
0
\n",
"
3322962084625312512
88.78286756313868
7.384598632215225
-0.09350717810996487
1.3495903680571226
1
0
\n",
"
3322962939322692608
88.73289357818679
7.407688975612043
-0.11002934783569704
1.002126813991455
1
0
\n",
"
3322963459013111808
88.80348931842845
7.438699901204871
0.800833828337078
-3.3780655466364626
1
0
\n",
"
3322962015904143872
88.74740822271643
7.387057037713974
-0.7201178533250112
0.5565841272341593
1
0
\n",
"
"
],
"text/plain": [
"
\n",
" source_id ra ... number_of_mates\n",
" deg ... \n",
" int64 float64 ... int16 \n",
"------------------- ----------------- ... ---------------\n",
"3322773965056065536 88.78178020183375 ... 0\n",
"3322774068134271104 88.8206092188033 ... 0\n",
"3322773930696320512 88.80843339290348 ... 0\n",
"3322774377374425728 88.86806108182265 ... 0\n",
"3322773724537891456 88.81308602813434 ... 0\n",
"3322773724537891328 88.81570329208743 ... 0\n",
"3322773930696321792 88.8050736770331 ... 0\n",
"3322773724537890944 88.81241651540533 ... 0\n",
"3322773930696322176 88.80128682574824 ... 0\n",
" ... ... ... ...\n",
"3322962359501481088 88.85037722908271 ... 0\n",
"3322962393861228544 88.82108234976155 ... 0\n",
"3322955831151254912 88.74620347799508 ... 0\n",
"3322962118983356032 88.76109637722949 ... 0\n",
"3322963527732585984 88.78813701704823 ... 0\n",
"3322961775385969024 88.79723215862369 ... 0\n",
"3322962084625312512 88.78286756313868 ... 0\n",
"3322962939322692608 88.73289357818679 ... 0\n",
"3322963459013111808 88.80348931842845 ... 0\n",
"3322962015904143872 88.74740822271643 ... 0"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results = job.get_results()\n",
"results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that this result has fewer rows than the previous result.\n",
"That's because there are sources in the Gaia table with no corresponding source in the Pan-STARRS table.\n",
"\n",
"By default, the result of the join only includes rows where the same `source_id` appears in both tables.\n",
"This default is called an \"inner\" join because the results include only the intersection of the two tables.\n",
"[You can read about the other kinds of join here](https://www.geeksforgeeks.org/sql-join-set-1-inner-left-right-and-full-joins/)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adding the Pan-STARRS table\n",
"\n",
"### Exercise\n",
"\n",
"Now we're ready to bring in the Pan-STARRS table. Starting with the previous query, add a second `JOIN` clause that joins with `gaiadr2.panstarrs1_original_valid`, gives it the abbreviated name `ps`, and matches `original_ext_source_id` from the best neighbor table with `obj_id` from the Pan-STARRS table.\n",
"\n",
"Add `g_mean_psf_mag` and `i_mean_psf_mag` to the column list, and run the query.\n",
"The result should contain 490 rows and 9 columns."
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"tags": [
"hide-cell"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SELECT \n",
"gaia.source_id, gaia.ra, gaia.dec, gaia.pmra, gaia.pmdec, best.best_neighbour_multiplicity, best.number_of_mates, ps.g_mean_psf_mag, ps.i_mean_psf_mag\n",
"FROM gaiadr2.gaia_source as gaia\n",
"JOIN gaiadr2.panstarrs1_best_neighbour as best\n",
" ON gaia.source_id = best.source_id\n",
"JOIN gaiadr2.panstarrs1_original_valid as ps\n",
" ON best.original_ext_source_id = ps.obj_id\n",
"WHERE 1=CONTAINS(\n",
" POINT(gaia.ra, gaia.dec),\n",
" CIRCLE(88.8, 7.4, 0.08333333))\n",
"\n",
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
},
{
"data": {
"text/html": [
"Table length=490\n",
"
\n",
" source_id ra ... pmdec \n",
" deg ... mas / yr \n",
" int64 float64 ... float64 \n",
"------------------ ------------------ ... -------------------\n",
"635559124339440000 137.58671691646745 ... -12.490481778113859\n",
"635860218726658176 138.5187065217173 ... -11.346409129876392\n",
"635674126383965568 138.8428741026386 ... -12.702779525389634\n",
"635535454774983040 137.8377518255436 ... -14.492308604905652\n",
"635497276810313600 138.0445160213759 ... -12.291499169815987\n",
"635614168640132864 139.59219748145836 ... -13.708904908478631\n",
"635821843194387840 139.88094034815086 ... -12.55978220563274\n",
"635551706931167104 138.04665586038192 ... -12.224246333795001\n",
"635518889086133376 137.2374229207837 ... -12.710314902969365\n",
" ... ... ... ...\n",
"612282738058264960 134.0445768189235 ... -13.651740929272187\n",
"612485911486166656 134.96582769047063 ... -11.998725329569156\n",
"612386332668697600 135.45701048323093 ... -12.436641304786672\n",
"612296172717818624 133.80060286960668 ... -12.50750861370402\n",
"612250375480101760 134.64754712466774 ... -14.061676353845487\n",
"612394926899159168 135.51997060013844 ... -13.526821099431533\n",
"612288854091187712 134.07970733489358 ... -11.162471664228455\n",
"612428870024913152 134.8384242853297 ... -14.247379430659198\n",
"612256418500423168 134.90752972739924 ... -12.145212331165776\n",
"612429144902815104 134.77293979509543 ... -13.962312685889454"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results = job.get_results()\n",
"results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"Create a new query base called `query7_base` that combines the `WHERE` clauses from the previous query with the `JOIN` clauses for the best neighbor and Pan-STARRS tables.\n",
"Format the query base using the column names in `column_list`, and call the result `query7`.\n",
"\n",
"Hint: Make sure you use qualified column names everywhere!\n",
"\n",
"Run your query and download the results. The table you get should have 3725 rows and 9 columns."
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"tags": [
"hide-cell"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"SELECT \n",
"gaia.source_id, gaia.ra, gaia.dec, gaia.pmra, gaia.pmdec, best.best_neighbour_multiplicity, best.number_of_mates, ps.g_mean_psf_mag, ps.i_mean_psf_mag\n",
"FROM gaiadr2.gaia_source as gaia\n",
"JOIN gaiadr2.panstarrs1_best_neighbour as best\n",
" ON gaia.source_id = best.source_id\n",
"JOIN gaiadr2.panstarrs1_original_valid as ps\n",
" ON best.original_ext_source_id = ps.obj_id\n",
"WHERE parallax < 1\n",
" AND bp_rp BETWEEN -0.75 AND 2 \n",
" AND 1 = CONTAINS(POINT(gaia.ra, gaia.dec), \n",
" POLYGON(135.306, 8.39862, 126.51, 13.4449, 163.017, 54.2424, 172.933, 46.4726, 135.306, 8.39862))\n",
" AND 1 = CONTAINS(POINT(gaia.pmra, gaia.pmdec),\n",
" POLYGON( -4.05037121,-14.75623261, -3.41981085,-14.72365546, -3.03521988,-14.44357135, -2.26847919,-13.7140236 , -2.61172203,-13.24797471, -2.73471401,-13.09054471, -3.19923146,-12.5942653 , -3.34082546,-12.47611926, -5.67489413,-11.16083338, -5.95159272,-11.10547884, -6.42394023,-11.05981295, -7.09631023,-11.95187806, -7.30641519,-12.24559977, -7.04016696,-12.88580702, -6.00347705,-13.75912098, -4.42442296,-14.74641176))\n",
"\n",
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
},
{
"data": {
"text/html": [
"Table length=3725\n",
"
\n",
"
source_id
ra
dec
pmra
pmdec
best_neighbour_multiplicity
number_of_mates
g_mean_psf_mag
i_mean_psf_mag
\n",
"
deg
deg
mas / yr
mas / yr
mag
\n",
"
int64
float64
float64
float64
float64
int16
int16
float64
float64
\n",
"
635860218726658176
138.5187065217173
19.09233926905897
-5.941679495793577
-11.346409129876392
1
0
17.8978004455566
17.5174007415771
\n",
"
635674126383965568
138.8428741026386
19.031798198627634
-3.8970011609340207
-12.702779525389634
1
0
19.2873001098633
17.6781005859375
\n",
"
635535454774983040
137.8377518255436
18.864006786112604
-4.335040664412791
-14.492308604905652
1
0
16.9237995147705
16.478099822998
\n",
"
635497276810313600
138.0445160213759
19.00947118796605
-7.1729306406216615
-12.291499169815987
1
0
19.9242000579834
18.3339996337891
\n",
"
635614168640132864
139.59219748145836
18.807955539071433
-3.309602916796381
-13.708904908478631
1
0
16.1515998840332
14.6662998199463
\n",
"
635598607974369792
139.20920023089508
18.624132868942702
-6.124445176881091
-12.833824027100611
1
0
16.5223999023438
16.1375007629395
\n",
"
635737661835496576
139.93327552473934
19.167962454651423
-7.119403303682826
-12.687947497633793
1
0
14.5032997131348
13.9849004745483
\n",
"
635850945892748672
139.86542888472115
20.011312663154804
-3.786655365804428
-14.28415600718206
1
0
16.5174999237061
16.0450000762939
\n",
"
635600532119713664
139.22869949616816
18.685939084485494
-3.9742788217925122
-12.342426623384245
1
0
20.4505996704102
19.5177001953125
\n",
"
...
...
...
...
...
...
...
...
...
\n",
"
612241781249124608
134.3755835065194
18.129179169751275
-2.831807894848964
-13.902118573613597
1
0
20.2343997955322
18.6518001556396
\n",
"
612332147361443072
134.14584721363653
18.45685585044513
-6.234287981021865
-11.500464195695072
1
0
21.3848991394043
20.3076000213623
\n",
"
612426744016802432
134.68522805061076
18.77090626983678
-3.7691372464459554
-12.889167493118862
1
0
17.8281002044678
17.4281005859375
\n",
"
612331739340341760
134.12176196902254
18.42768872157865
-3.9894012386388735
-12.60504410507441
1
0
21.8656997680664
19.5223007202148
\n",
"
612282738058264960
134.0445768189235
18.11915820167003
-2.5972485319419127
-13.651740929272187
1
0
22.5151996612549
19.9743995666504
\n",
"
612386332668697600
135.45701048323093
18.63266345155342
-5.07684899854408
-12.436641304786672
1
0
19.3792991638184
17.9923000335693
\n",
"
612296172717818624
133.80060286960668
18.08186533343457
-6.112792578821885
-12.50750861370402
1
0
17.4944000244141
16.926700592041
\n",
"
612250375480101760
134.64754712466774
18.122419425065015
-2.8969262278467127
-14.061676353845487
1
0
15.3330001831055
14.6280002593994
\n",
"
612394926899159168
135.51997060013844
18.817675531233004
-3.9968965218753763
-13.526821099431533
1
0
16.4414005279541
15.8212003707886
\n",
"
612256418500423168
134.90752972739924
18.280596648172743
-6.109836304219565
-12.145212331165776
1
0
20.8715991973877
19.9612007141113
\n",
"
"
],
"text/plain": [
"
\n",
" source_id ra ... g_mean_psf_mag i_mean_psf_mag \n",
" deg ... mag \n",
" int64 float64 ... float64 float64 \n",
"------------------ ------------------ ... ---------------- ----------------\n",
"635860218726658176 138.5187065217173 ... 17.8978004455566 17.5174007415771\n",
"635674126383965568 138.8428741026386 ... 19.2873001098633 17.6781005859375\n",
"635535454774983040 137.8377518255436 ... 16.9237995147705 16.478099822998\n",
"635497276810313600 138.0445160213759 ... 19.9242000579834 18.3339996337891\n",
"635614168640132864 139.59219748145836 ... 16.1515998840332 14.6662998199463\n",
"635598607974369792 139.20920023089508 ... 16.5223999023438 16.1375007629395\n",
"635737661835496576 139.93327552473934 ... 14.5032997131348 13.9849004745483\n",
"635850945892748672 139.86542888472115 ... 16.5174999237061 16.0450000762939\n",
"635600532119713664 139.22869949616816 ... 20.4505996704102 19.5177001953125\n",
" ... ... ... ... ...\n",
"612241781249124608 134.3755835065194 ... 20.2343997955322 18.6518001556396\n",
"612332147361443072 134.14584721363653 ... 21.3848991394043 20.3076000213623\n",
"612426744016802432 134.68522805061076 ... 17.8281002044678 17.4281005859375\n",
"612331739340341760 134.12176196902254 ... 21.8656997680664 19.5223007202148\n",
"612282738058264960 134.0445768189235 ... 22.5151996612549 19.9743995666504\n",
"612386332668697600 135.45701048323093 ... 19.3792991638184 17.9923000335693\n",
"612296172717818624 133.80060286960668 ... 17.4944000244141 16.926700592041\n",
"612250375480101760 134.64754712466774 ... 15.3330001831055 14.6280002593994\n",
"612394926899159168 135.51997060013844 ... 16.4414005279541 15.8212003707886\n",
"612256418500423168 134.90752972739924 ... 20.8715991973877 19.9612007141113"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Solution\n",
"\n",
"query7_base = \"\"\"\n",
"SELECT \n",
"{columns}\n",
"FROM gaiadr2.gaia_source as gaia\n",
"JOIN gaiadr2.panstarrs1_best_neighbour as best\n",
" ON gaia.source_id = best.source_id\n",
"JOIN gaiadr2.panstarrs1_original_valid as ps\n",
" ON best.original_ext_source_id = ps.obj_id\n",
"WHERE parallax < 1\n",
" AND bp_rp BETWEEN -0.75 AND 2 \n",
" AND 1 = CONTAINS(POINT(gaia.ra, gaia.dec), \n",
" POLYGON({point_list}))\n",
" AND 1 = CONTAINS(POINT(gaia.pmra, gaia.pmdec),\n",
" POLYGON({pm_point_list}))\n",
"\"\"\"\n",
"\n",
"columns = ', '.join(column_list)\n",
"\n",
"query7 = query7_base.format(columns=columns,\n",
" point_list=point_series['point_list'],\n",
" pm_point_list=point_series['pm_point_list'])\n",
"print(query7)\n",
"\n",
"\n",
"job = Gaia.launch_job_async(query=query7)\n",
"results = job.get_results()\n",
"results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Checking the match\n",
"\n",
"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": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<MaskedColumn name='best_neighbour_multiplicity' dtype='int16' description='Number of neighbours with same probability as best neighbour' length=3725>\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
\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
\n",
"
"
],
"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": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results['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": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"count 3725.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": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pandas as pd\n",
"\n",
"multiplicity = pd.Series(results['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": 61,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"count 3725.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": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mates = pd.Series(results['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. But we are more interested in the final match, using both criteria."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Transforming coordinates\n",
"\n",
"Here's the function we've used to transform the results from ICRS to GD-1 coordinates."
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [],
"source": [
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord\n",
"from gala.coordinates import GD1Koposov10\n",
"from gala.coordinates import reflex_correct\n",
"\n",
"def make_dataframe(table):\n",
" \"\"\"Transform coordinates from ICRS to GD-1 frame.\n",
" \n",
" table: Astropy Table\n",
" \n",
" returns: Pandas DataFrame\n",
" \"\"\"\n",
" skycoord = SkyCoord(\n",
" ra=table['ra'], \n",
" dec=table['dec'],\n",
" pm_ra_cosdec=table['pmra'],\n",
" pm_dec=table['pmdec'], \n",
" distance=8*u.kpc, \n",
" radial_velocity=0*u.km/u.s)\n",
"\n",
" gd1_frame = GD1Koposov10()\n",
" transformed = skycoord.transform_to(gd1_frame)\n",
" skycoord_gd1 = reflex_correct(transformed)\n",
"\n",
" df = table.to_pandas()\n",
" df['phi1'] = skycoord_gd1.phi1\n",
" df['phi2'] = skycoord_gd1.phi2\n",
" df['pm_phi1'] = skycoord_gd1.pm_phi1_cosphi2\n",
" df['pm_phi2'] = skycoord_gd1.pm_phi2\n",
" return df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now can transform the result from the last query."
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [],
"source": [
"candidate_df = make_dataframe(results)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And see how it looks."
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEJCAYAAACdePCvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABmGklEQVR4nO29fZSdx1kn+Kt2AgHMx4Cl28oGJbCQDFhSyBe2d3GTBDnDYsvZQMJJgFlLChhwd9u7rKNWyMJykI0TyEyOlDADmZEUzhliYoVPywu2ukkCS5wQSfdeqU0CAyHBfe/t2/IZhoEwC+q+tX/cW1fV1U9VPfXx3tvqrt8575H6vm9VPfX1fNVTVUJKiYKCgoKC7YeJcRNQUFBQUDAeFAFQUFBQsE1RBEBBQUHBNkURAAUFBQXbFEUAFBQUFGxTFAFQUFBQsE0xVgEghPg6IcRHhRCfE0J8Vghx2zjpKSgoKNhOeN6Yyz8O4A+klG8WQnwZgK8cMz0FBQUF2wZiXBvBhBBfA6AJ4Jslk4ibbrpJvuQlL6mUroKCgoKthgsXLjwnpdxh/j5OC+CbAVwBcFoI8XIAFwA8IKX8kv6REOJeAPcCwO7du3H+/PmRE1pQUFBwPUMI8UXq93GuATwPwCsB/Hsp5SsAfAnAUfMjKeUHpZSvllK+eseODQKsoKCgoCAS4xQASwCWpJSfHvz9UfQFQkFBQUHBCDA2ASClXAbwrBDiZYOfvgfAn42LnoKCgoLthnFHAc0C+PVBBNDnARwaMz0FBQUF2wZjFQBSygaAV4+ThoKCgoLtirITuKCgoGCbogiAgoKCgm2KIgAKCiqGlBLLy8sot+8VbDYUAVCwLdDr9XDp0iX0er2Rl93tdnHy5ElcvnyZFAI5BUSOvIrA2j4oAqBgrFDMptfrsZhOKHOSUqLT6eBjH/sYHnzwQSwuLkbnxaXF/L1Wq+HAgQM4e/Ysut3uhm8vX76MU6dObXgXg263O8wrtn56HgVbG0UAFKzDqLU/xWwWFxdZTCeUOXW7XZw4cQJPPvkk7r33Xtx8883ReXFpMX8XQmDv3r04fPgwarXahm8ff/xx3HXXXRvexaBWq+HQoUPDftTp8PWter9z585hHsUK2OJQnXw9PK961atkQT70ej3Z6XRkr9cb/tbpdOTDDz8sO50O6/uUsvTf19bWnO/V76E09Ho92W63ZaPR2FCvlPq4aNfzdZWhaGu320H16XQ6cnV1VTYaDdlqtTakVX3YbrfXle3qWzNdo9GQR48ele12O6RJrPTGtjE3jxzlbGUAOC8JnlosgE0GOUINnNJga7XaUFPVaZEeV4WPbpu2LITA5OQkJiYmMDk5CSGENZ2UEt1ud6gpc9pJCIFdu3Zh3759GzRwVbZZJhcq/crKygaNX+Wr6Ff+f72dut0uTp8+DSEEmwaV3yc+8Qncf//9ePe7372hTVUfTk5OYnJyEgCGmj1lhZjpAODRRx/Fl770JfK7EJj9HjO+9Txsazkca04O3IGdTqdYNgqUVNisz/VsAXA1FF1LW1tbk81mU66trVVSlk87bTab8qGHHpKdTkd2Oh350EMPyWazSX7fbrfl0aNHZavVYmnyXOhtoLeNT5utAj4rxtWOx44dk81mU7bb7SHdMW2iWwD1el3W63Xv+AhpK2WVtFqtoWWSYoWZ3/rGkS+PZrMp77jjDtlsNp3lUOh0OvLo0aPy6NGj0e1/vQIWC2DsTD3kuZ4FgM0sN8EZ7NyyYpljr9eTjUZDzs3NyUaj4XVlSHlNADQajeAJ7oJeF9O1YnOf6EySK0C5DIRqV19a1Zeq/026YwU2t59DGB2Vp+pb5Q6ivgmpA+Ve4qaPVYpUGXrbc+dkKI2bEUUAjBlq8OgaoA9VWwA2dDodOTc3J2dmZtg+YN0frphdDu3cVRcbA1S/LywsyDvuuEM2Go0NecQwUxstoWnN77mM3Pwu1oJwrTtQ6xrtdnuoDKytrQ3T69+EaPa6QFFpW63WOqEwCmYbOie5/cxp31ELkSIARgSONriZtAjKvDddANR3LqRoaSF0cy2AVqu1wYowJ32KO8a2gG3LN8SlYlo9sa4XBdMNYoPpilxYWBgyeMoNp1s5IRaJysO0CmyWH5WXPg44fcJtq5jvfO07DtellEUAVA6TsXDNSldeoxASXG00ZOCOYpCHlEGtI6T0TwgNqW0Rm95mLegavMuVZjJpxfx1xmoyWx/TdQlDM1/KVUNZF+12W87Ozsq5ubl1Qim3FcqN1oqxsFx55eIDRQBUDHMRlGNW2jTlVN9uiBVi0mCbpCHhipyBmzq4OZq1AlebDC2T0/6pFmGMNUX1l00gqHHrsgh8NNqEqtkOXKZs5qfaQAUk6FD01+v1dW6pUKvMRw/HagrBKJQHHUUAaKhCwzYXyjhlqEVe3UcdwmxjtXX9fehAjDGxXXnmWiw26dSRq7+5NOe0lmKYAJXGbAOdycYuSpt5m0qPXoZi4FTb2QSrnp9NEKpv1RqCWqfgCCxX+5hlhO7X8GEUSpKOIgA0pEjWWK2bghrUuo+aS9va2tq6jUChWie1kMeps84EuYt+Nm3dptWFIKTeqej1rkVHuRbHQ+hIsRBSxqLLmuGMQY4FpjNyF/O3jQMbTVSft1qtoSuIsjR848Ql3KseV6NAEQAaUjo0p1lG0cM1+c0Q0U4nLL46th66BWBOtFAXTA43kE+IcBhVCExLj0tnFQyE0+fcsqnIHFeeet+7hImPBq5VZRMGlCVDKTW+8e4aSz7LkipvswmNIgAyoeoIlxALwPTdhyx85RigNrcCV5PLVZ6LeVDl6Ywz1JXFbbeUOnMZi860bHXhlq0EwNLSkvWICT3PY8eOyYWFhQ0b9ELKlNI+n7iWD9fa4bhwQi0qffxRY4yyQnK7krgoAiATqrAAdNgGiYspmN+4JkaIdhZSdmheMd/o2h6HHlteurDMGS1C0UppyC6YtJkavjn+dCGQssFKfddoNOTU1JScmZnxastcC4DTTqqOKm3IfhluGbn7OMQC6PX67sPp6emg/TW5UASAAyGSOUZzDk1j01rV5qaUTStcbY3SYG3MModmw9GYlba4tLSUfFAZV6hR9HHehzJFnfHpGj3HzaZ+MzdTxdTDXFuKzYcLtcdAPypDCTNzL0pKGaaVMWoXTafT31z5lre8RR48eFC2Wq2RlKtgEwDb8jA4KdcfSKWODD5x4oT3aGD9oC8p6cOlqPxDjh3WD2TTf7vrrrvwyU9+kjw6WJXZG5yrbzvS18ybKsukWX1z880346677sLOnTs3fMttPxuo8kyaVlZWcPbsWTz33HPr6uxqextch9BRefj6kDoCWs+Xm15KiQMHDmDPnj2YmJjYcIy0EAK1Wm14OJ5e1q5du/DAAw8MD4AzIWX/QL+TJ09a6ZiYmMDLX/5yvPCFLyQPqFNtAyDpID293k899RRuu+02nD17FgBw4MABPPHEE3juuedw+vRp9piy9b0aNysrK+vKdR1Sxx1HnO/UHHzrW9+KHTt24MYbb0xut2ygpMJmfaqKAur1erLVag1DyELyoeKDqfyVtpF6lolNa1F+bbVjs9PZGFEU6org+pN9bqvY8mzfKK2dchNwzH1fWaarieNrVm1g01pDLICU+HDfGEkJvdXdTbncKdQhgrHWma3tqPFptpPZNtxx5HIhmnRRIbejAooL6BooX3mIz1EfoKGML/aANw5Nps/YN8hDESq8cvheXW3pMu1dzMPXDrF+aJtCEIIYQcgR0ty8XVBKxvz8vNdFxAW3vq55Ztt9z3Er6nnpwo1Dl2oP2ziiFInQtsjlqioCwILQTtLTxEzyKqOIOFomJ47dhtB65xi8rjL1dxzrQKeLo8mGrm3kWAsJAdU2Vfq2Vbvl3hXrg02wUtaDmc63kUyHre1sjNl3FpXroiBbPWMsWg6KADBAdRJ38lQxyXLmSVk4Kt+QOHYz7agXznxlUhqefuwyx+Xiq0vqBKyaIY+jP1KEXAzNtjJtY5ma275xz3GdcRmzTankWp7FAqhAANjcIVVIW6o8H3KWredl5htCVyxN4xQUNrdcar4xLgopR3/6Yw4lpkpFaBRzjCrDJwBc8ySUMdtciFzLswpsewHAYYS5GJfq6FB3SS7G5Ru8ofmEWkmjZnrjLt/mopAybUyFpNWZjumXpvLhutZciGnnUVjZMXOboyCGllsljwnFphUAAG4AUAdw1vdtTgugSnQ64dfeqXRVM65UrU3/vy2vcQ1ybvk5BX0VVodCiPWmux1M5YPqM4pmSuBz6l9FP+s0h66b5aDLNbZzWHscYZETm1kA/BSAD1ctAGIRYyrHdm5VDFVPn6rZUHnFHKkwTnCEWGg+uUExGld5LqvP1v9mfjmt5FjrkSpHRc4tLCyw0lLt5JuTXNra7bacmZlJtvb0ueOal7mCRjalAADwIgALAF4/LgHgGwjjMpWrSG/GIecQLHroaWg0lS/vEG00Ju9YgUjlk4seBYquHMyAIyhi9lnY5o3tBjEqja099d3CXO3aJzx9f9vyabfb6+7K5rSz7b2aO679AbnCxjerAPgogFcBeK1NAAC4F8B5AOd3796d1AgUfAPBNelSNXsuclkQtvPaucyPyl+1z+rqqpVpxIBiICGIbbMqTHFfnhyN1fZdjvJt9FDC3JaXjbm7LAAOE6bK5SgHocLG1t762g4lEHxt6Oor27zUsWUtAAB3Afh3g/9bBYD+bDYLoGrETlwbrVxmZ/udiqTguA048GmjoW6I2H7L0d+hY4jbD7Ftq5efEg+vvwuNcXeVw71S0qxPVWdT6XTp+YSMjZD2oL7NrYhsRgHwCIAlAF8AsAzgHwH8J1eazbYGUHU5ocw1ddDZylWgBEDVZVLguCFyCaIY5BKKMRaab2d6rEvB1FZdh9SNSnGxCYxUIW67GyHn/PMht+K56QTAOiLGZAGMirlzwHEDuExzW5qQ8nx5pLaXq44hvn7F6Mxzd8ZprZm05Q7n5cB0WVCIdSmotl1aWpILCwvy6tWrVitHMVCfn1zBV88YQZgylnT6dfdMCJ1cGmzaf+6osiIACKQwjNzaAOd7Uwuj6ObWaRzCj9KmYtcMKGaXUqeU/tL/XyWz4uRV1VEUSnDU63W5f//+DRE5umBpt9vD6xlD3CU2JYCrJMS0paJ1enp6na9fX6Dl9m3MOKC+q0KR2dQCgPuM0gKI0UpSGEEOurlnuVcdV8353mfN+PrGdfJmDM2h/aV/H6P5ucrN5VbLBZ0htlqtdSfOKphtEBJtptL6Fvx9fRRyzIlCq9WSMzMz8uLFi+vopfZIhJ4LlWIBhFrFPhQB4IBiiFevXl0X0x7KoFIYgQ8ctwzXv0t9pw9wc+CFTiwOM/UJITMPm3C1tUvVmrjNAqDyClU0Qqy4lKMFuHXudK4dNb66umpdZ4jdIKUzPBeTtZWhfq/X6+yDDlWZ1CU6o9LKdTqocRQj0GzY1gLAN9AVQzxz5syQMfpcFC7J7fuNS5cO2wDUBwlXs6e+010qpibmO3WRW68QC8nMw8b0bfmMQjO2waQpVRjZxlWz2Vx3N29IniF06YLGdSMdpy982q1L8NvK0Meu8ttzhBo31DVU+7fB1wcmTZwb3rjY1gLAN9ApC0DBxki40tkVORPi/06lg5O/zQIw6TUnTYgmGWshubS/cTF6G3IzENVurVZrnYX60EMPyTNnzrAsAJurSe9zV/uqOaKubqT6ntMXio56vS5nZmY2XI3oY/i2tlXuQJfCZrqlcsb0c2DjBea8U3VQ/a322KSM8W0tAFKYRE4B4GOkMbT5BnGqBUJN8NgNZanMuiozPBYuS8e3EzUkf8UYGo3GOguVYsjcvEy6bOf8+4S2j0Hb6KjX63JqakrW63XW97bxa7rAbN93Ov07efUL2W3hniZybcay8QxbG6p51mg0kpW8bS0AUuAza0MiE0LMzljaOCaz+o6zy5ZTxqi0cJsrZLO4evTfc0QomfmbjChHnvo4jrGwuOPNhKpLq9WKtoJ14cURgCYj1Rmya266hFyIcLDNG9tmR9UfOVxBRQBEwjcBQjQ7V14pGmKMRqZcCL7JUxWDzZUvp924wjoUPgvAFaXEodW3IBpTl9R2DxUI3Pxsbkf9d87YVr+trq5uYMwms9b/tlm3VDm6sFhYWFgXUBGipHQ6neEajnLzrK6uykajIev1etaAkiIAEuDSeHN0js2VE5M3N43PtE5xG3HKyyHwbO1mQtcUcx/z4EJMHbm0xuRtIrSfuVZjTDk2t6JeHqftVTqTMevvzDJczN5Vn3a7vWERnsrTZYUrOtXi+sLCgpyampL33Xcfa3GfiyIAEqBry0pryOmLdrkSRunzdmlBqTRxLROKHkpoKGbEWQDNaQH46q8LpZjyfG6BnIqHqy42RsaNOgotxxZY4NscZvtdtwBcVgZ3Idil2VMH3nGPwTYPU6QsgBwoAkBDjJmqmL5+/HEOV4Jr81Yq4/JZFrYFQdfidBWuB5eJbBMayn9br9dHGhnky7vT8R/HEIucCoHPTcWxALjujpC+5/SdOh6au/8hVNApOvS54/qOUgi5feUTZLnGcBEAGjhaHEfap0I3AanNWxyN3AUbM7K5GUI0zFQ3kf6NTg+XmaiFsUajMbJNO5y6pGqVru9CtWEXzDGgaPZZLlS/5Wx708qj6t1sNuX+/fvlmTNnohZfOe/MuWMTmO02fT8Ad/5w5nYOYVAEgAaOFmdqwlVolZ3OtR2W1EDmaOQuhFoA+jvfwAzRqqi2Mye6z51j0qXaJPa46BzIIaA56SimaGvLUNr12HNzI6Bvv0eoth9Cl6kY6C7YVqsl5+fnh+MlxBKxCWeqbq1WSzYajXX+fVOhMkM7fUqUbRy72iuHQlMEgAGOVqBP7JBOCNXsYrTsqsDd9cvVqmzMPdT0TxWGVSCHgFZMmOuOocZhjnGkM0fbjWBV9gGllOi/6S5YkwauJWKziKnv1bg1w3BNd63Ztno+rjxD2rBYABUIAA5D90nzlLx9yK1dceuSa2exgq7FhdKth+mlxFvnQBXaroJiTLa7Zn10cPKPGY82Ic2xwkKhj5PQfTe+vjEFLccCUNa5Pm457UiFlrrypMqO5TsuFAFgIJeWniuNCa5Wwy3fpZnkoj2WSdq+0Q+tC2FiKe2UmmcMlPZNWQA5BE9ugWgKgpzhtVXcAR3TdxyhQgkTl8Jjy9NsP53eXApZEQCRqHLi22AbXNyJTJnqVWgVtnJD28qWbnV1dbhJhutyov7OQW9VbeZDiiJQNapk2jlh67uUPm236TsPFEMPiVAyrQKdriIAxiwARjXx9XJSJ7ipoaXkw41o0csNbStbOm475HJxxH7D+S5n2+QQ5uMSaLmQwzJKmWeKMavon9g+qWrcmCgCgMBmmgT6YMxFV+pAVL7pKmLaOXRSGmboxE9tS8qacn1nayff+1hXY+x9AKHWT66d6qEIUQ5y1Ikqm7PuQJW9mdx3RQAQcPnrpBytgMgxWFJgG8Apu1pTyra9D9XaQhgv9X+1o1Q/tIzqF98ita8vY9xR+l6IXBZAiCKQwxXla3/b2TxmiKarThQ4tOv8ISZWnyrDFTYaSycHRQAQ8GlQXO2vKlTp6zUH3qgGIocW1/tQl4xvU5Pe/5SgUSdIqrBY21EgscoEh07z25zHBJvwKQI+BSXU6tTb0jyd86GHHiIveEm1TDnav16XWOWHUgpMAeDy8XPp5KAIAAs47oNUX3osXbk6n0IoQ99M7jIXQvrMZNqUoDGP4lWMiQrl4ygT5ruQfjAVkiqsslwuM059zF20OjM0hYNJoz43QmmuUpnRmT7HLeTad5OTziIAEhDjY3b9HnKCper83Aw4VJPmTrRxC4oQq41DqxnXHXtEc+p4MPPgWjgh2nhMv7ksSV+eJsM324FrpYa6z1yb7lLQ660/2sXmqjJpUY/vCIwU2ATABAq8EEJgcnISQogN77rdLk6dOoVut8v6fXl5Gb/wC7+AEydObHino1ar4fDhw6jVas78YuGqkw69XA4NqXRKKbG8vNzXTiJQq9Vw6NCh4f9d9aPawCx/ZWUFZ8+excrKCrrdLk6fPg0hBJmvLb9Op4Pl5WWSHpXX6dOnrW2maAIwzF+18+LiorW9bX1BtXFov6k8lpeX16XT28CX5+TkJB544AFMTk6SbWu2py0/c664oMp57rnncOrUKVy+fDl6rFF5P/300zh69Cj27NkDIQSee+45vOMd78Di4uKG75eXl3H8+HGsrKysq0O328XJkyexuLg4rFPKnHCCkgrqAXAbgF8GcAnAFQB/A+D/ATAN4Gtdaat4xmUBuGDTUmzSm3t4VKqZmwsxFgA3uoJ6n8PsteXBod9lecX0Acdf7bMcXEdpxFgAHNcEp14+SyunBcKxvjjrFHq75T7WnRor6nhn6rRftZfAHBum+yvHnECoCwjA7wM4CeBuAC8E8DwANwJ4JYD/E8DHAdxtS1/FM24BwBnQvo0bVB7UJE9d6KoKMUzU9p5iHtyJHEKj+ntUJy+a+aWs5ZjrFLYyQmi2jcGQ9s29PuXrd87hhK5IJdtY87lvue1ifqfPAZdCwjmaIseYjBEAN9nehXyT8xn3GgDlp6OkemhkBjXJq5hkOcDRRjjanN6eJnOjykjRgnzaai5NleOzDmUoHD9wTo2e0776t7kEZgwD18ERINx7PEzm7UpnjmVqgT62jXIqI8ECYDM+oxQAeuNTgzDFzWAra9Rb6nNoOrbfOGl97o1YDZVDsw4u87PRa3tHKQPccRMitKj62RgXR+v1Qf82RuGx9S+15yQXE+z1Nl4Kz6FP9attfKh+UuGq3IhBTr1yuH4UbAIgahFYCHE5Jt31BH3BSS3QTE5ODhelbAtP3MVVKs3KygprIU7K/sJip9NRllgUqAVetbBny5eqH3cB0fxOCIE9e/bgwIED2LlzJ6TcuNipl6sWyEzaVDqKZinlsA9tfcJdROx2uzh79iwOHDiw4VvXO05ZZtuo7wB4gwyotqnVajhw4ADOnj27Lq2tr7jjltOePlA0CCGwa9cuTE5OrlvY9tHV6/XQbDbRbredc0EIgZ07dzrrpcaQXqYQAnv37iX7TDHRO++8E2fPnh2mM7+lxidnzoQsbkeDkgoDQr/f8vwAgCu2dNwHwDcC+BiAzwJ4BsADvjQ5LACuBpTbF5xCm4lc6wOUFhmz5yHFeuD6SjudjTH56l1I/H1VloRPq41pm5DfqbaznV+fsrvbLEctpl69ejXouG6bu6bZbMpjx46xjw9vNptyampKzszMkPs5zHJt9Q2xAl2eAQo+F90oeA0i1gCuAvgQgNPE8/e2dNwHwC4Arxz8/6sB/AWAb3elySEAbB3tGgBVdFCqK4WKmc6xZjBqwceZCLqZbXOpcHfgcidtLuQ0412wCVelKOjx5krIh9yra7opzUtjXFebUqDaxdeXFJSQq9fr5I5uW9twFT6XoPXNFbPNXALCdjNgrvkYIwAuANhjefesLV3sA+B3Adzh+maUFoD+e+5wMSnpgZSy4ziXVZALOQWJPpFsm4W49+Xa2jqnrznEkkxRBDi0qPZqtVpyZmZGzszMyFarte6aRU45JlM1Bak6M4ljAdjqmGKFuxQJl7DxzZUUCzLEqrAJ0FwKRIwAuB3Absu7V9vSxTwAXoL+HoOvId7dC+A8gPO7d+9OaoQYKOnsOi7ApXlyGJNPK+Vo9zktgBz5VKX5hlpwtgls9otvonMnvqsvbVo6lzmlMMRGoyHvu+8+OT09PezbECFjjosQazVWmIWMIc5cjLW2Y+kPSWs7THBsFsCoHvT3FlwA8P2+bzfDfQD6oHQNUN2cTQ2bHKV2n6usKlxJNteXydBDtUGKXp/G69PoKWuOKptqJxszCGGI1Dit1+vWYwlcUDRSazC+dL7xH8q4be+5bRuDGCuO+10V84RClAAA8DoAvzVYpH0GwEcBvNaVJuQB8HwATwL4Kc73VbqAXKDu+dTdEa7zwrkLbS66qMW80DxCBuk49x+46ORqyy5mwHVV2CwAl5vOFDy+s21ssIVW6nn6xhUlBGNdjLoAUS4kbjqfm0lXlCiBzR0POevrKof7nqI7JCw4N2wCwBoGKoS4E8ApAI8D+CEAP4z+MRCnhBDf5wgsYkH047pOAvislPLfpubHRcxZNYuLi3jwwQexuLg4DLlbXFwcnl1ChW+GhoO66FpZWcETTzyBiYkJZ36uPLj1VuF4u3btig7z0yFl2Nk+LjqpsDjubyrfZ555Zni2jwuq/wCsC02lwvwo2lUIoX7ejJmnmoSdQUjv2toams3msDwdUsphm5w+fXoYKtnpdHDp0iX0ej2Sfr0PY8MKVbqdO3fixhtvZI+LWq2Gt7/97di7d68zBPeuu+7C448/ToYjm2cNUXSpkFTzDCIA3vr6xqfqo0OHDqFWq234vtfrodvt4p577tEVW3IcSynxD//wD+vKGkmopwuUVBgQ+HEALyd+3wfgE7Z03AfAdwGQ6J8z1Bg83+dKM04LQNfATc0mVGul6EjV3n0a5yh8nlT6UA2nKp9rqEWmkOKHtrlAdHdSs9kcutwWFhbk1NSUnJ6etu6ONneaNhqNdYuHVBu46hDS3rZvueOTM75NdxN1THIOehRCtXuzH1T7LywskOGxrrsARglELAJ/LuZdlc84zwLSO951xKsJm0slhjG6XDMcVwhnElJb5kOElTnIYxl6TDpfm8b4YH03ffngEsyqrefn52Wr1SIPDvMJLpO+UF94TJvF5OE625/Kn8Msc7hPfGPC5wpcWlqSCwsL8urVq1bFJ1b5yIkYAXAh5l2VzzgFgN7xoYuktkkZ4mvvdNyLsxzm5vLJmlqpqfXo6Vz55NJyYia3j1lxhGQOOrh5KMboisWPtaByrD1xy/dp4q5zeLgWcgztsd+E1Nm2zqDnqfJrNpvB4bK5ECMA/iuA3yOexwH8rS1dlc+4TwOVMm6R1DYIQ90LqYuzXE3M58oI1ehitPkUl1WIa8DXrjlo57o/KJgaPpeeWMER6k7hWBeuyKkcbhwXXO2g+r7RaKwTwq65YMvDdaKn6kO12Stkw1wuxAiA73Y9tnRVPikCIGZApZr/sTSlDP6UCevSHkOYmgnKHM5Vtxit3kWfr7yU/Dj5+lw63LqF0h5r6egWX07Bl2p56WPZxqB1q9oWqcShw/YNNe59FoAuUEZxI9jImXjKkyIAdDOM26DNZnPk0prSrkOZ7bFjx8ht5Zy03L0LIaDM4Zi8qbQ5LA6b4Iul1adBuvI1x5ypYebSjF00h3yv7w3wadt6/j5XYQg91Ld6/i4Bz9lgGaow+YSPK892uy1nZmaGQinXPIyxAN4IYFr7+9MAPj943mxLV+WTagG4/K1Up4zDAqAYeKirKNbMdFkAOcAxl11mdNX7E8x2zuH+4QotBWrMpWrDIeCOAcq9w3EHqTooBu2K8jFpslkWVPu4LJMUUMqM6eLS6el0rq2XmecnUf3Zbl+7MXCsFgCAPwHwjdrfDQDfAGA3gAVbuiqf1DWAVH9lbtiYQ4oFIGUewVWFG8rFyHxmdBXnMZm0clxyoa4LM22MEMshmLmWktneruNJOIzb5s5Q6Tmbtcw25QhqV5ulCHadXttvpsWm2pJqU59wi6XXRIwA+Izx9we0/3/Klq7Kp8pF4JABGZJfiGZj/j4KrddGU8phaTZm7hJMnLZIcQnY6OBo16ZG59p96qKz00k7ZsPlyvAJGYpu1ccUQ1K+6larFWTB6OVTEXMmrZyQarNtQ+YF1WYxFpVKow7S0+mlrAK9vroQNOl20ZKTJ8UIgL90vPsr27sqn1FEAaUwbh2cQcbRJjmhn7lC/nS6XWatDzEWQCi41hzFKHQ/e2h/2yY7t79TBLpNOzQvEKfGjMmY9fYxo2D0NlQMb3V11cnk9HfUEdSUQDG/DZlTuqUSOldd/WB7x2XGZnk26yVkE2kOr0SMAPh1AD9G/P7jAB61pavyGXUUkA0xzN1lbrsGkzmBKN8w1z2SIpRytN+o+kCfrBQzTHGR2YSB7f+5QTEY3cesNiraol50hrK6uioXFhbWXcBiKhWK4ZlhwDoNpr/dJXQbjcY6LX5paWm4+U3lr+haXV3doEHrbezbQ+GCbfy4BKhN2fJZRCYfUHT73GvjXgPYCeCT6N/a9W8Gz8cBPA2gZktX5ZMjCihF+0zpGKr8EEZBRSRVYQFUiZzlczX33C402zhyuQFy0uzTvkO002azKffv3y8XFhbWMVVdSJgXwFA0cTb/ra6uyvn5ebm0tLSOVnXJT71eH9Z5YWFBfsu3fItcWFhwWsEcJUtKWuC7FB2fC83XJzbYhCslhHKvd0WHgQJ4PYDZwfN63/dVPrktgFCGFNrRMWXZBvW4dhDmANd8ptKk9FdqupC89LFhswY448fF8Hx1MZmXTXtWoPYcKGvCXB9wCQBO+5r1UmUvLS3JmZmZ4T0FUsoNFgD36BWbyyRHODe33X1rQJzFa1MA5Ri/234fAIVQScvtiBQJbqatQhsYNcyJ6bNY1tbW1l1XaAqQUJ/oKNrQNpl1NwXHTeSzAHx56HXVmS7lJnRp0Gabc92MNvrMeik6W62WnJ+fl3Nzc9aQUPUt199vKhpcl1+KgmFG+oQKbwVqrOYYv9teAFQlWW1lxbpkuGZtLC0xbpFcNHDioKW85po4c+bMukgdxfhjLibhaK+ctvd9Q2n8IWfg+GAyeNvBfeZ6ACV0OVaL2Xchyo+rjjqztlkdep1C/P2x4zWlT8xYf9uY8dEW4q4KwbYXAD7TOgdCOyp2wIWUQ1kUqYfZpU4wnwavLAC1OGlOKp/PObZ9ONaXyTSpg/Nsi6U+Gjl069+YlpKtXrb0lLDyMW4fQpidlPRlS1SsPUcAhQorX7oURU6HrW1dFlduHpUkAAC8GMD+wf+/AsBXc9LlfnJbANQ3KR0eOnFiywsppwoLgMvIffn4vrVpfa6JnqIthloALg3crCOn7r5+NWltNpvrInhc5al+p0I+fe0QUo8Ui4EqI6RN1LehmwZtZXAWtzkIEbpVbURNWQT+MQCfUbH/AL71et0J7AOXsYZqV7npis07F01VaytmORzNSv9NnbrIXTDP1Z4u4crpQ9/F4Hp7q/wajQZ7Q5Tv4DMOfBaOem9ekGIDV6DY3pshqCkWgPm9yptzXEUMquIXFFIEQAPAlwGoa79d9qWr4qlaAHC141BtLrVzcw0Ol4CLKaPqQRtavq7lhoRDmppebL10JmtqtJwdoLb+MbVDfUepTQhTbdNqtWS9Xh9eOOPT9F0bomzauc8ysSG2zW1aeqpg0b+pWtHh0pOCFAHw6cG/9cG/zwNwyZeuimcUO4Eps54Ljqln+17/jRsFEgpX2lD3FZV37pj7UHDNaLOuJhOJbQvz6lBbeVKGrQGYv3OUC8riaDab8siRI3JmZmbYRpRVoiyRubk564aoXq/ntFhs4ai2ulF+f8444tbfRIhVPQpFJ3X++ZAiAH4RwE8D+ByAOwD8NoCHfemqeKo8DE7/JvbgsRBGb/teaZFqkpp5VKGN5GDeNu13FAg1/U1BnSs6zNb/NsEQAi7jNS0EM/Jqfn7e6jIxQz7r9fqG+Hu9ji6mpQtV19ighHYOZkiNCZfLpWoG7ENV0T8KKQJgYrAOcAbARwf/F750VTypAoDLoLmamInQDrNpgY1GQ87Nza0TAFztNoauXBPOdQxBlXAxXl/8t0rLWcj1wSXodQYYov0rmAzVrO/q6qo8c+bM0PVihlaqzVRLS0sbwkIpF5NN2XAxUf0bfQy7xganLVIYIVdgjULLd8GmDOYSSqlRQF8B4GWcb6t8qrIAuA1dpZZg00o5W/FdeUkZFzaZWwvOhRBafTtAVTubkTGcOriYlL6TVb3Td7PGTHZqsVOv7/z8vHzxi18sH3vssaHg05m6Sq/y4GxQ87lFfAve5vlCrvDOkB20IaDmUg7B4iszh5K2WSyAuwH8OYC/Hvz9HQB+z5euiqeqNQBuQ8daBhyYg1yfRKG+UTOvmHC22ElXtSblsuLMIwN8rpdOp7/eY0bScOpg669Op7PuLBvqe0p4uFxwHBfS0tKSPHTokFxaWiLroJehu3oo5q2iqPQTQE202205Ozsr5+bmNlgItgVvys3j20Eb6j6j+s42Zlx9Ego9PXfujNLqSBEAFwB8rREFtGUXgWMQyixDzF7KVPctUJuba2L8+1xNcNSw0dJut+WhQ4fk1NTUOo3fx3h96z1coW9aAOfOnZMXL15kXeaumK6tX5Ulow5t49IZo7D0etdulHOFcrbb/Zur5ufnrZu4OLT5LIBGoyGnpqZko9Eg623CnIuu8a/TFRvz77JubO0bKijMdDGwCYAJ+LEqpfw7xnfXDaSUWF5eVsIsGbVaDYcPH0atVmN93+12cerUKXS73eFvQghMTk5CCLHu74mJieHvtVoNd911Fx5//PF1aU2srKzg7NmzWFlZQbfbxenTpyGEGObtgmobAMNyKXo5eahBlrOtXdixYwd+9md/Fnv27Bn+pvdNt9vFyZMncfnyZUgpIYTA3r17h+8pWpeXl3H8+PFhmyjo/SWlRLfbRa1WgxACExMTmJiYwJEjR7C4uLjhexM7d+7EbbfdNuxXk449e/Zgbm4OTz755AY6qP5SUP126dIldDqdYX4uWoQQeO1rX4v3vve9+K7v+i689KUvxU033bThu8nJSfzQD/0Qnn76aSwuLmJ5eRmnTp0CgA1zQUqJlZWVde1qjm99nJt9um/fPuzYsWMDDRRUf+/cuRPLy8tDuq5cubKh3fQ+o9Dr9XDp0iX0ej1reaqN9Xqb7av3pz6XTFpdcyR0DrJBSQX9AXASwA8BuIT+JrD3A/gVX7oqnlwWQNV+ah8o7T7VDWX7JlRzoNomlDb9TJeYM9tdmqtvd7DLjaPS//zP/7w8c+bM0Edv1l0/fsJ17pDSHM19B7obhbMZzbQAqD6waaiutQFVXyqckwPKlWXWnzo+2tS4U07jtI1l35hUbdhoNMh9IZSlYOZH3WEQY7nrfRRz53OsFa8DCS6grwTwMPq7gT8D4CEAL/Clq+LJIQA4jZnTH+iCznBiwjtzDAwqz9i6m75el6/Zl4/OQE3BojNom+vA5vrp9XryzJkzJGPTBYxibGp9wHad39GjR+Xc3FxSFBHHHWLLs9VqyZmZmeFuVUp4t9vXLl4JOU7cXMzW87MdMqfaRBc4iulRawomQ/S5ymxHJ1y9enXDwrvqy6WlJbIevj5S6w/1en3dmHAJDXNtQ41bVwQXZ60qVWGNEgAAbgAw7/pmlE8OAcCR/JwG5w5U1/f6YV6h4Z1qgFcVex/CxCgG5gsJ5fih9+/fPxSOc3Nz6y4N0Sca5T93+VoVY7t69aqVIan7F5S1oTM21Verq6vrQittPmaTsXOEFqdf1tbW1h2l7FIGOBFRNu3abMujR4/K6elpeejQIdlqtTbkY9tBrC5/0a0YU8t2KUJ63maoqi7UTSHVaDSiLRDKKtSFnW2e6wqDfpcH1c56hJat3ikKqZRpFsDvAfha33cxD4DvRT/C6C8BHPV9n8sC4Haiq8F9nU91MmWGHjt2LOjcGj0tFcHCqTM3f67WYX5LaYGh+euaYbvdljMzM+Tl4kp78oW5tttt+eyzz8pz587JZ599dngsgrnfQteozcmsGIpyCynXgop5bzQa69pYTex6vT78Vy2cUkKfGj+2dmu320MhubCw4Fx8XVtb816qoudr7kNRY00JxkajIZ966ik5NTU1DDv1CfRGoyHvu+8+OT09va6dlDBeWloalu8SZDYFTrcA1HhREUqUJWPSZxN4vnq55hVlDdjaPLclb8ImADiLwP8fgMtCiJNCiBPqSV17EELcAOCXAfwvAL4dwNuEEN+emi+j3HULNKELuApUOrVQs7KygsnJSQDA5cuXcfLkSQAbF8dqtRruvvtufOpTn8LKyoq3TCmvLSbVajW8/e1vx759+7Br1y5ywUn9Xy2EhSwgmfXT8/V9W6vVcP/99+P+++8n21XRdujQIWu7T0xMYN++fZiYmECtVsOP/uiPYnZ2drhgBvQXTxcXF8lFSr0sVe+f+ZmfwT333IMHH3wQjzzyCE6ePIkvfelL6xYor1y5gkuXLuHKlSsQQmDHjh34zd/8Tbzvfe/DysoKer0ennnmGfz93/89er0eDh06hJtuuglXrlzBo48+6m3jf/zHf8S5c+dw4MABSCnXLTCrsbmysmLtL9XWAPDJT34S73znO/Ha174WO3fuxC233IK1tbVhXdR4XFxcxIc+9CHUajVMTEyQfann++ijj+LKlSvD9yr44MMf/jA+/vGP44knnsDNN9+M2dlZ/Oqv/io+9rGPodFo4OGHH0an09lAc7fbxW/8xm9ACIE3velNeOKJJ4aL3c888ww+9alPQQiBu+66a7iISi1eqz6/88470ev10Ov1hvW44YYb8PrXvx433HADAODGG2/E2972NtRqNVy5cmU4x6gxbC6w6n/7Fsypd6p9d+7cicOHD2PPnj0kj1FtPjk5iV27dq2bxyMDJRX0B8A91ONLx8j3NgBPan+/E8A7XWlGFQYa63OzWRcuH7hL+4ihSy9TX4xNNSFT/JAcq8sFmw9VuQ/q9Tq5jb7dbg+171arJS9evCgPHjwoL1y4IFut1lArPnbs2NCSMn3JCwsL8sUvfrF8wxveIKenp+Xhw4flLbfcIt/61rcONcxWqyWnp6dlvV7f4PYw/eW6Jq6sDeVCUhogx2Izv2k2m3JqakoePHhwQ2imuX6h2tM85ZJqM/W+1WrJw4cPyyNHjgz91UtLS3J+fl4eOXJE3nPPPXLv3r3yqaeeIl1h1GYs3bJot9vrrIyFhQVy8brT6ci5uTl5+PDhIZ2mZWi6H00r3OUulJJ/g5itT2LX81LnqAvYbBfCAHgzgP+o/f2vAXyA+O5eAOcBnN+9e3cljWOCawJy0oZObBdz5DIG5bvWIzNSB1fKAKXM9hCT14yiUAxI+VZbrdaGia37qmdnZ8koGUWLvpai2k655VZXV+Vjjz0m3/GOd8h6vS6XlpaG/+r0mK4kqu6qLHNxULmwOBefm1Dv1aXrBw8elLfffvsGf7fehnq5FG2tVmtduypalftKucCU26vVasmnnnpKvvKVr5T33XffOgbPEWKrq6uyXq/Lc+fODYWAzb3Z6/Xk/Py8nJqaGq4J6fNLT2u6BfX35i5lar+NrR62uukunRABbo6LKhAtAABcRj8EVH/+GMD7AHyDL70j37cQAuD9rjTj2AgWqq3qWo05mHx5hTBGl784dp0jxd/JoVVPa2sLSvgqBqGOMp6bmxsydSnX3yBmhiMqTd/UPvVydKGi8lK+dfO9bdLbrDy9HHNs2DR0Pa3vaGWTWdl8/brWr2vjqo3UcRg2y9FULNRdvsp6WlpakhcvXpT1et163IUNzWZT3n777fLQoUNyfn5+Q7SQOS9sGrreF3r4LsVs9cV91457qh6UQqPWlJaWltaNOcoaMq0Ec1yY4yYHUgTALwJ4BMDewfMwgF8AMAfgcV96R75jcwGFNG6MBWAbTL68uIxaSklqcbY0XAbsmrTqHfcCEqpddIFli4gxBeaxY8fkY489Jg8ePDjcEaq0O+Wi0E17innpbWUyOUowm1ctKmvCpqHbhDbVD+a3rjGh6KJ25LoUAFf7m+4JVbe5ubkNYbe2+igGrCwlPX3MuF9dXZUf+chH5MGDB4cRTTqDN9vfVndbmRTDNr+zjUeOotTpdOTs7KycmpoatocZtGD2qSlcuQpSLFIEwJ/YfkPCxTDo3yvweQDfhP6FM00AN7vSVLURLHQy+ZDqQ/QNhF6vt4EJhg6eGAtAbSoyj6rmpFc0qbBFPcLENTkXFhbk1NSUfMtb3iI/8pGPyNXV1XXx2bbJpDM6XfvlbFyi6HCFteruFarOlJBV5fvCHm3jMkRZcPWRqpvtqkgKSmNVVpG+lkJFtJlWCtVXyl3XaDTkxYsX5a233irr9fo6GnVtXgkfzjgPccPGMF7d8lJRSDarM9TC59DMQYoAaAK4Rfv7OwE0B/+v+9J78v4+AH8B4K8AvMv3fVUWAGUGpuSdko9Ne7ctntmYtY+pucq0QcWdHzlyxOvvNunQLxnRBYheP92MVvHT9XpdPvXUU/K+++4balSdTn8h8NChQ/LixYtW7dt0AVEan40BU/1A7Tno9XqyXq8PQ0dtfWbSpSy3mAV6G1Nz7Vrm5MllOJSQ1ZUSs3zTctB3iVNhkPV6Xd56663DvqUYqW3/hwlzTPoYfA5mays/VqvPYQ2kCIDXDNYB/nrwXBoIga8C8IO+9Dmf3GsA5kTVN2zEwsecbTSobyltUh/w+kRy0UnlYxtI3AGma17UxRWmGa1oUDtldSZP3Zo1Ozsrb731Vnn48GF59OjR4aFk6mIS/SrDc+fOyVtuuWXdwqHZ3qbrgGpP20KjeeSETetUwkgXar721MeIz+q0KQAUU6vX62Q0Uihs9PvcJTZLSKXVmbZayzH3Tkgp5dWrV+Vjjz0mL1y4MIwE8u3C5So3uRm8D6aCY67PcKyCsVoAww/7J4J+Hff7Kp6cAoDS1LmM0OcqCeks0yVATQzdQtEZgMsfbxMkKRYAxRh1+sx21BcLTX86dXwCZQE0m025tLQk5+bm5MWLF9eFUM7OzsqlpSWrxWWbXIoBnzt3Tk5PT2/YzWoKOps1odOtT2x90lPl6wLc7FsTel+7zvxRY8fmew6BbTz4LGUqnW4t6IJaX5cxy1EhrTMzM+uO4nDNOY7VzRnnIW4aPU8OfdQ6hmoHV5+N2wVUQ/9AuN8f/P3tAN7uS1fFk1MAUJo6t6FDtXwXTHNaN41d2o5+yBflRqhK03FNcrMd9RBJNdHVBKNizW1ot/s7Ow8fPjw87jnEZWGrQ71eJ48aNhmKzoQprdUU4rrgpia3/r3ufzejg8z2UjH4Lk0xhHmZ7cFRAFyCjYK+/mJzqZmMz3c2jgnbXDZp5Ch4Zp/52sYnFHX69Cs2zTHAjRiKRYoA+H0AP6j5/Z+Xsvib8uS8EcylDbsW/NQ3IX5+F4M06eEsSqq/leapx2tT34bQlutbKWn3iAIVxWTLf21tTdbrdXnhwoVhmGEOuBgNNVZ0DVufuLb1BBtDNie2imwyQxNN99Njjz1m3W+QglCFxsW8zT70uYViBBaVjzluKOFrzm3O/PMJDddcNue7uWbiUqT0No5ZJzKRIgA+M/i3rv3W8KWr4qnqTmDzG1cYF8fk85Vrag3mBiaKRnMiuQaITeuy0Wb6ozkhmhy4JrhqR33dxZa/PpltkyGEmegTLaROehmUMHApDZSwMQ9zczGQZrMpjxw5Yo3AioE+lrkuFMpi0WG2p2u9yiVIUuqjW06U+0/f9OdboI2xpGPz0tPpPCKHFZ8iAD4O4BsAXBz8fSuAT/jSVfFUdSew+Y1tMnM1f5sGQFkA5gSwTRhTANgGSK937UJu3ffuolUXRCrunHv8bcwE0aEPep9Vph6q3mpiq7UF14I+VSZHoOv0tdvtoTuI48fV/doqD8XU9U1trrJDIrs4sLW9Lc92+9o1kJyNb2YZtvJDAyZssfouVwml4bu09hTE5qX3se/k2FCkCIBXAvgTAH83+PcvAOzzpaviGfedwPqk92nVse4hrvZgo1lp9EqzpGjlCCROhApFrynAfG1LMR7Ogpop2Hq93rrdqWfOnLEeAUyVw+kzXQjrdHN81sqNZe6Ujd1Yp8Ch2yVYQyw918JtaNm+d9R703I2LWqXqyTUeh0XbN6CVLqTooAGfv+bAewB8HxOmiqeqgQAt5Fd/kwdvglnahwc7YtbltIi1NHD1GR11Te1fHVAm2K8FNOgtDGdCXPiyV3CTx0L4Ttmm2N96KDccDbN0ydkUq0prrD0acZ6O/hcWTHCPEawmXQrxh6yW9dFmy1NTisgFma/js0CAPD9rseWrspn3BZA6gDRGYDt/1woTdtcAFblKEboOnc+RCtTgsXl59bpUoxXRfCokzP1drCFxtk2NPnaicMUOXXlfm8yTnNB2Nw34CvPrJ+rPpwxw9GM9XI4riwXUse0GmPmJjmbAONaxT5auYIyFbH05UCMADg9eJ4A8LcAfnPw/BcAv2VLV+UTKwCqluxcDUTXdlO1frVJSr+v1OY/9O0Q1cvXmbNOm2IOnDhzva5UVJPNAvC1p0sI5ejjFIFguifa7fY6dxqnn7kuDy6tShibxxO4yk+JyjH73WcFm+h01q/j+DZomkySa6GbdFCCsgqeEcLUY4+TsSFlDeAsgF3a37uuNwGQW5rG5K+iRlJ3aaryzB25avJQzNnmu6bo7/XWRxfp/6eYCMU0XLSkwpZ3jj5OyYPS1m0MnZs/16KxMSvT0pqent4QRZSL0enauyvwwGWVmGPJJQAp2kMEgC2fKq0Bl4BNHSs+pAiARePvCfO3UT1VWQAxk0CX0Jz0jUZD3nrrrXJ6eprlDuC4abgHTVHajkubNN1LrigNiiGnapIu2PIepQUQW1ao1cctRykE5oKyaWmpzXc6g4wVSpSw0cMrOe4m39jiCkAfbSHwuc3MMkLLtLW3+XuKh4BCigD4AIAnARxE/zaw34fn3P6qnlEsAnMbu9FokLtIbTAjRXwTzzcQTbq54NSPU7bPjTMuVGG6m/C1O4cGM6LMZ5nZ0Ov1o5/UsRZqoxjFPKgyUhlYyDgwaYnZA5DqGvEpJ772MOufS4D6XGNj2wncT4s3oX8BzPsAvImTpoqn6kVg7oYYKddfHp5SJldTosxfjl+3CtqkHG1YXYjWlZsuiuH42ofjilCL4/oJpzF+fqV5q+sq1S1Z1L0HsXAJEZfmStWZiqIKGbtmlBmHZh2UxRqC0DDnHBiLBQBA2N6FfJPzqfpGMGVKczampHYK18SzvdPDK0cd31ylG4ZCiNaVk4Ze79piu4/h6OAIALMNY+jWlQAV9qvyM+894ORr065d7W2mcX1ra5eQuitLWh0YyNXU9bJSLNZRKj45ESMAPg5gFsBu4/cvA/B6AL8G4KAtfRVPVcdBx/rzUsrhbEO3pZfSfdxD1VoJ149pQyh9VfQThxEohcC3n8BHb640Lu07po3Mb5R2rSLLYiwuFx0urdxcD9D7h9qoRq07pdafg6rnVqqAsiFGALwAwH2D3b9tAH82uA/giwD+A4DvsKWt6kkRAFTHjUqaUxpsSHhgqHsgV71C/ZW6Ky20PqMGxxUwKoVA0eNTAkyfOdeK5OYXsvjvqosOjoBot9vy4sWLw6O5FW3UmT261WNGnoXMKYq2cUGnNdVFZUPqGsDzB+GfX8f5vqonRQBwmGhu94E5OWN397kYq21C5dAgQieIosW3cSfX7kYfHT7tb1yL11S7+hY3bVqyzjhC+kp9z70QJ3ZuuKwW/e/5+flhUIVi5tRdAJ1Oh9zlrdOsTleNPSJDIXcsvg1639oujklFkgDYLE9uC8BEDo3At0M3NNqIok/365quINeiX07Xi0uLtfmdx2VxxSJFw3blF7Oz12c1xPatPl7M+uqMKLVNbXVXv+ubvXzjzrYuo4SAOjXVd3eCr724C86x0NtE3zRYxR6EbS8AOMgxyW1+VDP/mAlFTVpzMTh32KgtPbWg52u/UJdSaB62ndCx/anXN4dQcS0Ou9ogR7vZ4BqTOgOMoU9HSvuZaX1HTNsu7gmhxxSAoQhtE92CCVm456IIgBGAM2hS3UF6HpR7yef2SBlYevrYXZcUQvKyab2dzkbfaS6tNZcFENtmKcqCj15XHbkuEA59se6UGFedEgKukFNf+6SOb6pNfDRQSkwuIWATABMoyIZut4sPfehDqNVqmJjoN62UEsvLy31pO/jm1KlTWFlZweTkJFZWVnDq1Cl0u92ksoUQqNVq6Ha7w7KobyYnJyGEiC5DpZ+cnMQDDzyAycnJFLKDUavVcPjw4WFdVdvVajXcf//9uP/++1Gr1TZ8GwO9vqltByC6zWLqobcN9zuzjhMTE9i3b99wLPd6PVy6dAm9Xi+YvpWVFZw9exYrKysb3plzxKTv9OnTQ5o6nQ46nY51jEsp0e12sWPHDms9Ad58SYFqk507dw7rtry8jOPHj2N5eXlIg97eapydPn0aly9fXvdtZaCkwqBBvhHAbwD4YwA/De0YaAC/Y0tX5XM9WACmNm5qAhxNxKed2FwTuX3sMbTlKickHWdxfJT05GqbUdQh5LuFhQW5f/9+5x0LMT531xzxWXm2vEw3ClW+aw0lh7Wsjt9Q+zM4+yDU/30HOIYCEWGg5wD8BIDvAPB+AJ8E8A2Dd3VbuiqfzSIAQhh0zEDqdOgNaRTDc5mNqYOYoiOnkPG5Bbj0UwzEDB0M8aGr/KiIG5dpnmOnq6tevjxzMzGTFnXRzurq6oa8beOCIxiowwRjd7+HtI8+/nIrT0pYzczMDA/gs81f2wVHo3ABuQRAw/j7RwA8A+B/xOB6yFE/VUcBcb/nLCDFaF/6AKGOpLCV6/s9dlFpbW1teKeAzlxzhU9SURYxk5JiRrrgsjFS18Xm5jtdKHB3uoYyFQ4Dd+Vpvsu5TtMz/Oo+y5ZDr+19iIIVAlMx4ygJZhtw55FNWJlji1I2Usq1IUYAPAPgBcZv+wH8JYCOLV2VT+59ALHf55LOZhk+y8GnyXLvSg3VrKn7d2PyM2mmFsypyJPYfRMcRkpNPMqq4txG5RLoMZYMBV8EjF5OTgFg5her5FDvTYWCkyZWodGVi9BjX6jwzFBazHHnGj9mmhTrJEYA/B8Avpv4/RUAztnSVflsFgugqjKqcBlQ+cZq1jH5hQ5q6rAtl6bOoc9VN19MfowGS33DuaKSQzclfKtyAcbQFwPTUrXtTchVDpWvS7Ca9JljxRQkLppjPAtjsQA247NZ1gBsiDX5c4V85dLKcpabOqj1fQ5cIRLSDz4XRA5tdG1tbXhJPbU/JAQmbc1mU87NzQ03PnHj5VPLzglTw9YZbQ7t1yyHot+16ctmrap3tkvqQ2nOzQ90xFgARwb/vh/ACfOxpeM8AH4JwOcAXALw22AeMZEqAFJikVMYq08L5h5BHUrPZgBl4ocgZvNUSPtUaRnqjE1ZAK1Wi+168EFpn+fOnZMzMzNyaWlpHW2KqS0sLFSuQecAleeoxrqPN7jaclSWfApiBMCBwb/3UI8tHecB8AYAzxv8/z0A3sNJlyoAYrZ26xI+djHVpwWHSvwqBggXMYM7hd4Qd04Kk8i56cn81ryghXvFoQ/qexUyqO4C0N0aCwsL8tixY0EWjo2hucKbubSmCOsc6WLzt7XlKMrOgU3rAkL/splf53w7DgtAaVncA6YoxHSoa8LEnsWfY2BxJ36IGyUWpusgpF9MmprNpty/f79XW85hXfgswpjxpVxBZow8VZZeDqXVm+1JWWGhc8lmycW6CF1txRl7PsHncw2mwEfTpnAByWsM+qUAPgjgKQB/qB5fOu4D4HEAP+J4fy+A8wDO7969O7oBUjUN/ZalUWgdoROGwzxyaG3cdhyFleK6E4GbVrUvFfI6aqS0d4irzcYgdWVHz4cai6H9G3oZjC9/V7vYBJyejoqM04MNqnRx2RQWRWuoS5iDFAHQBPCTAL4TwKvUw0g3D2CReN6offOuwRoA62axUYaBmshlwnIHbuiEqcIC2GyuG5Nxqbj0VAFQheZlc6OEWkSjdPfpTNAX1aIsgNXVVWs9dYRu+kthwK62Nt1yqs8VU7YpeDmtZ5srebNaABd838Q8g7WEpwF8JTfNKMNAuQidoKmmK4WqtJWYxdtQ5saxqijNSU+b6gKqgsmaebq0Uoom129mGo5LkMtclBDgWJPqtjRqw5zpOvEJaxdzzDm+KYFAWT22+qZ4AEZ1vwCFmEXgrx88P4f+zWC7tN++3paO8wD4XvRvGNsRkm4zhoGGMmrfVvYYjIKB5UynMyMus6EWU2OEpYueKsx9igHbBGtMm3c69Pk4NgFknpLJoZ16t7q6OlwU9VkLnU5Hzs3NydnZWWuoqrIkKGHearXkzMzM8LawnH1lCjyXhc1pOxdS3FqpiBEAfw3g84N/zefztnScZ7Cb+FkAjcHzK5x0VQiAUUrlGB9fTlfKKKwKTjp9IuSoX5UTJzdSNF2bhUDFqNsE0NLS0pCZcsug6DddRS6BrNNouozMyDyq/Eajse62sNRd6GY6rjUYEpIcY82FWIah2LRRQCFPFQIg160/HJdJqN9Uyrjt9zbktBRSaMktiDqd9TttY91XowDlzuLC5oYI6Vff8RBcLZW6qD0kskiBo4Ctra3JRqMhW60WWa5yRal5xfGhu5itbey4xqLL5ceFLf8c8zZlDeAFAH4KwG8B+E0A/zuMM4JG9YzTAvAxFZspbn4TGsGjT9jUgZBTW8pxLgoHXJeSfk0gpy9yupBiECMEbS4z2xgO1UJjmZ9Jo55HaKQWhwmaZc3Pzw9DeLkWtqv/Y6wMjgUQU++YvCikCIDHAJwE8LrB80EAZ3zpqnhyCYCYBvUxlRALgDNRqXepAyE2PaXd2c5FidkpGjL4bd/qbevri15vvd83p2XkQ+zYU21O1cvGsEJDl13t4HpnulD0Ml1Hefjqaut3Pc+HH35YtlqtdZZ16O14PobvGvPU9zGIETohSAoD5fw2iieXAIg1z1LdCimMJkQQ5DYluWW7JpULVX3rykMxKdeibCqotkpxDVBn5VDl6Awr13oT1/1hClf1t3nXra7FU3VxuchS9n/YaLYJDn3e207VtfULFxRvya2UpAiADwG4Vfv7FgD/zpeuiieXADC1RX1QVhn6mCLVKS3cpo2lHDWQS/PIZQFUQZ8+4VM33XBOkdTzTqGfwxzNMkLbNmb8Uye46mUqH369Xt8whl1j1XfKaUq8vE84U23tsww46zocpSCHsmkiRQB8FkAPwBcGTw/9uwIuA7jkS5/zqcICMP/v8x378qsK+qD3rUWkxCvnslJiv41lXqHlS7nREqC0P19+riCCqiLMUoWgy0IMWTvR09kEnQoXnZubk3NzcxuihlxM3CckQsdqiOVsMndqXJi0q7qurq5aaeAoBVXwkxQB8GLX40uf80kVAFTH5bYAYjWuXJeGUOXHMuYYhprqyuFaOhTdLm3e1y5UWZy2brVasl6vD2PUOXVM7SMbuC4cF116SCk1FnR3h2vsqjLU4vz8/Hz0vLLNjVBLwKUcmXPfbEtb++m/K2XAdeT3KK1wHSUMVOaTrK4OCmWCSuPiuiJiBkcII6LScX2bHOHBndShAsg3uX3tG8OUfRqqS9MOFTYcmlzHh3DPubFFnan2VTtm1XeuoAhlAYQGPoTCp1Xr4dcua8JVH70M3bLT24tydW0WFAEg84VacczDGAugil3COk2u2G1XupCYdVWnRqNBHqwWI4RTtSaXcEgp07avg5O3aWly3ROcxVFdAPR6/R2sc3Nzw0gcrgAwhbG5uJvip7b1SYxgoGjRxxl1zajPAqDKoCwAs83b7fa6tt4sKALAg1B3Q6dz7Rz22K3hZp62I31dNIQwHVP7cZm5MWWpMnSrJjYvHTF9w2UqNmbMOR4gxaLkCl+TobkigCg6O51rxzCoceoqu9dz334Vq7VTY406eM7mXvHBrBNlAeRai9HzM/P2bbQbF4oA8CBmcId2NkfzMi/1sCF0onC04NTBa7oJUi0tThobo+QyZcr0N9sqxIXGtVZsTFanixMX7hsHMdaG0phD/Pa+PG1uGvO+Z5umHVt+itCywbU5kzMmXH/7BHosigCoAKEMwDag19bW5Pz8vJybm4sSJhzGzRlMZj4xVobPjE49UMsEl1GG0My1jDj0UOVxLL0QqyuUYfryVncjcK1RKf1WRcwO4xwM0CZ4UvK1uclseZs0uP62/T8VRQAw4GNgHLgYkm3gqduoYiIlpMxniZi/x1gELo1UuSMOHTpkPYwsFFT7VnGeegw9FFT7+I4frqJs9d7lutTbkJoLMZp2TkZGIVSIpNLjagPKdWjbx7AZLIAJFAzR7XZx4sQJnDhxAt1uNyqPWq2Gw4cPo1arDfM8deoUut0uhBCYnJyEEGL4vZQSN910E9773vfi277t23D69OngsicnJ/HAAw9gcnISUkosLy/3pbsFy8vLOH78OJaXl9f9TtFHodfr4dKlS+j1ehve6fXX667eve1tb8POnTu9ZXCh06zKW1xcXFduCDjtx6EHwDAfPU/VPvv27cOuXbsghICUEp1OB51Ox1ouhy5X/0kpcfnyZZw8eRJXrlyx5qHacGVlBbt27RrSaL4325aqt4I5J3JC1cvW32abqP44dOjQOnpC+t3Wzt1uF48//jjuvPPOYZ8eP34c3W533fdCiOH8kFKuy8/2/8pASYXN+lwPFgCVJ9c3WpXJa4K7ld5lsXBOUI3xh/rA0fZCLADKSsnRH3o+entTsebqvb5RypWfrx1c9Ng2ventwY2GoUCNLa6GHtPenQ7vQhedPipKJ8UqMMeditRyuTtzrG+EAMUFNHpwzOWcZp6rTOqbkPBOHTmjKkInXuzEsQke2/kuIT51X9mK6Zw7d05OT08P/eoqf8VwVegsJ3oqlLaQcWbLm3OYobq8hTr2weYWtIVUcugOHcvtdlvOzMxsWNvglkMpLjbh7BO0qf0RgiIAxgBbx3E7NFRQhDLmcV5Rp5DTAjChM1hTG9Pf2fLL0T5KAKhY/OnpaWs8fsjF9DGKAzeNzQqgLD/l91Y0+ywAilmqsqi1Ca7lw/091so36fBZivodBqlju1gA16kASO1Qn3ZkIvRymxjNIrfFUiV07XBubm54ExbXTZRD8zK1fBfj0d0ZVWwIDKkP9S0lEHXXEsfdY7rAlHDkbqwM7ZNcEXIuC4BCs9mUU1NTcnZ2lqSVUu50Kyb3PCsCYEzg+j9taV3akYlQjTXW55pjF/SooNpQ97+Pyvcaqm2a1kLooYSc/LkWZYi1EGOR6UqNufaRMmdMKOHv2l9DCZXUvlcWwNLSElmuXiY1p3IoHzqKAMiIXP47TidzXBW5aA3Jj2KkqYO2CgGitGqXhpqD0VPpKTeBz/JQ9M7Pz69zH1DpQ+j2adKhGm5oW+jv1NixWUUp44iqh7473aeN56BBhy0fUxiHbtYLRREADoQ2dsjgSNVmYgaCnia3JuGiK3XQVkGrYgAuDS+1XJsGqU9q3VXC2ShmrgNQ6VNcOpRA4Ox+jm0L2ymjOmLOVuK4iLiWmC5kc0UDcueE6UrLjSIADFAMgKtlb0ZXhw5K84wVILmR07znluPT8FK1XddxyCbDdTEXV35U+thwV5vGyTn/KKQcvS1mZmbkrbfeKhuNhjVt6BqWyjvFkqH4gE9hyA3VH7YDFHOgCAADFAOIDYvU4RsoMfcCp9AQk2+s9mcrK7dFwqmTr5wc7e0qQ01qajHXpnBw24bjVvClp1wiqePGhl6vJ+v1upyenpYXL160atYpFoBNGPraVL3XF945At0sP7addMEb6+bloAgAAzEujByMx9RycgqfGHoocE1mV1l6W6nFPuUCqJLxKtiu/cwJjpXhYrS2PQi+NrLVh9vXnP51WcWxViX3tNsYZhuaRm/rlLOZUhWaTse9kS12LpooAiADOJ0dagHYJlpqx3OZiS0dJZBizGmV1+zsrJybm8si4EIFsc9NwM3TBd19wO3D3Ewlh6AzLRffYj+nTI77S//OdomMz+IKqbtZh9i5lsMCcCkDamNdakRYEQAZkEuT5JjaNu2RixBtkHJDUJpfCGMy65j7iI3Q8n2ullRNTm3iOnbsWDLjrspi4UAfd7b9CHp/unZS63mGuLZs91bkPOTPJpRG1fbmnKDGqIoGO3LkSPIFM0UAOJBrRyyVT+yEyGUB+NK6tGMfY6pCII4iD64FENIHlAXgKtNmLZnfK6bnO1EyF1SdfXsRbIqCqy4+Wns9+o5lvW1jLVrb96Z1mqoIcOFT8JTLLPQmOxuKAHAgJvpAQR9gVD4pEyIXXOW53vkmg+89xUBTNEQXYvLw9UOIFcbJ6+GHH5atVks2m03ZarWcbhRTG9b/tQmNnLAdZeDTxlPHtu+4Ce7d2Qq+dtLXp3Slp+ojxfW5UVVb6tiUAgDAgwAkgJs4329GC0AfYFwLYNTgWhyhGnAMA61KIMZM3BgBlkpfo9GQd9xxh2w0Gs5Jb2r+amxxLnbX84ml29Y2qYLHZ0XqkUB627sihGx15Yxf6r0SOLkVCtv7UVgdm04AAPhGAE8C+OK4BQCFELM1J4NPzY+aYBwmFjLhuRaFMun1qw+rFog+V4yN1liE5uFTElzMwGxbjsWRi1H7fufmaTJW17jTj4nQdw9zLUjb76oOtsi7Tsd/fSoFX99RG/x0eqpUEjejAPgogJcD+MK4BUBOt0SK1ugaJFyYdNsYOHfrua1t1AQxNW5X+VUMdErghWjKqcjtvqKYfEhcukKIsIhFaL4UY/Vp7/pBcfrfHIXE1n6q7VToJbW+YrMMYq0us0xXm6nycx4IuKkEAIC7ARwf/N8pAADcC+A8gPO7d+9OaoQQEyzFnPNtsrGBo72G0mVj4CkRRromR0Vs2KIbuNptCMw8zfK5DDq2/KqsCF2Qpe5MrcrFoPzn6oBCH01cV6tNePk0c+78NpUE9d53GFvK/KSsDh8PyRH+qTByAQBgHsAi8bwRwKcBfK1kCAD9SbUAfCZhDnPPpj1wJmGIkIhlBi4aQ8DRTs3DvjgaUozF5bJAuG0TaqmMyrLgRr646OHSGlonUwC42lBXGnx9HDJPOQqGLx2XOXOsyxDLrNVqydnZ2Q33Y29pCwDAXgArA8b/BQCrAP4GwKQvbaoAyDFJqvg+JK+YQTdq6BMqNIomNfoiRZO3xbZTcAliTmy5TwiHtofJjEPrbjLAkHSKNl3DN9vHp73reYYoJ6GC21cfTr+5BA7Vb7Z86/W6nJqacp6PlAubRgBsIGATrAFIOV4mmmqZVOWGSM0nxtKooh+4zJRyBfjypASxefywzZ3AOaLYx5TVd7b7Ijj9qoSH745oX1vYNOiQcRCiQXPr56NbF2D6jWwcWmz10tOafa3eqZDgKo8rUSgCwIMcHRAbTmorO1UwhJSbi/Gm5lPFRNAZs8/E55ZNhSWaFoDrcDEfU+QKI993nP7wWQ++dqEYWigNnLJ0NxJHSPtg0tVoNOTtt98u5+fnZa/XY81nzhxV7h7XcdhVK6CbVgCEPOM+CsKHlA1lFEIFQygo89l34YiPGVStyfjgolmvG7cNbfVRfb2wsGDNJ9Utw9GafRoxR2PmMnifJlyv16OtEM63nU5HHjt2TC4sLKzbSJerLHXujvLJc8YIt7xQyyY3igCwIGfDuzSGnOVUYQEomKarLRKiSldZSv3a7f4VgL6zU1InLufo4hQBYKtLKrOm3ucQ6p1OZ3jvsl7f0L500U9ZAL4w5BAaQt1LPmvFdC9xNolWJQiKALCAM2BiXDtmHtyBaZqO1GDgfBMLX94cmlKRwjjb7XZy+Fyu9k3NR9VFvzM3VZOsyv1gujnUb6H7Wrj0K8bvOoxO/d+1iYxTpg2ucWq2KVeZqkrBKgLAAs6EiL2pyLyL1jTpKRdMu92WDz30kNPM9WnpMXWO+cZGUypG4TpxIUWLDMnTl7c6j6derztdddz8bEitm8loTaZXxTWHKm/XWoD5jW2ROnbsusap2U/UoX7FAtgEAoBCLgvAt+1dMSulOanBurCwIPfv3y/n5+dJRpaqWVJ0pGqFuTTl1LxyoAot2TWGbMzM/N08F8hGTy6Nnus2UbSYi8BV9iNHEJptrtOjM+9YOjluHZclEFJualsWAUCgakbDYWydzvozT3q9fvTB/Py8nJubi/Yfu+iJ8ZtW2VaxkyIFodoXh+G4wLEAXO4MPQ/zmGBOX8VadJwoJEULN2Q1Z9+66HS1eer6jC1/n3JFWR6ciKZUoV4EAIGq3CchaWwuixwD1IStvpw65XTxmGVzjm7ILZBsk9e2WcnGpEPqqQSIzUXlq5MtD07fxFp0IUzdJSSVpZD70nMXnb6gjFy74UMErdletvONOGWFoAgAAinukxC/5qgETZV5hqblfm9zi/mO0khtU9vktd1G5VpIDGkbZfH5Fql9zCXUanLV3ZY+lklSmm3MvOGCM2ao+scK8xC6fFaJS+nIiSIAMkHXZPQBznEfVDHwR+EuiQHXvOVMEuo7HzN30eRjvJR27aqHbyctpfX5mKrPvZDS92beZvy7iwYOKAZb5VjlWI36N+r/5smcuWnkjGOblZLDQtFRBEBGUAPcNVlcprH5DaezQzWYcQgKNbBtkUwuhJjRUvIZFVdIhy4Qqu9tpnwOC9BkCDmtuUajQZ5JkzpuRjXuOOXojFZX4mwCNpWGFEGv6Es5sddEEQAGcgxOrhaqm78xkRs2hsfRrrk+xtzg0piSt89tFJun0oiXlpZYtOsCnqIjh3WYwwVmQ647sRW4LqZRgjNmuHRS35lrdqbLy5Y31fbKLcS5O4CLIgAMpJq3rk6xmfy2C759+VKDl8s0zAij0LrEokq3F2eNILR8PV91sFqM4AwdV7HWy2bRrjkuPB/zjamLLU3I7xzXEQUqnSkATCvDVpbLTZTzXuIiAAzETiDOhKX8jSmmZswkVL+bwoczAVJQxeDVwZ20MXUzrRaXwLYhtP45XTtVwNeOMZq1UkwUw4y5fzeEodpga2ufNU7Vj5pnHA+B65uc87MIgEyIsQBCtP0YhAqlqujQtZ6YSc0tI9ZM930XwxBseVKBAjb4FIZRjCFXXUKUD05bmdYW5+wmDl026zCkPr5vfApdCtNOtWJdKAJgCyNkUMdq5r3exrNeTKgBzL3FajOAw3BDJqIpBH2bu8zfXK4KKqwyt4ANBdVWnPUYlc52j0FI+ao8dYyK6e7UaQzdX8Ox0DhKBLesquaMTQBMoGAskFJieXm5L4UTIYTA5OQkAFjzVN+srKzg1KlT6Ha7QfR1u1088sgjuP/++7G4uEimqdVqOHz4MPbu3YuJiQlMTk5CCMHKf1xQNNdqteFv3W53XRuptrPVRUe328Xjjz+OAwcOYM+ePQCAgwcPDiecmbeZP1WWohHAMC1FtwkpJTqdDjqdjrOd1Xftdtv7LVVfs62EEDh9+rRzjKl0Qgg88MADw/Gr6OGOjW63ixMnTuDEiRMAgAMHDuDs2bPryl5eXsbx48exvLzMrpee9pFHHsHx48fR7XbJ/vH1HxcpaaNBSYXN+mwlCyCnBqe0lEaj4d1padNWfKY8xwIIwWbRYCnk0uJ0iyiHZWT2l3KhtNttubq6SkaTcEIJ1Xezs7PBYYcumkL83SY93AgydWBeq9Wy5qvfvWt+74Me4rvZrVkXUFxAmws5zT19osfuUE71P1b9/bgRUz99YZCKO09pA5NpLywsbDixNsQdozNu34mVPprMOwH0966duabLjeMe4igSutuHKxT1trmexqkNRQBsYXAnOpWOO9F978el0XMnaKjmF7O4qUPX/ns9+nrB2DbTGbZ6lpaWWNYZZy0oZlFT0VSv162HGFLWkasMjr8+dFE3dq5c7ygCoCAJPkYwLk2Jy6CazaacmpqSMzMzLIabwyLyxYHHtpmZV0jfmG4pqn4xioHrrCQfTbYytor2vRlQBEBBEjbrZMxlAdjcELa/OTCPHwixsFzfU9+6tFrK9cI5EpyDEH99wfhgEwAlCqiAhbFEKDDApWtiYgIvf/nL8cIXvpD81hf5Q0XvAO5Im5WVFZw9exYrKyteOs38beVRtPkib/SIIZXWjNLiRBVRUOkmJyc35fgocEOYg3Yz49WvfrU8f/78uMko2IKQUg7DKykmZnuvQgwBkOGMrjxd+YekDS2rYPtBCHFBSvnqDb8XAVBQEA8p5TC+vGjABZsVNgHwvHEQU1CwVSCEwK5du8ZNRkFBFMoaQEFBQcE2RREABQUFBdsURQAUFBQUbFOMTQAIIWaFEH8uhHhGCPGL46KjoKCgYLtiLIvAQojXAXgjgH1Syn8SQuwcBx0FBQUF2xnjsgB+EsC7pZT/BABSypUx0VFQUFCwbTEuAfBSALcLIT4thPiEEOI1tg+FEPcKIc4LIc5fuXJlhCQWFBQUbG1U5gISQswDmCRevWtQ7r8AcCuA1wB4TAjxzZLYlSal/CCADw7yvCKE+GIkSTcBeC4y7fWKUuftgVLn7YGUOr+Y+nEsO4GFEH+Avgvo44O//wrArVLKylR8IcR5aifcVkap8/ZAqfP2QBV1HpcL6HcAvB4AhBAvBfBl2H7SvKCgoGCsGNdREKcAnBJCLAL4ZwD3UO6fgoKCgoLqMBYBIKX8ZwA/MuJiPzji8jYDSp23B0qdtwey1/m6Og20oKCgoCAfylEQBQUFBdsURQAUFBQUbFNseQEghPiIEKIxeL4ghGho794phPjLwZlE/2qMZGaH7aylrVhnIcTPCSFaWj9/n/Zuy9VXhxDiQSGEFELcpP22JesshDgmhLg06OOnhBAv1N5t1Tr/khDic4N6/7YQ4uu0d+l1pi4K3qoPgH8D4GcH//92AE0AXw7gmwD8FYAbxk1jpnq+DsA8gC8f/L1zK9cZwM8BeJD4fUvWV6vfNwJ4EsAXAdy01esM4Gu0/98P4Fe2QZ3fAOB5g/+/B8B7ctZ5y1sACqJ/V98PAnh08NMbAfyGlPKfpJR/DeAvAXznuOjLDNtZS1u5zhS2en3fB+AIAD2SY8vWWUr537Q/vwrX6r2V6/yUlHJ18OenALxo8P8sdd42AgDA7QC6Usr/PPj7fwDwrPZ+afDbVoDtrKWtXOeZgZl8SgjxLwa/bdn6CiHuBtCSUjaNV1u2zgAghHhYCPEsgB8G8LODn7d0nTUcBvD7g/9nqfOWuBPYde6QlPJ3B/9/G65p/wBA3d593cTExpy1hOu4zp76/nsAx9CvyzH0XX2HcR3XF/DW+afRdw9sSEb8tiXqLKX8XSnluwC8SwjxTgAzAP5vbPE6D755F4BVAL+ukhHfB9d5SwgAKeV+13shxPMAfD+AV2k/L6HvQ1V4EYB2fuqqgavOQoifBPBbsu8s/FMhRA/9g6Su2zr7+lhBCPEfAJwd/Hnd1hew11kIsRd9v2+z79nEiwBcFEJ8J7ZonQl8GMAT6AuALV1nIcQ9AO4C8D2DOQ1kqvN2cQHtB/A5KeWS9tvvAXirEOLLhRDfBOBbAfzpWKjLj98BfdbSlqyzEGKX9uebACwO/r8l6yulvCyl3CmlfImU8iXoM4NXSimXsUXrDABCiG/V/rwbwOcG/9/Kdf5eAHMA7pZS/qP2Kkudt4QFwMBbsd79AynlM0KIxwD8Gfqm1bSUcm0cxFUA21lLW7XOvyiE+A70TeAvAPhxYMv3MYktXud3CyFeBqCHfuTTTwBbvs4fQD/S59zA2vuUlPInctW5HAVRUFBQsE2xXVxABQUFBQUGigAoKCgo2KYoAqCgoKBgm6IIgIKCgoJtiiIACgoKCrYpigAouG4wOM31JuL3u4UQRwf/nxJCXBRCrAoh3uzI6ysGx2TcQLz7kCvtKCGEqAkhPiyE+LwQ4oIQ4mkhxJsG714rhPg7IUR9cCLkHwkh7tLSkm0hhNghhPiDcdSnYHOhCICC6x5Syt+TUr578OffADiI/k5RFw6jv1u6snjxwQ70lPQC/U19fySl/GYp5avQ39PyIu2zP5ZSvkJK+TL0T8j8gBDiewbvyLaQUl4B0BFC/M8p9BVc/ygCoGBTQQjxksH55782ONzto0KIr9Q+mR1otZeFEP9ykOagEOIDACCl/IKU8hL6m4Vc+GEA6pwVIYT4gBDiz4QQTwDYqdHzqoGlcEEI8aTadSyEeM2AvqdF/8z2RY2WM0KIxwE8JYT4qsEBdZ8ZaOpvHHx3wyDdZwb5/DhB4+sB/LOU8lfUD1LKL0op309VSErZAPDz6J+R42uL3xm0QcE2RhEABZsRLwPwQSnlPgD/DcB92rvnpJSvRP8AuAdjMhdCfBmAb5ZSfmHw05sGZe4F8GMA/qfBd88H8H4Abx5o36cAPDxIcxrAT0gpbwNgWhG3ob/7+vXoH9z2h1LK16B/T8MvCSG+CsDbAfzd4PfXAPixwZZ+HTcDuBhYvYsA/iXju/Pon5BbsI1RBEDBZsSzUso/Gfz/PwH4Lu3dbw3+vQDgJZH53wTgv2p/TwF4VEq5JqVsA/jDwe8vA7AH/W34DQD/F4AXif6tTF8tpfzk4DvT3XROSvlfBv9/A4Cjg/QfB/ACALsHv/9vg98/DeAb0D/PxQohxC8LIZpCiM+4PnPloWEFwAu9XxVsaWyXs4AKri+Y55Pof//T4N81xI/f/44+I3aVCfSZ6TMDLf/aj9fuG7DhS0YePyCl/HMjDwFgVkr5pCOfZwD8wJBAKacHi+DnHWleAeCzHvqAfv3/O+O7gi2MYgEUbEbsFkIopvs2AP9vzsyllH8L4AYhhBICf4T+yYo3DHz8rxv8/ucAdihahBDPF0LcPEj/90KIWwffvdVR3JPor1uIQR6v0H7/yYGbCUKIlw5cQzr+EMALRP94b4WvhAVCiH0AfgbALzvoUXgprp2aWrBNUQRAwWbEZwHcI4S4BODr0ff3szBYnF0C8BYAvyqEeMby6VO45lr6bQD/GcDlQVmfAAAp5T8DeDOA9wghmgAaGKwPoO/D/6AQ4mn0tfy/s5RzDMDzAVwaLBQfG/z+H9E/yfHi4PdfhWHRDE5w/V8BfLcQ4q+FEH8K4NfQPx5Y4XYVBoo+479fSrnAaIvXoX+efsE2RjkNtGBTQQjxEgBnpZR7Ki7nFQB+Skr5ryPT3yil/IfB/48C2CWlfCAnjVVCCPFHAN44sGYKtinKGkDBtoSUsi6E+JgQ4obIvQB3iv61hM9D/2z6g1kJrBBCiB0A/m1h/gXFAigoKCjYpihrAAUFBQXbFEUAFBQUFGxTFAFQUFBQsE1RBEBBQUHBNkURAAUFBQXbFP8/FX9mwh3ud+8AAAAASUVORK5CYII=\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",
"plt.plot(x, y, 'ko', markersize=0.5, alpha=0.5)\n",
"\n",
"plt.xlabel('phi1 (degree GD1)')\n",
"plt.ylabel('phi2 (degree GD1)');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is similar to what we saw in the previous lesson, except that have fewer stars now, because we did not find photometry data for all of the candidate sources."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Saving the DataFrame\n",
"\n",
"Let's save this `DataFrame` so we can pick up where we left off without running this query again.\n",
"The HDF file should already exist, so we'll add `candidate_df` to it."
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [],
"source": [
"filename = 'gd1_data.hdf'\n",
"\n",
"candidate_df.to_hdf(filename, 'candidate_df')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use `getsize` to confirm that the file exists and check the size:"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.5835609436035156"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from os.path import getsize\n",
"\n",
"MB = 1024 * 1024\n",
"getsize(filename) / MB"
]
},
{
"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.\n",
"\n",
"But before you go on, you might be interested in another file format, CSV."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## CSV\n",
"\n",
"Pandas can write a variety of other formats, [which you can read about here](https://pandas.pydata.org/pandas-docs/stable/user_guide/io.html).\n",
"We won't cover all of them, but one other important one is [CSV](https://en.wikipedia.org/wiki/Comma-separated_values), which stands for \"comma-separated values\".\n",
"\n",
"CSV is a plain-text format that can be read and written by pretty much any tool that works with data. In that sense, it is the \"least common denominator\" of data formats.\n",
"\n",
"However, it has an important limitation: some information about the data gets lost in translation, notably the data types. If you read a CSV file from someone else, you might need some additional information to make sure you are getting it right.\n",
"\n",
"Also, CSV files tend to be big, and slow to read and write.\n",
"\n",
"With those caveats, here's how to write one:"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [],
"source": [
"candidate_df.to_csv('gd1_data.csv')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can check the file size like this:"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.7606849670410156"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"getsize('gd1_data.csv') / MB"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see the first few lines like this:"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"def head(filename, n=3):\n",
" \"\"\"Print the first `n` lines of a file.\"\"\"\n",
" with open(filename) as fp:\n",
" for i in range(n):\n",
" print(next(fp))"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
",source_id,ra,dec,pmra,pmdec,best_neighbour_multiplicity,number_of_mates,g_mean_psf_mag,i_mean_psf_mag,phi1,phi2,pm_phi1,pm_phi2\n",
"\n",
"0,635860218726658176,138.5187065217173,19.09233926905897,-5.941679495793577,-11.346409129876392,1,0,17.8978004455566,17.5174007415771,-59.247329893833296,-2.016078400820631,-7.527126084640531,1.7487794924176672\n",
"\n",
"1,635674126383965568,138.8428741026386,19.031798198627634,-3.8970011609340207,-12.702779525389634,1,0,19.2873001098633,17.6781005859375,-59.13339098769217,-2.306900745179831,-7.560607655557415,-0.7417999555980248\n",
"\n"
]
}
],
"source": [
"head('gd1_data.csv')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The CSV file contains the names of the columns, but not the data types.\n",
"\n",
"We can read the CSV file back like this:"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [],
"source": [
"read_back_csv = pd.read_csv('gd1_data.csv')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's compare the first few rows of `candidate_df` and `read_back_csv`"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" Unnamed: 0 source_id ra dec pmra pmdec \\\n",
"0 0 635860218726658176 138.518707 19.092339 -5.941679 -11.346409 \n",
"1 1 635674126383965568 138.842874 19.031798 -3.897001 -12.702780 \n",
"2 2 635535454774983040 137.837752 18.864007 -4.335041 -14.492309 \n",
"\n",
" best_neighbour_multiplicity number_of_mates g_mean_psf_mag \\\n",
"0 1 0 17.8978 \n",
"1 1 0 19.2873 \n",
"2 1 0 16.9238 \n",
"\n",
" i_mean_psf_mag phi1 phi2 pm_phi1 pm_phi2 \n",
"0 17.517401 -59.247330 -2.016078 -7.527126 1.748779 \n",
"1 17.678101 -59.133391 -2.306901 -7.560608 -0.741800 \n",
"2 16.478100 -59.785300 -1.594569 -9.357536 -1.218492 "
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"read_back_csv.head(3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that the index in `candidate_df` has become an unnamed column in `read_back_csv`. The Pandas functions for writing and reading CSV files provide options to avoid that problem, but this is an example of the kind of thing that can go wrong with CSV files."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Best practices\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",
"* For most applications, saving data in FITS or HDF5 is better than CSV. FITS and HDF5 are binary formats, so the files are usually smaller, and they store metadata, so you don't lose anything when you read the file back.\n",
"\n",
"* On the other hand, CSV is a \"least common denominator\" format; that is, it can be read by practically any application that works with 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
}
![](\"https://github.com/datacarpentry/astronomy-python/raw/gh-pages/fig/gd1-3.png\")
\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": 28,
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Getting photometry data\n",
"\n",
"The Gaia dataset contains some photometry data, including the variable `bp_rp`, which contains BP-RP color (the difference in mean flux between the BP and RP bands).\n",
"We use this variable to select stars with `bp_rp` between -0.75 and 2, which excludes many class M dwarf stars.\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": [
"## The best neighbor table\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. Use the `JOIN` operator to look up each `source_id` in the `panstarrs1_best_neighbour` table, which contains the `obj_id` of the best match for each star in the Gaia catalog; then\n",
"\n",
"2. 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",
"Before we get to the `JOIN` operation, let's explore these tables.\n",
"Here's the metadata for `panstarrs1_best_neighbour`."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Retrieving table 'gaiadr2.panstarrs1_best_neighbour'\n",
"Parsing table 'gaiadr2.panstarrs1_best_neighbour'...\n",
"Done.\n"
]
}
],
"source": [
"from astroquery.gaia import Gaia\n",
"\n",
"meta = Gaia.load_table('gaiadr2.panstarrs1_best_neighbour')"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"TAP Table name: gaiadr2.gaiadr2.panstarrs1_best_neighbour\n",
"Description: Pan-STARRS1 BestNeighbour table lists each matched Gaia object with its\n",
"best neighbour in the external catalogue.\n",
"There are 1 327 157 objects in the filtered version of Pan-STARRS1 used\n",
"to compute this cross-match that have too early epochMean.\n",
"Num. columns: 7\n"
]
}
],
"source": [
"print(meta)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here are the columns."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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": [
"for column in meta.columns:\n",
" print(column.name)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the [documentation for these variables](https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_crossmatches/ssec_dm_panstarrs1_best_neighbour.html) .\n",
"\n",
"The ones we'll use are:\n",
"\n",
"* `source_id`, which we will match up with `source_id` in the Gaia table.\n",
"\n",
"* `number_of_neighbours`, which indicates how many sources in Pan-STARRS are matched with this source in Gaia.\n",
"\n",
"* `number_of_mates`, which indicates the number of *other* sources in Gaia that are matched with the same source in Pan-STARRS.\n",
"\n",
"* `original_ext_source_id`, which we will match up with `obj_id` in the Pan-STARRS table.\n",
"\n",
"Ideally, `number_of_neighbours` should be 1 and `number_of_mates` should be 0; in that case, there is a one-to-one match between the source in Gaia and the corresponding source in Pan-STARRS.\n",
"\n",
"Here's a query that selects these columns and returns the first 5 rows."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"query = \"\"\"SELECT \n",
"TOP 5\n",
"source_id, number_of_neighbours, number_of_mates, original_ext_source_id\n",
"FROM gaiadr2.panstarrs1_best_neighbour\n",
"\"\"\""
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
}
],
"source": [
"job = Gaia.launch_job_async(query=query)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"Table length=5\n",
"![](\"https://github.com/datacarpentry/astronomy-python/raw/gh-pages/fig/join.png\")
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Joining tables\n",
"\n",
"Now let's get to the details of performing a `JOIN` operation.\n",
"As a starting place, let's go all the way back to the cone search from Lesson 2."
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"query_cone = \"\"\"SELECT \n",
"TOP 10 \n",
"source_id\n",
"FROM gaiadr2.gaia_source\n",
"WHERE 1=CONTAINS(\n",
" POINT(ra, dec),\n",
" CIRCLE(88.8, 7.4, 0.08333333))\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And let's run it, to make sure we have a working query to build on."
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO: Query finished. [astroquery.utils.tap.core]\n"
]
}
],
"source": [
"from astroquery.gaia import Gaia\n",
"\n",
"job = Gaia.launch_job_async(query=query_cone)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"Table length=10\n",
"