From ef86df308452b4f47a51ee30c1d7217036c95792 Mon Sep 17 00:00:00 2001 From: Julien Date: Fri, 8 Feb 2019 13:23:20 +0100 Subject: [PATCH 1/3] Add supplementary notebook for 1807.03078 --- Spark4Physicists.ipynb | 745 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 745 insertions(+) create mode 100644 Spark4Physicists.ipynb diff --git a/Spark4Physicists.ipynb b/Spark4Physicists.ipynb new file mode 100644 index 0000000..4c79200 --- /dev/null +++ b/Spark4Physicists.ipynb @@ -0,0 +1,745 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Analysing billion-objects catalogue interactively: Spark for physicists\n", + "
Authors: **S. Plaszczynski, J. Peloton, C. Arnault and J.E. Campagne **\n", + "
***LAL, Univ. Paris-Sud, CNRS/IN2P3, Universite Paris-Saclay, Orsay, France***\n", + "
Last Verified to Run: **2019-02-07**\n", + "\n", + "\n", + "\n", + "This jupyter notebook describes the analysis presented in https://arxiv.org/abs/1807.03078\n", + "The results presented here were run on the Universite Paris-Sud Spark cluster: https://www.informatique-scientifique.u-psud.fr/services/spark.html\n", + "\n", + "In order to run this notebook you need:\n", + "- a cluster running `pypark` (spark version >=2.2) on which the following python libraries have been installed (on each executor):\n", + "`matplotlib`,`numpy`,`pandas`,`pyarrow`,`healpy`\n", + "- a spark jupyter kernel (NERSC users: https://github.com/astrolabsoftware/spark-kernel-nersc)\n", + "- the `SparkFITS` library : https://astrolabsoftware.github.io/spark-fits/\n", + "- a set of FITS files, as produced for instance by the `CoLoRe` simulation : https://github.com/damonge/CoLoRe\n", + "The directory that points to these files should be put in the FITSDIR environment variable. Any file(s) with (RA,DEC,redshift) coordinates may also be used with little adaptation to match the exact key names.\n", + "\n", + "For testing, Spark (and these dependencies) can also be installed on a personal computer https://spark.apache.org/docs/latest/.\n", + "\n", + "For LSST/DESC collaborators with a NERSC account see https://github.com/LSSTDESC/desc-spark (FITSDIR=/global/cscratch1/sd/plaszczy/LSST10Y)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Initialisation" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Spark session started\n" + ] + } + ], + "source": [ + "from pyspark.sql import SparkSession\n", + "from pyspark import StorageLevel\n", + "from pyspark.sql import functions as F\n", + "from pyspark.sql.functions import randn\n", + "from pyspark.sql.types import IntegerType,FloatType\n", + "from pyspark.sql.functions import pandas_udf, PandasUDFType\n", + "\n", + "spark = SparkSession.builder.getOrCreate()\n", + "sc=spark.sparkContext\n", + "print(\"Spark session started\")\n", + "\n", + "#a class for timing\n", + "from time import time\n", + "class Timer:\n", + " \"\"\"\n", + " a simple class for printing time (s) since last call\n", + " \"\"\"\n", + " def __init__(self):\n", + " self.t0=time()\n", + " \n", + " def start(self):\n", + " self.t0=time()\n", + " \n", + " def stop(self):\n", + " t1=time()\n", + " print(\"{:2.1f}s\".format(t1-self.t0))\n", + "\n", + "timer=Timer()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Reading the data \n", + "\n", + "all CoLoRe FITS files assumed to be in a directory pointed by the FITSDIR env variable\n", + "Le us build the **gal** datafreame selecting only the \"RA\" and \"dec\" columns and building a new redshift column called \"z\" on the flight." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "6.1s\n" + ] + } + ], + "source": [ + "import os\n", + "dirfits=os.environ['FITSDIR']\n", + "timer.start()\n", + "gal=spark.read.format(\"fits\").option(\"hdu\",1)\\\n", + " .load(dirfits)\\\n", + " .select(F.col(\"RA\"), F.col(\"Dec\"), (F.col(\"Z_COSMO\")+F.col(\"DZ_RSD\")).alias(\"z\"))\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Pretty fast isn't it? This is because of lazy evaluation : data have not been actaully loaded (yet). The pipeline (Direct Acyclic Graph, DAG) has been updated but its execution has not been triggered since read/load/slect are not actions but transformations.\n", + "At this level the schema of the data is however available so you can looks at it" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "root\n", + " |-- RA: float (nullable = true)\n", + " |-- Dec: float (nullable = true)\n", + " |-- z: float (nullable = true)\n", + "\n" + ] + } + ], + "source": [ + "gal.printSchema()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "you may want to look at a few values to check eveything is OK. Let's print the first 5" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "+---------+---------+---------+\n", + "| RA| Dec| z|\n", + "+---------+---------+---------+\n", + "|225.80168|18.519966|2.4199903|\n", + "|225.73839|18.588171|2.4056022|\n", + "|225.79999|18.635067| 2.396816|\n", + "|225.49783|18.570776|2.4139786|\n", + "|225.57983|18.638515|2.3995044|\n", + "+---------+---------+---------+\n", + "only showing top 5 rows\n", + "\n", + "14.1s\n" + ] + } + ], + "source": [ + "timer.start()\n", + "gal.show(5)\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Still pretty fast. Here show is an action and has triggered the data access.But not the full data set : since we required only the first 5 rows, the DAG was (automatically) analyzed and only the first block was actually read" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we add an extra column that performs gaussian smearing on the \"z\" column. We use the Spark \"rand\" function that is highly optimized" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "+---------+---------+---------+---------+\n", + "| RA| Dec| z| zrec|\n", + "+---------+---------+---------+---------+\n", + "|225.80168|18.519966|2.4199903|2.3029215|\n", + "|225.73839|18.588171|2.4056022|2.6602228|\n", + "|225.79999|18.635067| 2.396816|2.3686476|\n", + "|225.49783|18.570776|2.4139786|2.4345846|\n", + "|225.57983|18.638515|2.3995044|2.3599765|\n", + "+---------+---------+---------+---------+\n", + "only showing top 5 rows\n", + "\n", + "8.3s\n" + ] + } + ], + "source": [ + "timer.start()\n", + "gal=gal.withColumn(\"zrec\",(gal.z+0.03*(1+gal.z)*randn()).astype('float'))\n", + "gal.show(5)\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we have all our columns ready, we put data in cache. Recall that the amount of memory available for the cache is about 60% of your total executors memory. You may simply use the cache() function. Here in case you do not have enough memory we use a slighlty more evolved form that fills the cache and spills the rest over disk if necessary.\n", + "Remember to finalize your command with an actiom as count() here.\n", + "This is the lenghty part but you need to perform it only once." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "#galaxies=5926764520\n", + "97.6s\n" + ] + } + ], + "source": [ + "timer.start()\n", + "print(\"#galaxies={}\".format(gal.persist(StorageLevel.MEMORY_AND_DISK).count()))\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We put about 6 billions of galaxies in cache, just as if we had a hge RAM. All further computattions will be boosted." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Basic statistics" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Le us first look at some simple statistics on a single column" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "+-------+-------------------+\n", + "|summary| z|\n", + "+-------+-------------------+\n", + "| count| 5926764520|\n", + "| mean| 0.875229446764415|\n", + "| stddev|0.47360539212211805|\n", + "| min| -5.93947E-4|\n", + "| max| 2.4352543|\n", + "+-------+-------------------+\n", + "\n", + "5.6s\n" + ] + } + ], + "source": [ + "timer.start()\n", + "gal.describe(['z']).show()\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On all the columns:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "+-------+------------------+--------------------+------------------+-------------------+\n", + "|summary| RA| Dec| z| zrec|\n", + "+-------+------------------+--------------------+------------------+-------------------+\n", + "| count| 5926764520| 5926764520| 5926764520| 5926764520|\n", + "| mean|179.99419917017775|-0.00314190695924...|0.8752294467644149| 0.8752288274953041|\n", + "| stddev|103.92702755739798| 39.17648956108155|0.4736053921221175|0.47714711978268487|\n", + "| min| 4.9597503E-8| -89.997185| -5.93947E-4| -0.124362685|\n", + "| max| 360.0| 89.9986| 2.4352543| 2.95336|\n", + "+-------+------------------+--------------------+------------------+-------------------+\n", + "\n", + "9.2s\n" + ] + } + ], + "source": [ + "timer.start()\n", + "gal.describe().show()\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you just need some specific values it will be more efficient to call them directly" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "min=-0.0005939470138400793 max=2.4352543354034424\n", + "4.5s\n" + ] + } + ], + "source": [ + "timer.start()\n", + "minmax=gal.select(F.min(\"z\"),F.max(\"z\")).first()\n", + "print(\"min={} max={}\".format(minmax[0],minmax[1]))\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Histograms" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The distributed way of building histograms is:\n", + "- add a column of bin values\n", + "- groupBy this column and count the number of elements in each group\n", + "- you may then drop the bin and add the center of the bin locations." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## with Spark native functions\n", + "Here is how to do it using Spark dataframe functions" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "#histogram caracteristics\n", + "zmin=minmax[0]\n", + "zmax=minmax[1]\n", + "Nbins=100\n", + "dz=(zmax-zmin)/Nbins\n", + "\n", + "#create the bin number column\n", + "zbin=gal.select(gal.z,((gal['z']-zmin-dz/2)/dz).cast(IntegerType()).alias('bin'))\n", + "# groupBy it + count + sort by bin values\n", + "h=zbin.groupBy(\"bin\").count().orderBy(F.asc(\"bin\"))\n", + "# add the bin locations (zbin) and drop the bin numer\n", + "h=h.select(\"bin\",(zmin+dz/2+h['bin']*dz).alias('loc'),\"count\").drop(\"bin\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now data is reduced to a simple histogram, we go back to the standard python world by converting the dataframe to a pandas object (this is the action so triggers the DAG materialization)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Index(['loc', 'count'], dtype='object') size=100\n", + "13.0s\n" + ] + } + ], + "source": [ + "timer.start()\n", + "p=h.toPandas()\n", + "print(p.columns,\"size={}\".format(p.index.size))\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us plot the result with matplotlib" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAERCAYAAACAbee5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEuZJREFUeJzt3W2MZOdZ5vH/xcR2JCcEwQwi8oynDRjBECU4aTmBwGJIFtn+YH8goLEIYOQw4sXhJQjkBeRE5hNBu2jZdchOwDJEEBMCggFNMIIYBQGO3CHB2I4MgzONW05w581ZZHDw7s2HrrbL7eqpM91Vdaqe+v+kVqrqnKm5z9Tk8j33OeepVBWSpLZ8Ud8FSJImz3CXpAYZ7pLUIMNdkhpkuEtSgwx3SWpQr+Ge5I4kjyd5oMO+lya5J8lHktyf5NpZ1ChJi6jvzv1O4OqO+/488N6qugI4DrxjWkVJ0qLrNdyr6oPAZ4ZfS/JVSf4kyYeT/GWSr93eHfjiweOXAI/NsFRJWigv6LuAEU4CP1RV/5jk1Wx16N8OvA340yRvBi4GXt9fiZI03+Yq3JO8CPgm4HeTbL980eB/bwDurKr/nuQbgXcneVlV/f8eSpWkuTZX4c7WmOhzVfUNI7bdxGA+X1V/k+SFwEHg8RnWJ0kLoe8Tqs9RVZ8HPp7kuwCy5RWDzf8MvG7w+tcBLwQ2eylUkuZc+lwVMsl7gKvY6sD/BXgr8AHgV4GXAhcAd1XVbUmOAe8CXsTWydWfqao/7aNuSZp3vYa7JGk65mosI0majN5OqB48eLBWVlb6+u0laSF9+MMf/lRVHRq3X2/hvrKywtraWl+/vSQtpCTrXfZzLCNJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3BfcysoKSUjCxRdfPPLxzucu+yC1b96+rEMdrKyssL7+7B3I2yt7Jhn5eNQ2SW2zc19A6+vrVBV7Xa75oosusqOXGmfnvoSeeuqpXbt9SW2wc9czhjt6u3hpsdm5L4idc/Zp2NnRS1pcdu4LYr9zdknLxXCXpAYZ7nNs+Br2WXP+Li02Z+5zbHsUA7OfgTt/lxabnbskNchw11iOaKTF41hmzsziksfz5YhGWjx27nPGSx4lTYLhLkkNGhvuSe5I8niSB3bZ/j1J7h/8/HWSV0y+TM0L5+/SYujSud8JXH2O7R8HvrWqXg78AnByAnVpTm3P36tq7s4NSHrW2BOqVfXBJCvn2P7XQ0/vBQ7vvyxJ0n5MeuZ+E/D+3TYmOZFkLcna5ubmhH/rxdXnnaj74YhGml8TuxQyybexFe7fvNs+VXWSwdhmdXXVy0EG+rwTdT+8RFKaXxMJ9yQvB34NuKaqPj2J95Qk7d2+xzJJLgV+H/jeqvqH/ZckSdqvsZ17kvcAVwEHk2wAbwUuAKiqdwK3Al8GvGPwT/Onq2p1WgVLksbrcrXMDWO2vwl408QqWhLzuMzAfmyfXAU4evQoZ8+e7bcgacm5tkxPFvUk6m48uSrNF5cfkKQGGe6S1CDDXRPnzU1S/5y5a+Kcv0v9s3OfoUVdZkDS4jHcZ8gv4pA0K4a7psr5u9QPZ+6aKufvUj/s3CWpQYa7JDXIcNfMOH+XZseZu2bG+bs0O3buktQgw12SGmS4S1KDDPcpc8mB0Ty5Kk2XJ1SnrLUv5ZgUT65K02XnLkkNMtwlqUGGu3rn/F2aPGfu6p3zd2nyxnbuSe5I8niSB3bZniS/kuRMkvuTvHLyZUqSzkeXscydwNXn2H4NcPng5wTwq/svS5K0H2PDvao+CHzmHLtcD/xmbbkX+JIkL51UgYvIa9sl9W0SJ1QvAR4der4xeG1p+XV6kvo2iXAf1Z6OTLUkJ5KsJVnb3NycwG+t1njljDQZk7haZgM4MvT8MPDYqB2r6iRwEmB1ddW2Vs/jlTPSZEyicz8FfN/gqpnXAE9U1Scm8L6SpD0a27kneQ9wFXAwyQbwVuACgKp6J3AauBY4AzwJ/MC0ipUkdTM23KvqhjHbC/jRiVUkDWzP3wGOHj3K2bNn+y1IWiDeoaq55fxd2jvXlpGkBhnuktQgw12SGmS4ayF4c5N0fjyhOiErKyusr6/3XUazPLkqnR879wlxPRlJ88Rwl6QGGe6S1CDDXQvHk6vSeJ5Q1cLx5Ko0np27JDXIcJekBhnuktQgw10LzZOr0mieUN0H70rtnydXpdHs3PfBu1IlzSvDXZIaZLhLUoMMd0lqkOGuZnjljPQsr5ZRM7xyRnqWnbskNahTuCe5OsnDSc4kuWXE9kuT3JPkI0nuT3Lt5EuVJHU1NtyTHABuB64BjgE3JDm2Y7efB95bVVcAx4F3TLpQ6Xw4f9ey6zJzvxI4U1WPACS5C7geeGhonwK+ePD4JcBjkyxSOl/O37XsuoxlLgEeHXq+MXht2NuANybZAE4Dbx71RklOJFlLsra5ubmHcvu3srLyTEcoSfOqS7iPSrGd99vfANxZVYeBa4F3J3nee1fVyapararVQ4cOnX+1c8AlByQtgi7hvgEcGXp+mOePXW4C3gtQVX8DvBA4OIkCJUnnr0u43wdcnuSyJBeydcL01I59/hl4HUCSr2Mr3Bdz7qLmeHJVy2jsCdWqejrJzcDdwAHgjqp6MMltwFpVnQJ+CnhXkp9ka2RzYzm30Jzw5KqWUac7VKvqNFsnSodfu3Xo8UPAaydbmiRpr7xDVZIaZLhLUoMMdy2V4ZOrnmBVy1wVUktl+OQqeIJV7bJzl6QGGe6S1CDDXUvNG5zUKmfuWmre4KRW2bl34EqQkhaN4d6BK0FKWjSGuyQ1yHCXBjy5qpZ4QlUa8OSqWmLnLkkNMtwlqUGGuzSC83ctOmfu0gjO37Xo7NwlqUGGuyQ1yHCXpAYZ7tIYnlzVIvKEqjSGJ1e1iDp17kmuTvJwkjNJbtlln+9O8lCSB5P89mTLnD1XgpS0yMZ27kkOALcD/xXYAO5LcqqqHhra53LgvwGvrarPJvnyaRU8K9srQYLdmqTF06VzvxI4U1WPVNUXgLuA63fs84PA7VX1WYCqenyyZUrzwfm7FkWXmfslwKNDzzeAV+/Y52sAkvwVcAB4W1X9yUQqlOaI83ctii7hPupv8M5vrXgBcDlwFXAY+MskL6uqzz3njZITwAmASy+99LyLlSR102UsswEcGXp+GHhsxD5/WFX/UVUfBx5mK+yfo6pOVtVqVa0eOnRorzVLksboEu73AZcnuSzJhcBx4NSOff4A+DaAJAfZGtM8MslCpXnj/F3zbOxYpqqeTnIzcDdb8/Q7qurBJLcBa1V1arDtO5I8BPw/4Ker6tPTLFzqm/N3zbNONzFV1Wng9I7Xbh16XMBbBj+SpJ65/IAkNchwlybA+bvmjWvLSBPg/F3zxs5dkhpkuEtSgwz3Ia4EqUlw/q554Mx9iCtBahKcv2se2LlLUoMMd2mKHNGoL45lpClyRKO+2LlLUoMMd0lqkOEuzYjzd82SM3dpRpy/a5bs3CWpQYa7JDXIcJd64Pxd0+bMXeqB83dNm527JDVo6cPdlSDVN0c0moalH8u4EqT65ohG07D0nbsktchwl+aIIxpNSqdwT3J1koeTnElyyzn2e0OSSrI6uRKl5bE9oqkq1tfX+y5HC2xsuCc5ANwOXAMcA25IcmzEfi8Gfgz40KSLlCSdny6d+5XAmap6pKq+ANwFXD9iv18A3g78+wTrkyTtQZdwvwR4dOj5xuC1ZyS5AjhSVX98rjdKciLJWpK1zc3N8y5WWibO37UfXS6FHHVtVj2zMfki4JeBG8e9UVWdBE4CrK6u1pjdpaXmJZLajy6d+wZwZOj5YeCxoecvBl4G/EWSs8BrgFOeVJWk/nQJ9/uAy5NcluRC4DhwantjVT1RVQeraqWqVoB7geuqam0qFUtLyBGNztfYsUxVPZ3kZuBu4ABwR1U9mOQ2YK2qTp37HSTtlyMana9Oyw9U1Wng9I7Xbt1l36v2X5ak3Wx38QBHjx7l7Nmz/RakubT0a8tIi8YuXl0s5fIDrgQpqXVLGe7bK0Fudz/SovJEq3bjWEZaYI5otJul7NwlqXWGu9QIRzQa5lhGaoQjGg2zc5ekBhnuUoMc0cixjNQgRzSyc5caZxe/nOzcpcbZxS8nO3dJapDhLi0RRzTLw7GMtEQc0SyPpencXQlSei67+LYtTee+vRIk2LFIYBffuqXp3CXtbriLt5Nvw9J07pJ2N9zFg518C+zcJT2P8/jFZ+cu6Xmcxy8+O3dJ52QXv5g6hXuSq5M8nORMkltGbH9LkoeS3J/kz5McnXypkvqw3cVXFevr632Xo47GhnuSA8DtwDXAMeCGJMd27PYRYLWqXg68D3j7pAuV1D+7+MXRZeZ+JXCmqh4BSHIXcD3w0PYOVXXP0P73Am+cZJGS5oOz+MXRZSxzCfDo0PONwWu7uQl4/6gNSU4kWUuytrm52b1KSXPHLn6+dQn3Uf95rhGvkeSNwCrwS6O2V9XJqlqtqtVDhw51r1LS3BmexX/yk5806OdMl7HMBnBk6Plh4LGdOyV5PfBzwLdW1VOTKU/SInBcM3+6dO73AZcnuSzJhcBx4NTwDkmuAP4PcF1VPT75MvfGxcKk2XNcMx/GhntVPQ3cDNwNfAx4b1U9mOS2JNcNdvsl4EXA7yb5aJJTu7zdTG0vFjZ8W7Wk6XJcMx863aFaVaeB0zteu3Xo8esnXJekBjiu6Y93qEqaCcc1s2W4S5oJxzWz5cJhkmbOcc302blL6pXjmukw3CX1ynHNdBjukuaGQT85hrukuWTQ709z4e5dqVJ7DPrz11y4e1eq1DaDvpvmwl3S8jDod2e4S2qCQf9chruk5hj0hrukxi1r0BvukpbGMgW94S5pKe0W9BdffHEToW+4S1p6w0H/5JNPjgz9RQt7w12SdjEc+ovW4S98uA/fkepdqZKmqUuHPy+hv/DhPnxHqnelSupDl9CfddD7ZR2SNCV9finJwnfukqTnM9wlqUGGuyQ1qFO4J7k6ycNJziS5ZcT2i5L8zmD7h5KsTLpQSVJ3Y8M9yQHgduAa4BhwQ5JjO3a7CfhsVX018MvAL0660GF+IYcknVuXzv1K4ExVPVJVXwDuAq7fsc/1wG8MHr8PeF2mmLx+IYekRXPRRRfN9LLILpdCXgI8OvR8A3j1bvtU1dNJngC+DPjU8E5JTgAnBk//NcnDeykaOJjkmffe+d+R4ecNPj4IfOp8f+0c1L3fx7t+5vt5PCfHNpXPvIE/l6l85n3+uWxbX18fN3k4yI78HHL0XL9wW5dwH1XBzpa5yz5U1UngZIff89wFJWtVtbrf91lEy3rsy3rcsLzHvqzHDZM59i5jmQ3gyNDzw8Bju+2T5AXAS4DP7KcwSdLedQn3+4DLk1yW5ELgOHBqxz6ngO8fPH4D8IFyIC5JvRk7lhnM0G8G7gYOAHdU1YNJbgPWquoU8OvAu5OcYatjPz7NopnAaGeBLeuxL+txw/Ie+7IeN0xifG2DLUnt8Q5VSWqQ4S5JDZrrcM8SL3vQ4dhvTLKZ5KODnzf1UeekJbkjyeNJHthle5L8yuDP5f4kr5x1jdPQ4bivSvLE0Od966xrnIYkR5Lck+RjSR5M8uMj9mn1M+9y7Hv/3Ie/6GKeftg6eftPwFcCFwJ/Bxzbsc+PAO8cPD4O/E7fdc/w2G8E/nfftU7h2P8L8ErggV22Xwu8n617K14DfKjvmmd03FcBf9x3nVM47pcCrxw8fjHwDyP+rrf6mXc59j1/7vPcuc/dsgcz1OXYm1RVH+Tc90hcD/xmbbkX+JIkL51NddPT4bibVFWfqKq/HTz+v8DH2LrjfVirn3mXY9+zeQ73Ucse7Dzw5yx7AGwve7Douhw7wHcO/pn6viRHRmxvUdc/mxZ9Y5K/S/L+JF/fdzGTNhirXgF8aMem5j/zcxw77PFzn+dwn9iyBwuoy3H9EbBSVS8H/oxn/wXTulY/83H+FjhaVa8A/hfwBz3XM1FJXgT8HvATVfX5nZtH/JJmPvMxx77nz32ew32Zlz0Ye+xV9emqemrw9F3Aq2ZUW9+6/L1oTlV9vqr+dfD4NHBBkoM9lzURSS5gK9x+q6p+f8QuzX7m4459P5/7PIf7Mi97MPbYd8wcr2NrXrcMTgHfN7iC4jXAE1X1ib6LmrYkX7F9PinJlWz9f/fT/Va1f4Nj+nXgY1X1P3bZrcnPvMux7+dz77IqZC9qPpc9mImOx/5jSa4Dnmbr2G/sreAJSvIetq4QOJhkA3grcAFAVb0TOM3W1RNngCeBH+in0snqcNxvAH44ydPAvwHHG2lkXgt8L/D3ST46eO1ngUuh7c+cbse+58/d5QckqUHzPJaRJO2R4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLA0l+aGjd7I8nuafvmqS98iYmaYfBeh8fAN5eVX/Udz3SXti5S8/3P9lap8hg18Ka27VlpD4kuRE4CtzccynSvjiWkQaSvIqtdfG/pao+23c90n44lpGedTPwpcA9g5Oqv9Z3QdJe2blLUoPs3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatB/Aq3INiHIhMzLAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "plt.bar(p['loc'].values,p['count'].values,dz,color='white',edgecolor='black')\n", + "plt.xlabel(\"z\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## With a User Defined Function (udf)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Although absolutely unnecessary here, we illustrate how one could have plugged an external UDF (to compute the bin number. Performances are extrememly dowgraded so you should NOT do it this way (see next part)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "111.3s\n" + ] + } + ], + "source": [ + "#define your UDF (here just binning)\n", + "binNumber_udf=F.udf(lambda z: int((z-zmin)/dz))\n", + "timer.start()\n", + "#create the bin number column\n", + "p_udf=gal.select(gal.z,binNumber_udf(gal.z).alias('bin')).groupBy(\"bin\").count().orderBy(F.asc(\"bin\")).toPandas()\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "That's not great.Let see how to be more efficient" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## With pandas_udf" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The recommended way to bind Spark to external python functions is through pandas_udf which achieve some level of vectorization the previous method miss." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "46.1s\n" + ] + } + ], + "source": [ + "import pandas as pd\n", + "@pandas_udf(\"float\", PandasUDFType.SCALAR)\n", + "def binFloat(z):\n", + " return pd.Series((z-zmin)/dz)\n", + "timer.start()\n", + "p_udf=gal.select(gal.z,binFloat(\"z\").astype('int').alias('bin')).groupBy(\"bin\").count().orderBy(F.asc(\"bin\")).toPandas()\n", + "timer.stop()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tomography\n", + "This is a popular method in cosmology when one builds \"shells\" in some redshift range of over-densities and compute their 2D power-spectrum (using some spherical pixelization as HealPix). As the redshift decreases, matter clusters more and more and it is expected that these spectra increase in amplitude. Some cross-correlations between nearby bins also appear because of the smearing between the true redhifts (\"z\") and the observed ones (\"zrec\")." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# define a function that filters galaxies within a redshift range and prject them onto a healpix map\n", + "# that is returned\n", + "\n", + "import numpy as np\n", + "import healpy as hp\n", + "\n", + "nside=512\n", + "\n", + "def redshift_bin(gal,z1,z2): \n", + " timer.start()\n", + " #filtering\n", + " shell=gal.filter(gal['zrec'].between(z1,z2))\n", + " print(\"shell=[{},{}] N={}M\".format(z1,z2,shell.count()/1e6))\n", + " \n", + " #projection into healpixels: build a pandas_udf that transforms (ra,dec) into (theta,phi) \n", + " # and calls healpy ang2pix function \n", + " @pandas_udf('int', PandasUDFType.SCALAR)\n", + " def Ang2Pix(ra,dec):\n", + " return pd.Series(hp.ang2pix(nside,np.radians(90-dec),np.radians(ra)))\n", + " map=shell.select(Ang2Pix(\"RA\",\"Dec\").alias(\"ipix\")).groupBy(\"ipix\").count().toPandas()\n", + " \n", + " #back to python world\n", + " myMap = np.zeros(hp.nside2npix(nside))\n", + " myMap[map['ipix'].values]=map['count'].values\n", + " timer.stop()\n", + " return myMap" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Do a loop over redshift bins (choosen here so that there are some overlap between adajacent ones for the true redshift) and store the Healpix maps with the number of galaxies in each pixel in a dictionary for later use." + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "shell=[0.0,0.13] N=68.557105M\n", + "27.8s\n", + "shell=[0.13,0.27] N=322.103946M\n", + "24.7s\n", + "shell=[0.27,0.43] N=669.646368M\n", + "30.3s\n", + "shell=[0.43,0.63] N=1049.274954M\n", + "35.3s\n", + "shell=[0.63,0.82] N=980.216568M\n", + "37.5s\n", + "shell=[0.82,1.05] N=984.086242M\n", + "33.9s\n", + "shell=[1.05,1.32] N=812.257315M\n", + "30.7s\n", + "shell=[1.32,1.61] N=529.298439M\n", + "30.5s\n", + "shell=[1.61,1.95] N=325.212359M\n", + "27.8s\n", + "shell=[1.95,2.32] N=156.932408M\n", + "26.4s\n" + ] + }, + { + "data": { + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#redshift bins\n", + "zshell=[0.0,0.13,0.27,0.43,0.63,0.82,1.05,1.32,1.61,1.95,2.32]\n", + "tmap=dict()\n", + "\n", + "#loop over bins and store healpix maps in a dictionary indexed by z1\n", + "plt.figure(figsize=(15,10))\n", + "for (z1,z2) in zip(zshell[:-1],np.roll(zshell,-1)): \n", + " tmap[z1]=redshift_bin(gal,z1,z2) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now do a standard analysis as the density contrast power spectra with (standard) healpix functions" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "bin=[0.0,0.13] mean #gals=21.7937167485555\n", + "bin=[0.13,0.27] mean #gals=102.39408683776855\n", + "bin=[0.27,0.43] mean #gals=212.87484741210938\n", + "bin=[0.43,0.63] mean #gals=333.5555248260498\n", + "bin=[0.63,0.82] mean #gals=311.60245513916016\n", + "bin=[0.82,1.05] mean #gals=312.83259137471515\n", + "bin=[1.05,1.32] mean #gals=258.209646542867\n", + "bin=[1.32,1.61] mean #gals=168.25944232940674\n", + "bin=[1.61,1.95] mean #gals=103.38222471872966\n", + "bin=[1.95,2.32] mean #gals=49.88746897379557\n" + ] + }, + { + "data": { + "text/plain": [ + "Text(0,0.5,'power spectrum')" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEKCAYAAAAiizNaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8lOW1x78nkxWSsCZACPuWEIRIWNyKCkVxpRT0yqUtVVzqpa2t9V6r7a3XWlutF+tu61brilr0ShFxi4oiGIOAhNUIEcIalhBIyDZ57h/vzGQymeXNMpmQnO/nMx8mz/tsk2hOznnOc35ijEFRFEVR2oqoSG9AURRF6Vyo4VEURVHaFDU8iqIoSpuihkdRFEVpU9TwKIqiKG2KGh5FURSlTVHDoyiKorQpangURVGUNkUNj6IoitKmREd6A5Gkd+/eZvDgwZHehqIoyinF2rVrDxljUpo7vlMbnsGDB5Ofnx/pbSiKopxSiMi3LRkf1lCbiMwQkW0iUigiv/bzPE5EXnE9/1xEBns9u83Vvk1ELvRqf0ZEDopIgc9c94nIVhH5SkTeEJHu4fxsiqIoSvMIm+EREQfwKHARMBqYKyKjfbotAI4aY4YDfwHudY0dDVwFZAEzgMdc8wE862rz5T1gjDFmLLAduK1VP5CiKIrSKoTT45kEFBpjdhhjqoHFwEyfPjOBf7je/xOYJiLial9sjKkyxuwECl3zYYxZCRzxXcwY864xptb15RogvbU/kKIoitJywnnG0x/Y7fV1MTA5UB9jTK2IHAN6udrX+Izt34S1rwFe8fdARK4HrgcYOHBgE6ZUFKWp1NTUUFxcTGVlZaS3ojSD+Ph40tPTiYmJadV5w2l4xE+br/hPoD52xvpfVOQ3QC3wor/nxpgngCcAJkyYoGJEihJGiouLSUpKYvDgwVjBDOVUwRjD4cOHKS4uZsiQIa06dzhDbcXAAK+v04G9gfqISDTQDSuMZmdsI0RkPnApMM+owp2iRJzKykp69eqlRucURETo1atXWLzVcBqeL4ARIjJERGKxkgWW+vRZCsx3vZ8D5LoMxlLgKlfW2xBgBJAXbDERmQHcClxujKloxc+hKEoLUKNz6hKun13YDI/roP+nwDvAFuBVY8wmEfm9iFzu6vY00EtECoGbgV+7xm4CXgU2AyuAhcYYJ4CIvAysBkaJSLGILHDN9QiQBLwnIutF5K+h9lhZ42ylT6soiqLYRTpzRKr34ExzqGhLpLehKB2WLVu2kJmZGeltKC3A389QRNYaYyY0d85OXavN2YmNrqIoSqTo3IanTg2PoiihKSoqIiEhgezsbE/bihUrGDVqFMOHD+eee+4JONZuv6b0veaaa0hNTWXMmDGetsrKSiZNmsS4cePIysrijjvu8Dw7efIk2dnZxMbGcujQITsfObwYYzrtq0vaCKMoSvjYvHlzpLdgiyeffNJMmDDB5OTkeF6XXXaZ5/nOnTtNVlaW5+va2lozdOhQ880335iqqiozduxYs2nTpkbz2u3X1L4ff/yxWbt2bYM91dXVmePHjxtjjKmurjaTJk0yq1evbjBu0KBBpqSkxP43xvj/GQL5pgW/ezt1kVANtSlK52Dq1KkcOWIVPNm6dSvPP/88V1xxBQCLFi1i48aN5ObmkpSUZGu+vLw8hg8fztChQwG46qqrePPNNxk9enSz+jW175QpUygqKmrQJiIkJiYC1sXdmpqadptR2KkNjzFQVeskLtoRurOiKC3izn9tYvPesladc3RaMndclhWyX25uLgCPP/44H374Id///vcBcDqdvPTSS+Tl5eFw2P89sGfPHgYMqL9qmJ6ezueff97sfk3tGwin00lOTg6FhYUsXLiQyZN9i8W0Dzq14QE4XllLXKIaHkXp6Dz33HO8/fbbLFmyxGNkSkpKKCwsJCcnp1H/tLQ0li9f7ncu4yda4s+7sNuvqX0D4XA4WL9+PaWlpcyaNYuCgoIG50DthU5veMpO1tA7MS7S21CUDo8dzyRcvPbaa7z44ou8+eabDeqOpaSkMGTIENauXdskjyc9PZ3du+tLURYXF5OWltbsfk3tG4ru3btz3nnnsWLFinZpeDp1VhtYHo+iKB2XZcuW8dhjj/H6668THx/f4JnD4eCHP/wh1157LcePH7c958SJE/n666/ZuXMn1dXVLF68mMsvv7xJ/aZNm8aePXuaPGcgSkpKKC0tBawstvfff5+MjAzb49uSTm94yiprIr0FRVHCyPz58ykuLubss88mOzubp59+usHzX/3qV5x11llMnTqVCRMmeF7BfulHR0fzyCOPcOGFF5KZmcmVV15JVla9R3fxxRezd+/egP3q6uooLCykZ8+eTZ4TYO7cuZx55pls27aN9PR0nn76afbt28f555/P2LFjmThxItOnT+fSSy9trW9jq9KpKxfE9RthXn9nJZeM7RfprShKh6SjVC4oKiri0ksvpaCgIHRnGxQUFPDMM89w//33t8p8dhk8eDD5+fn07t3b9hitXBAG1ONRFCUUDoeDY8eONbhA2hLGjBnTpkbHfYG0pqaGqKjI/9rv9MkFx30Mz9cHjjM8NbHd5r8ritL2DBgwoMHB/6lGQkIC69evj/Q2PETe9EWYspP1yQU7D5Uz/S8rWfl1OygpoSiK0kFRw+Pl8RQftWR8SiuqI7UdRVGUDo8anpP1hufwCcvg1Dg7b8KFoihKuOn0hsf7Hs+hE1UA1DjrIrUdRVGUDk+nNzzeobbD5W6PRw2PoihKuFDD45VccOi42+PRUJuiKEq4UMOjHo+iKEqb0qkNj+CbXODyeGrV8CiKUk9LFEibqhYaiFDr7d69m/PPP5/MzEyysrJ48MEHPc+2bdtGdna255WcnMwDDzwQOWXSlqjIneqvrmkjzaBbl5maWqcxxpiz/vSBGXTrMrPo3W1BFfkURbFHZ1cgNab5aqHe2Flv7969Zu3atcYYY8rKysyIESMCqqL26dPHFBUVedqCKZOqAmkr44iyqhOcqKqlW0IMh8s1q01ROiKRUiCF1lELtbNev3796NfPqjuZlJREZmYme/bsabSnDz74gGHDhjFo0CBbnzUcdGrDE+2wftBlJ2uJcURRWWMZHA21KUrHIlIKpMFoilpoU9crKipi3bp1fudcvHgxc+fObdJeW5tObXjcHo9vodDaOs1qU5RW5+1fw/6NrTtn39PgosDnK95EQoE0GE1RC23KeidOnGD27Nk88MADJCcnN3hWXV3N0qVL+dOf/tSkvbY2YTU8IjIDeBBwAE8ZY+7xeR4HPAfkAIeBfzPGFLme3QYsAJzAz40x77janwEuBQ4aY8Z4zdUTeAUYDBQBVxpjjgbbn7fhqfLycqo11KYoHYpIKZDawY5aqN31ampqmD17NvPmzfN4dd68/fbbjB8/nj59+jRrr61F2AyPiDiAR4HpQDHwhYgsNcZs9uq2ADhqjBkuIlcB9wL/JiKjgauALCANeF9ERhpjnMCzwCNYBsubXwMfGGPuEZFfu76+Ndgeo6PqQ21RUn+fR0NtihIGbHomrY1bgXTZsmVBFUgfeugh22c83mqh/fv3Z/Hixbz00ku291RSUkJMTAzdu3f3qIXeemv9r6tp06bx3HPP0b9/f9vrGWNYsGABmZmZ3HzzzX7XffnllyMeZoPwplNPAgqNMTuMMdXAYmCmT5+ZwD9c7/8JTBPLf5wJLDbGVBljdgKFrvkwxqwEjvhZz3uufwDfC7VBb4/HfYcn1hGloTZF6UBEUoEUmq4W2lR1Uvdaq1at4vnnnyc3N9eTNu0dKqyoqOC9997z6wm1NeEMtfUHvAUsigHfky5PH2NMrYgcA3q52tf4jO0fYr0+xph9rrn2iUhqqA1GRwm1WHd5KmucAKQmx2moTVE6EIcPHw7Z57rrruO6665r0rwXX3wxF198sd9n3r/wX375Zb991q1b57d98+bNzJ49m4SEBFvruddKS0vzexbkpkuXLra+F21BOD0efydfvt+VQH3sjG0WInK9iOSLSP6J48cBq1DooRPVJMVFkxgXTa0aHkVRvGhtBdJgtKU6aaSUScPp8RQDA7y+Tgf2BuhTLCLRQDesMJqdsb4cEJF+Lm+nH3DQXydjzBPAEwAZp2WbmLhoT6itd1IcMY4ordWmKEoDTnUF0kBESpk0nCbuC2CEiAwRkVisZIGlPn2WAvNd7+cAua5bsUuBq0QkTkSGACOAvBDrec81H3jTziaTE2IoO1nL4RNV9OoaS7RD9AKpoihKGAmb4THG1AI/Bd4BtgCvGmM2icjvRcR9avc00EtECoGbsTLRMMZsAl4FNgMrgIWujDZE5GVgNTBKRIpFZIFrrnuA6SLyNVYmXcgUGgGS4l0ez4lqeiXGujweNTyKoijhIqz3eIwxy4HlPm2/83pfCVwRYOzdwN1+2v3mAhpjDgPTmrI/EUiOj+F4ZQ2HTlSRM7gHJ6pqPRUMFEVRlNanU1enBiE5IZrSihqOVFTTu6t6PIqiKOGmUxset8ez60gFxkCvRE0uUBRFCTed2/BgJRdUVFt3eHonxhGjyQWKoihhpVMbHhCS4uuPuTS5QFEUJfx0asPjDrW56e0yPLUaalMURQkbndvwAMkJXh5PVyvUpiVzFEXxprnS14HkqANJUQfDrtS20+nk9NNP99R+g8BS25GSvu7UejxQ7/E4ooRuCTEaalOUTshTTz3F3/72twa1ztLS0li6tP7O+7Bhwzy3/J1OJwsXLuS9994jPT2diRMncvnllzdS+4yOjmbRokWMHz+e48ePk5OTw/Tp0xk9enSDufr378+sWbMC7s/uegAPPvggmZmZlJWVedri4uLIzc0lMTGRmpoazjnnHC666CLOOOMM1q9fz+DBg5v1fWsundvjEUhyGZ6eXWOJihINtSlKB2Tq1Kke7yI+Pp7XXnvN82zRokV8+umn5Obmkp+f73l5Gx1fvKWoY2NjPVLUvvTr14/x48cDDeWovbEjRW13veLiYt566y2uvfbaBu1NldoON53a4xHXPR6AXl1jAUsOW0NtitL63Jt3L1uPbG3VOTN6ZnDrpKCyW0D7kL4OJEdtR4ra7nq/+MUv+POf/8xxVwFkb5oitR1uOrXhwSu5ICUpDrD0eDTUpigdj0hKXweSo7YrRW1nvWXLlpGamkpOTg4fffRRo/5NkdoON53a8Ljv8UC9xxPjiMIYcNYZj1Ccoigtx45nEi4iKX0dTI7arhS1nfVWrVrF0qVLWb58OZWVlZSVlfGDH/yAF154oUE/O1Lb4UbPeFz3eHolWh5PtMMyNur1KErHwC19/frrrweVvvYXngqEtxR1dXU1ixcv9qtYGkqOOpAU9bRp0xqcBdlZ709/+hPFxcUUFRWxePFipk6d6jE6JSUllJaWAniktjMyMmx/3tamUxsesJIJbp4+ku9lWwKnsQ7rW6LnPIrSMYik9HUwOepAUtTNlb4ORjCp7UggwaRSOzo5ORPM2rX5Ddr+8VkRdyzdxJf/PZ2ervCboijNY8uWLWRmZkZ6Gy2mqKiISy+9lIKCgrCvVVBQwDPPPNNmKqQAgwcPJj8/n969ezd65u9nKCJrjTETmrtep/Z4/J0FaqhNURRfVPq6denUyQX+iHGH2mrV8CiKYqHS161Lp/Z4/OE+46mt67whSEVRlHCihscHDbUpiqKEFzU8PmioTVEUJbyo4fFBQ22KoijhRQ2PDxpqUxRFCS9qeHxwh9rU8CiKooQHNTw+1BseDbUpiqKEAzU8PsS4Q22aXKAoihIW1PD4EONJLlDDoyiKRXOlr934ylEHkqIOhp31SktLmTNnDhkZGWRmZrJ69eqg60VK+jqshkdEZojINhEpFJFf+3keJyKvuJ5/LiKDvZ7d5mrfJiIXhppTRKaJyJcisl5EPhWR4c3ZsyedWkNtitJpeOqpp5g4cWLQIqH+pK/ffvttNm/ezMsvv8zmzZsDzu+Wo3bjlqLesGED69evZ8WKFaxZsybgeLvr3XTTTcyYMYOtW7eyYcMGz5qB1nNXLggk6RAubBseEUkWkZ7ul43+DuBR4CJgNDBXRHwFwhcAR40xw4G/APe6xo4GrgKygBnAYyLiCDHn48A8Y0w28BLwW7ufzRsNtSlKxyNS0tfgX466qVLUdtYrKytj5cqVLFiwAIDY2Fi6d+/erPXCTchabSJyA/B74CTgdgMMMDTE0ElAoTFmh2uexcBMwNtMzwT+x/X+n8AjYn03ZgKLjTFVwE4RKXTNR5A5DeCW9usGBK8THoDmhNqqa+u48YW1/OqCUYxOSw49QFE6Ifv/+EeqtrSu9HVcZgZ9b789ZL9ISl8HkqNuihS1nfV27NhBSkoKV199NRs2bCAnJ4cHH3yQrl27Nnm9cGPH47kFyDLGDDbGDHG9QhkdgP6Ad1W9Yleb3z7GmFrgGNAryNhgc14LLBeRYuCHQPCgawCaE2o7UFbJB1sP8mlhSXOWVBSlDXBLX7/44ot+pa/dHpH7dfHFFwecy670tbcctS9uKeri4mLy8vKCSi7YWa+2tpYvv/ySG2+8kXXr1tG1a9cGZ0FNWS/c2KlO/Q1Q0Yy5/flxvt+9QH0CtfszlO45fwlcbIz5XET+E7gfyxg1XFDkeuB6gIEDBzaarDmhtipX30Mnqm2PUZTOhh3PJFxESvrajhy1HSlqO+ulp6eTnp7u8WTmzJnjNwnhVJG+vg34TET+JiIPuV82xhUDA7y+Tqdx+MvTR0SisUJkR4KM9dsuIinAOGOM2/d8BTjL36aMMU8YYyYYYyakpKQ0et6cUFtVrROAkuNVtscoitI2RFL6OpAcdSgp6uZIX/ft25cBAwawbds2AD744ANGj7aOwNub9LUdj+dvQC6wEWjKifsXwAgRGQLswUoW+HefPkuB+cBqYA6Qa4wxIrIUeElE7gfSgBFAHpYn5G/Oo0A3ERlpjNkOTAe2NGGvHppzgdTt8ajhUZT2x/z58+nZsydnn302AD/72c88B/BgSV8/+eSTTJ06tUFIKy0tLWCCgbcUtdPp5Jprrmkkff3UU08FzBbbt28f8+fPx+l0UldXx5VXXulJtQ4lfe27nvdaDz/8MPPmzaO6upqhQ4fy97//PeR6EcEYE/QFfBaqT5CxFwPbscJ1v3G1/R643PU+HngNKMQyLEO9xv7GNW4bcFGwOV3ts7CM4wbgI++5Ar1ycnKML3V1dWbQrcvM/e9ua/QsEKu/OWQG3brMXPiXj22PUZTOwObNmyO9hVZh586dJisrq03W2rhxo/nlL3/ZJmu5GTRokCkpKfH7zN/PEMg3zbQLxhhbHs+HrnORfwGeP+mNMUdsGLXlwHKftt95va8Erggw9m7gbjtzutrfAN4ItadQiAjRUdLEUJt6PIrSkfGWvg63YmdbS1+feeaZ7VL62h0eu82rzU469SlLjCOqaaG2GuuM50hFNTXOOk+4TlGUjoFKX7cuIQ2PMWZIW2ykPRHtkCYJwVW7KlkbA0fKq+mTHB9ihKIoSufFzgXSH/lrN8Y81/rbaR/EOqKaFmqrqe9bcrxKDY+iKEoQ7ITaJnq9jwemAV8CHdbwxDiiqKltelYbQMkJPedRFEUJhp1Q28+8vxaRbsDzYdtROyDaIU0SgnPf4wFNMFAURQlFc07BK7Du1XRYYh1R1NTZ93i8z4MOqcejKIoSFDtnPP+ivixNFFZV6FfDualIY4Xamp5OnRDjUI9HURQlBHbOeP7X630t8K0xpjhM+2kXNCfU5ogS+naLV8OjKIoSgqCGx6V/89/GmO+20X7aBTE+obbSimr2llYGlDyorq0j1hFFSmKcGh5FUZQQBD3jMcY4gQpXQkGnIdYn1PbEyh1c/sinfH3AfxHBqto64mKiSEmK0zMeRemAtET62p8cdbikr//yl7+QlZXFmDFjmDt3LpWVlQDs3r2b888/n8zMTLKysnjwwQeB9i19XQlsFJGnm1id+pTFN9R2tKKa2jrD797c5FcXo6qmjrjoKHonxqrHoyinIOGUvvYnRx0O6es9e/bw0EMPkZ+fT0FBAU6nk8WLFwNWkdFFixaxZcsW1qxZw6OPPsrmzZsjJn1t54znLdfLG/spX6cgMY4oyqvrU6QrXO9X7zjM0g17mZndUM+uqtZJXLSDlKQ4yiprqaxxEh9jX9tDUZTwMnXqVI4cscpLbt26leeff54rrrDKRC5atIiNGzeSm5tLUlKSrfm8pagBjxS1W4bAjVuO+tlnnwUsOerY2FiAZktfB1uvtraWkydPEhMTQ0VFhceg9OvXj379+gGQlJREZmYme/bsaTS+rbBjeLobYx70bhCRm8K0n3ZBjEMahNoqqp2M7JNIXLSDu9/awtSMVJLi68Wkqp11xEZboTawUqrTe3Rp830rSnvmk1e3c2j3iVads/eARL5z5ciQ/SIlfR1Mjrq1pa/79+/PLbfcwsCBA0lISOCCCy7gggsuaDRXUVER69ata/fS1/P9tP24lffRrojxKZlTUV1LYlw0d31vDCUnqngkt7BB//pQm9vwqBKporQ3IiF9HUyOurWlr48ePcqbb77Jzp072bt3L+Xl5Q2UTgFOnDjB7NmzeeCBB0hO9p8s1RYE9HhEZC5WZeohLmE2N0nA4XBvLJL4VqeuqHaSGBdN9oDunDGkF2t2NlSEqKq1DI/b49FzHkVpjB3PJFxESvrajhx1a0lfv//++wwZMgS3svL3v/99PvvsM37wgx8AVkhv9uzZzJs3z+PxRYpgHs9nwCJgq+tf9+tXwIzwby1y+FanrqhykuA6s0mKj/bIILiprm0YalPDoyjth0hKXweSow6H9PXAgQNZs2YNFRUVGGP44IMPyMzMBCyPacGCBWRmZnLzzTfb/pzhIqDHY4z5FvhWROYBe12ibYhIApAOFLXJDiOAb3XqippausRahich1sFJH8NTVeuke5dYenWtP+NRFKV9EGnpa39y1Lt372516evJkyczZ84cxo8fT3R0NKeffjrXX389AKtWreL555/ntNNO86SE//GPfwwaTgwndpILXgXO8vraiSVXPdF/91Mf31DbyWonXeKsb1VCjIOT1b6Gxwq1xUZH0b1LjHo8itKOOHw49MnAddddx3XXXdekeS+++OKAv7iXL68XSc7OziY/P7/B8x49erBu3Tq/Yzdv3szs2bNJSEiwtZ73WnfeeSd33nlnoz7nnHOO33OiSGEnuSDaGOM5LXe9jw3fliJPtE9WW3mVky6uUFt8jIPKRh5PHXGu51q9QFE6Ht7S1+GmraWvs7Oz26X0dYmIXG6MWQogIjOBtrviGgGs6tSW4amrM5yscTYItVXWNKzj5i6ZA9A7MU41eRSlg6HS162LHcPzE+BFEXkU6+JoMeBXlbSj4B1qq3Rp7XiH2qqdddQ664h2GZuqWidxMdb7lKQ4NhSXRmDXiqIopwZ2hOC+Ac4QkURAjDH2Uz9OUaIdgrPOUFdnKK9yGZ5Yd6jNMjCVtXUkug2P6x4PWIZHQ22KoiiBCRnUE5E+IvI08Jox5riIjBaRBaHGncrEuAxKTV2dJ5GgS2y9xwM0SDCoclUuACvUVlHtpLyqlt1HKti091hbbl1RFKXdY+c06VngHcB9W2k78Itwbag94D6vqXEayqtrAW+Px/rXnWBgjKG6to64aFdygesuzyUPfcJ3/vwhsx77rIE0tqIoSmfHjuHpbYx5FagDMMbUYqVUd1iiHVYpipraOk+B0ASv5AKoNzxu9VF3qC0rLZkusQ76JMczNSOV6to6jp2sadP9K4qitGfsGJ5yEemFqyK1iJwB2IoficgMEdkmIoUi8ms/z+NE5BXX889FZLDXs9tc7dtE5MJQc4rF3SKyXUS2iMjP7ezRH/5CbV19Q20BDE9mv2Q23Xkhr9xwJjOzLSexTA2PoiiKBztZbTcDS4FhIrIKSAHmhBrkUi99FJiOlQn3hYgsNcZ4i0gsAI4aY4aLyFXAvcC/icho4CogCyvE976IuAs9BZrzx8AAIMMYUyciqTY+m1+Chdp8z3iqfQyP67MD0C3BqgmlHo+iKEo9drLavhSRc4FRgADbjDF2fpNOAgqNMTsARGQxMBPwNjwzgf9xvf8n8IhYv7VnAouNMVXAThEpdM1HkDlvBP7dGOMOCR60sUe/eIfa6pMLLIMT18jjsf51n/F4o4ZHURSlMXay2uKBnwN3AXcCC11toegPeN+4Kna1+e3jOjs6BvQKMjbYnMOwvKV8EXlbREbY2KNf3KG22rr6Mx7frLZGZzwxjb+VangUpWPQEulrf3LUgaSog3HNNdeQmpoasIJ1qD0NHjzYU6ttwoQJQPuWvn4OK+T1MPAIMBp43sY4f3J6vsWCAvVpajtAHFBpjJkAPAk843dTIte7jFN+SUmJ3427DU91raHCFWprnFxQ5+pj/esOz3njMTwVangUpT0TLunrQHLUgaSog/HjH/+YFStWBHxuZ08ffvgh69ev99SOa8/S16OMMeO8vv5QRDbYGFeMdebiJh3YG6BPsYhEA92AIyHGBmovBpa43r8B/N3fpowxTwBPAEyYMMFv1bwYd6jN6e3x+Jzx2PB4kj0eT62/ZRRFaSMiJX0N/uWomyNFPWXKFIqKilplT5HGjuFZJyJnGGPWAIjIZGCVjXFfACNEZAiwBytZ4N99+izFUjhdjZWwkGuMMS7huZdE5H6s5IIRQB6WxxNozv8DpmJ5Oudi3TdqFt6htvLqWmIdUZ423+QCtzaPvzOeGEcUXWMdGmpTFODDZ5/g4Lc7WnXO1EFDOf/H14fsFynpazty1K0lRR1qTyLCBRdcgIhwww03eCQTIoEdwzMZ+JGI7HJ9PRDYIiIbAWOMGetvkDGmVkR+inX51AE8Y4zZJCK/B/JdRUefBp53JQ8cwTIkuPq9ipU0UAssNMY4AfzN6VryHqyacr8ETgDXNuk74YV3qO1ktdMTXgOIj7We+Xo8sdH+o5bdEmLU8ChKO8Atfb1kyRK/0te+pKWlNZAc8Mau9LW3HHX37t254ooreOGFFzyqoK0pRR1qT6tWrSItLY2DBw8yffp0MjIymDJlSovWbC5awns+AAAgAElEQVR2DE+z1UaNMcuB5T5tv/N6XwlcEWDs3cDdduZ0tZcClzR3r974htq6ehmeWEcUIvXJBf7Sqb1JVsOjKAC2PJNwESnp62By1K0tRR1qT+73qampzJo1i7y8vIgZHlt6PMB+lyLpEKz05WPGmG9dbR2OhllttQ08HhFpIAZXf4HU/3+03RJi9AKpokSQSEpfB5KjDiVF7St93dI9lZeXez5feXk57777btDsuHBjx/AsAZwiMhwrNDYEeCmsu4owDbPanHSNa+gYJsQ4PHIJ7ns8GmpTlPbJ/PnzKS4u5uyzzyY7O5unn366wfNf/epXnHXWWUydOjVoVps33lLUmZmZXHnllY2kr/fu3dtAjvq0006jrq6O66+/3iNFnZubS3Z2NtnZ2Z6wnj/pa4C5c+dy5plnsm3bNtLT0z2fw71WsD0dOHCAc845h3HjxjFp0iQuueQSZsxodjCrxdgJtdW5zmu+DzxgjHlYRPxrtnYQGoTaqpyehAI38TEOTlY3TKcOFGpTw6MokSXS0tf+5KiDSVEHkr5++eWXQ64VaE9Dhw5lwwY7ychtgx2Pp0ZE5mKJvy1ztcUE6X/K0yDUVlPrSaV2Y6mQ+q/V5osaHkU59VHp69bFjsdzNZYK6d3GmJ2uVOYXwrutyBLjMiI1rlBbFz+htkYlc2ICn/GcrHFa8tgBjJOiKO0blb5uXezUatuMVTLH/fVOrNTlDktMlCvUVmeF2ro0CrVFNSoS6q9yAUC3LpZzWFZZQ+/EuHBtWVEU5ZRB/wT3g0cWodbKavMNtcU38HjqEKk/F/JF67UpiqI0RA2PHzyhNqfhZI3/UJv3GU9cdJTfy2PgXTZHDY+iKAqEMDwi4hCR+9pqM+2FaFeoraLaSY3TNAq1NUguqHEGDLOBejyKoii+BDU8rjI1ORLoz/kOijvU5jYWwZILqp11ARMLoN7w6CVSRVEUC1tFQoE3ReQ1oNzdaIx5PWy7ijCOKMERJfWGx98Zj6dIaF3AVGqA5Hj1eBRFUbyxY3h6AoexKj+7MUCHNTxghduCGR63Ho/7jCcQqsmjKIrSkJDJBcaYq/28rmmLzUWSWEeUJzzmVh91kxDjoNpZh7POUFVbR2yAOm1gldJJiFFpBEU5lfGnQGpHERQCq4L6UwQNRqj1KisrmTRpEuPGjSMrK4s77rjD8yyQ4mm7VSAVkZEi8oGIFLi+Hisivw3/1iJLTHRUQI8nwSWNUFnjpKrWGdTjAa1eoCjtnaYqkEJoRVAIrQrqqwgajFDrxcXFkZuby4YNG1i/fj0rVqxgzZo1AAEVT9uzAumTwH8CfwMwxnwlIi8BfwjnxiKNd6gtwdfweKmQhgq1gRoeRYk0ra1ACqEVQaF1VUFDrSciJCYmAlBTU0NNTY3nmkdzFE/DiZ17PF2MMXk+bR1eyznGUe/xdPUJtcV7qZBW2SiFo4ZHUSJLbm4u69ev54YbbuDyyy9vpED69NNPN8no2MWfKqhb7sCtCJqTk8MTTzzRKus5nU6ys7NJTU1l+vTpflVNW0vxtCXY8XgOicgwrIQCRGQOsC+su2oHxEZHeVKm/RUJBSvUVl1bF1CLx01yQgzFRyvCs1FFOUUo/dc3VO8tD92xCcSmdaX7ZcNs9W1NBVK7BFMFDYciqMPhYP369ZSWljJr1iwKCgoanAm1puJpS7BjeBYCTwAZIrIH2AnMC+uu2gHuS6TgJ6st2jvU5iQuJrTHs3mvejyKEilaW4HULsFUQcOpCNq9e3fOO+88VqxY4TE8ra142hLsFAndAXxXRLoCUcYY+zJ9pzAxXtUIGmW1eTyeOuseT5DKBaChNkUBbHsmrY1bgXTZsmVBFUgfeuihVg+3eauC9u/fn8WLF/PSSy9RXl5OXV0dSUlJHkXQ3/3ud55x06ZN47nnnqN///621yopKSEmJobu3btz8uRJ3n//fW699VaAkIqnbY2drLZvRORF4IfAgFD9Owruem0iVjVqb+K9kgusygWhDU95tZMaZ114NqsoSkDCoUAKgRVBIbQqaDBF0OYqkO7bt4/zzz+fsWPHMnHiRKZPn86ll14KEFTxNBLYCbWNBiYD3wH+V0QygA3GmFlh3VmEcUsjdIlxNCoAmuCdXFDjDHnG0y3B+jaXnayhV2Ictc46RKzqCIqihJdwKZAGUgSF0KqgwRRBm6tAmpaWxrp1/sWhgymeRgI7WW1OoMb1bx1wADgYzk21B9yhtoTYxrbZO7nAVjp1l4Zlc2b/dTW//b+NrbldRVHCiCqQti52PJ4yYCNwP/CkMSb0nw8dAHeozTexABrf47GTTg2W4dlTepINu0vZsq+M/7owgx5dY1t554qitDaqQNq62DFxc4GVwH8Ai0XkThGZFt5tRR5PqM2P4XGf+RyvtDwYOxdIwTI8n2wvASzl0iVfFrfafhVFUU4V7NRqe9MY85/ADcBy4MfAsjDvK+K4Q23+DY/VVlrhNjyhzni8DM/Xh+iTHMfpA7vzct6udhV3VRRFaQvsZLUtEZFvgAeBROBHQA87k4vIDBHZJiKFIvJrP8/jROQV1/PPRWSw17PbXO3bROTCJsz5sIicsLO/YLhDbV3jGkcjLcVRKHWd2YQKtblVSI+WV/Np4SGmjEjh3ycN5JuScvJ2HmnpVhVFUU4p7ITa7gFGGmMuNMbcZYz52BhTGWqQiDiAR4GLsDLj5oqIb2GgBcBRY8xw4C/Ava6xo4GrgCxgBvCYSw016JwiMgHobuMzhcQdakvwI/ImIg0qTtsNtX1aeJhjJ2v4zsgULh2bRlJ8NC/n7WqN7SqKopwy2DE864GFIvJP1+tnIhITchRMAgqNMTuMMdXAYmCmT5+ZwD9c7/8JTHOpnc4EFhtjqowxO4FC13wB53QZpfuA/7Kxt5AEC7WBZZDcGjuh7vHERTuIj4li5fYSROCc4b1JiHXw/dP7s7xgP0fLq1tjy4qiKKcEdgzP40AO8JjrNd7VFor+gHcaSLGrzW8fY0wtcAzoFWRssDl/Ciw1xgStIyci14tIvojkl5SUBOwX7XAlF/gJtYF1zlN60jIYoc54wPJ6qp11nNa/Gz1dmWxzJw+kuraOZV/tDTleURSlo2DH8Ew0xsw3xuS6XlcDE22M83c70vckPVCfJrWLSBpwBfBwqE0ZY54wxkwwxkxISUkJ2M/j8fgJtYGV2eYOtcWGKJkD9eG274zo7WnL6JtMSlIc63cfCzleURSlo2DrAqmrOjUAIjIU6zJpKIppWGInHfD9097TR0SigW7AkSBjA7WfDgwHCkWkCOgiIoU29hiQ2CD3eMC6RGo31Abehqehscvsl8yWfWUt2aqiKMophR3D85/AhyLykYh8DOQCv7Ix7gtghIgMEZFYrGSBpT59lgLzXe/nALnGyi9eClzlynobAowA8gLNaYx5yxjT1xgz2BgzGKhwJSw0G3d16kChtoQYB2WVliyR3VBb11gH4wc2TAjM7JdE4cETVNdqHTdFaa80V/o6kBx1ICnqYNhZr7S0lDlz5pCRkUFmZiarV68OOr7dSl8bYz7A+sX/c9drlDHmQxvjarHOXd4BtgCvGmM2icjvRcRdfe9poJfLO7kZ+LVr7CbgVWAzsAJYaIxxBpqzKR/YLqGSC+K9QnCh0qkBrjl7CHd9b0yjvqP7JVPtrOObkhZngCuK0kzCJX0dSI46kBR1MOysd9NNNzFjxgy2bt3Khg0byMzMDDq+3Upfi0g8VtWCc7DOWT4Rkb/aSak2xizHunTq3fY7r/eVWGcz/sbeDdxtZ04/fRJD7S0U9aG2wB6Pm1Dp1ABnDe/tt310P0uMacu+MjL7RU6YSVE6MpGSvg4kR90cKepQ65WVlbFy5UqeffZZAGJjY4mNrS/JZWe/bYWdWm3PAcepP7ifCzxPAIPRUYgOUjIH6guFgj3DE4ghvbsSGx2l5zxKh+ftt99m//79rTpn3759ueiii0L2y83NBeDxxx/nww8/bCR9nZeXFxYhOPcaOTk5FBYWsnDhwkaS060lRb1jxw5SUlK4+uqr2bBhAzk5OTz44IN07dq1RfOGAzu/MUcZYxYYYz50va4HRoZ7Y5Gmvjp1gFBbdNNCbYGIdkSR0TeJzWp4FCWsuKWvX3zxRb/S126dGvfLV8qgubjlqIuLi8nLy6OgoMDzrDWlqGtra/nyyy+58cYbWbduHV27duWee+5p6fbDgh2PZ52InGGMWQMgIpOBVeHdVuTxlMwJFGpr4PG07C+lzL7JvLflAMaYRto/itJRsOOZhItISV974ytH3dpS1Onp6aSnp3s8pzlz5rRbw2PnT/XJwGciUuRKVV4NnCsiG0Xkq7DuLoIEq04NDZML7KRTB2N0WjJHyqs5eLyqRfMoitIYt/T166+/HlT6+vjx462+dklJCaWlpQAeOeqMjIyQUtTTpk1jz549TVqrb9++DBgwgG3btgHwwQcfBD0ziiR2fmPOAIYA57peQ4CLgUuBy8K3tcjiTqNOjvdfHaipyQXBcCcVbN6r4TZFaW0iKX0dSI46mBR1c6WvAR5++GHmzZvH2LFjWb9+Pbfffrut/bY1IUNtxphv22Ij7Y0LRvfhiR/mMLBXF7/PE2LrjY2dygXByOhnZdJs3lfG+RmpLZpLUZSGRFL6OpAcdTAp6uZKXwNkZ2eTn5/f5P22NW2ndXqKER/j4IKsvgGfuz2e2OioFp/LJMfHMKBngiYYKEo7RaWvWxc7yQWKH+JchqelYTY3mX3rS+dU1To5WFbFgJ7+vS1FUdoWlb5uXdTjaSYJrWx4Rqcls/NQOU99soNz//wRUxd9xP5jIe/oNqDGWceqwrYre6EoitIc1PA0k3rD0zopmJn9kjEG/vDWFnp2jaXGaVi9o2lG5K2v9jHvqc8p2KPVrhVFab+o4Wkm7ns8reXxTBmRwg1ThrL4+jNY9rNz6JYQw+pvQh+KeuOu9/a5ymkritKOUcPTTOK9kgtag4RYB7ddnMkZQ3sRFSVMGtKTNTuaZkCKDlcAsPZbNTyKorRf1PA0k9Y+4/HlzKG92HWkgj2lJ22P2XW4HID8oqMBUzUVRVEijRqeZhLvqlbQWmc8vpwxtBcAa5oQbis6XEGXWAcHj1ex+4h9g6UoitKWqOFpJp4znhaWywlERt8kenSJYfUOe4antKKaYydruGiMVWo9X8NtiqK0U9TwNBPPBdIWVi0IRFSUMHlIL9bYNDzu850LsvqQFBdN/rdHw7IvRVGUlqKGp5m4kwvC5fEAnDG0J8VHT7L7SEXIvt+6zneG9u5K9sDurC1Sw6MorUVzpa8hsBy13fFNWe/BBx9kzJgxZGVl8cADDzR4NnjwYE477TSys7OZMGEC0I6lrxX/xEVHIRK+Mx6AM4dZqqV2vJ5vD1cgAgN6dmHi4J5sP3icYydrwrY3RelIhEv6GgLLUdsdb3e9goICnnzySfLy8tiwYQPLli3j66+/btDnww8/ZP369Z56bpGSvlbD00xEhIQYR9hCbQAjUhPp2TW2UVq1MYaFL37JioJ6Nceiw+X0TY4nPsbBhEE9MAa+3KVej6KAJX3trgIdHx/Pa6+95nm2aNEiPv30U3Jzc8nPz/e8li5dGnTOKVOmNKog7YtbjnrBggWAJUfdvXt32+Obst6WLVs444wz6NKlC9HR0Zx77rm88cYbtudvS7RWWwsY1KtrwOrVrUFUlHDG0J6NPJ79ZZW8tXEfJ6pqmTHGKmT67eEKBrn2kj2wO44oIb/oCOeP0mrXSvtg+/a7OH5iS6vOmZSYyciR/x2yX6Skr9tSjnrMmDH85je/4fDhwyQkJLB8+XJPSA2sP5YvuOACRIQbbriB66+/vtX3YBc1PC3gXz89m6gwK4aePqAHyzfu59CJKnonxgGwaY9VTDRv5xGqap3ERTv49nA5383sA0CX2GhG90smX895FMWDW/p6yZIlfqWvfUlLS2sgOdAc3HLUDz/8MJMnT+amm27innvu4a677mrRvP7IzMzk1ltvZfr06SQmJjJu3Diio+t/xa9atYq0tDQOHjzI9OnTycjIYMqUKa2+Dzuo4WkB0WEMs7kZm94NgK+KS5maYRmWgr1WLbaTNU7W7SplTP9uHDpR3cD7yhnUg8Vf7KLWWdcm+1SUUNjxTMJFpKSv21qOesGCBZ6w3u233056errnmfscJzU1lVmzZpGXlxcxw6O/kdo5Y/p3I0rgq+L6wp8Fe8ro1y2eKIHPCg95MtoG96p338cN6EZlTR3flJS3+Z4VpT0RSenr5spRN0f6GuDgwYMA7Nq1i9dff525c+cCUF5e7vl85eXlvPvuu7az6cKBGp52Tte4aIanJjYwPJv2HmPSkJ6MG9CdTwsP8a3rDs8gL49nTJrlKWmlaqWzE0npawgsRx1ofEukr2fPns3o0aO57LLLePTRR+nRowcABw4c4JxzzmHcuHFMmjSJSy65hBkzZjT3W9piNNR2CnBa/+58vP0gxhgOl1ez71glY9K6MaBHFx7/+BuPcRnk5fEMTUkkIcbBxj3HmJ2THmhqRenwRFL6GgLLUQca3xLp608++cRvn6FDh7Jhw4aA+21rwurxiMgMEdkmIoUi8ms/z+NE5BXX889FZLDXs9tc7dtE5MJQc4rIi672AhF5RkRiCDd5T8LuL8K+zLgB1hnO3mOVbNprJRZk9U/m7OG9cdYZlnxZTO/EWBLj6v+OcEQJo9OS2bRXPR5FaSkqfd26hG0lEXEAjwIXAaOBuSLiG9xcABw1xgwH/gLc6xo7GrgKyAJmAI+JiCPEnC8CGcBpQAJwbbg+GwA1J+HtW2HNo2FdBmBsupX3/9XuUo8hyUrrxvhB3YmPieJAWVUDb8fNmLRkNu0to65OK1UrSktwS19HQiY6nLgvkO7Zs6dJd4paSjhN3CSg0BizwxhTDSwGZvr0mQn8w/X+n8A0ERFX+2JjTJUxZidQ6Jov4JzGmOXGBZAHhDe+dHAzGCccbN17Cf7I7JdEjEP4as8xNu0pY2DPLnRLiCEu2sHEwdZ/LIN6Nr5PNKZ/Nyqqnew4pAkGiqK0H8JpePoDu72+Lna1+e1jjKkFjgG9gowNOacrxPZDwG9tCRG5XkTyRSS/pKSkiR/Ji31fWf8e+hpqq5o/jw3ioh2M6pvEV8WlFOw9RlZasufZOcOtsjp+PZ7+VoKBv3DbkfJqxt/1HisK9oVp14qiKP4Jp+Hxd7PSN+YTqE9T2715DFhpjPF7ymaMecIYM8EYMyElJcVfF3vscx3UGadlfMLM2PTurNtVyreHKzwGBeDcUdZnGNU3qdGY4amJxEZH+c1se2/zfo6UV/OvDWp4FEVpW8JpeIqBAV5fpwN7A/URkWigG3AkyNigc4rIHUAKcHOrfIJg7P8KEq0LnW0RbhuXboXNgAYeT0bfZN6/eQoXjO7TaEyMI4rMfsls9GN43tl0AIBPvi6h1lkXpl0riqI0JpyG5wtghIgMEZFYrGQB36p7S4H5rvdzgFzXGc1S4CpX1tsQYATWuU3AOUXkWuBCYK4xJry/SZ21cGATjP4eREXDwU1hXQ6slGo3WWndGjwbnppEVJT/0j1j0pLZtKdhgsHxyho+/foQg3p1oayylg3FpeHZtKIoih/CZnhcZzY/Bd4BtgCvGmM2icjvRcR9M+tpoJeIFGJ5Kb92jd0EvApsxjqrWWiMcQaa0zXXX4E+wGoRWS8ivwvXZ+PQdqithP450GtEm3g8I/skEh8TRd/keFKS4myPO61/N45X1bL7aL2mz0fbSqh21vHbS0YTJfDx9rbT4VAURQnrBVJjzHJguU/b77zeVwJXBBh7N3C3nTld7W13GXa/K7Gg3zjoMxqKA9zlqauDmnKIa3z+0lSiHVGcMzyF7l2adj3JfR60cc8xTwLCik376Z0Yy9SMVLIHdOfj7SXcPH1ki/eoKIpiBy2Z0xz2fQXRCdB7BKRmQukuqPJT5yn/abhnECz9ORxret0lX578UQ73zRnbpDEj+iQS4xAKXBWtK2ucfLT1INNH98ERJZw7MpWviks5Ul7dor2drHby8fYWZAkqitJpUMPTHPZtgD5ZEOWA1Cyr7eDWxv2+XQXRcbD+JXjodFjdssumIoI0UYbBnYr9zqb9bNt/nFWFhyivdnJBlqXjc+6oFIyxkgyC8f7mA6zfHfgs6KlPdjD/mTytDad0SFoifR1IjtqfFHUgdu/ezfnnn09mZiZZWVk8+OCDjfpUVlYyadIkxo0bR1ZWFnfccUfI8Sp9fapgDOzfaIXZwPJ4wLpQ6sv+Ahg2FX62FgadCe/fCZVt/4v55ukjOVpRzSUPfcJdyzaTFBfNWcN6AdYZUI8uMax0nfPsPFTOO5v2Y+V4WKz+5jDXP5/PD5/+3FMJ25e3Nlpp2f/6yjdxUVHaP+GSvg4lR+0rRR2I6OhoFi1axJYtW1izZg2PPvoomzc3/J0TFxdHbm4uGzZsYP369axYsYI1a9YEHa/S16cKR4ug6hj0c4W8ug+CmK6NDU91ORwuhD5joMcgOP+34KyCLcvafMtTM/qQ+6vzmJOTTtHhCqaP7kNctKU94ogSzhmRwsfbD3LzK+uZtugjbnh+LXe/tcUqSnqiil+8so6BPbsQJcLCl76kssbZYP6dh8rZuv84MQ7hra/2NTBaitIeiJT0dWvJUffr14/x48cDkJSURGZmZiPZBBEhMTERgJqaGmpqajwREjvj2xKtTm2HyjI4uhP6jq2/ONrXZXiioiA1o7HhObgFMNDX5YanT7CM1MbX4PR5bbZ1Nz27xnLP7LHccO4weiXGNnh27sgU/rVhL8sL9rHgnCFU1tTx1Kc7OVnjZN+xSo5W1PDMf0xk/7FKFvwjnz+8tZk/fO80z/i3XdUPfj51BIve28763aWcPrBHm34+pf3z318XU3DiZKvOOSYxgbtGhK6OFSnp62By1M2Voi4qKmLdunUecTlvnE4nOTk5FBYWsnDhQr99go1vK9TweFNdDivvgzMWQqJXVYMP7oQvnrKMTZdeIA5I9ap3mpoJ299pONf+jda/fV2/oEXgtCvg0/vhxEFITA3vZwnAkN6NS+vMzLbc7PNGpdA7MQ5jDF3jovnrx98A8PuZWWSldSMrrRs3TBnK31buYPKQXlw2zhr39sb9jBvQnflnD+bh3EKWfbVPDY/S7oiE9HUwOermSFGfOHGC2bNn88ADD5CcnNzoucPhYP369ZSWljJr1iwKCgoanEGFGt9WqOHx5rNH4NO/WKGzc//TajMGtr1thcxqq2DHh9DnNIjxUjJMzYJ1L8CJknqDdaAA4pItL8fNaVfAJ/8Lm96AyTe03ecKQYwjijlemj0iwq0zRpGSFMfBskp+eEb9Z7jlwlHkFR3hN29sZMLgHtQ6DRv3HOO2izJIjo9hysgU3vpqH7+5ODPgpValc2LHMwkXkZK+hsBy1E2Voq6pqWH27NnMmzfP47EFonv37px33nmsWLHCY3iaMj7c6BmPm+MHYJUrU2Sr1znMgQIo2wOTfwL/sQZ+9CbM+mvDsf4SDPYXWJlv3lloqRmW0dr4Go0wBoo+teQWAnFkJ5S1zeG9iLDgnCHcdnFmg0y6GEcUf7kym9o6wy2vbWC5K6ngojH9ALhsXD/2l1WydtfRNtmnooQiktLX4F+OOpQUta/0tTGGBQsWkJmZyc03+68IVlJSQmmplXl68uRJ3n//fTIyMmyPb0vU8Lj56E/W4f/4H8G+9VDqKoLtDqGNuMA6zxl6Xv25jZu+pwFipU+DdXH0wCbLS/LltDnWhdMjOxu2f/EUPHsJPDAWPnvYCvt5c/RbeOJceGYG1FTWtxtjidHVNTzwDyeDe3flvy8dzarCwzzw/tdkpSUz0CW7PS2zD3HRUSzboNltSvsg0tLX/uSog0lR+5O+XrVqFc8//zy5ubmeJAl3GNC91r59+zj//PMZO3YsEydOZPr06Vx66aUhx0cEY0ynfeXk5FgCPge3GvM/PYx56xZjDhUac0eyMasft549+V1j/jrFhOQflxvzlzHGOJ3GHP7GmiP/7437Hd1lPfvgrvq20mJj7u5vzFMXWPPckWzMfSOMKfrMel5TZcwTU425q4/17OM/149d81er7bNHG6+zdbkJF3V1dWbBs3lm0K3LzCO5Xzd49pPn8032ne+Yo+VVIec5Wl5lXsnbZb7YedhU1tSGa7tKhNi8eXOkt9Aq7Ny502RlZbXJWhs3bjS//OUv22QtN4MGDTIlJSV+n/n7GQL5pgW/ezu3x1NeAp/cD28uhNiucO6t0GsYpGRa4bbyw5Z3MvLC0HNlz7MqGHy7ygqzgRVW86X7AMi8DFb+L3z5nOWxLL8F6mqtEN6P3oRr3oXYRHjucuvy6Qd3wp58mPW4NfaT++FYMexdD+/+1pp37d+tudy8dTO8fBWUbGv598kPIsI9s8fybxMGcMWEhrH7n00dwfHKWn6/zM/dJi8KD57ge4+u4r+WfMWcv67mtP95l5+/vE4VU5V2h0pfty6dO7ngWLH1Sz0qGi6+D7paompkXIL59H7K1j9HNwyMsGF4Mi6F2CTY8DJ0SweJqj/78eX7T8IrP4ClP4Odn8C25TD9Lug5xHo+cDJc+z68Nh/+70arbeK1kDUL0sbD1+/B8v+Cki3QpbeVqPD+HbBrjXVR9eAW+Ppda9ynD1gGKwz0TozjXj8lfEanJfOTc4fxyIeFXDYujfNHpXKy2smLn39LVW0dw1MTqXHWcdvrG4mLjuLZqydSVVvHO5v28/qXe5iTk86UkS3QSlKUVsYtfd3RcF8gbWs6t+FJ7AO/3QmO2AZJAOv6jeR/+6aw5eun+WdyKkPTTg89V2wXGDMLNi6x7uz0Gm61+SMmAa56CV6dDxtftdK0z/iPhn269IQfvG55NEd2wAWueqk9BsHZv4CP77GM21dadu8AACAASURBVI/fsqoorPxfWPusZXhWP2LVkht9uTX/+bdB94HN+x41k59NG86KTfv5zesb+d1lWfxx+RZ2Halo0CejbxJPzZ9Aeg/r+3TeqBQ+3lbCc6u/VcOjKB2Yzh1qi0mwaqm5jE6Ns4ZbPr6FH33xB/bFxmKA1/sNs5IK7JA9z6pGvfNj/4kF3kTHwZXPwXfvhCueBYefvwEcMXDRvTDvtYbp22ffBIPOhgv/CIPOssKEY6+Azf9nhda+etW6pDrtd4BYyQotZfs7sHed7e5x0Q7unT2WfWWV/OSFtTiihJeum8zG/7mA/1t4Nn/9QQ5LbjzLY3TcY66aNIDcrQcoPloRZHb/GGN4d9N+lqwtbvJYRVHajs5teKITGny5et9q3il6h6uzrmZZ6gWcV3GSpXVHqXHW2JtvwGToOdR675v55nf9WDjnF9a5UlOI7QJXL4czbqxvy/mxpRH04hXgrLE8qG7pMPbfrLOk4wdg61vwj8tg9WNNW293nnVetHhew4y6EOQM6sFdM8fwnxeO4u2bvsNZw3qTFB9D9oDuzBjTl65xjY3t3EmWZ/bS57uatMXtB44z76nPuf75tfzqtQ0seneblu5pJ+jP4dQlXD+7zh1qi24oqJa3L4+YqBj+I/s/iD+6i1ll3/B+9Q4+Kv6I6YOmh55PBMb9O3z4B/+JBUHYVbaL+764j5vG38TwHsObNBawwm39sq1U8MzL6o3ZOb+A9S/CIxOtGnMxXaxzpZRRMHya1afOaSUq9B4B8T63mavL4Y0brMuwZXsg/xk40ycsGIQfeF0+Zf9Ga52TR6z7SuOugh6DG/RP79GFqRl9eOWL3dz03RFU19bx91VFrNt1lF1HKth3rJIx/btxyWn9OHNYL9Z+e5T3Nh/g4+0lJMZF8/uZWWzeW8bDuYVU1ji5/eJMDp2oZteRCkb1TSLRj7FTwkd8fDyHDx+mV69eTa6srkQWYwyHDx9udPepNejc/xf6/I+Qtz+PcSnjiI+Oh5SRnP1vr5O65AKWfL0kqOE5dPIQ1717HZcPu5wfT1yAOKthSPDSF95UO6u55eNb2HJkC98e/5bFlyymS0yA86FgTLwWlv4Uzrqpvq33CJhwjZWdd9bPYOQMeOZCWHIt3LDSknZYcq2VjRcVY50RjbzISmRI7gfv/rd15+jHy+DjP8Mni6y7TnGJDdfeuRI2L7UyAIdNteb1Ztfn8PeLwHjdN/rsESupY9xVDX4WPzpzEO9vOcBtr29k5fYSDp2oJrNfMiNSkzhrWG/W7DjMHUvr5cbTeySw4Jwh/OTcYfTsGktdnSEuOoonP9nJS5/vorzaWrNHlxh+cu4wfnTmYBJiw3NLXWlIeno6xcXFlJSoVtOpSHx8vKfSQmsindkNnjBhgnGXIz9WdYzvLP4ON467kRuz60NYD697mCe/epJ357xL3659/c5z3xf38dzm5wC4cuSV3Db5NqKj7Nv0e/Pu5YUtL3D1mKt5tuBZLht2GXef00h8NTTGWNWz3dlxgThUCE+cZ6V2H99vlQKa+hs4cQC2v2tly0kUpE+C3WvgzJ/ChXdDcT48NQ2m/hamuEoK1Zy05B4+f9waY+ogub9lBM/6mXVOVX4Y/vYd6/28JZDUByqOwBs/gV2fWUZuxr1WO1BXZ5i66COKDlcwaXBPfnNJJuMGdG/wEb7ed5SiDSsZMO5cRvXt1uivaWMMz35WxDclJxiWkkif5HgWf7GbldtLSEmK46IxfZkyIoWJg3uCQI2zjvgYRwOPyBhD0eEK+ndPIDa6c0elFcUbEVlrjAkuIhSEzu3xeJF/IB+DYVK/SQ3aZw2fxRNfPcEbhW9w47gbG407dPIQr257lUuHXkpql1SeKXiGfeX7uH3y7aQnhf5LYWXxSl7Y8gJzM+Zyc87NxDvieXzD40zqO4mZw2c27UOIhDY6AL2Hw/ceg1d/aCVBXPGs5RkBTP89HPraSlDY+Br0G4fzvNtZ+vUb9E/sz6RRl8Cqh6HbACt0tm25lXU3+Sdw3m2w4yP48h9WmvqWpTDrCXjndqrKS7hl/EVMPbqRWb2HW3LgP14Gqx6AD/8EhR9Y96gm30CUI4bH5uVw4Hgl541MaRyiqalkxEcLGbHtLUj4DfT7Lz/fCuHqsxt+Ly4+rR95O4/wxMod/HNtMc+t/rbRt29UnyQmDO5BZU0dH28voeR4FePSu/HCtZNJim+a7LiiKP5Rj8fl8dyTdw9Lti9h1dxVxDoaygZc++617Dy2k9cvf51ucd0aPLvvi/t4YcsLLP3eUgYlD+LVba/yp7w/UWfqmDZwGvMy5zE+dXyjX56HTh5i8dbFvLjlRfon9ufFS14kzhGHs87Jde9d9//tnXm4HUd14H+nuvvu7963a3uSJdnajG1Z2NjBYBsvgO0QyEAGjE3CBL4EAvnI8EHGQJgsJEwyM86MIUOAhGRCDIHEC2AHL3iceMU23hdZtiTb2pcnvf29u3R31Zk/+sp+1mIJZJ4sUb/v6+/e7q7qrtPV3afrVNU5PLjjQfqKfcytzKW/1E81V6Uj18Hi2mLOHjib3mIvQ40hrll7DbdtvI2VfSt5z5L3sKJnBfdvv5/vr/s+G8Y3cNqs0zhzzpmc3Hsy3YVuRISR5gi3b7qdBzfczruWX8pZC8494DVaPbSaL9z3BZ4eehojhiuW/waX3fSngEKQz/qWzv+DzJXQdJ7+Adz4e1ngO3Vcdfq7+buhh8iZHP/yK//C8Z3TBlTsXg+3fAbW35YNzjj9Q7Dy/S/Nq5pOcxy+exlsuDsbhr7zKfjgjbDwzS+lcQ4GV2d9WflKZl6c7g1clZZ1PLxxhKe2jhEYQy4QhqcSHto4zKObRgmMcPaSXpb0d/BX/7aOVQs6+eaHzqCU899qHs/htni84mkrnnff8G56Cj387dv+dp90jw0+xm/e+pus6l/F1y/8OlGQffnubuzm4usu5m0L3/Yy09jOqZ1855nvcM3aaxiPx5nfMZ93LH4HJ3SewHNjz/Hs8LPcteUuUpdy7vxzueINV7ysdTTaHOWatdeweWIz2ya3MdgYZDKeZDwep2VbACzrWsYLYy8Qu5hT+k5h7fBamrZJR9TBRDJBNVdladdSntz95It58kGe/lI/2ya3YdVSDIvENuYP3/iHvHtJ5q02tjGP73qcx3c9zmODj3H31rvpLnTzydM+yW0bbuOOLXdw+cBb+fTyXyecfVJmPtuL4eYwQ40hlgQVuOnTrC6UuHz8Qc5fcD4P7XiIuZW5XH3J1URmr7zP3pKFjdj8QNbfNHdVpjj29HfZOGuNjW2GX/0aLLsoMxm2JuF37oX6EDzw9Uzp1aeH8RUYeEM2UXh0Y2ZeLPVkpsa+5ZnZcFpL0TmFtInZ/QzsXM2TG3fyrQd3MDCrhxPOvIRx00mcOgJjiAKhryPPm07oJQpebo5TVaZiy8hUzJxagXDX6syL+VmfgNo8AFqp5YdPbKeZON57+gBh4E16ntc+XvEcBnsUz3BzmHP/+Vw+seoT/NYpv7XftDc+dyOfu+dzvL/rrXz6tE+R9ndx5UNXct26615s7exNPalz28bbuPH5G/nJ9p+gKIIwrzKPN817Ex9Y8QEW1hYecnlVlbUja7lzy53ct+0+ju88nstWXMbi2mLG43Fufv5mHtv1GOcOnMt5C84jH+Rp2RaPDz7OutF1bJ/czvap7SyoLuDtC9/OQGWAT9/5ae7ddi/vXfpexuIx7t5yN/U0m0OzsLqQcwbO4SMrP0I1V8U6y5UPXcm31nyL/mI/Fy26iEsWX8LyruUEJsA6yzVrr+HLj3yZiWSCdx7/Tn731N/lY7d/jPF4nO+/6/vct+0+PnXnp/j4qR/noys/+jL5Rpuj1NM6c+tj8MjVsPNJiOtokjlMlSCfzVl68ydhSXuwx/Yn4BsXQrELJndkrbAV74ATLswGeDRG4JmbstZUkMsm0nbMhqndmeeKLW0Hq+d9Dk69LFNaT/xz1p+l+zpeTTTgdvd6brWnM0928zqzgS6ZZF24hMoJZxEufCM/3hHw8MZhNg3XaSaOCnWuKHyPy7kZg2OqspAHzvkmT46Xufr+jeyezD4MThmo8Zf/cSVLZnWAs6SbHiDsPWG/sZtUlXvXD7FxeIp3rpzrzYCeGcUrnsNgj+K5dcOtfPrOT3P1xVdzav+BfTH907c/x7L/8T1sFPD7Hw4ZLlneveTd/MlZf3LQc+2c2snu5m4WVRf9TCPW1Fp2fenLmFKRno985FUbmpq4hC/e/0WuW3cd3YVuLph9DufJcl63/By6++fvN88dm+/g+nXXc/fWu0ldSjEssqxrGfW0ztqRtZw550xWdK/gW2u+hapi1fKVC77COQPZSL8r7rqCH234EW+Z/xa6C90YMTw6+ChrR9aiKOfPP5+Prvwo3YVurl13LdeuvZax1hid+U66Cl1ceNyFXLb8Mmr5Gs20yU13/leefe5WLjru7Zz6pv+CVPqygQHjGxCEBdUFGDlAS2JsK/zwU7D25pe29a2AZRdnZsTZJ2f9UWmTXYPbCFdfT3XttQSNoez61RZRDzooDa8hIpvv9QJz2Vg5lWoxx+x4E71T6wjTSa7lQm6KT+X/RH/FDu3m/fHned2yJXz4zYsYb6R8/ntPUE0G+Ujng5w3eRNz2EWTPBsWvY8Fv/JZSt1zaSaWRzYM8bXbHmf9pm1Upc7cfINfXV7mzIU1OjpqFMtVpH9F5v3ip8A59TGUPIeEVzyHwR7F82f3/xk3PHcD977/3n3NP23Gb7qJbVd8hvHeIsVdE4ycNJ/wyj/kjLln/lQj2ADs6ChSLGLy+YMnBjSO2XrFFUzcfAsAtV97D3P++I+RMCTZuZOJW26hcPLJFFetQkRoPfccu666isbq1VTOOYfqRRdhSiUm77iTyXvuwdWnCKo1gs5O8suWUly1iomuPObWuxm/7nrsSBZLx9RqBJ01SC1qLfmlS+i67DIqZ58N1rLr/rtZe9cP2OZGeD4YYbAYc8kZv855r38PQanEhrENXPXIVcwvz+NjpUuoP/ggrefW09jwPIObnuWZxTluXSU835tySu8pvGH2G0g15dtrvs1EPEEgAU4dZw+czfGdxzPWGmPrxFYe2PEA5ajMefPP456t9zDaGiWUkFRTlnQtYUnnEh7a8RCDjSwOSikssbRrKeWoTOISBGFl/0reMvAWlvcsZ+PYBlY/+W22Dz1L2H8iQXUeu5q7WTu8lnWj6ygEBQY6BhjoGKC70E1n1EF33GBe/ynM71lOb7EXsTHD63+C3XAfhd0PcPfux3AiHF9ZwHG9y9i14pdZE+V4dPtz9NVHOPWRv6Y/rJDrXw6FTgpJg8rWxwgndwCwtnwaq2ddRGHzj3lrfAcWw6agwIa8MBLA0iRmWZxQPMDzazGsK5zCtjkX0Fr0VnrmL6O/I89wPWZwrE7SrDO3v5eFPSU2jzT4/v1reOHJ+1jVMcoHXlegV8ag2IWrzWcbfQQ9i+mafRy5MGRwosWWkTr5MOCkeVU/P+cXkNe04hGRi4AvAQHwDVX9i73254F/BE4DhoD3qeqG9r7PAh8GLPAJVb31lY4pIouA7wLdwCPAr6tq/ErlO+200/W2e3/E5T+8nBVmLp8fP4fJu+4mt3AhpTPPIL94MY3HHmPqxz9m/KabKZ5+GvO/8hXGbriRnV/8IrM+/3m6P3A5LzyynQeuXUOHG2X+6MOUNz9OMGs29dnLSaqzKBcd1ZKlvmWQTWuG2Zn0EkSGgRP7OP6958H4CKN33cvkk89iymXCvj6kp59GqY+6dDD50MMUVt/DossuppCMMfTVr1K58AJc12w23fU0o+XjiHMd2GovVGqE29ZTsmN0zqsRPnEPxZFNGHVgDIWVK4l6e7Hj49jhIVrPPZ91xgMYQ8cF51O54ALs0DDxls248QkkDMEYJu+9l2TXboI5c9GxMXRqApWQiY75jNUWM1WaAzhEHREpFRmnGkwhO7dgG01UAqSrh2DOXEy5g9bjj2BaUxROPBE541zSBSfiyjUKXY77d/2QcNsWTtnVj31hiGB4J8WxzVCfwlYKbC45Xogcs3OwuDhAhymzsznI5vpWWmmLPlOjU8okPT08P1BjdZcyUhsniFJsHDO4dS21SUf3lFCuK5WmkgTw1HHCswNCmC9yUm4hpzb6MFMNRscGmZwYQsam6Gh2ErkqTyzYypMLUwq5Mid0ncDSrqWMNkZ4eP1dFCdjRsrQKBz8pSyqiIIxhu6wBGGBsWTyxb45g8GoksrLn9UA4YTiLE6uHc/cwgmk9Spxo0HSmKI2vIalk49xnN3G1ijkgbCXx8NO0miKZtjEiXJSw7KoUaSslq3FKR4p5IlF6LOWbid0pTEdzr24FKxhl3SwOh+ypiA0jLC8VeGs4iKW1HoZclMMa51Kq8H8xjhzGuOk3SuI558Fc05lpLmNbRPPM56O09HRT6HST3+uxuJWTGvn84S2QSE0Waur1Jt53uiY03ZtVcg8uI9tyZZkKpvUnK9mHjumdjEyuY18vkapNh/K/S/NJTNhZl6tzNpvnySQPQNpM+tHzFcP3U3WLyivWcUjIgGwFngrsAV4EHi/qj49Lc3HgFNU9aMicinwH1T1fSJyIvAd4AxgLvD/gKXtbPs9poj8C3C9qn5XRL4GPK6qr+iWefacZfre3/oCk7lRTl6/htOfvJfcnD4aE8pIeRlTHYvpSSboy0NxwTyG5nWxcevTOE2pDe6iZ+1TDK98H5vdIgqNYfLFThaXisyJYNxadjW3s6u5iXEtkZjZiOlGTEAhTHGpJWbfFo9qiqSbyJtRIqoYKaHqaKaDNNNt5IzQneuhGtSYsDE7WsM4EXJRhSgNCWJolUu00iZoEzCICFEg9OS6mRN1YbXJzuYGBhvrCAq9lCvHU4z6aUpKc3wrLhkn0QCRAkiIOoc6h7MTqBtC7ShiKkg4BxP0AwKaEEUJzk7hbAPnFDGdBKaGISFwQwQ6hiKoFNtLGaSMSpHU1bFuArTZ7guDQEqYsB8T9GF0isBuJdJBDDkwHRhTwjlLvjVMJdmJ2AbqWiSmRKO4iPH8UlqmiLopxE1iUDqaw1TrgzgXkyDEJiTJ1UjyPaBCZXw9leYmJFSmKFAvVFHTgTGdaNRLs9CPC0LQFElHKE2tQ8PdaJpA6gisEqVKZC05p0RdNZoDnRgJMRMtkkYLS4QzeazLkSaWtB0GIh+C7XIkJSFMQGKLqSeYqRTTtES0KOct5Uqe8VzIqBEm4xSpK1GSB5MjDULSMMSooxhPUmyM0gxH2F0eZqQ6jimViYo9OIo0xkbomEjIJYoLDR2FKq5UYrPGjEQt8pqnnNQopGXUDYPdAW4EF5aohl0YzTFuE/KJIzUNhku7Ga40SQ3kXIW8rYJO4mQcUci1AqqtLvK2gnEGo4YkSGnlpqhFY6SRMBR1MBaWiJylksaUXAsNW7hckzSf0Mgr9ZwgQGfsqLYc9TRgp41IE4MKiFGKRtGcw+UUCZXOWKm1HLU4pK+hdLUgBrZEwtZ8DiMxvZLSjcUYgyv04aJZVKIWnWETwbHZxWy1KVZD+rVGp/SSz9eIygVcLiQJLXGUYE2LqDFJvj5JVG9g0xibpuRNwPxqPwu6j6NQ7WM8DBg30HAJsU1INKVgclTCAqWgQDEsU4pKSFikUahSz5dBchQSIWo5Gnlll5lkd3OIlo1JURBDf3k2czoW0FOehbgUcSkOJQZaAuOtMUamdjLaGKIUlekp9dFV6qMQlckFeUITIgBJPZuXl6uwR0vsad2+lhXPG4E/VtW3t9c/C6Cqfz4tza3tNPeJSAjsAPqAz0xPuyddO9s+xwT+AtgFzFbVdO9zH4ho2Yna87V/enG9nDZYNbmevniSYirkrdBCmVQllZRc0CQKmuQ0JWyBNFOaQZ7JYgfj+TITYY2xsMpkWKRoE2qppStN6IrrVFuTBGmDySiiHoaEJqZDxukIxohSReI81kYMF4ShYp5mGNHdarBgCnrigNFiwnjBkRpH6BKMS0ldSKx5mhrRCvI0gwKpCSnaOqV4ElHLSLGLwUIPqQnp0jFqjDImVTYHA2wPe6nZKQbiQebGu4lcjEtjrArNsEgrKhObCMFi1BKqJecsBeeInBJYCFRoBiH1MEfLBISmQY4pIlqQBmgagAsQURAIsYQaE2mCC4U0EmxoCJySc5YARxoqaQhWBOdyqAsxzpFzMTkXk5CjIQVakgckay2oUrYxFdcilBZp3pJEjpbkaLoyTS2RuByphjgNCV1KziXkXULJxZRtQk4SXJhio4SmiZjQChNSJHBKSVuUXEJKQExIU/JMmiJTpkRqAorUKTNFJDFOAxwBkbOU0xaFNCExAfUgRzMIybmUgovJaYLVEKsRikEytUwoKbmwSS5okrgCU7ZGQ8vExmADRQXyNqGYJuSsfdHpg4pgJcASEAcBLROQmIC8tii6OgXXQI0hNSFWDOIcoorBEhlLSAoYUhvhXERKli7FkBKSSEgqAWSvJQwQYYlICTTFqGsv2b0hzqBicMaQGoMzihVFxJF3MQWNMU5xhDgCRCHAYnBYhBQhwYAIarKzRpqQcymRpgROs7QS0JIcsUSggnEQOM3KhAVRmlKgIQVSCcmpJacODDhRnMnun1xqyVmHOABB1dAKDc3QYMWQ05S8TQmdItaAFQTFqEWwOAMaZMcTNaAGcYIgiCooWFFsu8KMmkxmVcL2M6ZOsRacGlxo0DCAQAmcJXTZNc5OJG1loKgoKoKKvOQBZM97vX0vGjSrJ7EoitWsXq0YUglwIuStJZdawlRJg6zOrJGsTtQi4kAURPnib//Ra3YC6TxgegCLLcCZB0rTVhhjQE97+/175Z3X/r+/Y/YAo6qa7if9AZmt2/kL/ThVxnmGE7k3OJuf1N5ALIfW9zKdQFO6GaKbIRbyAg1KTFBhC52MMBeVn85FS6QxSUcO9jOV5UCIZi+B6eUXtcxmBwVtsl0WMUaNMlMscBtZro8yGnSyubSAx0tn7lPGik5QoIHDYAlIiEjIkUhu71OT1wZ5WsTkaFFAD9SZ/yoi6pB2++iVzifqyNMiJCXAIigxORoUDloveW1iCUjl5SaaQFPKTFJhghwxw8yiTpmEiACH4Ij3c63y2iQm91Nfn5y2iIgJsYDSpEhLDuxDK9IWeVrkSGiSpy4vd3EUaIqVQ3v8jVoiYnLEhKQvKkjXVkgxeSyZwjuQXKLZXWTa9bX39fQcbfzRYeX+eSqe/Rm3925eHSjNgbbv765+pfT7Fkrkt4HfBlgwMJ+zer5ElNvFm7o7+XjHHKKom0YaM9wYZ7w1QT4CXOaRWcIeJOomVWGiPsbY9kE6ImFOf5XeShG0G+cGsK6JOmg0WrRaLSSqM2pyxAi10FAJQMIqDdPLhBZo2oR6PEqcjjNQ6mJRxyxKQcjOVoPVQ1sYnJhgXnUWc6vd5MOQhnU0XPZ1maNFoA1qYUQ1ijASMp402N5s0HKOpR09VHIvTWB1mpmxRF7yJaeqWOeo25jxeJzUtujNlylGFYzJvZg3C1ubkrqYlm3RSmNaNqYjiigEEcbkCYLsBRerktiUejyMdU2MyRFISIqhYaGlQjkqUw1DCoEhcUrTKVaVghEKgSESwbTPHTvHpHVMWUfBCJUgoGDkxbJZVSZSy1hqiZ1SCgzFwFAyhvy0dADOtXAuRdXSsI66UyadkGhIIYiIjFA0hmoYELVHebWcYzJ1L+4LBVQTnMv6YYKgst9O9oZ1jCYpxcDQEQYEIqgqdetoOCVA21/5Nru9RXCEJAqxOgrGUDYGo1MYk8OYlz4qmmmTqbSZmfEkIJCASCAUxZgImaYEUqdMWUvOmJddN6dKokqqSuoUESHInpPs+qOYaX1LqtrOmy17y+xUadmYRjJFK62TCwKKQYFckCMICkhb0SdOqVtLohDJSw+2be8LjRCQtWhCExIZgyrUnaNuHS3ncAoWJRKhIELOZPdnbGNiG6MSkFiDWqUcOUomIRJLSkA9dkBAKVciFxWIbcpofYLxxhRqHCYUMErFRBQ1JLBCS5SGdTRti9S2iJM6DkUlQE1AoVAmnysQSgRisWlCmrZoJQlxHKPOUsoXyYcREghWLYnNLAyJGmKnkCa4pImNYyqlMqVymSCIiC00Y4t1FieK1Ti79lawqSWUgADAOiQIcKqkalGxpM5iVUnIkbogM2OnMWJb5KMIg0GdpelSJtOElrXkRAgBo+AwOCeoBKQEpDbrfD8cfp6KZwswfTzuALDtAGm2tE1tNWD4IHn3t3030CkiYbvVs79zAaCqfwP8DWSj2hat3Le1GEVQLR6kqVHtgtkLXzFJV9crH+Ll1PbZMrtQYva8pftJO53KPlt6oio9BxixbfbzchQRwiCgGhSp5or7yfVSOpGInInIhWU6XqFhmBchb3JUov37t9snvdmfJC+RM4ZuY+g+wIdyIEJnFNIZHfyWNib/4gs8iqB6kPRZ+Qz53Mu/e0RyGLNv6286xcBQ3MsThohQDgPKL245lK//jn22FMJC5tD2EAiNUNvP6EsjktUVZMN19uHl98vBBrAZEYphnmKYJxvns3+iA5TnFRGomoBq+Eqt1DxMu7IHorbX81GM8tQKB8/nyThcxfPztIc8CCwRkUUikgMuBW7YK80NwAfb/38N+DfNOp1uAC4VkXx7tNoS4CcHOmY7z7+3j0H7mD/4Ocrm8Xg8np+Rn1uLp91n87vArWTfUn+vqqtF5AvAQ6p6A/B3wNUisp6spXNpO+/q9ii1p4EU+LhqNo18f8dsn/IK4Lsi8mfAo+1jezwej+c1hp9A2vbV5vF4PJ5D43CHU/tZUh6Px+OZUbzi8Xg8Hs+M4hWPx+PxeGYUr3g8Ho/HM6N4xePxeDyeGeUXelSbiOwCNh7pcrzK9JJNqD2WOdZl9PId9e2o8AAABgFJREFU3Rzr8gEsU9V9ZzUfIr/QAeRVte9Il+HVRkQeOpxhjkcDx7qMXr6jm2NdPshkPJz83tTm8Xg8nhnFKx6Px+PxzChe8Rx7/M2RLsAMcKzL6OU7ujnW5YPDlPEXenCBx+PxeGYe3+LxeDwez4ziFc9RjIjMF5F/F5E1IrJaRH6vvb1bRG4TkXXt358qMtBrDREJRORREfnX9voiEXmgLd8/t0NkHJWISKeIXCsiz7Tr8Y3HYP19sn1/PiUi3xGRwtFchyLy9yIyKCJPTdu23zqTjC+LyHoReUJEXn/kSn5oHEC+/9m+R58Qke+JSOe0fZ9ty/esiLz9UM7hFc/RTQp8SlVXAL8EfFxETgQ+A9yuqkuA29vrRzO/B6yZtv7fgf/dlm+Ew49LdST5EnCLqi4HVpLJeczUn4jMAz4BnK6qJ5GFM7mUo7sO/wG4aK9tB6qzi8niiS0hi3z81Rkq4+HwD+wr323ASap6CrAW+CxA+31zKfC6dp6/FjlIPHm84jmqUdXtqvpI+/8E2UtrHvAu4JvtZN8EfvXIlPDwEZEB4JeBb7TXBTgfuLad5KiVT0SqwDm0Y0epaqyqoxxD9dcmBIrtKMMlYDtHcR2q6l1k8cOmc6A6exfwj5pxP1mk5DkzU9Kfjf3Jp6o/akd3BrifLMozZPJ9V1VbqvoCsB4442Dn8IrnGEFEFgKrgAeAWaq6HTLlBPQfuZIdNlcB/wVw7fUeYHTaQ7CFTNkejSwGdgH/t21K/IaIlDmG6k9VtwJXApvIFM4Y8DDHTh3u4UB1Ng/YPC3dsSDrh4Cb2/9/Jvm84jkGEJEKcB3wn1V1/EiX59VCRN4BDKrqw9M37yfp0To0MwReD3xVVVcBUxzFZrX90e7reBewCJgLlMnMT3tztNbhwTiW7ldE5A/ITPzf3rNpP8kOKp9XPEc5IhKRKZ1vq+r17c079zTn27+DR6p8h8mbgHeKyAbgu2TmmavIzBV73D0NANuOTPEOmy3AFlV9oL1+LZkiOlbqD+BC4AVV3aWqCXA9cBbHTh3u4UB1tgWYPy3dUSuriHwQeAdwub40D+dnks8rnqOYdn/H3wFrVPV/Tdt1A/DB9v8PAj+Y6bK9GqjqZ1V1QFUXknVg/puqXg78O/Br7WRHs3w7gM0isqy96QLgaY6R+muzCfglESm179c9Mh4TdTiNA9XZDcBvtEe3/RIwtsckdzQhIhcBVwDvVNX6tF03AJeKSF5EFpENovjJQQ+oqn45ShfgzWTN2ieAx9rLJWT9ILcD69q/3Ue6rK+CrG8B/rX9f3H75l4PXAPkj3T5DkOuU4GH2nX4faDrWKs/4E+AZ4CngKuB/NFch8B3yPqrErIv/g8fqM7ITFFfAZ4DniQb3XfEZfgZ5FtP1pez5z3ztWnp/6At37PAxYdyDu+5wOPxeDwzije1eTwej2dG8YrH4/F4PDOKVzwej8fjmVG84vF4PB7PjOIVj8fj8XhmFK94PB6PxzOjeMXj8RwBRGShiPyn9v8PiMgVh5hvclr+pw6W3uN5LeIVj8czw4jI7wC3An8qIncAdwKHFMfE4zkWCA+exOPxvFqISAfZTP5fAVYAdwCjQEFEqtp28ioivw+8l2yW//dU9Y+OTIk9nlcf3+LxeGYWB+SAKoCqbtAsltLtZH7MEJG3kfm8OoPMpc5pInLOkSmux/Pq4xWPxzODqOoU8BvAfyMztV0pIiWy+CZ7oj6+rb08CjwCLCdTRB7PMYE3tXk8M4yq3iAiT5CZ204HPkWmiP66nUSAP1fVrx+hIno8P1d8i8fjmUFEpCIix7VX94Qr71BVCzzbjmF/K/ChdoA/RGSeiBy1UUg9nr3xLR6PZ2aJgK8DvWSu9DcBl7X33ULmVv4vRWQFcF8WwoZJ4AMc3QHhPJ4X8WERPJ4jgIgsBN6iqv8wbdsc4Juq+rYjVCyPZ0bwisfjOQKISCewUFUfO9Jl8XhmGq94PB6PxzOj+MEFHo/H45lRvOLxeDwez4ziFY/H4/F4ZhSveDwej8czo3jF4/F4PJ4Z5f8DEb/N/XiJb+YAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "for (z1,z2) in zip(zshell[:-1],np.roll(zshell,-1)): \n", + " galmap=tmap[z1]\n", + " #compute mean\n", + " Nmean=np.mean(galmap)\n", + " print(\"bin=[{},{}] mean #gals={}\".format(z1,z2,Nmean))\n", + " #and do fractional difference\n", + " dens_map=galmap/Nmean-1\n", + " #compute power-spectrum on the spehere using a standrd healpix function\n", + " cl=hp.anafast(dens_map)\n", + " plt.plot(cl,label=r\"$z \\in [{},{}]$\".format(z1,z2))\n", + "\n", + "plt.xlim(1,120)\n", + "plt.legend()\n", + "plt.xlabel(r\"$\\ell$\")\n", + "plt.ylabel(\"power spectrum\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From c024c70f80d1994e86db998fad38c60fa00f4e50 Mon Sep 17 00:00:00 2001 From: Julien Date: Fri, 8 Feb 2019 13:31:14 +0100 Subject: [PATCH 2/3] Update the README with instructions --- README.md | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index d3a4d58..4a15ff9 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,28 @@ -# 1807.03078 -Supplementary material for arXiv:1807.03078 +# Analysing billion-objects catalogue interactively: Apach Spark for physicists + +This repository contains supplementary material for arXiv:[1807.03078](https://arxiv.org/abs/1807.03078). + +## How to run the notebook + +You must have [Apache Spark](http://spark.apache.org/) and [Jupyter notebook](https://jupyter.org/) installed on your machine or your cluster. + +### On a local machine + +```bash +PACK="com.github.astrolabsoftware:spark-fits_2.11:0.7.2" +PYSPARK_DRIVER_PYTHON_OPTS="jupyter-notebook" pyspark \ + --master local[*] \ + --packages $PACK +``` + +### On a cluster + +Standalone mode: + +```bash +PACK="com.github.astrolabsoftware:spark-fits_2.11:0.7.2" +PYSPARK_DRIVER_PYTHON_OPTS="jupyter-notebook --debug --no-browser --port=$PORT1" pyspark \ + --master $SPARKURL \ + --packages $PACK \ + --driver-memory $MEMDRIVER --executor-memory $MEMEXEC --executor-cores $EXECCORES --total-executor-cores $TOTALCORES +``` \ No newline at end of file From b72e69999cacd8450f02b55210a8aff2b490d771 Mon Sep 17 00:00:00 2001 From: Julien Date: Fri, 8 Feb 2019 13:34:04 +0100 Subject: [PATCH 3/3] Add instructions for NERSC users --- README.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 4a15ff9..1e4c25c 100644 --- a/README.md +++ b/README.md @@ -25,4 +25,8 @@ PYSPARK_DRIVER_PYTHON_OPTS="jupyter-notebook --debug --no-browser --port=$PORT1" --master $SPARKURL \ --packages $PACK \ --driver-memory $MEMDRIVER --executor-memory $MEMEXEC --executor-cores $EXECCORES --total-executor-cores $TOTALCORES -``` \ No newline at end of file +``` + +### DESC members: working at NERSC + +Source your DESC environement. Then go to the Jupyter Lab web [interface](https://jupyter-dev.nersc.gov/), and execute the notebook with the desc-pyspark kernel. \ No newline at end of file