Gen_inp_HSC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "declared-woman",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from astropy.table import Table\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "earned-confidentiality",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_dir='./temp/hsc_test/' #directory where files for corr_pc will be placed. I recommend using full path rather than relative path as here.\n",
    "file_name='hsc_test'  #prefix to be appended to the files being saved. Should be different for different datasets.\n",
    "\n",
    "os.makedirs(file_dir,exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e1664dd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.cosmology import Planck15 as cosmo\n",
    "cosmo_h=cosmo.clone(H0=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a58afb5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#this is cosmology table, passed inside inp file with distances keyword. order of columns matter\n",
    "z=np.linspace(0,5,5001)\n",
    "H_z=cosmo_h.H(z).value\n",
    "D_A=cosmo_h.angular_diameter_distance(z).value\n",
    "D_C=cosmo_h.comoving_transverse_distance(z).value\n",
    "D_M=cosmo_h.luminosity_distance(z).value\n",
    "D_Z = (z,H_z,D_A,D_C,D_M)\n",
    "D_Z=np.column_stack(D_Z)\n",
    "np.savetxt('./temp/Distances_planck15.dat', D_Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "plain-ecuador",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" this function generates the input file for corr_pc\n",
    " corr_pc takes two samples to compute cross correlations. \n",
    " Shape sample: Sample of galaxies with shapes. For density correlations it is not necessary to pass shapes. \n",
    "               Prefix of filenames for this sample should be passed vias file_nameS\n",
    " Density sample: Sample of galaxies/particles that trace density field. For density correlations it is not necessary to pass shapes. \n",
    "                 If shapes are passed for this sample, w++ can be computed. Prefix of filenames for this sample should be passed vias file_nameD.\n",
    "                 If auto correlation is to be computed, simply set file_nameD=file_nameS.\n",
    " file_name: This is the prefix corr_pc will use to save the output files containing two point correlations. Input file to be passed to corr_pc is \n",
    "            also saved with this name. input file is saved in file_dir and the outputs are saved in new directory, file_dir+corr_data_out.\n",
    " corr: two point correlation to be computed. This should be density-density (wgg,wmm), shape-density (wg+) or shape-shape(w++). \n",
    "        When computing wg+, corr_pc will also compute and save wgg by default since pair counting is done anyways. \n",
    "        When doing w++, wg+ and wgg are also computed and saved.\n",
    " n_shape, n_density: number of galaxies in the shape sample and density sample respectively.\n",
    " do_jk: Whether to do jackknifing or not.\n",
    " njk: number of jackknife regions. jackknife regions are assumed to be labelled from 0...njk-1.\n",
    "\"\"\"\n",
    "def inp_gen(file_dir,file_name,file_nameS,file_nameD,corr=1,n_shape=0,n_density=0,njk=100,do_jk=1):\n",
    "    s=\"\"\n",
    "    s=s+\"which_corr     \"+str(corr)+\"\\n\" #0: density-density, 1: ia or WL: shape-density, 2 ia:shape-shape\n",
    "    s=s+\"coordinates    3\"+\"\\n\"  #0: rp-pi, 3: rp\n",
    "    s=s+\"estimator      4\"+\"\\n\"  #0:landy-szalay (LS): cross. works for both auto and cross correlations. For auto, pass same filenames for shape and density samples \n",
    "                                 # 1: LS:auto, for auto correlation only\n",
    "                                # 4: LS like, for mean quantities, i.e. for SD/RR like terms, we instead have SD/SD like terms.\n",
    "    s=s+\"data_sorted     1\"+\"\\n\"  #0: no sorting, 1:sorted by z, 2:sorted by dec (y)\n",
    "    s=s+\"use_comoving    1\"+\"\\n\" #use comving distances from the table\n",
    "    s=s+\"do_jk           \"+str(do_jk)+\"\\n\" #0: no jackknife, 1: do jackknife. jk regions must be passed in a file. are assumed to be from 0..njk-1\n",
    "    s=s+\"sig_crit        0\"+\"\\n\" #for weak lensing only. 0: no sig crit; 1: calc sig crit for each pair, 2: use g1 wt as sig crit\n",
    "\n",
    "    #following passes information about shape sample\n",
    "    s=s+\"shape_pos       \"+file_dir+file_nameS+'_pos.dat'+\"\\n\"  # x, y coordinates\n",
    "    s=s+\"shape_z         \"+file_dir+file_nameS+'_z.dat'+\"\\n\"    # z coordinate\n",
    "    s=s+\"shape_e         \"+file_dir+file_nameS+'_e.dat'+\"\\n\"    # e1, e2\n",
    "#     s=s+\"shape_wt         \"+file_dir+file_nameS+'_wt.dat'+\"\\n\"  # weights on galaxies. If no weighting, pass 0\n",
    "    s=s+\"shape_wt         0\"+\"\\n\"\n",
    "    s=s+\"shape_jk        0\"+\"\\n\"   # jackknife regions\n",
    "\n",
    "    # following passes information about density sample\n",
    "    s=s+\"density_patch   0\"+\"\\n\" #for weak lensing, not important for ia, gg corrs.\n",
    "    s=s+\"density_pos     \"+file_dir+file_nameD+'_pos.dat'+\"\\n\"  # x, y coordinates\n",
    "    s=s+\"density_z       \"+file_dir+file_nameD+'_z.dat'+\"\\n\"    # z coordinate\n",
    "    s=s+\"density_wt       0\"+\"\\n\"  # weights on galaxies/particles. If no weighting, pass 0\n",
    "#     s=s+\"density_wt       \"+file_dir+file_nameD+'_wt.dat'+\"\\n\"  # weights on galaxies/particles. If no weighting, pass 0\n",
    "    s=s+\"density_jk      0\"+\"\\n\"   # jackknife regions    \n",
    "    s=s+\"density_e      0\"+\"\\n\"   #if doing shape-shape, this file should be passed. # e1, e2\n",
    "\n",
    "    # following passes information about randoms corresponding to shape sample\n",
    "    s=s+\"Srandom_pos      0\"+\"\\n\"\n",
    "    s=s+\"Srandom_z        0\"+\"\\n\"\n",
    "    s=s+\"Srandoms_wt      0\"+\"\\n\"\n",
    "    s=s+\"Srand_jk         0\"+\"\\n\"\n",
    "\n",
    "    # following passes information about randoms corresponding to density sample\n",
    "    s=s+\"Drandom_patch    0\"+\"\\n\"\n",
    "    s=s+\"Drandom_pos      0\"+\"\\n\"\n",
    "    s=s+\"Drandom_z        0\"+\"\\n\"\n",
    "    s=s+\"Drandoms_wt      0\"+\"\\n\"\n",
    "    s=s+\"drand_jk         0\"+\"\\n\"\n",
    "\n",
    "    s=s+\"distances        ./temp/Distances_planck15.dat\"+\"\\n\"  #not needed for PB\n",
    "    s=s+\"patch_file       0\"+\"\\n\"  #not needed for PB\n",
    "    s=s+\"out_file         \"+file_dir+\"corr_data_out/\"+file_name+\"_\"+\"\\n\"  # prefix for the output files.\n",
    "    \n",
    "    os.makedirs(file_dir+\"corr_data_out/\",exist_ok=True) #directory where outputs are saved.\n",
    "    \n",
    "    s=s+\"n_threads       0\"+\"\\n\"   # number of openmp threads to use. 0== use maximum available.\n",
    "    s=s+\"n_shape         \"+str(n_shape)+\"\\n\"  #number of galaxies in shape sample\n",
    "    s=s+\"n_density       \"+str(n_density)+\"\\n\"  #number of galaxies in density sample\n",
    "    s=s+\"n_Srand         0\"+\"\\n\"\n",
    "    s=s+\"n_Drand         0\"+\"\\n\"\n",
    "    s=s+\"rand_subsample  0\"+\"\\n\"  #subsample randoms to speed up calculations.\n",
    "    s=s+\"n_jk          \"+str(njk)+\"\\n\"  #number of jk regions. regions are assumed to be labelled from 0....njk-1\n",
    "    s=s+\"n_patch         0 \"+\"\\n\"#not needed for PB/HSC\n",
    "\n",
    "    s=s+\"binR_min        0.1\"+\"\\n\"  # rp_min for the bins\n",
    "    s=s+\"binR_max       10\"+\"\\n\" # rp_max for the bins\n",
    "    s=s+\"n_bins         10\"+\"\\n\" #number of rp bins to use\n",
    "    s=s+\"lin_bin         0\"+\"\\n\"  #linear or log bins. 0: log, 1:linear\n",
    "    s=s+\"n_p_bin         1\"+\"\\n\"  #number of pi bins to use.. for WL leave to 1.\n",
    "    s=s+\"pmin         0\"+\"\\n\"  # pi_min (line of sight or z separation)\n",
    "    s=s+\"pmax         0\"+\"\\n\"   # pi_max (line of sight or z separation)\n",
    "    s=s+\"z_min         0\"+\"\\n\"    \n",
    "    s=s+\"z_max         4\"+\"\\n\" \n",
    "    s=s+\"dz            0.001\"+\"\\n\" #resolution of cosmology table.\n",
    "    s=s+\"z_sep_min     -100\"+\"\\n\"  \n",
    "    s=s+\"z_sep_max     100\"+\"\\n\"   \n",
    "    s=s+\"periodic_box  0\"+\"\\n\"\n",
    "    s=s+\"box_size     0\"+\"\\n\"  #size of the box.\n",
    "\n",
    "    fil=open(file_dir+file_name+'.inp','w')\n",
    "    fil.write(s)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "focal-concert",
   "metadata": {},
   "outputs": [],
   "source": [
    "hsc_dir='/hildafs/projects/phy200017p/ztq1996/ztq1996/hsc/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "environmental-least",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_dat_name='r025_rot00_zbin2.fits'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "940adcbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "66c44000",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[( 8279,     2, 36411452835235108,  30.62002, -6.376178,  0.5943237 , -0.07674798, 0.00756096,  0.00193293, -0.00612652, 4.789851, 0.714 )\n",
      " ( 8279,     3, 36411452835235125,  30.61326, -6.375836,  0.89743423,  0.3828571 , 0.00789049, -0.00303394, -0.00149775, 3.258304, 0.3827)\n",
      " ( 8279,     4, 36411452835235146,  30.62137, -6.37493 ,  0.13051362,  0.80013806, 0.00694159,  0.0030251 , -0.00234604, 5.383868, 0.6425)\n",
      " ...\n",
      " (16011, 32831, 70418247970748448, 247.1476 , 44.33325 ,  0.7923786 ,  0.10393251, 0.00550792,  0.00336679,  0.00675744, 5.18838 , 0.7885)\n",
      " (16011, 32834, 70418247970748684, 247.1574 , 44.338   , -0.6247564 , -0.24049284, 0.01167254,  0.00215797,  0.00726533, 5.601186, 0.8664)\n",
      " (16011, 32835, 70418247970748695, 247.1489 , 44.33728 ,  0.32006502, -0.32866257, 0.00719187,  0.00081219,  0.00230956, 5.949475, 0.6425)]\n"
     ]
    }
   ],
   "source": [
    "with fits.open(hsc_dir+S_dat_name) as f:\n",
    "    data=f[1].data\n",
    "    print(data)\n",
    "    S_dat=Table(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8e7d8f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_dat['z']=S_dat['z_source_mock']\n",
    "S_dat['ra']=S_dat['ra_mock']\n",
    "S_dat['dec']=S_dat['dec_mock']\n",
    "S_dat['e1']=S_dat['e1_mock']\n",
    "S_dat['e2']=S_dat['e2_mock']\n",
    "\n",
    "S_dat.sort('z')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "300006b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_dat=S_dat[:int(1e4)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ffdec25e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABPwAAAO1CAYAAAD6zfvcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACSz0lEQVR4nOz9fZhV5X0v/n8YmQGUxMHIAElkMAgkTXUUmVqpyTSpE6/JaE7TOs1TT2qskZgSevg1sW1CPI16mpRjQ+shD0ZJ0hx7mRPSJ4+URhINfg05ZnxCbaooUYiKDkbGVGEYkPn9MbM2a++99jAzDMzsxet1XV5h33utey/CvvbDe3/u+zOhr6+vLwAAAACAXKgZ6wsAAAAAAEaPwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHJo71BVBu//798fjjjxeNnXTSSVFTI58FAAAAONYcOHAgXnzxxaKxefPmxcSJ2dGewG8cevzxx+NXfuVXxvoyAAAAABinfvrTn8Zb3vKWzPuUjAEAAABAjgj8AAAAACBHBH4AAAAAkCP28BuHTjrppLKxn/70p3HyySePwdUAAADjyXs+8L/i1ls+Eb/56a/GD//yY4Xxixv+ML7btSbznOS+93zgf0VERN0PHoyIiN7fOjMiIm695ROZ52eNJXPcessnRuOvA8AQvPDCC2X9HrLyo4TAbxzK6sZ78sknx/Tp08fgagAAgPGipX1lTL3jP+L8a26Of7/xs9Fa01G4b2PfrRXP29h3a7TWdMTUiNhwYG20tK+MiIgfr7syWms64vxrbo7NA+fPv3ZVRERsWbG8MGdrTUdsOLC2MEdE+H4CMMay8qPCfUfxOgAAABih1pqO2LjuyoiIaFi9KSL6w7sNB9YW7k+0tK+M1pqOwlgS2JUel/y5fuu+wnmNV22Kxqs2FULB5HEAqB4CPwAAgCrQ29ZcFNZFRCHUKw3zkmAwPZY+rm59Z9St7yzc3rjuymhpX1kY23BgbdSt76x4HVnXAsD4YUkvAABAFdj+gVfjtPVRWFqbDvAS6SCu9L6k6i85P621piOirbnovuTPpWMAjH8CPwAAgCqx7erFERHRtXRx1G/dVxTcpcO40mW7Gw6sLar6K60IzAr1suYurSYUAAKMTwI/AACAKvDm2c/F/g8/G7FieTSs3pQZ8CUhXKUqvqyArnTsUEFfsrdf3Sj8nQA4MuzhBwAAUAW2vTgtIqKomUa6MUd6LOLgXntp6XNLj0/uzzovCfu2Xb04Nq67Mjauu1J1H8A4JvADAACoAn/b9H+it625qNlGomvp4mhpX1k0ngRzEcXLeistA+5ta46N664sNPRI9gNMH79lxfLCn7PCQwDGB0t6AQAAqsDP972u0E23taajP6CrsJ9eadVf+r7S/fqS23XRH+JtLJmnNFxsWrYqIiIaKnTxBWDsCfwAAACqwK9PebLodlYTjpb2lYW99bK69EZE1K3vjK6l/c0/Svf8S89Zyebr+6v8WldvGtHfA4Ajz5JeAACAKpK1pDdiIOyrcF9yXrJUt37rvqjfuq8oKCxt9FHakTc9rkMvwPimwg8AAGAcO/Xmz0dExL3vnFBUwVeqbmCJbRLGpZftpv83K8iri+xlwFlLh5OGHk3LVsXm65cL/wDGIYEfAADAOPbtt30tIiI27pke2z/wasyO5orHZgV9EcXLeUsbdXTPrc0M7pJ5kn0Dk+W+3XNrIyKiYfWmiOuXC/sAxiGBHwAAwDj23x57X0RE1E/eE32/rI269fdkHlca1qWr+DYOBIGllX0bS/b+K+3gmw4KE/bwAxj/BH4AAADjWP3kPRERcdrUF6L7tClDOmewJb2lkg69Ef0Vf+mlwekGIOnlvIPNB8DYE/gBAACMY1vvmtP/h7dH1P3tSRGxteKx6VCudM++5L6sZb1JyJfu0psEhulOvqWPJfQDGJ8EfgAAAOPYvjk9ERHx6PaZUXtubTSuzz6udMluurJvsI67ddEf+m1cd2XZkt7BQj1hH8D4JfADAAAYx6ZM3RsREfsfOTH+7Pf+IS5dcXfmcVlhX+l9iazGHpVuV+r6m54rCQwBGB9qxvoCAAAAqKz2jhOj9o4TY8uK5XHLm2cVVeGltdZ0lIVzWZLluUOp0Ctd1pt+nPR8Oz+6e9DHBODoEvgBAABUicGCuuS+rAq/pmWrIiIOGcpl3d/b1nzI5b3Tbzx+KJcPwFFiSS8AAMA41rB6U0REtK7eNKSls1mhXMPqTRHXL4+upYsPeW4S+vW2NUdd9DfyaGlfWdSxt1T33NrYbE8/gHFjQl9fX99YXwTFdu7cGQ0NDUVjXV1dMX369DG6IgAAYLxraV857H30Sqv2KjX3GGx/wPSxGnkAHBnDzYos6QUAABjHWtpXRkv7yoiovCS3taZjRE0zKu3Ll9yXDvBa2lfq2AtQJVT4jUMq/AAAgJFKqvCSkPBQQWBL+8qoW99ZMbRTvQcw9oabFdnDDwAAYBzLWlo7mGQfvo0Vluqm59xwYG3Ure+s+FiDNeoAYPwS+AEAAIxjIwncSs+ptE9f+r6kMUfWHn5Z5wEwfgn8AAAAxqnhVvcNZ75E0uwja+lvb1tztNZ09HcHFvYBVA1NOwAAAMappHFGpaq7oUov4S0NDtNLetN625oLIeBIGoIAMHZU+AEAAFS5dKCXtfQ2a/++ZLzSfUlFn2W8ANVH4AcAAFAlKoVvWXvyJQpLctddOeg+fPboA8gPgR8AAEAVyArkWtpXFv5cadntYGFg1nHpakEAqpPADwAAoEqNdG+9wUI9QR9A9RP4AQAAVIHRDOKEegD5pksvAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5MnGsL2A09fX1xSOPPBIPPfRQ7NixI3p6emLq1KnR2NgYZ599dsyePXusLxEAAAAAjqijFvjdc889sXjx4jhw4EDZfU8++WTMmTNnxHN3d3fHF7/4xVizZk08++yzFY9ramqKT3ziE3HJJZfEcccdN+LHAwAAAIDx6qgs6d27d29ceumlmWHf4Vq3bl3Mnz8/rrnmmkHDvoiIzZs3x2WXXRbNzc3xxBNPjPq1AAAAAMBYOyqB3+c+97n46U9/OurzrlmzJi666KLYuXPnsM574IEH4pxzzokHH3xw1K8JAAAAAMbSEV/Se99998X//J//c9Tnvf322+Pyyy+Pvr6+ovFZs2bFkiVLYtGiRTFjxox4+umn4+677441a9ZEd3d34bgXX3wx3v3ud8f9998fM2fOHPXrAwAAAICxMKGvNDEbRfv27Yuzzz47Hn744YiIqK2tjaamprj33nuLjhvuHn4vvfRSLFiwIJ5//vmi8fe9733xzW9+MyZPnlx2zq5du+K9731vbNy4sWi8vb09brvttiE/9tGwc+fOaGhoKBrr6uqK6dOnj9EVAQAAADBWhpsVHdElvddee20h7IuI+LM/+7N461vfetjzXnPNNWVh30UXXRS33HJLZtgXETFt2rS4/fbbY+HChUXj69ati+9973uHfU0AAAAAMB4cscBv8+bN8fnPf75w+81vfnN85jOfOex5d+3aFV/96leLxk488cT42te+FhMmTBj03Lq6uvjmN78ZEycWr2T+y7/8y8O+LgAAAAAYD45I4Ld///74yEc+Evv27YuIiAkTJsRNN90UkyZNOuy5v/3tb8crr7xSNHbppZcOeR++008/PS688MKisbvuuiu2bNly2NcGAAAAAGPtiAR+X/jCF+KBBx4o3L7iiiviN37jN0Zl7rVr15aNLVmyZFhzZB2fNS8AAAAAVJtRD/z+/d//Pa655prC7Te+8Y3xhS98YVTm3r17d/zoRz8qGjv11FNjwYIFw5rnHe94R1m1oX38AAAAAMiDUQ38Xn311bj00kujt7e3MPaVr3wlXvOa14zK/Pfdd1/R3BER55133rDnmTRpUjQ3NxeNdXZ2xv79+w/r+gAAAABgrI1q4PfXf/3X8ZOf/KRw+/3vf3/ZfnmHY/PmzWVjZ5999ojmWrRoUdHtnp6eeOyxx0Y0FwAAAACMF6MW+D322GPx3//7fy/cPumkk+Jv//ZvR2v6iIh4/PHHy8bmzJkzorkaGxvLxp544okRzQUAAAAA48WoBH4HDhyISy+9NHp6egpjq1atioaGhtGYvmD79u1lY1nB3VDMnj27bGzbtm0jmgsAAAAAxotRCfz+9m//NjZt2lS4/a53vSs+/OEPj8bURV544YWysRkzZoxorqzzfvGLX4xoLgAAAAAYLyYe7gRbt26NFStWFG6fcMIJccMNNxzutJm6u7vLxqZMmTKiubLO27Vr1yHP2759e2al4WgaynUAAAAAQJbDCvz6+vri0ksvjd27dxfGrrnmmhHvq3co6SXDicmTJ49orqzAb+/evYc87+tf/3p87nOfG9FjAgAAAMCRdlhLer/85S/HXXfdVbj9a7/2a/HHf/zHh31Rlezfv79srK6ubkRzTZo0qWxs3759I5oLAAAAAMaLEQd+Tz31VPzZn/1Z4XZtbW3cdNNNUVMzao1/y0ycWF6Q2NvbO6K5sqr5amtrRzQXAAAAAIwXI17Se9lll8XLL79cuP2nf/qncfrpp4/KRVWStXy3p6dnRMt69+zZUzaWVfVX6tJLL43zzz9/2I83HLt27Yr3vOc9R/QxAAAAAMinEQV+X/va1+IHP/hB4faCBQuKGnccKfX19WVje/bsyRw/lKzAb9q0aYc8b/bs2TF79uxhP95w7Ny584jODwAAAEB+DXv97c9//vP41Kc+Vbg9YcKEuOmmm4ZUHXe4Tj755LKx559/fkRzZZ33ute9bkRzAQAAAMB4MezA7wtf+EL88pe/LNz+2Mc+Fuedd96oXlQlWZV127ZtG9Fc27dvLxtrbGwc0VwAAAAAMF4MO/B75ZVXim5/5StfiQkTJgz5v7/7u78rm/PUU08tO+7BBx8sO27evHllY0899dRw/woRkR0UnnbaaSOaCwAAAADGiyPXUvcIOPPMM8vG7rvvvhHNde+99xbdnjJlSixYsGBEcwEAAADAeFFVgd/ChQujrq6uaOzuu+8e9jx79+6Nzs7OorFFixbFxIkjbloMAAAAAOPCsBOu5ubmePnll0f8gPfee2/Zctq2trY4/vjji8ayOu8ef/zxcd5558Udd9xRGHvyySfjscceG1Z13p133hl79+4tGrvggguGfD4AAAAAjFfDDvz+6I/+KP7oj/5oxA94ySWXlO3j9+UvfznmzJkzpPMvvvjiosAvIuKGG26IL37xi0O+hhtuuKFsrKOjY8jnAwAAAMB4VVVLeiMi3v/+98cJJ5xQNPaNb3wjnnvuuSGd//DDD8dtt91WNPb2t7895s+fP2rXCAAAAABjpeoCv2nTpsXHPvaxorHu7u5YsmRJ9PX1DXpub29vXHLJJbF///6i8U9/+tOjfp0AAAAAMBaqLvCLiFixYkXMmDGjaOzWW2+ND37wg2V78yV27doVF1xwQdx///1F4+3t7fbvAwAAACA3qrItbX19fXzrW9+Ktra2OHDgQGH829/+dmzcuDGuuOKKWLRoUTQ0NMQzzzwTd999d9x4443R3d1dNM/MmTPjpptuOspXDwAAAABHTlUGfhER73rXu+KGG26Iyy+/vGgp744dO+Kqq6465PnTpk2Lf/3Xf42ZM2ceycsEAAAAgKOqKpf0Ji677LK49dZb4+STTx7WeWeeeWbcc889cdZZZx2hKwMAAACAsVHVgV9ExIUXXhiPP/54fPazn41Zs2YNeuwZZ5wRN954Y9x7770xb968o3SFAAAAAHD0TOg7VGvbKtLX1xcPPfRQPPzww7Fjx47o6emJqVOnxuzZs2PRokXR2Ng41pc4JDt37oyGhoaisa6urpg+ffoYXREAAAAAY2W4WVHV7uGXZcKECdHU1BRNTU1jfSkAAAAAMCaqfkkvAAAAAHCQwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAcmTiaE+4c+fO2Lp1azz99NPx/PPPxyuvvBK9vb3xmte8Jk488cSYM2dOnHXWWXHiiSeO9kNHRERfX1888sgj8dBDD8WOHTuip6cnpk6dGo2NjXH22WfH7Nmzj8jjAgAAAMB4cFiBX29vb9x5552xadOm+PGPfxwPPfRQPP/884c8b8KECXHGGWfEH/zBH8SHP/zheN3rXnc4lxEREd3d3fHFL34x1qxZE88++2zF45qamuITn/hEXHLJJXHccccd9uMCAAAAwHgyoa+vr2+kJz/yyCNx+umnH9YF1NfXx+c///lYsmRJTJgwYURzrFu3Lj7ykY/Ezp07h3zOWWedFd/5znfitNNOG9FjHkk7d+6MhoaGorGurq6YPn36GF0RAAAAAGNluFnRmO/h193dHVdccUV8+MMfjv379w/7/DVr1sRFF100rLAvIuKBBx6Ic845Jx588MFhPyYAAAAAjFejuoffzJkzY/HixfGrv/qrsWDBgpg1a1a89rWvjQkTJkR3d3c89thj8f/9f/9f/Mu//Evs3r276Nybb745jj/++LjhhhuG/Hi33357XH755VFapDhr1qxYsmRJLFq0KGbMmBFPP/103H333bFmzZro7u4uHPfiiy/Gu9/97rj//vtj5syZh/V3BwAAAIDx4LCW9D766KOxdOnSuOiii6K9vX3Iy2NfeOGFWL58edx8881l9/3gBz+Id77znYec46WXXooFCxaU7Rn4vve9L775zW/G5MmTy87ZtWtXvPe9742NGzcWjbe3t8dtt902pGs/GizpBQAAACAx3KzosAK/w7VkyZL42te+VjR20UUXxa233nrIcz/5yU/GX//1X5ed+y//8i+D7gXY29sb5557btx///1F4//2b/8WF1xwwTCu/sgR+AEAAACQqKrA7+WXX45TTjmlaJltXV1ddHd3x5QpUyqet2vXrjjllFPilVdeKYydeOKJ8eijjw5pae7DDz8cCxcuLNoz8O1vf3tZ5d9YEfgBAAAAkKiqph1Tp04tq6rr7e2Nn//854Oe9+1vf7so7IuIuPTSS4e8D9/pp58eF154YdHYXXfdFVu2bBnS+QAAAAAwXo15l945c+aUjb3wwguDnrN27dqysSVLlgzrcbOOz5oXAAAAAKrJmAd+PT09ZWMnnnhixeN3794dP/rRj4rGTj311FiwYMGwHvcd73hHTJo0qWjse9/73rDmAAAAAIDxZswDvwcffLDo9pQpU2LevHkVj7/vvvuit7e3aOy8884b9uNOmjQpmpubi8Y6OzuL9vUDAAAAgGozpoHfQw89FHfffXfR2IUXXhh1dXUVz9m8eXPZ2Nlnnz2ix1+0aFHR7Z6ennjsscdGNBcAAAAAjAdjFvg9++yz8f73vz9effXVwlhdXV2sWLFi0PMef/zxsrGsfQCHorGxsWzsiSeeGNFcAAAAADAeHPXA77nnnou//uu/jtNPPz3+4z/+o+i+66+/Ps4444xBz9++fXvZWFZwNxSzZ88uG9u2bduI5gIAAACA8WDiaE+4c+fOuOKKK4rGXn311fjlL38ZTz75ZDz55JNl55x44onxla98JT7wgQ8ccv6sDr4zZswY0bVmnfeLX/xiRHMBAAAAwHgw6oHfK6+8Ev/wD/8wpGMbGxvjQx/6UPzxH/9xNDQ0DOmc7u7usrEpU6YM5xIHPW/Xrl2DnrN9+/bMKsPRdKhrAAAAAIBKRj3wG/IDT5wYCxcujDe/+c1RX18/5PN6enrKxiZPnjyia8gK/Pbu3TvoOV//+tfjc5/73IgeDwAAAACOtDFr2rF///74p3/6p/jwhz8cp556avzjP/7jkM8rNVhX38FMmjSpbGzfvn0jmgsAAAAAxoNRD/zmzJkTfX19Rf/19PTEc889F3fddVf81V/9Vbz1rW8tOufZZ5+N3/3d341Pf/rTh5x/4sTyosTe3t4RXWtWNV9tbe2I5gIAAACA8eCoLOmdNGlSzJgxI2bMmBFve9vb4sorr4z//b//dyxbtqxoT77Pf/7zMWPGjPjjP/7jinNlLd/t6ekZ0bLePXv2ZF7rYC699NI4//zzh/1Yw7Fr1654z3vec0QfAwAAAIB8GrM9/P7rf/2v8Za3vCXe+c53xn/+538Wxv/sz/4s3vve98bs2bMzz8va72/Pnj3D2gcwfV6padOmDXrO7NmzK17baNm5c+cRnR8AAACA/BqzPfwiIhYtWhTXXHNN0VhPT09cf/31Fc85+eSTy8aef/75ET1+1nmve93rRjQXAAAAAIwHYxr4RUR89KMfLeuWu27duorHZ1XXbdu2bUSPvX379rKxxsbGEc0FAAAAAOPBmAd+xx9/fPzar/1a0dijjz6a2VAjImLevHllY0899dSIHjsrKDzttNNGNBcAAAAAjAdjHvhFRMyaNats7MUXX8w89swzzywbu++++0b0uPfee2/R7SlTpsSCBQtGNBcAAAAAjAfjIvDbvXt32djUqVMzj124cGHU1dUVjd19993Dfsy9e/dGZ2dn0diiRYti4sQx62MCAAAAAIdtXAR+P/vZz4puH3/88fGa17wm89jjjz8+zjvvvKKxJ598Mh577LFhPeadd95Ztmz4ggsuGNYcAAAAADDejHng9+STT8YjjzxSNParv/qrg55z8cUXl43dcMMNw3rcrOM7OjqGNQcAAAAAjDdjHvh99rOfLRt773vfO+g573//++OEE04oGvvGN74Rzz333JAe8+GHH47bbrutaOztb397zJ8/f0jnAwAAAMB4NeLA76mnnooHHnjgsB78uuuui7//+78vGpsyZUp86EMfGvS8adOmxcc+9rGise7u7liyZEn09fUNem5vb29ccsklsX///qLxT3/608O4cgAAAAAYnw4r8Fu4cGFcdNFF8c///M/R09Mz5HOffPLJ+L3f+7341Kc+VXbfZz/72TjllFMOOceKFStixowZRWO33nprfPCDHyzbmy+xa9euuOCCC+L+++8vGm9vb7d/HwAAAAC5MKHvUCVxFfzwhz+Md7zjHYXbU6dOjbe97W1x5plnxumnnx7Tp0+PE088MWpra+M///M/4/nnn4+HHnoo7rzzzvjRj36UWYn3X/7Lf4m1a9dGbW3tkK7h9ttvj7a2tjhw4EDR+KxZs+KKK66IRYsWRUNDQzzzzDNx9913x4033hjd3d1Fx86cOTMeeOCBmDlz5vD/TzhCdu7cGQ0NDUVjXV1dMX369DG6IgAAAADGynCzolEL/A7Xhz70ofjGN74x5LAvcdNNN8Xll19+yKW8WaZNmxY/+MEP4qyzzhr2uUeSwA8AAACAxHCzohEv6Z04ceJITy0yZ86cuPXWW+Pmm28edtgXEXHZZZfFrbfeGieffPKwzjvzzDPjnnvuGXdhHwAAAAAcjhEHfuedd15s27YtvvrVr8bv/d7vRWNj45DPnT59erzvfe+LdevWxRNPPBEXXXTRSC8jIiIuvPDCePzxx+Ozn/1szJo1a9BjzzjjjLjxxhvj3nvvjXnz5h3W4wIAAADAeDPiJb1ZXnjhhXj88cdj27ZtsXPnznjllVfi1Vdfjde85jXx2te+NhoaGuKMM86IN77xjaP1kGX6+vrioYceiocffjh27NgRPT09MXXq1Jg9e3YsWrRoWMHkWLGkFwAAAIDEUdvDjyNH4AcAAABA4qjt4QcAAAAAjD8CPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMTR3vC3t7eeOyxx+I//uM/YufOnfHSSy/FpEmTYtq0aTFjxoxobm6OhoaG0X7Ygr6+vnjkkUfioYceih07dkRPT09MnTo1Ghsb4+yzz47Zs2cfsccGAAAAgLF22IHf/v37Y+PGjfH9738/7rzzzrjvvvti//79g57zpje9KX7/938/Lr/88njDG95wuJcQERHd3d3xxS9+MdasWRPPPvtsxeOampriE5/4RFxyySVx3HHHjcpjAwAAAMB4MaGvr69vJCc+88wz8d//+3+Pf/7nf45f/OIXI3rw2tra+MxnPhOf+cxnYuLEkWeP69ati4985COxc+fOIZ9z1llnxXe+85047bTTRvy4R8rOnTvLqiC7urpi+vTpY3RFAAAAAIyV4WZFI97D7z/+4z9izZo1Iw77IiL27dsXf/EXfxFve9vb4pe//OWI5lizZk1cdNFFwwr7IiIeeOCBOOecc+LBBx8c0eMCAAAAwHg06nv4nXbaafG2t70tFi9eHK9//eujoaEh9u/fH08//XTcddddcfPNN8euXbuKzvl//+//xbvf/e74/ve/H5MnTx7yY91+++1x+eWXR2mR4qxZs2LJkiWxaNGimDFjRjz99NNx9913x5o1a6K7u7tw3Isvvhjvfve74/7774+ZM2ce1t8bAAAAAMaDES/p/f73vx+tra0REXHyySfHhz70obj00kvjjDPOGPS8l19+OT71qU/FV7/61bL7Pve5z8VVV101pMd/6aWXYsGCBfH8888Xjb/vfe+Lb37zm5nB4a5du+K9731vbNy4sWi8vb09brvttiE97tFgSS8AAAAAiaO2pDcioqGhIa677rrYvn17/M3f/M0hw76IiKlTp8ZXvvKVzGDvC1/4wpCX5l5zzTVlYd9FF10Ut9xyS8UqwWnTpsXtt98eCxcuLBpft25dfO973xvS4wIAAADAeDbiCr8XX3wxJk2aFCeccMKIHrivry/OPffcuOeee4rGb7zxxrjssssGPXfXrl1xyimnxCuvvFIYO/HEE+PRRx8d0tLchx9+OBYuXFjUTfjtb397WeXfWFHhBwAAAEDiqFX4nXTSSSMO+yIiJkyYEP/tv/23svENGzYc8txvf/vbRWFfRMSll1465H34Tj/99LjwwguLxu66667YsmXLkM4HAAAAgPHqsJb0Hq63ve1tZWM///nPD3ne2rVry8aWLFkyrMfOOj5rXgAAAACoJmMa+J188sllY7/4xS8GPWf37t3xox/9qGjs1FNPjQULFgzrsd/xjnfEpEmTisbs4wcAAABAtRvTwO+5554rGzvppJMGPee+++6L3t7eorHzzjtv2I89adKkaG5uLhrr7Ows2tcPAAAAAKrNmAZ+mzdvLht74xvfOOxzzj777BE9/qJFi4pu9/T0xGOPPTaiuQAAAABgPBjTwO/v/u7vysbOP//8Qc95/PHHy8bmzJkzosdvbGwsG3viiSdGNBcAAAAAjAdjFvht2bIl/u///b9FY8cdd1y85z3vGfS87du3l41lBXdDMXv27LKxbdu2jWguAAAAABgPxiTw6+vriyuuuCL27dtXNP6hD30oZs2aNei5L7zwQtnYjBkzRnQdWecdqmkIAAAAAIxnE8fiQa+77rq44447isamTp0af/mXf3nIc7u7u8vGpkyZMqLryDpv165dg56zffv2zCrD0XSoawAAAACASo564Pf9738//vzP/7xsfNWqVfGGN7zhkOf39PSUjU2ePHlE15IV+O3du3fQc77+9a/H5z73uRE9HgAAAAAcaUd1Se+///u/R0dHR7z66qtF4x/4wAfisssuG9Ic+/fvLxurq6sb0fVMmjSpbKx0mTEAAAAAVJOjFvg99dRTccEFF5QtyT3nnHPipptuGvI8EyeWFyX29vaO6Jqyqvlqa2tHNBcAAAAAjAdHZUnvs88+G7/1W78VzzzzTNF4U1NTrF+/Po4//vghz5W1fLenp2dEy3r37NlTNpZV9Zd26aWXxvnnnz/sxxqOXbt2HbJbMQAAAABkOeKB3/PPPx/vfOc742c/+1nR+K/8yq/Ehg0bYtq0acOar76+vmxsz549meOHkhX4Hep6Zs+eHbNnzx72Yw3Hzp07j+j8AAAAAOTXEV3S+8ILL8T5558fjz32WNH4ggUL4o477ojp06cPe86TTz65bOz5558f0fVlnfe6171uRHMBAAAAwHhwxAK/F198MVpbW+ORRx4pGj/ttNPijjvuiBkzZoxo3qzqum3bto1oru3bt5eNNTY2jmguAAAAABgPjkjg193dHe9617viwQcfLBo/9dRT44477ojXv/71I5573rx5ZWNPPfXUiObKCgpPO+20Ec0FAAAAAOPBqAd+v/zlL+OCCy6I++67r2i8sbEx7rzzzjjllFMOa/4zzzyzbKz0sYbq3nvvLbo9ZcqUWLBgwYjmAgAAAIDxYFQDv5dffjna2triJz/5SdH4KaecEnfeeeeoLJdduHBh1NXVFY3dfffdw55n79690dnZWTS2aNGimDjxqDQuBgAAAIAjYtQCv927d0d7e3ts2rSpaPwNb3hD3HnnnXHqqaeOyuMcf/zxcd555xWNPfnkk2WNQQ7lzjvvjL179xaNXXDBBYd9fQAAAAAwlkYl8NuzZ09cdNFFcddddxWNz5o1K+64446YO3fuaDxMwcUXX1w2dsMNNwxrjqzjOzo6RnxNAAAAADAeHHbgt3fv3njve98bd9xxR9H4zJkz484774z58+cf7kOUef/73x8nnHBC0dg3vvGNeO6554Z0/sMPPxy33XZb0djb3/72I3KtAAAAAHA0HVbg19vbG7/7u78b3/ve94rGGxoa4o477jhiDTCmTZsWH/vYx4rGuru7Y8mSJdHX1zfoub29vXHJJZfE/v37i8Y//elPj/p1AgAAAMDRNuLAb//+/fH+978/1q1bVzQ+ffr0uOOOO+Itb3nLYV/cYFasWBEzZswoGrv11lvjgx/8YNnefIldu3bFBRdcEPfff3/ReHt7u/37AAAAAMiFCX2HKomr4Itf/GL8yZ/8Sdn4mWeeedh79n33u98d0nG33357tLW1xYEDB4rGZ82aFVdccUUsWrQoGhoa4plnnom77747brzxxuju7i46dubMmfHAAw/EzJkzD+uaR9POnTujoaGhaKyrqyumT58+RlcEAAAAwFgZblY0caQP9Mtf/jJz/MEHH4wHH3xwpNMOy7ve9a644YYb4vLLLy9ayrtjx4646qqrDnn+tGnT4l//9V/HVdgHAAAAAIdjVLr0jqXLLrssbr311jj55JOHdd6ZZ54Z99xzT5x11llH6MoAAAAA4Oir+sAvIuLCCy+Mxx9/PD772c/GrFmzBj32jDPOiBtvvDHuvffemDdv3lG6QgAAAAA4Oka8h9941dfXFw899FA8/PDDsWPHjujp6YmpU6fG7NmzY9GiRdHY2DjWl3hI9vADAAAAIHHU9vAbryZMmBBNTU3R1NQ01pcCAAAAAEddLpb0AgAAAAD9BH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMTx/oCRlNfX1888sgj8dBDD8WOHTuip6cnpk6dGo2NjXH22WfH7Nmzx/oSAQAAAOCIGtXA78knn4x777037rvvvsJ/u3btKjrmD/7gD+Kb3/zmaD5sdHd3xxe/+MVYs2ZNPPvssxWPa2pqik984hNxySWXxHHHHTeq1wAAAAAA48FhBX4//OEP49/+7d8qhntHw7p16+IjH/lI7Ny585DHbt68OS677LL40pe+FN/5znfitNNOOwpXCAAAAABHz2EFfn/zN38T//Iv/zJa1zJsa9asiY9+9KPR19c3rPMeeOCBOOecc+IHP/hBnHnmmUfm4gAAAABgDFTtHn633357XH755WVh36xZs2LJkiWxaNGimDFjRjz99NNx9913x5o1a6K7u7tw3Isvvhjvfve74/7774+ZM2ce5asHAAAAgCNj1AO/mpqaePOb31xokvE//sf/GO2HiJdeeik+/OEPx4EDB4rG3/e+98U3v/nNmDx5cmFs0aJF8du//dvxmc98Jt773vfGxo0bC/ft2LEjLrvssrjttttG/RoBAAAAYCwcVuA3ceLEeOtb3xpnn3124b8zzzwzTjjhhIiIeOqpp45I4HfNNdfE888/XzR20UUXxS233BITJkzIPGfatGlx++23x7nnnhv3339/YXzdunXxve99Ly644IJRv04AAAAAONoOK/D7zne+EzU1NaN1LUOya9eu+OpXv1o0duKJJ8bXvva1imFfoq6uLr75zW/GwoULY//+/YXxv/zLvxT4AQAAAJALh5XWHe2wLyLi29/+drzyyitFY5deeumQ9+E7/fTT48ILLywau+uuu2LLli2jdo0AAAAAMFaOfmJ3mNauXVs2tmTJkmHNkXV81rwAAAAAUG2qKvDbvXt3/OhHPyoaO/XUU2PBggXDmucd73hHTJo0qWjse9/73mFfHwAAAACMtaoK/O67777o7e0tGjvvvPOGPc+kSZOiubm5aKyzs7NoXz8AAAAAqEZVFfht3ry5bOzss88e0VyLFi0qut3T0xOPPfbYiOYCAAAAgPGiqgK/xx9/vGxszpw5I5qrsbGxbOyJJ54Y0VwAAAAAMF5UVeC3ffv2srGs4G4oZs+eXTa2bdu2Ec0FAAAAAONFVQV+L7zwQtnYjBkzRjRX1nm/+MUvRjQXAAAAAIwXE8f6Aoaju7u7bGzKlCkjmivrvF27dh3yvO3bt2dWGo6moVwHAAAAAGSpqsCvp6enbGzy5Mkjmisr8Nu7d+8hz/v6178en/vc50b0mAAAAABwpFXVkt79+/eXjdXV1Y1orkmTJpWN7du3b0RzAQAAAMB4UVWB38SJ5QWJvb29I5orq5qvtrZ2RHMBAAAAwHhRVUt6s5bv9vT0jGhZ7549e8rGsqr+Sl166aVx/vnnD/vxhmPXrl3xnve854g+BgAAAAD5VFWBX319fdnYnj17MscPJSvwmzZt2iHPmz17dsyePXvYjzccO3fuPKLzAwAAAJBfVbWk9+STTy4be/7550c0V9Z5r3vd60Y0FwAAAACMF1UV+GVV1m3btm1Ec23fvr1srLGxcURzAQAAAMB4UVWB37x588rGnnrqqRHNlRUUnnbaaSOaCwAAAADGi6oK/M4888yysfvuu29Ec917771Ft6dMmRILFiwY0VwAAAAAMF5UVeC3cOHCqKurKxq7++67hz3P3r17o7Ozs2hs0aJFMXFiVfUwAQAAAIAyVRX4HX/88XHeeecVjT355JPx2GOPDWueO++8M/bu3Vs0dsEFFxz29QEAAADAWKuqwC8i4uKLLy4bu+GGG4Y1R9bxHR0dI74mAAAAABgvqi7we//73x8nnHBC0dg3vvGNeO6554Z0/sMPPxy33XZb0djb3/72mD9//qhdIwAAAACMlaoL/KZNmxYf+9jHisa6u7tjyZIl0dfXN+i5vb29cckll8T+/fuLxj/96U+P+nUCAAAAwFiousAvImLFihUxY8aMorFbb701PvjBD5btzZfYtWtXXHDBBXH//fcXjbe3t9u/DwAAAIDcmNB3qLK4Q8jaUy+xe/fuWL9+fdFYY2NjLFq0qOI5S5cujd/8zd885OPefvvt0dbWFgcOHCganzVrVlxxxRWxaNGiaGhoiGeeeSbuvvvuuPHGG6O7u7vo2JkzZ8YDDzwQM2fOPOTjHU07d+6MhoaGorGurq6YPn36GF0RAAAAAGNluFnRYQd+EyZMOJzTy3zjG9+ISy65ZEjH3nTTTXH55ZcfcilvlmnTpsUPfvCDOOuss4Z97pEm8AMAAAAgMdysqCqX9CYuu+yyuPXWW+Pkk08e1nlnnnlm3HPPPeMy7AMAAACAw1HVgV9ExIUXXhiPP/54fPazn41Zs2YNeuwZZ5wRN954Y9x7770xb968o3SFAAAAAHD0HPaS3vGkr68vHnrooXj44Ydjx44d0dPTE1OnTo3Zs2fHokWLorGxcawvcUgs6QUAAAAgMdysaOLRuKijZcKECdHU1BRNTU1jfSkAAAAAMCaqfkkvAAAAAHCQwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwAAAADIEYEfAAAAAOSIwA8AAAAAckTgBwAAAAA5IvADAAAAgBwR+AEAAABAjgj8AAAAACBHBH4AAAAAkCMCPwA4xrXWdIz1JQAAAKNI4AcAx5iW9pUV72ut6YiW9pWDHgMAAIxvE8f6AgCAo6dp2aqIubWF2601HbHhwNrCnyMiNq67ckyuDQAAGB0q/ADgGNFa0xF7GiJeevOrmfcnwZ8KPwAAqG4q/ADgGDKlK2JK13GFQG/jQMiX6Fq6ODZfv3wsLg0AABglAj8AOEb0tjVHw+pN8cS3Fsbm6/+86L6W9pWx49za6K1Q/QcAAFQPgR8AHCO659ZGQ0S89ieTo+WWg0t269Z3xsYDa6O1piMeX33O2F0gAAAwKgR+AHCMqN+6L3rbmuOlN78a+965OyIifvrbfxERBxt2nPjocWN1eQAAwCjRtAMAjhF16zujbn1nHP/6l6PxpF3ReNKuQtC34cDa6G1rjvqt+8b4KgEAgMMl8AOAY0TShXf3s1Pj0e0z49HtMwv3tbSvjLr1ndE9t3asLg8AABglAj8AOIZsOLA23vKrP4++X9ZG3y9rY8PA3n116ztjw4G18ZKmHQAAUPXs4QcAx5DWmo545h/fWjSWVP61tK+MunNV+AEAQLVT4QcAx5ANB9ZG7R0nxvGvfzmOf/3LEdEf9EVEbFx3ZTS3/nQsLw8AABgFAj8AOMbUb90Xe16eFHtenhQR/UFfRH/13wM73jCWlwYAAIwCgR8AHGPq1nfG4nk/i8XzfhYRUejU29vWXAgBAQCA6mUPPwA4BqUr+XrbmiOiPwh8ct3asbokAABglAj8AOAY07V0cdTeMXDjtw8u6e1aunjMrgkAABg9Aj8AOMY0rN4UE3/4+rLx+q37orWmo9C1FwAAqE728AOAY0xvW3Ps/qs3xO6/6l/Wm+zhV7e+M5741sKxvDQAAGAUqPADgGPMxnVXxqk3f75wO6no621rjr5f+i0QAACqnU/1AHCMaWlfmTm+49zamLf0nqN8NQAAwGgT+AHAMWbjuivjtT+ZHK/9yeTCWGtNR0zpinh89TljeGUAAMBoEPgBwDGmtaYj9r3zpdj3zpcK+/dtOLA2fvlrPSr8AAAgBwR+AJBTrTUdhUAvbcOBtbH72amx+9mpRR15X/uTybHt6sVH8xIBAIAjQOAHAFWqadmqEZ3XWtMRs057IWad9kLRWP3WfbFlxfLRujwAAGCMCPwAoErVb91X8b50Y47SKr8NB9bGL34yI37xkxlFY3XrO0f/IgEAgKNu4lhfAAAwMoMFdBvXXVk21lrTERsOrI3Wmo7YMrCUd/61q2LLiuWF+wAAgOqnwg8AqtS2qxdX3KcvLd2YozTYe99/uSualq2KrqX9e/elKwMBAIDqpMIPAKpU41WbysZa2lcWKv+SEK9hkDlu3vxrMXtgaXBL+8rMykAAAKC6CPwAoEr1tjUX/tzSvjK659bG5nVXli/Pvb5yI47P/tq6+Ltb3tN/7vXLhX4AAJADlvQCQJWqW98Z3XNro3tubWxcd2U0rC6v+Is4uEw3a+nv1d//7aJzd5xbGy3tKy3tBQCAKqbCDwCqUGtNR/S2NReCupat+6J76eLMUK9ufWfh+KR6Lwn0GiMiPt5fLdha0xGNcXApMAAAUJ1U+AFAFSrtqLtx3ZWx+frlg3baTS/VrVvfGXXrO+ONVz0eERHdc2tjw4G10dvWHPVb90X9wL5+AABA9VHhBwBVauPAfn0RBxtupCv8kiq+uije7y9Stx/YsTtaazqifuB23frOQUNDAABg/BP4AUAVS8K51pqO/tDvwNpC0Nc9tzYiIjanxhLJfdee/i/xOwceiIiIOV++LuYdrQsHAACOGIEfAFSp0m683XNro7WmIzYOjG0e6M6bHksk97W074vfWRfxoXs+GhNe+6boSu0DqNIPAACqk8APAKpYUrmXBHqtGZ16NwxU+KX38Etsa+/fzvdnL70u3jz7udh9yxsEfQAAUOU07QCAoyiri+5IpYO5rL370rLCvoiI41//crS0r4yXdk+OkybtLiz1BQAAqpcKPwA4CtINNIaidLlu6e1E3frOiIiipbjJWHovv4aByr/SOXY/OzXq1t8T06M5ftS+IOLNrxZdc6WgEAAAGL8EfgBwFAwnOMsK99LNOTYcWFsI99Ldd5Olu3XrO4v27asUFkZE1L14XOHPx7/+5ai948SY8+XrIiLiKWEfAABUJYEfABxF6cAuHeJF9Id3G9ddOaQ99EqPmX/tqmit6cisIBxsvt6TXi1cT+3cxbH5+uWF62lZp8IPAACqkT38AGAMpEO4DQfWxoYDawtLcYdyXhLKtdZ0RGtNR2xZsbxwX1L1l9w3mAmv3RdNy1ZFRET91n1Fxwv7AAA41iXb5Az2uTrrvkOdN//aVRXnSn/WT+Z6zwf+19AvOgR+AHDEpZtopJfgJuMt7SujpX1l0VLdQ81Vemwy3rV0cXTPrY0NB9ZGb1tzUfiXZfYtx8Xm6/vDwu65tdHb1lwIIAEA4FiUhG4t7SsLTe2G8vk4OSfZZqfSeckP9pU+o6fPmX/tqqhb3xl1P3hwWH+HCX19fX3DOoMjbufOndHQ0FA01tXVFdOnTx+jKwJgtGS9qVf6EJDsyRdx6Gq79PE7zq2Nxqs2DTp/Wrq5R/3WfYN+OAEAgLxLPrN3LV0cERGbr18+5HOT/bAjIhrXHYiIg0300qt1Hl99TjSuO1DUYC8iirb/SX936O3bG3fF/y16rMGyInv4AcARln6jLt2/r/S4ZDw5brCuvkUfANqaC4060lWApfOW3k5+fextay58CElXJAIAwLGm9LN6a01H9LY1R936zkF/FG9pXxl15/ZXBG5ZsTzi4/1j6bAvoj9IPPHRiO65/attWkvuL/1zS/vK6N37csQPigO/wQj8AOAoSAdsGyt8SKgUypVKgr7HV58T85beExFRtP9faTVg19LFZWFiIjmve25tbFbRBwAAhc/byY/w265eHI1XbRr083myJ3fj+v6xlh/vi4j+z+ale/JtLvncX+mxn/jWwmit6YjupYtj/57aYf0dBH4AcISVVs1lVdwlx5V27y39c+lYskdf99zaog676SUADRXmTh/XsHpTxMBSBc06AAA4liWfmZuWrYr6tubYsmJ5zD/EOcmx6R/Uk6W6yfeBZPVO6TY8WXt5Ny1bFbO3Hhe9bc2xpyFi/yvD+zvYw28csocfQL6ULulNa2lfWTFgG6zKr3Tu9IeS0j1GBpun9MNGb+pDij38AAA4FpX+ON61dHFReFd6bNaP+Yn0fUkRQPJ5u7etOTauu7LoM3x6b+5k652IiN69L8f/+8E1RXMPlhUJ/MYhgR9Afg0lxBvs3FKV9umrdG6lDyiV7gcAgGNV0qV38/XLB/2hPjm2bn1nPL76nIjob9iRDuzS+/81LVtVaJZXuhIn/TjJcRHDD/ws6QWAMVIppKvUmTer4UcS1CVNO5JfCQ/1OOnwsGnZqsI+IoNVIwIAwLFk47oro6V9Zcz58nXx1MBefJU+V29MBXYR/QFf19LFsadhoIFHyktvfjX2NNRGb/s5ZXt8pz/LJ6Fglz38AKC6ZH1oGOyXw3RFX9fSxRFLF0dEFJYADKXDbulefg2rN0Xr6k3RtXRx1A/8AgkAAMe6pGpv3vqIlnXZzfdKP8unl+sW9thecXC5bkTEiUsXx+brl0fTslVRt/6esgq/ZI4dVy+O3vZzou7F4e/hJ/ADgKOodAPg4TTISH8QSD5ANC1bVXRf8ivkUELDZJ7k2MNZbgwAAHm17er+H9lLV8Ok99uL6G/UUZ9qqpeWfP6P6K/ca63piD0D85Y2+Yvo/5zee9KrceKjx8We4l3fhsQefuOQPfwA8ie9R978a1fFrB/vKwraEsmGwI+vPiee+vgny85PK13au+3qxTGlK4bUtONQ8wEAwLEu/Zk5+aE8vaVO6Q/npZ+7W9pXxrb2msLn+nQwmDQBSeZKAsL0foHpeXv79sZd8X+Lrs8efgAwxtJv/I1XbSps3hsRhT8nb+q9bc1FYV9yTPpNv6V9ZVlo13jVpkOGdsmHh9KwMOsDCgAAHMu6li6O+q37Cp10s/bT3nBg7cFQrq25sAKnYfWm2JhU7n28/3P49g+8GhERT/7+lREDy317U1V/ERHzr10VcW5tWTVhVnHYYFT4jUMq/ADyK/2rXkT/B4F0+Nc9tzbqt+4r+9Uv/cGiadmq2NMQRVWCTctWFfYISe/Pl3V++r6sJh8AAEC/+deuisarsj+XJ5If1dNdddN7+SXb7qSVjqWPjyiuAqwU+KnwA4Bxom59Z2H/vdJuvK01HRFLF0fd+s7+rrklS3OTY7J+12tYvSmaIgpNPCKyQ77SSr7Wmo5D7vkHAADHkvSP443rOwsrbNJLetOfq+vWd0ZrTUfUtzUXgrvk/uT2xnVX9lfvRf8P93O+fF00RvEqn4iDQd/mA2ujdfWmwpZA+195eVh/h5rR+D8CABi6uvWdmSHbhgNr45e/1hMR2XvstbSvLGra0dvWHK01HYUPCJuvX174RXEk1wNHStbzmeoylH/D0mOG8++evJYd7nUk93vOAXA4ks/ayVY4869dVQj1kiW8EVH2Y3pyfHL7mX98a0T0vy81LVsVW1Ysjy0rlkf33NpoXHegUOVX+kN9et7Wmo5ovGpTnPKFzmH9HSzpHYcs6QXIj/QbdUQU/TmRvMGXNvMY6nzp8aZlq6Jh4JfAdAVhpeW99u8DBlPptab0zxExomrhrHOG8ro0WDOi0te4NK9zAAzXnC9fF/OW3lP2OT69v1+6IjBRN1AZmJyfHB8RsePc2ug96dVoXHegcFzy2X3HubUx68cHf8RP3ict6Y2ILVu2xIMPPhhPP/107N69O44//vh44xvfGGeddVbMmzdvrC8PgGNI6ZfLlvaVEQOl/r2pkv+IiC0r+pf5trSvLLzxJ7ICw6wvu8lHgEoVe77scrQMFgwJmKtHuklQ1v6gye0NB9ZmLnNK/lx6fPL6t7FkztKqicGuK0vXwLYGyTWVfvGKGFkwCcCxqbWmI2L1OYU/b7t6cSGM2zyw/U5rTUc88a2FcdqH74+I/sAuIqJxff8cJz56XOH9KTGlK2LWjw9ERBSFfRH9y33T3xGS983evcNb0pubwK+npye+9KUvxVe/+tV44oknKh532mmnxRVXXBEf//jHY/LkyUfxCgE4lpV+UU6/safHkqBvsKqUdKONwapvkg5h869dFY1RueqlpX1l4Us3HI5Kyyi3XV38IVfYV13Sm4pnBXnpf89KIW/aUP/9s4Li5PUvXc2cft4llRPbrl7cP57aD6k+1REdAIai8Ln60YMN8rKOaWlfGb1tzbGtvSbmLT3YaCP5Qb63rTl2nFsbDav7g7z6VCOPpmWromF9Z3EomFE9OFy5CPzuueee+MAHPhBPPvnkIY994okn4k/+5E/iS1/6Unz729+O5ubmQ54DAKMlXS0Tqb1BIoq/3B7qC3Fhue7AcUnFSnpD4dIuv+kvy6XL8EofHyqpVME1WFfo9HmV7md8Ke3kPdSQLHk9qYso2ty80lYCg0lX60UUh4lZjY26li4uVFvM+vG+soqJ9DUK/QAYqpb2lbF54H2j0mec7oH3oJb2lfHEtxYWxpOKwKwuvcmqnvq25njmH98ab/id/s/uT3xrYbz2J5ML+2zXRf8PV/v31A7ruqu+acf69eujpaVlSGFf2s9+9rN4+9vfHhs2bDhCVwYA/ZJuXBFR6NS14cDa6J5bm/ml81Cb16c7fSXSHyKymoEk86WXBySPlT4ODqV0yWWyd+Rgz58tK8o7TjO+Jf+2yQblEeXLbktfq1prOgrhYLLEt9LcWVWCpc+hSpUUyWMVfjwZsKchil4HkyrnRDJW+joIAJUk7zNJJV7y59L3oOQ9a+O6K6P2qclR+9Tk2HBgbUzp6r+/admq6J5bG11LFxd+0Nq47srCD2vTbzy+MNfsW46LhtWbYtvViw/rPauqm3Y8/PDD8eu//uuxe/fuovH6+vr46Ec/Gr/xG78Rb3jDG+K5556Le++9N2644YZ47rnnio494YQT4ic/+Un8yq/8ytG89EFp2gGQH1lVc8mveenqvmQ8Isr2wSqdJ2tT4MEacxQqbgY+MJRWFfZa5sYhZD3/BqvYynquJs+vrD0qGb+SpbMRB/8dS1+Dsl6z0ga7r5JKr03JXF1LFxeW9KavdU9Df3Vfshn6vKX3FK41/aXJ8w+AwZT+4FT6ub30uPT9yY/9jVdtyvxMnijdoieRvMclf04qB3v3vhz/7wfXFB87SFZUtYHf/v3746yzzopHHnmkaPw3f/M345/+6Z+ivr6+7Jyenp74gz/4g/jOd75TNN7U1BT33XdfHHfccUfykodM4AeQL8PpNFmpSi/r+KyAr1LDj8G6WSYfSlRhMVLpJgwR5R+Kk+dlhKCl2iSvH1n/hqXddFvaV0b33P7lRun99ZLnRtbrUHoZbuljlv650jFZP2pEZAeFla4FANJKf4iv1D0+/R6Ytcd28l6X/vEsmTf9WOnOvMl96f1qRxL4Ve0efl/96lfLwr6FCxfG9773vairq8s8Z/LkyfHtb387du/eHbfddlthfPPmzXHTTTfFkiVLjug1A3Bsygrfmpb1b/yb3mMqeXPP+pKaXvaW7j45lDL/0q5gpR9WGq/aVFQtCFkG67Q7WHhSGkKr8KtOgwW2yevSxtTrWelWApX+vUvDvoiDX4ayvlilXzuzJK+LXUsXF1UJRqrZR73XOwCGKL2kN9mfNusHsCSFamlfWWgYlQ77Sn+wT4d82z/wapz24eLPSt1za+OlN78aSSnYSPbwq8oKv1dffTXe9KY3xfbt2wtjEydOjPvvvz9OP/30Q56/Y8eOeMtb3hIvvfRSYezUU0+NJ554Impqxn5bQxV+APmV9Stg1v1DWQJXWuGSLBMu7fSb/gJcaU+sSssUIGJoDV2SICZ9fGnlX4QKv2pTqRthRMSpN38+Zt9yXMXGHGmlrzHpJkKVljJlPVfS1XxZS4x3nFtbqFYurbboWro49jSoZgZgZE69+fNx2ofvr7hyJqJ4OW5phV9vW3OhEj7pwNu0bFXUb93X301+oNN8RMS29pp46uOfLJxfv3XfsVHhd/vttxeFfRERF1100ZDCvoiIWbNmxUc+8pH4m7/5m8LYk08+Gd///vfjXe9612heKgCUSSphEskGwPVb90VdFP8amHVu+sNDWmmwUmlvvqwv5UnTD5VXZBnKc6Jh9aZC19R0JWo62IkY2X5ujJ109W8SsiWe/P0/j5Zbiiv6kvu7B16f0puYlyq87kT582H+wG/fpUFi8uWoIcqbEEX0VyzHiuVFy4sjivdA0qUXgJF47U8mR29bc9FnmeR9LHn/25Oq3UqOSz7Tb2uviboX+++r3zrwI1aqiV/pvPHxKJw/ktU4VRn4rV1b/gFxuMtxL7/88qLAL5lX4AfAkVBabZf+JbC+wl4eh5Ku1kuWESQfELK+0GZVaaUrAX0BppLSvWiSPw8mfX9pEO25Vj12nFvbH6JFeTOWrE68SchWeH26fnnhy1FWA6Os50JL+8rYUqFqsLetOWKgy2H91n3RWtMR264++ONH4/r+/93WXhMnPpqquhhY4lsX5dscAMBQpCv3Iop/xJ/z5evixEePK9oqJ/2DetfSxdG4bl9EHCjMl7xnFrb3GZiradmqqI/iH/o3rrtyYDVocYXfYKpySe8b3/jGeOaZZwq3J0+eHC+99FLFvfsqmTNnTmzbtq1w+5RTTimrHBwLlvQC5F/p8sdEeolupcCuUviS3lMkLakWzNocPzH/2lVlncSg1KGaKSTVqqXP7UM1kWF8S3fqTe87OlhFcGtNRzy++pzCcqSyCr5rV0XvSa9G47oDma91yWOVXsfm65fHnC9fV+i+W9qlNwkakx9Xki9dyWvjtqsXx5aBCkDBMwDDMf/aVTHrx/vKfqRPV+ilf9hM/4CfhH/p6vP6rfuKtqJI3v/S71Xp+4ebFVVdhd/Pf/7zorAvIqK5uXnYYV9ExNve9raiwC+Z+w1veMNhXycAJLK+7KaXPxaWwA18AGit6YgdV5dXoKSXDZQq/bCR/lK+p+HQ1VhbViyPsK8Vg2hatio2Z4TM6f3b0nvPlAbUCQ1iqkO6e3dpZ8HSHxaSf+PSrQga1x0oLEcqPX5Wsv/e0oObl0f0fyGqFCo3rN4Uras3xVMH1kbr0o6y51LD6k3RsnVfbExte1AaPidLftNfuABgKKZ09Ve+p6XDvdLP6ElTjuScdFCY/Ei6ZcXyg/vMDnz+b7xqUzR1RWweOH7+tatiSlfE/j0vD+t6qy7w27x5c9nY2WefPaK5Fi1aFDfffHPR2IMPPijwA2BUVarKSyRfnhsiomVgH79ZP95XdEzyxTfdkKOlfWWh9L+1pqNo2XBrTUehQ1h6rsE6rWbdhoj+58XmIXTkTT7IJnuspc8vHDNIhSDjx2CNLdL/bqWvQ4l0hcNgVYDJfVn7F5UGeslrX/L8yqrQ23FubVHY3DQwnoTVyXhDROFHFwAYivqt+6JhdWfZj+TpfbCTH8wiDoZ5U+b2B3lJcJfMldh29eKipcAR/T/YJz+sJvf17i3+fnAoVRf4Pf7442Vjc+bMGdFcjY2NZWNPPPHEiOYCgEoqdcVNvnxWWraWlvUFOx3ylR6TLNEdimS5Xta1QER2o5dk/7bSZb4R/eF1aUVq6XL1dGMPxrfk3yj5EpMV4JaqtNwpma/Sv31pg47Sx4s4uIdSMkdSJdEQ/VURT3xrYWFZ1S9/rSciDlaozvnydTFvaH9tACgyWAOy5P2ocfWmokDwmX98a+x/ZKCiL/Vr6Lb2mpi39J5oaV8ZUwaqzgsdfNuaC2Fh8v4VEcPew6/qAr+sPfaygruhmD17dtlYeokvAIyGdFVfS/vK2HH14pj1430VlzYmy32zqu/SS9S6li4udAJLPhQkIUtjam+tSnsERvQ3+0j/wgiHktX4JSLVUCEOfiBOhzylDWay5mJ8Sl5b0q8zpfdHFL8mpfcxyqrEy2oAkn6sLJU6j6eXEXctXRy1Tx183T1t/cFzW2s6onFgDo07ABiJrB/y69Z3Rkv7ymgYaM6RNv3G4yNiX3TPrY3ek16Nib/avyx33u/8eyEkrN/a/70gvYVG0VYqI6xIr7rA74UXXigbmzFjxojmyjrvF7/4xYjmAoChqFvfGbOiObrn1hZVqUQc/MJb+r/pY5I9qurWd0Z9W3Nsvv7gviFJ1d+OqxfHrGiuuIw43VmztaajaD8uG9mTJb0UslIYc6il4cntpBqL6lEaxqUDvs0DXXjr1ndG/dZ9hdek1qs2VawATc9b6XmSKH0upTdAz/rSVXhdvWpTPL76nP59BCMK15I8xlA6oQNAqSe+tTCe/P0/LxtP3mdKf0jfuO7KQqi3+for49SbPx8R/e+hrTUdUd/WXNjrLwkLG1Zvij1XL664d/dQVV3g193dXTY2ZcqUEc2Vdd6uXbsGPWf79u1HvJPvoa4BgOqTDjs2X788fuWf/yI2XH/ofa1Kzy8N6LKqpkqrcNJ7+5V+AY8o3mg/q8svZIXQEeWBy7aBzaZnpZaslwbaSRfVrPkYnwpdb6P/Nax+oEouqTZOXleSCoeI/tek9BeVrHAuvedRJUN5jlTaFqG0c3nyWtda03FYX6AAOHad9uH7o+WWg42mIg4u560f+AEsLV35N//aVfHarskREYWqvu65tYVmIMl2PKU/so50C5SqC/x6enrKxiZPnjyiubICv7179w56zte//vX43Oc+N6LHA+DYlXzZ3TzwgeANq/89WqPy0rhKYUn6z1mhXlYn36zwsGzegS/uAhgGk1UBmu7Wmz5uKCGOPfyqx8Z1V0ZL+8qo37qvqEK40Lk3yqs802FeVgVo+nmSVf1XSTpUrBQqJ2FfUi2RVFEkP3yoNAVgJNJNqdLvZ5uT95aByr3EttXnRGM0D1T+HezwW7e+Mx5ffU7MW9q/92ztU5OL9zvO+Ow/XFUX+O3fv79srK5uZL/RTZo0qWxs3z77GAEwukpDjYbVmwpL4NIfCLKWTabf+OvWd0ZvW3Phy27pB4HkvKyN7gdbNid0YSiSqqh0uFLpedM9t7aoOUfEwZAmXWXqeVc9KlXjpZfKRmTvyVcquS/ZnLyS0tem+deuilmpx21atqroS1H6Grvn1kZDHNyjNP0amVRb6NILwHClf+DcdvVAuDfw+XzzwI9j6aW585beU7TS56U3vxoR/Y055i29J3rbmuPJ37+yaKuJ5P3xcFdEVF3gN3Fi+SX39vaOaK6sar7a2sE/eADAcB1qj6pES/vKiuFb6ZfqJCzM+lI9lDAlqW6p37qvrMum5ZakVeqqm67qS56fSbVX8oE3rfTDa/Jnz7PxLVmilPxbz/nydTFv6T0RUb7kKB2+le4hmvXvnF7endhwYG3Mv3ZV5lyNUfwjSENJAJleRrWnIQYqJ+4pmjs5zvMOgJEo2i7i6sXxxLcWxuxbDt6XXpqbrgaM6P/cXb/14Fy9qSXA6UKAhoj+Pbvj8D6PV13gl7V8N2uZ71Ds2bOnbCyr6i/t0ksvjfPPP39EjzdUu3btive85z1H9DEAOLoqffkt7WzaWtNR2H+qUhewrA3sNw4sFR5qgJLs2Zf1ZdwXYRLpvc6SpeKl+0ZmNYdJh4HJB+OG1LyHagDC+NGwelNRyPbUuisjPl58TPJvmQ7ckn/j9DYFScVDsvdfuqp5sK65hWXDA3sEds+tjfrIXk5e6MYb/YFf1nNMh14ARir9ntP/A/rk2LhueeF2EvZFlO+Nnb6dNP9If95P3u+659YW3isPR9UFfvX19WVjWcHdUGSdN23atEHPmT17dsyePXtEjzdUO3fuPKLzA3D0JW/gEeVL3ZI3/8I+fxn76RU2zR84NmvPtM3rroy4fnlmkJIea1q2KjZX2CNECENauptpsodbUqm3p+HgMaXSH2hLQ+qhVH4xfiSvXU3LVkV9lFQ2xMFOvRHZHcEjDr5eFe4rOb70OdB41aai+0q/MCV79JXKet1rWVe8sXpyfvLDCgAMR+n71pSu1OqFpYuL3gt7U03Mkh9Ok/fV1/6kNlo/3H/fnC9fF/Pi4OenpDlW8vm/dP/koaq6wO/kk08uG3v++edHNFfWea973etGNBcADCar++2hgo7SX/yS4KS1pqMQ0pV+AEhXYaUrB9OPlf6Sm/4SXroMWBhDRBRV6qXN+vHBfY+zmnlElDehicgOvxm/kteIhihfmpQo6gSe8RqSLFNK/9snkufAnC9fFxER85beE11LF5c1H0of21rTEd0ZVXql2xFsu3pxTOkq/zsVfiABgGFKf3ZJ9olNusAXmkYN3J9uEpW8n+786O6IiJh+4/ER0f+elGxZUVoMcLiqLvDLqq7btm3biObavn172VhjY+OI5gKAwQwWbFTaM6/SOVnL2BJD2ZuqdElxpZBPGENE8RLerEq99HMxvfQyYmihtudZdUhXGKelqw/mX1vc+TazI3jJucl9Jz56XGE8aWyUnJMswS1tBBNR/qWoN10VsfqcaLxqYA+/1DLijYO8hgLAoSQ/gm5rr4nGdQeKtr/onlsbe64++KNUEgJuXHdl9LY1x09/e6C53u+U78Wdvt3SvjJi4L13pO9ZVRf4zZs3r2zsqaeeGtFcWUHhaaedNqK5AGCkRhJ4DOdNfyh79OmaSiXppeBZlXrpKtHkF+10JVfp8y+9hN3zrTqkl+Zm/WCQ/PtnVdMl96eX0A7lh4XSzcsTyZLiQgfEyO4EHBGF7ofp85Lr2XG1ffwAGJkd5/Zvb9K4bl9sa6+Jeev7x5NqvuT9MP35KNkXOQkL6waOT37kKv2MlV5hccws6T3zzDPLxu67774RzXXvvfeWjZ111lkjmgsAhmu4gcdQji/9wFBpjnSzjwjVfFSWfNBMfyCNiKLllll7sZWGyKW/XnvOVY+kqUZ99C+TTb9+pJcuNazeFK2rD25Wng4D67fuK1QAThkI6wp7Gq3vjBgY61q6OOq37susVi407GhrPlhNUdLgKNljMnnsStWmjVdtipaBZekq/QAYjt6TXo2IiJ0f3R3zfuffo7etudA9ftvVi2PLiv4fuFq27is04Chdrtvb1hx7Gsr39o44+L7UWtMRkbEdxlBN6Ovr6xvx2WPklFNOiaeffrpwe/LkyfHSSy9FXV3pTh+DmzNnTlGV3ymnnJK5zPdo27lzZzQ0NBSNdXV1xfTp08foigAYj0o3x09+JczqBJy1FDOifE+urA33IeJg8JeENkljhXRDhsH2gEw3/zjcTag5+tINW9JdviMOvm4kexil70tkdelNxist/c368aK0Uq/SY5VeW+nzzZJeAEYi/X6WNDBL72ucln6/SppcFSr81ncWOvVWepzS967hZkVVV+EXEfGud70rvv71rxdu9/T0xJ133hkXXHDBkOd49NFHy5b0Dud8ABhrpR0x0/uHRERZFV+lL9Slc0Kpoi6rVx+szko3kIkoD4zT3aVLf5YVtlSPZBlS8hUj63UiWTqbFb5F9IeB+975UkQUh3bp50mprMfpnltb6DKeVRWRHp9/7apIdudOLz3vWro4GkZpQ3QAji0bDqwt2rO28ar+H7rSXevTkvee5L6kwm/b1YvjtA9vivj97MeptHfucNQc1tlj5OKLLy4bu+GGG4Y1R9bxv/d7vzfiawKAsdDb1ly2X1pi47ori5YEbDiwtmhJZoSAj6FJnkcbDqyNxqs2FT7cllaRJv8lS1uSDq1dSxfbI7KKpfcUSgK60pBt47orM7/oJFUQ9Vv3xU9/+y/ip7/9F8VLlQaO6Z5bG91zawv7BCaSxyrs55daMlzaObr0x4xkSVXp/ZuvX+65CMCIbVmxvPBf8t6Y9R6YyPqBasuKwd+LhtKI71Cqcknvq6++Gm9605uKlt9OnDgxHnjggfjVX/3VQ56/Y8eOeMtb3hIvvfRSYezUU0+NJ554Impqxj4DtaQXgKFKV+1VWrqbtUSu0n3pXxN9ISaR3vuxtDKv0vLc9HOt0nPQc6w6pCvw0pVyWdV82wYqQGf9eF/RUu9kI/OIKPtSVLrNQOl+R+nGMcl48pxrWrYq6reWL6VKX2cyRzpgrN+6T5UpAEfUUD/rDPW4Y2JJ73HHHRef+tSn4hOf+ERhbP/+/XHJJZfEpk2bBt3Lr6+vL5YsWVIU9kVE/Omf/um4CPsAYDgqdd9Nf0FPd8csvS/9hVr4QiVJZWi6u1wSlqS7oFYKA9PPuaxOvox/lYK+9L5EXUsPblSePja9kXnEwedKulNhom59Z3/n3dQWBYXgbuniouPmfPm6aNx6oGgD9KQiomnZqmiI7OeZsA+A8eRIfQavygq/iP6A76yzzopHHnmkaPwd73hH/OM//mPU19eXndPT0xOXXHJJ/J//83+KxpuamuK+++6L44477khe8pCp8ANgJA61gX1pdV+pSt1UIeu5U1qlN//aVYVNq+vWd2rKkTNJJV3pfkKllX/p2xGpzrqpcwdrLJQOA9NBXkQUPXZWJWBaUQfgyP5BxPMTgGpyTFT4RfQv4f37v//7+PVf//XYs2dPYfzOO++MN73pTfHRj340zjvvvHj9618fzz//fHR2dsYNN9wQO3bsKJrnhBNOiL//+78fN2EfAIy25ItuVuiXFeoJ+ig12HMiCVUaU8eVVvtlVZl6nlWPlvaVsTkV6KVDtNJ/29J9Q9PBX9a/eUv7yqiLKCz5rY/iph5lxw6EyVmVgGmn3vz5OG39/UWhYm9bc+w4tza2eO4BcAyo2sAvIuKMM86I7373u/E7v/M7sXfv3sL4rl27YuXKlbFy5cpBzo6YNGlS/MM//EO89a1vPdKXCgBHVEv7ythYUoGVdBDbsmJ5WcCS9QVZCMOhDFblV9owJuJgtWgSzCTneJ5Vl+0feDV+5Z//In46UDWXDnYb1ndGy0D1XvrftXQpbWk1XbpqL2vZbVdq+W6ie25tNKTO6Z5bm3m9Le0r48l1fx7x+9ndylt+bEkvAPlX9ZvWvfvd744f/vCHMWfOnGGdd+qpp8Zdd90VF1xwwZG5MAA4ikqXykX0b5o/68f7Ki7hLWVPNUYiXUGamPPl6wpBYFZwM9TnJOND7VOT46MLNkXTslVFy2qT0CwJ+9Jdc0sDtaZlq4q67ab3BEyW/CbLfre110T91oOvXS3tK/vDxYFqvqSr756Gyg1AWms6omnZqkK38uSxLeUF4FhR1RV+iV//9V+Pn/70p/GlL30pvvKVr8TPfvazisfOnTs3rrjiivijP/qjmDx58lG8SgA4ctKVVknwkt7/KvmCm67MiSj+Ul66NxekVVqWW7oXW0REY1tzYWzz9csjUuc1LVsV9QONPFT6jX+tNR2x5cDaOPXmnqhtKP5RobetOTamgr7BXkM2X788WlLddJOlvKWvSV1LF8eJj0bUrb+ncGzpvEnw11rTEbFiedmcyTzJ9UeU7+HnuQdA3uUi8IuImDJlSnzyk5+MT37yk/HYY4/FAw88EM8880zs3r07jj/++HjDG94QZ511VixYsGCsLxUARl1pxVRWCJO+b2PGl19fgBlM0hwhq4FCYijPoT0NEQ2ry/dnY3xK9umb8s69Uds1uTAWcfC1ZVt7/6KhbVcvHnR/0OQ5k+zxmK7G6x4I6F5686vRuO7AoM+PQkOOQe6LGHjODiwlTyRdhQEg73IT+KUtWLBAsAfAMSVrb7WsP1f6khyhOy+DS547lfZiS0tC5KzKr1mp6j+qw753vhQR/WFcRPGPCK01HTEvIuLj/fuFzh/CfHXrO4vm6G1rjoaBYK4hyl+DkoAu3ZG3e+ni2FxyXLKnX936zthc8pqY/mGjtaajqOoUAPIol4EfAByLKnXfTVf7ZXW+LO3iC6XmX7sqGjOeOxHlDToq3Vep4pTxraV9ZdTOPbGwjDYtvRw3+XcdrANu6XMguZ2u9osoXj7eWtMRkVqem7xOJXv8pa9p8/XL+7v9ZlTwCZkBONZM6Ovr6xvri6DYzp07o6GhoWisq6srpk+fPkZXBEA1qLQv1WD7VVWqCoRSyXLedIOY0r0gVYnmT2tNR2y7enFM6SrZOy8lvU9opedAuqFH+oeH0tenZMltQ2opblJFmq4uzXpda2lfWXh+di1dXLbnZNZ1AUC1GG5WVPVdegGAfr1tzUVfqhNZX4qTjpXJEjdfgjmUrO6mSQfUpANr19LF0bV0cebzME2VX/XobWuO3pP6l/ImS2tLXzOSJbpJKJy1zLtufWfhv3Ro2LV0cVEH33TFcaXQML2XZNb1RkRRl9/0fOlOwQCQZ5b0AkBOlIYxlWQt64VDyer0nCzFrB8IWbK6+KZvJ+d7/lWPHefWxlMfX37IHxOSSsAtKyrvjZeEcS3tKwvdfetjYElvqqLvpTe/WvT8Sfbm23z98sJekjvOrS17jmVVDpZWoXruAXCsEPgBQM6k97lKJNVZ869dFY1jeG1Uv9JGHMNdRk51abxqU8RAiJdU8iX/tvOvXdV/fxw6SEu/JvW2NUdrTUfUpc5LV+w99fFPRnz84LnJ8t6WrfsKoV5L+8rMx0yC6ERWN2kAOBYI/AAg59JfiresWF748g4jkVVJWhrwZTWNEQBWpySo25jx7zfc15NKTTuGck5isOXiWY1hKu0VCAB5p2nHOKRpBwCQN5omVLeswOxwQrR0JXKlOUqr9Y7k9QDAeDfcrEjgNw4J/AAAyDPhHAAMjy69AADAuCbsA4AjS+AHAAAAADki8AMAAACAHBH4AQAAAECOCPwAAAAAIEcEfgAAAACQIwI/AAAAAMgRgR8AAAAA5IjA7xjVWtMRrTUdhzwGAAAAgOoycawvgLGx4cDaUTkGAAAAgPFFhd8xLquKL139p8oPAAAAoLoI/ChIwr0NB9YWVfcJ/QAAAACqh8CvihyJ4G2wZbvpABAAAACA6iDwqyJJ8Fap4UalQHCoQWFpsJdU+qnwAwAAAKgeAr8qlVV119vWHC3tK4d0flZomN63L33fE99aeBhXCgAAAMDRJPCrEhc3/GHhz71tzWX3t9Z0RN36zti47sqi8fnXrsqs0kuq99IBYXJcurKvadmqePL3/7zwGAAAAACMbwK/KvHdrjWFP5eGehGV99lrvGrToPMmc1Xar2/z9csLoaC9/AAAAADGv4ljfQEM3/xrV8WWFQeDuLr1nUX7+/W2NUfd+s6IOFgNeKiwLl3dlw7/Wms6ou5I/UUAAAAAGHUq/KpIEsQlVXt16zuLwr4NB9bGtqsXx8Z1V0ZvW3Mh+KvU5KNUVii44cDazCXEAAAAAIxPKvyqSBLIJc05Nqaq+pLxWQPHJhV+6cq9SlraVxYtE04HfK01HYXHAQAAAGD8m9DX19c31hdBsZ07d0ZDQ0PRWFdXV0yfPr1wu7Qbb+m+funlvhGHXtKbrgBMAr+kerA0EAQAAADg6BlKVpRmSW8VKg3gklAvuS+iPwBMjtlwYG00LVs16JzpZcHJucI+AAAAgOpjSW8V2rjuyqJmGtuuXhwRlZffttZ0xOYKFX7Jct+mZauiITUWcTD8O9SSYAAAAADGDxV+Vai0AceWFcsL3Xmz7q80FnGwsm/z9cuLOvSWHjOUph8AAAAAjD2BXxXqbWsuqrhLwr6N666MlvaVRfdtOLC28F+W0r0AS+dNjtGpFwAAAKA6WNJbZZKlvC3tKwtLeiMO7uNXugQ3XZmXFfol+/OlO/3GQLiXzBVtzUX7BAIAAAAwfqnwqzIbDqyNbVcvju65tUXVe4Mtx03/b5akKjCpEqxb31nW3be3rfmQjT8AAAAAGHsq/KrQrB/vK7qddNLNCvdKKwHT4xERO86tjS0DVX516zv7Q72liwvHJY1AWms6isYBAAAAGJ8EflUoqb5LKvrqongfv6HYcW5tREQ0XrUpYqDpR1mV4PXLC3v3lS4R1rUXAAAAYHyypLdKlTbj2HBgbVkQGBEV997rPenV/v8O0Yyjbn1nUWOPwRqAAAAAADD2JvT19fWN9UVQbOfOndHQ0FA01tXVFdOnT4+Ig9V8pWFeVofd3rbm2HFubTRetSkzqEvCvEqVgS3tK6NufWehelCFHwAAAMDRdaisqJQlvVUm3YE3Ceu2f+DVePL3/7zovq6B/fb2NAws2y2RnFu3vrNoKW+lZcEb113Zv1egoA8AAABgXBP4VaHWmo7oWro46gduz77luGi5pTiMa1jdH/L1tjWXhXhJI46IiPnXrirbvy+9JHjjgbXRtGxVoTEIAAAAAOObwK/KJBV8TctWFcay9ulLKvw2X7+8aHz+tatiVltzIdSbVbKHX1aoV7+1vytwspS40nEAAAAAjD2BX5VJKu1Kg7yi7rpRHvQlYd2WgaW5SXCXXtLb29YcTctWFaoD0/v0JQFf6eMAAAAAML4I/KpMOnhL79O3ZSC0S0K6wYK5dPONRLopR2lDjrr1nQer+w7R1RcAAACAsSXwq1JlHXlXLC9U6pUGdklH37r1nTH/2lWxZeDYwrnRH+q1tK+MSC33TR4n3RHYUl4AAACA8U3gV8WSYC7dcKO0S+/m65cXOvL2tjXHlhX9S30LXXoH/lwXByv/0st9W2s6om5gvtJlwgAAAACMPwK/KlJavVdapZcEda01HbE5dd/GA2sLAV+itGtvEhpuu3pxNF5VvIdfa01HbL5+eWb1IAAAAADji8AvB0pDuNLbLe0ro259Z6HqL610L79ZP96XOX9SBQgAAADA+CbwqyK9A/vrZVX2pRtulIZ4Gwc68yZLctP397Y1R/fc2ph/7aqYNbBXXzowTFcGqu4DAAAAGP8m9PX19Y31RVBs586d0dDQUDTW1dUV06dPLxor3cOvrJFHRNF+fOnj0kFe0pAjfXwynm7YkQSLAAAAABw9Q82KEir8cqBp2arYXNK4Iwn/mpativqt+4qCvIiDIV/6uJfe/Go0rjtQCPVaazqKAr6W9pVlzUEAAAAAGF8EflXi4oY/jI19t0bEwSYb267u35OvdN+99JLdhjhYAZgeTy8L3nBgbaEpR9fSxWWhXnJM6XJfAAAAAMYfgV+V+G7XmrKxdDfdpKnGYIFc19LFhX380pKgL+ns21Ryf1IdKOwDAAAAGP8EflWmtOFGWnK7pX1lYSlu0qG3kqT6r2H1pmiKiIbVm6J+YJ7ksTZmNAIR/gEAAACMTwK/KlPUoTcVzKWr+5KOvREHw7p0o46spb0R/WFfRPH+foM1BgEAAABg/BH4VZF05V5EFP052aOvNJhLlvq21nTExoHQb2NJaJfu5Lvt6sUxpevgnL1tzYUqQWEfAAAAwPhXM9YXwNClK/fSjTiSgC8J5LqWLi46Z8OBtYXgLh0SJurWd0b33NroWrq4qAFIumIQAAAAgOqgwq9KZHXpTVfcpZfpbi5ZiptIluqWLs0tm2cgMEyWCidVggAAAACMfyr8qsR3u9YUltiml+sm+/KVBoDJ7Zb2lf3LeVNVeunj0vv6JUFgw+pNhf38NhxYG91zay3nBQAAAKgSKvyqTN36zkIwV7qHX6VQbtvVizPHI4qX6xb26ys5pn7rvsJjJASAAAAAAOOTwK9KXNzwhxHvPq8Q0GUty0036EgacXTPrY0tK5YXjokob/6Rvi8d6iVBYffc2qJjAAAAABi/JvT19fWN9UVQbOfOndHQ0FA01tXVFR+c8fGIKA/esvbrS8Z625qLOuw2LVsVm69fXtb8Izk24mDVX9KdN+sxAQAAADg6KmVF06dPzzxehV+VuLjhD6NuwqSiSrxCQJdRnZcO79LjDas3RVy/PLNpR1I1mBxfF8XdgIV+AAAAAOOfwK+KJZV4TctWRcPqTUVhXbIMt1RWaJc07kiCw/TS3+R28mcAAAAAxjeBXxXpbWsuVNqlw71YurjQcGNjRqBXGvKl9/BrrekoNOlIgsOupf17923OOAYAAACA8a1mrC+Aoen9rTOL9taLiELw17B6U1nH3vR/iZb2ldFa01F0bDoMTLrxNqze1L/0FwAAAICqI/CrErfe8oloaV9Z2JtvsP36kgq9iOL9+dJVfWnJMXXrO4vObVq26gj8TQAAAAA4kizprSJJx9xEurnGnC9fF/OW3hO9bc1F1XnpbryVlvy2tK+MGOjmGxGx7er+0K/xqoMNPuzhBwAAAFAdVPhViYsb/rDw52Qpb/JfRETjugPROxDapZfpdg3s75dU+JVW97XWdETd+s6iJh9Tuvr/S+6P6K8eFPoBAAAAjH8T+vr6+sb6Iii2c+fOaGhoKBr79d/6bEy94z8iorwJR7qKL7mddUx66W+WSvNkzQcAAADA0ZGVFXV1dcX06dMzj7ekt0rcessn4uJLvlFUiVca0CVjGw6sjfnX9u+/13jVpvKuviXnZv05PW9SIVj6WAAAAACMPwK/KlK3vjMaIiKuXx4RB5t1RBRX9bXWdMSs1H1pyTHpY5PbpZV8ye3SZcIAAAAAjF/28KsSFzf8YWw4sLaoS+/GdVcWlucm4+l9/ZLbWX8uHUvCvfQ+fb0loWHp/n8AAAAAjD8q/KrEd7vWRGtNR3QvXVwI3pJwLgn+kvGupYujoaSjb3J8EhCWLtEtreyLiKKOvq01HdGVemwVfwAAAADjk6Yd49BgTTvSjTey9tsr1bRsVWweWAKcDvxKZTXoaFq2KhpWb9K4AwAAAGAMDbdphyW9VaLuBw/2/+/6zqgbqN5Llu8mS3GTADAdBDas3lT4c2lImCgN+1prOsoafAj7AAAAAKqDJb1V4ud/1hxz/+qhwu3SEC5prFHWpbfC8WlZjToSLVv3Rd36zmit6Sja069SpSAAAAAAY0uFX5W45//3sehta87ccy+5nQ7lGq/aFI1XbYretuaiar3S2xHl+/klFX5Ny1YVgr2kW2/33Nronlt7ZP6SAAAAABw2FX5VJr3sNq2oKq99ZaHhRvr4pMlHpYq+suq+kuN725oL+wECAAAAMD4J/KrExQ1/GHUTJkVE9vLcQ+2zl+zzV1fSvTddGbix5Ji6KO8EnL4NAAAAwPgj8KsS3+1aEx+c8fHobWsuC+0S6Wq+ZGlvXWo8vc9f4f71nYU5W2s6oi4132ABokYeAAAAAOOTwK/KJKFdqfSy3N625kIFXmmzjfQ8yXlNy1bF5tQy4dIlv6VjWbcBAAAAGB8EflUkaxlv1njpctvSisB0ANha0xENEdG6elNRaFg4buDYZKlv99LFh/8XAQAAAOCImdDX19c31hdBseeffz5mzpxZNPbTn/40Tj755Ijo38/vu11riu4vHbu44Q8jon8p8G9++qvxw7/8WOZ5ifd84H9FRMStt3yi6Nz/f3t3HhT1fT5w/AER5FAOBQSPBdR6IR5BUeMRj5R4ZlrrkdGQSRqjjmljtUnTQ22TqU5aJ6mJ02jVWI3WTLR1atUYUXAs1lA1Koq3IIo3WfDgPj6/Pzr6c/kusPtl2S+7+37N8AcP+3w+X9iHB/bZ3e/36bWejgEAAAAAAMA5CgoKpFevXhax27dvS2RkpNXbM/Brhs6dO6e5EwEAAAAAAIDHzp49Kz179rT6NW8nHwsAAAAAAACAJsTADwAAAAAAAHAjDPwAAAAAAAAAN8I5/JqhqqoquXTpkkUsLCxMvL2Zzzbkm2++kcmTJz/5fOfOnTJ48GADjwhwPOocnoA6hyegzuEJqHN4CmodTa2mpkbMZrNFrFu3buLj42P19tajMJSPj0+dJ11E/UJDQzWfh4eHG3Q0QNOgzuEJqHN4AuocnoA6h6eg1uEMdV2R1xpeMgYAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBvxMfoAAEfq3LmzLF261OJzwN1Q5/AE1Dk8AXUOT0Cdw1NQ62huvJRSyuiDAAAAAAAAAOAYvKUXAAAAAAAAcCMM/AAAAAAAAAA3wsAPAAAAAAAAcCMM/AAAAAAAAAA3wlV64fGUUnLmzBnJysqSW7duSVlZmQQFBYnJZJJnnnmGqyuhWTCyTqurq+XEiROSnZ0td+/elYqKCmndurXExcXJoEGDJCIiosn2BpyFOvccFRUVcuHCBTl37pzcu3dP7t+/L35+fhIaGiqRkZEycODAJr2/6edwFqNr3UjUuue4d++eXLlyRfLz8+XOnTtSXFz85P4ODg6WmJgY6d+/vwQHBzfJ/vR0NGsKcLKRI0cqEXHIx8OHD3UfR2FhoVq8eLGKjo6ud4++ffuqdevWqaqqKgf+FOBqcnJy1Jdffql+8YtfqLFjx6rQ0FBNrbzyyisO39fIOr1x44ZasGCBCgsLq3NfLy8vNWzYMLV9+3aH7QvjOKvOTSaTQ/4GtG3bttHHQp27v8rKSrV//3717rvvqqSkJOXj49NgbcXFxaklS5ao/Px8hx0H/RxNzahap6fDWcrLy9XevXvVkiVL1PPPP68iIyNtqi0vLy/Vt29f9eGHH6qCggKHHAs9Ha6AgR+crjkM/Hbt2qXCw8Pt2qt///7q0qVLDv5poLlKT0+vd+hh7cPRAz8j63T9+vUqKCjIrr3Hjh2r7ty544DvHM5iVJ03lweH1Ll7y8/PVz/+8Y9V27ZtdddYy5Yt1W9/+1tVWVnZqGOhn6MpGV3r9HQ4y+nTpxtdZyEhIerTTz9VNTU1uo+Dng5XwcAPTmf0wG/dunXKy8tL135hYWHqxIkTjv+hoNl58cUX7a4PRw78jKzTJUuW6P6djImJUXl5eQ77OaBpGVXnzeHBIXXu/lJTUx32/8bgwYPV/fv3dR0H/RxNzehap6fDWRwx8Hv8MWvWLF0Dbno6XAnn8IPhTCaTJCYm6sr18bGvhPft2ydvvPGGKKUs4lFRUTJnzhxJTEyUyMhIyc/Pl4yMDFm/fr0UFRU9uZ3ZbJbx48fLt99+K+3bt9d1zEBDjKzT9evXy3vvvaeJd+nSRebOnSt9+vSRsLAwycvLk9TUVPn888+ltLT0ye2uXr0q48ePl8zMTAkMDLTvG4fH6tmzp/Tq1cvuvDZt2ujajzr3bF27dpXhw4fL0KFDJTo6WiIiIqSqqkry8/Pl0KFDsnnzZiksLLTI+eabb2T8+PGyf/9+adWqlc170c9hJGfW+tPo6XCW9u3by9ChQyU+Pl66d+8uUVFR0qZNG/Hy8pKioiK5cOGC/Pvf/5Z//vOfUlJSYpG7efNmCQgIkDVr1ti8Hz0dLsfYeSM8Ue1X+DXFec+sKSoqsnqeh+nTp6vS0lKrOWaz2eorEidMmOCUY4ZxrL3yydvbW/Xq1Uu9/PLL6te//rXm646oZSPrNCcnR/n7+2vW+dnPfqaqq6ut5uTl5an4+HhNzvz58+3+3uF8RtV57VeDLF26tNFr2oo69xxPv+qpXbt26q233lKnTp1qMO/hw4dq7ty5Vl8h8bvf/c7m/enncBaja52eDmc5d+6cGjNmjPrTn/5k19tj7927p2bNmmW11g8cOGDTGvR0uCIGfnA6owZ+ixYt0jS8SZMmNXj+hvLycjVgwABN7t69e51y3DDGlClTVO/evVVKSopauXKlysjIUI8ePXry9dzc3CYZhBhZp1OmTNHkv/nmmw3mmc1m1bFjR83QKDs72+a9YQyj6tzIB4fUuedITU1VERERasWKFaqkpMTufGtvnfL391d37961KZ9+Dmcxutbp6XAVb7zxhtW+bAt6OlwRAz84nREDP7PZrAIDAy32DQ4OVrdu3bIpPysrS3OlsxEjRjTxUcNIdT1b9lhTDEKMrNPz589rzkcSGxuriouLbcrfs2eP5ueRkpJiUy6MY0SdK2Xcg0Pq3LN89913FgNse9XU1KikpCTNfb527doGc+nncCYja10pejpcx8OHD1VISIjFfe7r69vgoJyeDlflLYAH+OKLL6S4uNgi9tprr9l87oQ+ffrIxIkTLWKHDh2SixcvOuwY0bx4ezu/PRpZp5999pnmfCQLFiyQgIAAm/YeN26c9OvXzyL25Zdfyv37923KhzGMqHMjUeeeJSwsrFHnKfLy8pIFCxZo4qmpqQ3m0s/hTEbWupGoddgrKChIkpOTLWIVFRVy/fr1evPo6XBVnvWfPjzWtm3bNLE5c+bYtYa121tbF9DLyDrdvn27xeetWrWSlJSURu1dVlYmu3btsmsNoClR57DX8OHDNbGGHhiK0M/hevTWupGodegRExOjiRUUFNSbQ0+Hq2LgB7dXUlIihw8ftojFxsZK9+7d7Vpn1KhR4ufnZxH7+uuvG318gIixdXr58mXJycmxiA0ZMkRCQkLs2vuFF17QxPgdQXNBnUOPdu3aaWLfffddvTn0c7giPbVuJGodepWVlWliwcHBdd6eng5XxsAPbu/48eNSUVFhERs2bJjd6/j5+cnAgQMtYkePHpWqqqpGHR8gYmydHjlyRBPTs3dMTIx07NixwbUBI1Dn0OP27duaWFhYWL059HO4Ij21biRqHXqdPHnS4nN/f3/p1q1bnbenp8OVMfCD2zt16pQm9swzz+haKzEx0eLzsrIyuXDhgq61gKcZWadNufeVK1fk0aNHutYCHIk6hx7W6qb2gyZbcujnaO701LqRqHXokZWVJRkZGRaxiRMniq+vb5059HS4MgZ+MNzFixdlwYIFkpSUJNHR0eLn5yetW7eWmJgYGThwoMyfP1+2bt0qhYWFuta/dOmSJmbt3A22MJlMmtjly5d1rQU8zcg6bcq9lVJy5coVXWvBcxw9elTmzZsn/fv3l8jISPH19ZWQkBCJi4uToUOHysKFC2XHjh1SUlKiew/qHHps3LhRExs7dmy9OfRzuCI9tV4Xejqao5s3b8qMGTOkurr6SczX11d+85vf1JtHT4cr8zH6AIAjR45oXlJcUVEhjx49kry8PDl27Jj8+c9/lsDAQHnttdfknXfesesZx2vXrmli1pqtLTp37qyJ5eXl6VoLeJqRdeqMvfv27atrPXiGPXv2aGL379+X+/fvS25urhw5ckQ++ugjCQsLk3nz5snChQvtfqsZdQ57Xbx4Uf71r39ZxFq0aCGTJ0+uN49+Dlejt9brQk9Hc3L79m3ZsmWLLFu2TMxms8XXPv74Y0lISKg3n54OV8bADy6juLhYPvnkE/n8889l/fr18sMf/tCmPGtXXYqMjNR1DNbymvMJjeE6jKzT2ns/fibeGXsD9jCbzfL73/9eNm7cKFu2bJERI0bYnEudwx5KKZk3b55UVlZaxGfOnClRUVH15tLP4UoaU+uNRU+HI9y7d0/mzZtnEauurpYHDx5Ibm6u5ObmanKCg4Pl008/lZdeeqnB9enpcGUM/NAsdO7cWUwmkwQHB0t1dbUUFhbK+fPnpaioSHPboqIimTJliqxYsUIWLVrU4NrW1vD399d1nNby9L7VGHiakXVae2+9++rZGxAR8fb2lri4OOnQoYMEBwdLeXm5mM1myc7OtvqWr/z8fBk9erT87W9/k2nTptm0B3UOe6xYsULS0tIsYkFBQbJs2bIGc+nncCWNqfW60NPhTMXFxfL3v//dptuaTCaZOXOmvPXWWxIREWFTDj0droyBHwwRHh4u06ZNkxdeeEFGjBghbdq00dympqZGTp06JatXr5aNGzdKeXm5xdfffvtt6dSpU4P/GFi79HqrVq10Hbe1Rln7uAA9jKzT2nvr3VfP3vBcnTt3lunTp0tycrIMHTrUau1UVVVJZmamrFq1SrZt22Zx3p3q6mpJSUmR6Ohom65YR53DVvv375df/vKXmvhHH30kHTp0aDCffg5X0dhafxo9Hc2dj4+PDBgwQHr06GHXq+To6XBlXLQDTvfBBx9Ifn6+rFq1SiZOnGh12Cfyv2cH+/fvL2vWrJFjx45Jz549Lb6ulJI5c+Y0+MyEtUud13clpvr4+flpYrXfAgHoYWSd1t5b77569oZn2rJli+Tm5sof/vAHGTNmTJ3PWPv4+Mizzz4rW7dulbS0NM0D0PLycnn11VdtqjHqHLbIzs6WqVOnWgwiREReeuklef31121ag34OV+CIWn+Mng5XUFVVJTt27JCUlBSJjY2Vf/zjHzbn1UZPh6tg4AenS0pKsrtZxcfHy8GDByUuLs4iXlRUJMuXL68318dH+0LWiooKu/Z/zNqzIC1bttS1FvA0I+u09t5699WzNzzTs88+K97e9v0LMmLECElLS5O2bdtaxC9fvixr165tMJ86R0OuXr0qycnJmrdQJSUlybp162xeh36O5s5Rtf4YPR1GiomJEaWUxUdZWZncvn1bDh06JB988IH07t3bIufmzZsyZcoU+dWvftXg+vR0uDIGfnAZERERsnnzZk18w4YNUlNTU2eetZc+W3tpti1KS0s1MWvPlgD2MrJOa++td189ewP2+N73vierVq3SxG15gEqdoz43b96UMWPGyI0bNyziffv2la+++koCAgJsXot+jubMkbXeWPR0NBU/Pz+JjIyU4cOHyzvvvCNnzpyRTZs2ad7Ku3z5clm5cmW9a9HT4coY+MGlDBkyRJKTky1iBQUFcuzYsTpzrJ2jwVrDs4W1vNDQUF1rAU8zsk5r7613Xz17A/aaPn269OjRwyJ28uRJuXPnTr151DnqcufOHRk9erTk5ORYxHv16iWpqal237f0czRXjq51R6Cnw1lefvllSU1NldatW1vE3333Xbl27VqdefR0uDIGfnA5kydP1sQyMzPrvH27du00sYb+iaiLtbzab0UA9DCyTmvvXVFRYfWKZE2xN2AvLy8vmTRpkkVMKSX//e9/682jzmFNQUGBjB07Vi5cuGAR7969u6SlpUl4eLjda9LP0Rw1Ra07Aj0dzpSYmCjvv/++RaysrEw+/vjjOnPo6XBlDPzgchISEjSx+ppu586dNbG8vDxde1t79sdkMulaC3iakXXK7whcjb1/B0Soc2iZzWZ5/vnn5cyZMxbxrl27SlpamkRGRupal36O5qapat1R6OlwptmzZ2suLLN79+46b09Phytj4AeXY+0ZyIKCgjpv361bN03s6tWruva21mC7du2qay3gaUbWaVPu7eXlJV26dNG1FlAXe/8OiFDnsFRUVCTf//735eTJkxbx2NhYSUtLk+joaN1r08/RnDRlrTsKPR3OFBAQIIMGDbKInT9/3upFLUTo6XBtDPzgckpKSjSx2s/SPK1fv36a2PHjx3XtXftcgf7+/tK9e3ddawFPM7JOHbl37byuXbtKUFCQrrWAutj7d0CEOsf/e/DggSQnJ2vuR5PJJOnp6dKpU6dGrU8/R3PR1LXuKPR0OFtUVJQmZjabrd6Wng5XxsAPLic/P18Ti4iIqPP2AwYMEF9fX4tYRkaG3fuWl5fL0aNHLWKJiYlWL9UO2MvIOh0yZIgmpmfvq1evyvXr1y1igwcPtnsdoCH2/h0Qoc7xP48ePZJx48Zpzg/WqVMnSU9Pd8jbm+jnaA6cUeuOQk+Hs1kbMtc1/KKnw5Ux8IPLOXDggCZW3z8tAQEBMmzYMItYbm6u5qTFDUlPT9e81Lv2FYMBvYys0y5dukhcXJxF7MiRI3afFHjv3r2aGL8jaAr2/h0Qoc7xvwd4EyZMkP/85z8W8Q4dOkh6errExsY6ZB/6OYzmrFp3FHo6nK32laoDAgI0V+99+mv0dLgqBn5wKY8ePZIvvvjCIubl5SVjx46tN+9HP/qRJrZmzRq79rZ2+6lTp9q1BlAfI+u09t5lZWWyadMmu/b+y1/+YvF5q1atNFfeAxrr+vXr8vXXX1vEQkJCNOfjsYY691ylpaUyadIkOXTokEU8KipK0tLSHH4eI/o5jOLsWm8sejqcLTc3V3MBm/j4+Hpz6OlwWQpwIT//+c+ViFh8DB48uME8s9msAgMDLfJCQkLUrVu3bNo3KytL+fj4WOSPGDGisd8OXFhubq6mFl955ZVGrWlknZ4/f155eXlZ5MbGxqri4mKb8vfs2aP5eaSkpNiUi+arKeq8saZOnao5phkzZtiUS517prKyMpWcnKy579q3b6/Onz/fJHvSz2EEI2q9sejpcLaZM2dq7vfly5fXm0NPh6ti4Aenyc/PV5cuXdKd/8knn2ialYioPXv22JS/aNEiTe7kyZNVTU1NvXnl5eVqwIABmty9e/fq/l7g+ppqEGJknU6ZMkWT/5Of/KTBvMLCQtWpUyeLPG9vb5WdnW3z3mieHF3nFy9eVDdv3tSVW1NTo95++23N8Xh7e6vTp0/bvA517lnKy8vVhAkTNPd5RESEOnv2bJPuTT+HMxlR6/R0OFNubq769ttvG7XGH//4R029+Pv7q2vXrjWYS0+HK2LgB6dJT09XLVq0UNOnT1f79u1TlZWVNuVdv35dzZo1y+qwLzk52eb9CwsLVWRkpGaNGTNmqLKyMqs5ZrNZPffcc5qcCRMm2Lwv3FNTDfyMrNOcnBzl7++vWWfhwoWqurraas61a9dUfHy8Jmf+/Pl2f+9ofhxd5xs2bFB+fn5q9uzZKiMjo8F/kh87f/681VetiIiaM2eOXcdAnXuOyspK9YMf/EBzv4WHh6szZ840+f70cziLUbVOT4czpaenKxFREydOVDt27FClpaU25+bk5Fh9NamIqGXLltm0Bj0drshLKaUEcIKDBw/KqFGjnnweFhYmo0aNkoSEBOndu7eEhYVJcHCw1NTUiNlslrNnz8rBgwdl9+7dUlVVpVkvPj5eDh06JKGhoTYfw759+2TcuHFSU1NjEY+KipJ58+ZJYmKiREREyI0bNyQjI0PWrl2rOSlq+/bt5cSJE9K+fXv7fgBwOdbO1/FYSUmJfPXVVxYxk8kkiYmJdea8+eab8txzzzW4r5F1um7dOpk9e7Ym3rVrV5k7d64kJCRISEiIXL9+Xfbt2yebNm2S0tJSi9v27t1bMjMzJTAw0K69YQxn1vlf//pXefXVV598Hh0dLSNHjpSEhATp0aOHhIaGSps2baSiokLMZrOcOnVKDhw4IAcOHBBr/66MHj1a9uzZI35+fg18l5aoc8/w4YcfyqJFizTxfv36Nfo8Ztu3b7fpdvRzOINRtU5PhzPVfiwZFBQkw4cPl379+kmfPn0kPDxcgoODpWXLlvLw4UO5c+eOZGVlSXp6uhw+fNhqzb344ouybds2admypU3HQE+HyzF03AiP8vhZGUd8JCUlqRs3bug6jrVr12rOg2DrR2hoaKNfSg7X4ah6ffyxYcMGm/c2sk4XL16s+3s0mUzq6tWruveG8zmzzjds2OCwfSZMmKAePHig+/umzt3f0qVLHV7fjz/sQT9HUzOq1unpcCZHPpYUETVz5kxVUVFh93HQ0+FKuEovXEpgYKC8//77cvjwYYmOjta1xuuvvy47d+6Udu3a2ZXXr18/yczMlP79++vaF7CHkXX63nvvybp16yQoKMiuvDFjxkhmZqaYTCbdewMNCQsLk9WrV8uuXbukdevWutehzuEs9HOgbvR02MrHx8ch68TExMjOnTtl8+bNNr+y72n0dLgUoyeO8BzFxcVqx44dasGCBWrQoEEqICDApmcjfH191aBBg9TKlStVUVGRw46nsLBQLV68WEVFRdW7f0JCglq7dq2qqqpy2N5wDbbUpz0f9rzC7zEj6/TGjRvqpz/9qQoNDa1zXy8vLzV06FC1bds2h+0L53JmnRcWFqqtW7equXPnqn79+ilfX1+b1vT391cjR45Un332mV3n7LEFde6+mssr/B6jn6OpGFXr9HQ4W15enlq9erWaNm2aMplMNtdxeHi4mj59utq9e7fDeis9Ha6Ac/jBMDU1NXLlyhXJy8uT/Px8KSwslJKSEmnRooUEBwdLSEiIxMTEyIABA+w+l4c9lFKSlZUlp0+fllu3bklZWZkEBQVJ586dJTExkWdC0CwYWafV1dVy/Phxyc7Olrt370plZaW0bt1aYmNjJSkpSSIjI5tsb7i3yspKuXz58pO/Aw8ePJDS0lLx8fGRkJAQCQ0Nla5du0pCQoLDntmvC3UOZ6Gfw13R0+FsBQUFcunSJcnLy5N79+5JcXGxVFdXS+vWraVNmzYSEREhCQkJ0rFjxyY7Bno6mjMGfgAAAAAAAIAb4Rx+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4EQZ+AAAAAAAAgBth4AcAAAAAAAC4kf8DlMV8E3qEL+UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1500x1125 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_=hist2d(S_dat['ra_mock'],S_dat['dec_mock'],1000,norm=LogNorm())"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d637dd82",
   "metadata": {},
   "source": [
    "# read in density sample.\n",
    "with fits.open(hsc_dir+D_dat_name) as f:\n",
    "    data=f[1].data\n",
    "    print(data)\n",
    "    D_dat=Table(data)\n",
    "D_dat.sort('z')\n",
    "\n",
    "D_dat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "69b054ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "D_dat=S_dat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f4e4c745",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<TableColumns names=('tract_id','irow_in_original_fits_file','object_id','ra_mock','dec_mock','e1_mock','e2_mock','shear1_sim','shear2_sim','kappa','weight','z_source_mock','z','ra','dec','e1','e2')>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_dat.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "atlantic-italic",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'e1' in S_dat.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "another-reunion",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function to write data files in format that can be read by corr_pc\n",
    "\"\"\"\n",
    "def write_data(data,file_name,file_dir,wt_col=None):\n",
    "    file_name=file_dir+file_name\n",
    "    np.savetxt(file_name+'_pos.dat',np.column_stack((data['ra'],data['dec']))) #ra,dec position of galaxies\n",
    "    np.savetxt(file_name+'_z.dat',data['z']) # z position of galaxies\n",
    "    if 'jk' in data.columns:\n",
    "        np.savetxt(file_name+'_jk.dat',data['jk'],fmt='%i') #jackknife is read as integer\n",
    "    if 'e1' in data.columns: #ellipticities. corr_pc takes e1, e2\n",
    "        np.savetxt(file_name+'_e.dat',np.column_stack((data['e1'],data['e2']*-1)))\n",
    "    if wt_col is not None: #wt to be applied to each galaxy/particle.\n",
    "        wt=data[wt_col]*1.\n",
    "        wt/=wt.mean()\n",
    "        np.savetxt(file_name+'_wt.dat',wt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "infectious-highlight",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_nameS=file_name+'_shape'\n",
    "write_data(data=S_dat,file_name=file_nameS,file_dir=file_dir,)#wt_col='mass')\n",
    "\n",
    "# file_nameD=file_name+'_density'\n",
    "# write_data(data=D_dat,file_name=file_nameD,file_dir=file_dir,wt_col='mass')\n",
    "file_nameD=file_nameS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "efficient-anniversary",
   "metadata": {},
   "outputs": [],
   "source": [
    "inp_gen(file_dir=file_dir,file_name=file_name,file_nameS=file_nameS,file_nameD=file_nameD,corr=1,n_shape=len(S_dat),n_density=len(D_dat),\n",
    "        njk=0,do_jk=0) \n",
    "#corr=0 will only do wgg like correlation, no shapes. set corr=1 for wg+ and wgg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fc726ce",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6875828a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "524d0a5e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "intel",
   "language": "python",
   "name": "intel"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}