DES_Comparison_check.ipynb

{
 "cells": [
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "from obscond import SkyCalculations",
   "execution_count": 1,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "from lsst.sims.photUtils import BandpassDict, Bandpass",
   "execution_count": 2,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import numpy as np",
   "execution_count": 3,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import pandas as pd",
   "execution_count": 4,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "%matplotlib inline\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set_style('whitegrid')",
   "execution_count": 5,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "class SkyMagCalcs(SkyCalculations):\n    def __init__(self, deshwfname=None):\n        self.deshwfname = deshwfname\n        self._deshw = None\n        super(SkyMagCalcs, self).__init__(hwBandpassDict=BandpassDict.loadBandpassesFromFiles()[1])\n\n    @property\n    def deshw(self):\n        if self._deshw is not None:\n            return self._deshw\n        else:\n            deshw = self.desHWBandpassDict(self.deshwfname)\n            self._deshw = deshw\n            return self._deshw\n\n    @staticmethod\n    def simlib(fname):\n        df = pd.read_csv(fname, skiprows=2,\n                         delim_whitespace=True,\n                         names=['trash', 'ROW','LIBID','RA','DEC','MJD','BAND','ZP_pe','SKYMAG','PSF','M5SIG'],\n                         index_col='ROW')\n        del df['trash']\n        return df\n\n    def desHWBandpassDict(self, fname=None):\n        if fname is None:\n            # private bandpass currently\n            fname = '/Users/rbiswas/data/DES/Kessler_DES_Cadence/STD_BANDPASSES_Y3A1_FGCM_20170320.dat-2.txt'\n        desdf = pd.read_csv(fname, delim_whitespace=True)\n        desdf.rename(columns={'#lambda':'wave'}, inplace=True)\n        desdf.wave = desdf.wave / 10.0\n        df = desdf[['g', 'r', 'i', 'z', 'Y']].divide(desdf['atm'], axis=0).join(desdf['wave'])\n        bandpassNameList= ['g','r','i','z','y']\n        bandpassList = list(Bandpass(df.wave.values, sb=df[bandName].values) for bandName in ['g','r','i','z','Y'])\n        hwBandpassDict = BandpassDict(bandpassList, bandpassNameList)\n        return hwBandpassDict\n    \n    def desSkyMags(self, pointings, raCol='RA', decCol='DEC', mjdCol='MJD', bandCol='BAND'):\n        \n        num = len(pointings)\n        idxs = np.zeros(num, dtype=int)\n        des_mags = np.zeros(num)\n        lsst_mags = np.zeros(num)\n        sm = self.sm\n        \n        count = 0\n        for obsHistID, row in pointings.iterrows():\n            ra = row[raCol]\n            dec = row[decCol]\n            bandName = row[bandCol]\n            mjd = row[mjdCol]\n            sm.setRaDecMjd(lon=ra, lat=dec,\n                           filterNames=bandName, mjd=mjd,\n                           degrees=False, azAlt=False)\n            lsst_mags[count] = self.skymag(bandName, sm=sm,\n                                           hwBandPassDict=self.adb.hwbandpassDict)\n            des_mags[count] = self.skymag(bandName, sm=sm, hwBandPassDict=self.deshw)\n            idxs[count] = obsHistID\n            count += 1\n        return pd.DataFrame(dict(obsHistID=idxs, lsst_mags=lsst_mags, des_mags=des_mags)).set_index('obsHistID')",
   "execution_count": 6,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "smc = SkyMagCalcs(deshwfname='/Users/rbiswas/data/DES/Kessler_DES_Cadence/STD_BANDPASSES_Y3A1_FGCM_20170320.dat-2.txt')\nsmc.adb.hwbandpassDict",
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 7,
     "data": {
      "text/plain": "<lsst.sims.photUtils.BandpassDict.BandpassDict at 0x1a21bb1a90>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "plt.plot(smc.deshw['g'].wavelen, smc.deshw['g'].sb, color='r')\nplt.plot(smc.adb.hwbandpassDict['g'].wavelen, smc.adb.hwbandpassDict['g'].sb, color='g')",
   "execution_count": 8,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 8,
     "data": {
      "text/plain": "[<matplotlib.lines.Line2D at 0x1a21bcfbe0>]"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD7CAYAAACVMATUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8XHWd//HXzCSdJG3SNJmUFgoqSL8ElCKgtlAEFooWFetlXVwvD2GrywNcFPixy7qywC6I61JdRauiD+THroj4w3KpUMDlIrR0wSACdvqtZSstTdMm6SW3ziQzc35/nJkw5DozmcnMOX0/Hw+YnDknZz7nm/bdb75z5vsNOI6DiIj4Q7DcBYiISPEo1EVEfEShLiLiIwp1EREfUaiLiPiIQl1ExEcU6iIiPqJQFxHxEYW6iIiPVE33C7744otOOBye7pctiXg8jl+upRTUPuNT20xM7TPawMBA1ymnnNIy2XHTHurhcJjW1tbpftmSiEajvrmWUlD7jE9tMzG1z2htbW2v5XKchl9ERHxEoS4i4iMKdRERH1Goi4j4iEJdRMRHFOoiIj6iUBcR8RGFurh+/nPYu7fcVYjIFE364SNjTBBYDSwC4sBKa+3WrP3/B/gUkAK+bq1dU6JapVT+8Af46792//vZz8pdjYhMQS499RVAjbV2CXANsCqzwxjTCFwOLAHOA/6jFEVKiW3Z4j4ODJS3DhGZslxCfSmwDsBauxE4NWtfP/AaMDP9X6rYBco0eC396eOjjipvHSIyZbnM/dIAHMjaThpjqqy1ifT2DmATEAJunuxk8XicaDSad6GVKBaLefZaAoODHH711fSffjrhLVtoAjqHhugq4vV4uX1KTW0zMbVP4XIJ9R6gPms7mBXoy4H5wNvS248YY9Zba58b72Sa0KtCXHstPPYYDY89BvPmAdDS0EBLEa/H0+1TYmqbial9Rmtra8vpuFyGX9YD5wMYYxYDL2ft2wccBOLW2hiwH2jMq1KZfr//Pdx8M1xwgbvd0eE+xuPlq0lEiiKXUF8DxIwxG4BvA1cYY640xlxgrX0aeB7YaIx5FtgCPFa6cmVKHAeefx4+8xloaYGf/hTe9a439ivURTxv0uEXa20KuGTE05uz9l8HXFfkuqTYHAe+/GW49VaorYV774WmJrj/fnjhBbjsMoW6iA/ow0eHivvvdwP97/4Odu6E5cvd5488Ej7yEQiHFeoiPjDtKx9JmdxyCxx9NHzrW1A1xo9doS7iC+qpHwpeeAHWr4cvfWnsQAeFuohPKNQPBf/6rzBrFlx00fjHKNRFfEGh7ne33w733Qdf/So0TnC3qUJdxBcU6n722mtwySVw3nlw5ZUTH6tQF/EFhbofPfMMLFoE73sfBALwk5+4oT0RhbqIL+juF7+JxdwpdHfscLdvuMG9bXEyCnURX1Co+80tt7iB/uCDYAwce2xu36dQF/EFhbqf/OAHcP318MlPwoc+lN/3KtRFfEFj6n6xdi1ceiksWwY//nH+369QF/EF9dT9IBZz53VpbYUHHoDq6vzPoVAX8QWFuh/89Kfwv/8Ljz1WWKCDQl3EJzT84gd33QUnnADnnlv4OTKh7jjFq0tEpp1C3asSCXe63NNOc+9Lv/DCqZ0vcx/70NDUaxORstHwi1f95CewerX79YIF8PnPT+18mVCPx2HGjKmdS0TKZtJQN8YEgdXAIiAOrLTWbk3vOwn4j6zDFwMrrLXrSlCrZAwNwTe+AUuWwKOPuoFc6Fh6Rnao19dPfKyIVKxceuorgBpr7ZL0GqWrgI8AWGtfBM4CMMb8JdCuQC+x116Df/s39/HWW93ZF4shO9RFxLNyCfWlwDoAa+1GY8ypIw8wxswEbgDeV9zy5E3+9Cd3TdH+fnjb2+CDHyzeuTNDLoODxTuniEy7XN4obQAOZG0njTEj/zH4G+CX1tquolUmb9beDhdf7N6dcuON8KtfQbCI73OHQu5jKlW8c4rItMulp94DZA+yBq21iRHHfBr4RC4vGI/HiUajOZZX2WKx2LRcS82mTbzlU58iODREx7XXsu9jH3N3FPG1G3bt4ghg65YtDBWptz5d7eNFapuJqX0Kl0uorwc+DNyTHlN/OXunMWY2ELbW7sjlBcPhMK2trXkXWomi0ej0XMvXvw61tfDHPzLv2GOZV4rXePFFAN5+9NHuRGBFMG3t40Fqm4mpfUZra2vL6bhcQn0NsMwYswEIABcZY64EtlprHwAWAn8usE6ZTFcX3HMPfPGLuc+4WIjMUI6GX0Q8bdJQt9amgEtGPL05a//zuHfISCn89Kfum5eXjPwRFJlCXcQX9InSStbXB6tWwdlnu9MAlJJCXcQX9InSSvXb38J3vgO7d7sLR5eaQl3EFxTqleg3v3HnRQf42tdg8eLSv6ZCXcQXFOqVJpGAyy931xV9+ml4y1um53UV6iK+oFCvNA8/7N5//otfTF+gg0JdxCf0Rmmlue02mDcPPvrR6X1dhbqILyjUK8mOHfDQQ+50AFOddTFfCnURX1CoV5Lbb3dDdeXK6X9thbqILyjUK8WDD8JNN8F557kzME43hbqILyjUK8Hrr8Nf/RW8853uikbloFAX8QWFeiW47TaIxeCXv3RvZSwHhbqILyjUyy0ehx/9CM4/H44+unx1KNRFfEGhXm6/+AXs2QNf/nJ561Coi/iCQr2cHAe++11obYVzzy1vLQp1EV/QJ0rL6cknoa0NfvADCATKW0sm1JPJ8tYhIlOinnq5OA7ccAPMnw+f/3y5q1FPXcQn1FMvl8cfh6eecodfamrKXY0WnhbxiUlD3RgTBFYDi4A4sNJauzVr/3LguvTmC8Bl1lqnBLX6y3XXwYIF8IUvlLsSl3rqIr6Qy/DLCqDGWrsEuAZYldlhjKkH/h34kLV2Me5apZES1Okvzz0H69fD1VdXRi8dFOoiPpFLqC8F1gFYazcCp2btOw14GVhljHka2G2t7Sx6lX5z661QX18ZY+kZCnURX8hlTL0BOJC1nTTGVFlrE7i98rOBk4A+4GljzLPW2i3jnSwejxONRqdSc8WIxWJ5X0to3z6Ovftu9l14Ibt37oSdO0tUXX5mbNvGMcDr27fTW6SfTyHtc6hQ20xM7VO4XEK9B6jP2g6mAx2gG3jeWtsBYIz5LW7Ajxvq4XCY1tbWAsutLNFoNP9rueceSCRo+tKXaKqkdkjfUrng8MPd++aLoKD2OUSobSam9hmtra0tp+NyGX5ZD5wPYIxZjDvcMvw6wDuMMRFjTBWwGNiUX6mHmCefhFmz4JRTyl3Jm2n4RcQXcumprwGWGWM2AAHgImPMlcBWa+0Dxph/BB5JH3uPtfaVEtXqD088AWecAVUVdjepQl3EFyZNFmttCrhkxNObs/bfDdxd5Lr8qb0dNm+Gv/mbclcymkJdxBf0idLp9Ej6F5ply8pbx1gU6iK+oFCfTg89BEccASeeWO5KRlOoi/iCQn26DA3Bo4+686aXe/KusSjURXxBoT5d1q+Hnh5YvrzclYxNoS7iCwr16fLww1BdXf5508ejUBfxBYX6dHAcePBBWLrUnR6gEinURXxBoT4d2togGoULLyx3JeNTqIv4gkJ9Otxxhzsb4yc/We5KxqdQF/EFhfp0uPdeuOACaGwsdyXjU6iL+IJCvdQOHICOjsqb62UkhbqILyjUS+0//9N9PPbY8tYxGYW6iC8o1Evt0Ufdx9NOK28dk1Goi/iCQr3UtmyBj38cDjus3JVMTKEu4gsK9VI6cABefRWOO67clUxOoS7iCwr1UvrVryCRcO98qXQKdRFfUKiX0l13wTHHwLvfXe5KJpcJ9WSyvHWIyJRMukiGMSYIrAYWAXFgpbV2a9b+7wKnA73ppz5irT0w6kSHmr4+ePxxuOaaypyVcaRMjeqpi3haLmuqrQBqrLVL0muUrgI+krX/ZOD91tquUhToWTt3ugHplcVzAwG3t65QF/G0XIZflgLrAKy1G4FTMzvSvfhjgduMMeuNMReXpEov6uhwH+fPL28d+VCoi3heLj31BiB7OCVpjKmy1iaAmcCtwLeAEPCEMeZ31tqXxjtZPB4nGo1OpeaKEYvFxr2WOb/5DfOAV/v7GfTI9ZpAgL179tBZpHonap9DndpmYmqfwuUS6j1A9nyxwXSgAwwA37HWDgAYYx7HHXsfN9TD4TCtXhmSmEQ0Gh3/Wn79a6iq4pizz67c6XZHCoWINDURKdLPZ8L2OcSpbSam9hmtra0tp+NyGX5ZD5wPkB5Tfzlr30LgGWNMyBhTjTtU80J+pfqM48Att8Dvfw//9E/eCXTQ8IuID+TSU18DLDPGbAACwEXGmCuBrdbaB4wxPwM2AkPAndbaP5auXA94/nm4+mo46yy4/PKSvETKSTEwNEA4FCaWiLFt/zZe2fMKfYN9xBIxQoEQVcEqUk6KfbF97D24l70H99I10IWDw96De9nTv4fO/k6Obzme9RevJ6A3SkV8YdJQt9amgEtGPL05a/83gW8WuS7v2rDBfbzrLmhqKtppX9nzCmuia7jP3scLu9xfhgIEcHAm/d666jrm1MyhqbaJUDBEc20zpx5+Kk/++Umeff1Z2nvbOaLhCIW6iA/k0lOXfDz/PBxxRFHueuns72T186v5+Ss/x3ZbAgQ4ef7JfHXpV2msaaQn3kN9uJ7D6w/n5Pkn01jTSDgUJuWkSKQSBAIB5tTMIVwVHvP8T2x7gr+48y/Y1LlJoS7iEwr1Ytu0CU48cUqneHn3y3zvue9x50t3EkvEOOdt53D5ey/n460f57BZxZsY7PiW4wHY1LmJZccsU6iL+IBCvZgcx53A64wzCvr2Xb27uPShS7lv832EQ2E+t+hzXLH4ClpbSnMXwNyZc5kdns2f9v7JfUKhLuJ5CvVi6uqC3l53vpc8rd2yls+t+RwHEwe58ewb+eIpX6RlZksJinxDIBBgYfNCtnRvcZ9QqIt4nkK9mLamp8TJI9RTToqbfnsT1z15HSfNO4m7P3E3C5sXlqjA0RY2L+SZ7c+4Gwp1Ec/TLI3FtHatG4w5rkeaSCX41L2f4p+f/Gc+feKnWX/x+mkNdHBDffuB7cQSMYW6iA8o1IslmYQ77oDly3O688VxHL6y7ivc88d7+MY53+DOFXdSW11b+jpHWNi8EAeHV/e+qlAX8QGFerE8+ii0t8PFuc1p9p3/+Q7ff/77XLXkKv5h6T+4H/4pg8xvBlu6tyjURXxAoV4st98OkQh86EOTHvr0a09z1aNX8dHjPso3l5X3c1vHNh0L4N4Bo1AX8TyFejHEYu54+oUXwowZEx7a2d/JhfdeyDFzjuGOFXcQDJT3R1Afrmdm9Uz29O9RqIv4gO5+KYYHH3SDffnyCQ9LOSk+u+azdA908+uVv6Yh3DBNBU6sua6ZroEuhbqIDyjUi2HVKmhogDPPnPCwb67/Jo+8+gg//OAPOWneSdNU3OQidRG6D3Yr1EV8QMMvxWAtfOYzMHPmuIe8tPslrn3iWj55wif54ilfnMbiJtdcq566iF8o1Kdq3z7Yvx+OPnrcQxzH4bKHLmN2eDarz19dtjtdxhOpiyjURXxCwy9TtXGj+/jOd457yN2v3M0z25/htg/dRnNd8zQVlrtIXYTugW4IHqFQF/E49dSnasMGCIVg6dIxd/cP9nP1Y1dz8vyTufhdlbkud3NtM/ti+0iEAgp1EY+btKdujAkCq3HXHo0DK621W8c45tfA/dbaH5ai0Ir1+OPwjndAXd2Yu7/+9NfZ2buTX3ziF4SCoWkuLjeRuggAe8Mp5irURTwtl576CqDGWrsEuAZYNcYxNwLFW+bHI6r//Ge3p75s2Zj7t+3bxi3P3sJnT/wspx91+vQWl4dMqHfXpNzpDkTEs3IJ9aXAOgBr7Ubg1OydxphPACng4aJXV+Ea1q1zv/jyl8fcf/1T1xMMBLn5nJunsar8Zcb5u2pSGn4R8bhcQr0BOJC1nTTGVAEYY94B/DXwzyWoreI1rFvnLoixYMGofdHOKP/10n9x2bsvc5eKq2CZnrpCXcT7crn7pQeoz9oOWmsT6a8/BxwBPA68FRg0xvzZWrtuvJPF43Gi0WiB5VaQVIrjXn2VrrPOonOM67liwxXUhGr4aMtHK/569w3sA2BXIEZfby87ilRvLBar+GsvF7XNxNQ+hcsl1NcDHwbuMcYsBl7O7LDW/n3ma2PM9UDHRIEOEA6HaW0tzfJs06qrC5JJIiecQGTE9bzY8SKPvP4IXzvja5z+rsodS89469BbYS30zAoyq6+2aD+faDTqj591CahtJqb2Ga2trS2n43IJ9TXAMmPMBiAAXGSMuRLYaq19oPASPW73bvfxsNELQV/7xLU01jRy1WlXTXNRhamtrqWuuo6uGQkNv4h43KShbq1NAZeMeHrzGMddX6SavGGcUN/4+kbWblnLTX9xE401jWUorDDNtc10zxhUqIt4nD58VKhxQv36J6+npa6Fy997eRmKKlykLqKeuogPKNQLNUaoJ1NJnvzzk3zmxM8wa8asMhVWmEhdhK7qIYW6iMcp1Au1z71jhMY3hli2H9hOPBnn+Jbjy1RU4ZrrmumuVk9dxOsU6oXq6yNVW+vObJi2uct9q+G4yHHlqqpgkdoIXdUaUxfxOoV6oXp7SY2YP912WwBMsylHRVMSqYuwrypBwtE0ASJeplAvVG8vqRGTeNkuS1Nt0/AnNL0kM1XAvtBgmSsRkalQqBeqr29UT33L3i0sbF5YcYtg5GJ4qoCqoTJXIiJToVAvVG8vyRGh3t7bzoKG0fPAeMFwqFerpy7iZQr1QvX14YwYfuno62DezHllKmhqmmvd4ZfuavXURbxMoV6oET31WCLG/th+5s3yZqi/0VNPTHKkiFQyhXqhenvf1FPf3ed+GMmroT48p7p66iKeplAvRDIJu3a9qafe0dcBwPz6+eWqakrqquuoTYXonqGeuoiXKdQLcf/9AIT6+oaf2tW3C/BuTx0gkgy787+IiGcp1AvR1QXAgfPPH34q01P3cqg3J8N0hfXhIxEvU6hPweDb3z78dUdfBwECtNS1lLGiqYmkwnTPUKiLeJlCvRAJd4jCyZr3paOvg0hdhOpQdbmqmrJIqkY9dRGPU6gXIpkOvqo31hjZ1bfLs2+SZjSnatzFp0XEsyZd+cgYEwRWA4uAOLDSWrs1a/9lwOcBB/gXa+3a0pRaQTI99VBo+KmOvg5Pj6cDRJxa9s9IkUglqArmstKhiFSaXHrqK4Aaa+0S4BpgVWaHMSYCXAqcBpwD/MAY472JT/I1zvCLH0LdCcC+g/vKXYqIFCiXUF8KrAOw1m4ETs3ssNZ2AYustUPAPGC/tdYpRaEVZcTwi+M4np4iIKPZqQWg+2B3mSsRkULl8jt2A3AgaztpjKmy1iYArLUJY8yXgBuA7052sng8TjQaLajYStG8axdzgYNDQ0SjUfbH9zOYHCTQH/D0tYV7kzAHfrfpdzgtU/+3ORaLebo9SkltMzG1T+FyCfUeoD5rO5gJ9Axr7feMMbcBDxtjzrbWPjHeycLhMK2trYVVWynmzAGgZuZMWltb2dS5CYCT3n6Sp68tVu+utzpr7ixaj5v6dUSjUU+3RympbSam9hmtra0tp+NyCfX1wIeBe4wxi4GXMzuMMQa4Gfg4MIT7Rqr/b59IJt1l7NLzpvvhg0cAzQF3LpvuAQ2/iHhVLqG+BlhmjNkABICLjDFXAluttQ8YY/4APIt798vD1tqnSlduhUgk3nQ7o19CPYIb6l0DXWWuREQKNWmoW2tTwCUjnt6ctf8G3PH0Q0ciAVm3M+7q9f68LwB1wTC1Qwp1ES/Th48KkUyO6qmHQ2Fmh2eXsagiCAZpPhjQ3S8iHqZQL8TI4Zf+DubXz/fk2qRvEgwSOaieuoiXKdQLMWL4xQ8fPALcnvqAQl3EyxTqhRgx/LKrd5dvQj0yoA8fiXiZQr0QY/XUPf5pUsAN9X5HPXURD1OoFyJrTH0wOUj3wW7f9NSbB9y5X5IpTcEr4kUK9UJkDb/s6d8DeP92RmB4+MXBYV9Mk3qJeJFCvRBZwy9++eARMBzqoDdLRbxKoV6IrOEXv4V680H3S4W6iDcp1AuRNfyyu2834J9Qz/TUNf+LiDcp1AsxxvDLYbMOK2dFxREKafhFxOMU6oUYMfzSWNNITVVNmYsqgqoqmhXqIp6mUC9E1vBLR38Hh830QS8dIBSibghqQjX6AJKIRynUC5E1/LK7b7c/xtMBQiECQHPNHPXURTxKoV6I7J66X+Z9geF/qCI1TQp1EY9SqBdixJi670I9PEfDLyIeNekiGcaYILAaWIS7XN1Ka+3WrP1XABemNx9KL5rhb+nhl4OJg/QO9vpqTB2geUYjL/ZuKXMxIlKIXHrqK4Aaa+0S4BpgVWaHMeZo4NPAacAS4DxjzImlKLSipIdfumNub9Z/PfVGDb+IeFQuob4UWAdgrd0InJq1bwfwAWttMr3sXTUQK3qVlSY9/NIVc4PPd6E+o1GTeol4VC4LTzcAB7K2k8aYKmttwlo7BHQZYwLAvwO/t9ZO+Ht7PB4nGo0WXnEFeFt/P4MDA7T3tgMw0DlANOHtawJo6OjgCMA5MISDw/+89D/MCc8p+HyxWMzzP+tSUdtMTO1TuFxCvQeoz9oOWmsTmQ1jTA1wO9ALXDrZycLhMK2trfnWWVmqqqiZM4eeZA8AS965hMPrDy9zUUXw0ksAHDf/aNgOTQuaOC5yXMGni0aj3v9Zl4jaZmJqn9Ha2tpyOi6X4Zf1wPkAxpjFwMuZHeke+v3AH6y1f2utPTR+X88afgkQoKWupdwVFUdm+KXKXUBb87+IeE8uPfU1wDJjzAYgAFxkjLkS2AqEgDOBsDFmefr4f7TWPluSaitF+u6XrlgXkboI1aHqcldUHJm7X6rcX8z0ZqmI90wa6uk3QC8Z8fTmrK99MOlJntJ3v3TFdvrnTVJ4o6ceagAU6iJepA8fFSI9/NId98kydhnDoa6euohXKdQLkTX84ospdzPSoV7nVFFTVaNQF/EghXohkkmcKjfU5830UU89PfVBIJWipa6FzoHOMhckIvlSqBcikaCnKkU8Gffl8AvJJJG6iHrqIh6kUC9EIkF7tbuYpy/uT8/ICvWWmeqpi3iRQr0QySQ7q9xQP6LhiDIXU0TqqYt4nkK9EIkE7SF33Tff9tTrWujsV09dxGsU6vlyHEgmaQ/2AzB/1vwyF1REI3rqvYO9xBPx8tYkInlRqOcrlQKgPdhPfXU9M2fMLHNBRTSipw66V13EaxTq+Uq4c5m1B/qYWzu3zMUU2Yg3SgG9WSriMQr1fKVDfSe9vg71SF0EUE9dxGsU6vlKuhNRttPL3Br/hnpm+EVvlop4i0I9X4kEqQDscnpoqfXJlLsZ6qmLeJ5CPV+JBN21MESSw2p9NO8LvCnUm2qbCBDQmLqIxyjU85VM0p5eB8rPY+qhYIim2ib11EU8RqGer0SCne504/4N9fSbwZoqQMR7Jl0kwxgTBFYDi4A4sNJau3XEMS3ABuCd1tpYKQqtGInEcE+9pcZnY+rpWRozbwZrqgAR78mlp74CqLHWLgGuAVZl7zTGvB94FPDZAPM4soZf/PxGKaCpAkQ8KJdQXwqsA7DWbgROHbE/BZwL7C1uaRUq3VOfG2qgOuiTtUkzRoS6euoi3pNLqDcAB7K2k8aY4WEba+1j1tpDZ9n5RIKd9XB4dVO5Kym+MXrqXQNdpJxUGYsSkXxMOqYO9AD1WdtBa22i0BeMx+NEo9FCv73swn/6E+310JSsIxaLefpaRgr292OA3e3t7I1GSfYmSTpJnnvpOWbPmJ33+fzWPsWktpmY2qdwuYT6euDDwD3GmMXAy1N5wXA4TGtr61ROUV4DA7TXwykNC6ipqfH2tYw04E4nfFgkwmGtrbxj6B3wB2ha0MTC5oV5ny4ajfqrfYpIbTMxtc9obW1tOR2Xy/DLGiBmjNkAfBu4whhzpTHmginU51lDgzF2z4L54eZyl1J8Y4ypgz5VKuIlk/bUrbUp4JIRT28e47i3FqmmitY+sBsnAEfV+Ght0owxxtRB87+IeIk+fJSnHf3tABzptykCYNyeuj6AJOIdCvU8bT/YAcBRdT5a8SgjEHD/y/TUZ2qhDBGvUajnafvALgCOnOWjtUmzhULDoV5XXUdtVa2GX0Q8RKGepx2x3TQNwMwZs8pdSmlkhTq4vfWug+qpi3iFQj1P2+N7OOoAb8yT4jcjQ11TBYh4ikI9T9vjnRzZg39DvaoKhoaGNzVVgIi3KNTztGOo0+2pZ+4U8Zv6eujtHd7U9Lsi3qJQz0NvvJd9yX5/D7/MmQP79w9vRmrVUxfxEoV6Hnb07ADgSL+H+r59w5stM1voG+wjlvD3NPkifqFQz8OOA26o+3r4pbHxzT31zAeQ9GapiCco1POw/cB2AP8Pv2T11I+afRQA2/ZvK1dFIpIHhXoeth/YTpAA8/s4ZEL9uMhxAGzuGjXdj4hUIIV6Hrbt38aCYCNVKfw9/NLTM3yv+lGzj6K2qpZop+a2FvEChXoeXux4keND6dkZ/dxTBzjgLnYVDAQxEcPmbvXURbxAoZ6j3X27+WPnHznz1SQEgxAOl7uk0mhsdB+z3ixtjbSqpy7iEQr1HD2+7XEAzvnt67BiBczy6dwvmZ76iHH11w68xsDQQJmKEpFcKdRz9N/b/pvG6gZOfnUAPvzhcpdTOk3pBbX37Bl+qjXiLitmu2w5KhKRPEw6MGyMCQKrgUVAHFhprd2atf8LwN8CCeBGa+3aEtVaVo9ve5yzqt9OyHkB3vOecpdTOieeCNXV8NRTsHw5AK0tbqhv7trMu+a/q5zVicgkcumprwBqrLVLgGuAVZkdxph5wOXA6cD7gZuNMf4ZbHYc2LuXbc88yLb92zhn7SZoaQFjyl1Z6TQ0wBlnwJo10NkJjsOxTccSDASJdmlcvWCOA4kExGLQ10ewpwe6u93fiLq63OGunh7o73ePGRx070BynHJXLh6Tyy0cS4F1ANbajcaYU7P2vQdYb62NA3FjzFbgROD5Yhf6u/bfsfL4VKnXAAAE50lEQVSBlQwmB3FwcNJ/2B3SjxNs53LsmPv278cZGGBXvVvDuQeaYdXN/r2dMWPlSvj0p2HuXKipIRwOc/RFcOujN3LvfTdDIPvgwOjvTz/lOA6BsfYPmySwJtw9wc5Jc3C6v9fJ4byTCAz/b4yvJ9jO6bxTKqok3+o4DoHAeAcV+LqFlhsIwPz5UFdX4AlcNVU13LniTk6Ye8KUzjOZXEK9ATiQtZ00xlRZaxNj7OsFZk90sng8TjSaf4+vq7eLBeEFJFIJAAKBwHBgDD+O+EOQ/Xwux2aeyxxTFdhLyNkL1bXMO/x4UnddRTQQgHT9sVisoGupeCedRPhXv2LmM89QtXcvgaEhrji4lUeCr0MIAo7DqJRysr5If/3GX8yRgTOe0cc5I/eNcfhk58j1tcY88YTnGX2cM9b+zDKBQffRCQTctgmF0sc4bo/ccdsu8KbtdHumOxuB4fYd4xHeOH5MY1Y30SETPZmb8c432b/njkPAGavCAmuZynUFAsSrDycVntrNEeFQmI7tHQS7S/tWZi6h3gPUZ20H04E+1r56YD8TCIfDtLa25lUkQCutfOA9H8j7+0opGo0WdC2e0Nrq3uWTdmn6v3z4un2mSG0zMbXPaG1tbTkdl8s/GeuB8wGMMYuBl7P2PQecYYypMcbMBlqBV/IrVUREiiWXnvoaYJkxZgPub2wXGWOuBLZaax8wxnwXeBr3H4h/stZqjlYRkTKZNNSttSngkhFPb87a/2Pgx0WuS0RECqAPH4mI+IhCXUTERxTqIiI+olAXEfERhbqIiI8EnGmeW6Ktra0TeG1aX1RExPvecsopp7RMdtC0h7qIiJSOhl9ERHxEoS4i4iMKdRERH1Goi4j4iEJdRMRHcpml8ZBnjJkLtAHLcNdivQN3hv1XgMustSljzHXAB9P7v2Ktfa5M5U4rY8w/AhcAM3DXsn0KtQ/GmGrg/wJvBZLAF9CfHQCMMe8F/s1ae5Yx5u3k2CbjHVuOa6hk6qlPIv2X80fAwfRT3wK+Zq09A3cq4o8YY04GzgTeC1wIfL8ctU43Y8xZwGm4a9SeCRyJ2ifjfKDKWnsa8C/ATahtMMb8PfAToCb9VD5tMurY6azdKxTqk7sF+CHQnt4+Bbc3CvAwcC7uOq6PWmsda+12oMoYM+mHBHzg/biLpqwBHgTWovbJ2IJ7nUHcZR+HUNsAvAp8LGs7nzYZ61gZQaE+AWPM54FOa+0jWU8HrLWZT2xl1mTNe61Wn4gApwJ/iTvn/s9wlztU+0Af7tDLZtz1Br6L/uxgrb0X9x+4jHzaZKxjZQSF+sQuxl316UngJOBOYG7W/syarHmv1eoT3cAj1tpBa60FYrz5L9qh3D5X4LbNQmAR7vj6jKz9h3LbZMseE5+sTcY6VkZQqE/AWvs+a+2Z1tqzgBeBzwEPp8eSAZbjLuW3Hni/MSZojDkKt7faVY6ap9kzwAeMMQFjzOHATOC/1T4A7OON3uZeoBr4vdpmlHzaZKxjZQTd/ZK/q4AfG2NmAFHg/1lrk8aYp4Fncf+hvKycBU4Xa+1aY8z7cBcgz1z3NtQ+AN8Gbk9f9wzgq8DvUNuMlM/fp1HHlqPgSqcJvUREfETDLyIiPqJQFxHxEYW6iIiPKNRFRHxEoS4i4iMKdRERH1Goi4j4iEJdRMRH/j+RgHwZWzTPRwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "pointings = SkyMagCalcs.simlib('/Users/rbiswas/data/DES/Kessler_DES_Cadence/SIMLIB_DUMP_DES-griz.DAT')\npointings[['RA', 'DEC']] = pointings[['RA', 'DEC']].apply(np.radians)",
   "execution_count": 9,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "from joblib import delayed, Parallel",
   "execution_count": 10,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "pts = np.array_split(pointings, 40)",
   "execution_count": 11,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "# This dataframe was calculated using the class above. Once done, we can use the saved values for making plots\n#dff = smc.desSkyMags(pointings)\n#dff.to_hdf('DES_pointings.hdf', key='0')\ndff = pd.read_hdf('DES_pointings.hdf')",
   "execution_count": 13,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df = dff.join(pointings)",
   "execution_count": 14,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df['diff_des'] = - df.des_mags + df.SKYMAG\ndf['diff_lsst'] = - df.lsst_mags + df.SKYMAG",
   "execution_count": 15,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.diff_des.mean()",
   "execution_count": 16,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 16,
     "data": {
      "text/plain": "-0.06865641764850226"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, ax  = plt.subplots()\ndf.diff_lsst.hist(histtype='step', bins=np.arange(-8, 8, 0.05), ax=ax, color='k', alpha=1, lw=1, label='LSST: median {0:0.2f}'.format(np.median(df.diff_lsst.dropna())))\ndf.diff_des.hist(histtype='step', bins=np.arange(-8, 8, 0.05), ax=ax, color='r', alpha=1, lw=1, label='DES: median {0:0.2f}'.format(np.median(df.diff_des.dropna())))\nplt.legend(loc='best')\nax.set_xlabel('Model - DES')",
   "execution_count": 17,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 17,
     "data": {
      "text/plain": "Text(0.5,0,'Model - DES')"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8VOW97/HPkJAECBdBCkpkBsrhJypoAVsEBKy3im2puq3UFsTL7tF6qRR2FV+g8dQLKqJWRauiuKn1AurpPlWUg8cLUtHdVLy0w4NVMxjZWwELCpqEJHP+WDNxEibJJJnJTLK+79crLzNrnjXrN4l858mznvWsQDQaRURE/KFbtgsQEZGOo9AXEfERhb6IiI8o9EVEfEShLyLiIwp9EREfUeiLiPiIQl9ExEcU+iIiPpKf7QIa27RpU7SwsLDN+1dVVdGe/TNFdbWO6mqdXK0Lcre2rlbXl19+uWPcuHEDW2qXc6FfWFjIqFGj2rx/OBxu1/6ZorpaR3W1Tq7WBblbW1erq6ysLJJKOw3viIj4iEJfRMRHFPoiIj6Sc2P6ItI6+/bto6amhnA4nO1Sktq3b19O1tZZ6yoqKqKkpITu3bu36fUV+iKdXEVFBQMGDOCggw4iEAhku5z9fPXVV/To0SPbZeynM9YVjUbZuXMnFRUVDBs2rE2vr+EdkU6usrKSfv365WTgS3oFAgEGDBhAZWVlm19DoS/SBSjw/aO9v2uFvkgXEwqFCAQCafsKhULZfkuSRhrTF+liIpEI6bz3dUs9y9dff53HHnuM2267bb86rr/+eqqrq4lGoxxxxBHMmzePbt268fTTT/P000+Tl5dHNBrlggsuYPLkyZxzzjnU1dXxwQcf0L9/f/r168fEiRO56KKL0vZ+Gjv99NNZtmwZ69evZ+DAgUybNi0tr/vXv/6VG2+8kby8PKZMmcIvfvGLBs/v3LmTefPmUV1dzeDBg7nhhhsoKirigQce4IUXXqBXr178/Oc/Z+rUqWmpJ06hL5KqUIhv7tsHH3+c7Uo6haVLl/Kzn/2Mo48+mqKiIi655BJeeOEFJkyYwLJly3jmmWcoKCjgk08+4cwzz+Sll17i4YcfBuDKK69k+vTpTJkypcPqPfPMM9P6etdccw333HMPQ4YM4fzzz2fz5s0ceuih9c/feeednHbaacyYMYNly5axatUqxo8fz7p161i1ahV1dXXMnDmTCRMmpHW5CIW+SKoiETj4YAiFoLw829XkvIMPPpinn36a/Px8jj76aG6//Xby8/Opq6ujtraWRx99lOOOO46hQ4eybt06unVrerTZOccTTzzBokWL6rfV1NRw6qmnMnr0aLZu3cqkSZPYtWsXb7/9NiNHjuTGG2/k448/5uqrr6a6upqioiKuu+46Bg0axJIlS9iwYQNDhgxh9+7dANx2222UlJRw2mmnsWjRIj755BN2797NtGnTuPTSS5k/fz49e/bk448/Zvv27dx8880NQjzRrl27iEajlJSUADB58mQ2btzYoH1ZWRmXXXYZAFOmTOHuu++mX79+HH300RQUFABwyCGHsGXLFkaPHt2+X0YCjemLtMIw8MJfWjR37lyOPPJI7rzzTiZOnMiCBQv44osvyMvL46GHHiISiXDBBRdw3HHHsXr16mZfy8waBH5cRUUF8+fP5/e//z0PPvggs2fPZtWqVbz22mvs3buXxYsXc95557Fy5Upmz57N0qVLCYfDvPXWWzzyyCPccMMN7Nmzp8Fr/td//Rfjxo3jwQcf5JFHHuGRRx6pf+6QQw5h+fLlzJw5kyeeeKLJevfs2UNxcXH94169evHFF180aLN371569+5d//yePXswM/7yl7/w5Zdf8tlnn7Fp0ya++uqrZn82raWevkgrbNu2LdsldBobN25kzpw5nHXWWdTV1XHTTTexbNkyzj33XCorK7n66qsB+PDDD7ngggsYN24cZtaqY/Tv35/BgwcDUFxcXD93vbi4mKqqKrZs2cKyZcu49957qauro6ioiPfee48jjjiCQCBAnz59GDlyZIPX7NevH5s2beK1116jd+/e7Nu3r/65ww47DICDDjqIv/3tbw32e/jhh1m3bh0AS5YsYe/evfXPJQZ8XK9evdi7dy/9+vWrf37kyJGcccYZnH/++YRCIcaMGcMBBxzQqp9JS9TTF5GMuOWWW9iwYQPgBdywYcMoKChgx44dzJ8/v35YZciQIRxwwAFtusK0pZPMw4YN44orrmDlypWUlpZy8sknM2zYMN555x3q6urYu3cvH3zwQYN9Vq9ezYABA7j11luZPXt2g552c8c755xzWLlyJStXrmTQoEEAfPTRR0SjUV599VXGjx/foP3YsWN5+eWXAXjllVcYP348O3bsYM+ePTz66KMsWLCATz/9lG9+85ut+pm0RD19kS4mGAymdd5+MBhssc2GDRs4/fTT6x/feuut3H777Vx33XXs2rWLwsJCSkpKKC0tpbi4mNmzZ3POOedQVFREbW0tZ555JsOHD2/y9ZON6adiwYIFXHvttVRVVVFdXc2iRYsYPXo03/3ud/npT3/K4MGD6d+/f4N9jjnmGObPn8/rr79Oz549OeSQQ9ixY0erjgtw7bXXMm/ePGpra5kyZQpHHHEEn332GaWlpfz2t7/lkksu4corr+TRRx9lwIABLFmyhMLCQiKRCGeccQYFBQVcccUVzZ7raItAOqd2pUM4HI5qPf2Oo7paIRAgAEQBcujfTTgcJhQK5eSSAtA5lzvIplTqSvbvo6ysrGzcuHHjm9ilnoZ3RFKgC5Skq0hpeMfM3gR2xx5+CPwOuAOoAdY65641s27AMuBIoAq4wDn3DzOb0Lhtmt+DSMbEwz6SOGMnGNS0Tem0Wgx9MysCcM5NS9i2CTgD+AB4xszGAiGgyDl3TCzobwVmAPc2buuc+2ua34dIRkSSTM8MAeWatimdVCo9/SOBnma2Nta+FCh0zr0PYGbPA8cDBwHPATjnNprZeDPr00Rbhb50Wsk+CEQ6i1RC/0tgCfAA8D+ANcCuhOe/AIYDffh6CAigNrbt8yRtm1RVVdWuGxtUVlbm5I0RVFfr5Gpd9WvWb9uWM/Xt27ePaDSa9ot40iVXa+vMdbXnBjCphP4W4B/OuSiwxcx2A4lznHrjfQj0jH0f1w0v8HsnadukwsLCds3ayMlZH6iu1sq1uj4ECAZ5Yc0ar65AIGfqC4fDBAKBr2d8hELpvWo4GGzX+YvOPEsmG1Kpq3v37slm76T0+qmE/nnAaOAXZnYwXrjvNbNv4o3TnwxcC5QAPwCeiI3pv+Oc+9zMqpO0FelUQuAFX4707psViaR3SmkKq2xefvnljBgxgmg0Sk1NDbNnz2b69OlUVFTwwx/+kMMPP7zBPitWrGDfvn2Ulpby6aefEggEKC4uprS0NO1XoCZ65ZVXePbZZ1m8eDG/+tWvuOeee9LyunV1dZSWluKco6CggOuuu26/6xueeOIJHnvsMfLz87nooos47rjj6p9bsWJF/UVrmZZK6C8HVpjZq3hTlM8D6oBHgDy8GTmvm9l/Aiea2Z+BAHBubP8LG7dN83sQkSybMGFC/dLKe/fuZdasWQwbNozevXszfPhwVq5cud8+jz32GAceeCCLFy8GvOC7++67WbhwYYfUvHTp0rS91rp166iurubxxx9n06ZNLF68uMEHyvbt21m5ciVPPvkkVVVVnH322UyaNIm6ujoWLlzI22+/zUknnZS2eprTYug756qBs5M8NaFRuzq8gG+8/8bGbUWk6+rVqxdnnXUWzz33XLPLFQ8ZMoTVq1czduxYvv3tbzNr1qz6+wDcfPPNfO9732PMmDH17e+8804ikQj//Oc/2b17N2effTZr167lww8/5KabbuKoo45i5cqV/OlPfyIQCDB9+nRmz57N+++/z1VXXUWPHj3o0aMHffv2BeD444/nz3/+M2+88QZ33XUX4J1Luummm+jevTvz5s1j8ODBfPTRR4wePZprr216kKKsrIxjjz0WgKOOOop33323wfNvv/023/rWtygoKKCgoIChQ4eyefNmgsEgP/rRj5g4ceJ+y0FkipZhEGmDcmJz+DVXP6kBAwbUL0j2wQcfMGvWrPrnDj/8cK688kqmTZtGdXU1q1evZsGCBYwcOZKFCxdiZvz6179O+rpFRUUsX76c++67j5dffpl7772XJ598kmeeeYbi4mKeffZZ/vCHPxAIBJgzZw6TJ0/mjjvu4LLLLmPSpEncd999+4Xre++9xy233MKgQYO49957ee655/jBD35AeXk5y5cvp0ePHpxwwgls376dgQMHJq2r8aqaeXl51NTUkJ+fX/984oJr8VU1+/bty+TJk3nqqafa9oNuA4W+SAuCweB+J0aHAVFN3WzStm3b6le/bGp458033+SYY47hpJNOora2lj/+8Y8sWLCg2QCMr3LZu3dvRowYAUDfvn3rV9Tctm0bc+bMAWD37t1s3bqV9957r/4vhrFjx+4X+oMGDeL666+nZ8+efPLJJ4wdOxaAoUOH1gf5wIEDqaqqqt9n7969XHihN7AxceJEiouLG6yqWVdXVx/4wH7PJ1t1s6Mo9EVaUF5e3uLJTPnanj17WLVqFXfccUez7Z555hl69erF3LlzycvLw8zqbx7SlOYWkhs+fDgjRozggQceIBAIsGLFCkaOHMnw4cN58803mTJlyn7DLgALFy5k3bp1FBcXc8UVV9QPMTV3rF69ejX4IHv++ed58cUXmT59Ops2bdpvueYxY8Zw++231y/89v777+/XpqMo9EW6mmAwvR9SKayyuXHjRmbNmkW3bt2ora3l0ksvZfjw4VRUVOw3vANwww03cPnll/Ob3/yGGTNm0KNHD3r27Mn1118PJB/Tb8mhhx7KMcccw09+8hOqq6sZM2YMgwYN4pprrmHu3LksX76c/v3773frwRkzZvDjH/+YPn36cOCBB/Lpp5+mfMy4E088kQ0bNjBz5kyi0Sg33HADAA899BBDhw7l+OOPZ9asWZx99tlEo1Hmzp2b1lsgtoZW2ewgqqt1cqWuQCDg9fwCAYhG6+sKBAI5s9qmVtlsm85cl1bZFBGRlCj0RUR8RKEv0gXk2jCtZE57f9cKfZFOrqioiF27din4fSAajbJz506Kiora/BqavSPSklAopRks2VJSUsLmzZvrbzSea/bt29emm55nWmetq6ioiJKSkja/vkJfpCXpXsAszbp3705+fn5OzHZKJldmYjXm17o0vCPSjPiSyiJdhUJfpBkh0Po60qUo9EVEfEShLyLiIwp9EREfUeiLiPiIQl9ExEcU+iIiPqLQFxHxEYW+iIiPKPRFmhAKhbJdgkjaKfRFmhDRjc+lC1Loi4j4iEJfRMRHFPoiIj6i0BcR8RGFvkgTtJa+dEUKfZEmhEBr6UuXo9AXaYOg/gKQTkqhL9IG5foLQDqplG6MbmbfAMqAE4EaYAUQBd4FLnbO1ZnZNcCpsecvd869YWYjkrVN95sQEZHUtNjTN7PuwO+Ar2KblgILnXPHAgFghpmNBaYC3wFmAnc31Ta95YuISGukMryzBLgX2BZ7PA54Ofb9GuAEYDKw1jkXdc5tBfLNbGATbUVEJEuaHd4xsznAdufc82a2ILY54JyLxr7/AugL9AF2Juwa356sbbOqqqoIh8Opv4NGKisr27V/pqiu1smFukbBfjUk1pXs+WzJhZ9XU3K1Nr/W1dKY/nlA1MxOAI4C/h34RsLzvYFdwOex7xtvr0uyrVmFhYWMGjWq5cqbEA6H27V/pqiu1smVuhrX0LiuXKgRcufnlUyu1tbV6iorK0upXbPDO865Kc65qc65acAmYDawxsymxZqcAqwHNgAnm1k3MxsKdHPO7QDeTNJWRESyJKXZO43MA+43swIgDKx2ztWa2XrgNbwPkoubapuGmkVEpI1SDv1Ybz9uapLnS4HSRtu2JGsrIiLZoYuzRER8RKEv0kblALqlonQyCn2RNhoGoFsqSiej0BcR8RGFvoiIjyj0RUR8RKEv0kZaU186I4W+SBtpTX3pjBT6IiI+otAXSSKk+ffSRSn0RZKIpDj/vhx0gZZ0Kgp9kXbQBVrS2Sj0RUR8RKEvIuIjCn0RER9R6IuI+IhCX0TERxT6IiI+otAXEfERhb6IiI8o9EVEfEShLyLiIwp9EREfUeiLiPiIQl9ExEcU+iIiPqLQFxHxEYW+iIiPKPRF2iEYDGa7BJFWUeiLtEN5eXm2SxBpFYW+iIiP5LfUwMzygPsBA2qBc4EAsAKIAu8CFzvn6szsGuBUoAa43Dn3hpmNSNY2/W9FRERakkpP/wcAzrlJwNXA0tjXQufcsXgfADPMbCwwFfgOMBO4O7b/fm3T+g5ERCRlLYa+c+5/Az+PPQwCnwDjgJdj29YAJwCTgbXOuahzbiuQb2YDm2grIiJZ0OLwDoBzrsbMHgZOA/4F+L5zLhp7+gugL9AH2JmwW3x7IEnbJlVVVREOh1N/B41UVla2a/9MUV2tkyt1Na4hWV09gIOHDOH9des6sLKGcuXnlUyu1ubXulIKfQDn3DlmdgXwOt7/53G9gV3A57HvG2+vS7KtSYWFhYwaNSrVsvYTDofbtX+mqK7WyXZdHwIEg/vVkKyuABDdti2r9Wb759WcXK2tq9VVVlaWUrsWh3fMbJaZLYg9/BIvxP9iZtNi204B1gMbgJPNrJuZDQW6Oed2AG8maSuS00IAmo4pXVAqPf2ngIfM7BWgO3A5EAbuN7OC2PernXO1ZrYeeA3vw+Ti2P7zGrdN83sQEZEUtRj6zrm9wI+TPDU1SdtSoLTRti3J2oqISMfTxVkiIj6i0BcR8RGFvoiIjyj0RUR8RKEv0kgoFMp2CSIZo9AXaSQSiWS7BJGMUeiLiPiIQl9ExEcU+iIiPqLQFxHxEYW+iIiPKPRFRHxEoS8i4iMKfZF0CAZBF3VJJ5DynbNEpBnl5RAIZLsKkRappy8i4iMKfRERH1Hoi4j4iEJfRMRHFPoiIj6i0BcR8RGFvoiIjyj0RdJFF2hJJ6CLs0TSRRdoSSegnr6IiI8o9EVEfEShLyLiIwp9kQShUIhgMJjtMkQyRidyRRJEIhGi0ahOyEqXpZ6+SGOhkDf9MkXBYJCQpmpKJ6GevkhjkQhEoyk3Ly8vJ6C/DKSTaDb0zaw78CAQAgqB64C/AyuAKPAucLFzrs7MrgFOBWqAy51zb5jZiGRtM/JORESkRS0N7/wM2OmcOxY4BbgLWAosjG0LADPMbCwwFfgOMBO4O7b/fm3T/xZERCRVLYX+KmBRwuMaYBzwcuzxGuAEYDKw1jkXdc5tBfLNbGATbUW6Li3FIDmu2eEd59weADPrDawGFgJLnHPxAc8vgL5AH2Bnwq7x7YEkbZtVVVVFOBxuzXtooLKysl37Z4rqap1s19XUsZurKxwOw5o1jDrssA6vPds/r+bkam1+ravFE7lmdgjwNLDMOfcHM7s54enewC7g89j3jbfXJdnWrMLCQkaNGpVC6cmFw+F27Z8pqqt1sl1XU8durq7E7R1de7Z/Xs3J1dq6Wl1lZWUptWt2eMfMBgFrgSuccw/GNr9pZtNi358CrAc2ACebWTczGwp0c87taKKtSE5qz4VZmrYpnUVLPf2rgAOARWYWH9v/JfBbMysAwsBq51ytma0HXsP7ILk41nYecH9i23S/AZF0ac+FWZq2KZ1FS2P6v8QL+camJmlbCpQ22rYlWVsREckOXZErIuIjCn0RER9R6IuI+IhCX0TERxT6IiI+otAXEfERhb4IumOW+IfW0xch4cKsdIgvulZenp7XE0kjhb5IupWX63aLkrM0vCMi4iMKfRERH1Hoi4j4iEJfRMRHFPoiIj6i0BcR8RGFvoiIjyj0RUR8RKEvIuIjCn3xPa27I36i0Bffi0QilKdhnZxgMEgoFGr364hkktbeEUmT8vJyAlpzR3KcevoiIj6i0BcR8RGFvkgmxNfUF8kxGtMXSRQKeYHdXlpTX3KUQl8kUSQC6bqDlkgO0vCOiIiPKPRF4tI1tCOSwzS8IxKnoR3xAfX0RUR8RKEvvpbudXe0FIPkupSGd8zsO8BNzrlpZjYCWAFEgXeBi51zdWZ2DXAqUANc7px7o6m26X8bIm0TiUSIpnFIR0sxSK5rsadvZr8GHgCKYpuWAgudc8cCAWCGmY0FpgLfAWYCdzfVNr3li4hIa6QyvPM+cHrC43HAy7Hv1wAnAJOBtc65qHNuK5BvZgObaCsiIlnSYug7554E9iVsCjjn4n8PfwH0BfoAuxPaxLcnayuSezIxXVNLMUgOasuUzcQx+d7ALuDz2PeNtydr26yqqirC4XAbyvJUVla2a/9MUV2t05F1hcNhRkUihP/+d2jhmKnWFQ6HYc0aRh12WIe8j1z9PULu1ubXutoS+m+a2TTn3EvAKcCLwD+Am81sCVACdHPO7TCzZG2bVVhYyKhRo9pQliccDrdr/0xRXa3TkXXFj5PK8VKpKxgMcsopp9TfmKUj3keu/h4hd2vranWVlZWl1K4toT8PuN/MCoAwsNo5V2tm64HX8IaMLm6qbRuOJ9KpaAaP5LKUQt85Vw5MiH2/BW+mTuM2pUBpo21J24qISHbo4izxLd0QXfxIa++Ib6X7wiyRzkA9fRERH1Hoi4j4iEJfRMRHFPoimbx5iq7KlRyjE7kimbx5im6QLjlGPX2RDNC6+pKr1NMXX8r0HH1dlSu5SqEvvqQ5+uJXCn3xLw2/iA8p9MW/IpGMvnx8XL88o0cRaR2FvkiG1I/rx6dtxpZaFskmzd4RybR42Gs4SXKAevoiHUHz9SVHqKcvvhIKhQgEAlpSWXxLPX3xlQZTNePj7SI+otAXf9PJVfEZDe+IiPiIQl9ExEcU+iIZpIXXJNco9MWfMrmGfoLy2DmDUPx4+gCQLFPoi698CF7wRiIddhK3Pvgh40s/iLREs3ekS4kPpZQ3EeghyErwaqllyRUKfelSIo0CPfFDoCI/H/LyKCkpyUJlIrlBoS9dWuKHQEltbeZuiyjSSWhMX7qUD2Nf0PDuWBX5+VTk5WWtLpFcoZ6+dCmh2H8DgQAf5eVRUltLRey/6uWLKPSlk2tw4jZ2w5JQMEgUvBO20SglmiYpUk/DO9I5hULeiVkSxu0jEaYFg19PxYzPwy8v1xo7IjEKfel0QrF59iW1tZSD16sPBKjIy/t6qmauBr0u0JIsU+hL7giFmgzEUChUP5TzUiTihWe8Jx+NesM4NTUdUma76C5akmUZH9M3s27AMuBIoAq4wDn3j0wfV3JHRX4+NbW1TAsGG1w0VT8eDxCJeOPxse3rKyqoAEpqaqjIz6e8ttZrGwh4J2tzsRefKt1FS7KoI3r6PwKKnHPHAFcCt3bAMSWTQiFvOCV//z5DRX7+19tj7cA7uVoeiUAgQHn8KxKp/woFg4SiUYi3i4uHY6w3H0rs4XcyDRZf0zCPZElHhP5k4DkA59xGYHwHHFNSEQpRPHp0g1UgE4dR4iry873wjX3FQxposL1B7zWhXUlNjde7TQjuUDD49eNotOFYfDTKnnfe8fZrPGyTq2P1KUhcfC0E9R+C9V/6EJAOEIhmeO6ymT0APOmcWxN7vBUY7pxLOgBbVla2HdCqVCIirRMcN27cwJYadcQ8/c+B3gmPuzUV+ACpFC0iIm3TEcM7G4DpAGY2AXinA44pIiJJdERP/2ngRDP7MxAAzu2AY4qISBIZH9MXEZHcoYuzRER8RKEvIuIjXWqVTTPrCzwG9AKqgZ855/47u1WBmeUBS/GuUSgESp1zf8puVQ2Z2aHA68Ag51xlDtTTF/g90AcoAH7lnHsti/Xk5JXlZtYdeBDvYuZC4Drn3H9ktagEZvYNoAw40Tm3Odv1AJjZAuCHeP9fLXPOLc9ySfHf48N4v8da4F8z9fPqaj39OcA7zrkpwOPAv2W3nHqzgO7OuUnADGBElutpwMz64F0pXZXtWhL8CnjBOTcV7/d6d3bLydkry38G7HTOHQucAtyV5XrqxYLsd8BX2a4lzsymAROBScBU4JCsFvS16UC+c24i8L+A6zN1oK4W+u/w9TUBfYB9Wawl0clAhZk9A9wP/J8s11PPzALAfcBVwJdZLifRbXiBAd5fpNn+6yNXryxfBSxKeJxLq84tAe4FtmW7kAQn4+XE03j/DnPlL+4tQH7sL8qMZlenHd4xs/OBuY02XwycZGZ/B/oDx+ZIXdvxQuv7wBTgodh/O1QTtUWAx5xzb5lZR5cENFnXuc65/zSzwXjDPJd3fGUN9AF2JzyuNbP85i407AjOuT0AZtYbWA0szGY9cWY2B9junHs+NpySKw4Egnj/FocB/2Fmhzrnsj2NcQ/e0M5mvBq/n6kDdakpm2b2FPC8c+53ZjYG+L1zbkwO1PUYsMo592Ts8X875wZnuSwAzOwfQEXs4QTgjdjwWNaZ2Wi8czTz48t4ZLGWpcBG59wTsccVzrmSbNYUZ2aH4PVclznnHsx2PQBm9grerQ6iwFF4PdkfZvscm5ktxvswujX2+C288w2fZrmupUCVc25B7Pf5/4DRmTi/1ml7+k34J1/3xj7F653lglfxxuyeNLMjga1Zrqeec67+/IKZlQMnZa2YBGZ2GN7QxVnOubeyXQ/eleU/AJ7IpSvLzWwQsBa4xDn3QrbriUvsOJjZS8CF2Q78mFeBX8ZC9iC8SR87s1sS4GVXfEjnM6A7kJeJA3W10F8EPGBmv8D7of1rluuJux+4x8w24l2VfGGW6+kMbgSKgDtiw067nXMzslhPrl5ZfhVwALDIzOJj+6c453Lm5Gkucc79ycymAG/gndO82DlXm+WywDuH9aCZrcebVXSVc25vJg7UpYZ3RESkeV1t9o6IiDRDoS8i4iMKfRERH1Hoi4j4iEJfRMRHutqUTfGR2DoqLwIznXOPJ2x/G/irc25OCq9RBGx2zoWaOcaFzrmZKdZUDfw59rAH8DzeAnt1sesgtgJ1CbvMc86VmdmVwAmx56J4U/bKUjmmSGso9KWz2wz8BG+BvfhVvL2yWM9nzrlpsVoCeGvR7NVFAAACEklEQVTPXAzcGXv+pMZXWcYuRPshMMk5FzWzo/BWXDyyw6oW31DoS2f3FjDSzPo553bhrTr5CDAUwMx+irduTxXwHvBzvCWIH8G7qKl+eeTYB8Zv8S6+2gmc157CYgF+K97Sx3c20/TTWL3nmdlzzrlNZvbt9hxbpCka05eu4CngtFjP+tvEhlfMbABwLfBd59xkYBfwP/GWan43tlTA7xJe5368KzSnAc8Cv05DbZ/gLaAVt9bMXop9vQDgnNtBrKcPvGZmm8ngglvib+rpS1fwB+Ae4ANgfcL24cDfnHNfxB6/wtdrC8WXSX7dzOJrnowClsWWfeiOt0jYfszsEuBfYg9/6pz7uJnagny9oB0kH94ZAXzunDsv9ng88KyZveic+6yZ1xZpNfX0pdNzzn2AN45/Gd4yzHEfAoeZWXyMfypekG8GjgEws2/hBTyAA2bHevq/Bp5p4nh3Oeemxb6aDPzY2ujz8VYKbc4YvLWZimKPt+AtHJgLa8JIF6OevnQVjwOznHNbzGw4eMMmZnYN8KKZ1eGN31+Jd6ORh8zsVbwPgPgdwy4C/j12e0uA84GDW1lH/9iqknV4Hyb/F0i8Hd/aWC1xdzjnnjKzUcDrZrYHrzP2b865xPX7RdJCC66JiPiIhndERHxEoS8i4iMKfRERH1Hoi4j4iEJfRMRHFPoiIj6i0BcR8RGFvoiIj/x/OgOwo5DY4MYAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, ax = plt.subplots(2, 2, sharey=True, sharex=True)\ndf.hist('diff_des', by='BAND', bins=np.arange(-6., 2., 0.1), histtype='step', alpha=1, color='k', normed=1, \n        cumulative=0, lw=1, sharey=True, ax=ax)\nfor axx, band in zip(ax.flatten(), df.BAND.unique()):\n    xx = df.dropna().query('BAND == @band').diff_des\n    axx.axvline(np.median(xx), ls='dashed', color='k')\n    axx.text(-5, 1.5, 'median = {0:0.2f}'.format(np.median(xx)))\n    axx.set_xlabel('model - DES')",
   "execution_count": 18,
   "outputs": [
    {
     "output_type": "stream",
     "text": "/Users/rbiswas/soft/lsst_stack/python/miniconda3-4.5.4/envs/lsst-scipipe-fcd27eb/lib/python3.6/site-packages/pandas/plotting/_core.py:2396: UserWarning: When passing multiple axes, sharex and sharey are ignored. These settings must be specified when creating axes\n  yrot=yrot, **kwds)\n/Users/rbiswas/soft/lsst_stack/python/miniconda3-4.5.4/envs/lsst-scipipe-fcd27eb/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n/Users/rbiswas/soft/lsst_stack/python/miniconda3-4.5.4/envs/lsst-scipipe-fcd27eb/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n/Users/rbiswas/soft/lsst_stack/python/miniconda3-4.5.4/envs/lsst-scipipe-fcd27eb/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n/Users/rbiswas/soft/lsst_stack/python/miniconda3-4.5.4/envs/lsst-scipipe-fcd27eb/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n",
     "name": "stderr"
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 4 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAETCAYAAAA/NdFSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X14VPWZ//F3UAiEp1KRbaiYWVnzJayLliBStWK9IlYri0UrCouLYH2ou3uhIloXt27Xx66KSsUKlQXZiqJiFQXtgigrCP0RxaIMt6ImRHEVRAQSCQ+Z3x8zkw5kJplkJnPOzHxe15WLzJyZMx8m5557vuexIBQKISIi+aeD1wFERMQbagAiInlKDUBEJE+pAYiI5Ck1ABGRPKUGICKSp9QARMT3nHNDnHPPeJ0j1xToOAARkfx0pNcBpCnn3M3AJGA3sBK4wMwCnoYS8ZBz7kzgN2Z2gtdZcolWAfmMc+4cYAJwMlAOdPc0kIjkLDUA/zkPeNrMdppZCHjY60AikpvUAPznAFAQc/ugV0FEJLepAfjPS8CFzrmekduTAG2pF5G0UwPwGTN7FZgNvOmcWwf0BOq8TSUiuUh7AfmMc24IcMDMBkZuXw909jaViLfM7DVAewClmRqA/7wP3OScu5Lwqp8twJXeRhKRXKQDwURE8pS2AYiI5Ck1ABGRPJXxbQDr168PFRYWxp1WX19PommZ5JccoCzJ5qirq9teXl5+tEeR2iQbagH8k8UvOcD/WZKth4w3gMLCQsrKyuJOCwaDCadlkl9ygLdZ6urCe58WFRV5niVWvByVlZXVHsVps2yoBfBPFq9zxNaD11lipVIP2gtIEjrvvPMAeO2117wNIuIDuVgP2gYgIpKnNAKQhNasWUN9fT2BQICqqiqv44hImmkEIAnV19czfPhwqquzbvW6iCRBDUBEJE9pFZA0a8KECbz++utexxDx3IQJE7yOkHYaAUizcnGhF2mLCRMm5Fw9qAFIs7Zv3+51BBFf2L59e87Vg1YBSbMuuugiryOI+EK0FnQcgIiIZD2NAKRFJSUlFBQU0LdvXz799FOv44hImqgBSIuiB4EVFBQ0/0ARySpaBRTHJ598wvjx4wG47rrr2LdvX8Zee8eOHUycOJGxY8cyefJk6uvr4z7unXfeacwI8N577/GDH/yA8ePHM378eJYsWZKpyJLj/FIP//mf/8k333wT93GH10MwGGTs2LGMHz+eSZMm5dzG23TRCKAF06dPz+jrzZw5k/PPP5/Ro0cza9YsXnnlFU466aRDHjN79mxeeOEFunTp0njfxo0bufzyy5k4cWJa81xzzTVpnZ9kNy/r4fbbb+epp55qsitmvHq44447uPXWWykrK+PJJ59k9uzZ/OIXv0gpSy7WQs41gEWLFrFixQr27t3Ltm3buOyyy1i+fDkffPABU6dOpaKigqVLlzJ37lw6dOhAeXk5U6ZM4YsvvmDKlCmEQqFDzq191llnsXTpUqqrq7n77rtpaGhg165dTJs2jcGDBzNixAgGDx7Mxx9/zFFHHcWMGTM44ogjGp8/ffp03nrrrUMyPvbYY3Tq1Clu/srKSq666ioAzjjjDH71q181ecyxxx7LjBkzmDp1auN97777Lh9//DHLly+npKSEW265hW7duqX0XgKMGTMm5XmId3KpHgYPHswf/vCHJg0gXj3cf//99OnTB4CDBw+m5dz9uVgLSTUA59wpwD1mduZh918PTAK2Re66yswsrQnboLa2ljlz5vDSSy8xd+5cFi5cyNq1a3n88ccZMmQIM2bM4Nlnn6VLly7ceOONrFq1ijfeeIPzzz+fiy++mNmzZ7Ny5cpD5rl582ZuuukmnHMsXryYRYsWMXjwYGpqapg3bx7FxcVccsklbNiw4ZBv7Nddd12rsu/Zs4fu3bsD0LVr18ZzkMc655xz+OSTTw65b9CgQfz0pz/lhBNO4JFHHuHhhx/mpptuatVrx1NTU0O/fv1Sno94J1fqoUuXLuzevbvJY+LVQ/TD/6233uK///u/+f3vf9+q142npqYGIKfqocUG4JybCowHauNMHgxcZmaV6Q6WiujFEbp3707//v0pKCigZ8+e1NfXs2XLFnbs2MGVV14JhIujpqaGDz74gFGjRgEwYMCAJgt8nz59mDlzJp07d6a2trbx23WvXr0oLi4GoLi4uMk6+5a+8axbt44HH3wQgEmTJtGtWzdqa2sbX6dr165J/Z/PPvtsevTo0fj7f/zHfyT3ZrVg/PjxObXfcz7KlXr45ptvGpfxZCxZsoRHHnmEWbNm8e1vfzvp5yUS3caQS/WQzAjgQ2A0MD/OtHLgF8657wAvmdld6QzXVs3trXLMMcdQXFzMnDlz6NixI4sWLaKsrIyPPvqIt99+mwEDBrB58+Ymz7vjjju499576d+/Pw899FDj7pAt7RnT0jeeIUOGMH/+X97a//3f/+X1119n9OjRrFy5koEDBzb7/KhJkyZx6623MmjQIN58803+9m//NqnnSe7LlXp46623KC8vb/b5Uc8//zxPPfUU8+fP51vf+lZSz8lHLTYAM3vWORdIMPlJ4GFgF/Ccc+58M3uxufnV19cTDAbjTtu7d2/CacnaunUrX375JcFgkJqaGnbu3EkwGOSjjz6itraWzz//nBEjRnDhhRfS0NBAnz59OO644/jRj37Evffey9NPP03v3r2pra0lGAyyb98+Nm3axCmnnMIVV1xBz5496d27N7t27SIYDHLgwIHGzLt27aK6urpV31IOV1FRwYMPPsi8efPo0aMHP//5zwkGgzzwwAOMGzeOo48OX+bz888/55tvvml87X/8x39k2rRpHHnkkfTq1avxeamqq6s7ZD7pmGeq0rGc+EF71wLkVj107dqVH//4xy3Ww8GDB/nVr35F7969G3eKOOGEE7j00ktTei+jq2ODwaCvlsGUsoRCoRZ/SktLA6WlpWsOu6+gtLS0Z8ztn5eWlt7a0rw2btwYSqS5aZnklxyhkLdZgNDw4cMPue0H8d6TdevWrQslsSz76ScbaiEU8k8Wr3MMHz68sR68zhIrlXpIZS+gHsC7zrkywtsHzgLmpDA/ERHJoFY3AOfcWKCbmc1yzt0CrADqgeVmpqOPckAgEKC6upqjjz6aG264wes4Ir6Qi7WQVAMwsypgWOT3J2Lun0/8jcOSxaqrqwmFQl7HEPGVkSNHeh0h7XQqCEnIzPDBYR0ivpCL9ZBzRwJL+kSPwMyl/Z5F2ioX60EjABGRPKUGICKSp7QKSEQkCWvWrPE6QtqpAYiIJCHRtTmymRqAJDRt2jSvI4j4Sp8+fXLq8qhqAJJQRUWF1xFEfOXzzz8H4Lvf/S4FBQWUlJQ0XjI1G2kjsCS0fv161q9f73UMEd+I1sOyZcsIhUJUV1d7nCg1agCS0OTJk5k8ebLXMUR8I9fqQQ1ARCRPqQGIiOQpNQARkWYEAgEKCgrScmF5v1EDEBFpRvTsuMOGDfM6StppN1BJ6M477/Q6gohv5GI9JNUAnHOnAPeY2ZmH3T8S+DfgADDHzGanPaF45tRTT/U6gohv5GI9tLgKyDk3Ffgd0Pmw+zsC04ERwHDgSufcd9ojpHhj9erVrF692usYIr6Qi/WQzAjgQ2A0Ta/8VQZsNrOvAJxzbwA/AJ5ubmb19fUJr2Cf0tXt08gvOcC7LMFgsHGf53nz5h1yv9f89PdJRTbUAvgni5c5Dq+H2CxevzepvC8tNgAze9Y5F4gzqQfwdczt3UDPluZXWFhIWVlZ3GnBYDDhtEzySw7wLktZWRlFRUWNvwP07duXc8891/ND3+O9J5WVlR6labtsqAXwTxYvcxxeD7FZvH5vUqmHVPYC2gV0j7ndHdiZwvzE55YtW5b1h76LyF+k0gCCwPHOuW875zoBZwBvpieWiIi3ovv/l5SUeB2l3bR6N1Dn3Figm5nNcs5dD7xCuJHMMbPsPz+qiAh/2f8/lyXVAMysChgW+f2JmPsXA4vbJZl47oEHHvA6gohv5GI96EAwOUQgEGgc8p500kkepxHxj1ysBzUAOUTssHfZsmWALgwjArlZD2oAktDtt98O5NYCL9JWuVgPOhmciEgblZSUEAgEvI7RZhoBiIi0UVVVFQUFBV7HaDONAERE8pQagIhIntIqIEno0Ucf9TqCiG/kYj2oAUhCzjmvI4j4Ri7Wg1YBSUKLFy9m8WId6C0CuVkPGgFIQvfddx8AI0eO9DiJiPcS1UNJSUnjSeO8PlV6a6kBiIikIPqhn427g2oVkIhInlIDEBHJUy2uAnLOdQBmAicC9cAVZrY5ZvpDwGmELwkJMMrMvm4yI/G1QCBAdXV1Tl/8QkQOlcw2gAuAzmb2fefcMOA+YFTM9MHAOWa2vT0CSmbEu/jF/Pnzmzwumzd4iaQiXj1ku4KWrnjjnLsf+JOZPRm5/amZfTfyewfgM2AV8FfAY2Y2p7n5rV+/PlRYWBh32t69e+ncuXOr/xPp5pcckLksAwcOZOPGjUlnSebx7SXee1JXV1dZXl4+xJNAbZQNtQD+yZLpHM0t4/GyeFUTqdRDMiOAHkDsKp2DzrkjzewA0BWYAdwPHAGscM6tM7M/J5pZYWFhkyvYR8W7ur0X/JIDMpvl8Nd56qmnABgzZkzcLF69R/Hek8rKSk+ypCIbagH8k8WLHLGvF1sPibJ48T6lUg/JbATeBXSPfU7kwx+gDnjQzOrMbDfwKuFtBZIDHnnkER555BGvY4hkVKKLwediPSQzAlgFjAQWRrYBbIiZVgo86ZwbTLiZnA7MS3tKEZEMiJ7bP9cvBh+VTAN4DjjbObcaKAAud85dD2w2sxecc78H1gD7gcfN7L32iysi0n7i7QyRy1psAGbWAFx92N2bYqb/Gvh1mnOJz2lvIJHspwPBpE2qqqoIhUJUV1d7HUXEF6JfirLpEpE6F5Ak9Mwzz3gdQcQ3WqqHbDwnkBqAJNS7d2+vI4j4Ri7Wg1YBSUJz585l7ty5XscQ8YVcrAc1gDyXaJ9nyM0FXqStcrEetAooz+Xbbm8i8hcaAYiI5Ck1ABHJe82tCs1lagCSkui+z9m2/7NIrOiq0HQc1JhNxwNoG4AktGTJkhYfE1sw2bT/s0hrJVMPkF3HA6gBSEJFRUVeRxDxjVysB60CkoRmzpzJzJkzk358Ng19RVqrLfXg91pQA8hjgUCg2Y1eCxcuZOHChUnPT+cHkmwR3egb/Ulm429b6sHvtaAGkGOiC3a8bx6HL/SAzuQpeSm60Tf601514PdRcYvbACLX/Z1J+Epf9cAVZrY5ZvrPgKuAA8DtZvZiO2WVwwQCAaqrqw85JXN0wY5+wEcfA+GFUQd9ST6JfvBWVVUdUi+Z4vcNwsmMAC4AOpvZ94GbgfuiE5xz3wH+BTgNOAe4yzkX/yrXkrRAIJDUN4bYo3gPH8pGv3kA7f4tR8QLsaPdw0e3saNc4JBaAPJuf/9EkmkApwMvA5jZGiD2SvNDgVVmVm9mXwObgUFpT+lDiRa4dPxEFRQUUFFRkfD1ogtxdN177Id89D4vPvTjDXsPz+/XIbG0XXN/43TUy8CBA5vUSeyINnaVTmw9HF4L0fsyya+rggpaWiXgnPsd8KyZLY3c3gIcZ2YHnHP/APydmd0UmfY44ctCLks0v8rKym2Av7eMSDYqKS8vP9rrEK2hWpB2lFQ9JHMcwC6ge8ztDmZ2IMG07sDO5maWbUUq0l5UC+K1ZFYBrQLOA3DODQM2xEz7E/AD51xn51xPoAx4N+0pRUQk7ZJZBRTdC2gQUABcTrghbDazFyJ7AV1JuJncaWbPtm9kERFJhxYbgIiI5CadCyjDnHOdgU1mFkgw/UzgajO7JMn57QNWR252AV4BbjOzBudcFbAFaIh5yg1mVumcuxmoiEwLAbeYWWWr/0MibaRa8J4aQPbbYWZnAjjnCoDfAtcCMyLTR5jZ3tgnOOcGAn8PnGZmIefcScA8wgf7iWQr1UIrqQG0gnNuAjCS8LeLYuBBYBRwAjDFzJ53zo0DJhM+avoDwttHCoHfA70IHysRnd/fAQ8R3rbyJTAxlXyRBfg+YA5/Wejj+QI4FpjonHvZzNY754am8tqSX1QLuUHnAmq97mZ2HnAPcA0wmvCCfblz7ijg34GzzOx0wrvEXgVMAN41szOAR2PmNRu4NvKtZQkwNQ35Pgd6x9z+o3PutcjPcgAz207kWw/wpnNuE3B+Gl5b8otqIctpBNB6b0f+3QkEI980vgI6A8cB75nZ7shjVgIjIr9Hj6Ze65zbH7mvDJjpnAPoCLwf7wWdc/8EXBS5Oc7MPm0mXwnwSczteMPevwF2mdnEyO0hwBLn3Aoz29HMvEViqRaynBpA6zW329THwEDnXFczqwWGE16QG4DvA887575HeAEHMOAyM9vinDuN8FC6CTP7DfCbloJFdtmdAjzZwkMHAdc450ZGCuJ94GvgYEuvIRJDtZDl1ADSyMy2O+d+CaxwzjUQXsd5M+Ezpf6Xc+4NYBPhdaIQHjY/7pw7InJ7EtC3lS/7befca4QLqyPwP8BjMdP/GMkS9aCZLXLOlQFrnXN7CK8KvDFyPieRlKkWsoOOAxARyVPaCCwikqfUAERE8pQagIhInlIDEBHJU2oAIiJ5Sg1ARCRPqQGIiOQpNQARkTylBiAikqfUAERE8pQagIhInlIDEBHJU2oAIiJ5Sg1ARCRPqQGIiOQpNQARkTylBiAikqfUAERE8pQagIhIntJF4X3KOXcm8CBQC3QDTjaz+mafJJJjnHPXAj+LuWsgcI+Z3epRpJyii8L7VKQBLAeOM7Nqj+OIeM45dw0wCRhuZrVe58kFGgH4W40+/EXAOfcTYApwmj7800cNwN/2eB1AxGvOudOAh4EKM/s/r/PkEm0EFhHfcs6VAU8DY81so9d5co1GACLiZw8AnYB7nXPRz6t1ZnaFh5lyhjYCi4jkKa0CEhHJU2oAIiJ5Sg1ARCRPqQGIiOSpjO8FtH79+lBhYWHcafX19SSalkl+yQH+zNLQ0ABAhw7efH+I957U1dVtLy8vP9qTQG2UDbUA/snilxwAe/fupVOnTp7VQKxU6iHjDaCwsJCysrK404LBYMJpmeSXHODPLGeeeSYAr732mqc5YlVWVmbdEdPZUAvgnyx+yQEwdOhQioqKPKuBWKnUg/ftS0REPKEGICKSp9QARETylBqAiEie0rmApNUmTJjgdQQRT11wwQX07dvX6xgpUwOQVlMDkHz3k5/8xDd7JKVCq4AkaRUVFQQCAbZv38727du9jiPima+++ionakAjAEna1q1bAbjooosA744DEPHa5MmTfXMcQCo0AhARyVNqACIieUoNQEQkT6kBxPHJJ58wfvx4AK677jr27dvXrq/36quvcuGFFzJmzBgWLlyY8HF33nknCxYsaLw9a9YsRo0axbhx41ixYkW7ZowqKSnh9ddfZ82aNRl5PfGe3+qhurqaSy+9lLFjx/LLX/6y8eSEixYt4qc//SmjR4/m4YcfbteMuUINoAXTp0+nU6dO7Tb//fv3c9dddzFnzhzmz5/PU089xbZt2w55zI4dO7jiiit49dVXG+8zM1588UUWLlzInDlzeOihh/jmm2/aLWdUVVUVTz75JPX19e3+WuI/fqiHu+66i8mTJ/PEE08QCoVYvnw5W7ZsYcGCBcyfP59nnnmG/fv3s3///nbLeckll3DNNde02/wzJef2Alq0aBErVqxg7969bNu2jcsuu4zly5fzwQcfMHXqVCoqKli6dClz586lQ4cOlJeXM2XKFL744gumTJlCKBQ65NSqZ511FkuXLqW6upq7776bhoYGdu3axbRp0xg8eDAjRoxg8ODBfPzxxxx11FHMmDGDI444ovH506dP56233jok42OPPdZYRB9++CHHHnssPXv2BKC8vJx169Zx7rnnNj6+traWf/7nf2blypWN93344YcMHTq0MWtJSQlmxkknnZT+N/UwY8aM4ZJLLmn315HU5WI9vPfeewwdOhSAM844g1WrVvHll19ywgkncNNNN7Ft2zauvvpqOnbs2G7v67nnnpsTxwHkXAOA8AfmnDlzeOmll5g7dy4LFy5k7dq1PP744wwZMoQZM2bw7LPP0qVLF2688UZWrVrFG2+8wfnnn8/FF1/M7NmzD/mwBdi8eTM33XQTzjkWL17MokWLGDx4MDU1NcybN4/i4mIuueQSNmzYcMiH8HXXXdds1j179tC9e/fG2127dmXPnj2HPKZfv37069fvkEzOOWbNmsWePXvYv38/b7/9NmPGjEnlbUtaTU1NRl5H0iPX6iEUClFQUNA4fffu3Xz11VesW7eOBQsWUF9fz6WXXsozzzxDjx49Un374vrss8/o1q0b/fr1a5f5Z0pSDcA5dwpwj5mdedj91wOTgOgY7Sozs7QmbINoZ+7evTv9+/enoKCAnj17Ul9fz5YtW9ixYwdXXnklEC6OmpoaPvjgA0aNGgXAgAEDmizwffr0YebMmXTu3Jna2lq6desGQK9evSguLgaguLi4yaqRlr7xdOvWjdra2sZptbW1hxRAIv3792fcuHH87Gc/o6SkhBNPPJFevXol/R6lIro+WLJDrtVD7EVYamtr6dGjB9/61rcYOnQo3bp1o1u3bvTv35+qqioGDRrUtjetBTfffHNOHAfQYgNwzk0FxgO1cSYPBi4zs8p0B0tF9NtBPMcccwzFxcXMmTOHjh07smjRIsrKyvjoo494++23GTBgAJs3b27yvDvuuIN7772X/v3789BDD/Hpp5+2+FrQ8jee/v37U11dzc6dOykqKmLdunVMmjSpxf/jjh07+Oqrr1iwYAG7d+9m4sSJHH/88S0+T/JPrtXDwIEDWbt2LaeccgorV65k2LBh9O/fnyeeeIL6+noOHjzYuCpJmpfMCOBDYDQwP860cuAXzrnvAC+Z2V0tzay+vp5gMBh32t69exNOS9bWrVv58ssvCQaD1NTUsHPnToLBIB999BG1tbV8/vnnjBgxggsvvJCGhgb69OnDcccdx49+9CPuvfdenn76aXr37k1tbS3BYJB9+/axadMmTjnlFK644gp69uxJ79692bVrF8FgkAMHDjRm3rVrF9XV1a0edo4bN45x48bR0NBARUUFO3bs4J133uGll15iwoQJjfPftm1b4+uFQiE2bNjAj3/8Y4488kjGjx/P+++/n9J7l4xgMEhdXV3j715Ix3LiB+1dC5Bb9fDCCy9w7bXXctFFF3HPPfdw4MABjjnmGI499lgaGho4/fTTueCCC4Dwydo+++wzPvvss5Tfw3gaGhqoq6vzxXKY0rISCoVa/CktLQ2UlpauiXP/L0tLS3uXlpZ2Ki0tfam0tPT8lua1cePGUCLNTcskv+QIhfyVJby4hELDhw9v/N0L8d6TdevWrQslsSz76ScbaiEU8k8Wv+QIhUKhk08+OTR8+HCvY4RCodTqoc27gTrnCoAHzGy7me0DXgK+19b5iYhIZqWyF1AP4F3nXBnh7QNnAXPSkkp87YYbbuD111/3OoaIZyZMmJD1ewBBGxqAc24s0M3MZjnnbgFWAPXAcjNbku6A4j8jR470OoKIp374wx/mz3EAZlYFDIv8/kTM/fOJv3FYcpgP9vQV8dTHH39Mhw4dcM55HSUlOXkgmLSvq666yusIIp667bbbcuI4AJ0LSEQkT6kBiIjkKTUAEZE8pQYgItJGgUCAQCDgdYw200ZgabVp06bpOADJa1dddRUlJSWcffbZXkdJiUYA0moVFRVeRxDx1KmnntpYByUlJRQUFGTlSEAjAGm19evXex1BxFPBYLDxVNdVVVVAy2dC9SM1AGm1yZMnex1BxFN33303RUVFXsdImVYBSYsCgQAFBQX07dvX6ygikkYaAUiLqqurCYVCvjj3uYikj0YAIiJ5Sg1ARCQNsnFvIK0Ckla78847Oe2007yOIeKZyZMnEwgEDqmDbNwbKKkG4Jw7BbjHzM487P6RwL8BB4A5ZjY77QnFd0499VSvI4h46sYbb2Tr1q2UlJR4HSUlLa4Ccs5NBX4HdD7s/o7AdGAEMBy4MnJxeMlxq1ev9jqCiKe2bt3KqlWrGr/1Z6tktgF8CIyOc38ZsNnMvopcE/gN4AfpDCf+dMstt3gdQcRzuVAHLa4CMrNnnXOBOJN6AF/H3N4N9GxpfvX19Ql3J9y7d68vdjX0Sw7wT5ZgMNiYpa6urvE+L/jlPUlVNtQC+CeLX3JE1dXVxc3Tt2/fxuNmli1b1u45UnlfUtkIvAvoHnO7O7CzpScVFhYmvJZmMBj0xXU2/ZID/JOlrKysMUv0CEivcsV7TyorKz3JkopsqAXwTxa/5IgqKiqKm+fTTz8FwhuDM5E3lXpIZTfQIHC8c+7bzrlOwBnAmynMT7JIYWFh1u3yJpKq6FHxnTp18jpKWrS6ATjnxjrnrjSz/cD1wCuEP/jnmNmn6Q4o/jRs2DBCoRDV1dVeRxHJmOhR8SeeeKLXUdIiqVVAZlYFDIv8/kTM/YuBxe2STHzrgQce8DqCiKduvvlmjjvuOK9jpEwHgkmrnXTSSV5HEPFUWVmZr7ZHtJVOBSGttmzZsozs3SDiV6tXr86JGtAIQFrt9ttvB3RlMMlfjz76KEVFRVlfAxoBSELRPR6y/XB3EYlPIwBJKLrHg4jkJo0ARETylBqAiEie0iogabVHH33U6wginrrtttvo37+/1zFSpgYgreac8zqCiKf++q//OifqQKuApNUWL17M4sXhA8Cz8TJ4IqlasWJFYw1kM40ApNXuu+8+AEaOHJmVl8ETSdXcuXMpKipi5MiRXkdJiUYAIiJ5Sg1AmtABYCJNBQKBnKsJrQKSJnQAmEhTra2LkpISAoGAr68b3GIDcM51AGYCJwL1wBVmtjlm+kPAaYQvCQkwysy+bjIjEZE8UlVV5fttY8mMAC4AOpvZ951zw4D7gFEx0wcD55jZ9vYIKP4zf/58ryOIeOruu+/m+OOP9zpGypLZBnA68DKAma0BhkQnREYHxwOznHOrnHMT2yWl+Eq/fv3o16+f1zFEPFNcXJyOSkOpAAAI/ElEQVQTNZDMCKAHELtK56Bz7kgzOwB0BWYA9wNHACucc+vM7M+JZlZfX5/wCvapXN0+nfySA7zLEu81o1mWLl0KwLnnntvic9qDn/4+qciGWgD/ZPFDjujrP//88yxdurRJDTT3nPaS0vsSCoWa/SktLb2/tLT04pjbn8T8fkRpaWn3mNu/Li0tHd/c/DZu3BhKpLlpmeSXHKGQN1nCi0VT0SzDhw8PDR8+PKnntId478m6devWhVpYlv32kw21EAr5J4vXOWKX8ZNPPrlJDbT0nPaSSj0kswpoFXAeQGQbwIaYaaXAG865I5xzHQmvLnqrba1IREQyKZlVQM8BZzvnVgMFwOXOueuBzWb2gnPu98AaYD/wuJm9135xRUQkXVpsAGbWAFx92N2bYqb/Gvh1mnNJlsmGfZ5F2iIQCFBdXZ1zB4GBDgSTNMmGfZ5F2iKXD4xUA5BWe+aZZ7yOIOKpBx54gNLSUq9jpEznApJW6927N7179/Y6hohnevXqlVQN+P106WoA0mpz585l7ty5XscQ8cxzzz2XVA1UVVU1rj7yYyPQKiBpteiCP2HCBE9ziHjlD3/4A0VFRUnXgF+vm6ERgBwilVPe+n24KyKH0ghADpHKHg9+/ZYjIvFpBCAikqfUACTttCpIJD6/1YZWAUmrLVmypNnpWhUkue63v/0tAwYMaPXz/FYbGgFIqxUVFVFUVOR1DJF21dy1sbt06ZITNaAGkMcCgUCbhqIzZ85k5syZLT4uen4gkWwU3SEi3vmtFixYkFQN+J0aQB6rrq6muroaaP7bzuEWLlzIwoULW3xcVVVV4/xFcsnLL7+cVA34nRpAHol+yEd/SkpKGjdKAQm/7Yjkumht5NuIVRuB80Ds6Wxz9ayGIqmIru6JNoFcPf3z4VpsAJELv88ETgTqgSvMbHPM9J8BVwEHgNvN7MV2yiqtEF2IAX3wixwm9kM+dtQb/dBv73rxy/UzklkFdAHQ2cy+D9wM3Bed4Jz7DvAvwGnAOcBdzrnC9giaa5LdAFtRUZFwaHr4Kp3YH/jL9Z69XMj8tt+zZJfWrpo5/PGJagQ45CRt0Q/+qqqqjNRL7PYxL1c/FbTU6Zxz9wN/MrMnI7c/NbPvRn7/e+A8M7s6cvs54E4z+3+J5ldZWbkN0JZBSbeS8vLyo70O0RqqBWlHSdVDMtsAegBfx9w+6Jw70swOxJm2G+jZ3MyyrUhF2otqQbyWzCqgXUD32OdEPvzjTesO7ExTNhERaUfJNIBVwHkAzrlhwIaYaX8CfuCc6+yc6wmUAe+mPaWIiKRdMtsAonsBDQIKgMsJN4TNZvZCZC+gKwk3kzvN7Nn2jSwiIunQYgMQEZHcpCOBRUTylBqAiEieUgMQEclTagDSosiOAL6iI87FC7lWCzoZnMTlnDsOuB8YAhyILPgbgOvM7P0M5hgJ/AbYD/yrmT0VmbQUOCtTOSR/5XIteNYAnHPnAX8DLAbmAqWED4u/2szWZzhLp8Pu+iNwNlBgZvsymGMY8DDwDXCzmb0Ruf85M/tJpnJE/A74hZmtPSzffxE+91Om/CvwPcK7ID/tnOtsZvMit3OCaiFuDtVCU2mvBS+HM7cBzwAzgFvNrJjwWUUf8SDLF8AWYBNgwCnA+5HbmXQfcCnh9+Eh59yIyP3fynAOCJ8AcG3sHWa2xoMc+8xsh5l9CYwC/sk590Mgl/Zfvg3VwuFUC02lvRa8XAVUb2ZbnXOY2UoAM3vHOedFlmHAvYS7/Abn3Aoz+6EHOfZHh5SRb4X/45wbizcfdu845+YALxM+31N3wgcA/jnDOaoiJyS81cx2O+dGA6/gzQdBe1EtNKVaaCrtteBlA6h0zv0GWOWcewx4kfCbujHTQcxsk3PuUmCWc+5FvPt2ucs59y/Ao2b2f5EFfiHgxQbPnxM+FfjphE/6t4vw3+i5DOeYCPwDkb+JmdVEvvX8IsM52pNqoSnVQlNprwXPjgSObEgZT/g6An9FeOj5BjA7k+sa4+R6ESg1s1IPXrsHcD3hDU7Hm1mlc24g4VNsXJDpPJIZqoW4r61ayABfnArCOfdnMxvkdQ7wTxa/5JDM8tPf3S9Z/JIjF/lun1ZplDN7uYikSLXQTtQA/Mv7oZmIP6gW2olfGoCfOrxfsvglh2SWn/7ufsnilxw5xy8NYLrXAWL4JYtfckhm+env7pcsfsmRc3yxEVhERDLPLyMAERHJMDUAEZE8pbOBZphzrjOwycwCCaafSfgkYJckOb99wOrIzS6EDw2/zcwanHNVhM/r0hDzlBsiB9XcDFREpoWAW8ysstX/IZE2Ui14Tw0g++0wszMBnHMFwG+BawmfWAxghJntjX1C5IjKvwdOM7OQc+4kYB5wYsZSi6SfaqGV1ABawTk3ARhJ+NtFMfAg4bPynQBMMbPnnXPjgMlAPfABcCXh85f8HugFbI6Z398BDxHeze1Lwuf6aLPIAnwfMIe/LPTxfAEcC0x0zr1sZuudc0NTeW3JL6qF3KBtAK3X3czOA+4BrgFGE16wL3fOHQX8O3CWmZ0O7CR8OtsJwLtmdgbwaMy8ZgPXRr61LAGmpiHf50DvmNt/dM69FvlZDmBm24l86wHedM5tAs5Pw2tLflEtZDmNAFrv7ci/O4Fg5JvGV0Bn4DjgPTPbHXnMSiB6HvOXAcxsrXNuf+S+MmBm5LS/HQmfd70J59w/ARdFbo4zs0+byVcCfBJzO96w92+AXWY2MXJ7CLAkcurfHc3MWySWaiHLqQG0XnMHTnwMDHTOdTWzWmA44QW5Afg+8Lxz7nuEF3AIX3DjMjPb4pw7jfBQugkz+w3hS8E1K3JWySnAky08dBBwjXNuZKQg3id8nvODLb2GSAzVQpZTA0gjM9vunPslsMI510B4HefNwAHgv5xzbxC+slJ95CnXAI87546I3J4E9G3ly37bOfca4cLqCPwP8FjM9D9GskQ9aGaLnHNlwFrn3B7CqwJvNLOvW/naInGpFrKDjgQWEclT2ggsIpKn1ABERPKUGoCISJ5SAxARyVNqACIieUoNQEQkT6kBiIjkqf8P3wcFFx9M2x0AAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "So on the average, the model works quite well in each band. However, in individual measurements there are large variations including a tail which indicates that the model can be extremely bright compared to DES measurements on a handful of events. "
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.groupby('BAND').agg(dict(diff_des=lambda x:np.median(x), diff_lsst=lambda x:np.median(x)))",
   "execution_count": 19,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 19,
     "data": {
      "text/plain": "      diff_des  diff_lsst\nBAND                     \ng    -0.146185  -0.095597\ni    -0.120146  -0.108651\nr     0.063058   0.097259\nz     0.191831  -0.266670",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>diff_des</th>\n      <th>diff_lsst</th>\n    </tr>\n    <tr>\n      <th>BAND</th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>g</th>\n      <td>-0.146185</td>\n      <td>-0.095597</td>\n    </tr>\n    <tr>\n      <th>i</th>\n      <td>-0.120146</td>\n      <td>-0.108651</td>\n    </tr>\n    <tr>\n      <th>r</th>\n      <td>0.063058</td>\n      <td>0.097259</td>\n    </tr>\n    <tr>\n      <th>z</th>\n      <td>0.191831</td>\n      <td>-0.266670</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)\nfor (band, group), ax in zip(df.groupby('BAND'), ax.flatten()):\n    group.plot(x='SKYMAG', y='diff_des', kind='scatter', ax=ax, title=band, marker='.')\n    group.plot(x='SKYMAG', y='diff_lsst', kind='scatter', ax=ax, title=band, color='r', marker='.')",
   "execution_count": 20,
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 4 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAESCAYAAADnvkIDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd8FEX7wL9X0hMChJJAKIFwS+gtFAW7othAX2zYULG/FlCsr/raRQUbvir2imL5qdhFUDokdBI2BRI6hISQcrm73O3+/riSu8tecoH0zPfzuU9uZ2Znn83N7LPzPM/M6FRVRSAQCAQCf/RNLYBAIBAImidCQQgEAoFAE6EgBAKBQKCJUBACgUAg0EQoCIFAIBBoIhSEQCAQCDQRCkIgELRoJEkaJUnS100tR2tEJ+ZBCAQCgUALY1MLIGh4JEl6ELgRKAX+ASbLsty7SYUSCOoJSZJOA96QZXlQU8vS2hAmplaOJEkTgeuBVGAkENOkAgkEghaDUBCtn0nAIlmWi2VZVoH5TS2QQCBoGQgF0fqxAzqvY0dTCSIQCFoWQkG0fn4CLpUkKdZ1fCMgIhMEAkGtCAXRypFl+S9gAbBakqQ0IBYwN61UAoGgJSCimFo5kiSNAuyyLA9wHc8EwptWKoGg/pBleRkgIpgaAKEgWj9ZwAOSJN2M07S0G7i5aUUSCAQtATFRTiAQCASaCB+EQCAQCDQRCkIgEAgEmrRoH8SmTZvUsLAwrFYrYWFhTS1Os5BDyFB/MpjN5iMjR47sXE8i1Ynm1LaFDM1HhvqSI9i23aIVRFhYGCkpKWRmZpKSktLU4jQLOYQM9SdDenp6fj2JU2eaU9sWMjQfGepLjmDbtjAxCQQCgUAToSAEAoFAoIlQEAKBQCDQRCgIQYtHsTs4krsbVVGaWhSBoFUhFISgRZP113LKQ4zEJvcio+8QFLtYrFYgqC+EghC0WHb8sYx+Z55CNBACDMjbzoHtclOLJRC0GoSCELRI7BYbyeecDvhudnFwyaqmEUggaIUIBdGGqKutvjnb9nev3oCBKuWguj5Db7+66YQStHgUi5Win39nT/oWNvQdyvp57+Kw2XwL2Wzw3Xdw3XVgbt0r5wsF0UZQ7A52DBhFbL8k5L5DcNgqA5a1lVWw+aNFZPUZSIfk3mT2G07Wb/9wKCObguw8HLbKJlccvSeMojQ00qMYFMB69CjGcLGSucBJcXYeW4aMY/kpF7Hvr7/ZE9+TA1ddi+J+qNtssGwZOJx+q71rN0BEOB3On0jiqKEM37mFUTNnoIaFcciUQs8BA/hx0uUoYWGol1yC+vHHqFFRFMk1mDUtFsqfe5EDiX2wHjjQ8Ddd36iq2mI/GRkZqqqqqvtvU9Mc5Agkw96NGWolqCqoCqiZSQNVR6VdtZaa1a2ffqfays1q3up09Y/zrlDtrjKKV3nvj9312dZ7oFpZYVULcvJVxeFQHZV29c8X56uH9WHq6lfeVg/vyFUPbt2hbho0Wi07cKDe77WywqpmLP5L/eP1T1VbRUVQ/4e6kJaWlqaKtt3yZNi1Sy3rm1yt3Xp/CkF1eB3v+HxRtTavBmj//nkHQ6PU7P0F1cQ4tHFbtfOObdzWuP+LAATbtlv0UhuC4Ni7OYPuwwcCzrdtHWDalcGmT79h6PQrGIiKCvR0fcDXru9/rHcdD8jbTllEGO2BvUAH4AxXmU733OIp3wUgIYHD27PoMqBfvd2XMTyUlPNPp+kXPxA0GcXF8Pjj0L07XHstthm3ELL4ByJd2d4mSO/vHbzyVcB01VSf8u5yeJXTQgd0sZXzyFMf8e6bszzpP8+Zx3kPzKwmQ+kFF9Ju78463mTTIRREK0SxOyjK30dcUiKH5J10H+ZUDu7O4OwsKsOnX+5J9/7rT03KItp1nBigrLtjqsD+y66my7a1db8hF5biMtY88woRPy+m+wv/pft5Z6EzGHwLKQoUFDjf77RQFDh0CHQ66NrV+VfQMikuRu3Qoer4gQcIcX31f9B7P6T983V+ZfzzvfFuz95lXnrvUVYflel17+3kffkV5819TrOe5699nNdruqdmhlAQrQhLcRk75i2g4ycL6GgpIzuuB0mFe4CaG79/Bwr2kRlMh/JOtx4uYM3ES3EMGcbo/84iJDKy6oHepUuND2tLcRmhHWI41Z1w4bmYjaFElJehD3U9FhQFx4RT0K1eRY9hwyAtDfQuN5vdDsuXo0yZgu7YMQDMAwYRsTEdfWhokHcsaDYoCqUz7/O8oEDgtut+oHsfu3GfE2z7DzTCiLVZGLtwASxcQEKAsuXAcw9dUkPtzQ/hpG4FlB0q4q/rbyesQwyTF8yji6WMEKBf4R6M+HYQLeXgTUO+T48u2MWY37/lpJcewxgVxR6dkdKYWNT4eNSRI50PcX/sdti0iU233+3pyO5PpN1GwZo0UBSUffs5uHQV+lUr0asqURs3YlmxBkaMoFynwxESgnrGGeiOHas6P2Mb5REx2C226tcVNFvKcndTGh5F5AfvVctTA3z3VgL4pWsdu0fagfBvi1off340nUx0TEwNtTY/hIJowZiPHGPJTbOIio/j9I/+B1Q1TlXjO1Q1fO8O0FhGFv/Ok4iDaHOZ83jjRipDQljbrT+l+/eDopD71bfYQ0JQhw9nzBfva9a59+BRypNNkNidrmdN8MkLO/Vk1I0biaTKb+I/6olWbOxevaF+b1TQYBTLO4lK7kV0pQW3cTFQWw5GIQQi0EP+eHDLds6ShfVUY+PRqCYmSZL0wJvAUMAK3CTLco5X/gzgFsAOPC3L8uLGlK8lYT5yjIjO7T1OYX8bK37H/grD+5zGoCY/hvvYCIw+IEP37hwA+miU9b4/MzDi8vN9yvi/+dV0jypgwRkyK2j+2C02IlKcQQ7efjP/UUNNJtTGxC3XfiAsP59OiYk1FW+WNPYIYjIQLsvyOOBB4GV3hiRJ8cBdwMnAROA5SZKafvumZsqWee8A1e2nACtTz2ZbfF+8DTY1OelqGko3Jt73Ek/1tzi3nBWuv1qRKt71BPNQWJs8Er1RuOJaArtXbyBUVXyCLYL5zZtSOaTGT6azxUKnnj1rLN9caeyeMR74FUCW5TWSJHm/uo0GVsqybAWskiTlAEOA9YEqs1qtZGZmYrFYyMzMbEi5g6Ix5Yi85Ax4tvrDfdXvvxOXmIilrAz96NFAzcrBP7+pYnr8O7o7TWtkFMGJPwg8ymbu082i7fjTnNp2c5FB6RRFaWgkMTYzNpwmiBiars0Gwt227rzmWd6/73xyd9ZvWGtj/h6NrSDaAce8jh2SJBllWbZr5JUCsTVV1py2ZYTG35LQXFBM2lMv48jfQ0hyX8Y8OZOO+fmkpKSw7bP/q1N0BkHknyg1jV5qe+DXtz0Y4OQ7PmHlBRPRBYieSk9Pr6er1p3m1LabkwyWAwfY289EYtEhaoo9a4qXHe+2tWvdOuanpjbIdeppy9GgyjW2gijBqfTd6F3KQSsvBihuLMFaIpGdYjnl1Sc186QpE3GAx5FXhnPOgj9N6YeoqRM3RAd3d+BZQO83vmPl7RcHVA6C5odiq8TRLYFEq7lW31Jjj4rdbesN4LK8PPr06tUIV214GtsHsRKYBCBJ0lhgq1feOmCCJEnhkiTFAinAtkaWr8Vjt9nI+XMlWUPHgd7AzoQ+7E/fimqocudo2WtrC+urCzXVE6wDuT47tfe9JY2bzdRVq7jrjslCObQQFFsl5Ru2kvn2R0T6KQettlbbyLM+fW7evpDv77mff6sqXVuJcoDGH0F8B5wtSdIqnL/ddEmSZgI5siz/IEnSa8BynIrrEVmWLY0sX4vmSO5ukoeNIBznIno6oNehfDYtW003h7XGN/j6iPyoLTqqsZ3h3te7OW48WdffR/ZzkwgJCQl4jqB5YTdbsHaMY6S1atXUmtqtfxn/dlyfIwp3+5pIGIuKDjDZe1Z3K6FRFYQsywpwq1/yDq/8BcCCxpSptVCUt4+4ZOebi7tj2HV6zCFhDJ1183HVWZdw2GDK+ockNhTe15hFJ76duYCMp84hMjIy4DmC5odid7Cz/1D6eY0aanrga+VppZ8o3u1r5euv8/udd9Zj7c0LMVGuFaDYHVgGDQF8O8KaWx4gwmYhBG0Tktaooa6mpqq1nWpXAIGikoK9Tl3KbcVI0oyPuSNnDXkvTxbKoSVx5AjKddeT99vfJO3J8glrrWk+TaC2VV/+CO++kQsc2LSJ8a1YOYBQEK2Covx9dC4vqjZZbMzLj5KdPAS73kBGjxSfPK2PAmwM7+Spt7a5BFqT8rw7s3e6f57WzG6tugM9GPzL+vgZLnudXb/+xq63r6Zv3741nClobtj37kft3Bndxx+RdMGZWDDUGP0W6MXB28+m1eZqOlcLd9nPgWP79pGsqnQbOrQONbRMhIJoBcQlJZKVPBQFZ0PO65SIo7yckMhIEpf9zqqpN1F+2+2U7jvIslsfoMAYAa6yWZ16YNE7LY16oLuhyn/hj79CceNdthLQ2kaoNod4TZ0/GCX1Yp/TyIzpysk3fUTGBzdg6pkgnNAtDLvFRmUP52xj98M9EgcK2nN5tF4canr419RuA+F9zqtTH2aaqtK+W7cgzmwdiCmkrQCdXk9KZjqFO3eza1ceqWefik6v50jubuKSezEB4EtQH/43p7nPcX1MR3xXe40rPwr4vvG7OQr8ddENDLjuX/S6bDKRDlu1Mn/P/A9nzn2qxslu/iOTQB21tg7sGTHc+hnzbhiLaXhPVolZ0S2W3as3kITq06Y2durDiCM7q70s1OZfCCaSKZhgCjPOWPvKVauY2L59bbfQ6hAjiFaC3migsymJmJ7x6PR6LMVlPk5r74dyoLcx9/HhjCz2h0RhB4qANaYRVOLcZOW8nz4l5dJJHuXg5o0B55K3No0z5z5VTTbvumubvKdSXYHUVF/SHQvJfm0qU1L7YBTKoUXjvY2sBR2p17+F7vlnNc1FOREdqThSyI6ozrWaivzNlRB45OFmL3DvU59BaSndVZXe48Yd9321ZISCaKVsf/9LIPAyG/7Kwf3dgIpt5266lhWz4+/1xAJjszYQ4sqPdNg067kz41e6jUmtlmcGcuOTqqVrdVp/efxxl60EHu83ke9+XsbO16aKsNWWiHvjJq+NnfRGI5HHjrLxp6U8/Nk6/pl/DcOumYJDp/f89mbginvfp0/JISLjOmIq2scmY8eAwRX+baumeRHedYRlZ/PKo1cRFa01vbTtIBREK8BusZH27hesHXwyu39ZhmK3M/T2aZr+Au+0jN4DcbjSvTtO5pffs3L6bQw4NRUD2iYhrVFBqCvFO+/IyvX0Orzbc44C7A2N9qnLWykEGvq788dc+yaDZv/ITX++wyXnnYpeL5pwS0OxVWIbexJqYiLqaac5lYULY3goIyadxryrRhEZGYkhNBSlpIw/XniD7F+WUVZUxsKXr0dnd1D08+/o9TqGVhzm3Tv+7XONQOZMf/xHrBs69ORITg5dkpPr7X5bMmJM3sKxW2yYI8IY6U6YtYrS++4m0lyG7Wgp6+e+Cbl5jP7f81RW2Nn4n+cJTerByHtnkGIMQU7PoP/YIT5K5OxPqjZFrCkstaYZre6HfmThYbL7DsaUu5WcHhK73/kQ6aopUFgWlB/Cu+6kaxeQ/e51YsTQQrGbLeS/9wkd7ruLDjYLOsC+YiX6Q4fRJ8QHPC80OoIeF56BKSUFLBbM731I2Iwb6AAoOj1bv/yGm+Y722ywEz69FYO7fZUaQhl+OFes7uuF+E+0cHL/WokJ304Qo1aya/UG+pw+lglPzvakh7WDCe/MwW6xkfvXSiwFhcSnDuOax77hoycv9RlOaimGusacK8CRYWPpn5FGUf4+pKREOu3aS/vCfTU6prU69UsvfEbOzMuEn6ElYrdT9PMfxF48yWePDxXYEp9Mj6j2dA6mHosFNSKCCK869KpC4pVXeI79qS0Awj3a3frxIgZfOVkoBz/Ef6OFYy086vnueRPShQTcBMdSXIa9UxwmR9U2m2+j9+y45l+XFloPcK35D6ff/A7LusVWJe7fT9T0azEEqMP/Gm7WxvVh5szLMIjO2/Kw21FiY+lgdi6V4f+bv3v1A8yPCW7bF8v/LSZMo47vzrqcm377GKg5NDpQm7YWFjK0Y8egZGhrBGXAlSTpAr/jyxpGHEFdGXj5BR4nngL88vxrRNvMmm9CdosNe+dORLkikKpizZVqb1QVQHZkHOA7FPd36K1t39sTeVKqD+WtIROZP2wSpUNH8Pf7t6M7/XR2pIykfXIv6NGD8OV/e+rUwl3/jNNvYVVXE6nXzGdA7gahHFoo9m0Z6MxmTV+TGXj96atrnK+i2B0UyLso3X2Apf1SNefiuJUDaPscvIMb9hLKZmMH/jz7UlbcNIvK0lIihHIISI29zqUYTgaulCTpJFeyAbgI+KqBZRMEgSE0FFtJGTu+/42UqZNIyg1sQ929egNJdmtA85D7ePuiH+l85umsXryKftee41Pm1xvvw/TITLYsT0fp2IkLzhrG6GeX0HVfHru79OSW0/pxpmEvkRPPRac4cKxchUlxYCS42HMr0P/2L3h1xlg+3XIDf186kKioqLr+WwTNhJzOPelmCCPGYfWkuX9rc14eUQH8SYrdwf4tmdgnnk+vI7uJAzL7DkGa+RVnZa/jqhWLGH80v8b25J23/IX/MX7mTfQwGulRD/fVVqjttWwzEIfzhVJ2pSnAFw0plKBuhEZHMGja5FrLOePMI4ixOTftzO3Si+g/f6Pokssx5Wxmf4d4Ou/MZJBrQtCVl53C+vuTST2Ugwps7zWAc995AZ1eT98k16BSUUi/ZQQ5uvGY4mPQ6/WsSreR3j2FEfsy2dC9P8P3bK81Oglg4hUvwsjh5NxzKkajkYuH9T6Rf4ugGWBKaMeI+z7jvPW/0v1oAQNeup/9j7/GeQvn0rl7d81zdq9JJ3bcaLp7rUqsA6Sd2/h7xkg2l43mpNSXaryuf6j1+Fkz0BsMNZwh0KJGBSHL8h7gI0mSPnEl6YFxQEZDCyaof5xx5sXkrlhPbI8E+vbrjU6vp0tmOkX5++iRlIjOK2w0NCyMAdlb+ez3jZxu6kh8uJ+tWFFQTz8dw6pVSCedhG7pUgDaRxh58qG32LUjj2ExKm89c01A5aAAb/Y6mZ5vP88nY0fSpV24WCKjtaAolK/bSPoL/0IPOHR6lNFfYPnyDbY88BSp858mrF07n1N2r0mnxzin/8zfJ7aze1/6mXpT8eIr1Xxm/ni/eBSuW0dPoRyOi2ANu88BO4FewAjgEHBdQwklqD/sFhu7V2+g94RR6I1GjOGhJJ91sk8ZvdFAp77am6pHx4Rz1cWj2TFgFKbcrWT2HUz/jDT0RgPKocMoK1ZiVBw+4Yo6nY7PbzmJwvJRRD6nveOd2y7c995vWX7vKHr0EAP/VoXFgn3QYKJzcwDXJExVYcVTL3HqnMc4GeDT1ynZd5h23ZwxTIrdQcfxJ3vKu3GPAvpmb6Fy126MTz7lk6flczgClIdEY1zzNz1HjGiIO2wTBDvLaLwsy28D42RZPhdIbECZBCeIpaiEzf+dy8GtmVgjwkg6YxxlUe2wW2zaJ2jMavWmKH8fptytGBUHptytFOXvA6Awqj3p3VOo1BtI755CYVTVWjV6vY7OUSEc/fYnoLpzWzbGkJOzm10vTxbKobVhNqNGRGDIzanmND51zmPglbblAefDPueXPykMjybKtbGVd3sB2P7771QMGERIchI9zcWePP/gChXYEN6ZDlYrvW2lJArlcEIEqyAMkiSNBvIkSQqF4MKWBQ1Pcf5+Vow9l4z4PqwYcRr7t2whLC6WIU/MouuQAUTi7EQxtgryVqyvdr7dbKF8yHCUbt1QTjnFZ1arm7ikRLL6DsauN5DVdzBxSc73g04xYbx6/xtce8XTzJs9n07mYh8loxw6THz2Nk+HdwBHjBGkXvc/uh7Kw9S3hzAntTbsdpS4qug3N/4zm93HqfOfJueXP+k76Ww6Oao2kPQus2rOfEaecw5Ru3J9ou/8gy3e730SK35ewbDS/RhDQ+v7ztokwZqYPgZeB24A5gCvNphEgqA5vCOXzinJeAxGh3bBFVVhpP7D79geCT7nl+wvwNA93hPmqq5YgTlrF5H9ffdQ0On1nsluKW4/haJQuWU7H997JgaHA740wD2gnnwyvPkm4Bxh7HQ5q9O79efY+++z3hLBirP7ER4e3hD/EkETY9+WgcFiqfbw3tKxJ8YO7RmQuwVw+p6sBQVEtmtHzOSLgeoPfKveSKVOx0mz7/DJ98Yzy/7GD8l680pChWKoV4JSELIsvwm8CSBJ0ssu57WgCbGVVdApxblejP+bmveqre4OtK1HCoP69faUMx85Rkz3LtXOT//uDyY8VH2THR8/hc2GfchQQuQdVddxODxLJ+gKCwHnCOPfLmd175QkFp6VykQxYmjV5HTuSTdjODH2qtHAmJveZfnr06g8ZmbFY88TkzqCwddeQqTrYb68xxAuzV0DVLXX7MhOJJmPaE6Mc6MChSHhrPx7MzliNd8GIaj/qCRJd+EMdW0PTJck6VdZlmfW5UKSJMUCnwLtgFBgpizLq/3KvIZz3kWpK+liWZaP1eU6bYWs739jIIGVg7c34dqZ7zH/rolsmXINXf/7EF2HDGDLvHcYQ/XJS2PvqSX2wG5Had8eQ0VFtWspQHq3/jhcvghvZ3Wn6FBhTmoDmBLaMfyhr+meL9OzpIBbX/w3y9qFs23mI4z431zGA7wDFZOLiezkVBBnLf8WunXzabP9zEdqXdkXIPLAXi52mbQE9U+wKvdK4FTgV2AgsOQ4rjUTWCLL8iuSJEk451L4e5BGABNlWT5yHPW3KVKmTsJxtXPWIlS3yR7Rh5M5/mySPnmHVyrMxPROZAjAD5+zo+cAhm9dA8/O9nS0ZU/MY/zsWwiJiCAQit1B0cdfEKehHHQ4wxjfvfM57o2oalZ6vY7OQS6lIGj56PV60h87j5yCCZi6RmMrKiGscwdPR/eswTTvHcY+cz8A0R3iyI2IpW/FsaBeeMAZ1WTfv59YoRwalGCd1CqQABySZVkFjmdu+jzgbdd3I2DxzpQkSQ/0A96RJGmlJEk3HMc12gyG0FCUUjPfPPYK4BvF4QA6lB/llL9/oGvHWNq7fApuJWLanUFpwVEqCopZ+/Ac7OXlnP74PbUqh6yuScTdeK3nem42dOzljGRKHMAzt54lRgqtFLvFxs6la1Ds9hrLGY16+ie0Q6/Xk/WC013p75we+YjTr6DYHeT0NFVTDlr8Pm4SOSvX8e17PxFSUUFsQkItZwhOFJ0aILTRG0mSngGuxjmSmAoUybJcfeuwqvI3Avf6JU+XZXm9JEnxwC/APbIs/+11TgxwNzAX54vxUuAGWZa3BLrOpk2b1LCwMCwWS7NwejaFHA6Hg0s+kjlv81LOkVdhPWkMyffdREh4OJbcXApfXsCZ/yz26ZybE1MI/XWRz6S42ihauY6TZlxf7Q0vM64n9876H0cPFNGldzxzzu2G1Wpt8t+jPn4Ls9mcPnLkSO1VDxuY5tS2iw4d4uhnP3D6B/OJdNgoDYskf+2qoCKFcvYVc8HZJ3mOLXoj8qoVnglyJfn7GX3eWZpzGbxZ2zGJmL9/aNLZ0M3ht6gvOYJu26qq1uljMplC6nqO17mDTSbTdpPJdJ5GnsFkMsV4Hc8xmUzX1FRfRkaGqqqq6v7b1DSVHAU5+aoCns8xQ6hqTR3tk6aAagd11PS31MqSEtX6yefq1g+/Vu1Wa82VOxyqIy9ftbvqUL0+CqiP3j5Xtdsd6uESi6ooiqqqzeP3qA8Z0tLS0tTjbOsn+mkubbv0YKFPG3L/7rl/ra79ZKtVLft0oaftVIJaWVrqU+Rw1i6fduV9Dfm3perKh55TM9ZuUbdt3dpAdxg8Tf1buGnMtl3bYn2r0fATSZKELMsnaZxSU10DgEXA5bIsb9YoYgIWSpI0AqfpazzwUV2u0RaxW2xEJycBVfbaGIcN1q/zeStT0PHXH+tZM64/+uhoDDidSY7rddhKywmN1jAv2e2oEyagW7++2oQn9zIZ5/37SgwGvfAztELsFhtHUwYRhd/aRsawgMvJe7DZUCMiiFSq1lMyADtffoOkR+5DsSvsXLoK84FDhBtCifZaYVgFdke0p99Zp2A65zQAMjMzG+AOBbVRm5P6ipoyJUkaI8vy2iCv9RwQDrzq9FFzTJbliyVJmgnkyLL8gyRJnwFrcK7M+7Esy9uDrLvNsnv1BpL8lus2oyfCPbfBlZ7ZdzCnjO5P9j0PIVHV4Q2o7Pj+t+qL/SkKyvgJ6NauqRZN4v7e984vye0n5ky2RhS7g5yUEUhHD/i0LQcQeqyo1o11yn/9k0hFqWY66vPEQ5Q+8Qgheh39FIcn3b9cyMoVdTKBChqG2hbry6/l/OeAM4K5kCzLFwdIn+v1fQ7OiXgCnHv3Hly2kp1/rCA6xcTgqydj8LP7+q/Q6gB0RwtJHzqBYbu3sS+mM6H//EWfnj0xxMYgeT3q3aOAlKmTql/70GGUdes8y3T78+6lt5L7yqUYxCJorZKi/H30zc/02QckNy6RpL3ZhARh/16aXcj5ru/+m/jEoIBSPULJM2enWz8GDU058ZsQnDAnOrNEhKs0EJYjxTg6dyQBFXeshuNGsJWafcxB7hVaM39fQU72PnqdPpCh7dszPHsju3bk0XdgEqiwq6eJPqjVOqWlsJAoDWdjYVR7crv1Z8w+58K9/uedP/8ZoRxaMXFJiWQmD6Ffzhaye5jo/PtPJLtW/w2Gc27/F477rvUJw/Z+0fA+9n6ILL3nv5z+8qNi9NBMONFfofYQKEGdsRSVENq5A5GuB7r7Y8A5Qc58+Chr75jN2hfmY7dYMIaHMuCiM7ho1jWEukJVjaFG+g1JRm8wUJS/j54H8/zeBntgOXaMqAC7aXWKCePTq+8DfHeTM6PjwmcXk9ClQwP/FwRNiXt5lWPZu+Dnr+ks9anbQ7tSYdsLb/LNnU8ANSsH7xDtCc/eL5RDM0L8Es0Mxe7gYMrgahPf3Pbfnv2TiOhH1FjYAAAgAElEQVTakdFvvsjoB+9EHxGBpbisxjrjkhLJTh6CXadnZ0IfDm/LolfONgre+oCc3/7RjGvX6XS88uTVlIY5txMtM4ZzxW1vUniggB8fnCTmOrQB3Mur1OWBbbfY2PL5D4TERjP0gdu59I0n2Na1D7+fdI6PKcmiM+Cgqm3fM+Y67OXlNc7FETQ+wsTUgCh2B0X5+4jz24inJory95F4eI+PclCALRdPJfmn74kZNQzw/cdvf/9LRs68MWCd3ovt9ekRj2N9GoZBJs/Wi6WhEUQeK8YY7mtqMoaGEFl8lNxV6cQMH8LC9pFCMQgCYrfYqIjtwGCbGahqo/0P7WTQoZ2eNBWQ3/2cAZddSMb3v9Ll3LN4pWO0aFvNkBqfWpIkfeD6e0uAIp/Xu0StALvFRub3f5Cf0IfY5N5k9h+JYnfUfiLOt3257+Aqk44+BD0w7PtFRNttmiOLobdPq7VevdFApx7x0LEjBq9NWdxLge9evUHzPGN4KMlnjKNrhyjRgQU1snv1BqJt5mpt1Ow68m49UYkJhEZHMGTaFOLjYkTbaqbUNoIYJknSi8BUSZJ6eWfIsvywLMsLGk60lknJ/gJCEhPorzoVgg7ol72Zovx9AXdt80an15OyYwOH5F0cyNuLqfgAXH2V5tr6ZkMIxqOFhAUzq1JRcIwZi768vFpdZSHhtce1CwS14IyoiyTGNYIA2NixFz1X/Qb9+3vabak+hL5njGsaIQV1oja7xxRgM86VXGW/j8CP4oMHienehXDV4TPpZ2+XHp5NdoJBbzQQPzCZ4eefxj8p432ceCqwIbobG374i/CKcsJiYoKqUzl0GHXz5uqmq449CQsirl0gqA1nRN1Rdnz3BytnP8OONVsYdjiXTpLErtVp7OjSmy0fLSLaahbtrYVQ26/0lizL50qSdIYsy2JWcw3YyioYfcaZQHUTUOKuzOOOzDh7SAKmWV9z4eZf+dfmv5l7+Sy+mncjhjp2MPfmPal7t6MDNiSYePaGp1j42GRCxCYrgnrCGB5KyuSzYPJZPul9xo50bmglaFHU9pRpJ0nSImCCJEk+/gZZlq9qOLFaHs79GVSfEL6dsfH03J9LSGTkcddrNBrZ9vTF/JU1jj49Yvn6OB3FnWLC+PfDb7MrI5ekztHMu/McvhZ+BYFAUAO1KYjzgMFAMlVLdQs0SJk6Ccc1egyqgoKOg2mb6TN8YL3EdIeHG5k0pNsJ1SE27xEcN4qCRc5l+09LGXr71RibwYqmgsahNgURI8vyCkmSrgWsjSFQS8UQGsrmtPUYMneTMnUS3Zuh2UZs3iOoK4qtki4XXERY3k5GAOqsGViOlhLePrqpRRM0ArW93rq3FX0VeMvvI/AjNCKCQdOqr5ckELREFLuD7L6D6Ji302dG//b3v2xiyQSNRW0jiCmSJE2myu9aCYTgtxucQCBofRTl7yN5b1a1dbiCmXcjaB3UNoLoDwwA/sK5j4MJuARY3tCCCQSCpiUuKZEdSYM80XgA5QcPCh9EG6JGBSHLslWWZQvQV5blda60jTgVh0AgaMXo9HpSsjbx1+eLWHrlbViLi4nu2rWpxRI0IsEG0xdLkvQUsA4YB+Q1mEQCgaDZoDca6DZsIClX/qupRRE0AcHGYE4DDuIMez0ITG8wiQQCgUDQLAhqBCHLcjkwv4FlEQgEAkEzQuwHIRAIBAJNhIIQCAQCgSaNtqSiJEk6YC+Q7UpaLcvyQ35lHgfOB+zAPe7IKYFAIBA0Po255m5fYIMsyxdqZUqSNAI4FRgD9AC+AVIbTzyBQCAQeNOYJqaRQHdJkpZKkvSzJEmSX/544HdZllVZlncDRkmSOjeifAKBQCDwokFGEJIk3Qjc65d8B/CcLMuLJEkaD3yK7wihHVDodVwKxAIFga5jtVrJzMzEYrGQmZlZP8KfAM1BDiFD85HhRGhObVvI0HxkaGw5GkRByLL8HvCed5okSZE4fQu4VojtLkmSTpZl9yz+EsB7e7QYoLim64SFhZGSkkJmZiYpKSn1dwPHSXOQQ8hQfzKkp6fXkzR1pzm1bSFD85GhvuQItm03ponpceAeAEmShgK7vZQDwEpgoiRJekmSegJ6WZaPNKJ8AoFAIPCiMZ3UzwOfSpLkjlK6HkCSpDnA17Isr5MkaTmwGqfiuqMRZRMIBAKBH42mIGRZPoozhNU/fbbX9yeAJxpLJoFAIBAERkyUEwgEAoEmQkEIBAKBQBOhIAQCgUCgiVAQAoFAINBEKAiBQCAQaCIUhEAgEAg0EQpCIBAIBJoIBSEQCAQCTYSCEAgEAoEmQkEIBAKBQBOhIAQCgUCgiVAQAoFAINBEKAiBQCAQaCIUhEAgEAg0EQpCIBAIBJoIBSEQCAQCTYSCEAgEAoEmQkEIBAKBQBOhIAQCgUCgiU5V1aaW4bhJT08vAPKbWg5Bq6XXyJEjOzfFhUXbFjQwQbXtFq0gBAKBQNBwCBOTQCAQCDQRCkIgEAgEmggFIRAIBAJNhIIQCAQCgSZCQQgEAoFAE6EgBAKBQKCJUBACgUAg0MTY1AIIGgdJkk4DXgXKgWggVZZla5MKJRCcAJIk3QHM8EoaALwgy/J/mkikVoeYKNdGcCmIJUAfWZbFDF1Bq0KSpNuAG4FTZVkub2p5WgtiBNG22COUg6C1IUnSFOA+4GShHOoXoSDaFmVNLYBAUJ9IknQyMB84S5blg00tT2tDOKkFAkGLRJKkFGARcJUsyxlNLU9rRIwgBAJBS+UVIBR4SZIk97MsTZblm5pQplaFcFILBAKBQBNhYhIIBAKBJkJBCAQCgUAToSAEAoFAoIlQEAKBQCDQpEVHMW3atEkNCwvDarUSFhbW1OI0CzmEDPUng9lsPtJUe1I3p7YtZGg+MtSXHMG27RatIMLCwkhJSSEzM5OUlJSmFqdZyCFkqD8Z0tPTm2zWeXNq20KG5iNDfckRbNsWJiaBQCAQaCIUhEAgEAg0adEmJoFAIKhvikstzP87l1ln9aXE4uChbzehR8ctQ0KbWrRGRygIgUAgcLFl7yEueiMNgAUr8nzy/siCdcnldGkf1QSSNQ3CxCQQCFoVxaUWnlm8HYvFUi1PUVQKSq2szdvLgEd/4rJXfqG83LlC+O4jxzzKIRBPLM5sEJmbK2IE0dYpKYF582D6dOjRA3S6JhHDUlzGpvnvoZSUkfr4vYRERjaJHIKWTV5BMae9vBJwjgA2PXIm7WPCsdsVMg+W8OA3m9l+oGrV+3UHFQY+tSzo+l+c0vRRTI2JGEG0ZUpKUGNjUZ94ArVXL5TRo+HAAQi0gKOiwJ498PPP8OabmB98lPKQMPImXkR5aBj7r7kJpbwc3cGDFK5JQ921CyU1lZK7Z6JaLOSt38T2Hv0pzsvzqdZSXEZohxjGPHoP4+Y8ijEqCvORYw1//4IWhcVi5/tNe1m8aQ92u71afnGphbOf/4Pp/3zCb69fD0eOcN3bS8g5eITkR3/hwjdWkrG/hE7lR0FV0amK53ttjO4extf/6kpUVNsxL4EYQbRdFIXiWfcRC3jGDGlpKN26AaDPz2fh/Le4cM5z/DFjBuPf/4hwVSVKqfRUEeH62+v3HwGI/PQ9+PQ9JK/L6ICYtDR4bR693IlJSeSt20jv1GEAbHrzA8Z4yaECW+a9w9hn7q/nmxa0VCx2O/2f+M0rZQs7nphIeHjVI2z+L1vJmnuZ53jXe9cz+5TrOOugnRDVxuRty7g67WcGFOWTHhnPwIojRKiVpCUO5MqrnkPVab8vGxw2vhwcxY7w8Ia6vWaLUBBtgf37UaZPp+iFucSGGig/5xzYt4dojaJ6nA9otVcvLnelXbxggSe/NgOUiu+D3vu7+3wVKL/kCtizw3nNkBCf8wFGPnIHRXn72DztBiKLihm49Hui4+ODuVtBK+SHbQeqpT343WamDEtgvKkrBoOBWQXrAN82N+efj3jhn4982q0OGGM+UPV973Y6lx7hcLsuRJUf5unF81HQM/u8uzh1TzrvLn4F3UvQH8BqhdC2E80kFERrRlFg0ybUkSPRAXG/DwIg1pXtfli7H+Q6v3Tv78F6Jvw7oj/ueosGDabiww850D2ZYbddjTr7Nk8Zy5EjVBw+SoekRE5zJyYkcDAjm/iUZN8KjxyB+++H+fNB+C1aLYvlSp9jvWJH/m05878v5/oeg/nyjpG8sbyQj1z5Wm1W62XFzYL/3YklJJzRlUWetEv+d011QX79FS666PhvpIUhFEQrRbFYqRw+nNAdzqgLrbf6mvAvc7wKw7s+b0vv+F+/hl+/Jgkw60PQAwZXXviUKaRHdOVkP7nzr5pO/MblzoSiIsqnXkHkX384jz/8EF1BAXTqVEfJBC2BuWe245qfStCpCt2O7uHP9+8m3OH0QzjQIfENdBsA+LYzf2WgAgrVna9DMEOluVrZau38n3+EghC0bEry9mFMSiSCwA9y/8av9VaFV5p3+eONcwp0XqTLr+HOdyxfzuDPvoDfv/bp7P03reCoTkcoEOn6eMtXcuddtFv4+XFKJ2jOGKOjSTiwjD8/fszndwcwoDJy3w5CbeU+6VqKAqorh0AjXc32+vTTdZK7pSOimFoZZYcPE6OhHNympHICD79revDXV/Cr6vXRqlvF2SjbTbuymlztgPb4Kgbvc9+1itFDa6SkzMp9c/+PVS7l4P27u0cEBx2lfPTNMz7nBXqhqamte5tcfbj6anZs2ABtzFEtFEQrw/LSO4CvScf9N2Phd2gF6dXkK6hvdFTv4P7KQqtj6zQ+/lz5z1f1Kqug6Tl4tIxRTyxm+adPAtqj3t9ie/P3V8+h18j3bkvBmlY1y338cZtTDtDIJiZJkvTAm8BQwArcJMtyjlf+DOAWwA48Lcvy4saUrzVQdMetsPhzT+exAznhcXSSN9G/Vw+gdn/E8fgY6oq/Y9yfulzffa/33fsWn52QVILmgqKorJHzeevht1m++JVqLxXenHcsDwisHLw5rnbduXOTTSBtahp7BDEZCJdleRzwIPCyO0OSpHjgLuBkYCLwnCRJTb87RwsjuUdHRt3/BW+kTiHlnq+Ys+BPehXkUfDyG9XesPwdx97pDU19mqwAkqa+xlt3n1lPtQqaEpvNwaTHFjFmQBIf/d/zdLFbApqJAo0o67UNb99en7W1KBrbST0e+BVAluU1kiSN8sobDayUZdkKWCVJygGGAOsbWcYWjV6vZ+2zl5FTMIk7ukZTlLubiq7d6W8u0Sx/vFFJTY1bMRwFUm/9jK1zLiQmJqYpRRLUA3a7wuj/fMuGOZdXGzVomU0bvN0eOOAcQbRRGltBtAO811BwSJJklGXZrpFXSlXIviZWq5XMzEwsFguZmU2/iFZzkMNisZCdLQPw4/d/c9FtzvkFgUYL7jx/auqAtSkU/5DYQPUcD973kAqc/cCPLL6uL3v37vUp1xx+ixOhObXtxpQhr6CclS9PqzYq0Br5NoRycLev/AcfpGLaNDh61Pmh+bSpxpSjsRVECeD9mqd3KQetvBiguKbKmtO2jNA8tiR0y7B/w1Yf5aCFv8KoaSKR9zlacyS0FFB9d2J3/WZg4Gn/YfHrlzFo0CDNsvW0LeMJnX8iNKe23VgyHDxaxgefL+JcR2VAf8OJhlrXhg5AlultMlXLaw6/RX3JEWzbbmwFsRK4EPhKkqSxwFavvHXAM5IkhQNhQAqwrZHlaxWUHSoiYeQQILgIJa0yO+L7IB3cGXTkh0rDvdm55bUB0s2f8M3NKeSNHFnPVxE0JQePljH2hb/R29vzmM5IlFp9Mb5GMYMeOABiSRcPja0gvgPOliRpFc7fe7okSTOBHFmWf5Ak6TVgOU7n+SOyLFdf0F0QELvFxr5f/qbnLN+RQ6ARQaDZ0psN7RnsUg4OqmY4a5UNZhJdoDWZgsFdfmKvKXDZdHY+fz56vYjObm3c/81mjHYLm+Ze4VEOjeofO3y4TfsaAtGoCkKWZQW41S95h1f+AmABgjpjt9ioiGnHmXYrUN3kE+gB7T+a+OPxZzn7vw97yhj8yntTF/9EoIlxgfAubzrnKbZ+fx/hbTAOvS1QVm5j9Y6DbH75X56JcG4aXDmMHQtLl7bJOQ7BIJbaaCXsXr6OJLs14CghkLJQgY/jh3Fm8U523ng7Y72UQ01vbzU5vbVGDP7XDoR3+cfDB3Lxqi/JGjYAXRuNQ2/tKIrKhXOXsPnlqT7KoVFGDv/8AxMmNPRVWjRCQbQSEjtUrWQaaA0a0H47u+7gJmcd85/3pGs5nP3ztGZDB5qsFEyHd9exH1jxv/f47y3ThWJo5RSWWnj1tduIRG3ckUNREXTo0NBXafEIY24rYWdCH0pDI1EBq87A5Rc/SsnBw2z43yeseeBZiiO1dn+oItCEI//lwP1nQAfyKfgrqWCVw5gHFhF+5AiX33qDUA6tGEVROXTMQsb2bQw8kh9wxFnvfPihcwc5oRyCQowgWjJmMyxcCNdcgymhHSNmL8SUu41NvQax5aFTMa5cTcf20USdewqxLzwccNJRTZ2zpuUL/I/rMjPbv5wCSHd/TdazU4QTupWTdbCQ819ZQyWAqrKoSx9GHd7pyQ/UPk8Y4YiuM0JBtETsdsxLlhFx7tnO4xtvRJ+dzYbfn4G09Zi7diP0+b3ogN5+p7o7nQM4GhpJqM3smXyiZQ6qadXV2t76auvYZuBwWDseOGUG/37hZrKGJgvl0MrJOVTEOa+sweCwMWZvBrkdEjka5myB3u1JwTkJqgP14I/YvBmGDDmRGtosQkG0NOx2lA4diCgrA6retNR+/Tymn6gDez15blTAjo68Dt0497pXuWFgO2J7JjBt0ihwWKt1Tu/HdKCJcVrUNCnPX6EMnPg0oy6YwPvXDqNdu3Y13LSgtXDDB2l0LshgzfuzPW3M4fpb9fKi45xpT7Pks0eAE1AOkyfDF1+ICKUTQCiIFkb5qvVElpXVGLaq9UDPjU+i8udfUDt3IatbLHq9niO5u4lyKQetDVU8yofqpiQIPLrwx3/9HICkEXez5sMbiReTklo9FoudpVkFjOsdS9nOTNLfnw1UtRmD13c7cNKMd1ixYIZPmTpxzjnw7bcQpbW4vaAuCAXRwli5dAMuw1KNPgRvO64C9Nq5nZCICJ+64pISyeg3DCl7E2UYWdJjCJfs2VBjNIn3wz6YZQ/c5cfd8wXtD2znjcVvcf7lL7LmmXOFcmjFKIpKbmEJn6/M44M1Vetk/bDoRSCwOdMA/LVgBiEcp3KwWiE09HjFFvghFEQL47T7p6M+cScQnJ/AbAgl5FgRoX7KAUCn15OSkUbBzj0URbVHUhRKe3YhxjXo1/QxhIVRsSaNbf99gdT/+7RWB7QCnD57ESufvYQicyWdop9EFtFJrRpFUTn/5T/ILKz0SdcrdnqXFgA1j3zr/N4/ejQ8+ihMmgQGrXn/guNFeARbGKGRkaxb+JPn2LMUxZUvsb3PEI9JSAWyf/6LCEs5oa6htq2sgm2f/R8Om81zvt5ooKupNynd29MpwohB7+ymWqMIG4DVSsTwwaT+36c1yrm+fQ8uv+BhsrLz+Pv5SzEY9HSOCROhq22A/KIyj3LQK3YGHNjB2TtWIO3LIsZh8wmp9jZhBppPE5APP4TKSli7Fi68UCiHBkCMIFogqZecTWlIBDGVFZh1Rgbd/SXyCxdg1N/DysW/EZ67j2F3XIPJ5ZyzFJex8Y33SP3PPQwEHNfosZWUERrtO6ow5mQRqdgDmpVCqX2GtQpsCOnIis9/5ZMz+hIWJvZ8amv8sHEPYdYSbl7+Fbem/x+RGmXcbUqhyrQEdTArlZSA2P+jwQlqBCFJ0gV+x5c1jDiCYDCEhBBefJS0H5eyMXMPOS9dTGhoKHqjgTgpiVGzZmB0KYfDWbsI6xDD2P/c43EGGlSFrO9/q1ZvXOpwzOFRPqMQLfNVILOSCmwI64zjn5+Yed4AoRzaIIqi8tVfm9nxylXMdCkH/4mV7s8VFz7EpRfOrKYctMKmVXCOEKZNg/JyoRwaiRpHEC7FcDJwpSRJJ7mSDcBFgNghvgkJjQwj9YLTtDMVBQoKKLHr6Cz1AXw7nwNImTqp2mk6g4GIY0cpWptO+/Jj5FcobPz6N/r++CVDSw/WOmv6trvf4OWnr2P3nj0ndG+ClkthuY2Llju3kvd/4KtA/39/yrUbfyHUamP2svcZUXrIp6z/dzcWk4mIjAxhRmpkajMxbQbigApAdqUpwBcNKZTgOLHbMWzdSsXlVxGeuQ2lQxfAt6Na9UYMpccwBoj00IeGEDdhLNjtJK5Lp6BzJENLD9bojFaAAlnmLY1NVgRti07RoaRfMh1Wfu4zErh14l0s6X8SW+dfS7i9ygcWlEkpO5s8m40UoRwanRoVhCzLe4CPJEn6xJWkB8YBGQ0tmKCO2O0o0TH0szq30NAB0QUHAL95CPv2UvLWe+yPi2fAlRdj0FIUNhtK+/aEVFQwOsDl3HVefvZ9fPjNY8SLIb8A0KkqXy58uOoYZ1vJT+jDpteuJly112hKqsaSJZCcDM1gq8+2SLBO6ueAnUAvYARwCLiuoYQS1A27xUbu/AWYrBafzpcR35foZX+x579PEzZoIIOmXUJUQjxhQEfAcb0OW2m5x1mds+Qfwi68mHbhBtpVVAScX+Gmzx0LyZp3CSEhIQ19i4IWgKKoZG/LInn9eo9fQQUsgOpQiVLttZqSfNDr4dRTG0JUQZAEG+Y6Xpblt4FxsiyfCyQ2oEyCOmC32DC374Dpvjs9aSpgxsDTj75PkimRUz9/mxF3zSCidw+gymloQPU4q3OW/EPfs04lsaKYdkcLNe3Hznrh5gtmszlrH7mvTRXKoQ2gqCoFpVZUNfA7f1FJBSOe+JmJn+WwvlOyT4BDOPD15w/W/cIVFcLn0MQEO4IwSJI0GsiTJCkUEEsi1oTNBqtWYRs6gtwvvsV045UY6juix2bDtuQvNvyziTFWc7U48sL8PXzZIx6dqsLhw2T9sYqBVF8Qze2sDpt8CVDdsWg2hDD6rk8YvD8bfVwcs2ddwtu9O4tF9doIZnMlN329iwPmXYzo0Y6vbxuPXl/17m+zOfhx+25mfemyOut0XHnDXEZt/4uvfnrFmQREO6woQ4Zg2LKl9ot+/z2cf75QDs2AYBXEx8DrwA3AHODVBpOoJWOzofz8M7rzzwcgBOgPqHdMx3asjNB29bM2jL24BH3HDoSoCmNcaSpQGhrB71Ou44J3XiDRVknJ9Tdh+HsZUbt3MWDcuGr7Sx/auJFuoaGgKKQ9+yqJd13tE6lkB3av3cSWYf1ds6BDxUS3NoTFYmfAk797jjfsKWHnoRKSE2I9+f2fqB4urer0zF3yPuBrRjLcfLNzAb0ePZx7Mrg5eBAiIuCVV+CBB0CERzcbgnoNlGX5TVmWx8iyvB14WZbl9xpYrhaHYq5AGj7coxy85wzogJxPTywq2FZ0jK2338uKcedCh1h0quITX7724TlEl5cw8PG7CLVVouvciXYfv09U/k50qoqyahX2rJ38POtxstoncHj7droNGwaKgnLqaUy+62ofed2Oxf7D+4tZ0G2UpVkF6FSF7oV5zFryHlgsLExzrqukKCoXv7k84LkTLp8L+M2MvvNO+Ne/nPuYvP02dOniVA5du0K7dvDYY0I5NDOCGkFIknQXzlDX9sB0SZJ+lWV5Zl0uJElSLPAp0A7npNyZsiyv9ivzGs55F6WupItlWT5Wl+s0BYrdQX7yIHqrqo+Jxtvs0/GyS4+7/rL8/UT17s4grzR/U1Hq43ejNzp/zqJ77yfOr4zcqTdduiZw/ktPwEtPeOqx5+/BsGJ5NdOSAx29dm5HJ0xJbZaz+8fx5cezSD2YDcCdad9hfdAZGVdYbkM+bK52TkRFEbOXfsbefgMADUf0mjVw7BjcfLPzI2jWBNv7rwQ+As4DBgLDjuNaM4ElsiyfClwPzNcoMwKYKMvyaa5Ps1cOAEX5+0g8sEvTsQsw+ekf6BpX9zBQxWJl/yeLiOrdHQg8I9VaVITBaMSSIXPo6bl0+PQDnzIO4PH/fECnGL+3M7sd+kvVrnvrpJk4yss8s7EFbRPjoUOkHsz2aXeO515g7m872F9YwsgeVXt46FSFngXZZLx2Lddv/Y1Hv52nHaXUr59z5CBoEQTrg1CBBOCQLMuqJEkdj+Na8wCr13Ut3pmSJOmBfsA7kiR1Bd6TZfn947hOoxOXlEhG38Gk5G5BB2zp2pfZNz7DI/kr6PPCY3zfrVOdzTOWg0cISehMguvYf2RiB56ZfDcPfPI0keHh2EeOImzTRk71Kq8Cj590Ddd/MYf5ZWXk/LmS9j270alfb3R6PeZ16UTYrNU68s23XkBYpNYKOoK2ROGmrfh39FPNgziyNJfXluYSGaJj+QOncbSwhD5nnkzU7l3aSsGb7dtBmCpbDMEqiKXAPziX3JgHfFNTYUmSbgTu9UueLsvyekmS4nGamu7xy4/C6Qifi9OXulSSpDRZlgOGPVitVjIzM7FYLGQ29USa//uMpdl57CvX07tfAvNiwtDphlBacoQdJUfqVJWluJhhJzlXNvHvSh7/QPtunPXQ9exasYLuN91EzL591Zb6VtAx5uJU9IOG0KW0APd727akQeh/+IKMJev5l1d5gNKQCMJ6xZ3Q/7M5/B7NQYYToanbtsVu57J/7GTq9BhUBYC1nZM50rGrp0yFzUHWT79yxjXOVqS1BMu+adPo8PXXFN5zD+XTpkFOzvHJ0wx+z+YgQ6PLoapqnT4mkymkrud4nTvYZDJtN5lM52nkGUwmU4zX8RyTyXRNTfVlZGSoqqqq7r9NTX3I4ai0q3kdu6mKM85DVUFVXB/v46se+ES1Dx6imaeA6gC1RONcFdRKdGpBTr5qq6hQK3V6VQHVjk694c431UqrtVn8H5qDDGlpaWnqceR9o6kAAA+pSURBVLb1E/00Vdu2Wu3qsh0H1V4PLFZ7PbBY7XPft+plVzytjrztQ7XX7B896b1n/6CuTZCqtS2fz5Ahqqoo9SJXa2lT9UFjtu3aFutbjcaMeEmSkGX5JI1TaqprALAIuFyW5c0aRUzAQkmSRuD0jYzH6fdoUxTl76N70X7NLUXd39d16ccnr9+CzmyuVi4/uhPPTriO+7b8QL992ZpvdTm9UpCSEp1mpqISVn3+IwOmnse7ndqJSKU2it2ukL6nkCveXufT4UMqzUzITSd3qO/uf6lhVlIPyIFNSkYjbNggzEktnNpMTFfUlClJ0hhZltcGea3ncE6qfFWSJIBjsixfLEnSTCBHluUfJEn6DFgDVAIfu8Jq2xROf8YQBuQ6LWtmIK/PYKRd2ymMiOGG6+fwzfUj0Y0eVc1hPfXBL5h+XjJPd2pPx4Fzq5mcADK69yMlZ7MnOimyfRRn3V7jzyxo5djtCgP+8ws2v1fBiIoiMl67FoA71n+Hkjqa3975huHdY4g/kI/u9zHOzXr8WbECxo1zLpUhaNHUtlhffi3nPwecEcyFZFm+OED6XK/vc3BOxGuz6PR6UnZsYP+OnezcsZtxF00gRW+gcNdedF26sLhdOPL+Y3RzbRgEzm1Fx8z6ko23DGfjsuVEnXkTBtRqG/ucO+1lfvrg355wWIEAYG3u4WrKQacq/PHu3c7vrjTD+nVMykuD1KnOCLjYWMjKgl9+gQMHYNgwuPRS5+hB0Co40V9SjB8bAL3RQPdB/eg+qJ8nrXO/Xp7vpoR2DH9wEQl7cjDo9cz5z+VsvnIy+jmrSXWV8Tc9rY2X+OHd2zGKtZMELhRFZf9RM9M+SK+WF2c+RjfzUd8OrtPBlClVx8eOgdUKd93V4LIKmoYTVRBBrdgrqCOuDX/o0kXThqvX60l/7DxyCsowdY2mfEsm+jWrtU1KMQk8duc8Fj4xhZAAe0AI2h6KonLlgjWs21Wkmb/i0bPRf9UFDh+uGon6L9YXEgIDBjS0qIImRBgJ6xObDfsvv5G/6EeUyso6n16y9xBrplxLaVIyavfuqKed5lQW/igKxsICTJ0iObBFJmp41Rxr/13fIh+dzdfPXCaUg8CHwnIb6flHfd7wDA4bJ2ev4ck1HxEWHw+HDwM1mAmKi4WfoZUjTEz1hGKuQBoxAp2i0BPnUhVKuRljZHCzkYvz9xPbu7tn8T0dYF++Av2hw+gTvCJIFAX19NNh5UoqQsLpain3WT/J37S048wpJIlIEoEfnaJDGdmrA2l5RThUiC47yOb5N2m+MXraVVSUcz9oo9FpXhKTKVs9Nap/SZI+cP29JUCRz+tdohaI3WJjT98UdErVAnoGVPb8vCSo8xW7g/KBgwHfxfKy4npSGNXet+yhwzhWrETncBBpKfdoeP99GwCSZnzM2UO7n8itCVopOp2OL2aMZfVDZzLBWMJWl3LQeX18OHbMOWLYutXpdxDKoU1Q2whimCRJLwJTJUnq5Z0hy/LDsiwvaDjRWgb/3969B0dVWHEc/+4SDA+DUAORBlCeh0Sx1tCiCSgyvlD7cqYjouNo1ToWsVKnKL6gD1vrjK0gjrQqdXxQKo/ptDO1tR3B8hQawdYCBxBFoRIwgCIImmz6x92UJdyQgNm7D36fGWZ2783u/jKczdl7995zE3X1bCw7G9u2+ZBzFeqJ0WlU666GtXPzVnrs3XnIp/96YPIDz/Bik/lJtZ27sqm0jLO3rOGT9h3oWHeA9X3K2Dt+PD9r6MfuN1dz7ysvsO2x6Wy8bBgFOqJEmhGngR67anj2wbHA4eNcKC6G885j3aRJlHVJzl0644yQZ5J81dJfj28RnLB2BeDpj5N7dm7eSv/Na1M+9ffi/lHfpX74COZ2a931H07u24s1/YZQvunfAOyLF7B7aw0vlnQ77MS14qJCxk+awdvr3uG0QX2Yfkkfyvv1JhaPMz/RwPLVJZw7806d8CZHlkjABRfA4sWh41zYsAH69w8OksiC8RKSGS01iBnufqmZjXL34+6s5tY4uW8v1g44k4Eb/8X63oPY+fICppV0oeSkjq3+Ix2LxynzVWxbv4ltb2/lzIurKG3mcNRYLMasWyqp3Tv0sAv4xOMxunUsUHOQlu3YAUuXhh8EUVUF/frpLGhpsUF0MbM5wAgzO+T7Bncfm75YuSMWjzN4zT/ZuXkrsf17GGGntPygEPGCdvQsH0jP8oEt/2w8Rvemo7tFjkaPHlBZCUuWBF86HzhwcN3SpUEDKSlp/vFyXGipQYwGhgADgF+nP05uihe0o7h/H3ZoU1xyRSwGCxbAmjUwZMih6zp21DUbBGj5PIgid18MXAe83+SfiOSyeBxOPz2Ym5Tqvfe0e0mAlrcgfpD8NzVkXatmMIlIFovFguF6K1cGlwBduBC6dct0KskSLR7FZGbf5OARcJ8B7WlyNTgRyWHxOAwbBm+ETeGX41lLu5gGA+XAKwTXcRgEXAksSncwERHJrCM2CHc/4O77gf7uviK5bBVB45CjlUhATc3hQ89ERLJQa0+z3W1mPwFWAOcC76QtUb5qnKG0dClUVhJbsECDzkQkq7X2L9Q1wDaCw163ATekLVGe+v8Mpbo66hcvIVGzPdORRESOqFVbEO6+F3g8zVny2v9nKG1dy+ulZfTr3JXumQ4lInIEmuQWkUNmKJX1ZbbOhBaRLKcGEZEjzVASEclGkTUIM4sBW4ANyUXL3H1Sk5+ZDFwO1AF3NB45lS80Q0lEckmUWxD9gdfd/WthK83sbOB8YBjQG5gHfCW6eCIikirWENEx+WZ2FXAX8CHwCTDB3T1l/e1AJ3d/KHl/FXCxu+9o7jlXr17dUFhYyP79++nQoXWX9kynbMihDG2XYd++fdUVFRVD2yjSUcmm2laG7MnQVjlaW9tp2YIwsxuBCU0WjwN+7u5zzGw48DyHbiF0AWpT7u8BTgKabRCFhYWUlZWxdu1aysrK2ib855ANOZSh7TJUV1e3UZqjl021rQzZk6GtcrS2ttPSINz9aeDp1GVm1onguwXcfbGZlZpZzN0bN2E+AopSHlIE7E5HPhERaVmUp/JOBu4AMLMvAe+mNAeAJcAlZhY3sz5A3N0/iDCfiIikiPJL6oeA582s8Sil6wHM7GFgrruvMLNFwDKCxjUuwmwiItJEZA3C3XcRHMLadPnElNtTgClRZRIRkeZpWpyIiIRSgxARkVBqECIiEkoNQkREQqlBiIhIKDUIEREJpQYhIiKh1CBERCSUGoSIiIRSgxARkVD53SASCaipgYiueSEikk/yt0HU1dFQVUVDaSkNI0cGzUJERFotPxtEIkFi+AhYvpxYfT31ixaTqNme6VQiIjklLxtEomY7iZUriQENwBs9B1HbuWumY4mI5JS8bBC1nbtSXVrGZ/F2vP7FwTw46TcUFxVmOpaISE6J8oJBkSkuKmT8pBlsWvcOpYNOZd6tlcRisUzHEhHJKXnZIGKxGLNuqaR271CKTzxBzUFE5BjkZYMAiMdjdNduJRGRY5aX30GIiMjnF2vI4ZPIqqurdwCbM51D8tapFRUV3TPxwqptSbNW1XZONwgREUkf7WISEZFQahAiIhJKDUJEREKpQYiISCg1CBERCaUGISIioXL2TGozGwb8wt1HpiwbC4x393OjzmBmPYAngW5AO+A6d38r4gxnATOAOmA9cJO7p/VCGGbWHpgJnAYUAj8F1gDPEAzTfRMYl84czWR4F3gMqAcOEPx/1KQrQ1tSbYdmiLS2VdeBnNyCMLOJwFNAh5RlZwE3ApEMXgrJ8DDwgrufB9wHDM5AhsnAj919OEFBXZ7uDMC1QK27jwBGA9OBXwL3JZfFgG9kIMNUgj+oI4H5wF1pztAmVNvNZoi6tlXX5GiDAN4Crmy8Y2YnAw8Bd2QqA1AF9DKzvwPXAAszkGEV8AUziwFFwGcRZJgD3J9yvw6oAF5N3n8JuDADGca4++rk/QJgf5oztBXVdniGqGtbdU2ONgh3n0eyQMysHfA0MAHYk4kMSacBu9z9QoLNwLR/Yg3JsAGYBqwFSojgjezuH7v7HjMrAuYSfMKMuXvjKfp7gJOizuDu7wOYWSVwG/CrdGZoK6rtZjNEWtuq60BONogmKoCBwBPAbKDczB7NQI5a4I/J238ChmYgw1RghLsPBp4FHoniRc2sN7AAeM7dZwGp+2WLgN0ZyICZXUWw3/pyd9+R7gxpoNo+KPLaVl3nQYNw9xXufnpyn9wYYI27R7k53mgxcFny9nnAfzKQYSfwUfL2fwm+VEwrMysBXgbucveZycWrzGxk8vZoYFHUGczsWoJPWCPdfVM6Xz9dVNuHiLS2VdeBnD2KKQvdCTxlZrcCHwJjM5DhJmC2mdUBnwI3R/Ca9xC8We83s8b9pd8HppnZCQS7BOZGnKEdcAbBNNT5ZgbwqrtPTnOOfHU81rbqGk1zFRGRZuT8LiYREUkPNQgREQmlBiEiIqHUIEREJJQahIiIhNJhrlnMzO4mOJ0/QTAg7B5gPDDb3f9iZgXALOADgsPuLnP30SmPnwf8jeB0/N8C57j7a8l17YH3genuPiW5bBjBsd1V7r4y5Xm+DDwIdE0+1y7gdnffmr7fXvKZajs3aAsiS5lZOfB14CJ3v5hgvMHMlPXtgReBTe7+PYJBXgVmdmNy/RigvbvPSD5kHXB1yktcSnBMe6qbCM5QHZfyOj2BF4AJ7l7p7qOA5wgGuIkcNdV27lCDyF7bgT7Ad8ysNDmg66vJdYUEkxxXu/vdAMkZMTcQnFRTDtxLMAG00UvARWbW+H9+NfC7xpVmdiIwCvgRUGVmxclV1wFPubs3/qy7/4Fg0qTIsVBt5wg1iCzl7h8QfMqqApaZ2TrgiuTqaUBnoFeTx2wBHgCWARObzGn5NLn8/OTwry7AlpT1Y4D57r4f+D0H34B9gY0AZtbRzBaa2cLGZSJHS7WdO3QmdZYyswEA7t5YwEOBPwPLgX8AjxLsU33c3Z9v8tht7n5Kyv3rCWb4v0zwZllCcNr+CcAp7j7FzJYTjBP+GOhE8AYdAEwCPnb3qUd6DZHWUm3nDm1BZK8zgSfMrPGCKesJ9qvWA2+6ex3BbP6Hzayslc+5EDgH+DYpc2TMbAjQzt2Hu/ulyQvDvEXwqe5Z4GYzG5Ty8xXAiZ/nl5Pjmmo7R6hBZCl3n09Q9K+Z2RLgr8APSfnyLTnNcSIwx8w6teI5EwRHfuxz949SVt1M8OVcqieB29z9PYI36yPJTfDlBJc+vOhYfzc5vqm2c4d2MYmISChtQYiISCg1CBERCaUGISIiodQgREQklBqEiIiEUoMQEZFQahAiIhLqf65l0lc3Bnx+AAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "There is not a particular time of night where the variations are large. The two populations are very similar."
  },
  {
   "metadata": {
    "collapsed": true,
    "trusted": false
   },
   "cell_type": "code",
   "source": "df['time'] = df.MJD % 1\ndf.time = df.time.astype(np.float) * 24.0",
   "execution_count": 19,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "df.query('diff_des > 0.5 or diff_des < - 0.5').time.hist(histtype='step', lw=2, alpha=1, normed=1)\ndf.time.hist(histtype='step', lw=2, alpha=1, normed=1)",
   "execution_count": 22,
   "outputs": [
    {
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x127ad2550>"
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAAHcCAYAAAATCPhsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X9wlOW9//9XkiXZRYzRBaKiNaC12TaWJWm0VjJojSIz\n2sCxOIczYzN8TcFSK9ZRkJASFSQSdao9cbRyOttzgj9qU9Qz+AE5tJ+PHntaaQohqclqCaYVjSEr\nWYxmswvZ/f6Rk9VtAuSGJHvt7vMxg9n7uq/7zvv2zQ2v3Nx7b1okEokIAAAAMFR6vAsAAAAAToTA\nCgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYzXJgDYVC\nqqysVHFxsUpKSuTxeE66TWNjo0pLS4eN79ixQ/Pnz9ecOXN022236cMPP7RaDgAAAJKc5cC6adMm\ntba2qr6+XtXV1aqrq9POnTuPO/+dd97RXXfdpX/8BNg9e/bonnvuUUVFhV566SVNmjRJd999t/Uj\nAAAAQFKzFFgDgYAaGhpUVVWl/Px8lZaWqqKiQlu2bBlx/gsvvKAlS5Zo6tSpw9Z5PB6VlZVp8eLF\nysvLU1VVlbq7u+X3+0/tSAAAAJCULAVWr9ergYEBud3u6FhRUZGam5tHnP/mm2+qtrZW5eXlw9bt\n3r1b1113XXT5ggsu0G9/+1vl5ORYKQkAAABJzlJg7e7uVk5Ojmw2W3TM6XQqGAyqp6dn2Py6uroR\n713t7e3VkSNHdOzYMd12222aO3euVqxYoa6urlM4BAAAACQz28mnfC4QCCgzMzNmbGg5FAqNej99\nfX2SpIceekh33323Zs6cqccff1y33367XnrppVHt49ixYzpy5IiysrKUns7DDgAAAEwTDocVDAZ1\n1llnxVzwtMrSlllZWcOC6dCyw+EY9X4yMjIkSYsXL9ZNN90kSXr00Ud11VVXqampKeaWg+M5cuSI\nOjo6Rv09AQAAEB95eXlyOp2nvL2lwJqbmyu/369wOBy9qunz+WS325WdnT3q/Zx99tmy2WyaOXNm\ndCwnJ0c5OTnq7OwcVWDNysqSJE2dOlVTpkyxchhIQMFgUJ2dnTrvvPOivUfyot+phX6nFvqdWj79\n9FP5fL7T7rWlwOpyuWSz2dTU1KTCwkJJg89YLSgosPRNMzIyVFBQIK/XqwULFkiSDh8+rJ6eHs2Y\nMWNU+xgKzFOmTDmtxI7E0NfXp87OTuXk5Gjy5MnxLgfjjH6nFvqdWuh36vH5fKd9+6alre12u8rK\nylRdXa2Wlhbt2rVLHo8n+hQAn8+nYDA4qn0tXbpU9fX12rFjh9rb21VZWamvfvWr+vrXv279KAAA\nAJC0LN/9umbNGj3wwAMqLy/XmWeeqZUrV0afBDB37lw9/PDDWrhw4Un3M3/+fH3yySeqra1VT0+P\nrrjiCj355JPWjwAAAABJLS3yjx9BlSD6+vrU1tZ22jfxIjEM9dvlcvFPSCmAfqcW+p1a6Hdq+fjj\nj9XR0XHa/eZ5UAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAF\nAACA0QisAAAAMBqBFQAAAEYjsAIAAMBotngXcLoe+n2dukOH413GhPnJ1St16dRZ8S4DAABgwiR8\nYD0aPqrgQCjeZUyYcCQS7xIAAAAmVMIH1oy0DH3prBnxLmNc9QT86g19Fu8yAAAA4iLhA+s59hyt\nm//jeJcxrv6j6Tfa9s6ueJcBAAAQF7zpCgAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQC\nKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAA\njEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAF\nAACA0QisAAAAMBqBFQAAAEYjsAIAAMBolgNrKBRSZWWliouLVVJSIo/Hc9JtGhsbVVpaetz127dv\nV35+vtVSAAAAkAJsVjfYtGmTWltbVV9fr4MHD2r16tWaMWOGrr/++hHnv/POO7rrrruUlZU14vre\n3l499NBDSktLs1oKAAAAUoClK6yBQEANDQ2qqqpSfn6+SktLVVFRoS1btow4/4UXXtCSJUs0derU\n4+6ztrZWF110kbWqAQAAkDIsBVav16uBgQG53e7oWFFRkZqbm0ec/+abb6q2tlbl5eUjrt+9e7d2\n796t22+/3UoZAAAASCGWAmt3d7dycnJks31+J4HT6VQwGFRPT8+w+XV1dce9dzUUCmndunW6//77\nj3u7AAAAAGDpHtZAIKDMzMyYsaHlUChk6Rs/+eSTKigo0JVXXqndu3db2vaLIpGI+vr6Tnn7RHDs\n6NHo62CwP+mPdySBQCDmK5Ib/U4t9Du10O/UEgwGx2Q/lgJrVlbWsGA6tOxwOEa9n3fffVe//vWv\ntW3bNkmDofNUHT12VG1tbae8fSL4+PDh6OuOjg4d6+qPYzXx1dHREe8SMIHod2qh36mFfsMKS4E1\nNzdXfr9f4XBY6emDdxP4fD7Z7XZlZ2ePej87d+7UJ598omuvvVaSFA6HFYlEVFhYqAcffFA33njj\nqPc1yTZJLpfLymEknH1v/1XyD77Oy8vTl8+ZGd+C4iAQCKijo0N5eXmWfjhCYqLfqYV+pxb6nVr8\nfr86OztPez+WAqvL5ZLNZlNTU5MKCwslDT5jtaCgwNI3/d73vqeysrLoclNTk1atWqVXXnlFTqfT\n0r7S0tI0efJkS9skGtukSdHXWVn2pD/eE3E4HCl9/KmGfqcW+p1a6HdqGKtbPywFVrvdrrKyMlVX\nV2vjxo3q6uqSx+PRww8/LGnwauuZZ5550jdRZWdnx1yRHUreF154odX6AQAAkOQsf9LVmjVrVFBQ\noPLycq1fv14rV66MPglg7ty52r59+5gXCQAAgNRl+ZOu7Ha7ampqVFNTM2yd1+sdcZtFixZp0aJF\nx93n5ZdfnvRvnAIAAMCpsXyFFQAAAJhIBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDR\nCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAA\nADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0WzxLuB0BULH\n9F9v/S3eZYyrvx3+JPo6EDwax0oAAAAmXsIH1k8+O6qfvdgU7zLGle3CQ5p03uDrI5+G4lsMAADA\nBOOWAAAAABgt4a+wZqSn6Y7Fs+Ndxrj6dVuHPjn5NAAAgKSU8IE1LU2a/828eJcxrnb8LZPACgAA\nUha3BAAAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARkv4\nDw5Acml777B++vyeYePhSFihUEiZOz5Welry/Zw17WyHHvrBVfEuAwAAIxFYYZTQ0QF1fvzZCWYE\nJqyWCZUW7wIAADAXgRXGypyUIXtmhiQpEoloYGBAGRkZSktLnnT3yWeheJcAAIDxCKww1sJ5F+vW\nBS5JUl9fn9ra2uRyuTR58uQ4VzZ2/uUn29XbR2gFAOBEku9mQAAAACQVAisAAACMRmAFAACA0SwH\n1lAopMrKShUXF6ukpEQej+ek2zQ2Nqq0tHTY+DPPPKNrr71WRUVFWrp0qdrb262WAwAAgCRnObBu\n2rRJra2tqq+vV3V1terq6rRz587jzn/nnXd01113KRKJxIw///zz+uUvf6l169Zp69atmjFjhr7/\n/e8rGAxaPwoAAAAkLUuBNRAIqKGhQVVVVcrPz1dpaakqKiq0ZcuWEee/8MILWrJkiaZOnTps3csv\nv6zbbrtN8+bN00UXXaT7779fPT092rNn+EPjAQAAkLosBVav16uBgQG53e7oWFFRkZqbm0ec/+ab\nb6q2tlbl5eXD1q1evVo33nhjdHno2Zq9vb1WSgIAAECSsxRYu7u7lZOTI5vt88e3Op1OBYNB9fT0\nDJtfV1c34r2rklRYWKjc3Nzo8osvvqiBgQEVFRVZKQkAAABJztIHBwQCAWVmZsaMDS2HQqf+8PN9\n+/aptrZWFRUVcjqdlraNaPCh8sksHA5Hf7QIhkJJfbz9wf7o66NHj0aPNRAIxHxNHoP3dkfC4aTu\nq1XJ22+MhH6nFvqdWsbqvUmWAmtWVtawYDq07HA4TqmAvXv3atmyZZo3b57uvPNO6zuIRNTW1nZK\n3ztRBINB6X//93744YdqO5o8H036j/7+0eeB1efzqa3taMz6jo6OCa5ofB0bGJAkhUJHk/738alI\ntn7jxOh3aqHfsMJSYM3NzZXf71c4HFZ6+uAlP5/PJ7vdruzsbMvf/K233tLtt9+ukpISPfbYY5a3\nlySlpcnlcp3atgki690/Rl+ff/75cl2SvMd7LPNjST5J0tSpU+VyXSJp8Cfxjo4O5eXlnfIPRyay\nZXRJCiszc1LS/z62Iln7jZHR79RCv1OL3+9XZ2fnae/HUmB1uVyy2WxqampSYWGhpMFnrBYUFFj+\nxu+++65WrFihq6++Wo899lg0AFuVJiXVZ8uP5Iv/b7IyM5P6eO1Zn0VfT5o0adixOhyOJDv+wavl\naenpSXZcYyP5+o0Tod+phX6nhrG69cNSSrTb7SorK1N1dbVaWlq0a9cueTye6FMAfD7fqO9VWLdu\nnc4//3zdd999Onz4sHw+n6XtAQAAkBosX9Zcs2aNCgoKVF5ervXr12vlypXRJwHMnTtX27dvP+k+\nfD6f9u3bp/379+vqq69WSUlJ9NdotgcAAEDqsHRLgDR4lbWmpkY1NTXD1nm93hG3WbRokRYtWhRd\nnjp1Km8wAQAAwKhYDqyIr798/BeF9vvjXca4Oej7VBnT/i5J8h1zSOKNSAAApDoCa4J54+B/642D\n8a5ifGXOHPz6/tGwpOvjWgsAAIi/U3trPgAAADBBuMKaAKbpEr13YPATkW69wSVnjj3OFY2ftoOd\n+r8f/le8ywAAAAYhsCaAM9OmacA3+LGd3zi3SBedZ/1DGhJFpNdLYAUAADG4JQAAAABGI7ACAADA\naARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYA\nAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiN\nwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAA\nAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAo1kO\nrKFQSJWVlSouLlZJSYk8Hs9Jt2lsbFRpaemw8W3btum6666T2+3WHXfcoZ6eHqvlAAAAIMlZDqyb\nNm1Sa2ur6uvrVV1drbq6Ou3cufO489955x3dddddikQiMePNzc2qqqrSj370I7344os6cuSI1qxZ\nY/0IAAAAkNQsBdZAIKCGhgZVVVUpPz9fpaWlqqio0JYtW0ac/8ILL2jJkiWaOnXqsHXPPvusFixY\noO985zu69NJL9cgjj+j111/XBx98cGpHAgAAgKRkKbB6vV4NDAzI7XZHx4qKitTc3Dzi/DfffFO1\ntbUqLy8ftq6pqUnFxcXR5XPPPVfnnXee9u3bZ6UkAAAAJDlLgbW7u1s5OTmy2WzRMafTqWAwOOL9\np3V1dSPeuzq0r+nTp8eMTZ06VR999JGVkgAAAJDkbCef8rlAIKDMzMyYsaHlUChk6Rv39/ePuC+r\n+4lI6uvrs7RNojl27Fj0daC/X319ltqWUIKhYPR1OByO9jYQCMR8TR6D93ZHvnCsSOZ+YyT0O7XQ\n79QSDAZPPmkULCWfrKysYYFyaNnhcFj6xsfbl91ut7QfRSJqa2uztk2C8fv90dcHDhzQZ4cnxbGa\n8fVh56Ho62B/cFhvOzo6Jrii8XVsYECSFAodTfrfx6ci2fqNE6PfqYV+wwpLgTU3N1d+v1/hcFjp\n6YN3E/h8PtntdmVnZ1v6xtOnT5fP54sZ8/l8w24TOKm0NLlcLmvbJJg33mmV9JkkadasWfpS7pT4\nFjSOuiM2qX3wdZY9K9rbQCCgjo4O5eXlWf7hyGS2jC5JYWVmTkr638dWJGu/MTL6nVrod2rx+/3q\n7Ow87f1YCqwul0s2m01NTU0qLCyUNPiM1YKCAsvf2O12689//rMWLlwoSers7NRHH32k2bNnW9pP\nmqTJkydb/v6J5Iv3DDvs9qQ+3qzMrOjr9PT0YcfqcDiS7PjTBv87wrEiGfuNE6HfqYV+p4axuvXD\n0puu7Ha7ysrKVF1drZaWFu3atUsejyf6FACfzzfqexWWLFmiV155RQ0NDfJ6vVq9erWuueYazZgx\nw/pRAAAAIGlZ/uCANWvWqKCgQOXl5Vq/fr1WrlwZfRLA3LlztX379lHtx+1268EHH9STTz6pf/mX\nf1FOTo42btxotRwAAAAkOctvN7fb7aqpqVFNTc2wdV6vd8RtFi1apEWLFg0bX7hwYfSWAAAAAGAk\nlq+wAgAAABOJwAoAAACjJe8T6JHwPgv79X8P/I+kwWf0fvhJp3x/7x32gROJLJzzd2VMPqbglCwN\nhAeUkZ4R75IAADAOgRXG6hn4SE/9qT528NDIcxPWBVKmpICkgfD3CKwAAIyAwJpg9u3vVufHn8W7\njHGz/wP/yScBAICUQmBNMJtf/ku8SxhfGUeVcc7XJEnfcOXqmwXnShq8JaCz8yOdd965SXVLwDO/\nf0XKSt4fQAAAGAsEVphlYJIGui+UJF00+1KVXjz4caV9fX1q62uT6yJXUn0yyub/t1MRAisAACdE\nYE0AV152ns51Jk9IGy1X3jnxLgEAABiAwJoAvuHK1TdcufEuAwAAIC54DisAAACMRmAFAACA0Qis\nAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAw\nGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUA\nAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYj\nsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDTLgTUUCqmyslLFxcUqKSmRx+M57tzW1lbdcsstcrvd\nWrx4sd5+++2Y9f/6r/+qefPm6fLLL9ePf/xjHT582PoRAAAAIKlZDqybNm1Sa2ur6uvrVV1drbq6\nOu3cuXPYvEAgoGXLlqm4uFhbt26V2+3W8uXL1d/fL0l64YUXtHXrVj322GN67rnndOjQIf3kJz85\n/SMCAABAUrEUWAOBgBoaGlRVVaX8/HyVlpaqoqJCW7ZsGTb31VdflcPh0L333qtZs2Zp7dq1OuOM\nM7Rjxw5J0htvvKEFCxboG9/4hi655BJVVFToD3/4w9gcFQAAAJKGpcDq9Xo1MDAgt9sdHSsqKlJz\nc/Owuc3NzSoqKooZKyws1N69eyVJOTk5ev3119XV1aX+/n5t27ZNX/va107lGAAAAJDELAXW7u5u\n5eTkyGazRcecTqeCwaB6enpi5h46dEjTp0+PGXM6nerq6pIk/fCHP1R6errmzZunoqIi7dmzR48+\n+uipHgcAAACSlO3kUz4XCASUmZkZMza0HAqFYsb7+/tHnDs07+DBg5o8ebJ+/vOfKzs7W5s2bVJl\nZaV+8YtfWDqAiKS+vj5L2yDxBAKBmK/JIxJ91RcI6FjGsTjWYo7k7TdGQr9TC/1OLcFgcEz2Yymw\nZmVlDQumQ8sOh2NUc+12uyTpvvvu0+rVqzVv3jxJ0uOPP65rrrlGzc3N+vrXvz76oiIRtbW1WTkM\nJLCOjo54lzCmIp/nVXm9Xk1Kt3RKJr1k6zdOjH6nFvoNKyz97Zibmyu/369wOKz09MG7CXw+n+x2\nu7Kzs4fN7e7ujhnz+XyaNm2aDh8+rM7OTn3lK1+Jrjv33HN19tln68MPP7QWWNPS5HK5rBwGElAg\nEFBHR4fy8vKG/XCUyNKaXopeY83Pz1dmxqS41mOKZO03Rka/Uwv9Ti1+v1+dnZ2nvR9LgdXlcslm\ns6mpqUmFhYWSpMbGRhUUFAybO3v2bG3evDlmbM+ePVqxYoXOOussZWZmqr29XTNnzpQkHT58WH6/\nXxdccIGlA0iTNHnyZEvbIHE5HI4k63da9NVkh0OZtswTzE09yddvnAj9Ti30OzWM1a0flt50Zbfb\nVVZWpurqarW0tGjXrl3yeDwqLy+XNHgFdehehfnz56u3t1cbN25Ue3u7NmzYoEAgoBtuuEEZGRn6\np3/6J23atEmNjY169913tWrVKs2ZM2fE8AsAAIDUZfmDA9asWaOCggKVl5dr/fr1WrlypUpLSyVJ\nc+fO1fbt2yVJU6ZM0dNPP63GxkbdfPPNamlp0ebNm6P3sFZWVuq6667TPffco+9973s666yzVFdX\nN4aHBgAAgGRg+R0edrtdNTU1qqmpGbbO6/XGLF922WXaunXriPvJzMzUqlWrtGrVKqslAAAAIIVY\nvsIKAAAATCQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgB\nAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0\nAisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAA\nAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjGaL\ndwEABv1/G/5LaZGMeJcxIR5cfqVmnn9WvMsAACQIAitgiCOfBqUUCawD4Ui8SwAAJBACKxBHkyal\nK/S/r6efMzmpr7B+8llIgeCxeJcBAEhABFYgjr58YY7ePnRIkvTUqm8r05YZ54rGz1O/2af/8z8d\n8S4DAJCAeNMVAAAAjEZgBQAAgNEIrAAAADCa5cAaCoVUWVmp4uJilZSUyOPxHHdua2urbrnlFrnd\nbi1evFhvv/12zPodO3Zo/vz5mjNnjm677TZ9+OGH1o8AAAAASc1yYN20aZNaW1tVX1+v6upq1dXV\naefOncPmBQIBLVu2TMXFxdq6davcbreWL1+u/v5+SdKePXt0zz33qKKiQi+99JImTZqku++++/SP\nCAAAAEnFUmANBAJqaGhQVVWV8vPzVVpaqoqKCm3ZsmXY3FdffVUOh0P33nuvZs2apbVr1+qMM87Q\njh07JEkej0dlZWVavHix8vLyVFVVpe7ubvn9/rE5MgAAACQFS4HV6/VqYGBAbrc7OlZUVKTm5uZh\nc5ubm1VUVBQzVlhYqL1790qSdu/ereuuuy667oILLtBvf/tb5eTkWDoAAAAAJDdLgbW7u1s5OTmy\n2T5/fKvT6VQwGFRPT0/M3EOHDmn69OkxY06nU11dXert7dWRI0d07Ngx3XbbbZo7d65WrFihrq6u\n0zgUAAAAJCNLHxwQCASUmRn7YPOh5VAoFDPe398/4txQKKS+vj5J0kMPPaS7775bM2fO1OOPP67b\nb79dL730kqUDiEjR/SF5BQKBmK/JIjwQjr7uCwR0LCN5Pwnq2LHPj62/v/+E522y9hsjo9+phX6n\nlmAwOCb7sRRYs7KyhgXToWWHwzGquXa7XRkZgx8/uXjxYt10002SpEcffVRXXXWVmpqaYm45OKlI\nRG1tbVYOAwmso6Mj3iWMqc/6Pou+9nq9mpSevB8+98V/hXnvvfcUPHLyT/VKtn7jxOh3aqHfsMLS\n3465ubny+/0Kh8NKTx+8m8Dn88lutys7O3vY3O7u7pgxn8+nadOm6eyzz5bNZtPMmTOj63JycpST\nk6POzk5rgTUtTS6Xy8phIAEFAgF1dHQoLy9v2A9HieyMnv8n/e9Fhvz8fGVmTIprPePpf/a3SRoM\n6DNnztSs87OPOzdZ+42R0e/UQr9Ti9/vV2dn52nvx1JgdblcstlsampqUmFhoSSpsbFRBQUFw+bO\nnj1bmzdvjhnbs2ePVqxYoYyMDBUUFMjr9WrBggWSpMOHD6unp0czZsywdABpkiZPnmxpGyQuh8OR\nVP1Oz/j8NvLJDocybSe/6piovnjvu91uH1Ufk63fODH6nVrod2oYq1s/LL3pym63q6ysTNXV1Wpp\nadGuXbvk8XhUXl4uafAK6tC9CvPnz1dvb682btyo9vZ2bdiwQYFAQDfccIMkaenSpaqvr9eOHTvU\n3t6uyspKffWrX9XXv/71MTkwAAAAJAfLHxywZs0aFRQUqLy8XOvXr9fKlStVWloqSZo7d662b98u\nSZoyZYqefvppNTY26uabb1ZLS4s2b94su90uaTDQrlmzRrW1tfrud78rSXryySfH6rgAAACQJCy/\nw8Nut6umpkY1NTXD1nm93pjlyy67TFu3bj3uvhYvXqzFixdbLQEAAAApxPIVVgAAAGAiEVgBAABg\nNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisA\nAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxG\nYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNFu8CwAwaPcHTbKl\nJ+8peSj8N6WffUiS9HHAp0uUE+eKAACJInn/dgQSzM/+6Il3CeMu68uDX709F+gKXRLfYgAACYNb\nAgAAAGC4vm+IAAAVFklEQVQ0rrACcXTdxSWac15BvMuYEP/V3KyuyP54lwEASEAEViCOvvWlb8S7\nhAmzp+UTAisA4JRwSwAAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAA\noxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwmuXAGgqFVFlZqeLiYpWUlMjj8Rx3bmtrq2655Ra5\n3W4tXrxYb7/99ojztm/frvz8fKulAAAAIAVYDqybNm1Sa2ur6uvrVV1drbq6Ou3cuXPYvEAgoGXL\nlqm4uFhbt26V2+3W8uXL1d/fHzOvt7dXDz30kNLS0k79KAAAAJC0LAXWQCCghoYGVVVVKT8/X6Wl\npaqoqNCWLVuGzX311VflcDh07733atasWVq7dq3OOOMM7dixI2ZebW2tLrrootM7CgAAACQtS4HV\n6/VqYGBAbrc7OlZUVKTm5uZhc5ubm1VUVBQzVlhYqL1790aXd+/erd27d+v222+3WjcAAABShKXA\n2t3drZycHNlstuiY0+lUMBhUT09PzNxDhw5p+vTpMWNOp1NdXV2SBu+FXbdune6//35lZWWdav0A\nAABIcraTT/lcIBBQZmZmzNjQcigUihnv7+8fce7QvCeffFIFBQW68sortXv3bsuFD4lI6uvrO+Xt\nkRgCgUDMVySegfBA9PWxo0dPeN7S79RCv1ML/U4twWBwTPZjKbBmZWUNC6ZDyw6HY1Rz7Xa7/vrX\nv+rXv/61tm3bJkmKRCKWC4+KRNTW1nbq2yOhdHR0xLsEnKLPPv1MmjL4+uPDh0d13tLv1EK/Uwv9\nhhWWAmtubq78fr/C4bDS0wfvJvD5fLLb7crOzh42t7u7O2bM5/Np2rRpeu211/TJJ5/o2muvlSSF\nw2FFIhEVFhbqwQcf1I033jj6otLS5HK5rBwGElAgEFBHR4fy8vKG/XCExHDG/vboa+c555zwvKXf\nqYV+pxb6nVr8fr86OztPez+WAqvL5ZLNZlNTU5MKCwslSY2NjSooKBg2d/bs2dq8eXPM2J49e/SD\nH/xA1157rcrKyqLjTU1NWrVqlV555RU5nU5LB5AmafLkyZa2QeJyOBz0O0FlpGdI4cHXtkmTRtVH\n+p1a6Hdqod+pYaxu/bD0piu73a6ysjJVV1erpaVFu3btksfjUXl5uaTBK6hD9yrMnz9fvb292rhx\no9rb27VhwwYFAgEtWLBA2dnZuvDCC6O/cnNzJUkXXnghv3kBAAAQw/IHB6xZs0YFBQUqLy/X+vXr\ntXLlSpWWlkqS5s6dq+3bt0uSpkyZoqefflqNjY26+eab1dLSos2bN8tut4/tEQAAACCpWbolQBq8\nylpTU6Oampph67xeb8zyZZddpq1bt550n5dffjlvnAIAAMCILF9hBQAAACYSgRUAAABGI7ACAADA\naARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYA\nAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiN\nwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAA\nAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFY\nAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwmuXAGgqFVFlZqeLiYpWUlMjj8Rx3bmtrq2655Ra53W4t\nXrxYb7/9dsz6Z555Rtdee62Kioq0dOlStbe3Wz8CAAAAJDXLgXXTpk1qbW1VfX29qqurVVdXp507\ndw6bFwgEtGzZMhUXF2vr1q1yu91avny5+vv7JUnPP/+8fvnLX2rdunXaunWrZsyYoe9///sKBoOn\nf1QAAABIGpYCayAQUENDg6qqqpSfn6/S0lJVVFRoy5Ytw+a++uqrcjgcuvfeezVr1iytXbtWZ5xx\nhnbs2CFJevnll3Xbbbdp3rx5uuiii3T//ferp6dHe/bsGZsjAwAAQFKwFFi9Xq8GBgbkdrujY0VF\nRWpubh42t7m5WUVFRTFjhYWF2rt3ryRp9erVuvHGG6Pr0tLSJEm9vb1WSgIAAECSsxRYu7u7lZOT\nI5vNFh1zOp0KBoPq6emJmXvo0CFNnz49ZszpdKqrq0vSYHjNzc2NrnvxxRc1MDAwLOQCAAAgtdlO\nPuVzgUBAmZmZMWNDy6FQKGa8v79/xLn/OE+S9u3bp9raWlVUVMjpdFopSRFJfX19lrZB4gkEAjFf\nkXgGwgPR18eOHj3heUu/Uwv9Ti30O7WM1XuTLAXWrKysYYFzaNnhcIxqrt1ujxnbu3evli1bpnnz\n5unOO++0Us6gSERtbW3Wt0NC6ujoiHcJOEWfffqZNGXw9ceHD4/qvKXfqYV+pxb6DSssBdbc3Fz5\n/X6Fw2Glpw/eTeDz+WS325WdnT1sbnd3d8yYz+fTtGnTostvvfWWbr/9dpWUlOixxx47tSNIS5PL\n5Tq1bZEwAoGAOjo6lJeXN+yHIySGM/Z//tg65znnnPC8pd+phX6nFvqdWvx+vzo7O097P5YCq8vl\nks1mU1NTkwoLCyVJjY2NKigoGDZ39uzZ2rx5c8zYnj179IMf/ECS9O6772rFihW6+uqr9dhjj0UD\nsFVpkiZPnnxK2yLxOBwO+p2gMtIzpPDga9ukSaPqI/1OLfQ7tdDv1DBWt35YSol2u11lZWWqrq5W\nS0uLdu3aJY/Ho/LyckmDV1CH7lWYP3++ent7tXHjRrW3t2vDhg0KBAJasGCBJGndunU6//zzdd99\n9+nw4cPy+Xwx2wMAAADSKXxwwJo1a1RQUKDy8nKtX79eK1euVGlpqSRp7ty52r59uyRpypQpevrp\np9XY2Kibb75ZLS0t2rx5s+x2u3w+n/bt26f9+/fr6quvVklJSfTX0PYAAACAZPGWAGnwKmtNTY1q\namqGrfN6vTHLl112mbZu3Tps3tSpU3mjFAAAAEbl1G4cBQAAACYIgRUAAABGI7ACAADAaARWAAAA\nGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAK\nAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACj\nEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEA\nAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBotngXACD1/O5Pf9efXv/dcdeHwxEF\ng/3K2uVXenraBFY29kov/5IWzrsk3mUAQEIjsAKYcH3Z7+gzx3vxLmNC/LZ7mhbq3niXAQAJjcAK\nYEJkZKRJ4cHXaZlBpWUG41vQBPk0zB+zAHC6+JMUwIS49hsXyfvW6P7IiUiKRCJKS0tTIt4QEI5I\nA5Fj8S4DAJIGgRXAhPjWl4r0rS8VjWpuX1+f2tra5HK5NHny5HGubOy9feBj3f+H+5VmI7QCwFjg\nKQEAAAAwmuXAGgqFVFlZqeLiYpWUlMjj8Rx3bmtrq2655Ra53W4tXrxYb7/9dsz6bdu26brrrpPb\n7dYdd9yhnp4e60cAAACApGY5sG7atEmtra2qr69XdXW16urqtHPnzmHzAoGAli1bpuLiYm3dulVu\nt1vLly9Xf3+/JKm5uVlVVVX60Y9+pBdffFFHjhzRmjVrTv+IAAAAkFQs3cMaCATU0NCgX/ziF8rP\nz1d+fr4qKiq0ZcsWXX/99TFzX331VTkcDt177+DjXNauXas33nhDO3bs0MKFC/Xss89qwYIF+s53\nviNJeuSRR3TNNdfogw8+0IwZM8bo8AAAAE6utr5R/t7UeHrJF1WUFWjWjLPiXcZJWQqsXq9XAwMD\ncrvd0bGioiL9/Oc/Hza3ublZRUWxb7AoLCzU3r17tXDhQjU1NWn58uXRdeeee67OO+887du3j8AK\nAAAmlPdvPnUf6Yt3GRPO3zdLUpIF1u7ubuXk5Mhm+3wzp9OpYDConp4enX322dHxQ4cO6dJLL43Z\n3ul0av/+/dF9TZ8+PWb91KlT9dFHH1k+CAAwVTgyoL9+1BnvMozW39+vD3qPKPNQl+x2e7zLwTgz\ntd/Bs9rl+PJf4l3GhNvculdPfXl9vMs4Kcu3BGRmZsaMDS2HQqGY8f7+/hHnDs072fqTCYcHn0B+\nduaZ+vjjj0d/EEhIweDgP9P4/X4FAoE4V4Pxluj97u/rVW7WOUqzDUiSfvbHJ+NcUYJoincBmFCG\n9fucCyXJGe8yJly2LXtcc9Snn34q6fPcdqosBdasrKxhgXJo2eFwjGru0E9TJ1t/MkN/od10wZXq\n6OgY9TEgsXV2cqUqlSRyv5fO/E68SwCAUZmIHBUMBjVlypRT3t5SYM3NzZXf71c4HFZ6+uADBnw+\nn+x2u7Kzs4fN7e7ujhnz+XyaNm2aJGn69Ony+XzD1v/jbQLHc9ZZZykvL09ZWVnRWgAAAGCOcDis\nYDCos846vftkLQVWl8slm82mpqYmFRYWSpIaGxtVUFAwbO7s2bO1efPmmLE9e/ZoxYoVkiS3260/\n//nPWrhwoaTBKykfffSRZs+ePbrCbTY5nal36R4AACCRnM6V1SGWLk3a7XaVlZWpurpaLS0t2rVr\nlzwej8rLyyUNXiEd+qf6+fPnq7e3Vxs3blR7e7s2bNigQCCgG264QZK0ZMkSvfLKK2poaJDX69Xq\n1at1zTXX8IQAAAAAxEiLRCIRKxv09/frgQce0GuvvaYzzzxTFRUVuvXWWyVJ+fn5evjhh6NXTVta\nWlRdXa0DBw7oK1/5ih544AHl5+dH9/Xyyy/riSee0JEjRzR37lytX7/+tC8ZAwAAILlYDqwAAADA\nROLdSgAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMlZGANhUKqrKxUcXGxSkpK5PF44l0S\nxtGuXbuUn58vl8sV/bpy5cp4l4UxFgqFdNNNN+lPf/pTdOzgwYNaunSp5syZoxtvvFG///3v41gh\nxtJI/d6wYcOwc/3ZZ5+NY5U4XV1dXbrzzjt1xRVXaN68eXr44YejH8vO+Z18TtTv0z2/LX3SlSk2\nbdqk1tZW1dfX6+DBg1q9erVmzJih66+/Pt6lYRzs379f3/72t7VhwwYNPYUtKysrzlVhLIVCId19\n993av39/zPgPf/hD5efn6ze/+Y127dqlO+64Q9u3b9e5554bp0oxFo7X7wMHDuiee+7RokWLomNj\n8Qk5iJ8777xTOTk5eu655+T3+1VZWamMjAzde++9WrFihVwuF+d3EjlRv0/3/E64K6yBQEANDQ2q\nqqpSfn6+SktLVVFRoS1btsS7NIyT9vZ2ffnLX9Y555wjp9Mpp9PJX2JJpL29XbfccosOHjwYM/6H\nP/xB77//vh588EHNmjVLy5Ytk9vtVkNDQ5wqxVg4Xr+H1n31q1+NnudOp5MfThPYgQMH1NzcrJqa\nGl188cUqKirSnXfeqW3btumPf/yjDh48yPmdRE7Ub+n0z++EC6xer1cDAwNyu93RsaKiIjU3N8ex\nKoyn9vZ2zZw5M95lYJzs3r1bV155pX71q1/pi59j0tzcrK997Wsxf6AVFRWpqakpHmVijByv359+\n+qm6urqUl5cXv+IwpqZNm6Z/+7d/0znnnBMz3tvbq3379nF+J5mR+h2JRNTb2zsm53fC3RLQ3d2t\nnJwc2Wyfl+50OhUMBtXT06Ozzz47jtVhPLz33nv67//+bz311FMKh8O64YYbdOedd2rSpEnxLg1j\nYMmSJSOOd3d3a/r06TFjTqdTXV1dE1EWxsnx+n3gwAGlpaXpqaee0htvvKGcnBwtXbo0+lHfSDxn\nnnmmrrrqquhyJBLRli1bdOWVV3J+J6Hj9ftb3/rWmJzfCRdYA4GAMjMzY8aGlodu7EXy+PDDD9Xf\n36+srCw98cQTOnjwoDZs2KBgMKjKysp4l4dxdLxznfM8OR04cEDp6em6+OKLdeutt2r37t36yU9+\noilTpqi0tDTe5WEM1NbWqq2tTQ0NDfJ4PJzfSa62tlZer1cNDQ36y1/+ctrnd8IF1qysrGG/oYeW\nHQ5HPErCODr//PP11ltvKTs7W5KUn5+vcDisVatWac2aNUpLS4tzhRgvWVlZOnLkSMxYKBSS3W6P\nU0UYTwsXLtS3v/3t6Ll+6aWXqqOjQ88//zyBNQk88sgjqq+v1+OPP65LLrmE8zvJ/WO/L7nkktM+\nvxPuHtbc3Fz5/X6Fw+HomM/nk91uj/6PQHL5x75efPHFCgaD8vv9caoIEyE3N1fd3d0xYz6fT9Om\nTYtTRRhv/3iuz5o1S4cOHYpTNRgr69ev17//+7/rkUceiYYTzu/kNVK/pdM/vxMusLpcLtlstpgb\nsxsbG1VQUBDHqjBe3nzzTV1xxRUKBoPRsdbWVuXk5HC/cpKbPXu2WltbY/5F5c9//nPMGy6RPH72\ns59p6dKlMWNtbW284TLB1dXV6Ve/+pV++tOfasGCBdFxzu/kdLx+j8X5nXCB1W63q6ysTNXV1Wpp\nadGuXbvk8XhUXl4e79IwDubMmSOHw6G1a9fqvffe0+uvv65HHnlE3//+9+NdGsbZ5ZdfrvPOO0/3\n3Xef9u/fr2eeeUYtLS367ne/G+/SMA6uueYa/elPf5LH49H777+v5557Tv/5n/+pioqKeJeGU9Te\n3q6nnnpKy5Yt05w5c+Tz+aK/OL+Tz4n6PRbnd1rki88VSRD9/f164IEH9Nprr+nMM89URUWFbr31\n1niXhXHS3t6ujRs3qqmpSWeccYb++Z//WStWrIh3WRgHLpdL//Ef/6Hi4mJJ0vvvv6/Kyko1Nzfr\nS1/6ktauXatvfvObca4SY+Uf+/273/1OTzzxhP72t79pxowZ+vGPf8z9qwnsmWee0U9/+tOYsUgk\norS0NLW1tenvf/+71q5dy/mdJE7W79M9vxMysAIAACB1JNwtAQAAAEgtBFYAAAAYjcAKAAAAoxFY\nAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACj/f82zQHv\noRmqEwAAAABJRU5ErkJggg==\n",
      "text/plain": "<matplotlib.figure.Figure at 0x127adf5d0>"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ]
  },
  {
   "metadata": {
    "collapsed": true,
    "trusted": false
   },
   "cell_type": "code",
   "source": "",
   "execution_count": null,
   "outputs": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python",
   "version": "3.6.6",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}