From 1d0a7fc31f7ad8cff6deaf44808fe1ad94671b79 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 25 Oct 2021 16:42:48 +0100 Subject: [PATCH] Initial commit of xcorr (#22) * initial commit of xcorr * checks against CCL, fixed units in calculation of mag kernel * added sacc file handling * exposed defaults in xcorr initialize * removed dummy_pk * updated code style * added proposal for xcorr b parm * added docstring, made sacc load default behaviour * whitespace * added other parms to docstring, smaller b1 proposal * added s1 proposal * removed redundant chi setup option * changed to use soliket constants * changed theory to provider * rename cosmo to provider for clarity * changed chi result to dict * removed zeff and autoCMB for now * tidy up codestyle Co-authored-by: Ian Harrison --- .../dev/cross_correlation_simulatedata.ipynb | 486 +++++++++ soliket/__init__.py | 2 + soliket/constants.py | 1 + .../data/unwise_g-so_kappa.sim.sacc.fits | 932 ++++++++++++++++++ soliket/tests/test_runs.py | 6 +- soliket/tests/test_xcorr.py | 171 ++++ soliket/tests/test_xcorr.yaml | 33 + soliket/xcorr/XcorrLikelihood.yaml | 41 + soliket/xcorr/__init__.py | 1 + soliket/xcorr/data/clgg.txt | 3 + soliket/xcorr/data/clkg.txt | 3 + soliket/xcorr/data/dndz.txt | 401 ++++++++ soliket/xcorr/limber.py | 121 +++ soliket/xcorr/xcorr.py | 264 +++++ 14 files changed, 2463 insertions(+), 2 deletions(-) create mode 100644 notebooks/dev/cross_correlation_simulatedata.ipynb create mode 100644 soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits create mode 100644 soliket/tests/test_xcorr.py create mode 100644 soliket/tests/test_xcorr.yaml create mode 100644 soliket/xcorr/XcorrLikelihood.yaml create mode 100644 soliket/xcorr/__init__.py create mode 100755 soliket/xcorr/data/clgg.txt create mode 100755 soliket/xcorr/data/clkg.txt create mode 100755 soliket/xcorr/data/dndz.txt create mode 100644 soliket/xcorr/limber.py create mode 100644 soliket/xcorr/xcorr.py diff --git a/notebooks/dev/cross_correlation_simulatedata.ipynb b/notebooks/dev/cross_correlation_simulatedata.ipynb new file mode 100644 index 00000000..205217b5 --- /dev/null +++ b/notebooks/dev/cross_correlation_simulatedata.ipynb @@ -0,0 +1,486 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "296ce09a", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from astropy import units\n", + "import pyccl as ccl\n", + "import sacc\n", + "\n", + "import sys\n", + "sys.path.append('../../')\n", + "\n", + "from matplotlib import pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "id": "88cc94b4", + "metadata": {}, + "source": [ + "# Set up cosmology" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "5db24a17", + "metadata": {}, + "outputs": [], + "source": [ + "cosmo = ccl.Cosmology(Omega_c=0.25,\n", + " Omega_b=0.05,\n", + " h=0.7,\n", + " n_s=0.965,\n", + " A_s=2.11e-9,\n", + " Omega_k=0.0,\n", + " Neff=3.046,\n", + " matter_power_spectrum='linear')\n", + "b1_unwise = 1.0\n", + "s1_unwise = 0.4" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "c5d1b487", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "pyccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.7, n_s=0.965, sigma8=None, A_s=2.11e-09, Omega_k=0.0, Omega_g=None, Neff=3.046, w0=-1.0, wa=0.0, T_CMB=None, bcm_log10Mc=14.079181246047625, bcm_etab=0.5, bcm_ks=55.0, mu_0=0.0, sigma_0=0.0, c1_mg=1.0, c2_mg=1.0, lambda_mg=0.0, m_nu=0.0, m_nu_type=None, z_mg=None, df_mg=None, transfer_function='boltzmann_camb', matter_power_spectrum='linear', baryons_power_spectrum='nobaryons', mass_function='tinker10', halo_concentration='duffy2008', emulator_neutrinos='strict')" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cosmo" + ] + }, + { + "cell_type": "markdown", + "id": "418ed394", + "metadata": {}, + "source": [ + "# Set up binning" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "39d96e49", + "metadata": {}, + "outputs": [], + "source": [ + "ell_max = 600\n", + "n_ell = 20\n", + "delta_ell = ell_max // n_ell\n", + "\n", + "ells = (np.arange(n_ell) + 0.5) * delta_ell\n", + "\n", + "ells_win = np.arange(ell_max + 1)\n", + "wins = np.zeros([n_ell, len(ells_win)])\n", + "\n", + "for i in range(n_ell):\n", + " wins[i, i * delta_ell : (i + 1) * delta_ell] = 1.0\n", + " \n", + "Well = sacc.BandpowerWindow(ells_win, wins.T)" + ] + }, + { + "cell_type": "markdown", + "id": "88c43116", + "metadata": {}, + "source": [ + "# Set up unWISE tracer" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "c8e82dfc", + "metadata": {}, + "outputs": [], + "source": [ + "dndz_unwise = np.loadtxt('../../soliket/xcorr/data/dndz.txt')\n", + "ngal_unwise = 1.\n", + "fsky_unwise = 0.4\n", + "\n", + "tracer_unwise_g = ccl.NumberCountsTracer(cosmo,\n", + " has_rsd=False,\n", + " dndz=dndz_unwise.T,\n", + " bias=(dndz_unwise[:,0], b1_unwise * np.ones(len(dndz_unwise[:,0]))), \n", + " mag_bias=(dndz_unwise[:,0], s1_unwise * np.ones(len(dndz_unwise[:,0])))\n", + " )\n", + "\n", + "Nell_unwise_g = np.ones(n_ell) / (ngal_unwise * (60 * 180 / np.pi)**2)" + ] + }, + { + "cell_type": "markdown", + "id": "8595cbf0", + "metadata": {}, + "source": [ + "# Set up SO CMB lensing tracer" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "b41fbc8d", + "metadata": {}, + "outputs": [], + "source": [ + "zstar = 1086\n", + "fsky_solensing = 0.4\n", + "\n", + "tracer_so_k = ccl.CMBLensingTracer(cosmo, z_source=zstar)\n", + "\n", + "# Approximation to SO LAT beam\n", + "fwhm_so_k = 1. * units.arcmin\n", + "sigma_so_k = (fwhm_so_k.to(units.rad).value / 2.355)\n", + "ell_beam = np.arange(3000)\n", + "beam_so_k = np.exp(-ell_beam * (ell_beam + 1) * sigma_so_k**2)" + ] + }, + { + "cell_type": "markdown", + "id": "2da86e56", + "metadata": {}, + "source": [ + "# Calculate power spectra" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "75519761", + "metadata": {}, + "outputs": [], + "source": [ + "n_maps = 2\n", + "cls = np.zeros([n_maps, n_maps, n_ell])\n", + "\n", + "cls[0, 0, :] = ccl.angular_cl(cosmo, tracer_unwise_g, tracer_unwise_g, ells) + Nell_unwise_g\n", + "\n", + "cls[0, 1, :] = ccl.angular_cl(cosmo, tracer_unwise_g, tracer_so_k, ells)\n", + "cls[1, 0, :] = cls[0, 1, :]\n", + "\n", + "cls[1, 1, :] = ccl.angular_cl(cosmo, tracer_so_k, tracer_so_k, ells)" + ] + }, + { + "cell_type": "markdown", + "id": "4113c57a", + "metadata": {}, + "source": [ + "# Set up covariance" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "620c0df5", + "metadata": {}, + "outputs": [], + "source": [ + "n_cross = (n_maps * (n_maps + 1)) // 2\n", + "covar = np.zeros([n_cross, n_ell, n_cross, n_ell])\n", + "\n", + "id_i = 0\n", + "for i1 in range(n_maps):\n", + " for i2 in range(i1, n_maps):\n", + " id_j = 0\n", + " for j1 in range(n_maps):\n", + " for j2 in range(j1, n_maps):\n", + " cl_i1j1 = cls[i1, j1, :]\n", + " cl_i1j2 = cls[i1, j2, :]\n", + " cl_i2j1 = cls[i2, j1, :]\n", + " cl_i2j2 = cls[i2, j2, :]\n", + " # Knox formula\n", + " cov = (cl_i1j1 * cl_i2j2 + cl_i1j2 * cl_i2j1) / (delta_ell * fsky_solensing * (2 * ells + 1))\n", + " covar[id_i, :, id_j, :] = np.diag(cov)\n", + " id_j += 1\n", + " id_i += 1\n", + " \n", + "covar = covar.reshape([n_cross * n_ell, n_cross * n_ell])" + ] + }, + { + "cell_type": "markdown", + "id": "9024625e", + "metadata": {}, + "source": [ + "# Construct sacc file" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "a00dfdbc", + "metadata": {}, + "outputs": [], + "source": [ + "s = sacc.Sacc()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "b91478f0", + "metadata": {}, + "outputs": [], + "source": [ + "s.add_tracer('NZ', 'gc_unwise',\n", + " quantity='galaxy_density',\n", + " spin=0,\n", + " z=dndz_unwise[:,0],\n", + " nz=dndz_unwise[:,1],\n", + " metadata={'ngal': ngal_unwise})\n", + "\n", + "s.add_tracer('Map', 'ck_so',\n", + " quantity='cmb_convergence',\n", + " spin=0,\n", + " ell=ell_beam,\n", + " beam=beam_so_k)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "58a67a06", + "metadata": {}, + "outputs": [], + "source": [ + "s.add_ell_cl('cl_00',\n", + " 'gc_unwise',\n", + " 'gc_unwise',\n", + " ells, cls[0, 0, :],\n", + " window=Well)\n", + "\n", + "s.add_ell_cl('cl_00',\n", + " 'gc_unwise',\n", + " 'ck_so',\n", + " ells, cls[0, 1, :],\n", + " window=Well)\n", + "\n", + "s.add_ell_cl('cl_00',\n", + " 'ck_so',\n", + " 'ck_so',\n", + " ells, cls[1, 1, :],\n", + " window=Well)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "f76675aa", + "metadata": {}, + "outputs": [], + "source": [ + "s.add_covariance(covar)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "71c72b97", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: VerifyWarning: Keyword name 'META_ngal' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]\n" + ] + } + ], + "source": [ + "s.save_fits('../../soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits', overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "id": "d98c9b8e", + "metadata": {}, + "source": [ + "# Read sacc file and compare to 'theory'" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "6b6831ff", + "metadata": {}, + "outputs": [], + "source": [ + "s_load = sacc.Sacc.load_fits('../../soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits')" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "f29b20ef", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ck_so cmb_convergence \n", + "gc_unwise galaxy_density \n" + ] + } + ], + "source": [ + "for n, t in s_load.tracers.items():\n", + " print(t.name, t.quantity, type(t))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "d1150f7a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Data types: ['cl_00']\n", + "Tracer combinations: [('ck_so', 'ck_so'), ('gc_unwise', 'ck_so'), ('gc_unwise', 'gc_unwise')]\n", + "Nell: 60\n" + ] + } + ], + "source": [ + "\n", + "# Type of power spectra\n", + "data_types = np.unique([d.data_type for d in s_load.data])\n", + "print(\"Data types: \", data_types)\n", + "\n", + "# Tracer combinations\n", + "print(\"Tracer combinations: \", s_load.get_tracer_combinations())\n", + "\n", + "# Data size\n", + "print(\"Nell: \", s_load.mean.size)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "165c4abd", + "metadata": {}, + "outputs": [], + "source": [ + "ell_theory = np.linspace(1,ell_max,ell_max)\n", + "\n", + "z_unwise = s_load.tracers['gc_unwise'].z\n", + "nz_unwise = s_load.tracers['gc_unwise'].nz\n", + "\n", + "tracer_gc_unwise = ccl.NumberCountsTracer(cosmo, has_rsd=False,\n", + " dndz=[z_unwise, nz_unwise],\n", + " bias=(z_unwise, b1_unwise * np.ones(len(z_unwise))), \n", + " mag_bias=(z_unwise, s1_unwise * np.ones(len(z_unwise)))\n", + " )\n", + "\n", + "tracer_ck_so = ccl.CMBLensingTracer(cosmo, z_source=zstar)\n", + "\n", + "cl_gg_theory = ccl.angular_cl(cosmo, tracer_gc_unwise, tracer_gc_unwise, ell_theory)\n", + "cl_gk_theory = ccl.angular_cl(cosmo, tracer_gc_unwise, tracer_ck_so, ell_theory)\n", + "cl_kk_theory = ccl.angular_cl(cosmo, tracer_ck_so, tracer_ck_so, ell_theory)\n", + "\n", + "Nell_unwise_g = np.ones_like(ell_theory) / (s_load.tracers['gc_unwise'].metadata['ngal'] * (60 * 180 / np.pi)**2)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "7a14780c", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAEWCAYAAABPK/eBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABjDklEQVR4nO3dd3zV5dnH8c+VDUmAAGFvWQpKhMhwW6virlrrXq1b+9jaVjueWmuXtUPbakXcPnVvVBxIVVRACQoIygh7JqyEkUDGuZ8/fifhEDJOkpPzS875vl+v88o5v3WuE+DmXL/7vq/bnHOIiIiIiIjEgwS/AxAREREREYkWJUAiIiIiIhI3lACJiIiIiEjcUAIkIiIiIiJxQwmQiIiIiIjEDSVAIiIiIiISN5QAiYiIiIhI3FACJCIiAJjZ0WY208yKzWybmX1qZkeE7L/SzL4ysxIz22RmD5pZJx9D3o+Z/cLMptbYtqyObRcGn68ys2+H7KvzdxA8ttTMdoU87o/GZxMRkchRAiQiIphZB+BN4F9AZ6A38Ftgb3D/T4A/Az8DOgLjgf7ANDNLqeOaZmaH17J9lJkltsDHmAEcVXVtM+sBJAOja2wbHDy2Zlz1/g6CznTOZYQ8bm6BzyEiIi1ICZCIiAAMBXDOPeucq3TOlTrn3nPOLQgmBr8Ffuice8c5V+6cWwV8Dy8JurSOaw4A3jOziVUbzOwYYDpwcF2BmJkzs8Ehr58ws9+HvF5lZj81swXBnprnzSwNmIOX8OQEDz0W+ABYUmPbcufchsb8DuqKVURE2h4lQCIiArAUqDSzJ83sVDPLCtl3JJAGvBJ6gnNuF/A2cFJtF3TOrQTOA542sxPMbGzwGpc65xY2M97vAROBgcBhwJXOuTLgM7wkh+DPj4FPamw7oPcnqL7fgYiIxAglQCIignNuB3A04ICHgc1mNsXMugNdgS3OuYpaTt0Y3F/XdWcAFwMv4Q0vu9Y5904EQv6nc26Dc24b8Ab7eng+Yl+ycwxeAvRxjW0f1RFrfb+DKq+ZWVHI45oIfBYREYkiJUAiIgKAc+4b59yVzrk+wEigF3AfsAXoamZJtZzWM7i/PmuACsCAVREKd1PI8xIgI/h8BnB0sPcm2zm3DJgJHBncNpK6e4Dq+x1U+Y5zrlPI4+EIfR4REYkSJUAiInIA59xi4Am8JGAWXiGAc0OPMbN04FS8OT21MrODgGnA7cD1wFQzG9HA25cA7UNe92hE6LPwijRcC3wK1T07G4LbNgSH5jWoxu9ARERihBIgERHBzIab2U/MrE/wdV/gImC2c64YrwjCv8xsopklm9kA4EVgHfB/dVyzF15y9Afn3BPOuZeBn+IVRhhUTzjzgIvNLDFYQOG4cD+Hc64UyANuxRv6VuWT4LY6e3/q+x2E+/4iItL6KQESERGAncA44DMz2433pX8h8BMA59w9wC+BvwI78IoNrAVOdM7trfWKsBX4iXPuwaoNzrmngRuBwnpiuQU4EygCLgFea+Rn+Qjohpf0VPk4uK3OBIgGfgdBb9RYB+jVRsYmIiI+M+ec3zGIiIiIiIhEhXqAREREREQkbigBEhERERGRuKEESERERERE4oYSIGlVzCzBzO4ws7VmtsHMzjSzMjPLqm+f33GLiD8aaDO+a2azg8d1M7NZZnar3zGLiH/qaxfUZsSP2ha1E/HTncDxwHhgN/A2UOCc225md9W1z5dIRaQ1uJO624zDgPlmNhx4Cfhf59xrPsUpIq1Dfe2C2ow4oQRIWg0zywZ+DBzunFsf3PYecER9+/yKV0T8FUa7MApvAdeZwBnOuZm+BCoirUl97YLajDihIXDSmpwI5Dvn8kO2dQa+amCfiMSnhtqFUcAQvPV9hkc5NhFpneprF9RmxAklQNKadAU2VL0ws0TgVGBBA/tEJD7V2S6YWUegH16S9CfgV2amUQ8icay+dkFtRnxRAiStyTfAkWY22Mw6AP8EDsK7m1vfPhGJT/W1C6OAlc65bcBjQApwhW+RikhrUF+7oDYjjigBklbDOTcdeBb4EsgD5gElwOL69vkRq4j4r4F2YRTBHmLnXBlwN7qjKxLv6msX1GbEEXPO+R2DSK3M7HrgdOfcmY3ZJyLxSe2CiIiEQ1mttBpmNh7YCKzFG4N7F3BWQ/tEJD6pXRARkaZQAiStyeHAW0AysBS40jk3O4x9IhKf1C6IiEijaQiciIiIiIjEDRVBEBERERGRuKEESERERERE4kbU5wCZ2UTgH0Ai8Ihz7u4a+y24/zS8cqZXOue+MLO+wFNADyAATHbO/SN4zl+AM4EyYDlwlXOuqKFYunbt6gYMGBChTyYiVebOnbvFOZftdxyRpjZDpGWozRCRxmhumxHVBCi4SvcDwEnAOmCOmU1xzn0dctipwJDgYxzwYPBnBfCTYDKUCcw1s2nBc6cBv3DOVZjZn4FfALc3FM+AAQPIy8uL4CcUEQAzW+13DC1BbYZIy1CbISKN0dw2I9pD4MYC+c65FcFFpp4Dzq5xzNnAU84zG+hkZj2dcxudc18AOOd24q0A3jv4+j3nXEXw/NlAn2h8GBERERERaVuinQD1xluvocq64LZGHWNmA/DKn35Wy3t8H3i7uYGKiIiIiEjsiXYCZLVsq1mHu95jzCwDeBn4kXNux34nmv0Kb6jc03UGYHatmeWZWd7mzZvDDlxERERERNq+aCdA64C+Ia/7ABvCPcbMkvGSn6edc6+EnmRmVwBnAJe4ehY3cs5Nds7lOudys7Njbr6liIiIiIjUI9oJ0BxgiJkNNLMU4EJgSo1jpgCXm2c8UOyc2xisDvco8I1z7u+hJwQry90OnOWcK2n5jyEiIiIiIm1RVKvABau03Qy8i1cG+zHn3CIzuz64fxIwFa8Edj5eGeyrgqcfBVwGfGVm84LbfumcmwrcD6QC07w8idnOueuj86lERERERKStiPo6QMGEZWqNbZNCnjvgplrO+4Ta5wfhnBsc4TBFRERERCQGRXsIXEy44KFZXPDQLL/DEBFpGY+f7j1ERMJw1TtXcdU7V0XsOJGWpgRIRCQetERSo0RJRETaICVAIiLS8pQsicS9XWW72Lh7I/MK5zV4rHqLpCVFfQ6QiIiIiMSOXWW7KCgp4I5P76DSVbJ1z1ZKyksYmjWU/x3/vwBc+talLN6+GIAr3rmC4/scz/F9j+ecIecAsGbHGrq260r75PbV19xZvpN5hfPI6Zbjy+eS2KUESERE9rd3B+wphrWfQ9+xfkcjIj6o6n15fOLj+22vDFSyYMsCZqybwZbSLZw35DyWbF+Cw/Fq/qtkpWbRK6MXGckZpCenV58XILDvuQswa+MsUpNSOWfIOTjnuPDNC9lZvpNu7bqRlZbFku1LALjmvWt4+OSH602C6opVpC5KgERE4kHNpKa0CHZuhIzu0L4z7NgAS9+Fbcth01eAg8dPg3Mmw6Hnwo6NsOoTSEqBpDRI6wS7C6F8T+QTpaqhcle9FblrikizrNu5jpeXvcxr+a+xpXQLiZbI2B5jmbNpDg5v/flES+TyEZdz9aFXH3D+bUfcxuVvX47DkZaYxuSTJjMqexQADscdR97B2h1rWb1jNZ9v+rz6vPJAOTM3zOQ3M3/DyK4jGd1tNId3P5yBHQYSXPpEpNGUAImItFXhJgprP4eCr8A5ePRkSG4H5cE1o78zCXIugqI18OaP9j8vUA5L3/YSoE0L4JUDv9QA8ORZcPztMPcJ6NAHOvaGjn28x8FnQXpX9SqJ78zsMeAMoNA5N7KW/T8DLgm+TAIOBrKdc9vMbBWwE6gEKpxzudGJuvWYvmY6jy18jGN7H8vpg07nyN5H0iGlA/MK52EYDkdyQjK53Wv/1eR0y2FY1jB2lu/k7mPu3q9HJ8ESmDhgYvXreYXzqpOl5IRkDu58MN9s+4aP133MlOVTAMhKzeLOI+/kW/2+xY6yHewq26XhchI2JUAiIrHGOVgzC758GvoeASVbvW3eTug+wktMOvaGPkd4m3uNhh9/DQUL4dkLwQW8np7c4CTkAUfDzXlQsRcq9kDeYzDvaW9fZZmXQPUeA8XrYfVMr0fJVUL/o2HbipBepVPhhF9BziWQ2T3avxmJb0/gLZz+VG07nXN/Af4CYGZnAj92zm0LOeQE59yWlg6ytSjaW8SGXRt4cN6D3JBzA+cMOYdTBpxCj/Qe+x1XX2JTU0ZKBhkpGQ0mKbVd84R+J+CcY9WOVXxZ+CVfFHxB74zezCucx7Lty3A4Ln/7cs4efDbnDjmXQ7seSlKCvuZK7fQ3Q0QkVgQqYeErMOMvsGUJpGRC18Ew4BiwhGBS0w5O+eOBvTBJKcGem97QfaTXW3PeI/uOS0mHrkP2He8CMP9Z72diCoy6aP9rBiq9JKhDL/j0PggOkSFQAdN/6z2yBsD1n0BqppegaTiLtCDn3AwzGxDm4RcBz7ZgOK2Wc477v7yf/KJ8ACYvmMyEXhPI6ZZDh5QOtZ4TbmLTmDk6tV3TzBjYcSADOw7k3CHnAvDIV49UD8FzOF7Pf53X8l8jMzmT177zGt3ad6OssoyUxBRA84XEowRIRCRWvHglfDMFuo2As/8NI77jJS5Qe1JTl9QO3qO+4/qOrf+aCYnQqa/3fL8ELA0m/hnKdnq9RamZwdivgN1b4ZCzvB6r8hINlxNfmFl7YCJwc8hmB7xnZg54yDk3uZ7zrwWuBejXr19Lhhpxu8t38/vZv+fNFW9Wb3M48gry6k1uWiKZCPeaud1zq4fgpSWmcd8J97G7fDcLNi8gu102AHfOvJOvtnzF0b2PpqCkgMpApYbLxTklQCIibZlzUFkOickw/kYv6TnkHEioscxbOElNY4V7zYaSJYDuh8KiV+Dt2/Zte+J0uPKt+q+vggkSeWcCn9YY/naUc26DmXUDppnZYufcjNpODiZHkwFyc3Ndbce0Vh+u/ZCpK6dy3pDzeGXZKw3O62kN6hqCd/KAk6uPGdtzLNv2buO5xc9R4SoAuPKdK3li4hNKguKUFkIVEWmrSrfDhrnw0g+81/0nwMjzDkx+WoPUDtCxb93JzHE/gxtnwbgb9m2rLIdVH0NlBRQujk6cInAhNYa/Oec2BH8WAq8CMdU1uadiDwCnDzqdN7/zJnceeSfDsobRO6N3gyWoW4MXz3qRd857p844vzP4O0z69iSuOeya6m0BFyCvII/KQCW3fXQbLy19iS2l+6Z4aSHW2NYK/5ds/XaUlrN+ewlzV2/3OxQRiVfL3ofCRV5BgiVTveFikXLVW+H1qIR7XGOMPNcbLgeQlOoNn1syFf49zuvtWTwVAoH6ryHSRGbWETgOeD1kW7qZZVY9B04GFvoTYeQt2rqI0185nYVbvI/Ut4M3dDUjJYOe6T1bffLTGEf2OhLDm2uYmphKbvdcNpVsYsGWBfx21m/51gvf4rKpl/Gfr/9D0d4iNu7eyLzCef4GLS1CQ+Aa6d5pS/hm004Azp80kwcvHcMpI3o0cJaISASVFsHrN+177QJeT0kszJepbbjc7q1w0l3w+SPw3EWQfTAc/WM49Lt+RyttiJk9CxwPdDWzdcBvgGQA59yk4GHnAO8553aHnNodeDW45kwS8Ixz7p1oxd2S8rfnc/2060lPTqdru6777YvFIgF1DZd7+9y3Wbp9KR+s/YD3V7/Pn+f8uXpe0dXvXc3fj/s7x/Y91t/gJaLMuTY1PDWicnNzXV5eXtjHf7xsM5c9uv9d1o7tkpjzq5NISVJnmkgVM5sbi+tkNLbNaBHOwTMXQP77XuKD8yq7XTElNhIgqHteT2UFLHwZPrnXK7198xyYfHz4xR2k1VKbEX0bd23kkqneskdPTnyyuudH4J7P7+H/vvk/gOoeo3E9xzFxwES+3f/bdEzt6Gd4QvPbDH1rb4QJg7pw7bGDqiu1piYl8LOThyv5EZHoMYNx18HZ90OPQ6FT/9hKfuqTmASjLoAbZsIVb8D6ud66RUWrvYIJaz7zO0KRVq1qXktpRSm3fHALpRWlTD5pspKfGk4ecHJ14pOckMyZg85k/a713DnrTo5//nhun3E7oR0Imi/U9mgIXCMkJSbwy9MO5pNlm9mxp4J/XHg4Y/pnAbBuewl9str7HKGIxLSy3V5Z68Eneq+/fDryld1ag4bmFSUkeOsVLXgu2AuG1yP00pVw3qPQ/8gWD1GkLUuwBIZ3Hs7Nh9/M4KzBfofT6tQ2VM45x9fbvubtFW/jcASHRPLQ/IfYXLKZ8kC5Smu3IUqAwlBWEeD8h2Zx/bGDOPXQnmSmJZOZllyd/Dw/Zw2/fm0R7996HP26KAkSkRZQsg0mHQ3H3ApHXO13NK1D6PpCCclQUQaPn+pVwjvrfkgJaY9VLlsE8BY6TU1M5a6j7vI7lFbtxbNe3O+1mTGiywhGdBlRva1gdwEPLXiI8kA54PUE3XPcPZzU/6SoxiqNp7FbYfjv4kLmry0iLSWx1v3HD+uGGdw3fWmUIxOR2pjZY2ZWaGa1Vmoyzz/NLN/MFpjZ6GjH2Ghv3w67CqBPjPX2NEdVwYRO/eGqqfCjr+C422HvLkhu5x0Tx/NcRWraVrqNBVsW8NYK3QiIhO7p3fnByB9Uv65wFdz64a3MWFfrElHSiigBCsP73xTQsV0yxwzuWuv+7h3SuOLIAbz25XryC3dGOToRqcUTeCu51+VUYEjwcS3wYBRiaro1s+GrF7zKZz0P27e9JcpQtzWh6wultIcTfgkXP+/NlSpeBw8e5RWMEIlzcwvmsmLHCsoD5fxm5m9U3jlCjup91H6ltb879LuM6T4GgFeWvcLvZ/+exdu0jllroyFwDQgEHB8uKeS4odkkJdadL15/3EH8Z/Zq/v3Bcv5+QU70AhSRAzjnZpjZgHoOORt4ynmzWGebWScz6+mc2xidCBshEIB3fwmZPb0ESBpWVammZKs3N+g/50FqJ0jN8NZLirU5UyJhmLxgcvXzikAFeQV5mq8SATndcnjq1KfIK8gjt3vufr/T9bvW8+qyV3l+yfOM6DKCc4ecy/OLn2d3xe79ynBL9KkHqAEL1hezZVcZJx7crd7jOqen8L3cvkxfXEhJWUWUohORJuoNrA15vS64rfUp+Ao2LoAT7/AKIMj+6usF6zkKbvgUci6DvUWwYx08flpkF40VaQM27d5EXsG+ctzJCcnkdo+5quO+yemWw9WHXn1AQvPDw3/If7/3X34+9ueUBcr43ezfsbRoKet3reea965RL5yPop4AmdlEM1sSHHv/81r21zo238z6mtkHZvaNmS0ys1tCzulsZtPMbFnwZ1ak4k0wOOfw3hw7JLvBY3/4rcHM+NkJtE9Rx5pIK2e1bKt1soiZXWtmeWaWt3nz5hYOqxY9R8EP58JhF0b/vWNBUip0GbjvdaDCWzRWJI48880zGMagjoPondGbh09+WL0PUdIxtSOXHHwJL5/5MhcMu6B6e1mgjFs/vJVnFz/LrrJdPkYYn6KaAJlZIvAA3vj7Q4CLzOyQGofVNTa/AviJc+5gYDxwU8i5PwemO+eGANODryPisD6duPeCHLLSUxo8tktGKh3bJwMQzwvMirQB64DQhS/6ABtqO9A5N9k5l+ucy83ObvhGSETtDc4pzOrvlX6WpqmqFgdeQjTgGNgwD2Y9AIFKX0MTiYb/Gf0/PH7K43RO60zP9J5KfnxgZpwx6Izq+UJJlkRGcgZ//OyPfOvFb/G7Wb9jybYlPkcZP6L9P+pYIN85t8I5VwY8hzcWP1T12Hzn3Gygemy+c+4LAOfcTuAb9g1ZORt4Mvj8SeA7kQjWOcfabSWNSmY2Fe/hjH99zFtfeVMJLnhoFhc8NCsS4YhI5EwBLg/2OI8Hilvd/B/nvOFab9zS8LFSv9BqcVe84b3+6kVvbtUTp8PW5fsf//jp+8pmi7RhARegpLyEpIQkDs0+lMcnPs7jEx/3O6y4VTVf6JbRt/DoKY8y5ZwpPHv6s5wy4BReX/46579xPoUlhdXHnz/lfCa+PFFD5VpAtBOgcMbdN3hMcHLz4UDVst/dq768BH/WOWGnMcNZ1m4r5Zh7PuDFvHX1HhcqOzOV7bvLeX7O2oYPFpEWYWbPArOAYWa2zsx+YGbXm9n1wUOmAiuAfOBh4EafQj1Q1Zfv5dNh0wLorXH6ERFaLQ7g5N/DdyZBwddepbjPHvIKTojEkOlrpnPqK6eyomiF36FIUM35QiO7juR3R/2O6edP52/H/41u7b2vsDe+fyOLty/WfKEWEu3JKuGMu6/3GDPLAF4GfuSc29HYAJxzk4HJALm5ufV27Xy5djsAI3p3CPv6iQnG+bl9uO/9ZazdVtLY8EQkApxzFzWw3wE3RSmcpvn8EUjvBodd0PCx0nhmkHMRDDoOpvwPvH2bt3bQ6Mv9jkwkIioDlfx73r/pmNqR/h36+x2ONKBjasfqBVQrA5Vs2LVvVPaeyj28uPRFDss+jATTcOhIiPZvMZxx93UeY2bJeMnP0865V0KOKTCznsFjegKFRMBX64pJTUpgWPfM/bY/f90Enr9uQp3nfS+3L2bwQp56gUSkCSr2wrJ34fBLIanh+YcShrqqxXXoBZe8COc/AaOCefPuLVC0VtXipE17Z9U75Bflc2POjSQm1L6Qu7ROiQmJ3HnkndXzhQCmLJ/C09887WNUsSXaCdAcYIiZDTSzFOBCvLH4oWodm29mBjwKfOOc+3st51wRfH4F8Hokgl1auIvB3TLqXf+nNr06teOYIdm88sV6FUMQkcbbVeDNARpzRcPHSvOZwYhzIDEZln8AWxZD8Wp48kwlQdImBVyAhxY8xJCsIZzc/2S/w5EmyOmWw7CsYfTO6M3jpzzOPcfew2kDTwPgv2v+yx8/+yMri1f6HGXbFdUhcM65CjO7GXgXSAQec84tqhqX75ybhDc2/zS8sfklwFXB048CLgO+MrN5wW2/dM5NBe4GXjCzHwBrgPMjEe+ygp2MH9SlSeded+wg1m8v5cW5a2sd0yciUqu9OwAHJ/4Gsgb4HU38Wf/FvucVe+DT++CCp/ctrirSBswtmMvK4pXcfczdGjLVhr141ou1bl9RvIKXlr7Es4uf5ajeR3HpwZdy39z72FW+SwushsniuYciNzfX5eXl1brPOccbCzbSPTOVcU1MgoDqCnD1DZkTiTVmNtc5F3Oz9+trMyJi7efw2CngApDUDq6Ysm/SvkRH6J+BJXg/Bx0PF/wHUjMbPF2aRm1GZDnn+LLwSw7NPpTkhOSov7+0vC2lW3hx6Yu8sOQFtpRuqd6elpgWF+s8NbfN0G2BOpgZZ43q1azkZ+uuvWzasYdAHCeZItIIqz72vnADVJZpwU4/hJbMvuodOONeSOsEKRm1H6+S2W2GmT1mZoVmtrCO/cebWbGZzQs+7gjZV+8i7q2NmTG6+2glPzGsa7uu3DDqBt477z1OGXBK9fbyQDn3zb2PNTvW+Bhd66cEqA7LN+9i7urtVAaanrzMX1fE6q0lrNy8m7mrt0cwOhGJSX1CensSU7wFOyX6qkpm9xsHud+H7z3pDYErWgPPXeIVSJC26AlgYgPHfOycywk+7oKwF3FvNe6adRf3zr3X7zAkSpITk7n04EurCyYkWALzNs/jjFfP4KbpNzFzw0zNR6+FEqA6PPvZGi55ZHaz5u9kpHpTrLbsLuOSR2YrCRKR+u0p9n6md9Pwt9ao4GuvSMK/J8DcJ7xCFdJmOOdmANuacGo4i7i3CltLt/Jq/qvsrdzrdygSRaELrD52ymNM++40rht1HQu3LOS6adfxnde/Q/HeYr/DbFWivQ5Qm7F6Wwn9O6eTkND0FGjOqn0JT3lFgNkrtjKmf1YkwhORWLTwZUhIgi5DlPz4qbZy2QDDJsKNs2DKzfDGLfD1617J7IpSb+6Q/sxiwQQzm4+3/MZPnXOLqH2B9nF+BNeQN1e8SUWggvOHRqQWlLQhOd1y9pv3c1POTVxz6DW8u+pd8gry6JjaEYBTXjqFSlfJX4/7a8zPE6qPeoDqsHZbCX07t2/WNcYP6lLdg5SQYE2uKCcicaC8FJa+A+27quJYa5bVHy57HU77K6ya6ZXMLloNT56lktlt3xdAf+fcKOBfwGvB7eEs4l7NzK41szwzy9u8eXPko6yDc46Xl71MTnYOB3U6KGrvK61XSmIKZx50Jr898rcAfLbxMzbs3kBBSQFXvHMFz3zzTNwOj1MCVAvnHGu2ldCvmQnQmP5ZDO+RSaIZpx/aS70/IlK3nZug5ygvAZLWLSEBxl4D467dt61yr9cjJG2Wc26Hc25X8PlUINnMuhLeIu6h15nsnMt1zuVmZ2e3aMyh5m2ex8rilZw75Nyovae0LV9t+ar6ecAF+NPnf+LcKefyRcEX9ZwVm5QA1WLr7jJKyirp17lds6/VoV0yh/frxH0X5jQ/MBGJXZ0HwvffgXad/I5EwnXwmV6p7CqfT4aP/wYVZf7FJE1mZj2Ci65jZmPxviNtJbxF3H3XJa0LFw+/eL+KYCKhcrvnVhdLSE1M5dpDryU5IZmsNO8G/eodq1m/a72fIUaN5gDVIjMtiReum0CfrOYnQACJwXlEzjlMQ1tEpCbnYO9OSOtQ9/wTaX2qSmbvKYaJf4L5z8H0u2DBi3DmfdBv/P7HV5XL1p+xL8zsWeB4oKuZrQN+AyRD9ULs3wVuMLMKoBS40Hnjg2pdxN2Hj1Crq97x1ot/fOLj/GLcL3yORlqzqmIJeQV55HbPJadbDj8c/cPq/f/44h9MXzOdE/qewCUHX8KfP/9zzC6uqgSoFqlJiYwd2Dmi1/zFKwvYtruMhy6LuXXeRKS5ChbB5OPgwmdgqO7etimpHbzH8NO9x5K34a2feoupXvIyDPm23xFKkHPuogb23w/cX8e+qcDUlogrEnaV7yJvUx5juo/RjVapV81iCaFuO+I2+nfoz0tLX2L6munV269575qYW1xVQ+BqMXf1Nt6YvyGiE8MyUpP47+JCikvLI3ZNEYkRy96FQIU3B0jatmGnwk2fwYm/gUHHedu2r4JAwNewJLZt3LWRX37yS1zdtRlEGtQjvQe3jL6Fad+dxon9TqzeXh4o5/ONn7Np9yYfo4ssJUC1eH7OWu568+uI3kU5/bBelFc6pn1dELFrikiMWPoe9MyBzB5+RyKNddVbBw5pS82AY26FxGTYuwsePQUePckrdFG8VtXiJKKK9xZTXFbM6G6jSTB9rZPmS0tK48oRV1bPF0pOSCYxIZGJL0/k1g9vZW7B3DZfPU7/UmqxsXgPvTpFZv5PlVF9OtK7UzveWlBn4RgRiUd7d8L6PBh8YsPHStuT3B6+/RvYuhy25QdLZp+pJEgiYl7hPPKL8gGYtnoa8wrn+RuQxIzQxVUfPvlhTh90OpePuJzPNn7Gle9cyQVvXsBr+a9REajwO9QmUQJUi03Fe+jRITWi1zQzTj+sJ5/kb6G4RMPgRCRozWxv+NvAY/2ORFpCQgLkXLx/yeyKPfB1qysiJm1QXkFe9bC3SldJXkGezxFJLMnplsPVh15NTrcceqT34NYxt/L++e9zx4Q7KA+UM2n+pOpeovOmnMfElye2mSRcRRBqsWnHHo48KDKLlj5/3YTq5+eO7k3HdskRua6IxIjs4XDy76Fvq1xYXiJl8Ldhxl/ABSAhySuhDbBxAXQfAQmJ+x+vinEShkO7Hlr9PDkhmdzuKrQkLatdUjvOH3o+3x3yXQpLCklMSGTOpjks3b4UgB+8+wMePeXRVl8wQT1ANezeW8HOPRV075gW8WsP79GBm04YTMf2SoJEJKhTXzjyh5Ac2WG30spUlczu1B+uehv6jfPmBD16Ejx0LCyb5pVDF2mEcT3HMTRrKD3Te8ZclS5p3cyM7undAfhs42fV28sCZdz64a28vPRlSspL/AqvQUqAamiXnMjHt53ABbl9Gz64CUrLKpn61UYNgxMRb/2YRa96PyX2pXaAjn29ZAggvRt8599Qtgue/m5wbtAcf2OUNqdDSgd6Z/RW8iO+Obr30dVD4ZISkkhJTOHOWXeyeNtigFY5T0hD4GpISDD6dm7fYtdfVriTG5/+gnu+exjfa6EkS0TaiJUz4MUrvR6B/kf6HY20tJrD2RISYOR5MPxMmPsEfPRneOxkuGUB7N3hJcZrP9+XMImEKK0o5dKpl1IZqCQrLcvvcCSO1VxgdVT2KBZuWcjIriMBuPvzu8kvyueCYRfw7X7fZtHWRfstxuoHJUA1zF9bxKwVW7lsfH/SUyP/6zm0d0f6dm7HWws2KgESiXerPoGkdtBb4/bjWlKKVyQh52JY+RHs3AgFC735Qk+cAVe+qSRIDjB7w2yWbl/KQyc9xJG9dANF/FVzgdVDs/fNTxvUcRCfrP+E22bcRmZKJjvLdgKQlpjm29BNDYGr4dPlW7j77cW01ELKZsbph/bi0/wtbN9d1jJvIiJtw5rZ0CfX+wIskpoBw0+HVR97yQ9A5V545VpY+XHtc4QeP31fwQSJKx+t+4j05HSO6H6E36GI1Ovigy9m6rlTefDbD9K1Xdfq7eWBcuZsmkNloDLqMSkBCnHBQ7N48tNVZKYl0T6l5TrHzjisJxUBx7uLYmdFXRFppLLdsOkrVX+TAw04BqoWtExIgtLt8OQZ8PAJUPC1v7FJqxBwAT5a9xFH9TqK5EQVVpLWL8ESOLr30dx15F37LbDaKbUTp71yGpMXTKZgd0H04onaO7URZZUBenSIfAW4UCN6dWBAl/Z8vmpbi76PSLwys4lmtsTM8s3s57Xs72hmb5jZfDNbZGZXRT3IjfPBVWpokxyoZsW4nyyG0/8OleWQ4VVdYtsKL4mWuLRoyyK2lG7h+L7H+x2KSKPUXGB1WOdh9M3sy7++/Bcnv3wyN0+/mf+u+W+LF06I+hwgM5sI/ANIBB5xzt1dY78F958GlABXOue+CO57DDgDKHTOjQw5JweYBKQBFcCNzrkmLbNdXhmgW4QXQa3JzHjphiPpkq5hLyKRZmaJwAPAScA6YI6ZTXHOhd46vwn42jl3ppllA0vM7GnnXPTGpfabAP/z5b4vtCKhUjt4j6oE+YgfQO73wcwbCvfyNbBtOViSN4RSxRLiSnJiMqcPOp2jex/tdygijVZzvtAjpzzCmh1reDX/VV7Pf515m+cx/fzpXDTlIorLirnn2HsiPk8oqj1AIV9MTgUOAS4ys0NqHHYqMCT4uBZ4MGTfE8DEWi59D/Bb51wOcEfwdZOUVzq6ZrRsAgTQNSMVa6mJRiLxbSyQ75xbEUxongPOrnGMAzKDN1wygG14N0+ixww6D4KU9Ki+rbQRV711YNW4qv8zzOCUP0D2wVCyGXash8dOgblPRj9O8cXwzsO5+5i7Vf1NYka/Dv24ZfQtvPfd93jilCf4Zus3LNm+hI27N3L525fzwJcPUFpRGrH3i/YQuHC+mJwNPOU8s4FOZtYTwDk3A++LSk0O6BB83hHY0NQAD+vTkd99Z2TDB0bAP6cv4+on86LyXiJxpDewNuT1uuC2UPcDB+O1FV8BtzhXNes8CpyDN37kTWwXaYp+42HIt/e9dgFY+o73vLwU9u468BwVTIgJpRWlrCpehdPCuRKDkhKSGJw1mLyCPBze33GHY9KCSZz4won8fvbvWVG8otnvE+0EKJwvJuEcU9OPgL+Y2Vrgr8Avmhpgghkd0qIzodCA978pYGNx5DJaEaG2rtWa3xROAeYBvYAc4H4z60AtzOxaM8szs7zNmzdHJsKty2Hu4948DpGmCi2WkJQG42/0ns97Gv5+MLx5qzfXTGLKZxs/48zXzmRuwVy/QxFpMbndc0lLTCPREklNTOWX437JcX2P47X811hZtLLZ1492AhTOF5NwjqnpBuDHzrm+wI+BR+sMoJ4vM+WVAVZt3c2iDdFZlf2MUb0AeGvBxqi8n0icWAeELrLVhwN7ha8CXgn2NOcDK4HhtV3MOTfZOZfrnMvNzs6OUITBKYqasyHNEVos4Yo3YOAx3vY+R8CwU+HL/8BDx3qPOY96C6sWr/XmC8UhM3vMzArNbGEd+y8xswXBx0wzGxWyb5WZfWVm88zM16EbszfOJi0xbb91VkRiTU63HB4++WFuPvxmHjn5ES4afhF/OuZPTD9/Osf2PbbZ1492AhTOF5NwjqnpCuCV4PMX8Yba1aq+LzN7KwIU7NjLxqI9DbxdZAzsms6hvTvyxvwmj9gTkQPNAYaY2UAzSwEuBKbUOGYNcCKAmXUHhgHR647Z8CWkZEDXYVF7S4lRqR2gY9/9k+meo+DcyfDTJXDqPVBZAXMegcJFULTaK6m95jP/YvbPE9Q+j7jKSuA459xhwO+AyTX2n+Ccy3HO+bpy8ewNsxndfTSpiS0/X1nETzndcrj60Kv3K4DQMbUjyQnNH6kV7QQonC8mU4DLzTMeKHbONdRFsgE4Lvj8W8CypgRXUelNAeiSEb3qbGeO6sn8dcWs2VoStfcUiWXOuQrgZuBd4BvgBefcIjO73syuDx72O+BIM/sKmA7c7pzbErUgN8yDHodBglYikGaqrVhClXZZMO46uOFTGH7GvsVVK/bCcxfDzPthVx3DOmNwvlA984ir9s90zm0PvpyNdwO2VSnYXcDy4uVM6DnB71BE2rSolsF2zlWYWdUXk0TgsaovJsH9k4CpeCWw8/HKYFevz2FmzwLHA13NbB3wG+fco8A1wD/MLAnYg1c9rtHKK72RdtGoAlfl9MN6saFoDwkJ3kKsAM9fp4ZNpDmcc1Px2pLQbZNCnm8ATo52XME3h8q9Gv4m0WMGQ06Cj//qJUEJyZCeDe/9Ct7/DQw5GU68A7od7HekrckPgLdDXjvgPTNzwEPOuZq9Q1Hx2Sav5258r/F+vL1IzIj6OkBhfDFxeGt01HbuRXVs/wQY09zYyn3oAerdqR13njUiau8nIj4zg+tmeImQSLRUzRfaUwznPeK9LlwM85+BBS9AYvD/vU0LvQV69xTD3h1xub6QmZ2AlwCFLrJzlHNug5l1A6aZ2eJgj1Jt519L8EZsv379IhrbcX2O4+/H/52hWUMjel2ReKPxFyEqnSPBoH1KdPPCQMAxZ9U29lZURvV9RcRHWgdMoq3mfKFuw+Gku+DHX0OXg7xtM+7xiiYUfBWcL3RmXBVNMLPDgEeAs51zW6u2B3uNcc4VAq/SxLnGzdUxtSMn9T+JBNPXN5Hm0L+gEH2z2pPbP/qLim3ZtZfvPTSLzTv3Rv29RSTKPvgTvHCF31FIPKprvlDoXLQz7oOhp+57XbEH3vll/deNkflCZtYPr6DSZc65pSHb080ss+o53vDZWivJtaStpVt5ctGTFOwuiPZbi8QcJUA1mA93Zbt1SGP8wC5s3V2mhc1EYt2KD2CXvsBIK9W+Mxxz6771hRKSvLLaABVl8N6v2+zaQsF5xLOAYWa2zsx+UKM4yh1AF+DfNcpddwc+MbP5wOfAW865d6Id/5xNc/hr3l/ZXBqh9chE4ljU5wC1Zqu3lpCW7E9OeHZOL2at2MqqrSXMXb2dMT70RIlICwtUwqavYLR6gKQVq22+EHh/dz+bBDP/6e3PuRgO/Z43V2hPcaufL1TXPOKQ/VcDV9eyfQUw6sAzouuLwi9ol9SO4Z1rXbJMRBpBPUAhtu0uY9feCl/eu3dWOwAKd+7lkkdmM3f19gbOEJE2Z8tSKC/x1mkRac1qW1+ozxj4yRI47a9e0YR3fwl/G+YlRkWr4cmz4mq+ULR9UfAFo7JHkZSge9cizaUEKERFIEBSgj8TkxesK65+Xl4RYPaKrfUcLSJtUsEi72cPreAurVxd84Xad4ax18C1H8CNs2HgsXgVovHmC037NRR+U/s1Y2SukB92lO1g6faljO4+2u9QRGKCEqCgsooAAQclZZW+9L6MH9SluihUclIC4wd1iXoMItLCUjNh4HHQdYjfkYg0X7eD4YRf7psvZAmwdg78ezw8fCLkPe4NjZNmy9+eT6IlMqZbs1f8EBGUAFX7eJk3qXDHngpfhqCN6Z/FwT0y6ZPVjqevHq85QCKxaOgpcMUUSIreYssiLapqvlCn/vD9d7whcif/Acp2wZs/gml3eMc55yVDxWs1TK4JRncfzacXfcrh3Q/3OxSRmKAEKOizFduqn/s1BC0zLZn0lERue2k+23eXRf39RaSFVZb7HYFI5IXOF8rIhiNv9obHXf1fmHCzd8yX/xe3awtFSvvk9iQnJPsdhkhMUAIUdMrIHiS0giFoyYkJLN+8mze/2ujL+4tIC9mzA/7YC/Ie8zsSkciqbb6QmVc0oWq454Yv9+2r2AOv3gCLXtNNgTBUBCq4btp1fLj2Q79DEYkZSoCCxvTPYngrGILWPiWR4T0yeeWLdb68v4i0kMJvoLIMMnv6HYlI9I26aP+1hfYUwZT/gYA/lVfbkuVFy5m5YSa7y3f7HYpIzFAtxaAPlhSyoWgPg7LTfZ1/Y2acN7oPf5j6Dcs37+Kg7AzfYhGRCHn8dNgZ7NXtdoi/sYj4oebaQr3HwJZlkNzO78havQVbFgBwaFdVjxSJFPUABS0v3EVRaTn+FMH2PH/dBJ6/bgLfObw3SQnG83PW+hiNSPSY2bFmNsjM/mNmL5jZsX7HFHHlJZCSCZ36+R2JiD9C5wolJEK3+hf0jIt2IQwLtyykY2pH+mb29TsUkZihHqCg4lJvHHKiT+sAhcrOTOWnpwxjRK8OfociEi0XAanArUAR8CQww8+AIq5st1c22PxvY0R8Udu6QvWL/XYhDAs2L2Bk15GY2g6RiGkwAQrecVkH3AWkAPc752KuASoqKScpwVpNA3P9cQf5HYJINI0AdjrnCgHMLPYWD0nPhjFX+B2FSFsS++1CAwIuQM/0nozrOc7vUERiSjg9QHFxB6ao1EuAWpO120r4NH8LF47VkBmJeb+mejl5AN7xK5AWsXeHN9m761C/IxFpS2K7XQhDgiXw72//2+8wRGJOOAlQXNyByWqfTLuURL/D2M+rX67n79OWMuGgLvTvku53OCItxjn3EYCZDQfOBnqb2QnABmCKc+4bP+NrlrWfQ8FCcAF48ixvIdS+Y/2OSqTVi+l2IUwBFyDBNF1bJNLC+Vf1a+DPIa/fbaFYfHXX2SMZ2j3T7zD2c35uHxIMFUOQuGBmtwHPAQZ8DswJPn/WzH7uZ2zNsupjL/kBrwz2qo/9jUekDYnZdiFMt824jRvfv9HvMERiToM9QFV3YEJev9py4Uionh3bccKwbrw4dx0/PmkoyYm6CyQx7WpghHNuv5URzezvwCLgbl+iaq4Bx+x7npiy/2sRaUhstgthmr95PjnZOX6HIRJzwv5GbWbDzex2M/unmf0j+Pzglgwumi5+eDabivf4HcYBLhnfj8079/L2wk1+hyLS0gJAr1q29wzua5v6joXk9pCYquFvIo0Xm+1CGLbv2c6m3Zs4pIvWDhOJtLDKYJvZ7XjFEJ7D64IG6IPXBf2cc65N34FxzvHZym10z0z1O5QDHD+0G0O6ZbB6i1aAlpj3I2C6mS0DqsZ99gMGAzf7FVREVJZD+85KfkQa70dEqF0ws8eAM4BC59zIWvYb8A/gNKAEuNI590Vw38TgvkTgkWh871m8bTEAwzvXv16SiDReuOsA/YAIdUE31Ig00ADV2XiZ2Q/xGsMK4C3n3G3hxlRaXkllwLWKNYBqSkgwpt5yjIa/Scxzzr1jZkOBsUBvvHH+64A5zrlKX4NrjpJtECjXivciTRDhduEJ4H7gqTr2nwoMCT7GAQ8C48wsEXgAOKnqvc1sinPu60a+f6NUJUAHd46ZwTYirUa4CVBVF/TqGtsb1QUdZiNSawMU3PcEtTRewaowZwOHOef2mlm3cGMC2LmnAmgdi6DWpir5Kdixh+4d0nyORqTlOOcCwGy/44ioxBToMhRSM/yORKRNilS74JybYWYD6jnkbOAp55wDZptZJzPrCQwA8p1zKwDM7LngsS2WAF31zlUU7y3m0oMvpVNap5Z6G5G4FW63wo/wuqDfNrPJwcc7wPTgvnCNJdiIOOfK8IbUnV3jmOoGyDk3G6hqgAguwLqtluveANztnNsbPK6wETGxc4/XsdXa1gEK9egnKznmzx+weedev0MRiTozu6qRx080syVmll9XpSgzO97M5pnZIjP7qLZjIiI1A344B679oMXeQiQeNbZdCENv9g2zA+9Gbe96treojqkduX3s7S39NiJxKawEyDn3DjAU+C1eGez3gDuBYc65txvxfuE0Ik1paIYCx5jZZ2b2kZkd0YiYMDNy+2eRnNR6h5kdPyybssoAz36+xu9QRPzw23APDOlpPhU4BLjIzA6pcUwn4N/AWc65EcD5kQu1hg1fwvovWuzyInEs7HYhTLXdBXX1bK/9ImbXmlmemeVt3ry5SYHsKNvB2p1rmVc4r0nni0j9wv7G75wLOOdmO+deds69FHxe2cg7MOE0Io1qaIKSgCxgPPAz4IXgXKIDA6ilYTooO4OXbjiSDmnJDbyNfw7KzuC4odn83+zV7K2o5IKHZnHBQ7P8DkskYsxsQR2Pr4DujbhUOD3NFwOvOOfWQON7jRvlwz/Da1rHQ6QpItguhGMd0DfkdR+8RVfr2l4r59xk51yucy43Ozu70UHMK5zHsu3LKCgp4Pvvfl9JkEgLiESXR2PuwITTiDSqoQk555XgsLnP8eYlda3twOY2TH665phBbN65l1e/WO93KCItoTtwOXBmLY+tjbhOOL3IQ4EsM/vQzOaa2eV1XazZd3O3LIWuQxp/nohA5NqFcEwBLjfPeKDYObcRb/HVIWY20MxSgAuDx7aIvII8XPC+b2WgkryCvJZ6K5G4FW4Z7AV17aJxd2CqGxFgPV4jcnGNY6YANwcnGY5jXwNUn9eAbwEfBqvFpABbwg1qyvwN/Gv6Mp69djxdM1pfKewqRw3uwqG9O/LM52tol5zodzgikfYmkOGcm1dzh5l92IjrhNOLnASMAU4E2gGzzGy2c27pASc6NxmYDJCbm9tQb/T+Kspg+yoYeW6jThORapFqFzCzZ4Hjga5mtg74DZAM4JybBEzFq0Cbj1eF9qrgvgozuxlvCkAi8JhzblHTPk7DcrvnVj9PSUzZ77WIREa4VeC6A6cA22tsN2BmuG9WVyNiZtcH99fZAEHtjZdz7lHgMeAxM1sIlAFXBKu4hKWgeA/LCneR2ornAIE3V+neC0aRnZnGtU/pjpDEFufcD+rZV/NGSX3C7Wne4pzbDew2sxnAKOCABKhZtq8EVwldh0b0siLxIoLtAs65ixrY74Cb6tg3Fe/7SYvL6ZZDamIqDsfDJz9MTrecaLytSFwJNwGK2B2Y2hqRYOJT9by+BqjWxis4zv/SxsQRaueecswgPSXcX4d/BnfLBGBHaTk795Qzd/V2xvTP8jkqkVYlnJ7m14H7zSwJr8d4HHBvxCPZEsynNARORMJUHiinrLKMbu27KfkRaSFhfeOP5B2Y1mjHngoyUpNIaMVlsEO9tWAD32zaCcAlj8zm6avHKwkSCQqnp9k5902wlP8CvDmDjzjnFkY8mIHHwpVvQbZWcheR8DjnuOfYexjQcYDfoYjErNbf5REFO/dUtOoKcDUt37yr+nl5RYDZK7YqARIJ0VBPc/D1X4C/tGggaR1hwNEt+hYiEltSElOYOHCi32GIxLTWPeklSoZ2z+DYobUWjWuVjhqcXT3LOzEhgfGDuvgaj0ikmNlBZvammbUL2XaXmdXZC92qzXsWVs7wOwqRNi3m2oUGLNqyiPmb5/sdhkhMUwIEXHfcQfzp3MP8DiNsY/pnMbxHJgkGg7LTGd2vk98hiUSEc2453vyc982si5n9CxgEPOFrYE31/p0w/3m/oxBp02KuXWjA5AWT+fWnv/Y7DJGYFnYCFG93YFq7Du2S6de5PYs37WTm8kgvhSDiH+fcw8ADwHIgA7jMOVfpb1RNULYbdm2CzgP9jkSkzYuZdiEM+UX5DO402O8wRGJa2AlQLN+BOfv+T7hzSouV9G8x2Zmp/P17oxg3sLPfoYhETHChwfOBt/HW6unvb0RNtG2l97PzIH/jEIkBMdMuNKC0opS1O9cyJEuVI0VaUqOGwMXqHZh120sprwz4HUajJZhx7ug+JCVqJKPEBjPLwPuC82mw7P1NwFQzG+FvZE2wbYX3UwmQSLPEVLvQgBVFK3A4hnRSAiTSkhpVBa6OOzCrIh9WdO3cU0FmG6oCB/D8dROqn78+bz3Pz1nLU98fq2RI2rp2wIPOuZcAnHMfm9nFQAd/w2qC6gRIQ+BEmil22oUGLCtaBqAhcCItLOwEKHgH5nXgbefcX83sGLw7MOc759re+LGgsooAZZUBMlIT/Q6lydKSE5m5fCvPzlnLZeNjclSAxAnn3GbgpRrb5vkTTTNNuBlGnOOVwhaRJoupdqEBJ/U/if4d+tM3s6/foYjEtMZ0F1TdgfkreHdg8FZXb9N3YErLvBF87VPa7pJIJx/SnfGDOnPvtKUUl5b7HY6IACQmQZZuSIhI+NKT0zm82+EkJrTdm7IibUFjiiBsrup+Dtk2zzk3K/JhRdc5h/dmaPdMv8NoMjPjf08/hO0lZfz7g3y/wxERgPd/C8ve9zsKEWlDHlv4GAs2L/A7DJGYF/cTRjq2T+beC3I4ekjbWQi1NiN7d+S7o/vw+KerKNyxx+9wRJrNzI42s5vMbFDItrYxoaZ8D3xyL6yf63ckIjGlTbcLDSjeW8y9c+8lryDP71BEYl7bHfclB/jpKcM4K6cX3Tqk+R2KSCRkA2OBsWa2FXgG+Alwka9RhaNoNeBUAEEk8tpuu9CA/CJvBIcqwIm0vLjvAfo0fwvD/vdt5q7e5ncozda9QxrHDMkGYG9Fm69OLnHOOfcq8H280vsLgWOApb4GFS6VwBZpEW26XWhA/vZgAqQ1gERaXJN6gMzsaGAUXkW4FcFtA51zKyMZXDTs3lvB3ooAqUmxM+HwP7NXM3nGCqbecgwZqerkk7bBzH4NlDjn/la1LbjO2OfBR9uhBEgkImKqXWhAflE+6cnpdG/f3e9QRGJeU78dx0wXdEl1FbjYSYAO7tmBtdtL+Pt7S1m0oRjYf90gkVbqMiCn5kYzuxrIds79KeoRNVXJVmiX5T1EpDlip11owLpd6xjUcRBm5ncoIjEvrCFwZvZrM/tJ1etY6oLetbcCgPQY6ikZ0z+Li8f244mZK9kd/HwibUCpc66klu3/B1wa7WCa5cQ74KfLQF9kRJorou2CmU00syVmlm9mP69l/8/MbF7wsdDMKs2sc3DfKjP7Krgv4pUKHjjxAR789oORvqyI1CLcOUCXAfv9qwx2QR8GdHfO3euc+02kg4uGkrLYS4AAbps4nK4ZqSzfvJuAc36HIxKOUjPrWXOjc24v0PYy+cRkvyMQiQURaxfMLBHvxu2pwCHARWZ2SI3r/sU5l+OcywF+AXzknAudJHxCcH9uIz9HgxIsgY6pWjhZJBrCTYDqugPzFG3tzmwNw3t04KKx/WiXHDtD4AA6tkvmz+cdRml5Jau27Gbu6u1+hyTSkL8Br5vZfquHmlk3IOBPSE1QWQ7PXAhL3/M7EpFYEMl2YSyQ75xb4ZwrA54Dzq7n+IuAZxv5Hk2yomgFv5n5G9bsWBONtxOJe2EnQHXcgSmjLd6ZDXHs0Gz+dO6hJCbE3lCVDu2SMYPNu8q45JHZSoKkVXPOvYh3d3aumb1pZr83sz8CnwJ/9Te6RiheB0vfhl0Ffkci0uZFuF3oDawNeb0uuO0AZtYemAi8HBoO8J6ZzTWza+t6EzO71szyzCxv8+bNYQW2aOsiXln2ChWuTX+lEmkzwk2AYuPObC32VlQSCMTmELHZK7ZSNfqtrCLAx8vCa4hF/OKcexIYCLwAJAN7gIucc0/7GlhjFK32fmb1r/84EQlLBNuF2u501vUF4Ezg0xrD345yzo3GG0J3k5kdW0e8k51zuc653Ozs7LACW1m8kiRLom9m37COF5HmCWvii3PuxeDdkLlmNhuYh5c8nQ/c2WLRRcHtLy1g3toiPvzZCX6HEnHjB3UhwSDgvMeCtUV+hyRSJzMbCxyN1748HZxniJl1MrMM59wuP+MLW1FwCEunfv7GIRIDItwurANCM4w+wIY6jr2QGsPfnHMbgj8LzexVvCF1Mxrx/nVatWMVfTL7kJyguYMi0RD2QqiRugMTRgUWM7N/BvcvMLPRIfseM7NCM1tYx7V/ambOzLqGG8+uvZW0T4mtAghVxvTPYniPTPpkteO7Y3rz3yWbeeWLdX6HJVKXy4BhwOHAU2b2v2aWiXfX9hFfI2uM7avBEqFDH78jEYkFkWwX5gBDzGygmaXgJTlTah5kZh2B44DXQ7alB98XM0sHTsarghsRK4pWMLDjwEhdTkQa0Khv/s65nXiFD5okpALLSXh3YuaY2RTn3Nchh50KDAk+xuFVnxsX3PcEcH9tMZhZ3+B1GzWDsKSsgvTU2CqAECozLZnMtGTuPvcw1mwr5X9fW8hhfToxuFuG36GJ1HQ7cD3euPtEvLH5N+CNwQ/7Zo3vEpKg1+GQGJs3VkSiLGLtgnOuwsxuBt4NXusx59wiM7s+uH9S8NBzgPecc7tDTu8OvBpcoycJeMY5906TP9X+cWFmDO40OBKXE5EwRPt/6OoKLABmVlWBJTQBOht4yjnngNnBbu6ezrmNzrkZZjagjmvfC9xGyB2bcOwuq6RTu9jvck5KTOCfFx7Oaf/8mJuf+YI3fng0yYlt5zulxL5gpcm/A383s1S8myDZQC9gdH3ntion/MJ7iEizRbpdcM5NBabW2Dapxusn8G64hm5bAYxq7PuFw8x49exXcVqyQiRqop0A1VaBZVwYx/QGNtZ1UTM7C1jvnJvf0ArKwcot1wL069ePTnsr6N0pLewP0NY8f92E6uc9OqZx7wU5bCwqVfIjrVpwjY/q4SXBmyVhM7OJwD/w7vI+4py7u47jjgBmAxc4515qesQi0tKa2y60dg19fxGRyIn2t+BwKrA0pkpLVanKXwF3hBNAzeosFx7Rl4kjD6jwHbOOG5rNhWO9ydnbdpf5HI1IeKomPocjnMUOQ477M95wmMgo3wOTjoGFr0TskiJSu8a0C63Zy0tf5vpp11NeWe53KCJxI9oJUDgVWBpTpQXgILziDPPNbFXw+C/MrEc4AV19zCDOGtUrnENjSt6qbRzz5/8y/RutVSIxJ9zFDn+IN4+gMGLvXLwWNi2ASt1cEJHwzN88n8XbFpOcGPvD8UVai2gnQOFUYJkCXB6sBjceKHbO1Tn8zTn3lXOum3NugHNuAF4CNdo5tymcgAp37mFPeUzcRGqUEb06MjA7nVuem8eygp1+hyMSSQ0udmhmvfEmOu839r/ZtgfXAOqkNYBEJDwri1eqApxIlEU1AXLOVQBVFVi+AV6oqsBSVYUFb3LiCiAfeBi4sep8M3sWmAUMM7N1ZvaD5sUDY/8wnUc+XtGcy7RJ7VISmXxZLu1SErny8TkU7NjDBQ/N4oKHZvkdmkhzhTOM9j7g9nCG0DRqVXctgioijeCcY0WxSmCLRFvU67Q2VIElWP3tpjrOvSiM6w8IN5ZAsOJKrK4D1JBendrx+JVHcMFDs7jisc/JSE0iMUGTMKXNC2cYbS7wXHDScVfgNDOrcM69VvNizrnJwGSA3Nzc+ss0Fa2GxBTICGsErojEue17t7OjbAcDOgzwOxSRuBLXpcCqEqB2KbG7DlBDRvbuyKTLxjDhoC7s3lvB+qJS5q7e7ndYIs3R4FBb59zAkGGzLwE31pb8NFpGDxh6CiTEddMqImHaXb6bsT3GMrzzcL9DEYkrcf2/dCB4L7d9HCdAAMcMyeaMw3qxpGAn67aXcskjs5UESZsV5lDbljHhRrjgPy36FiISO/pm9uXRUx5lbM+xfociElfic+xXUFUPUFpyfCdAALNXbK1OCPeUB5i1fAtj+mf5G5RIE4Wz2GHI9iujEZOIiIi0DnHdA5SckMDtE4czrHum36H4bvygLoRO//lyTRGBgFalFgnb3l1wzyD4Uj1AIhKeX33yK65/v2U7pkXkQHGdACUlGjccfxADuqb7HYrvxvTPYniPTPpkteP83D5MX1zIHVMW4pySIJGwFK2Bkq2QlOZ3JCLSRizbvqyepd5FpKXE/RC4NVtL6NkpjeTEuM4FAchMSyYzLZl7zjuMzukpPPTRCgZ0SefqYwb5HZpI61e0xvuZNcDXMESkbXDOsWbnGnK65fgdikjciesEaEdpBcf+5QOm/+Q4DsrO8Dsc3z1/3YTq5z+fOJw+ndpx9uG96zlDRKoVaRFUEQnf1j1b2V2+m/4d1GaIRFtcd3tUDe9qpyIIBzAzLpswgA5pyZSWVfL3aUvZU97gmpEi8atoDSS3h/SufkciIm3Amh1er3G/zH4+RyISf+I6AQooAQrLJ/lb+Of0ZVz9ZB6791b4HY5I69R9BIy+HEyLCYtIw9KT0znroLMY3Gmw36GIxJ04T4C8n/G8EGo4TjqkO389fxSzVmzlgsmzKNyxhwsemsUFD83yOzSR1iPnYjj1z35HISJtxLDOw/jD0X+gZ0ZPv0MRiTtxngB5GVBqUlz/GsLy3TF9eOTyXFZs3s05/57J5p17WV9UqgVTRaqUlfgdgYg0wMwmmtkSM8s3s5/Xsv94Mys2s3nBxx3hnttYu8t3q9KqiE/i+pt/Zloyv/vOSExDVsJywvBuvHDdBBLNWLl1N+u2l3LJI7OVBImUFsEfe8Jnk/2ORETqYGaJwAPAqcAhwEVmdkgth37snMsJPu5q5Llhu/KdK/nxhz9uziVEpIniOgFqn5LIZeNVfaUxRvbuyPeO6EPVTauyigCzV2z1NygRv1VVgMvs7m8cIlKfsUC+c26Fc64MeA44OwrnHsA5x+odq+mR3qOplxCRZojrBKisIkB+4S6/w2hzJhzUlYRgp1nAwdKCnVQG1I0vcWy7SmCLtAG9gbUhr9cFt9U0wczmm9nbZjaikeeGZUvpFkorSlUBTsQncZ0Abdqxh+v/M9fvMNqcMf2zGN4jk96d0jh1ZA9en7eBKx//nOKScr9DE/FH1SKonfRlRqQVq228e827d18A/Z1zo4B/Aa814lzvQLNrzSzPzPI2b95cayCrd3g3TbQGkIg/4joBCjinEthNlJmWTJ+s9jx46Rj+dO6hzF6xlbMf+IT8wp1+hyYSfUWrIbUDtMvyOxIRqds6oG/I6z7AhtADnHM7nHO7gs+nAslm1jWcc0OuMdk5l+ucy83Ozq41kDU7g2sAddBNExE/JPkdgJ8CAa0B1FTPXzeh+vlFY/sxuFsGNz/zBYU79jK4W6aPkYn4YNAJ0KGX1gASad3mAEPMbCCwHrgQuDj0ADPrARQ455yZjcW7UbwVKGro3MYY1nkY1xx6DT3TVQJbxA/xnQA5R5rWAIqIIwZ05qOfnUBaMKH8YHEhRw3uSopKjEs8GH4acJrfUYhIPZxzFWZ2M/AukAg85pxbZGbXB/dPAr4L3GBmFUApcKHzalXXem5TYxnRZQQjuoxo+EARaRFxnQA5B+3VAxQxVclPfuEuvv/kHEb16cQDl4zm1ufnAfv3GonEDOdg63Lo1BeSUv2ORkTqERzWNrXGtkkhz+8H7g/33KbK355Pj/QeZKRkROJyItJIcX17vkfHNK46aoDfYcScwd0yeODi0eQX7uL0f37MhqJSLZoqsatkK9w/BvIe8zsSEWkDnHNcPPViHpj3gN+hiMStuE6AMtOSGDeoi99hxKTTDu3Jmz88mk7tklm7vVSLpkrsUglsEWmEwpJCSitKVQFOxEdRT4DMbKKZLTGzfDP7eS37zcz+Gdy/wMxGh+x7zMwKzWxhjXP+YmaLg8e/amadwoll994K1m4rafZnktoN6JrOOYfvWyahXIumSiwqWuX9zNKXGRFpmCrAifgvqgmQmSUCDwCnAocAF5nZITUOOxUYEnxcCzwYsu8JYGItl54GjHTOHQYsBX4RTjwrt+zmuTlrGvMRpJGOHpJdvWhqclICgYDjqVmr8OaUisSA6h4gfZkRkYZpDSAR/0W7B2gskO+cW+GcKwOeA86ucczZwFPOMxvoZGY9AZxzM4BtNS/qnHvPOVcRfDkbrz5/gxwqg93SqhZN7ZPVjqevHs+ywl3c8foivv/EHDbv3Ot3eCLNV7Qa2nWGVJV/F5GGrdmxhuSEZHq07+F3KCJxK9pV4HoDa0NerwPGhXFMb2BjmO/xfeD5cANKUwLU4qbecmz189H9OjGmfxZ/nPoNE++bwR/OGcnEkVoHQdqwwy6EfqpwKCLhOX3Q6QzvPJzEBH3/EPFLtHuAalslsOZYqHCOqf3iZr8CKoCn6znmWjPLM7M8gHZaByiqzIwrjhzAGz88mh4d07j+P1/w5Zp9hREueGgWFzw0y8cIRRqp/wQYdaHfUYhIGzGs8zBOG6R1w0T8FO0EaB3QN+R1H2BDE445gJldAZwBXOLqmWDinJvsnMt1zuUCpCUpAfLD0O6ZvHbTUTx4yWgO75cFwPLNu9hRWq6S2dJ2BAKw4kPYvcXvSESkDXDOMX3NdAp2F/gdikhci3YCNAcYYmYDzSwFuBCYUuOYKcDlwWpw44Fi51y9w9/MbCJwO3CWcy7ssm4DurRnwkEqg+2X5MQETj3UG/62cstuTrl3Bt9s2qmS2dJsYVSbvCRYNXKBmc00s1FNeqOdG+Gps+Hr15sds4jEvsKSQn70wY/4cO2HfociEteimgAFCxXcDLwLfAO84JxbZGbXm9n1wcOmAiuAfOBh4Maq883sWWAWMMzM1pnZD4K77gcygWlmNs/Mqld1rk9mWjK9OrWLxEeTZuqb1Y6jBnetfr23PMCs5bqrLo0XZrXJlcBxwcqRvwMmN+nNioIV4FQCW0TCUFUCu2+Hvg0cKSItKdpFEHDOTcVLckK3TQp57oCb6jj3ojq2D25KLDv2lLN9dxlZ6SlNOV0iKCkxgf85cQgzlm3GOW/S11tfbeT64w4iKTGu1+uVxquuNglgZlXVJr+uOsA5NzPk+LArRx6gugT2gCadLiLxZe1Or8ZTv0yVzRfxU9QToNZk9dYS8jfv4oj0zn6HInglsw/ukcmO0nK+M7oPqYkJ1clPWUWAlCQlQhKWcKpNhvoB8HaT3qloNWDQSXdzRaRhq3esJikhiZ7pqn4q4qe4ToAAUvWlulUJLZldZebyLfz0hfn8+oxDmDiyB2a1FQoUqRZ2JUkzOwEvATq6zouZXYu3KDP9+tW4a7t9NWT2hKTUpsYqInFk7c619MnooxLYIj6L+wRI6wC1fukpSXRol8wNT3/BuIGd+fUZh/C7N73RTM9fp/VX5ABhVZI0s8OAR4BTnXNb67qYc24ywTlCubm5+ydSR/9IJbBFJGy3jrmV7XtU4EfEb3Hf/aEeoNZvVN9OvPnDo/nd2SNYWrCTM+//hGUFu1QuW+rSYLVJM+sHvAJc5pxb2uR3evNW+Oie5sQqInHkxx/8mJ/N+BnzCuf5HYpIXIv7b//qAWobkhITuGzCAD786QmccWhPtpeUVZfLVrU4CRVmtck7gC7Av4OVI/Ma/UYVZbB9FWxbCWs/j1D0IhKrPl3/KUu2L2H9rvVc8941SoJEfBTXCdCgrulktVcFuLakY/tkhvfsUD2hY295gGufmsvzc9ZQURnwNTZpPZxzU51zQ51zBznn/hDcNqmq4qRz7mrnXJZzLif4yG30myx+E3asg53r4cmzlASJtAHNWSPMzFaZ2VdNvWkybfU0XPB/r/JAOXkFjb/vIiKREdcJUHpqkiqLtUHjB3UhITjNPTkxge4d07j95a84+d4ZvLlgA4FArfPdRSIrf/q+55VlsOpj/2IRkQZFaI2wE5p60yQrLav6eXJCMrndG3/fRUQiI66LIBSVlPsdgjTBmP5ZDO+RyY49FfzjwsMZ3a8T074u4K/vLeHmZ77krFEF/POiw/0OU2Jdepd9zxNTYMAx/sUiIuGI3hphtXDOkUgiN+TcwLie48jplhOpS4tII8V1ArRue4nfIUgT1SyXffKIHpx4cHemzF9Pl3SvJHFxaTmzV2zlvmlL2bnXS5bG9M+q7XIijWcJYIlwwi9h4LHQd6zfEYlI/Zq7RpgD3jMzBzwUrBAZtjU719CnQx+uG3VdY04TkRYQ1wmQ1pOJLYkJxjmH77tZ99LcddXlsgEueWQ2T189XkmQRMb21dCpHxz7U78jEZHwNHeNsKOccxvMrBswzcwWO+dm1HJurWuHrdmxhn6Z/WoeLiI+iOsJMAnKf2LaFRP6c8Zh+1bb3lMeYPKM5TinOUISAaf9BS56zu8oRCR8jV0j7OzQNcKccxuCPwuBV/GG1B3AOTfZOZfrnMvNzs6u3v7UqU/x2yN/G4nPISLNFNcJkHqAYltSYgJXHTWwOtE1g4Ide6v/3ItLNQdMmiG9K3Qb7ncUIhK+Jq8RZmbpZpZZ9Rw4GVjYmDdvn9ye7PbZDR8oIi0urofAqQco9u1XMOGCHIb0yAS8+V/f+ttHnDKiB1cdNYDR/TQsThphTzHMngQjvgPZw/yORkTC4JyrMLOqNcISgceq1ggL7p/E/muEAVQEK751B14NbksCnnHOvRPuey8vWs7r+a9z8cEX0yO9R0Q/l4g0XlwnQP06p/sdgkRBzYIJAClJCVw+vj/P563ljfkbOLR3Ry4c25f/m7WaXSqYIA3Zkg8f/hF6HKoESKQNcc5NBabW2DYp5PnVwNW1nLcCGFVze7gWbF7A44se5/yh5zf1EiISQXE9BC4tOa4/flzrlpnG/55xCLN/cSJ3nT2C8soAd7y2kCUFO1m3vZSLH57N3FXb/A5TWqtty72fXQ7yNw4RaRPW7FxDkiXRI0O9PyKtQVxnADv3aA5IvEtPTeLyCQN4+5ZjuOrogVTVR9hbEeCa/5vLgx8uV7l0OdDW5YBBp/5+RyIibcCq4lX0yexDckKy36GICHGeAG3dVeZ3CNJKmBmnjuxZPS8sOdHompHCn99ZzNF//oDvPjiTdxZu9DdIaT22LYeOfSE5ze9IRKQNWFm8koEdB/odhogExfUcIFWBk1Bj+mfx4vVHMnvFVsYP6sKY/lms2VrCGws28NqX61lftAeA4pJy3lm0kUc/WUlJWaXmC8WjojXQZZDfUYhIGxBwAXaU7VACJNKKWDyvidJ90CGuYMXXDR8occ85R8B5i62+uWADNz/zZfW+pATjHxfmcPphvXyMsHUxs7nBykkxJTc31+Xl5UGgEvbuhHad/A5JJCbEepvhnKPCVWgInEiENLfNiOshcOoBknCZGYnB8XGnH9qTKyYMqN5XEXDc9MyXrN3mzRXavHMvZRUBP8KUaElIVPIjImEzMyU/Iq1IXCdAWgdImsLMOCunV/Xfn9SkBG46/iD6dm4PwG/fWETOXe/xgyfm8OTMVXz7bx9y9N3Tmbt6u49RS8QUfA1T/ge2r/I7EhFpA95Y/ga3zbiN8oAKL4m0FlFPgMxsopktMbN8M/t5LfvNzP4Z3L/AzEaH7HvMzArNbGGNczqb2TQzWxb8GdaEjG4dNIFZmqZqvtDPThnGM9eM52cTh1fvu+CIvpw3ug/5m3fxmymLyN+8m3VFe7jkkdnMXb2dLbv2Es9DT9u89XnwxZPg1MsnIg2bs2kOczbNUQ+QSCsS1SIIZpYIPACcBKwD5pjZFOdc6EScU4Ehwcc44MHgT4AngPuBp2pc+ufAdOfc3cGk6ufA7Q3Fk6QuIGmGMf2zai1+cMyQbI4Zkg3AH976moc/XglAeUWAmcu3cPmjy+nYLplxg7owbmBnxg7szMCu6ZgZp/1jBjv2aCHWVq3ga0huD50G+B2JiLQBqgAn0vpEuwdoLJDvnFvhnCsDngPOrnHM2cBTzjMb6GRmPQGcczOA2lanPBt4Mvj8SeA74QRTUlbZ+E8g0ggTQ0trJyWQ2z+L2yYO5/B+WXy8bDM/f+UrvvW3j/j3h8uZu3o732wKWYhVQ+Zap8KvIXs4JMT1CGIRCdOK4hUM6DDA7zBEJES0y2D3BtaGvF7Hvt6d+o7pDdS3CEt359xGAOfcRjPrFk4wJWUV4Rwm0mS1ldaecFBXrjhyAM45lm/ezecrt3F4v078d3Hh/guxPpXHt4Z3I6dvJyaO7EHXjFR/P4x4Cr+BoSf7HYWItAHlgXL2lu1lSNYQv0MRkRDRToBqG3NWczJEOMc0PQCza4FrAbr2Vpe0tLy6hsqZGYO7ZTC4Wwbg9UiagXPe8MxBXdP5YHEhL81dx6g+neiakcoHSwp5c/5GRvTqwJMzV1EeCPCvi0ZruFy0uACkpEP3kX5HIiJtQGWgkmFZwzi488F+hyIiIaKdAK0D+oa87gNsaMIxNRWYWc9g709PoLCuA51zk4HJAAcdfJhmokurMaZ/Fgf3yNxvDpBzjrXbSunR0SvYsbFoDx8t3czLX6yrPu+7D87kye+P5dih2azcspsEg75Z7UnQHLfIswS45QtQEQsRCUNaUhovn/Wy32GISA3RToDmAEPMbCCwHrgQuLjGMVOAm83sObzhccVVw9vqMQW4Arg7+PP1cILROkDS2ky95dj9XpsZ/bq0r3598bh+XDyuH/e8s5h/f7gc8LpHv1pfxLFDs7l32lKmzN9Au+REBnfLYGj3TEb27sCLeWtVXCGS1HaIiIi0WVFNgJxzFWZ2M/AukAg85pxbZGbXB/dPAqYCpwH5QAlwVdX5ZvYscDzQ1czWAb9xzj2Kl/i8YGY/ANYA54cTj26QS1t14sHdmfTRcgIO0pITGD+oKwA3nnAQRw/uypKCnSwt2Mkn+ZvJW7WNtdtLCDg4f9JMBnfL4PC+WQzKTmdQdgbDumful2SpEl09itfCGz+CM+/zOxIRaQNWFq/kgXkPcFPOTX6HIiIhot0DhHNuKl6SE7ptUshzB9TaUjjnLqpj+1bgxMbGkpmmmvzSNtVWXAFgeI8ODO/RYb9j//H+Uu59fxngjdwqLavk/W8K2JpXBsC3hnfjsSuPAOD7T3zO1xt3AnDR5Nn8+9LRnDi8m3pLq+zZAXt3+B2FiLQRpRWlJFqi32GISA1RT4BaE32nk7asruIKNR09JJsHP1pOeUWA5KQE7gv27BSXlLNiyy4Sg12hlQHHV+v3fbkvqwxw9ZN5XDGhP789eySVAcc/pi+jX+f29O/Snl+98hUlZRX8oxUWYTCzicA/8HqaH3HO3V1jvwX3n4bX03ylc+6LBi9cWQZ9joh8wCISFc1pGxo6tzYOx+HdDo/shxCRZovrBEgkHozpn8XTV48/oLeoY/tkDu+3L3FJTDAmXTqG8yfNJOAgJTGBS8f349sHdwegYMce7v/vMgI15v9fOHkWz107gYOy05kyfwNDumVG7bPVJgILLtdviEpgi7RFzWkbwjz3wPfEOCz7sMh/GBFpFiVAInEg3N6iuobWAfTq1I7FvzuVddtLeODDfF6eux7weo5mr9hKRWWAO15fRE7fTi31McJVveAyQLCgytlA6BeV6gWXgdlm1qmqkmS9V05Kgy4HtVDYItLCmtw2AAPCOPcA7ZPa0y6pXcQ/iIg0j5YyF5H9jOmfxU0nDK41YUpJSmBQdgYXj+1fXUQkJSmB8YO6MHZgZz775Yn87XujohzxAepaTLmxxxwoI7u5sYmIf5rTNjSpzchKa13Dg0XEox4gEWm0unqKundIo7vPsRHhBZdDF0/u169f8yITET81p21QmyESQ5QAiUiThDuszgcRXXA5dPHk3NxcrYAq0nY1p21ICeNcQG2GSFugIXAiEmuqF1w2sxS8BZen1DhmCnC5ecYT3oLLItK2NadtCOdcEWkj1AMkIjGluQsui0hsak7bUNe5PnwMEYkAJUAiEnOas+CyiMSuZi7GfsC5ItI2aQiciIiIiIjEDSVAIiIiIiISN5QAiYiIiIhI3FACJCIiIiIiccO8+X7xycw2A6v9jiOoK7DF7yBCKJ76KZ76DXPOZfodRKSpzaiX4qmf4qmf2oyW19r+zBVP/RRP/ZrVZsR1FTjnXLbfMVQxszznXK7fcVRRPPVTPPUzszy/Y2gJajPqpnjqp3jqpzaj5bXGP3PFUzfFU7/mthkaAiciIiIiInFDCZCIiIiIiMQNJUCtx2S/A6hB8dRP8dSvtcUTi1rb71jx1E/x1K+1xROLWtvvWPHUT/HUr1nxxHURBBERERERiS/qARIRERERkbihBChKzOwxMys0s4Uh2zqb2TQzWxb8mRWy7xdmlm9mS8zslAjH0tfMPjCzb8xskZnd4nM8aWb2uZnND8bzWz/jCXmPRDP70sze9DseM1tlZl+Z2byqyic+x9PJzF4ys8XBv0cT/P7zijVqM+qNR21Gw7GozYgzajPqjUdtRsOxxFeb4ZzTIwoP4FhgNLAwZNs9wM+Dz38O/Dn4/BBgPpAKDASWA4kRjKUnMDr4PBNYGnxPv+IxICP4PBn4DBjvVzwhcd0KPAO86eefV/A9VgFda2zzM54ngauDz1OATn7/ecXaQ21GvfGozWg4FrUZcfZQm1FvPGozGo4lrtoM3/6hxuMDGFCjYVoC9Aw+7wksCT7/BfCLkOPeBSa0YFyvAye1hniA9sAXwDg/4wH6ANOBb4U0TH7GU1vD5Es8QAdgJcE5hH7HE8sPtRlhxaI2o/Z41GbE4UNtRlixqM2oPZ64ajM0BM5f3Z1zGwGCP7sFt/cG1oYcty64LeLMbABwON7dEN/iCXYDzwMKgWnOOV/jAe4DbgMCIdv8jMcB75nZXDO71ud4BgGbgceDXfePmFm6j/HEE99/x2oz6nQfajPqojbDP77/jtVm1Ok+1GbUpcXbDCVArZPVss1F/E3MMoCXgR8553b4GY9zrtI5l4N3R2SsmY30Kx4zOwModM7NDfeUlown6Cjn3GjgVOAmMzvWx3iS8IZZPOicOxzYjdcV7Vc8ojZDbcaB1GZIfdRmqM2oKa7aDCVA/iows54AwZ+Fwe3rgL4hx/UBNkTyjc0sGa9Reto594rf8VRxzhUBHwITfYznKOAsM1sFPAd8y8z+42M8OOc2BH8WAq8CY32MZx2wLnj3DOAlvIbK978/cUBtRg1qM2qnNkOC1GbUoDajdvHWZigB8tcU4Irg8yvwxshWbb/QzFLNbCAwBPg8Um9qZgY8CnzjnPt7K4gn28w6BZ+3A74NLPYrHufcL5xzfZxzA4ALgf865y71Kx4zSzezzKrnwMnAQr/icc5tAtaa2bDgphOBr/2KJ86ozUBtRkPUZkgItRmozWhIXLYZkZqwpEeDE7qeBTYC5XiZ6g+ALngT4JYFf3YOOf5XeFUslgCnRjiWo/G6BhcA84KP03yM5zDgy2A8C4E7gtt9iadGbMezb3KiX7+fQXjVTeYDi4Bf+f37AXKAvOCf2WtAVmv484qlh9qMeuNRm1F/DGoz4vChNqPeeNRm1B9D3LUZFjxJREREREQk5mkInIiIiIiIxA0lQCIiIiIiEjeUAImIiIiISNxQAiQiIiIiInFDCZCIiIiIiMQNJUAiIiIiIhI3lACJiIiIiEjcUAIkrV5wheL7zWy837GISOunNkNEGkNtRvxRAiRtwfVAKt7K0iIiDVGbISKNoTYjzigBkrZgIrAUmOdzHCLSNqjNEJHGUJsRZ5QASatmZmlAIjAa+MjncESklVObISKNoTYjPikBktZuCF7DtNg5V+53MCLS6qnNEJHGUJsRh5L8DkCkAdnAUOBsvwMRkTZBbYaINIbajDikHiBp7XoBLwMJZpbldzAi0uqpzRCRxlCbEYeUAEmrZWZJeGNyewCTgEp/IxKR1kxthog0htqM+GXOOb9jEBERERERiQr1AImIiIiISNxQAiQiIiIiInFDCZCIiIiIiMQNJUAiIiIiIhI3lACJiIiIiEjcUAIkIiIiIiJxQwmQiIiIiIjEDSVAIiIiIiISN/4fkuPOmKgRZBgAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(1, figsize=(3.*4.5, 3.75))\n", + "\n", + "plt.suptitle(r'SO $\\times$ unWISE')\n", + "\n", + "plt.subplot(131)\n", + "plt.title(r'$gg$')\n", + "ell, cl, cov = s_load.get_ell_cl('cl_00', 'gc_unwise', 'gc_unwise', return_cov=True)\n", + "cl_err = np.sqrt(np.diag(cov))\n", + "plt.plot(ell_theory, 1.e5*(cl_gg_theory + Nell_unwise_g), '--', c='C0')\n", + "plt.plot(ell, 1.e5*cl, 'o', ms=3, c='C0')\n", + "plt.errorbar(ell, 1.e5*cl, yerr=1.e5*cl_err, fmt='none', c='C0')\n", + "plt.xlim([1,ell_max])\n", + "plt.xlabel(r'$\\ell$')\n", + "plt.ylabel(r'$C_\\ell \\times 10^5$')\n", + "\n", + "plt.subplot(132)\n", + "plt.title(r'$g\\kappa$')\n", + "ell, cl, cov = s_load.get_ell_cl('cl_00', 'gc_unwise', 'ck_so', return_cov=True)\n", + "cl_err = np.sqrt(np.diag(cov))\n", + "plt.plot(ell_theory, 1.e5*ell_theory*cl_gk_theory , '--', c='C1')\n", + "plt.plot(ell, 1.e5*ell*cl, 'o', ms=3, c='C1')\n", + "plt.errorbar(ell, 1.e5*ell*cl, yerr=1.e5*ell*cl_err, fmt='none', c='C1')\n", + "plt.xlim([1,ell_max])\n", + "plt.xlabel(r'$\\ell$')\n", + "plt.ylabel(r'$\\ell C_\\ell \\times 10^5$')\n", + "\n", + "plt.subplot(133)\n", + "plt.title(r'$\\kappa\\kappa$')\n", + "ell, cl, cov = s_load.get_ell_cl('cl_00', 'ck_so', 'ck_so', return_cov=True)\n", + "cl_err = np.sqrt(np.diag(cov))\n", + "plt.plot(ell_theory, 1.e5*ell_theory*cl_kk_theory , '--', c='C2')\n", + "plt.plot(ell, 1.e5*ell*cl, 'o', ms=3, c='C2')\n", + "plt.errorbar(ell, 1.e5*ell*cl, yerr=1.e5*ell*cl_err, fmt='none', c='C2')\n", + "plt.xlim([1,ell_max])\n", + "plt.xlabel(r'$\\ell$')\n", + "plt.ylabel(r'$\\ell C_\\ell \\times 10^5$')\n", + "plt.subplots_adjust(wspace=0.25);" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/soliket/__init__.py b/soliket/__init__.py index d5d524ce..0e54220c 100644 --- a/soliket/__init__.py +++ b/soliket/__init__.py @@ -3,6 +3,8 @@ from .ps import PSLikelihood, BinnedPSLikelihood # noqa: F401 from .clusters import ClusterLikelihood # noqa: F401 from .mflike import MFLike # noqa: F401 +from .xcorr import XcorrLikelihood # noqa: F401 + try: import pyccl as ccl # noqa: F401 from .ccl import CCL # noqa: F401 diff --git a/soliket/constants.py b/soliket/constants.py index 7b049cf5..0b8858c6 100644 --- a/soliket/constants.py +++ b/soliket/constants.py @@ -2,6 +2,7 @@ C_M_S = constants.c C_KM_S = constants.c * 1.e-3 +C_HMPC = constants.c * 1.e-5 h_Planck = constants.h k_Boltzmann = constants.k elementary_charge = constants.e diff --git a/soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits b/soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits new file mode 100644 index 00000000..275bd27f --- /dev/null +++ b/soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits @@ -0,0 +1,932 @@ +SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T NMETA = 0 END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 24 / length of dimension 1 NAXIS2 = 0 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 3 / number of table fields TTYPE1 = 'id ' TFORM1 = 'D ' TTYPE2 = 'min ' TFORM2 = 'D ' TTYPE3 = 'max ' TFORM3 = 'D ' SACCTYPE= 'window ' SACCCLSS= 'TopHat ' EXTNAME = 'window:TopHat' END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 24 / length of dimension 1 NAXIS2 = 0 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 3 / number of table fields TTYPE1 = 'id ' TFORM1 = 'D ' TTYPE2 = 'min ' TFORM2 = 'D ' TTYPE3 = 'max ' TFORM3 = 'D ' SACCTYPE= 'window ' SACCCLSS= 'LogTopHat' EXTNAME = 'window:LogTopHat' END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 168 / length of dimension 1 NAXIS2 = 601 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 2 / number of table fields TTYPE1 = 'values ' TFORM1 = 'K ' TTYPE2 = 'weight ' TFORM2 = '20D ' TDIM2 = '(20) ' SACCTYPE= 'window ' SACCCLSS= 'Bandpower' SACCNAME= 5828661008 EXTNAME = 'window:Bandpower' END ????????? ? +? ? ? ??????????????????? ?!?"?#?$?%?&?'?(?)?*?+?,?-?.?/?0?1?2?3?4?5?6?7?8?9?:?;?<?=?>???@?A?B?C?D?E?F?G?H?I?J?K?L?M?N?O?P?Q?R?S?T?U?V?W?X?Y?Z?[?\?]?^?_?`?a?b?c?d?e?f?g?h?i?j?k?l?m?n?o?p?q?r?s?t?u?v?w?x?y?z?{?|?}?~??????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????? ? +? ? ? ??????????????????? ?!?"?#?$?%?&?'?(?)?*?+?,?-?.?/?0?1?2?3?4?5?6?7?8?9?:?;?<?=?>???@?A?B?C?D?E?F?G?H?I?J?K?L?M?N?O?P?Q?R?S?T?U?V?W?X?Y?Z?[?\?]?^?_?`?a?b?c?d?e?f?g?h?i?j?k?l?m?n?o?p?q?r?s?t?u?v?w?x?y?z?{?|?}?~??????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????? ? +? ? ? ??????????????????? ?!?"?#?$?%?&?'?(?)?*?+?,?-?.?/?0?1?2?3?4?5?6?7?8?9?:?;?<?=?>???@?A?B?C?D?E?F?G?H?I?J?K?L?M?N?O?P?Q?R?S?T?U?V?W?XXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 16 / length of dimension 1 NAXIS2 = 0 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 2 / number of table fields TTYPE1 = 'name ' TFORM1 = 'D ' TTYPE2 = 'quantity' TFORM2 = 'D ' SACCTYPE= 'tracer ' SACCCLSS= 'Misc ' EXTNAME = 'tracer:Misc' END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 16 / length of dimension 1 NAXIS2 = 0 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 2 / number of table fields TTYPE1 = 'name ' TFORM1 = 'D ' TTYPE2 = 'quantity' TFORM2 = 'D ' SACCTYPE= 'tracer ' SACCCLSS= 'Misc ' EXTNAME = 'tracer:Misc' END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 16 / length of dimension 1 NAXIS2 = 3000 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 2 / number of table fields TTYPE1 = 'ell ' TFORM1 = 'K ' TTYPE2 = 'beam ' TFORM2 = 'D ' SACCTYPE= 'tracer ' SACCCLSS= 'Map ' SACCNAME= 'ck_so ' SACCQTTY= 'cmb_convergence' EXTNAME = 'tracer:Map:ck_so:beam' MAP_UNIT= 'none ' SPIN = 0 END ??,?چ?>?\-? +D??5L:?>g& ? +?z` ?݋ ?2' ?-:d?G#?R%?L ?5[??g[?ð>?7?V?VS?bl?3d??ϗ?Қ?)yl?>]w?B% ?6N +!?"?#?ׯ[$?aM %?} &?Д8;'?ܠ(?ˆ)?,3*?7$W+?vV/,?-?=.?6/?Gp0?R1?E2?i1.3?%J 4?j15?m>6?7?t8?+59?9Y:?;?k0???~@?zAA?v|sB?rVC?n:GD?iE?ev"F?`^G?\p H?WAI?S)DJ?NmcK?I +.L?D xM??ֺ2[N?: +@O?5P?0xQ?+}R?&>]S? T?0U?hV?3;W? X?klY?bZ?a[?$\?A]?N^?Jl_?6T`? k?Hjgl?_um?f!n?y\nkgo?rB[hp?kpq?cr?\s?U6Uft?Melu?FNv?>hy\w?7$[-x?/vy?'#bEz?{? o|?_}?@F~?_??2\?}pY'?,Om?ЍW?X?ֶ?D1#?#LO?s?*?V/?q<?|]h?wt?vb?m<D?dd?Z0B?Qh-G?HMN?> C?5u?+hw?!y?k?<?b"?x?}1?rS?Wd?+uu?杖?/:?FX> +?#?[?ͥ6i?/Z5?~>?sy?hQ?^z?S"}?H" ?c?=8[?1 b?&B?|y?*_f?l?T˼?xh{?=Ɂl?֙p?YBf? H?K{&?f+?p0?j?SzF?w,-?jWm<?^'?RUɢ!?Ez?9tG*?,Sk? Rn ???FZ?$6?In?^M?bC?V?:,?E?c1?'Q?&j?v1?i:[U?[ƿ?N H?@]k?2$?$͹?[W?F??,;|?kl ?ЖR3?S'?b??*z ??z~N?kmI?\#܎?Mu??j3?08$/?!BR?;n?$!?&O? )?~"C?%?HF?DY?^l?!wH?xwW?h<"c?XA?INW?9,]_?)0?$? y??^?QV?zC?,?M?y W??u*?dk+?S6?B ?1P< +? ^@ ? ?~U ?dH?"N?Ϲ+?lк??vt?\?s=C?a]_a?OĀ?=O ?, >?? ΅?S?І.?љf5?Q?0 ?!?"?ud#?b3v$?PNz%?=UZ&?* '?<(?)? +_*?D+?&%,?B)-?| u.?/Ҍ/?~Ұ +0?keu1?W&.2?DZM3?07ڹ4? 5? O6?7?=8?ͱ 9?.:?;?\V?T}??@[@?+,pA?fA&B?SqC?0тD?}E?ĹF?G?ިH? 7|I?qc3J?[CnK?FL?13M?RtN?yO?"P?6Q?ƖR?y_S?W +/T?񅟎dU?o'%V?YtlW?Dv]cX?.-Y?Z?u[?ΏZJ\?Քa]?I[AMϠ?##?%? p?)Q?;N&?<~?-v?gQnJ?K?0?O[5?u<?~j?U?lo?ˋ?p?TXld ?8y?W8??M]?Ȝm?yP) ?E ?t`'?W?;J?7Y?Qc???b?揝K?K?rǞ'V?Up ?8E?7Լ?ΡI??pƗ)?)Z?E?ljO$?N_A?1j7K?<?)<?pj?产`x?#?~䤪?`&?B +. ?$?J?b^??f?㫼C?Q?nթ?PJ g?17w?.u?E|?y{>?⶜p?◯?x?YJ?:{?Zm?O?qݯ?p&?D?~?^)??T ??1??69?F X?FI+Q?`5?@+?R>?Hz2?R ?߾?ߞe?};U?]k@mr?<!?B?T09?ځv?޹**?ޘvc?w1) ?V?5q ?=J0 +?L ?ѥ 6 ?ݰ@! ?ݎ@?mG=U?K ?* @z?Xg)?b?ĽpH?ܢ׳j?܀V?^۶/y? .?Y=;/?0?C1?ذZvn82?،s]3?hS#~4?Dʉ}5? ϛ`W6?Ĉ[ 7?ةP8?״}K9?אBuSQ:?kѶ;?G +~?&`<??ֵ@?֐pA?l!.B?GTKKC?"wD?<E?؍F?ճdG?Վc*'H?i5I?C0J?B,K?L_LL?\M?Ԯ`9N?ԈO?c3RXP?=*jQ?wER?>S?-T?Ӧ)mU?Ӏ*V?Z[W?3EX? b.Y?HZ?=ko[?Қ٠\?tm]?M^?']"_?i\`? Y\a?ѳLVb?ь|ggc?e0vd?> +|e?48f?E(Ug?x;˚h?ТGi?{&j?Sqk?,S(l?m?_Q~n?ϵ>^o?ώ,pNp?fzTq?> r?xs?M\t? |pu?Ο:v?vDձw?Nx?&&y?ew4z?%{?ͭ%|?ͅ<}?\M~?4&o? ?h?̺ ?̑9M?hV8dH??bX2?_ɭ?K_q?' R?˚ʃ8?q_&?H[7?~?@?x?ʢkmfE?x$H?O`T?%NbE?y,?ѕdaO?ɧk?}*?Sx:?)b7?-p? ?Ȫ$?Ȁ/.~h?U+<?+5+@?xv??ǫEE?ǀfH?U +?*K?f?Z?Ʃ1?~?SZq?(R?S?YR?ť%?z^N_?NȐ?#"˻I?l+?˧*b?ğO?smv?G„??٠U?ó?×}F?k7<?>j?{W7?q`}?¹Y2,?Œ=P?`C? +?3?ӟA?Y}? +~u?P8?S?%š?Ѻ5?˛&?UO3?p?C_??"?rƸ?;ޖ?ￍ` +?_2?1? i3?'8#?ﾨ1 M?z*B?L?i?f(?rE'd?ｓ'?dW)?6??MS?#DH?Ｊ}B?{Ӎ?Mn?+,?E?P]?ﻑJ$^?b4H\:?3s?ئi?Ԓ?ﺥ=#?und?Fa¯?o?Fm?ﹷ?﹇nD]?X%?(Pz?k?uVg?︘pi?h[$"?85?*?׻?﷧f?w}?F?B?r?ﶴt?ﶄ"?STn?"q?3?U?﵏6}?^Ьp?-izݔ?55 ?’ +?ﴙ  ?h8 ?6 ?@Ѭ?ӭ[?ﳢ ??pU"e?>? \{G?}?ﲨ}]?v-%#?De"??z?ﱮ7i?{+?I +MW?R?䌰9?ﰱ&?W ?LZ"!?C"? +R#?ﯴ.W$?﯁<%?N;&?)A'?#(?﮴h[)?ﮁ*?NC5+?,?pыf-?ﭳ(a.?נּ^r/?L B0? Z1?L"@2?בּ|K93?}4?IK5?"6?ᛕ%47?﫭{8?yK>9?E O":?͆;?]kC?>?? +A*!@?Ւme A?褐Ҵ B?lZ_C?7'D?8 E?:®vF?礼,G?c!/H?-kI?얆J?V|K?留]L?XN7M?#dN?킡NoO?醴hHP?廬9NQ?L}?QR?S?T?索@BgU?tV?> +X"W?ƬX?қxY?卵`nWZ?f_[?/n\?PK]?e^?K-0_?U`?6a?Nb?c?zAd?CtIe? z^f?Jjg?Zh?g7Vi?0F^j?Vk?D?l??Flm?RNn?@n'go?KVp?Y;q?tKr?< +us?4t?Ȅ|{u?ҍ}v?\8w?$9x?ݖYy?\z?|{?CWK|? ^y}?%~?f?a ?).nӼ?zo?C{?~䥧j?Fo? 6x.? Ӓf??aظܪ?(?e?A?|9M?CBk? >?10?i?\?#!?R%V?s^m?ut%?;n?w\j?Y.V?+?=?R폒??BX??iYR?.ʌ?0??~'H?C? "?7q?=%e?X2??Td?bx?kl"?0 b?hD?10{?}?B?jJ%?ʲ?H?S%?.?8x߅?2?c)?&Ϣ?•?}נ~?r)d1?5>?QcL?/?:?CQ??"?PjK?Pnٗx?}l?|OT?k?\K-???o]?g,59?)N?>s?uS?q1-?3dӶ?4?ڬ+?y ?<z9?ˁ? ++??C쌖?O??B?ʔ?J{D? C9?y]8q?ϰ?P8r??ҹ@J?!j?TN?ހv?N?j?Xôg??J??[e?Z]?ܛ/ǚ?c?]eŻ??y?v]oy?v?vsr?v1"Xg?u7 ?uR$!?uiIq"?u'D|#?t"k$?t:%?t_:&?tmd'?sّ6Į(?sU$)?sSx{*?sH1+?r̓@,?rYIg-?rGz.?r"9/?q}A0?q}G=1?q9mU2?plU3?p~k4?pn5?p+!Ծ,6?o\?7?o$i8?o_f9?o_w:?nתK;?n?m[??m*@?m>NA?lޱ8B?l^C?lpϯkD?l,0E?k炬CF?kOG?k]BH?kEI?j.^J?j2GK?jJ'gL?j zM?iDžTN?izyO?i5_[P?h@$Q?h%R?he'}eS?hdVAT?g -U?gcV?gN=rW?g?&X?fǻY?f}0~Z?f7=q|[?e:\?e(?2S]?eetH^?e3}_?dؔ};`?dDQ a?dKmb?duԎc?c8d?cxiMe?c1ˬ=f?beg?bbxm%h?b]qi?b}j?ak?aֿl?aA8Xm?`Aln?`Po?`lTpp?`% 8q?_ݶr?_Ophs?_Ngt?_UoYIu?^ +v?^x:V w?^0i"x?]Wy?]Fz?]X[{?] |?\}?\q~?\8?[9l?[gJT?[`?[71?ZQ?ZՊ?Z>J&?YI'r?Y1?YdMuw?YXk?XҭP?Xuy?X@$4U?W?W;L?We`?WZ ?V&9?V8-?V@a?U??Uo5d?Uc?U?K?TЏ?Tc?T=U?S&OU?S9v?S_=?S2?R)D?R?R6\?Qkݕ?Q ?QWv?Q 5?P¯?PxL_?P-u ?Ou?Ou?OM+`ڤ?OIQi?NW?NlW7?N!G.>?M'J?M#U?M?#?Lml-?L.?L]<?L)6?KƟfւ?K{bh?K/\UM?J?J ?JL?J5?I(?Ih#ڣ?I?HT?H?H7v_?G%~K?GI4?GRT?Go ?FF\?Fl>?F&?E?m?Es5?E9Ƃ7?D쮋Q?D +?DRF\'?D>?Cl|<?Ck6c;?C +~?BЛ?B7+<?B5Ö*?A@ .?A?AM So?@\?@}?@cZ?@4<????z0,??+)76?>=S?>gW?>Ayg?=0Ӽ?=؋f?=Vq#D?=??!rc x??!z_z@?!'+LA? q+>B? C? +Ͼ,D?s:E?72F?/츣KG?IR*H?I?3IsEJ?@:K? hL?6}M?!ҤN?zO?9c(P?{IQ?ד\R?;S?N2eT?vWU?zq?a]?၈ +?#, ???fH? ?ߩ*n^?J%#N?0?ލפ?.=E ?]P:?pn?qj +?ܲe{,?SJ?!5#?۔G?54~; ?N +?v-} ?w, ?ٷ ?XfO?n~W?ؙgV?9_q?ٕu?y?!?ֹإ?Yh?Ʃ?ՙ?9t/?4&?x?v?Ӹf 3?W9?v7?Җ ?5"M&!?-"?tjmp#?0$?в]%?Q&?DE'?ϏbE(?.:#)?;*?lH1+? Q,?ͩP-?HJ.?G/?̅Y#W0?#op1?-,r2?`[3@3?Z4?ʜ5?;&M6?@J7?wL8?H9?ȳ7=:?Q 0;?SK? +??e@?$RA?Š.KB?> +moC?g~sD?x E?MF?ó'oG?PJ + +H?_!.I?Še""J?'\bK?ERL?a #SM?%%N?O?7XbXP?6Q?pB +R? XS?|T?E2|U?BV?~FiW?sX?':2Y?ScqZ?"(O[?#vh4\?'NRS]?^?^Н_?U`?Qxa?1b?̮͗c?i%yd?.e?Õf?;wgg?ćMh?r?i? HLj?p|k?Cl?ޖWm?y9n?&lo?cxp?J58q?m;]r?l9s?Ut?a1u?Ox Tv?a#w?` x?(Vy?=z?SRLtT{?|?@}?!%a~?ʨ?U>q?ﰊF?A?#\?.?V?sD? <5f?#f4? 5G?Vs?Ͽ?}jy?"\?z?TX?/F?j?¥??<?Qa?SS?`?_(?č?L +?|x?}Y?[7?P?G?>?wmu???,??pi?כi?o?\w;?'̑?6^.l?Β2?f3 `???HĚB?,,xP?%Nu|?[~+-?<?>?!3'?R'?Od?g?}\߃?Cd#?%"?An?آyv?oPFv?֔?)m?3@a?y?_߼?8"?O?"CHw?Q?O ?:?{?ë?n?9I?94?hp?۷?6.?'fH??Q\ǻ?0?zA?h̐6?0g?8o+3?^x?a]})?MK??@Y?oX?F5S?ڤ[ߦ?nb $?H[?T?*c?0Eg?Q /?~G?~yA5?xX1ɦ?x' ?w❖?w>"?v?O;B?vc틸?u^|?u\?uEp"%?tթ ?t?. +?s6{ ?sdj' ?r ?r ?rD^?q_!?q>q?p}|}?pbO?o~?ov\?oqj8y?n YC?n:D&?m,?m]?l!?l:9|?l?ke?k3F?jhZ ?jV Oz!?i 9wM"?iw&#?i"$?h)%?h* n&?g: +'?gK(?f܄!)?fm9*?eX+?e +x(,?ep!b-?dF.?d?9`/?cOk0?c_}1?b[P2?b3?b`x4?aF5?a/6?`oZG7?`O<#8?_&9?_n:?^N;?^{?]8{B?[yq'C?[TD?Z#WE?Z&~F?Y+G?YDԛH?Xӄ"TI?XbV+J?W(K?WѺ]L?Wz~}M?VbN?V+fuO?U!@P?UHУQ?T8TjR?TeK¢S?Sp5T?SAU?S7wV?RS"W?R,'_yX?Q)+Y?QH1/Z?PB[?Pc!\?O*-]?Of\^?O _?NAa`?N(J=a?Mֈb?MCHBc?Lгd?L^ޭe?K^#Df?KxJg?K嚎h?J@i?J l^j?Ik?I:l?Hm?HS3n?G=do?Gmxҩup?F/dq?FءAr?FsTs?EO }t?E,Iu?DVv?DEWw?Cѭ]x?C]"y?B0z?Bv]{?B}(J|?A}?Am~?@?@2pۮ??Kip??J9c?>K?>a?=,:\?=x%?=J;C?<ģ?<1SM%?;I?;1K)?:%F?:HZ{T?9ӂ|?9^^A?8-?8t +?7Rsv?7}(Z?7SϬ?6}?6*?5Q76?5@"sk?4ʳv?4U74?3߭ĹT?3j0Pe?2p ?2~K?2n ?1/-?1SDu?0i{]?01rx?/m?/E[Z?.:⊿?.Y  ?-ј?-l~?,1 ?,_?, []uC?+ۺ9U?+Nv?*H?*/ ?)TD?)A1?(ʿ֧F?(S?'K?'eD?&U?&w?&?%f!?%K;N?$zB?$##a?#?6<?#4ʵQO?"HC?"EJ?!?!Vp,? ޸q? f? ?w>k[?Oun?SJ?I{?2N2? ~t?"?.=?Mͻ?=d?ŊT?ML|b?ԐN?[?a ?j3`?IIO?y4!t?_wζ?}K +?y?o?(?mC?*G?_5 ??k闒?/;k?xhDX? j? Ѫ? +X? ">? y? ^? "{3 ? +HL2? +.h? ]? 9_B(?`?Dd?i?Ok{?m)?Z!]?h^r?dv~?ͬU?n7?\u?y ?tx?i3??'?Jƴ?\j?db* ?= GA +?#T ?F= ?,ԙ ?Um?5o!?;l?=T?!2?FAw?T'?NZ5?RT&?V=P?A7>?]똻E?ᮛ>?edIxL? L?lh?5aK?sŽ?(8t ?zh!?߅"?26y#?p7X$?򇠘J%? +&?kW'?p(?)?&*?b+?/>,?Fj6 -?!^.? @/?'Iv0?ڛB%1?,^|+2?3?1>uf4?곚(5?5gp6?+7?:_Y8?輆sa9?>P5,:?;?BQf:?[1??J%̉@?A?MΐB?5UC?P4wD?SRE?S*VF?L-`?U6a?h`e>hb?]c?hM6[d?1e?hf?%g?gJh?:j{i?f!+j?pmk?eRl?rm?d1rn??Cbo?c>Gp?#,q?b;r?<s?`[U3t?lq(u?^qk\v?h۲w?\S;/x?0oy?Zʑz?k{?Wzi|?#5a}?T?~?N=?Q1.?E/?N?Я?KU?ɖgB@?G(l?}?D +7?LT7?@? Z?;Pl??7 ?RH?3^?T(?.@I\??)HE?fZ?$$w%?sD???*aL?ME?c3?l,? h03?W@?9]g?z?k?} t?@aS?vʖ?vC?o{?xn6%?h ?F-I?a8?}?ZZ?H'?Ri?|d.?J5?|R ?Bh F?H_ ?:K ?z?1?F43?(;?xL?3?vY?„?A? j??m?<?[?umZ?r!?kjβX?V*t?a4*?c?VBxE?уN`?L.&??A^פ?휻Ro?6[L C?훰wq?+$`?횥vf?+ ?홙&R?UR?혎<?NX?헂S-{?K;?v60??i?㫷?]cd?]#?P !??m?C(6w?푽=?6D?퐰5?)Zh؂?폢c?ؚ?펕*?#ޗ?퍇6jW?<?!y]?kg'N?^ $?\m?T-?N?ơ//m??4G?퇷i F?05|W?톨Gy?!?텙W0??턉٥?D ?z(?=dJ ?jE] ?Aw ?Z/O ?44 ?I* ?># ?9l ?~y ?~(LI ?}T} +?} FY ?|ZI ?|G ?{~. ?zx ?zlw*+ ?yK ?y[? ?xd p ?xI|" ?w ?w7Y ?vx ?v%^.~ ?u76`a ?u ?tÄs ?tvGd ?swV ?rf ?rdE0AF ?q- ?qQ:pj !?pǢwT "?p=Z #?oL $?o* %?n[ &?nt0 '?m + (?m )?ly *?kq +?kdڷ ,?jܼ| -?jPQ .?irKN /?i<+= 0?h +X 1?h'vSk 2?g ( 3?g^9 4?f "> 5?eveM 6?er) 7?d+n3 8?d]s6XD 9?cҮ :?cGN ;?b" ?a? ??` 'x @?` A?_zǙ B?^ C?^dQ4C D?]R` E?]MӐ F?\B4g G?\6. H?[ON I?[1_ J?Z*T K?Z L?Y|AI M?XkQ N?XeL(T O?Wuy P?WM_ Q?V R?V5 S?U T?U U?Th/ V?T;+3 W?Sy:Ǣ X?R嶄 Y?R`h,2 Z?Q Ф [?QG \?P% ]?P.rt ^?O+M _?Os `?Nu a?M" b?MoRU c?L.A d?LU_V{ e?K6 f?K;γS g?JB h?J!¤cT i?Iϫ j?I8/ k?HzTl l?Gw5 m?G_̘ + n?Fve< o?FE +O p?E +? q?E*)] r?Dl s?D t?Cn8 u?B0* v?Bf V w?AD$ x?AJs@h y?@B z?@.]W {?? |??I }?>G^ ~?= ?=hbae ?<. ?@` ? +՜R ? +^a ? ŋ ?Mk ?a ?  ? +#?*ps +$?BS +%?ެVI +&?Q` +'?݀ +(? +)?Sހ +*?۽'W ++?'G] +,?ڑ=p +-? +.?dÀ +/?0 +0?7O: +1?סBfo +2? +M +3?t%~s +4?݅ +5?F +6?԰"t +7?_i[ +8?ӂM +9?J +:?T9 +;?ѽ|z4 +?MV +??aD +@?vl. +A?3?9o +B?͛_ +C? +D?mRX+o +E? +F?>z>V% +G?ʦ +H?s + +I?w4R +J?> +K?Hc +L?ǰ +M? +N?ƁHo +O?n +P?Qtp +Q?ĹE +R?!L^ +S?É9 +T?|H +U?Y\S +V?1+1` +W?( +X?e9 +Y?i(? +Z?`Y +[?ǩ +\?/9JC +]?콖@ +^?44 +_?e> +`? +a?4W' +b?캛  +c?1 +d?j +e?9 +f?8T8 +g?췟d! +h?i +i?mb +j?P? +k?;2V(^ +l?촢~ +m?x* +n?o +o?G2 +p? +s?pN: +t?3 +u?=M +v?쮣w +w? +J- +x?pe +y? +z?= +{?쫣99 +|? X*E +}?okd +~?s? +?;pP +?쨡a{ +?G6ȣ +?m!` +?b&c +?8 +?쥞k٬ +?t$ +?i7 +?Oj +?4ȣ +?좚X +?M$ +?e4v} +?ʑ: +?/i +?쟕( +?c/ +?_fm2 +?Ķ +?t:^8 +?s;Q5 +?r +?r`܈rJ +?qÉ +?q$Kj +?p{E +?o6 +?oJL| +?nN9 +?n A +?mopC +?l +?l2} +?k0ƂP +?jE +?jW$ +?iy9 +?i + +?h{=1 +?g܅D +?g=< +?f +?f_+= +?ea5c +?dF=t2 +?d#K= +?cEޠ +?b5~ +?bFIL +?am, +?a +?`hM +?_;T +?_)[ +?^D +?]!| + +?]Ke: +?\.J +?\ +Ӷ +?[m +?Zn ?Z- ?Y ^ ?XLvm ?XNRp ?W1 ?W ?Vn ?Uˌ ?U.է ?T<^ +?Sk) ?SN ?Rd ?R+h ?Qm&H ?P͖YB ?P-;k ?O嶳 ?NeA ?NK ?Mc(W ?M +ѵ ?Lj5T ?KɎ_ ?K(ͻ ?J ?IV ?IF} ?H ?H ?GcH ?Fʳ ?F!9 !?E "?DߍG #?D>cw $?C.r8 %?Bs &?BZl '?AN$ (?AwW )?@vC *?? J] +??3|d ,?>cP -?=h9 .?=NƗݺ /?<8J 0?< b 1?;i 2?:Ӭh 3?:%%r 4?9: 5?8+Q 6?8@3. 7?70i1 8?6"W 9?6Z :?5K` ;?5U ?3.G=" ??2̖ @?1.s- A?1G~ B?0EE C?0 D?/`1B E?. F?.N G?-xCB H?,ՈΈ I?,2Âmp J?+| K?*= L?*J3 M?)Cr=Q N?)H O?(aB P?'2!8 Q?'Lg R?&w S?%dP= T?%12B U?$@K V?# W?#G` X?".]A Y?"@ Z?!]BBH [? +l \? +ec| ]?r< ^?Α _?+8> `?} a?^ b?? c?  d?$h e?T4, c f?8 g? 2.b h?h!˃ i? j? k?{ l?t m?3.= n?F o?% p?Fgd q? l` r?26 s?X[ t?Y u?6 v?j8s w?,D x?!q(] y? | z? c {? 35o |? ]k }? +.[ ~? +D4_L ? / ? j ?U ?} ? +n ?ez ?6 ? ?u ?,*+ ?*/B ?E ?ª+E ?:5 ?f ?P ?IMs ? ? ?XA3/ ?03| ? NϹ ?fb ?l 3J ?k$ ?t_a ?I ?()[R' ?K ?-q ?5gI ??, ?3= ?B_: ?!5 ?& ?O/]Y ?. ?2~ ?[uz` ?U ? ?g]f ?j ?Ġ8" ?rn ? ?% ?~xD ?0| ?/ ?y ? ?:ֹ ?W$. ? ?DĹsXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 16 / length of dimension 1 NAXIS2 = 3000 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 2 / number of table fields TTYPE1 = 'ell ' TFORM1 = 'K ' TTYPE2 = 'beam ' TFORM2 = 'D ' SACCTYPE= 'tracer ' SACCCLSS= 'Map ' SACCNAME= 'ck_so ' SACCQTTY= 'cmb_convergence' EXTNAME = 'tracer:Map:ck_so:beam' MAP_UNIT= 'none ' SPIN = 0 END ??,?چ?>?\-? +D??5L:?>g& ? +?z` ?݋ ?2' ?-:d?G#?R%?L ?5[??g[?ð>?7?V?VS?bl?3d??ϗ?Қ?)yl?>]w?B% ?6N +!?"?#?ׯ[$?aM %?} &?Д8;'?ܠ(?ˆ)?,3*?7$W+?vV/,?-?=.?6/?Gp0?R1?E2?i1.3?%J 4?j15?m>6?7?t8?+59?9Y:?;?k0???~@?zAA?v|sB?rVC?n:GD?iE?ev"F?`^G?\p H?WAI?S)DJ?NmcK?I +.L?D xM??ֺ2[N?: +@O?5P?0xQ?+}R?&>]S? T?0U?hV?3;W? X?klY?bZ?a[?$\?A]?N^?Jl_?6T`? k?Hjgl?_um?f!n?y\nkgo?rB[hp?kpq?cr?\s?U6Uft?Melu?FNv?>hy\w?7$[-x?/vy?'#bEz?{? o|?_}?@F~?_??2\?}pY'?,Om?ЍW?X?ֶ?D1#?#LO?s?*?V/?q<?|]h?wt?vb?m<D?dd?Z0B?Qh-G?HMN?> C?5u?+hw?!y?k?<?b"?x?}1?rS?Wd?+uu?杖?/:?FX> +?#?[?ͥ6i?/Z5?~>?sy?hQ?^z?S"}?H" ?c?=8[?1 b?&B?|y?*_f?l?T˼?xh{?=Ɂl?֙p?YBf? H?K{&?f+?p0?j?SzF?w,-?jWm<?^'?RUɢ!?Ez?9tG*?,Sk? Rn ???FZ?$6?In?^M?bC?V?:,?E?c1?'Q?&j?v1?i:[U?[ƿ?N H?@]k?2$?$͹?[W?F??,;|?kl ?ЖR3?S'?b??*z ??z~N?kmI?\#܎?Mu??j3?08$/?!BR?;n?$!?&O? )?~"C?%?HF?DY?^l?!wH?xwW?h<"c?XA?INW?9,]_?)0?$? y??^?QV?zC?,?M?y W??u*?dk+?S6?B ?1P< +? ^@ ? ?~U ?dH?"N?Ϲ+?lк??vt?\?s=C?a]_a?OĀ?=O ?, >?? ΅?S?І.?љf5?Q?0 ?!?"?ud#?b3v$?PNz%?=UZ&?* '?<(?)? +_*?D+?&%,?B)-?| u.?/Ҍ/?~Ұ +0?keu1?W&.2?DZM3?07ڹ4? 5? O6?7?=8?ͱ 9?.:?;?\V?T}??@[@?+,pA?fA&B?SqC?0тD?}E?ĹF?G?ިH? 7|I?qc3J?[CnK?FL?13M?RtN?yO?"P?6Q?ƖR?y_S?W +/T?񅟎dU?o'%V?YtlW?Dv]cX?.-Y?Z?u[?ΏZJ\?Քa]?I[AMϠ?##?%? p?)Q?;N&?<~?-v?gQnJ?K?0?O[5?u<?~j?U?lo?ˋ?p?TXld ?8y?W8??M]?Ȝm?yP) ?E ?t`'?W?;J?7Y?Qc???b?揝K?K?rǞ'V?Up ?8E?7Լ?ΡI??pƗ)?)Z?E?ljO$?N_A?1j7K?<?)<?pj?产`x?#?~䤪?`&?B +. ?$?J?b^??f?㫼C?Q?nթ?PJ g?17w?.u?E|?y{>?⶜p?◯?x?YJ?:{?Zm?O?qݯ?p&?D?~?^)??T ??1??69?F X?FI+Q?`5?@+?R>?Hz2?R ?߾?ߞe?};U?]k@mr?<!?B?T09?ځv?޹**?ޘvc?w1) ?V?5q ?=J0 +?L ?ѥ 6 ?ݰ@! ?ݎ@?mG=U?K ?* @z?Xg)?b?ĽpH?ܢ׳j?܀V?^۶/y? .?Y=;/?0?C1?ذZvn82?،s]3?hS#~4?Dʉ}5? ϛ`W6?Ĉ[ 7?ةP8?״}K9?אBuSQ:?kѶ;?G +~?&`<??ֵ@?֐pA?l!.B?GTKKC?"wD?<E?؍F?ճdG?Վc*'H?i5I?C0J?B,K?L_LL?\M?Ԯ`9N?ԈO?c3RXP?=*jQ?wER?>S?-T?Ӧ)mU?Ӏ*V?Z[W?3EX? b.Y?HZ?=ko[?Қ٠\?tm]?M^?']"_?i\`? Y\a?ѳLVb?ь|ggc?e0vd?> +|e?48f?E(Ug?x;˚h?ТGi?{&j?Sqk?,S(l?m?_Q~n?ϵ>^o?ώ,pNp?fzTq?> r?xs?M\t? |pu?Ο:v?vDձw?Nx?&&y?ew4z?%{?ͭ%|?ͅ<}?\M~?4&o? ?h?̺ ?̑9M?hV8dH??bX2?_ɭ?K_q?' R?˚ʃ8?q_&?H[7?~?@?x?ʢkmfE?x$H?O`T?%NbE?y,?ѕdaO?ɧk?}*?Sx:?)b7?-p? ?Ȫ$?Ȁ/.~h?U+<?+5+@?xv??ǫEE?ǀfH?U +?*K?f?Z?Ʃ1?~?SZq?(R?S?YR?ť%?z^N_?NȐ?#"˻I?l+?˧*b?ğO?smv?G„??٠U?ó?×}F?k7<?>j?{W7?q`}?¹Y2,?Œ=P?`C? +?3?ӟA?Y}? +~u?P8?S?%š?Ѻ5?˛&?UO3?p?C_??"?rƸ?;ޖ?ￍ` +?_2?1? i3?'8#?ﾨ1 M?z*B?L?i?f(?rE'd?ｓ'?dW)?6??MS?#DH?Ｊ}B?{Ӎ?Mn?+,?E?P]?ﻑJ$^?b4H\:?3s?ئi?Ԓ?ﺥ=#?und?Fa¯?o?Fm?ﹷ?﹇nD]?X%?(Pz?k?uVg?︘pi?h[$"?85?*?׻?﷧f?w}?F?B?r?ﶴt?ﶄ"?STn?"q?3?U?﵏6}?^Ьp?-izݔ?55 ?’ +?ﴙ  ?h8 ?6 ?@Ѭ?ӭ[?ﳢ ??pU"e?>? \{G?}?ﲨ}]?v-%#?De"??z?ﱮ7i?{+?I +MW?R?䌰9?ﰱ&?W ?LZ"!?C"? +R#?ﯴ.W$?﯁<%?N;&?)A'?#(?﮴h[)?ﮁ*?NC5+?,?pыf-?ﭳ(a.?נּ^r/?L B0? Z1?L"@2?בּ|K93?}4?IK5?"6?ᛕ%47?﫭{8?yK>9?E O":?͆;?]kC?>?? +A*!@?Ւme A?褐Ҵ B?lZ_C?7'D?8 E?:®vF?礼,G?c!/H?-kI?얆J?V|K?留]L?XN7M?#dN?킡NoO?醴hHP?廬9NQ?L}?QR?S?T?索@BgU?tV?> +X"W?ƬX?қxY?卵`nWZ?f_[?/n\?PK]?e^?K-0_?U`?6a?Nb?c?zAd?CtIe? z^f?Jjg?Zh?g7Vi?0F^j?Vk?D?l??Flm?RNn?@n'go?KVp?Y;q?tKr?< +us?4t?Ȅ|{u?ҍ}v?\8w?$9x?ݖYy?\z?|{?CWK|? ^y}?%~?f?a ?).nӼ?zo?C{?~䥧j?Fo? 6x.? Ӓf??aظܪ?(?e?A?|9M?CBk? >?10?i?\?#!?R%V?s^m?ut%?;n?w\j?Y.V?+?=?R폒??BX??iYR?.ʌ?0??~'H?C? "?7q?=%e?X2??Td?bx?kl"?0 b?hD?10{?}?B?jJ%?ʲ?H?S%?.?8x߅?2?c)?&Ϣ?•?}נ~?r)d1?5>?QcL?/?:?CQ??"?PjK?Pnٗx?}l?|OT?k?\K-???o]?g,59?)N?>s?uS?q1-?3dӶ?4?ڬ+?y ?<z9?ˁ? ++??C쌖?O??B?ʔ?J{D? C9?y]8q?ϰ?P8r??ҹ@J?!j?TN?ހv?N?j?Xôg??J??[e?Z]?ܛ/ǚ?c?]eŻ??y?v]oy?v?vsr?v1"Xg?u7 ?uR$!?uiIq"?u'D|#?t"k$?t:%?t_:&?tmd'?sّ6Į(?sU$)?sSx{*?sH1+?r̓@,?rYIg-?rGz.?r"9/?q}A0?q}G=1?q9mU2?plU3?p~k4?pn5?p+!Ծ,6?o\?7?o$i8?o_f9?o_w:?nתK;?n?m[??m*@?m>NA?lޱ8B?l^C?lpϯkD?l,0E?k炬CF?kOG?k]BH?kEI?j.^J?j2GK?jJ'gL?j zM?iDžTN?izyO?i5_[P?h@$Q?h%R?he'}eS?hdVAT?g -U?gcV?gN=rW?g?&X?fǻY?f}0~Z?f7=q|[?e:\?e(?2S]?eetH^?e3}_?dؔ};`?dDQ a?dKmb?duԎc?c8d?cxiMe?c1ˬ=f?beg?bbxm%h?b]qi?b}j?ak?aֿl?aA8Xm?`Aln?`Po?`lTpp?`% 8q?_ݶr?_Ophs?_Ngt?_UoYIu?^ +v?^x:V w?^0i"x?]Wy?]Fz?]X[{?] |?\}?\q~?\8?[9l?[gJT?[`?[71?ZQ?ZՊ?Z>J&?YI'r?Y1?YdMuw?YXk?XҭP?Xuy?X@$4U?W?W;L?We`?WZ ?V&9?V8-?V@a?U??Uo5d?Uc?U?K?TЏ?Tc?T=U?S&OU?S9v?S_=?S2?R)D?R?R6\?Qkݕ?Q ?QWv?Q 5?P¯?PxL_?P-u ?Ou?Ou?OM+`ڤ?OIQi?NW?NlW7?N!G.>?M'J?M#U?M?#?Lml-?L.?L]<?L)6?KƟfւ?K{bh?K/\UM?J?J ?JL?J5?I(?Ih#ڣ?I?HT?H?H7v_?G%~K?GI4?GRT?Go ?FF\?Fl>?F&?E?m?Es5?E9Ƃ7?D쮋Q?D +?DRF\'?D>?Cl|<?Ck6c;?C +~?BЛ?B7+<?B5Ö*?A@ .?A?AM So?@\?@}?@cZ?@4<????z0,??+)76?>=S?>gW?>Ayg?=0Ӽ?=؋f?=Vq#D?=??!rc x??!z_z@?!'+LA? q+>B? C? +Ͼ,D?s:E?72F?/츣KG?IR*H?I?3IsEJ?@:K? hL?6}M?!ҤN?zO?9c(P?{IQ?ד\R?;S?N2eT?vWU?zq?a]?၈ +?#, ???fH? ?ߩ*n^?J%#N?0?ލפ?.=E ?]P:?pn?qj +?ܲe{,?SJ?!5#?۔G?54~; ?N +?v-} ?w, ?ٷ ?XfO?n~W?ؙgV?9_q?ٕu?y?!?ֹإ?Yh?Ʃ?ՙ?9t/?4&?x?v?Ӹf 3?W9?v7?Җ ?5"M&!?-"?tjmp#?0$?в]%?Q&?DE'?ϏbE(?.:#)?;*?lH1+? Q,?ͩP-?HJ.?G/?̅Y#W0?#op1?-,r2?`[3@3?Z4?ʜ5?;&M6?@J7?wL8?H9?ȳ7=:?Q 0;?SK? +??e@?$RA?Š.KB?> +moC?g~sD?x E?MF?ó'oG?PJ + +H?_!.I?Še""J?'\bK?ERL?a #SM?%%N?O?7XbXP?6Q?pB +R? XS?|T?E2|U?BV?~FiW?sX?':2Y?ScqZ?"(O[?#vh4\?'NRS]?^?^Н_?U`?Qxa?1b?̮͗c?i%yd?.e?Õf?;wgg?ćMh?r?i? HLj?p|k?Cl?ޖWm?y9n?&lo?cxp?J58q?m;]r?l9s?Ut?a1u?Ox Tv?a#w?` x?(Vy?=z?SRLtT{?|?@}?!%a~?ʨ?U>q?ﰊF?A?#\?.?V?sD? <5f?#f4? 5G?Vs?Ͽ?}jy?"\?z?TX?/F?j?¥??<?Qa?SS?`?_(?č?L +?|x?}Y?[7?P?G?>?wmu???,??pi?כi?o?\w;?'̑?6^.l?Β2?f3 `???HĚB?,,xP?%Nu|?[~+-?<?>?!3'?R'?Od?g?}\߃?Cd#?%"?An?آyv?oPFv?֔?)m?3@a?y?_߼?8"?O?"CHw?Q?O ?:?{?ë?n?9I?94?hp?۷?6.?'fH??Q\ǻ?0?zA?h̐6?0g?8o+3?^x?a]})?MK??@Y?oX?F5S?ڤ[ߦ?nb $?H[?T?*c?0Eg?Q /?~G?~yA5?xX1ɦ?x' ?w❖?w>"?v?O;B?vc틸?u^|?u\?uEp"%?tթ ?t?. +?s6{ ?sdj' ?r ?r ?rD^?q_!?q>q?p}|}?pbO?o~?ov\?oqj8y?n YC?n:D&?m,?m]?l!?l:9|?l?ke?k3F?jhZ ?jV Oz!?i 9wM"?iw&#?i"$?h)%?h* n&?g: +'?gK(?f܄!)?fm9*?eX+?e +x(,?ep!b-?dF.?d?9`/?cOk0?c_}1?b[P2?b3?b`x4?aF5?a/6?`oZG7?`O<#8?_&9?_n:?^N;?^{?]8{B?[yq'C?[TD?Z#WE?Z&~F?Y+G?YDԛH?Xӄ"TI?XbV+J?W(K?WѺ]L?Wz~}M?VbN?V+fuO?U!@P?UHУQ?T8TjR?TeK¢S?Sp5T?SAU?S7wV?RS"W?R,'_yX?Q)+Y?QH1/Z?PB[?Pc!\?O*-]?Of\^?O _?NAa`?N(J=a?Mֈb?MCHBc?Lгd?L^ޭe?K^#Df?KxJg?K嚎h?J@i?J l^j?Ik?I:l?Hm?HS3n?G=do?Gmxҩup?F/dq?FءAr?FsTs?EO }t?E,Iu?DVv?DEWw?Cѭ]x?C]"y?B0z?Bv]{?B}(J|?A}?Am~?@?@2pۮ??Kip??J9c?>K?>a?=,:\?=x%?=J;C?<ģ?<1SM%?;I?;1K)?:%F?:HZ{T?9ӂ|?9^^A?8-?8t +?7Rsv?7}(Z?7SϬ?6}?6*?5Q76?5@"sk?4ʳv?4U74?3߭ĹT?3j0Pe?2p ?2~K?2n ?1/-?1SDu?0i{]?01rx?/m?/E[Z?.:⊿?.Y  ?-ј?-l~?,1 ?,_?, []uC?+ۺ9U?+Nv?*H?*/ ?)TD?)A1?(ʿ֧F?(S?'K?'eD?&U?&w?&?%f!?%K;N?$zB?$##a?#?6<?#4ʵQO?"HC?"EJ?!?!Vp,? ޸q? f? ?w>k[?Oun?SJ?I{?2N2? ~t?"?.=?Mͻ?=d?ŊT?ML|b?ԐN?[?a ?j3`?IIO?y4!t?_wζ?}K +?y?o?(?mC?*G?_5 ??k闒?/;k?xhDX? j? Ѫ? +X? ">? y? ^? "{3 ? +HL2? +.h? ]? 9_B(?`?Dd?i?Ok{?m)?Z!]?h^r?dv~?ͬU?n7?\u?y ?tx?i3??'?Jƴ?\j?db* ?= GA +?#T ?F= ?,ԙ ?Um?5o!?;l?=T?!2?FAw?T'?NZ5?RT&?V=P?A7>?]똻E?ᮛ>?edIxL? L?lh?5aK?sŽ?(8t ?zh!?߅"?26y#?p7X$?򇠘J%? +&?kW'?p(?)?&*?b+?/>,?Fj6 -?!^.? @/?'Iv0?ڛB%1?,^|+2?3?1>uf4?곚(5?5gp6?+7?:_Y8?輆sa9?>P5,:?;?BQf:?[1??J%̉@?A?MΐB?5UC?P4wD?SRE?S*VF?L-`?U6a?h`e>hb?]c?hM6[d?1e?hf?%g?gJh?:j{i?f!+j?pmk?eRl?rm?d1rn??Cbo?c>Gp?#,q?b;r?<s?`[U3t?lq(u?^qk\v?h۲w?\S;/x?0oy?Zʑz?k{?Wzi|?#5a}?T?~?N=?Q1.?E/?N?Я?KU?ɖgB@?G(l?}?D +7?LT7?@? Z?;Pl??7 ?RH?3^?T(?.@I\??)HE?fZ?$$w%?sD???*aL?ME?c3?l,? h03?W@?9]g?z?k?} t?@aS?vʖ?vC?o{?xn6%?h ?F-I?a8?}?ZZ?H'?Ri?|d.?J5?|R ?Bh F?H_ ?:K ?z?1?F43?(;?xL?3?vY?„?A? j??m?<?[?umZ?r!?kjβX?V*t?a4*?c?VBxE?уN`?L.&??A^פ?휻Ro?6[L C?훰wq?+$`?횥vf?+ ?홙&R?UR?혎<?NX?헂S-{?K;?v60??i?㫷?]cd?]#?P !??m?C(6w?푽=?6D?퐰5?)Zh؂?폢c?ؚ?펕*?#ޗ?퍇6jW?<?!y]?kg'N?^ $?\m?T-?N?ơ//m??4G?퇷i F?05|W?톨Gy?!?텙W0??턉٥?D ?z(?=dJ ?jE] ?Aw ?Z/O ?44 ?I* ?># ?9l ?~y ?~(LI ?}T} +?} FY ?|ZI ?|G ?{~. ?zx ?zlw*+ ?yK ?y[? ?xd p ?xI|" ?w ?w7Y ?vx ?v%^.~ ?u76`a ?u ?tÄs ?tvGd ?swV ?rf ?rdE0AF ?q- ?qQ:pj !?pǢwT "?p=Z #?oL $?o* %?n[ &?nt0 '?m + (?m )?ly *?kq +?kdڷ ,?jܼ| -?jPQ .?irKN /?i<+= 0?h +X 1?h'vSk 2?g ( 3?g^9 4?f "> 5?eveM 6?er) 7?d+n3 8?d]s6XD 9?cҮ :?cGN ;?b" ?a? ??` 'x @?` A?_zǙ B?^ C?^dQ4C D?]R` E?]MӐ F?\B4g G?\6. H?[ON I?[1_ J?Z*T K?Z L?Y|AI M?XkQ N?XeL(T O?Wuy P?WM_ Q?V R?V5 S?U T?U U?Th/ V?T;+3 W?Sy:Ǣ X?R嶄 Y?R`h,2 Z?Q Ф [?QG \?P% ]?P.rt ^?O+M _?Os `?Nu a?M" b?MoRU c?L.A d?LU_V{ e?K6 f?K;γS g?JB h?J!¤cT i?Iϫ j?I8/ k?HzTl l?Gw5 m?G_̘ + n?Fve< o?FE +O p?E +? q?E*)] r?Dl s?D t?Cn8 u?B0* v?Bf V w?AD$ x?AJs@h y?@B z?@.]W {?? |??I }?>G^ ~?= ?=hbae ?<. ?@` ? +՜R ? +^a ? ŋ ?Mk ?a ?  ? +#?*ps +$?BS +%?ެVI +&?Q` +'?݀ +(? +)?Sހ +*?۽'W ++?'G] +,?ڑ=p +-? +.?dÀ +/?0 +0?7O: +1?סBfo +2? +M +3?t%~s +4?݅ +5?F +6?԰"t +7?_i[ +8?ӂM +9?J +:?T9 +;?ѽ|z4 +?MV +??aD +@?vl. +A?3?9o +B?͛_ +C? +D?mRX+o +E? +F?>z>V% +G?ʦ +H?s + +I?w4R +J?> +K?Hc +L?ǰ +M? +N?ƁHo +O?n +P?Qtp +Q?ĹE +R?!L^ +S?É9 +T?|H +U?Y\S +V?1+1` +W?( +X?e9 +Y?i(? +Z?`Y +[?ǩ +\?/9JC +]?콖@ +^?44 +_?e> +`? +a?4W' +b?캛  +c?1 +d?j +e?9 +f?8T8 +g?췟d! +h?i +i?mb +j?P? +k?;2V(^ +l?촢~ +m?x* +n?o +o?G2 +p? +s?pN: +t?3 +u?=M +v?쮣w +w? +J- +x?pe +y? +z?= +{?쫣99 +|? X*E +}?okd +~?s? +?;pP +?쨡a{ +?G6ȣ +?m!` +?b&c +?8 +?쥞k٬ +?t$ +?i7 +?Oj +?4ȣ +?좚X +?M$ +?e4v} +?ʑ: +?/i +?쟕( +?c/ +?_fm2 +?Ķ +?t:^8 +?s;Q5 +?r +?r`܈rJ +?qÉ +?q$Kj +?p{E +?o6 +?oJL| +?nN9 +?n A +?mopC +?l +?l2} +?k0ƂP +?jE +?jW$ +?iy9 +?i + +?h{=1 +?g܅D +?g=< +?f +?f_+= +?ea5c +?dF=t2 +?d#K= +?cEޠ +?b5~ +?bFIL +?am, +?a +?`hM +?_;T +?_)[ +?^D +?]!| + +?]Ke: +?\.J +?\ +Ӷ +?[m +?Zn ?Z- ?Y ^ ?XLvm ?XNRp ?W1 ?W ?Vn ?Uˌ ?U.է ?T<^ +?Sk) ?SN ?Rd ?R+h ?Qm&H ?P͖YB ?P-;k ?O嶳 ?NeA ?NK ?Mc(W ?M +ѵ ?Lj5T ?KɎ_ ?K(ͻ ?J ?IV ?IF} ?H ?H ?GcH ?Fʳ ?F!9 !?E "?DߍG #?D>cw $?C.r8 %?Bs &?BZl '?AN$ (?AwW )?@vC *?? J] +??3|d ,?>cP -?=h9 .?=NƗݺ /?<8J 0?< b 1?;i 2?:Ӭh 3?:%%r 4?9: 5?8+Q 6?8@3. 7?70i1 8?6"W 9?6Z :?5K` ;?5U ?3.G=" ??2̖ @?1.s- A?1G~ B?0EE C?0 D?/`1B E?. F?.N G?-xCB H?,ՈΈ I?,2Âmp J?+| K?*= L?*J3 M?)Cr=Q N?)H O?(aB P?'2!8 Q?'Lg R?&w S?%dP= T?%12B U?$@K V?# W?#G` X?".]A Y?"@ Z?!]BBH [? +l \? +ec| ]?r< ^?Α _?+8> `?} a?^ b?? c?  d?$h e?T4, c f?8 g? 2.b h?h!˃ i? j? k?{ l?t m?3.= n?F o?% p?Fgd q? l` r?26 s?X[ t?Y u?6 v?j8s w?,D x?!q(] y? | z? c {? 35o |? ]k }? +.[ ~? +D4_L ? / ? j ?U ?} ? +n ?ez ?6 ? ?u ?,*+ ?*/B ?E ?ª+E ?:5 ?f ?P ?IMs ? ? ?XA3/ ?03| ? NϹ ?fb ?l 3J ?k$ ?t_a ?I ?()[R' ?K ?-q ?5gI ??, ?3= ?B_: ?!5 ?& ?O/]Y ?. ?2~ ?[uz` ?U ? ?g]f ?j ?Ġ8" ?rn ? ?% ?~xD ?0| ?/ ?y ? ?:ֹ ?W$. ? ?DĹsXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 16 / length of dimension 1 NAXIS2 = 401 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 2 / number of table fields TTYPE1 = 'z ' TFORM1 = 'D ' TTYPE2 = 'nz ' TFORM2 = 'D ' SACCTYPE= 'tracer ' SACCCLSS= 'NZ ' SACCNAME= 'gc_unwise' SACCQTTY= 'galaxy_density' EXTNAME = 'tracer:NZ:gc_unwise' HIERARCH META_ngal = 1.0 END ?zG{? >?zG{?)U}?Q??+M?Gz?ҧ6?ۅQ?ʜğ?(\)?b??-z?p +=q?#X1?zG? e-?޸Q??~v?\(\???ho +?QR?tڿ? +=p?()?\(?)22?GzH?? .?ᙙ?V@?Q?&?=p +>??s?\(?غ.t?Gz?.pu?333333?J )?Q?Pe? +=p?w]퉚?(\)? B?zG{?2-(??;s T?Q?f4S?p +=q?gmOn?\(?8E?zG?(A?fffffg?%U!ӎ?Q?&|? +=p +?&6?\(\?oG?zG? ??T?QR?I)? +=p?6?\(?g!0Q?GzH? ?`?\(?9<:?Gz?~ϯ?333333? +?Q?(S_? +=p?Kp^b ?(\)?ĎC2?zG{?d]??*?Q?UHl?p +=q?Ŗn#pO?\(?(˙?zG?H?fffffg?VƊ+?Q?ƕ}.JD? +=p +?c?\(\?Ȥv?zG?NW C??NjV6,8?(\)?N4t?QR??zG{?=kLM? +=p?wkq??ȱ 1?\(?7[?Q?"+?GzH?YńE?p +=q?ɐ;?񙙙? A?\(?qǥ?Q?1Hk?zG?ff[?=p +>?ʙ*?fffffg?չ0J?\(?6?Q?0yL?Gz?a? +=p +?˒ZY?333333? ]ݼ?\(\?ڶd?Q?w!?zG?M(ӆJ? +=p?z->??̦2)?(\)?1H?QR?/X?zG{?'ǯB? +=p?Q&ل O??z5?\(?͢j?Q? +E%P}?GzH?KHq?p +=q?DA???̜!/=?fffffg?(C?\(?|?Q?1nYH?Gz?Qi? +=p +?qY>?333333?ϐS?\(\?Ϯ?Q?SR?zG?Ye6? +=p?ݦ].??0Jc?(\)?HC+?QR?+C?zG{?8mRKA? +=p?E~X??QOF{?\(?]Hǩ?Q?he6?GzH?tKs>?p +=q?Vx ??Њm?\(?ДY ?Q?ОG?zG?Шx%?=p +>?вnI?fffffg?лF'??y?(\)? /$?QR?Dz?zG{?hS]? +=p?^O8Ψ??#`V/?\(?(kZ?Q?,X?GzH?08:q?p +=q?4o??8z0[y?\(?;?Q??+?zG?Bk?ي?=p +>?D¬*N?fffffg?G5 +?\(?Ijo ?Q?Ka?Qf@QR?Pht@ffffff?Ov%@zG{?NQ@\(?M\\٢@ +=p?K› @Q?Ic@?G}<@Gz?E5@\(?CV@ +=p +?@ @Q?={@333333?:W0@GzH?7uy1@\(\?4o@p +=q?0*@Q?-+5pT@?)CH@zG?%.x#@\(? @ +=p?~*ʡ@Q?m@?/@zG?-2g,@(\)? >Uo@=p +>?͊Q@QR?_8@fffffg?ђI@zG{? +;!Ls@\(?#jX@ +=p?O޵\@Q?6k@?ڈ +=@Gz?RĈ)@\(?c7@ +=p +?ƚIA@Q?п/@333333?и.@GzH?бh^@\(\?Ъ @p +=q?Т a@Q?Кs@?ГE1v@zG?ЋiR:@\(?Ѓn7o@ +=p?{S:c@Q?s)H@?j~t^@zG?bIzd@(\)?YA('@=p +>?Qjn@QR?H5E@fffffg??Jy@zG{?6DC2/@\(?-"؃@ +=p?#}#@Q?e@?K@Gz?* @\(?@ +=p +?f+y@Q?ԼT@333333?$ @GzH?Ϭܕ @\(\?Ϙ@p +=q?τDQq@Q?o@?ZU@zG?FZ@\(?1o@ +=p?F @Q?qŬ@?^ +@zG?=@(\)?jiC?@=p +>?ίrIdׁ@QR?ΙUx^@fffffg?΃+F@@zG{?lꐲ@\(?V/JӦ@ +=p??@Q?(ä#@?~@Gz?u@\(?Q3@ +=p +?uvN@Q?͵@333333?͝a@GzH?͆ Lp@\(\?nYpHC@p +=q?V{uZ@Q?>@?&:N@zG?w@\(?p,@ +=p?-D@Q?j@?̭fh@zG?̔k{<@(\)?|Ls@=p +>?c4"@QR?JW/]@fffffg?2r:\ @zG{?1Bn@\(?=ՌfC@ +=p?9±@Q?%c@?˵&@Gz?˛i@\(?˂zY&@ +=p +?i?h@ Q?O̜@ 333333?6{G'@ GzH? OM@ \(\?$Ž@ p +=q?@ Q?c@ ?ʶTz>@ zG?ʝj @ \(?ʃfmki@ +=p?iƙU@ Q?O9M@ +?6P[@ +zG?I@ +(\)?p=@ +=p +>?wHH@ +QR?άx@ +fffffg?ɴw@ +zG{?ɚӌ4@ +\(?ɀGc@ + +=p?f/pm@ +Q?Lv-@ +?2PG@ +Gz?{@ +\(?9@ +=p +?黙W@ Q?\@ 333333?Ȱ݃?@ GzH?Ȗ]@ \(\?| M@ p +=q?bĥ@ Q?Hz,@ ?.|@ zG?cU@ \(?wD@ +=p?,'S@ Q?ƷW@ ?ǬSS@ zG?ǒyץ@ (\)?x"$@ =p +>?^ߥu'@ QR?DYw@ fffffg?+  @ zG{?()@ \(?Jd@ +=p?r @ Q?à*@ ?ƩѕU@ Gz?Ɛ 5@ \(?vNB@ +=p +?\yZ@ Q?B唦G@ 333333?)?v4@QR?ŮOcI^@fffffg?Ĭ9!@zG{?ēv@\(?z6p@ +=p?aW@Q?H"C@?/ҹ@Gz?@\(?B3@"@ +=p +?{H*@Q? c@333333?ôL0b@GzH?ÛӦ8@\(\?ÃAfJ;@p +=q?jTKU@Q?Rn"S|@?:@zG?!O9@\(? HN"@ +=p?wwv@Q?_@?U!<XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 16 / length of dimension 1 NAXIS2 = 401 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 2 / number of table fields TTYPE1 = 'z ' TFORM1 = 'D ' TTYPE2 = 'nz ' TFORM2 = 'D ' SACCTYPE= 'tracer ' SACCCLSS= 'NZ ' SACCNAME= 'gc_unwise' SACCQTTY= 'galaxy_density' EXTNAME = 'tracer:NZ:gc_unwise' HIERARCH META_ngal = 1.0 END ?zG{? >?zG{?)U}?Q??+M?Gz?ҧ6?ۅQ?ʜğ?(\)?b??-z?p +=q?#X1?zG? e-?޸Q??~v?\(\???ho +?QR?tڿ? +=p?()?\(?)22?GzH?? .?ᙙ?V@?Q?&?=p +>??s?\(?غ.t?Gz?.pu?333333?J )?Q?Pe? +=p?w]퉚?(\)? B?zG{?2-(??;s T?Q?f4S?p +=q?gmOn?\(?8E?zG?(A?fffffg?%U!ӎ?Q?&|? +=p +?&6?\(\?oG?zG? ??T?QR?I)? +=p?6?\(?g!0Q?GzH? ?`?\(?9<:?Gz?~ϯ?333333? +?Q?(S_? +=p?Kp^b ?(\)?ĎC2?zG{?d]??*?Q?UHl?p +=q?Ŗn#pO?\(?(˙?zG?H?fffffg?VƊ+?Q?ƕ}.JD? +=p +?c?\(\?Ȥv?zG?NW C??NjV6,8?(\)?N4t?QR??zG{?=kLM? +=p?wkq??ȱ 1?\(?7[?Q?"+?GzH?YńE?p +=q?ɐ;?񙙙? A?\(?qǥ?Q?1Hk?zG?ff[?=p +>?ʙ*?fffffg?չ0J?\(?6?Q?0yL?Gz?a? +=p +?˒ZY?333333? ]ݼ?\(\?ڶd?Q?w!?zG?M(ӆJ? +=p?z->??̦2)?(\)?1H?QR?/X?zG{?'ǯB? +=p?Q&ل O??z5?\(?͢j?Q? +E%P}?GzH?KHq?p +=q?DA???̜!/=?fffffg?(C?\(?|?Q?1nYH?Gz?Qi? +=p +?qY>?333333?ϐS?\(\?Ϯ?Q?SR?zG?Ye6? +=p?ݦ].??0Jc?(\)?HC+?QR?+C?zG{?8mRKA? +=p?E~X??QOF{?\(?]Hǩ?Q?he6?GzH?tKs>?p +=q?Vx ??Њm?\(?ДY ?Q?ОG?zG?Шx%?=p +>?вnI?fffffg?лF'??y?(\)? /$?QR?Dz?zG{?hS]? +=p?^O8Ψ??#`V/?\(?(kZ?Q?,X?GzH?08:q?p +=q?4o??8z0[y?\(?;?Q??+?zG?Bk?ي?=p +>?D¬*N?fffffg?G5 +?\(?Ijo ?Q?Ka?Qf@QR?Pht@ffffff?Ov%@zG{?NQ@\(?M\\٢@ +=p?K› @Q?Ic@?G}<@Gz?E5@\(?CV@ +=p +?@ @Q?={@333333?:W0@GzH?7uy1@\(\?4o@p +=q?0*@Q?-+5pT@?)CH@zG?%.x#@\(? @ +=p?~*ʡ@Q?m@?/@zG?-2g,@(\)? >Uo@=p +>?͊Q@QR?_8@fffffg?ђI@zG{? +;!Ls@\(?#jX@ +=p?O޵\@Q?6k@?ڈ +=@Gz?RĈ)@\(?c7@ +=p +?ƚIA@Q?п/@333333?и.@GzH?бh^@\(\?Ъ @p +=q?Т a@Q?Кs@?ГE1v@zG?ЋiR:@\(?Ѓn7o@ +=p?{S:c@Q?s)H@?j~t^@zG?bIzd@(\)?YA('@=p +>?Qjn@QR?H5E@fffffg??Jy@zG{?6DC2/@\(?-"؃@ +=p?#}#@Q?e@?K@Gz?* @\(?@ +=p +?f+y@Q?ԼT@333333?$ @GzH?Ϭܕ @\(\?Ϙ@p +=q?τDQq@Q?o@?ZU@zG?FZ@\(?1o@ +=p?F @Q?qŬ@?^ +@zG?=@(\)?jiC?@=p +>?ίrIdׁ@QR?ΙUx^@fffffg?΃+F@@zG{?lꐲ@\(?V/JӦ@ +=p??@Q?(ä#@?~@Gz?u@\(?Q3@ +=p +?uvN@Q?͵@333333?͝a@GzH?͆ Lp@\(\?nYpHC@p +=q?V{uZ@Q?>@?&:N@zG?w@\(?p,@ +=p?-D@Q?j@?̭fh@zG?̔k{<@(\)?|Ls@=p +>?c4"@QR?JW/]@fffffg?2r:\ @zG{?1Bn@\(?=ՌfC@ +=p?9±@Q?%c@?˵&@Gz?˛i@\(?˂zY&@ +=p +?i?h@ Q?O̜@ 333333?6{G'@ GzH? OM@ \(\?$Ž@ p +=q?@ Q?c@ ?ʶTz>@ zG?ʝj @ \(?ʃfmki@ +=p?iƙU@ Q?O9M@ +?6P[@ +zG?I@ +(\)?p=@ +=p +>?wHH@ +QR?άx@ +fffffg?ɴw@ +zG{?ɚӌ4@ +\(?ɀGc@ + +=p?f/pm@ +Q?Lv-@ +?2PG@ +Gz?{@ +\(?9@ +=p +?黙W@ Q?\@ 333333?Ȱ݃?@ GzH?Ȗ]@ \(\?| M@ p +=q?bĥ@ Q?Hz,@ ?.|@ zG?cU@ \(?wD@ +=p?,'S@ Q?ƷW@ ?ǬSS@ zG?ǒyץ@ (\)?x"$@ =p +>?^ߥu'@ QR?DYw@ fffffg?+  @ zG{?()@ \(?Jd@ +=p?r @ Q?à*@ ?ƩѕU@ Gz?Ɛ 5@ \(?vNB@ +=p +?\yZ@ Q?B唦G@ 333333?)?v4@QR?ŮOcI^@fffffg?Ĭ9!@zG{?ēv@\(?z6p@ +=p?aW@Q?H"C@?/ҹ@Gz?@\(?B3@"@ +=p +?{H*@Q? c@333333?ôL0b@GzH?ÛӦ8@\(\?ÃAfJ;@p +=q?jTKU@Q?Rn"S|@?:@zG?!O9@\(? HN"@ +=p?wwv@Q?_@?U!<XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 50 / length of dimension 1 NAXIS2 = 60 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 6 / number of table fields TTYPE1 = 'tracer_0' TFORM1 = '9A ' TTYPE2 = 'tracer_1' TFORM2 = '9A ' TTYPE3 = 'value ' TFORM3 = 'D ' TTYPE4 = 'window_ind' TFORM4 = 'K ' TTYPE5 = 'ell ' TFORM5 = 'D ' TTYPE6 = 'window ' TFORM6 = 'K ' NTRACER = 2 SACCTYPE= 'data ' SACCNAME= 'cl_00 ' EXTNAME = 'data:cl_00' END ck_sock_so>4;Y@.[jOck_sock_so>E߼@F[jOck_sock_so>]K2>@R[jOck_sock_so>?`~n@Z@[jOck_sock_so>xנA<@`[jOck_sock_so>}B-N^J@d[jOck_sock_so>x^C@h`[jOck_sock_so>u^Y@l [jOck_sock_so>r @o[jOck_sock_so>p`BG @q[jOck_sock_so>m l@k +@s[jOck_sock_so>iﯗ @u[jOck_sock_so>gE^. @wp[jOck_sock_so>ds @yP[jOck_sock_so>bFizB@{0[jOck_sock_so>a4@}[jOck_sock_so>_P@~[jOck_sock_so>\U@h[jOck_sock_so>Z0yY@X[jOck_sock_so>XE@H[jOgc_unwiseck_so>~WT@.[jOgc_unwiseck_so>\@F[jOgc_unwiseck_so>F^]@R[jOgc_unwiseck_so>zcB@Z@[jOgc_unwiseck_so>uk<@`[jOgc_unwiseck_so>qʸW9@d[jOgc_unwiseck_so>nb[Q@h`[jOgc_unwiseck_so>iۛӢN*@l [jOgc_unwiseck_so>fr?rD@o[jOgc_unwiseck_so>c+| @q[jOgc_unwiseck_so>aX[ +@s[jOgc_unwiseck_so>]&{ @u[jOgc_unwiseck_so>ZP @wp[jOgc_unwiseck_so>Wm @yP[jOgc_unwiseck_so>TyJ@{0[jOgc_unwiseck_so>RO!@}[jOgc_unwiseck_so>P|O@~[jOgc_unwiseck_so>NMU@h[jOgc_unwiseck_so>L4Kv@X[jOgc_unwiseck_so>I7E@H[jOgc_unwisegc_unwise>@.[jOgc_unwisegc_unwise>-j"@F[jOgc_unwisegc_unwise> 3^_@R[jOgc_unwisegc_unwise>@qۦ @Z@[jOgc_unwisegc_unwise>[$@`[jOgc_unwisegc_unwise>Iv5@d[jOgc_unwisegc_unwise>@h`[jOgc_unwisegc_unwise>3˯@l [jOgc_unwisegc_unwise>s1@o[jOgc_unwisegc_unwise>cDJ @q[jOgc_unwisegc_unwise>~ry6l +@s[jOgc_unwisegc_unwise>}v= @u[jOgc_unwisegc_unwise>|=ޠ @wp[jOgc_unwisegc_unwise>{w @yP[jOgc_unwisegc_unwise>{j䦡@{0[jOgc_unwisegc_unwise>z\@}[jOgc_unwisegc_unwise>zGq/b @~[jOgc_unwisegc_unwise>z>8 E@h[jOgc_unwisegc_unwise>y+s@X[jOgc_unwisegc_unwise>yY".;@H[jOXTENSION= 'IMAGE ' / Image extension BITPIX = -64 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 60 NAXIS2 = 60 PCOUNT = 0 / number of parameters GCOUNT = 1 / number of groups EXTNAME = 'covariance' SACCTYPE= 'cov ' SACCCLSS= 'full ' SIZE = 60 END :<[G9}F<>ԍ <'-PA<6gՙP<78< ]Mm<1gV؃& <13K<@<+JzZ<+ }&b+<Ab8 <%ZVs<%<|kY9< yPe:<~ɼn<^\"?ԍ jrlo;HQ <k +ܷ^<;bg;ލ2M< [ +s<9oPS;T<*/<86 \ No newline at end of file diff --git a/soliket/tests/test_runs.py b/soliket/tests/test_runs.py index 239a7f61..b95bcd7d 100644 --- a/soliket/tests/test_runs.py +++ b/soliket/tests/test_runs.py @@ -10,7 +10,8 @@ "lensing", "lensing_lite", "multi", - "cross_correlation"]) + "cross_correlation", + "xcorr"]) def test_evaluate(lhood): info = yaml_load(pkgutil.get_data("soliket", f"tests/test_{lhood}.yaml")) info["force"] = True @@ -24,7 +25,8 @@ def test_evaluate(lhood): "lensing", "lensing_lite", "multi", - "cross_correlation"]) + "cross_correlation", + "xcorr"]) def test_mcmc(lhood): info = yaml_load(pkgutil.get_data("soliket", f"tests/test_{lhood}.yaml")) info["force"] = True diff --git a/soliket/tests/test_xcorr.py b/soliket/tests/test_xcorr.py new file mode 100644 index 00000000..7877c866 --- /dev/null +++ b/soliket/tests/test_xcorr.py @@ -0,0 +1,171 @@ +# pytest -k xcorr -v --pdb . + +import pytest +import numpy as np + +from cobaya.yaml import yaml_load +from cobaya.model import get_model + +import os +import pdb + + +def get_demo_xcorr_model(theory): + if theory == "camb": + info_yaml = r""" + likelihood: + soliket.XcorrLikelihood: + stop_at_error: True + datapath: soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits + k_tracer_name: ck_so + gc_tracer_name: gc_unwise + + theory: + camb: + extra_args: + lens_potential_accuracy: 1 + + params: + tau: 0.05 + mnu: 0.0 + nnu: 3.046 + b1: + prior: + min: 0. + max: 10. + ref: + min: 1. + max: 4. + proposal: 0.1 + s1: + prior: + min: 0.1 + max: 1.0 + proposal: 0.1 + """ + elif theory == "classy": + info_yaml = r""" + likelihood: + soliket.XcorrLikelihood: + stop_at_error: True + datapath: soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits + k_tracer_name: ck_so + gc_tracer_name: gc_unwise + + theory: + classy: + extra_args: + output: lCl, tCl + path: global + + params: + b1: + prior: + min: 0. + max: 10. + ref: + min: 1. + max: 4. + proposal: 0.1 + s1: + prior: + min: 0.1 + max: 1.0 + proposal: 0.1 + + """ + + info = yaml_load(info_yaml) + model = get_model(info) + return model + + +@pytest.mark.parametrize("theory", ["camb"])#, "classy"]) +def test_xcorr(theory): + + params = {'b1': 1.0, 's1': 0.4} + + model = get_demo_xcorr_model(theory) + + lnl = model.loglike(params)[0] + assert np.isfinite(lnl) + + xcorr_lhood = model.likelihood['soliket.XcorrLikelihood'] + + setup_chi_out = xcorr_lhood._setup_chi() + + Pk_interpolator = xcorr_lhood.theory.get_Pk_interpolator(("delta_nonu", "delta_nonu"), + extrap_kmax=1.e8, + nonlinear=False).P + + from soliket.xcorr.limber import do_limber + + cl_gg, cl_kappag = do_limber(xcorr_lhood.ell_range, + xcorr_lhood.provider, + xcorr_lhood.dndz, + xcorr_lhood.dndz, + params['s1'], + params['s1'], + Pk_interpolator, + params['b1'], + params['b1'], + xcorr_lhood.alpha_auto, + xcorr_lhood.alpha_cross, + setup_chi_out, + Nchi=xcorr_lhood.Nchi, + dndz1_mag=xcorr_lhood.dndz, + dndz2_mag=xcorr_lhood.dndz) + + ell_load = xcorr_lhood.data.x + cl_load = xcorr_lhood.data.y + # cov_load = xcorr_lhood.data.cov + # cl_err_load = np.sqrt(np.diag(cov_load)) + n_ell = len(ell_load) // 2 + + ell_obs_gg = ell_load[n_ell:] + ell_obs_kappag = ell_load[:n_ell] + + cl_obs_gg = cl_load[:n_ell] + cl_obs_kappag = cl_load[n_ell:] + + # Nell_unwise_g = np.ones_like(cl_gg) \ + # / (xcorr_lhood.ngal * (60 * 180 / np.pi)**2) + Nell_obs_unwise_g = np.ones_like(cl_obs_gg) \ + / (xcorr_lhood.ngal * (60 * 180 / np.pi)**2) + + import pyccl as ccl + h2 = (xcorr_lhood.provider.get_param('H0') / 100)**2 + + cosmo = ccl.Cosmology(Omega_c=xcorr_lhood.provider.get_param('omch2') / h2, + Omega_b=xcorr_lhood.provider.get_param('ombh2') / h2, + h=xcorr_lhood.provider.get_param('H0') / 100, + n_s=xcorr_lhood.provider.get_param('ns'), + A_s=xcorr_lhood.provider.get_param('As'), + Omega_k=xcorr_lhood.provider.get_param('omk'), + Neff=xcorr_lhood.provider.get_param('nnu'), + matter_power_spectrum='linear') + + g_bias_zbz = (xcorr_lhood.dndz[:, 0], + params['b1'] * np.ones(len(xcorr_lhood.dndz[:, 0]))) + mag_bias_zbz = (xcorr_lhood.dndz[:, 0], + params['s1'] * np.ones(len(xcorr_lhood.dndz[:, 0]))) + + tracer_g = ccl.NumberCountsTracer(cosmo, + has_rsd=False, + dndz=xcorr_lhood.dndz.T, + bias=g_bias_zbz, + mag_bias=mag_bias_zbz) + + tracer_k = ccl.CMBLensingTracer(cosmo, z_source=1100) + + cl_gg_ccl = ccl.cls.angular_cl(cosmo, tracer_g, tracer_g, xcorr_lhood.ell_range) + cl_kappag_ccl = ccl.cls.angular_cl(cosmo, tracer_k, tracer_g, xcorr_lhood.ell_range) + + assert np.allclose(cl_gg_ccl, cl_gg) + assert np.allclose(cl_kappag_ccl, cl_kappag) + + cl_obs_gg_ccl = ccl.cls.angular_cl(cosmo, tracer_g, tracer_g, ell_obs_gg) + cl_obs_kappag_ccl = ccl.cls.angular_cl(cosmo, tracer_k, tracer_g, ell_obs_kappag) + + assert np.allclose(cl_obs_gg_ccl + Nell_obs_unwise_g, cl_obs_gg) + assert np.allclose(cl_obs_kappag_ccl, cl_obs_kappag) diff --git a/soliket/tests/test_xcorr.yaml b/soliket/tests/test_xcorr.yaml new file mode 100644 index 00000000..d666a950 --- /dev/null +++ b/soliket/tests/test_xcorr.yaml @@ -0,0 +1,33 @@ +# A simple cobaya likelihood for SO + +debug: True + +likelihood: + soliket.XcorrLikelihood: + stop_at_error: True + +params: + b1: + prior: + min: 0. + max: 10. + ref: + min: 1. + max: 4. + proposal: 0.1 + latex: b_1 + s1: + value: 0.4 + latex: s_1 + +theory: + camb: + stop_at_error: False + extra_args: + lens_potential_accuracy: 1 + +sampler: + evaluate: + + +output: chains/test_xcorr diff --git a/soliket/xcorr/XcorrLikelihood.yaml b/soliket/xcorr/XcorrLikelihood.yaml new file mode 100644 index 00000000..b20a2972 --- /dev/null +++ b/soliket/xcorr/XcorrLikelihood.yaml @@ -0,0 +1,41 @@ +auto_file: +cross_file: +dndz_file: + +datapath: soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits +k_tracer_name: ck_so +gc_tracer_name: gc_unwise +# tracer_combinations: gc_unwise-ck_so, gc_unwise-gc_unwise + +high_ell: 600 +nz: 149 +Nchi: 149 +Nchi_mag: 149 + +Pk_interp_kmax: 10.0 + +params: + b1: + prior: + min: 0. + max: 10. + ref: + min: 1. + max: 4. + proposal: 0.1 + s1: + value: 0.4 + H0: + value: 70. + omegam: + value: 0.3 + ombh2: + value: 0.0245 + omch2: + value: 0.1225 + omk: + value: 0.0 + ns: + value: 0.965 + As: + value: 2.11e-9 diff --git a/soliket/xcorr/__init__.py b/soliket/xcorr/__init__.py new file mode 100644 index 00000000..ebb75aa1 --- /dev/null +++ b/soliket/xcorr/__init__.py @@ -0,0 +1 @@ +from .xcorr import XcorrLikelihood diff --git a/soliket/xcorr/data/clgg.txt b/soliket/xcorr/data/clgg.txt new file mode 100755 index 00000000..cb9fe448 --- /dev/null +++ b/soliket/xcorr/data/clgg.txt @@ -0,0 +1,3 @@ +3.550000000000000000e+01 7.650000000000000000e+01 1.265000000000000000e+02 1.765000000000000000e+02 2.265000000000000000e+02 2.765000000000000000e+02 3.265000000000000000e+02 3.765000000000000000e+02 4.265000000000000000e+02 4.765000000000000000e+02 5.265000000000000000e+02 5.765000000000000000e+02 +1.171259386682930336e-07 1.074925971098873865e-07 7.605910992260147274e-08 5.561247886809266997e-08 4.345656443980961222e-08 3.471274000295990805e-08 2.776643347807964378e-08 2.248227907027175782e-08 1.867015082977402487e-08 1.590664246359614681e-08 1.377366655041607371e-08 1.201950138681249373e-08 +1.171386555992278504e-08 1.074818774537805139e-08 7.605695767423656527e-09 5.562085700838220625e-09 4.345907898321461207e-09 3.470700836078987746e-09 2.777317313635696080e-09 2.248324340497516305e-09 1.866778760685389341e-09 1.590802127607215577e-09 1.377460168904915491e-09 1.201961522334897607e-09 diff --git a/soliket/xcorr/data/clkg.txt b/soliket/xcorr/data/clkg.txt new file mode 100755 index 00000000..c962e2ae --- /dev/null +++ b/soliket/xcorr/data/clkg.txt @@ -0,0 +1,3 @@ +3.550000000000000000e+01 7.650000000000000000e+01 1.265000000000000000e+02 1.765000000000000000e+02 2.265000000000000000e+02 2.765000000000000000e+02 3.265000000000000000e+02 3.765000000000000000e+02 4.265000000000000000e+02 4.765000000000000000e+02 5.265000000000000000e+02 5.765000000000000000e+02 +1.425194692235914338e-07 1.234653760415391458e-07 8.611126066544479263e-08 6.248167114060786722e-08 4.838074702560995877e-08 3.853147673160977486e-08 3.090768970482272427e-08 2.506062365326854316e-08 2.078165440429319254e-08 1.768399884093924078e-08 1.532187351346465496e-08 1.340877933987993602e-08 +1.424982171873059966e-08 1.234732647730792068e-08 8.611508401982062142e-09 6.248438515784782318e-09 4.837007692483073648e-09 3.852305256406458467e-09 3.090458671774709973e-09 2.505965588344152520e-09 2.078012214338600843e-09 1.768136302352318975e-09 1.532330752048001851e-09 1.340720288202413845e-09 diff --git a/soliket/xcorr/data/dndz.txt b/soliket/xcorr/data/dndz.txt new file mode 100755 index 00000000..c59e742f --- /dev/null +++ b/soliket/xcorr/data/dndz.txt @@ -0,0 +1,401 @@ +0.000000000000000000e+00 0.000000000000000000e+00 +1.000000000000000021e-02 4.950249168745840913e-05 +2.000000000000000042e-02 1.960397346613510525e-04 +2.999999999999999889e-02 4.367004900968286698e-04 +4.000000000000000083e-02 7.686315513218585889e-04 +5.000000000000000278e-02 1.189036780625892789e-03 +5.999999999999999778e-02 1.695176160451647689e-03 +7.000000000000000666e-02 2.284364858769573781e-03 +8.000000000000000167e-02 2.953972308437234542e-03 +8.999999999999999667e-02 3.701421300348473834e-03 +1.000000000000000056e-01 4.524187090179798351e-03 +1.100000000000000006e-01 5.419796518543995231e-03 +1.199999999999999956e-01 6.385827144363533450e-03 +1.300000000000000044e-01 7.419906391278744198e-03 +1.400000000000000133e-01 8.519710706908299042e-03 +1.499999999999999944e-01 9.682964734781900756e-03 +1.600000000000000033e-01 1.090744049876750561e-02 +1.700000000000000122e-01 1.219095659981774687e-02 +1.799999999999999933e-01 1.353137742486260585e-02 +1.900000000000000022e-01 1.492661236767768920e-02 +2.000000000000000111e-01 1.637461506155963795e-02 +2.099999999999999922e-01 1.787338262364262381e-02 +2.200000000000000011e-01 1.942095491069197738e-02 +2.300000000000000100e-01 2.101541378621318498e-02 +2.399999999999999911e-01 2.265488239871673790e-02 +2.500000000000000000e-01 2.433752447098140245e-02 +2.600000000000000089e-01 2.606154360016054231e-02 +2.700000000000000178e-01 2.782518256857829941e-02 +2.800000000000000266e-01 2.962672266506444219e-02 +2.899999999999999800e-01 3.146448301667866682e-02 +2.999999999999999889e-01 3.333681993067730276e-02 +3.099999999999999978e-01 3.524212624657709764e-02 +3.200000000000000067e-01 3.717883069817298075e-02 +3.300000000000000155e-01 3.914539728536838709e-02 +3.400000000000000244e-01 4.114032465567884350e-02 +3.500000000000000333e-01 4.316214549527120498e-02 +3.599999999999999867e-01 4.520942592940280919e-02 +3.699999999999999956e-01 4.728076493212692716e-02 +3.800000000000000044e-01 4.937479374513208813e-02 +3.900000000000000133e-01 5.149017530558542416e-02 +4.000000000000000222e-01 5.362560368285115842e-02 +4.100000000000000311e-01 5.577980352395765090e-02 +4.199999999999999845e-01 5.795152950768799743e-02 +4.299999999999999933e-01 6.013956580717060713e-02 +4.400000000000000022e-01 6.234272556084808486e-02 +4.500000000000000111e-01 6.455985035170455633e-02 +4.600000000000000200e-01 6.678980969463277351e-02 +4.700000000000000289e-01 6.903150053182431634e-02 +4.799999999999999822e-01 7.128384673606742716e-02 +4.899999999999999911e-01 7.354579862183914518e-02 +5.000000000000000000e-01 7.581633246407917803e-02 +5.100000000000000089e-01 7.809445002453517526e-02 +5.200000000000000178e-01 8.037917808557028254e-02 +5.300000000000000266e-01 8.266956799132499367e-02 +5.400000000000000355e-01 8.496469519612769028e-02 +5.500000000000000444e-01 8.726365882004862018e-02 +5.600000000000000533e-01 8.956558121149418850e-02 +5.700000000000000622e-01 9.186960751673982351e-02 +5.799999999999999600e-01 9.417490525630062281e-02 +5.899999999999999689e-01 9.648066390804095616e-02 +5.999999999999999778e-01 9.878609449692474231e-02 +6.099999999999999867e-01 1.010904291913106851e-01 +6.199999999999999956e-01 1.033929209056964338e-01 +6.300000000000000044e-01 1.056928429098187572e-01 +6.400000000000000133e-01 1.079894884440163461e-01 +6.500000000000000222e-01 1.102821703407646514e-01 +6.600000000000000311e-01 1.125702206522920940e-01 +6.700000000000000400e-01 1.148529902841894657e-01 +6.800000000000000488e-01 1.171298486349243301e-01 +6.900000000000000577e-01 1.194001832411745179e-01 +7.000000000000000666e-01 1.216633994288953413e-01 +7.099999999999999645e-01 1.239189199700362481e-01 +7.199999999999999734e-01 1.261661847448246609e-01 +7.299999999999999822e-01 1.284046504095344154e-01 +7.399999999999999911e-01 1.306337900696591947e-01 +7.500000000000000000e-01 1.328530929584103848e-01 +7.600000000000000089e-01 1.350620641204617889e-01 +7.700000000000000178e-01 1.372602241008635537e-01 +7.800000000000000266e-01 1.394471086390489944e-01 +7.900000000000000355e-01 1.416222683678591554e-01 +8.000000000000000444e-01 1.437852685175109169e-01 +8.100000000000000533e-01 1.459356886244358720e-01 +8.200000000000000622e-01 1.480731222449169771e-01 +8.300000000000000711e-01 1.501971766734529479e-01 +8.399999999999999689e-01 1.523074726657792921e-01 +8.499999999999999778e-01 1.544036441664775050e-01 +8.599999999999999867e-01 1.564853380411035111e-01 +8.699999999999999956e-01 1.585522138127689884e-01 +8.800000000000000044e-01 1.606039434031082958e-01 +8.900000000000000133e-01 1.626402108775664435e-01 +9.000000000000000222e-01 1.646607121949426600e-01 +9.100000000000000311e-01 1.666651549611269933e-01 +9.200000000000000400e-01 1.686532581869663872e-01 +9.300000000000000488e-01 1.706247520501989190e-01 +9.400000000000000577e-01 1.725793776613946395e-01 +9.500000000000000666e-01 1.745168868338436707e-01 +9.599999999999999645e-01 1.764370418573316390e-01 +9.699999999999999734e-01 1.783396152757439623e-01 +9.799999999999999822e-01 1.802243896684420643e-01 +9.899999999999999911e-01 1.820911574353535067e-01 +1.000000000000000000e+00 1.839397205857211670e-01 +1.010000000000000009e+00 1.857698905304554782e-01 +1.020000000000000018e+00 1.875814878780353301e-01 +1.030000000000000027e+00 1.893743422339042304e-01 +1.040000000000000036e+00 1.911482920033083233e-01 +1.050000000000000044e+00 1.929031841975243911e-01 +1.060000000000000053e+00 1.946388742434262686e-01 +1.070000000000000062e+00 1.963552257963386283e-01 +1.080000000000000071e+00 1.980521105561285078e-01 +1.090000000000000080e+00 1.997294080864849553e-01 +1.100000000000000089e+00 2.013870056373381623e-01 +1.110000000000000098e+00 2.030247979703702632e-01 +1.120000000000000107e+00 2.046426871875703934e-01 +1.130000000000000115e+00 2.062405825627874056e-01 +1.140000000000000124e+00 2.078184003762342802e-01 +1.150000000000000133e+00 2.093760637518989332e-01 +1.159999999999999920e+00 2.109135024978168460e-01 +1.169999999999999929e+00 2.124306529491611917e-01 +1.179999999999999938e+00 2.139274578141075467e-01 +1.189999999999999947e+00 2.154038660224298074e-01 +1.199999999999999956e+00 2.168598325767855384e-01 +1.209999999999999964e+00 2.182953184066490637e-01 +1.219999999999999973e+00 2.197102902248513623e-01 +1.229999999999999982e+00 2.211047203866861044e-01 +1.239999999999999991e+00 2.224785867515421123e-01 +1.250000000000000000e+00 2.238318725470234971e-01 +1.260000000000000009e+00 2.251645662355177546e-01 +1.270000000000000018e+00 2.264766613831743769e-01 +1.280000000000000027e+00 2.277681565312566492e-01 +1.290000000000000036e+00 2.290390550698284400e-01 +1.300000000000000044e+00 2.302893651137406517e-01 +1.310000000000000053e+00 2.315190993808805553e-01 +1.320000000000000062e+00 2.327282750726488436e-01 +1.330000000000000071e+00 2.339169137566293188e-01 +1.340000000000000080e+00 2.350850412514167154e-01 +1.350000000000000089e+00 2.362326875135686566e-01 +1.360000000000000098e+00 2.373598865266485225e-01 +1.370000000000000107e+00 2.384666761923254497e-01 +1.380000000000000115e+00 2.395530982235001272e-01 +1.390000000000000124e+00 2.406191980394230823e-01 +1.400000000000000133e+00 2.416650246627743692e-01 +1.409999999999999920e+00 2.426906306186741302e-01 +1.419999999999999929e+00 2.436960718355921374e-01 +1.429999999999999938e+00 2.446814075481268114e-01 +1.439999999999999947e+00 2.456467002016238244e-01 +1.449999999999999956e+00 2.465920153586047825e-01 +1.459999999999999964e+00 2.475174216069769284e-01 +1.469999999999999973e+00 2.484229904699957470e-01 +1.479999999999999982e+00 2.493087963179516919e-01 +1.489999999999999991e+00 2.501749162815539718e-01 +1.500000000000000000e+00 2.510214301669835280e-01 +1.510000000000000009e+00 2.518484203725890880e-01 +1.520000000000000018e+00 2.526559718071984539e-01 +1.530000000000000027e+00 2.534441718100204910e-01 +1.540000000000000036e+00 2.542131100721104486e-01 +1.550000000000000044e+00 2.549628785593751057e-01 +1.560000000000000053e+00 2.556935714370905388e-01 +1.570000000000000062e+00 2.564052849959098657e-01 +1.580000000000000071e+00 2.570981175793355389e-01 +1.590000000000000080e+00 2.577721695126324875e-01 +1.600000000000000089e+00 2.584275430331589574e-01 +1.610000000000000098e+00 2.590643422220910708e-01 +1.620000000000000107e+00 2.596826729375191767e-01 +1.630000000000000115e+00 2.602826427488927341e-01 +1.640000000000000124e+00 2.608643608727914676e-01 +1.650000000000000133e+00 2.614279381100015343e-01 +1.660000000000000142e+00 2.619734867838749426e-01 +1.669999999999999929e+00 2.625011206799504615e-01 +1.679999999999999938e+00 2.630109549868153151e-01 +1.689999999999999947e+00 2.635031062381882894e-01 +1.699999999999999956e+00 2.639776922562015460e-01 +1.709999999999999964e+00 2.644348320958633258e-01 +1.719999999999999973e+00 2.648746459906807238e-01 +1.729999999999999982e+00 2.652972552994236621e-01 +1.739999999999999991e+00 2.657027824540098560e-01 +1.750000000000000000e+00 2.660913509084941175e-01 +1.760000000000000009e+00 2.664630850891406832e-01 +1.770000000000000018e+00 2.668181103455626313e-01 +1.780000000000000027e+00 2.671565529029090169e-01 +1.790000000000000036e+00 2.674785398150825166e-01 +1.800000000000000044e+00 2.677841989189702065e-01 +1.810000000000000053e+00 2.680736587896693779e-01 +1.820000000000000062e+00 2.683470486966932911e-01 +1.830000000000000071e+00 2.686044985611380498e-01 +1.840000000000000080e+00 2.688461389137953295e-01 +1.850000000000000089e+00 2.690721008541952508e-01 +1.860000000000000098e+00 2.692825160105618010e-01 +1.870000000000000107e+00 2.694775165006665363e-01 +1.880000000000000115e+00 2.696572348935651897e-01 +1.890000000000000124e+00 2.698218041722001415e-01 +1.900000000000000133e+00 2.699713576968562623e-01 +1.910000000000000142e+00 2.701060291694533100e-01 +1.919999999999999929e+00 2.702259525986613897e-01 +1.929999999999999938e+00 2.703312622658249897e-01 +1.939999999999999947e+00 2.704220926916813816e-01 +1.949999999999999956e+00 2.704985786038589524e-01 +1.959999999999999964e+00 2.705608549051431999e-01 +1.969999999999999973e+00 2.706090566424955712e-01 +1.979999999999999982e+00 2.706433189768120973e-01 +1.989999999999999991e+00 2.706637771534095016e-01 +2.000000000000000000e+00 2.706705664732254046e-01 +2.010000000000000231e+00 2.706638222647194669e-01 +2.020000000000000018e+00 2.706436798564643142e-01 +2.030000000000000249e+00 2.706102745504125329e-01 +2.040000000000000036e+00 2.705637415958286329e-01 +2.049999999999999822e+00 2.705042161638736009e-01 +2.060000000000000053e+00 2.704318333228306614e-01 +2.069999999999999840e+00 2.703467280139603135e-01 +2.080000000000000071e+00 2.702490350279734832e-01 +2.089999999999999858e+00 2.701388889821121908e-01 +2.100000000000000089e+00 2.700164242978251306e-01 +2.109999999999999876e+00 2.698817751790292818e-01 +2.120000000000000107e+00 2.697350755909456166e-01 +2.129999999999999893e+00 2.695764592394986336e-01 +2.140000000000000124e+00 2.694060595512700051e-01 +2.149999999999999911e+00 2.692240096539948446e-01 +2.160000000000000142e+00 2.690304423575922699e-01 +2.169999999999999929e+00 2.688254901357189919e-01 +2.180000000000000160e+00 2.686092851078369370e-01 +2.189999999999999947e+00 2.683819590217855211e-01 +2.200000000000000178e+00 2.681436432368480283e-01 +2.209999999999999964e+00 2.678944687073046560e-01 +2.220000000000000195e+00 2.676345659664614574e-01 +2.229999999999999982e+00 2.673640651111473532e-01 +2.240000000000000213e+00 2.670830957866694999e-01 +2.250000000000000000e+00 2.667917871722190748e-01 +2.260000000000000231e+00 2.664902679667189300e-01 +2.270000000000000018e+00 2.661786663751035120e-01 +2.280000000000000249e+00 2.658571100950249955e-01 +2.290000000000000036e+00 2.655257263039749738e-01 +2.300000000000000266e+00 2.651846416468158218e-01 +2.310000000000000053e+00 2.648339822237126273e-01 +2.319999999999999840e+00 2.644738735784577521e-01 +2.330000000000000071e+00 2.641044406871821959e-01 +2.339999999999999858e+00 2.637258079474437689e-01 +2.350000000000000089e+00 2.633380991676864125e-01 +2.359999999999999876e+00 2.629414375570627294e-01 +2.370000000000000107e+00 2.625359457156132836e-01 +2.379999999999999893e+00 2.621217456247942335e-01 +2.390000000000000124e+00 2.616989586383476896e-01 +2.399999999999999911e+00 2.612677054735080273e-01 +2.410000000000000142e+00 2.608281062025357699e-01 +2.419999999999999929e+00 2.603802802445751019e-01 +2.430000000000000160e+00 2.599243463578258528e-01 +2.439999999999999947e+00 2.594604226320258999e-01 +2.450000000000000178e+00 2.589886264812357730e-01 +2.459999999999999964e+00 2.585090746369205772e-01 +2.470000000000000195e+00 2.580218831413230718e-01 +2.479999999999999982e+00 2.575271673411210216e-01 +2.490000000000000213e+00 2.570250418813640469e-01 +2.500000000000000000e+00 2.565156206996837551e-01 +2.510000000000000231e+00 2.559990170207708804e-01 +2.520000000000000018e+00 2.554753433511154359e-01 +2.530000000000000249e+00 2.549447114740023279e-01 +2.540000000000000036e+00 2.544072324447592681e-01 +2.550000000000000266e+00 2.538630165862491572e-01 +2.560000000000000053e+00 2.533121734846046080e-01 +2.569999999999999840e+00 2.527548119851965591e-01 +2.580000000000000071e+00 2.521910401888343700e-01 +2.589999999999999858e+00 2.516209654481909030e-01 +2.600000000000000089e+00 2.510446943644485396e-01 +2.609999999999999876e+00 2.504623327841605795e-01 +2.620000000000000107e+00 2.498739857963244437e-01 +2.629999999999999893e+00 2.492797577296606004e-01 +2.640000000000000124e+00 2.486797521500941355e-01 +2.649999999999999911e+00 2.480740718584334148e-01 +2.660000000000000142e+00 2.474628188882417312e-01 +2.669999999999999929e+00 2.468460945038983279e-01 +2.680000000000000160e+00 2.462239991988428855e-01 +2.689999999999999947e+00 2.455966326940014766e-01 +2.700000000000000178e+00 2.449640939363879055e-01 +2.709999999999999964e+00 2.443264810978776347e-01 +2.720000000000000195e+00 2.436838915741498246e-01 +2.729999999999999982e+00 2.430364219837937989e-01 +2.740000000000000213e+00 2.423841681675758253e-01 +2.750000000000000000e+00 2.417272251878629930e-01 +2.760000000000000231e+00 2.410656873281995516e-01 +2.770000000000000018e+00 2.403996480930334634e-01 +2.780000000000000249e+00 2.397292002075879502e-01 +2.790000000000000036e+00 2.390544356178762597e-01 +2.800000000000000266e+00 2.383754454908544318e-01 +2.810000000000000053e+00 2.376923202147103731e-01 +2.819999999999999840e+00 2.370051493992846869e-01 +2.830000000000000071e+00 2.363140218766204281e-01 +2.839999999999999858e+00 2.356190257016389233e-01 +2.850000000000000089e+00 2.349202481529376885e-01 +2.859999999999999876e+00 2.342177757337083610e-01 +2.870000000000000107e+00 2.335116941727704842e-01 +2.879999999999999893e+00 2.328020884257193668e-01 +2.890000000000000124e+00 2.320890426761838266e-01 +2.899999999999999911e+00 2.313726403371924034e-01 +2.910000000000000142e+00 2.306529640526433778e-01 +2.919999999999999929e+00 2.299300956988777622e-01 +2.930000000000000160e+00 2.292041163863512421e-01 +2.939999999999999947e+00 2.284751064614025939e-01 +2.950000000000000178e+00 2.277431455081164180e-01 +2.959999999999999964e+00 2.270083123502773248e-01 +2.970000000000000195e+00 2.262706850534126624e-01 +2.979999999999999982e+00 2.255303409269223403e-01 +2.990000000000000213e+00 2.247873565262918372e-01 +3.000000000000000000e+00 2.240418076553877536e-01 +3.010000000000000231e+00 2.232937693688316072e-01 +3.020000000000000018e+00 2.225433159744515110e-01 +3.030000000000000249e+00 2.217905210358078194e-01 +3.040000000000000036e+00 2.210354573747918272e-01 +3.050000000000000266e+00 2.202781970742941908e-01 +3.060000000000000053e+00 2.195188114809420998e-01 +3.070000000000000284e+00 2.187573712079019239e-01 +3.080000000000000071e+00 2.179939461377464749e-01 +3.089999999999999858e+00 2.172286054253836085e-01 +3.100000000000000089e+00 2.164614175010452501e-01 +3.109999999999999876e+00 2.156924500733346239e-01 +3.120000000000000107e+00 2.149217701323291041e-01 +3.129999999999999893e+00 2.141494439527378280e-01 +3.140000000000000124e+00 2.133755370971114895e-01 +3.149999999999999911e+00 2.126001144191031200e-01 +3.160000000000000142e+00 2.118232400667772741e-01 +3.169999999999999929e+00 2.110449774859669836e-01 +3.180000000000000160e+00 2.102653894236757581e-01 +3.189999999999999947e+00 2.094845379315236067e-01 +3.200000000000000178e+00 2.087024843692350529e-01 +3.209999999999999964e+00 2.079192894081680620e-01 +3.220000000000000195e+00 2.071350130348819374e-01 +3.229999999999999982e+00 2.063497145547425193e-01 +3.240000000000000213e+00 2.055634525955633563e-01 +3.250000000000000000e+00 2.047762851112818472e-01 +3.260000000000000231e+00 2.039882693856681362e-01 +3.270000000000000018e+00 2.031994620360657033e-01 +3.280000000000000249e+00 2.024099190171622376e-01 +3.290000000000000036e+00 2.016196956247897631e-01 +3.300000000000000266e+00 2.008288464997518008e-01 +3.310000000000000053e+00 2.000374256316776178e-01 +3.320000000000000284e+00 1.992454863629005979e-01 +3.330000000000000071e+00 1.984530813923614245e-01 +3.339999999999999858e+00 1.976602627795328859e-01 +3.350000000000000089e+00 1.968670819483666623e-01 +3.359999999999999876e+00 1.960735896912602083e-01 +3.370000000000000107e+00 1.952798361730426202e-01 +3.379999999999999893e+00 1.944858709349783499e-01 +3.390000000000000124e+00 1.936917428987879886e-01 +3.399999999999999911e+00 1.928975003706847324e-01 +3.410000000000000142e+00 1.921031910454254477e-01 +3.419999999999999929e+00 1.913088620103756965e-01 +3.430000000000000160e+00 1.905145597495870313e-01 +3.439999999999999947e+00 1.897203301478865844e-01 +3.450000000000000178e+00 1.889262184949768442e-01 +3.459999999999999964e+00 1.881322694895454506e-01 +3.470000000000000195e+00 1.873385272433838455e-01 +3.479999999999999982e+00 1.865450352855143323e-01 +3.490000000000000213e+00 1.857518365663230475e-01 +3.500000000000000000e+00 1.849589734617008152e-01 +3.510000000000000231e+00 1.841664877771879427e-01 +3.520000000000000018e+00 1.833744207521245118e-01 +3.530000000000000249e+00 1.825828130638040558e-01 +3.540000000000000036e+00 1.817917048316303186e-01 +3.550000000000000266e+00 1.810011356212762046e-01 +3.560000000000000053e+00 1.802111444488443681e-01 +3.570000000000000284e+00 1.794217697850283566e-01 +3.580000000000000071e+00 1.786330495592739209e-01 +3.589999999999999858e+00 1.778450211639395762e-01 +3.600000000000000089e+00 1.770577214584558023e-01 +3.609999999999999876e+00 1.762711867734826354e-01 +3.620000000000000107e+00 1.754854529150638442e-01 +3.629999999999999893e+00 1.747005551687788594e-01 +3.640000000000000124e+00 1.739165283038896226e-01 +3.649999999999999911e+00 1.731334065774840503e-01 +3.660000000000000142e+00 1.723512237386136403e-01 +3.669999999999999929e+00 1.715700130324259720e-01 +3.680000000000000160e+00 1.707898072042905169e-01 +3.689999999999999947e+00 1.700106385039181767e-01 +3.700000000000000178e+00 1.692325386894731332e-01 +3.709999999999999964e+00 1.684555390316772594e-01 +3.720000000000000195e+00 1.676796703179059822e-01 +3.729999999999999982e+00 1.669049628562755405e-01 +3.740000000000000213e+00 1.661314464797209178e-01 +3.750000000000000000e+00 1.653591505500640324e-01 +3.760000000000000231e+00 1.645881039620719355e-01 +3.770000000000000018e+00 1.638183351475041571e-01 +3.780000000000000249e+00 1.630498720791491996e-01 +3.790000000000000036e+00 1.622827422748497894e-01 +3.800000000000000266e+00 1.615169728015155848e-01 +3.810000000000000053e+00 1.607525902791246686e-01 +3.820000000000000284e+00 1.599896208847115808e-01 +3.830000000000000071e+00 1.592280903563428607e-01 +3.839999999999999858e+00 1.584680239970789883e-01 +3.850000000000000089e+00 1.577094466789228100e-01 +3.859999999999999876e+00 1.569523828467540572e-01 +3.870000000000000107e+00 1.561968565222491556e-01 +3.879999999999999893e+00 1.554428913077868235e-01 +3.890000000000000124e+00 1.546905103903381817e-01 +3.899999999999999911e+00 1.539397365453423761e-01 +3.910000000000000142e+00 1.531905921405657123e-01 +3.919999999999999929e+00 1.524430991399457758e-01 +3.930000000000000160e+00 1.516972791074190641e-01 +3.939999999999999947e+00 1.509531532107323271e-01 +3.950000000000000178e+00 1.502107422252372260e-01 +3.959999999999999964e+00 1.494700665376685056e-01 +3.970000000000000195e+00 1.487311461499045417e-01 +3.979999999999999982e+00 1.479940006827111243e-01 +3.990000000000000213e+00 1.472586493794671714e-01 +4.000000000000000000e+00 1.465251111098734293e-01 diff --git a/soliket/xcorr/limber.py b/soliket/xcorr/limber.py new file mode 100644 index 00000000..153f263d --- /dev/null +++ b/soliket/xcorr/limber.py @@ -0,0 +1,121 @@ +import numpy as np +import pdb +from ..constants import C_HMPC + + +oneover_chmpc = 1. / C_HMPC + + +def mag_bias_kernel(provider, dndz, s1, zatchi, chi_arr, chiprime_arr, zprime_arr): + '''Calculates magnification bias kernel. + + ''' + dndzprime = np.interp(zprime_arr, dndz[:, 0], dndz[:, 1], left=0, right=0) + norm = np.trapz(dndz[:, 1], x=dndz[:, 0]) + dndzprime = dndzprime / norm #TODO check this norm is right + + g_integrand = (chiprime_arr - chi_arr[np.newaxis, :]) / chiprime_arr \ + * (oneover_chmpc * provider.get_param('H0') / 100) \ + * np.sqrt(provider.get_param('omegam') * (1 + zprime_arr)**3. + + 1 - provider.get_param('omegam')) \ + * dndzprime + + g = chi_arr * np.trapz(g_integrand, x=chiprime_arr, axis=0) + + W_mu = (5. * s1 - 2.) * 1.5 * provider.get_param('omegam') \ + * (provider.get_param('H0') / 100)**2 * (oneover_chmpc)**2 \ + * (1. + zatchi(chi_arr)) * g + + return W_mu + + +def do_limber(ell_arr, provider, dndz1, dndz2, s1, s2, pk, b1_HF, b2_HF, + alpha_auto, alpha_cross, + chi_grids, + #use_zeff=True, + Nchi=50, dndz1_mag=None, dndz2_mag=None, normed=False): + + zatchi = chi_grids['zatchi'] + # chiatz = chi_grids['chiatz'] + chi_arr = chi_grids['chival'] + # z_arr = chi_grids['zval'] + chiprime_arr = chi_grids['chivalp'] + zprime_arr = chi_grids['zvalp'] + + chistar = provider.get_comoving_radial_distance(provider.get_param('zstar')) + + # Galaxy kernels, assumed to be b(z) * dN/dz + W_g1 = np.interp(zatchi(chi_arr), dndz1[:, 0], dndz1[:, 1] \ + * provider.get_Hubble(dndz1[:, 0], units='1/Mpc'), left=0, right=0) + if not normed: + W_g1 /= np.trapz(W_g1, x=chi_arr) + + W_g2 = np.interp(zatchi(chi_arr), dndz2[:, 0], dndz2[:, 1] \ + * provider.get_Hubble(dndz2[:, 0], units='1/Mpc'), left=0, right=0) + if not normed: + W_g2 /= np.trapz(W_g2, x=chi_arr) + + W_kappa = (oneover_chmpc)**2. * 1.5 * provider.get_param('omegam') \ + * (provider.get_param('H0') / 100)**2. * (1. + zatchi(chi_arr)) \ + * chi_arr * (chistar - chi_arr) / chistar + + # Get effective redshift + # if use_zeff: + # kern = W_g1 * W_g2 / chi_arr**2 + # zeff = np.trapz(kern * z_arr,x=chi_arr) / np.trapz(kern, x=chi_arr) + # else: + # zeff = -1.0 + + # set up magnification bias kernels + W_mu1 = mag_bias_kernel(provider, dndz1, s1, + zatchi, chi_arr, chiprime_arr, zprime_arr) + + c_ell_g1g1 = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) + c_ell_g1kappa = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) + c_ell_kappakappa = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) + + c_ell_g1mu1 = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) + c_ell_mu1mu1 = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) + c_ell_mu1kappa = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) + + for i_chi, chi in enumerate(chi_arr): + + k_arr = (ell_arr + 0.5) / chi + + p_mm_hf = pk(zatchi(chi), k_arr) + p_mm = p_mm_hf + p_gg = b1_HF * b1_HF * p_mm_hf # lets just stay at constant linear bias for now + p_gm = b1_HF * p_mm_hf + + W_g1g1 = W_g1[i_chi] * W_g1[i_chi] / (chi**2) * p_gg + c_ell_g1g1[:, :, i_chi] = W_g1g1.T + + W_g1kappa = W_g1[i_chi] * W_kappa[i_chi] / (chi**2) * p_gm + c_ell_g1kappa[:, :, i_chi] = W_g1kappa.T + + # W_kappakappa = W_kappa[i_chi] * W_kappa[i_chi] / (chi**2) * p_mm + # c_ell_kappakappa[:,:,i_chi] = W_kappakappa.T + + W_g1mu1 = W_g1[i_chi] * W_mu1[i_chi] / (chi**2) * p_gm + c_ell_g1mu1[:, :, i_chi] = W_g1mu1.T + + W_mu1mu1 = W_mu1[i_chi] * W_mu1[i_chi] / (chi**2) * p_mm + c_ell_mu1mu1[:, :, i_chi] = W_mu1mu1.T + + W_mu1kappa = W_kappa[i_chi] * W_mu1[i_chi] / (chi**2) * p_mm + c_ell_mu1kappa[:, :, i_chi] = W_mu1kappa.T + + + c_ell_g1g1 = np.trapz(c_ell_g1g1, x=chi_arr, axis=-1) + c_ell_g1kappa = np.trapz(c_ell_g1kappa, x=chi_arr, axis=-1) + c_ell_kappakappa = np.trapz(c_ell_kappakappa, x=chi_arr, axis=-1) + + c_ell_g1mu1 = np.trapz(c_ell_g1mu1, x=chi_arr, axis=-1) + c_ell_mu1mu1 = np.trapz(c_ell_mu1mu1, x=chi_arr, axis=-1) + c_ell_mu1kappa = np.trapz(c_ell_mu1kappa, x=chi_arr, axis=-1) + + clobs_gg = c_ell_g1g1 + 2. * c_ell_g1mu1 + c_ell_mu1mu1 + clobs_kappag = c_ell_g1kappa + c_ell_mu1kappa + # clobs_kappakappa = c_ell_kappakappa + + return clobs_gg.flatten(), clobs_kappag.flatten()#, clobs_kappakappa.flatten() diff --git a/soliket/xcorr/xcorr.py b/soliket/xcorr/xcorr.py new file mode 100644 index 00000000..6c918ad9 --- /dev/null +++ b/soliket/xcorr/xcorr.py @@ -0,0 +1,264 @@ +r""" Likelihood for cross-correlation of CMB lensing and galaxy clustering probes. + +""" + +import numpy as np +import sacc +import pdb + +from scipy.interpolate import InterpolatedUnivariateSpline as Spline + +from ..gaussian import GaussianData, GaussianLikelihood +from .. import utils +from .limber import do_limber + + +class XcorrLikelihood(GaussianLikelihood): + '''Cross-correlation Likelihood for CMB lensing and galaxy clustering probes. + + Based on the original xcorr code [1]_ used in [2]_. + + Accepts data files containing the two spectra from either text files or a sacc file. + + Parameters + ---------- + + datapath : str, optional + sacc file containing the redshift distribtion, galaxy-galaxy and galaxy-kappa + observed spectra. Default: soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits + k_tracer_name : str, optional + sacc file tracer name for kappa. Default: ck_so + gc_tracer_name : str, optional + sacc file tracer name for galaxy clustering. Default: gc_unwise + + dndz_file : str, optional + Text file containing the redshift distribution. + auto_file : str, optional + Text file containing the galaxy-galaxy observed spectra. + cross_file : str, optional + Text file containing the galaxy-kappa observed spectra. + + high_ell : int + Maximum multipole to be computed for all spectra. Default: 600 + nz : int + Resolution of redshift grid used for Limber computations. Default: 149 + Nchi : int + Resolution of Chi grid used for lensing kernel computations. Default: 149 + Nchi_mag : int + Resolution of Chi grid used for magnification kernel computations. Default: 149 + + Pk_interp_kmax : float + Maximum k value for the Pk interpolator, units Mpc^-1. Default: 10.0 + + b1 : float + Linear galaxy bias value for the galaxy sample. + s1 : float + Magnification bias slope for the galaxy sample. + + + References + ---------- + .. [1] https://github.com/simonsobs/xcorr + .. [2] Krolewski, Ferraro and White, 2021, arXiv:2105.03421 + + ''' + + def initialize(self): + name: str = "Xcorr" # noqa F841 + self.log.info('Initialising.') + + if self.datapath is None: + + dndz_file: Optional[str] # noqa F821 + auto_file: Optional[str] # noqa F821 + cross_file: Optional[str] # noqa F821 + + self.dndz = np.loadtxt(self.dndz_file) + + self.x, self.y, self.dy = self._get_data() + if self.covpath is None: + self.log.info('No xcorr covariance specified. Using diag(dy^2).') + self.cov = np.diag(self.dy**2) + else: + self.cov = self._get_cov() + + else: + + self.k_tracer_name: Optional[str] # noqa F821 + self.gc_tracer_name: Optional[str] # noqa F821 + # tracer_combinations: Optional[str] # TODO: implement with keep_selection + + self.sacc_data = self._get_sacc_data() + + self.x = self.sacc_data['x'] + self.y = self.sacc_data['y'] + self.cov = self.sacc_data['cov'] + self.dndz = self.sacc_data['dndz'] + self.ngal = self.sacc_data['ngal'] + + # TODO is this resolution limit on zarray a CAMB problem? + self.nz: Optional[int] # noqa F821 + assert self.nz <= 149, "CAMB limitations requires nz <= 149" + self.zarray = np.linspace(self.dndz[:, 0].min(), self.dndz[:, 0].max(), self.nz) + self.zbgdarray = np.concatenate([self.zarray, [1100]]) # TODO: unfix zstar + self.Nchi: Optional[int] # noqa F821 + self.Nchi_mag: Optional[int] # noqa F821 + + #self.use_zeff: Optional[bool] # noqa F821 + + self.Pk_interp_kmax: Optional[float] # noqa F821 + + self.high_ell: Optional[float] # noqa F821 + self.ell_range = np.linspace(1, self.high_ell, int(self.high_ell + 1)) + + # TODO expose these defaults + self.alpha_auto = 0.9981 + self.alpha_cross = 0.9977 + + self.data = GaussianData(self.name, self.x, self.y, self.cov) + + + def get_requirements(self): + return { + 'Cl': {'lmax': self.high_ell, + 'pp': self.high_ell}, + "Pk_interpolator": { + "z": self.zarray[:-1], + "k_max": self.Pk_interp_kmax, + #"extrap_kmax": 20.0, + "nonlinear": False, + "hubble_units": False, # cobaya told me to + "k_hunit": False, # cobaya told me to + "vars_pairs": [["delta_nonu", "delta_nonu"]], + }, + "Hubble": {"z": self.zarray}, + "angular_diameter_distance": {"z": self.zbgdarray}, + "comoving_radial_distance": {"z": self.zbgdarray}, + 'H0': None, + 'ombh2': None, + 'omch2': None, + 'omk': None, + 'omegam': None, + 'zstar': None, + 'As': None, + 'ns': None + } + + def _bin(self, theory_cl, lmin, lmax): + binned_theory_cl = np.zeros_like(lmin) + for i in range(len(lmin)): + binned_theory_cl[i] = np.mean(theory_cl[(self.ell_range >= lmin[i]) + & (self.ell_range < lmax[i])]) + return binned_theory_cl + + def _get_sacc_data(self, **params_values): + + data_sacc = sacc.Sacc.load_fits(self.datapath) + + # TODO: would be better to use keep_selection + data_sacc.remove_selection(tracers=(self.k_tracer_name, self.k_tracer_name)) + + ell_auto, cl_auto = data_sacc.get_ell_cl('cl_00', + self.gc_tracer_name, + self.gc_tracer_name) + ell_cross, cl_cross = data_sacc.get_ell_cl('cl_00', + self.gc_tracer_name, + self.k_tracer_name) #TODO: check order + cov = data_sacc.covariance.covmat + + x = np.concatenate([ell_auto, ell_cross]) + y = np.concatenate([cl_auto, cl_cross]) + + dndz = np.column_stack([data_sacc.tracers[self.gc_tracer_name].z, + data_sacc.tracers[self.gc_tracer_name].nz]) + ngal = data_sacc.tracers[self.gc_tracer_name].metadata['ngal'] + + data = {'x': x, + 'y': y, + 'cov': cov, + 'dndz': dndz, + 'ngal': ngal} + + return data + + + def _get_data(self, **params_values): + + data_auto = np.loadtxt(self.auto_file) + data_cross = np.loadtxt(self.cross_file) + + # Get data + self.ell_auto = data_auto[0] + cl_auto = data_auto[1] + cl_auto_err = data_auto[2] + + self.ell_cross = data_cross[0] + cl_cross = data_cross[1] + cl_cross_err = data_cross[2] + + x = np.concatenate([self.ell_auto, self.ell_cross]) + y = np.concatenate([cl_auto, cl_cross]) + dy = np.concatenate([cl_auto_err, cl_cross_err]) + + return x, y, dy + + def _setup_chi(self): + + chival = self.provider.get_comoving_radial_distance(self.zarray) + zatchi = Spline(chival, self.zarray) + chiatz = Spline(self.zarray, chival) + + chimin = np.min(chival) + 1.e-5 + chimax = np.max(chival) + chival = np.linspace(chimin, chimax, self.Nchi) + zval = zatchi(chival) + chistar = \ + self.provider.get_comoving_radial_distance(self.provider.get_param('zstar')) + chivalp = \ + np.array(list(map(lambda x: np.linspace(x, chistar, self.Nchi_mag), chival))) + chivalp = chivalp.transpose()[0] + zvalp = zatchi(chivalp) + + chi_result = {'zatchi': zatchi, + 'chiatz': chiatz, + 'chival': chival, + 'zval': zval, + 'chivalp': chivalp, + 'zvalp': zvalp} + + return chi_result + + def _get_theory(self, **params_values): + + setup_chi_out = self._setup_chi() + + Pk_interpolator = self.provider.get_Pk_interpolator(("delta_nonu", "delta_nonu"), + extrap_kmax=1.e8, + nonlinear=False).P + + cl_gg, cl_kappag = do_limber(self.ell_range, + self.provider, + self.dndz, + self.dndz, + params_values['s1'], + params_values['s1'], + Pk_interpolator, + params_values['b1'], + params_values['b1'], + self.alpha_auto, + self.alpha_cross, + setup_chi_out, + Nchi=self.Nchi, + #use_zeff=self.use_zeff, + dndz1_mag=self.dndz, + dndz2_mag=self.dndz) + + # TODO: this is not the correct binning, + # but there needs to be a consistent way to specify it + bin_edges = np.linspace(20, self.high_ell, self.data.x.shape[0] // 2 + 1) + + ell_gg, clobs_gg = utils.binner(self.ell_range, cl_gg, bin_edges) + ell_kappag, clobs_kappag = utils.binner(self.ell_range, cl_kappag, bin_edges) + #ell_kappakappa, clobs_kappakappa = utils.binner(self.ell_range, cl_kappakappa, bin_edges) # noqa E501 + + return np.concatenate([clobs_gg, clobs_kappag])