diff --git a/examples/cluster_number_counts/cluster_SDSS_counts_mean_mass_redshift_richness.ini b/examples/cluster_number_counts/cluster_SDSS_counts_mean_mass_redshift_richness.ini new file mode 100644 index 00000000..edfcd655 --- /dev/null +++ b/examples/cluster_number_counts/cluster_SDSS_counts_mean_mass_redshift_richness.ini @@ -0,0 +1,58 @@ +[runtime] +sampler = test +root = ${PWD} + +[default] +fatal_errors = T + +[output] +filename = output_rp/number_counts_samples.txt +format = text +verbosity = 0 + +[pipeline] +modules = consistency camb firecrown_likelihood +values = ${FIRECROWN_DIR}/examples/cluster_number_counts/cluster_richness_values.ini +likelihoods = firecrown +quiet = T +debug = T +timing = T + +[consistency] +file = ${CSL_DIR}/utility/consistency/consistency_interface.py + +[camb] +file = ${CSL_DIR}/boltzmann/camb/camb_interface.py + +mode = all +lmax = 2500 +feedback = 0 +zmin = 0.0 +zmax = 1.0 +nz = 100 +kmin = 1e-4 +kmax = 50.0 +nk = 1000 + +[firecrown_likelihood] +;; Fix this to use an environment variable to find the files. +;; Set FIRECROWN_DIR to the base of the firecrown installation (or build, if you haven't +;; installed it) +file = ${FIRECROWN_DIR}/firecrown/connector/cosmosis/likelihood.py +likelihood_source = ${FIRECROWN_DIR}/examples/cluster_number_counts/cluster_SDSS_redshift_richness.py +sampling_parameters_sections = firecrown_number_counts +use_cluster_counts = True +use_mean_log_mass = True + +[test] +fatal_errors = T +save_dir = output_counts_mean_mass + +[metropolis] +samples = 1000 +nsteps = 1 + +[emcee] +walkers = 20 +samples = 4000 +nsteps = 10 diff --git a/examples/cluster_number_counts/cluster_SDSS_redshift_richness.py b/examples/cluster_number_counts/cluster_SDSS_redshift_richness.py new file mode 100644 index 00000000..67e9607b --- /dev/null +++ b/examples/cluster_number_counts/cluster_SDSS_redshift_richness.py @@ -0,0 +1,58 @@ +"""Likelihood factory function for cluster number counts.""" + +import os + +import pyccl as ccl +import sacc + +from firecrown.likelihood.gaussian import ConstGaussian +from firecrown.likelihood.binned_cluster_number_counts import ( + BinnedClusterNumberCounts, +) +from firecrown.likelihood.likelihood import Likelihood, NamedParameters +from firecrown.modeling_tools import ModelingTools +from firecrown.models.cluster.abundance import ClusterAbundance +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.murata_binned_spec_z import ( + MurataBinnedSpecZRecipe, +) + + +def get_cluster_abundance() -> ClusterAbundance: + """Creates and returns a ClusterAbundance object.""" + hmf = ccl.halos.MassFuncBocquet16() + min_mass, max_mass = 13.0, 16.0 + min_z, max_z = 0.2, 0.8 + cluster_abundance = ClusterAbundance(min_mass, max_mass, min_z, max_z, hmf) + + return cluster_abundance + + +def build_likelihood( + build_parameters: NamedParameters, +) -> tuple[Likelihood, ModelingTools]: + """Builds the likelihood for Firecrown.""" + # Pull params for the likelihood from build_parameters + average_on = ClusterProperty.NONE + if build_parameters.get_bool("use_cluster_counts", True): + average_on |= ClusterProperty.COUNTS + if build_parameters.get_bool("use_mean_log_mass", False): + average_on |= ClusterProperty.MASS + + survey_name = "SDSSCluster_redshift_richness" + likelihood = ConstGaussian( + [BinnedClusterNumberCounts(average_on, survey_name, MurataBinnedSpecZRecipe())] + ) + + # Read in sacc data + sacc_file_nm = "SDSSTEST_redshift_richness_sacc_data.fits" + sacc_path = os.path.expanduser( + os.path.expandvars("${FIRECROWN_DIR}/examples/cluster_number_counts/") + ) + sacc_data = sacc.Sacc.load_fits(os.path.join(sacc_path, sacc_file_nm)) + likelihood.read(sacc_data) + + cluster_abundance = get_cluster_abundance() + modeling_tools = ModelingTools(cluster_abundance=cluster_abundance) + + return likelihood, modeling_tools diff --git a/examples/cluster_number_counts/demo_sacc_for_clusters.ipynb b/examples/cluster_number_counts/demo_sacc_for_clusters.ipynb index d671abd8..56eeaed7 100644 --- a/examples/cluster_number_counts/demo_sacc_for_clusters.ipynb +++ b/examples/cluster_number_counts/demo_sacc_for_clusters.ipynb @@ -26,7 +26,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "7f10789d-5ee5-4ee4-ba3c-b4605909b22b", "metadata": {}, "outputs": [], @@ -44,7 +44,7 @@ "from numcosmo_py import Nc\n", "from numcosmo_py import Ncm\n", "\n", - "from firecrown.sacc_support import sacc\n", + "import sacc\n", "\n", "os.environ[\"CLMM_MODELING_BACKEND\"] = (\n", " \"nc\" # Need to use NumCosmo as CLMM's backend as well.\n", @@ -63,7 +63,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "8fad6f65-7a43-4965-8e27-01fbc4f7afe5", "metadata": {}, "outputs": [], @@ -85,7 +85,8 @@ "\n", "dist = Nc.Distance.new(2.0)\n", "\n", - "tf = Nc.TransferFunc.new_from_name(\"NcTransferFuncEH\")\n", + "# tf = Nc.TransferFunc.new_from_name(\"NcTransferFuncEH\")\n", + "tf = Nc.TransferFuncEH.new()\n", "\n", "psml = Nc.PowspecMLTransfer.new(tf)\n", "\n", @@ -125,7 +126,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "472f1c1a-2aa4-4f01-8988-45f958f5aed4", "metadata": {}, "outputs": [], @@ -137,13 +138,18 @@ "zl = 0.2\n", "zu = 0.65\n", "\n", - "cluster_z = Nc.ClusterRedshift.new_from_name(\n", - " f\"NcClusterRedshiftNodist{{'z-min': <{zl:22.15e}>, 'z-max':<{zu:22.15e}>}}\"\n", - ")\n", + "# cluster_z = Nc.ClusterRedshift.new_from_name(\n", + "# f\"NcClusterRedshiftNodist{{'z-min': <{zl:22.15e}>, 'z-max':<{zu:22.15e}>}}\"\n", + "# )\n", + "\n", + "# cluster_m = Nc.ClusterMass.new_from_name(\n", + "# f\"NcClusterMassAscaso{{'M0':<{3.0e14 / 0.71:22.15e}>,'z0':<0.6>, \"\n", + "# f\"'lnRichness-min':<{lnRl:22.15e}>, 'lnRichness-max':<{lnRu:22.15e}>}}\"\n", + "# )\n", "\n", - "cluster_m = Nc.ClusterMass.new_from_name(\n", - " f\"NcClusterMassAscaso{{'M0':<{3.0e14 / 0.71:22.15e}>,'z0':<0.6>, \"\n", - " f\"'lnRichness-min':<{lnRl:22.15e}>, 'lnRichness-max':<{lnRu:22.15e}>}}\"\n", + "cluster_z = Nc.ClusterRedshiftNodist(z_max=zu, z_min=zl)\n", + "cluster_m = Nc.ClusterMassAscaso(\n", + " M0=3.0e14 / 0.71, z0=0.6, lnRichness_min=lnRl, lnRichness_max=lnRu\n", ")\n", "\n", "# mean richness-mass relation parameters\n", @@ -212,7 +218,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "id": "b0db89c8-f549-497f-bce8-f33b79166701", "metadata": {}, "outputs": [], @@ -259,7 +265,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "b3144e1f-5b5c-420c-83d5-2e7f34d447c6", "metadata": {}, "outputs": [], @@ -269,7 +275,7 @@ "bin_richness_labels = []\n", "\n", "survey_name = \"NC_mock_redshift_richness\"\n", - "s_count.add_tracer(\"cluster_survey\", survey_name, area)\n", + "s_count.add_tracer(\"survey\", survey_name, area)\n", "\n", "for i, z_bin in enumerate(zip(z_edges[:-1], z_edges[1:])):\n", " lower, upper = z_bin\n", @@ -327,7 +333,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "49dfb4d9-ab69-48f0-9005-f0363263d2d7", "metadata": {}, "outputs": [], @@ -337,20 +343,65 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "id": "34090216-8645-4b69-b85d-4a5498497a9e", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "['cluster_counts', 'cluster_mean_log_mass']" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "s2.get_data_types()" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "id": "0aa582db-74f9-4b85-bc25-711644109c3c", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0'), value=74639, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1'), value=26398, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2'), value=3158, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3'), value=105, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0'), value=135594, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1'), value=45119, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2'), value=4894, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3'), value=162, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0'), value=185542, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1'), value=58654, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2'), value=5866, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3'), value=113, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0'), value=13.377663350048888, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1'), value=13.639310773443059, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2'), value=14.01797746367139, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3'), value=14.426967540211269, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0'), value=13.369226957014275, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1'), value=13.610514626747987, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2'), value=13.983144611997671, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3'), value=14.332596440985231, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0'), value=13.359366691889559, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1'), value=13.586210260775394, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2'), value=13.91311887488307, ),\n", + " DataPoint(data_type='cluster_mean_log_mass', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3'), value=14.204315477919724, )]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "s2.data" ] @@ -367,10 +418,39 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "14b1d729-d470-42b0-a679-98180c60708c", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/f9/69jj5zld1gz1dpwjv9vmd77c0000gn/T/ipykernel_31895/1676763377.py:1: RuntimeWarning: divide by zero encountered in log10\n", + " plt.imshow(np.log10(s2.covariance.covmat))\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAGdCAYAAADdSjBDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAilUlEQVR4nO3dcXCU1f3v8c+CZgmaZEwj2eSyhGixWkNpGyzCRQlWUjIdrgr1QnUY8IeMDISRyfVaI52breOQ3ygy3JZKizM3xTtSuZ0qOsUq6UWCXi79BZSRMpYf1HCTCmkKtVmgkpTsuX9Q1i6JIbtnN/scnverc2bcZ/c8z8nOTr98v+c8zwkYY4wAAIAnjcj2AAAAwOcjUAMA4GEEagAAPIxADQCAhxGoAQDwMAI1AAAeRqAGAMDDCNQAAHjYVdkewKVisZiOHz+uvLw8BQKBbA8HAJAkY4xOnz6t0tJSjRiRuXzw3Llz6u3ttT5PTk6ORo0alYYRZYbnAvXx48cVDoezPQwAgKWOjg6NHTs2I+c+d+6cysuuVWdXn/W5QqGQ2traPBusPReo8/LyJEnvt16vvGuT/5fY9aX/nu4hAQCSEI1GFQ6H4/9/ngm9vb3q7OpT2/4y5eelnrVHT8dUXvn/1NvbS6Aeqovl7rxrRygvhS8/Pz8/3UMCAKRgOKYv8/NGWAVqF2Tsr3v++edVXl6uUaNGqbKyUu+8806mLgUA8Kk+E7NuXpeRQL1161atWrVKq1ev1vvvv6877rhDNTU1am9vz8TlAAA+FZOxbl6XkUC9bt06LVmyRA8//LBuueUWrV+/XuFwWBs3bszE5QAAPhVLw/+8Lu2Bure3V/v371d1dXXC8erqau3Zs6ff53t6ehSNRhMaAAC4IO2B+uTJk+rr61NxcXHC8eLiYnV2dvb7fGNjowoKCuKNW7MAAEPVZ4x187qMLSa7dLWfMWbAFYD19fXq7u6Ot46OjkwNCQBwhfHDHHXab88qKirSyJEj+2XPXV1d/bJsSQoGgwoGg+keBgAAV4S0Z9Q5OTmqrKxUc3NzwvHm5mZNmzYt3ZcDAPhYTEZ9Fs2XGbUk1dXVaeHChZo8ebKmTp2qTZs2qb29XcuWLcvE5QAAPmVbvvZtoJ4/f75OnTqlp556SidOnFBFRYXeeOMNlZWVZeJyAABcsTL2CNHly5dr+fLlmTo9AADWK7ddWPXtuWd9X3R96b+n9NzuA+12t3d9dRyrzgHAFbF/NJv+XndlP8kcAADHeTajBgDgci6u3rbp73UEagCAs/rMhWbT3+sI1AAAZzFHDQAAsoqMGgDgrJgC6lP/fSSS6e91BGoAgLNi5kKz6e91lL4BAPAwMmoAgLP6LEvfNn2HC4EaAOAsPwRqSt8AAHgYGTUAwFkxE1DMWKz6tug7XAjUAABnUfoGAABZRUYNAHBWn0aozyLn7EvjWDLligvUtvtJv3D4jpT7Lv3SO1bXBgAkx1jOURvmqAEAyBzmqAEAQFaRUQMAnNVnRqjPWMxRO/CsbwI1AMBZMQUUsygOx+T9SE3pGwAADyNQAwCcdXExmU2z0djYqEAgoFWrVqXnDxoApW8AgLPs56hTL323trZq06ZN+spXvpLyOYaCjBoAgCSdOXNGDz74oF544QVdd911Gb0WgRoA4KwLi8nsmiRFo9GE1tPTM+h1V6xYoW9/+9u6++67M/43UvoGADgrZvkI0YurvsPhcMLxhoYGRSKRAfu8/PLLeu+999Ta2prydZNBoAYA+F5HR4fy8/Pjr4PB4Od+7tFHH9WOHTs0atSoYRkbgRoA4Kx0LSbLz89PCNSfZ//+/erq6lJlZeVn5+jr0+7du7Vhwwb19PRo5MiRKY9nIARqAICzYhoxrA88+eY3v6mDBw8mHHvooYd0880363vf+17ag7REoAYAOKzPBNRnsQNWsn3z8vJUUVGRcOyaa67RF77whX7H04VAfQmbrSr/0zu1Vtd+/Y4NVv0BAFceAjUAwFl9lqu++9LwrO9du3ZZn2MwBGoAgLNiZoRiFovJYhZPJhsuPPAEAAAPI6MGADjLC6XvTCNQAwCcFVPyK7cv7e91lL4BAPAwMmoAgLPsH3ji/XyVQA0AcJb9I0S9H6i9P0IAAHyMjBoA4Kx/3lM61f5eR6AGADjLD6VvAjUAwFn291F7P1B7f4QAAPgYGTUAwFkxE1DM5oEnFn2HC4EaAOCsmGXpm/uofcZ2P+meEzek3DdY8pHVtQEA3kSgBgA4y36bSzJqAAAypk8B9VncC23Td7h4/58SAAD4GBk1AMBZlL4BAPCwPtmVr/vSN5SM8f4/JQAA8DEyagCAsyh9AwDgYWzKAQCAhxnLbS4Nt2cBAAAbZNQAAGdR+gYAwMP8sHuW9/8pAQCAj5FRAwCc1We5zaVN3+FCoPYQm60q/9BRknLfG8MnUu4LANlE6RsAAGQVGTUAwFkxjVDMIue06TtcCNQAAGf1mYD6LMrXNn2HS9r/KRGJRBQIBBJaKBRK92UAAPCFjGTUt956q37zm9/EX48cOTITlwEA+JwfFpNlJFBfddVVZNEAgIwzlrtnGQeeTJaRER45ckSlpaUqLy/XggUL9NFHn3/bUU9Pj6LRaEIDAGAo+hSwbl6X9kA9ZcoUvfjii3rrrbf0wgsvqLOzU9OmTdOpU6cG/HxjY6MKCgriLRwOp3tIAAA4K+2BuqamRvPmzdPEiRN19913a/v27ZKkzZs3D/j5+vp6dXd3x1tHR0e6hwQAuELFzGfz1Km1bP8Fl5fx27OuueYaTZw4UUeOHBnw/WAwqGAwmOlhAACuQDHLOWqbvsMl4yPs6enRhx9+qJKS1B9xCQCAX6U9UD/22GNqaWlRW1ubfvvb3+o73/mOotGoFi1alO5LAQB8LqaAdfO6tJe+//jHP+q73/2uTp48qeuvv16333679u7dq7KysnRfCgDgc354MlnaA/XLL7+c7lMCAOBbPOv7CmGzVeWWI9+wuvYDE/7Nqj8ApMoPi8kI1AAAZ8Vk+QhRB+aovf9PCQAAfIyMGgDgLGO5cts4kFETqAEAzmL3LAAAPMwPi8m8P0IAAHyMjBoA4Cw/lL7JqAEAzhruR4g2NjbqtttuU15ensaMGaN7771Xhw8fztBfdwGBGgCAIWppadGKFSu0d+9eNTc36/z586qurtbZs2czdk1K3wAAZw136fvNN99MeN3U1KQxY8Zo//79uvPOO1Mex2AI1AAAZ6UrUEej0YTjwWBQwWDwsv27u7slSYWFhSmP4XIofQMAfC8cDqugoCDeGhsbL9vHGKO6ujpNnz5dFRUVGRsbGTUAwFnpyqg7OjqUn58fPz6UbLq2tlYffPCB3n333ZSvPxQEagCAs9IVqPPz8xMC9eWsXLlSr7/+unbv3q2xY8emfP2hIFADADBExhitXLlSr776qnbt2qXy8vKMX5NADev9pP/z/30k5b7/a+pPra4NwN+M7LaqNEl+fsWKFdqyZYtee+015eXlqbOzU5JUUFCg3NzclMcxGAI1AMBZw3171saNGyVJVVVVCcebmpq0ePHilMcxGAI1AMBZwx2ojUk2B7fH7VkAAHgYGTUAwFl+2JSDQA0AcJYfAjWlbwAAPIyMGgDgLGMCMhZZsU3f4UKgBgA4K5U9pS/t73WUvgEA8DAyagCAs/ywmIxADQBwlh/mqCl9AwDgYWTUAABnUfoGAMDD/FD6JlDDms1WlbHOCVbXHhE6YtUfgNuMZUbtQqBmjhoAAA8jowYAOMtIstl5cvg3rUwegRoA4KyYAgrwZDIAAJAtZNQAAGex6hsAAA+LmYACV/h91JS+AQDwMDJqAICzjLFc9e3Asm8CNQDAWX6Yo6b0DQCAh5FRAwCc5YeMmkANAHCWH1Z9E6gBAM7yw2Iy5qgBAPAwMmoAgLMuZNQ2c9RpHEyGEKiRVbb7SX/8x5KU+/6HsSesrg0g+/ywmIzSNwAAHkZGDQBwlpHdntIOVL4J1AAAd1H6BgAAWUVGDQBwlw9q3wRqAIC7LEvfcqD0TaAGADiLJ5MBAICsIqMGADjLD6u+CdQAAHeZgN08swOBmtI3AAAeRkYNAHCWHxaTEagBAO7ywX3UlL4BAPAwMmo4zWaryl8crbS69v1f3G/VH4A9Vn0DAOB1DpSvbVD6BgDAw8ioAQDO8kPpO+mMevfu3ZozZ45KS0sVCAS0bdu2hPeNMYpEIiotLVVubq6qqqp06NChdI0XAIDPmDQ0j0s6UJ89e1aTJk3Shg0bBnz/mWee0bp167Rhwwa1trYqFApp1qxZOn36tPVgAQBIFEhD87akS981NTWqqakZ8D1jjNavX6/Vq1dr7ty5kqTNmzeruLhYW7Zs0SOPPGI3WgAAfCati8na2trU2dmp6urq+LFgMKgZM2Zoz549A/bp6elRNBpNaAAADAml7+R0dnZKkoqLixOOFxcXx9+7VGNjowoKCuItHA6nc0gAgCsZgTo1gUBizd8Y0+/YRfX19eru7o63jo6OTAwJAIC0ef7551VeXq5Ro0apsrJS77zzTsauldZAHQqFJKlf9tzV1dUvy74oGAwqPz8/oQEAMCQXt7m0aUnaunWrVq1apdWrV+v999/XHXfcoZqaGrW3t2fgD0xzoC4vL1coFFJzc3P8WG9vr1paWjRt2rR0XgoAgPjuWTYtWevWrdOSJUv08MMP65ZbbtH69esVDoe1cePG9P+BSmHV95kzZ3T06NH467a2Nh04cECFhYUaN26cVq1apTVr1mjChAmaMGGC1qxZo9GjR+uBBx5I68ABAEiXSxcyB4NBBYPBfp/r7e3V/v379cQTTyQcr66u/txF07aSDtT79u3TzJkz46/r6uokSYsWLdLPfvYzPf744/r000+1fPlyffLJJ5oyZYp27NihvLy89I0aAAApbdtcXrqQuaGhQZFIpN/HT548qb6+vqQWTdtKOlBXVVXJDFIrCAQCikQiA/6BAACkVYrzzAn9JXV0dCSskRoom/5nySyatsWzvgEAvjfUxcxFRUUaOXJkUoumbRGo4Vu2+0nf++7ylPtum/681bUBXBAwF5pN/2Tk5OSosrJSzc3Nuu++++LHm5ubdc8996Q+kEEQqAEA7krTHHUy6urqtHDhQk2ePFlTp07Vpk2b1N7ermXLllkM5PMRqAEA7krTHHUy5s+fr1OnTumpp57SiRMnVFFRoTfeeENlZWWpj2MQBGoAAJK0fPlyLV+e+vRXMgjUAAB3ZaH0PdwI1AAAd/kgUGdkUw4AAJAeZNQAAHf5IKMmUAMA3JWFVd/DjdI3AAAeRkYNAHDWcD+ZLBsI1AAAd/lgjprSNwAAHkagBgDAwyh9AwCcFZDlHHXaRpI5BGogRTZbVX6r4F9S7vtW9/9IuS9wxeH2LAAAkE1k1AAAd/lg1TeBGgDgLh8EakrfAAB4GBk1AMBZPJkMAAAvo/QNAACyiYwaAOAuH2TUBGoAgLP8MEdN6RsAAA8jowYAuMsHjxAlUAMA3MUcNQAA3sUcNQAAyCoyaiALbLaqrPnif7W69q+PPmvVH/AUSt8AAHiYZenbhUBN6RsAAA8jowYAuIvSNwAAHuaDQE3pGwAADyOjBgA4i/uoAQBAVhGoAQDwMErfAAB3+WAxGYEaAOAsP8xRE6gBAG5zINjaYI4aAAAPI6MGALiLOWoAALzLD3PUlL4BAPAwMmrAMbb7SX9z5pqU+/7vt5+0ujaQdpS+AQDwLkrfAAAgq8ioAQDuovQNAICH+SBQU/oGAMDDyKgBAM7yw2IyAjUAwF0+KH0TqAEA7vJBoGaOGgAADyOjBgA4yw9z1GTUAAB3mTS0DDl27JiWLFmi8vJy5ebm6sYbb1RDQ4N6e3uTOg8ZNQAAGfD73/9esVhMP/3pT/XFL35Rv/vd77R06VKdPXtWa9euHfJ5CNQAAGd5ufQ9e/ZszZ49O/76hhtu0OHDh7Vx40YCNQDAJ9K06jsajSYcDgaDCgaDFiceWHd3twoLC5PqQ6AGfMZmq8rbFq2zunbr5jqr/kCmhMPhhNcNDQ2KRCJpvcYf/vAH/ehHP9Jzzz2XVD8WkwEA3JWmxWQdHR3q7u6Ot/r6+s+9ZCQSUSAQGLTt27cvoc/x48c1e/Zs3X///Xr44YeT+hPJqAEAzgr8o9n0l6T8/Hzl5+cPqU9tba0WLFgw6GfGjx8f/+/jx49r5syZmjp1qjZt2pT0GAnUAAAkoaioSEVFRUP67Mcff6yZM2eqsrJSTU1NGjEi+UJ20j12796tOXPmqLS0VIFAQNu2bUt4f/Hixf1KALfffnvSAwMA4LI8fB/18ePHVVVVpXA4rLVr1+rPf/6zOjs71dnZmdR5ks6oz549q0mTJumhhx7SvHnzBvzM7Nmz1dTUFH+dk5OT7GUAALgsL9+etWPHDh09elRHjx7V2LFjE94zZugXTjpQ19TUqKamZtDPBINBhUKhZE8NAEByPLwpx+LFi7V48WLr82Rk1feuXbs0ZswY3XTTTVq6dKm6uro+97M9PT2KRqMJDQAAXJD2QF1TU6OXXnpJO3fu1HPPPafW1lbddddd6unpGfDzjY2NKigoiLdL72UDAGBQHpyfTqe0r/qeP39+/L8rKio0efJklZWVafv27Zo7d26/z9fX16uu7rOHIESjUYI1AGBIvDxHnS4Zvz2rpKREZWVlOnLkyIDvZ+oxbQAAXAkyHqhPnTqljo4OlZSUZPpSAAC/8fBisnRJOlCfOXNGR48ejb9ua2vTgQMHVFhYqMLCQkUiEc2bN08lJSU6duyYnnzySRUVFem+++5L68ABAKD0PYB9+/Zp5syZ8dcX55cXLVqkjRs36uDBg3rxxRf117/+VSUlJZo5c6a2bt2qvLy89I0aAACfSDpQV1VVDXqj9ltvvWU1IAAAhozSNwAA3kXpGwD+ie1+0tU53025747en1tdG3AVgRoA4C5K3wAAeBiBGgAA7/LDHHVGNuUAAADpQUYNAHAXpW8AALwrYIwCgzzbYyj9vY7SNwAAHkZGDQBwF6VvAAC8i1XfAAAgq8ioAQDuovQNAIB3UfoGAABZRUYNAHAXpW8ASB+brSprwo+m3PfXHf895b7wNj+UvgnUAAB3+SCjZo4aAAAPI6MGADjNhfK1DQI1AMBdxlxoNv09jtI3AAAeRkYNAHAWq74BAPAyVn0DAIBsIqMGADgrELvQbPp7HYEaAOAuSt8AACCbyKgBAM5i1TcAAF7mgweeEKgBAM4iowYAj7DZqvLu//i01bV/83++b9UfsEGgBgC4ywervgnUAABn+aH0ze1ZAAB4GBk1AMBdrPoGAMC7KH0DAICsIqMGALiLVd8AAHgXpW8AAJBVZNQAAHfFzIVm09/jCNQAAHcxRw0AgHcFZDlHnbaRZA5z1AAAZFhPT4+++tWvKhAI6MCBA0n1JVADANx18clkNm0YPP744yotLU2pL4EaAOCsi7dn2bRM+/Wvf60dO3Zo7dq1KfVnjhrAFc92P+lvLHwu5b7/9j//i9W1MTyi0WjC62AwqGAwaH3eP/3pT1q6dKm2bdum0aNHp3QOMmoAgLtMGpqkcDisgoKCeGtsbLQfmjFavHixli1bpsmTJ6d8HjJqAICzAsYoYDHPfLFvR0eH8vPz48cHy6YjkYh+8IMfDHre1tZW7dmzR9FoVPX19SmPTyJQAwCg/Pz8hEA9mNraWi1YsGDQz4wfP15PP/209u7d2y/oT548WQ8++KA2b948pOsRqAEA7or9o9n0T1JRUZGKioou+7kf/vCHevrpp+Ovjx8/rm9961vaunWrpkyZMuTrEagBAM5KV+k7E8aNG5fw+tprr5Uk3XjjjRo7duyQz8NiMgAAPIyMGgDgLoee9T1+/HiZFDJ4AjUAwF22TxcbpieT2SBQAwCcZft0seF4Mpkt5qgBAPAwMmoAgLsofQMA4F2B2IVm09/rKH0DAOBhZNQAAHdR+gYA2GxVOWvE/VbXbo79wqr/Fc+h+6hTRekbAAAPI6MGADjLy8/6ThcCNQDAXT6Yo06q9N3Y2KjbbrtNeXl5GjNmjO69914dPnw44TPGGEUiEZWWlio3N1dVVVU6dOhQWgcNAIBfJBWoW1patGLFCu3du1fNzc06f/68qqurdfbs2fhnnnnmGa1bt04bNmxQa2urQqGQZs2apdOnT6d98AAAnzP6bE/qVJr3E+rkSt9vvvlmwuumpiaNGTNG+/fv15133iljjNavX6/Vq1dr7ty5kqTNmzeruLhYW7Zs0SOPPJK+kQMAfM8Pc9RWq767u7slSYWFhZKktrY2dXZ2qrq6Ov6ZYDCoGTNmaM+ePQOeo6enR9FoNKEBADAkRp/NU6fUsv0HXF7KgdoYo7q6Ok2fPl0VFRWSpM7OTklScXFxwmeLi4vj712qsbFRBQUF8RYOh1MdEgAAV5yUA3Vtba0++OAD/fznP+/3XiAQSHhtjOl37KL6+np1d3fHW0dHR6pDAgD4jVU2bblifJikdHvWypUr9frrr2v37t0aO3Zs/HgoFJJ0IbMuKSmJH+/q6uqXZV8UDAYVDAZTGQYAwO9ikgbOA4fe3+OSyqiNMaqtrdUrr7yinTt3qry8POH98vJyhUIhNTc3x4/19vaqpaVF06ZNS8+IAQDwkaQy6hUrVmjLli167bXXlJeXF593LigoUG5urgKBgFatWqU1a9ZowoQJmjBhgtasWaPRo0frgQceyMgfAADwLz+s+k4qUG/cuFGSVFVVlXC8qalJixcvliQ9/vjj+vTTT7V8+XJ98sknmjJlinbs2KG8vLy0DBgAgDgfPJksqUBthvAHBQIBRSIRRSKRVMcEAAD+gWd9AwDcRUYNALBhu5/07JIVKfd988SPra7tBB8EavajBgDAw8ioAQDu8sF91ARqAICzuD0LAAAvY44aAABkExk1AMBdMSMFLLLimPczagI1AMBdlL4BAEA2kVEDABxmu6e09zNqAjUAwF2UvgEAQDaRUQMA3BUzsipfs+obAIAMMrELzaa/x1H6BgDAw8ioAcDDbLaqnHX7U1bXbt7736z6DwsfLCYjUAMA3MUcNQAAHuaDjJo5agAAPIyMGgDgLiPLjDptI8kYAjUAwF2UvgEAQDaRUQMA3BWLSbJ4aEnM+w88IVADANxF6RsAAGQTGTUAwF0+yKgJ1AAAd/ngyWSUvgEA8DAyagCAs4yJyVhsVWnTd7iQUQMA3GXMhfJ1qm0Y5qi3b9+uKVOmKDc3V0VFRZo7d25S/cmoAQDuMpZz1BkO1L/85S+1dOlSrVmzRnfddZeMMTp48GBS5yBQA8AVynY/6cn/si6lfn2956yue6U4f/68Hn30UT377LNasmRJ/PiXvvSlpM5D6RsA4K5YzL5JikajCa2np8d6aO+9954+/vhjjRgxQl/72tdUUlKimpoaHTp0KKnzEKgBAO66eB+1TZMUDodVUFAQb42NjdZD++ijjyRJkUhE3//+9/WrX/1K1113nWbMmKG//OUvQz4PgRoA4HsdHR3q7u6Ot/r6+s/9bCQSUSAQGLTt27dPsX9k66tXr9a8efNUWVmppqYmBQIB/eIXvxjy2JijBgA4y8RiMgH727Py8/OVn58/pD61tbVasGDBoJ8ZP368Tp8+LUn68pe/HD8eDAZ1ww03qL29fchjJFADANyVhVXfRUVFKioquuznKisrFQwGdfjwYU2fPl2S9Pe//13Hjh1TWVnZkK9HoAYAIAPy8/O1bNkyNTQ0KBwOq6ysTM8++6wk6f777x/yeQjUAAB3xYwU8O591M8++6yuuuoqLVy4UJ9++qmmTJminTt36rrrrhvyOQjUAAB3GSPJ4jGgGQ7UV199tdauXau1a9emfA5WfQMA4GFk1AAAZ5mYkbEofRv2owYAIINMTHalb+/vnkWgBgA4yw8ZNXPUAAB4mOcy6ov/uolGo1keCQD4W6q7YPX9/UK/4chWz5seq/L1ef09jaPJDM8F6ouPXAuHw1keCQDAxunTp1VQUJCRc+fk5CgUCundzjeszxUKhZSTk5OGUWVGwHisQB+LxXT8+HHl5eUpEAj0ez8ajSocDqujo2PIz2X1O76z5PGdJY/vLHlX6ndmjNHp06dVWlqqESMyN8N67tw59fb2Wp8nJydHo0aNSsOIMsNzGfWIESM0duzYy34umQeo4wK+s+TxnSWP7yx5V+J3lqlM+p+NGjXK0wE2XVhMBgCAhxGoAQDwMOcCdTAYVENDg4LBYLaH4gy+s+TxnSWP7yx5fGcYCs8tJgMAAJ9xLqMGAMBPCNQAAHgYgRoAAA8jUAMA4GHOBernn39e5eXlGjVqlCorK/XOO+9ke0ieFYlEFAgEElooFMr2sDxl9+7dmjNnjkpLSxUIBLRt27aE940xikQiKi0tVW5urqqqqnTo0KHsDNYjLvedLV68uN/v7vbbb8/OYD2gsbFRt912m/Ly8jRmzBjde++9Onz4cMJn+J1hME4F6q1bt2rVqlVavXq13n//fd1xxx2qqalRe3t7tofmWbfeeqtOnDgRbwcPHsz2kDzl7NmzmjRpkjZs2DDg+88884zWrVunDRs2qLW1VaFQSLNmzYo/k96PLvedSdLs2bMTfndvvGH/PGZXtbS0aMWKFdq7d6+am5t1/vx5VVdX6+zZs/HP8DvDoIxDvvGNb5hly5YlHLv55pvNE088kaUReVtDQ4OZNGlStofhDEnm1Vdfjb+OxWImFAqZf/3Xf40fO3funCkoKDA/+clPsjBC77n0OzPGmEWLFpl77rknK+NxQVdXl5FkWlpajDH8znB5zmTUvb292r9/v6qrqxOOV1dXa8+ePVkalfcdOXJEpaWlKi8v14IFC/TRRx9le0jOaGtrU2dnZ8JvLhgMasaMGfzmLmPXrl0aM2aMbrrpJi1dulRdXV3ZHpJndHd3S5IKCwsl8TvD5TkTqE+ePKm+vj4VFxcnHC8uLlZnZ2eWRuVtU6ZM0Ysvvqi33npLL7zwgjo7OzVt2jSdOnUq20NzwsXfFb+55NTU1Oill17Szp079dxzz6m1tVV33XWXenp6sj20rDPGqK6uTtOnT1dFRYUkfme4PM/tnnU5l259aYwZcDtMXPg/zIsmTpyoqVOn6sYbb9TmzZtVV1eXxZG5hd9ccubPnx//74qKCk2ePFllZWXavn275s6dm8WRZV9tba0++OADvfvuu/3e43eGz+NMRl1UVKSRI0f2+xdmV1dXv3+JYmDXXHONJk6cqCNHjmR7KE64uEKe35ydkpISlZWV+f53t3LlSr3++ut6++23E7by5XeGy3EmUOfk5KiyslLNzc0Jx5ubmzVt2rQsjcotPT09+vDDD1VSUpLtoTihvLxcoVAo4TfX29urlpYWfnNJOHXqlDo6Onz7uzPGqLa2Vq+88op27typ8vLyhPf5neFynCp919XVaeHChZo8ebKmTp2qTZs2qb29XcuWLcv20Dzpscce05w5czRu3Dh1dXXp6aefVjQa1aJFi7I9NM84c+aMjh49Gn/d1tamAwcOqLCwUOPGjdOqVau0Zs0aTZgwQRMmTNCaNWs0evRoPfDAA1kcdXYN9p0VFhYqEolo3rx5Kikp0bFjx/Tkk0+qqKhI9913XxZHnT0rVqzQli1b9NprrykvLy+eORcUFCg3N1eBQIDfGQaX1TXnKfjxj39sysrKTE5Ojvn6178ev8UB/c2fP9+UlJSYq6++2pSWlpq5c+eaQ4cOZXtYnvL2228bSf3aokWLjDEXbp1paGgwoVDIBINBc+edd5qDBw9md9BZNth39re//c1UV1eb66+/3lx99dVm3LhxZtGiRaa9vT3bw86agb4rSaapqSn+GX5nGAzbXAIA4GHOzFEDAOBHBGoAADyMQA0AgIcRqAEA8DACNQAAHkagBgDAwwjUAAB4GIEaAAAPI1ADAOBhBGoAADyMQA0AgIcRqAEA8LD/DyGIzDbGagX4AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "plt.imshow(np.log10(s2.covariance.covmat))\n", "plt.colorbar()" @@ -397,7 +477,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "id": "6226ebf5-5ced-4886-bbed-c04284043f75", "metadata": {}, "outputs": [], @@ -407,7 +487,7 @@ "bin_richness_labels = []\n", "\n", "survey_name = \"NC_mock_redshift_richness\"\n", - "s_count2.add_tracer(\"cluster_survey\", survey_name, area)\n", + "s_count2.add_tracer(\"survey\", survey_name, area)\n", "\n", "for i, z_bin in enumerate(zip(z_edges[:-1], z_edges[1:])):\n", " lower, upper = z_bin\n", @@ -441,7 +521,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "id": "cd235d47-bbc4-4b74-9276-56aa5fac14b0", "metadata": {}, "outputs": [], @@ -480,7 +560,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "id": "c912d717-eb26-4072-b9ba-281befa7ff9b", "metadata": {}, "outputs": [], @@ -510,7 +590,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "id": "433d8c0e-4c4c-4214-a72d-0b011f2541fc", "metadata": {}, "outputs": [], @@ -544,10 +624,104 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "id": "04db6a66-6b33-45d7-8a7f-5ab00890523c", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0'), value=74639, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1'), value=26398, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2'), value=3158, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3'), value=105, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0'), value=135594, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1'), value=45119, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2'), value=4894, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3'), value=162, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0'), value=185542, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1'), value=58654, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2'), value=5866, ),\n", + " DataPoint(data_type='cluster_counts', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3'), value=113, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0', 'bin_radius_0'), value=29050865783927.395, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0', 'bin_radius_1'), value=18588353432314.89, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0', 'bin_radius_2'), value=10796854943136.59, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0', 'bin_radius_3'), value=5772635608161.832, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0', 'bin_radius_4'), value=2886119336410.716, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_0', 'bin_radius_5'), value=1368962867902.2593, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1', 'bin_radius_0'), value=41233092797499.516, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1', 'bin_radius_1'), value=27505847912994.324, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1', 'bin_radius_2'), value=16589396009095.674, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1', 'bin_radius_3'), value=9154038784864.775, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1', 'bin_radius_4'), value=4693620704875.559, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_1', 'bin_radius_5'), value=2270526335611.045, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2', 'bin_radius_0'), value=66416557623257.34, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2', 'bin_radius_1'), value=47096483269269.5, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2', 'bin_radius_2'), value=30097538584127.676, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2', 'bin_radius_3'), value=17462598870593.883, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2', 'bin_radius_4'), value=9328022913001.797, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_2', 'bin_radius_5'), value=4660315570393.327, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3', 'bin_radius_0'), value=106986213336953.16, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3', 'bin_radius_1'), value=80919708596517.55, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3', 'bin_radius_2'), value=55202529359557.2, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3', 'bin_radius_3'), value=33994259233850.05, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3', 'bin_radius_4'), value=19094956517762.902, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_0', 'bin_rich_3', 'bin_radius_5'), value=9932690237459.707, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0', 'bin_radius_0'), value=32587950929002.02, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0', 'bin_radius_1'), value=20379038976194.39, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0', 'bin_radius_2'), value=11600721702709.613, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0', 'bin_radius_3'), value=6100334051829.585, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0', 'bin_radius_4'), value=3010165023263.342, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_0', 'bin_radius_5'), value=1413317712470.7285, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1', 'bin_radius_0'), value=45458814464648.16, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1', 'bin_radius_1'), value=29521288295870.605, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1', 'bin_radius_2'), value=17375548318583.297, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1', 'bin_radius_3'), value=9392796640628.133, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1', 'bin_radius_4'), value=4737377269273.788, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_1', 'bin_radius_5'), value=2262425615341.657, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2', 'bin_radius_0'), value=73842163412903.06, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2', 'bin_radius_1'), value=50910928102274.35, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2', 'bin_radius_2'), value=31666725548146.758, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2', 'bin_radius_3'), value=17942129951519.47, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2', 'bin_radius_4'), value=9399024516477.635, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_2', 'bin_radius_5'), value=4623929348356.979, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3', 'bin_radius_0'), value=112247981630090.03, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3', 'bin_radius_1'), value=81914840579736.38, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3', 'bin_radius_2'), value=53861228335403.6, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3', 'bin_radius_3'), value=32059337825829.258, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3', 'bin_radius_4'), value=17494323338239.246, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_1', 'bin_rich_3', 'bin_radius_5'), value=8890086692097.94, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0', 'bin_radius_0'), value=35936619495183.12, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0', 'bin_radius_1'), value=22007299985791.734, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0', 'bin_radius_2'), value=12302736751882.656, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0', 'bin_radius_3'), value=6374754560854.587, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0', 'bin_radius_4'), value=3109279660830.6567, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_0', 'bin_radius_5'), value=1446785207664.3884, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1', 'bin_radius_0'), value=49526225007016.42, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1', 'bin_radius_1'), value=31404690109286.375, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1', 'bin_radius_2'), value=18095981017891.367, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1', 'bin_radius_3'), value=9611537780060.18, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1', 'bin_radius_4'), value=4780362356676.691, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_1', 'bin_radius_5'), value=2258254452608.0923, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2', 'bin_radius_0'), value=76891264521828.5, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2', 'bin_radius_1'), value=51353282953337.69, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2', 'bin_radius_2'), value=31005956516000.42, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2', 'bin_radius_3'), value=17124881536532.178, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2', 'bin_radius_4'), value=8787092289236.374, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_2', 'bin_radius_5'), value=4253201308523.909, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3', 'bin_radius_0'), value=111603221923155.64, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3', 'bin_radius_1'), value=78083691672770.1, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3', 'bin_radius_2'), value=49254937057929.83, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3', 'bin_radius_3'), value=28251504646716.047, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3', 'bin_radius_4'), value=14948905345881.504, ),\n", + " DataPoint(data_type='cluster_shear', tracers=('NC_mock_redshift_richness', 'bin_z_2', 'bin_rich_3', 'bin_radius_5'), value=7412718656160.18, )]" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "s_count2.data" ] @@ -562,7 +736,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "id": "ff31f15c-6808-4870-bc92-665d6f1299a2", "metadata": {}, "outputs": [], @@ -599,10 +773,31 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "id": "209ca0b4-a743-40da-a956-1688291682cc", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAG1CAYAAAAC+gv1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABffUlEQVR4nO3dd3hU1drG4d/MpEEaJQk1QOiE3nvvoKIePSqgIIiiVBHx82Dl2A5KEQyeozQ7diwgSO+991BCJ4SahJI2M98fmyDdJGSyZ5Lnvq65cK+Z2fMikDx519prW5xOpxMRERERucpqdgEiIiIi7kYBSUREROQGCkgiIiIiN1BAEhEREbmBApKIiIjIDRSQRERERG6ggCQiIiJyAwUkERERkRt4mV2AJ3I4HBw/fpzAwEAsFovZ5YiIiEgGOJ1OEhMTKV68OFbrnXtECkhZcPz4ccLDw80uQ0RERLLgyJEjlCxZ8o6vUUDKgsDAQMD4HxwUFGRyNSIiIpIRCQkJhIeHX/0+ficKSFmQPq0WFBSkgCQiIuJhMrI8Rou0RURERG6ggCQiIiJyA02xiYhInmW320lNTTW7DMkm3t7e2Gy2bDmXApKIiOQ5TqeT2NhYzp8/b3Ypks0KFChA0aJF73obHgUkERHJc9LDUVhYGPnz59eedrmA0+nk0qVLxMXFAVCsWLG7Op8CkoiI5Cl2u/1qOCpcuLDZ5Ug2ypcvHwBxcXGEhYXd1XSbFmmLiEiekr7mKH/+/CZXIq6Q/ud6t2vLFJAyISoqisjISOrXr292KSIicpc0rZY7ZdefqwJSJgwYMICdO3eybt06s0sRERERF1JAcgeL3oUlo2/93JLRxvMiIiKSYxSQ3IHVBovevjkkLRltjFuzZ08HERHxbK1atWLo0KG3fb5MmTKMHz/epTW88cYb1KpVK8Ovt1gszJw502X1uIoCkjtoOQJajzTC0OL/GGPp4aj1SON5ERFxC+PmRTNhwd5bPjdhwV7GzYvO4Yr+sm7dOp5++mmXfsbw4cNZsGCBSz/jWocPH+bee+/F39+fkJAQBg8eTEpKiss/V5f5u4uWI+DyOVj8Diz5DzjtCkciIm7IZrUw9koIGty2wtXxCQv2MnZeNMPaVzSrNEJDQ112bqfTid1uJyAggICAAJd9zrXsdjtdu3YlNDSU5cuXc+bMGXr16oXT6WTixIku/Wx1kNyJf4jxq9Nu/Lr7d1g5EeKPmVeTiEgu53Q6uZSSluHHU80jGNSmPGPnRTPmzz1cSkljzJ97GDsvmkFtyvNU84gMn8vpdGa63rS0NAYOHEiBAgUoXLgwr7zyytXz3DjFZrFYmDx5Mg888AD58+enQoUK/Prrrxn6nMWLF2OxWJg7dy716tXD19eXZcuW3XKKberUqVStWhVfX1+KFSvGwIEDr3v+9OnTWarhzz//ZOfOnXz55ZfUrl2bdu3aMWbMGD799FMSEhIydI6sUgfJnaRcMn61WMDphBNbjMefr0LpplD9IYjsBvkLmVuniEgucjnVTuRrc7P03okL9zFx4b7bHv+dnaM6kt8nc9+KP/vsM/r27cuaNWtYv349Tz/9NKVLl6Zfv363fP2bb77J6NGjef/995k4cSI9evTg0KFDFCqUse8lI0aM4IMPPqBs2bIUKFCAJUuWXPf8xx9/zLBhw3jvvffo3Lkz8fHxrFixIltqWLVqFdWqVaN48eJXxzp27EhycjIbNmygdevWGfo9ZIU6SO5iyWhY9oExrfb6eWj6vDEeHA444dBy+H0ofFABvn4Etv0AKRdNLFhERMwQHh7OuHHjqFSpEj169GDQoEGMGzfutq/v3bs3jz32GOXLl+edd97h4sWLrF27NsOfN2rUKNq3b0+5cuVuufP4W2+9xQsvvMCQIUOoWLEi9evXv2kheVZriI2NpUiRIteNFSxYEB8fH2JjYzP8e8gKdZDcwa0WZLd/A3zyG+ONB0JAGGz7HmK3QfQc4+GdHyp1geoPQ7k24OVj6m9DRMQT5fO2sXNUx0y/7+PF+5m4cB/eNgupdieD2pTn2VblMv3ZmdWoUaPrNkNs3LgxY8aMwW633/L1NWrUuPrf/v7+BAYGXr1fWUbUq1fvts/FxcVx/Phx2rZte8dz3E0Nt9r40el0unyjTwUkd+C4zYLs9GOHHZoOMR6n9hjdo23fw7kY2P6D8chX0Jh+q/4wlGoCVjUHRUQywmKxZHqaa8KCvUxcuI9h7SsyuG2Fqwu0vW3W6xZuuwNvb+/rji0WCw6HI8Pv9/f3v+1z6fc+c1UNRYsWZc2aNdeNnTt3jtTU1Js6S9lNASkToqKiiIqKum1Kz7LWL9/+uRtDU2glaDMSWv8Ljm28EpB+hAsnYcN04xFYHKo9aISlYjWNNU0iIpItrr1aLT0Mpf96q6vbstvq1atvOq5QocJd3Zg1qwIDAylTpgwLFixwyXqgxo0b8/bbb3PixAmKFSsGGAu3fX19qVu3brZ/3rUUkDJhwIABDBgwgISEBIKDg80txmKBknWNR4e34OAyo7O081dIPA6rPjIehcsbQanaQxBS3tyaRURyAbvDeV04Spd+bHdk/sq0zDhy5AjDhg3jmWeeYePGjUycOJExY8a49DPv5I033qB///6EhYXRuXNnEhMTWbFiBYMGDbrrc3fo0IHIyEgef/xx3n//fc6ePcvw4cPp168fQUFB2VD97Skg5QZWG5RtZTy6joG984zO0p4/4Mw+WPyu8ShW60pYehCCiv/NSUVE5Faev8M+RzkxvfbEE09w+fJlGjRogM1mY9CgQS7fHPJOevXqRVJSEuPGjWP48OGEhITw0EMPZcu5bTYbs2bN4rnnnqNp06bky5eP7t2788EHH2TL+e/E4szKJgx5XHoHKT4+3uUJ9q4kJ8Lu2cZ6pf0L/9pfCQuUaWZsG1DlPm0bICJ5SlJSEjExMURERODn52d2OZLN7vTnm5nv3+og5Wa+gVDzEeNx8TTsnGlMwx1eZUzJHVwGs4ZD+XZGWKrUGXxuvxhPREQkr9ClTnmFfwjUfwr6zIGh26Ddm1CkOjhSIfoP+LEvvF8BfnwKoueCPdXsikVExIX69+9/9bYhNz769++fZ2q4HU2xZYHHTLFlRNxuY73Stu/h3MG/xvMVhMj7r2wb0FjbBohIrqEpNkNcXNxtb9cRFBREWFiYR9aQXVNsCkhZkKsCUjqnE45tMKbgtv8IF6/ZwCuoxF/bBhStoW0DRMSjKSDlblqDJNnLYoGS9YxHx7chZqnRWdr5GyQcM26au3IiFK5gBKXqD0HhzO0YKyIi4ik0byI3s9qgXGvoFgXDo+GRL43pNi8/OLMXFr8DE+vAJ61gVRQknDC7YhERkWylDpLcmbcfVLnXeCQlwO5ZRmdp/yI4vsl4zB15ZduAhyHyPmP9koiIiAdTQJKM8wuCWo8Zjwun/to24Mjqa7YNeAEqtDem4Cp2Nm64KyIi4mEUkCRrAkKhQT/jce6QsbB7+49wcjvsmW08vP2hclejs1SuNdi8//68IiIibkBrkDIhKiqKyMhI6tevb3Yp7qVgaWg+DJ5dAc+uguYvQIHSkHoRtn0HXz8MH1SE35+HQyshE3eRFhGRv7Rq1YqhQ4fe9vkyZcowfvx4l9bwxhtvUKtWrQy/3mKxMHPmTJfV4yoKSJkwYMAAdu7cybp168wuxX0ViYS2r8GQLdB3PjR4BvxD4fJZWD8VpnWG8dXhz1fhxFZjewEREU+y6F1YMvrWzy0ZbTxvknXr1rn8vmzDhw9nwYIFLv2MdFu2bOGxxx4jPDycfPnyUaVKFT788MMc+WxNsYlrWCwQXt94dHwHDi411ivt+g0SjsLKCcYjpOKVG+j+Q9sGiIhnsNpg0dvGf7cc8df4ktHGeOuR5tQFhIaGuuzcTqcTu91+dafrnLBhwwZCQ0P58ssvCQ8PZ+XKlTz99NPYbDYGDhzo0s9WB0lcz+YF5drA/ZNg+F745xfGTXJtvnA62viCMrEOfNIaVk2CxFizKxaRvMTphJSLGX80HgAtXjS+di18yxhb+JZx3OJF4/mMnisLXfS0tDQGDhxIgQIFKFy4MK+88grpez7fOMVmsViYPHkyDzzwAPnz56dChQr8+uuvGfqcxYsXY7FYmDt3LvXq1cPX15dly5bdcopt6tSpVK1aFV9fX4oVK3ZTeDl9+nSWaujTpw8TJkygZcuWlC1blp49e/Lkk0/y008/Zej9d0MdJMlZ3n7GVgCR90FSvLFtwLbv4cBiOL7ReMz9F0Q0NzpLVe7VtgEi4lqpl+Cd4ll779L3jcftjv/Ov45n+ibhn332GX379mXNmjWsX7+ep59+mtKlS9OvX79bvv7NN99k9OjRvP/++0ycOJEePXpw6NAhChUqlKHPGzFiBB988AFly5alQIECLFmy5LrnP/74Y4YNG8Z7771H586diY+PZ8WKFdlaw7Xi4+Oz9L7MUkAS8/gFQ63uxuNCHOyYaYSlo2uNnbxjlhrbBpRP3zagk7YNEJE8Lzw8nHHjxmGxWKhUqRLbtm1j3Lhxtw1IvXv35rHHHgPgnXfeYeLEiaxdu5ZOnTpl6PNGjRpF+/btb/v8W2+9xQsvvMCQIUOujt14MdPd1pBu1apVfPfdd8yaNStT78sKBSRxDwFh0PBp43HuoLFlwLYfIW4H7JllPHwC/to2oGwrbRsgItnDO7/Rycms5eOMbpHNB+wpxvRas+cz/9mZ1KhRIyzX3BOzcePGjBkzBrvdfsvX16hR4+p/+/v7ExgYSFxc3C1feyv16tW77XNxcXEcP36ctm3b3vEcd1sDwI4dO+jWrRuvvfbaHQNbdlFAEvdTsIyxVUDzF+Dkjis30P0Bzh+Grd8aj/yFjdufVH8YwhuCVcvpRCSLLJZMT3OxZLQRjlqPNBZqpy/Qtvlcv3DbDXh7X//DpMViwZGJ7Vb8/W///yZfvnw5UsPOnTtp06YN/fr145VXXsnw++6GApK4tyJVjUfb1+DoOmMKbsfPcPEUrJ9iPIJKQvV/GGGpSDXji52IiKtce7VaehhK//VWV7dls9WrV990XKFCBWw2m8s+83YCAwMpU6YMCxYsoHXr1i75jB07dtCmTRt69erF22+/7ZLPuBUFJPEMFguENzAeHd+FmCXXbxuw4kPjEVLJCErV/wGFyppdtYjkRg779eEoXfqx49ZTXdnlyJEjDBs2jGeeeYaNGzcyceJExowZ49LPvJM33niD/v37ExYWRufOnUlMTGTFihUMGjTors+9Y8cOWrduTYcOHRg2bBixscZVzjabzaVbGoACkngimxeUb2s87hkLe/80OkvRf8LpPbDoLeNRoq4Rlqo+AIFFjc3brLZb/2S3ZPSVL3ov5/zvR0Q8y52+TuTA9NoTTzzB5cuXadCgATabjUGDBrl8c8g76dWrF0lJSYwbN47hw4cTEhLCQw89lC3n/v777zl16hRfffUVX3311dXx0qVLc/DgwWz5jNuxOJ3ayjizEhISCA4OJj4+nqCgILPLkXRJ8bDrdyMsxSwB55X5bYsVyjQ3FnnvmXXzT363apeLSK6VlJRETEwMERER+Pn5mV2OZLM7/flm5vu3OkiSe/gFQ+0exuNCnLFWadv3xtqlmCv7dliu7IAbtxPu/xhWTlQ4EhGRm+jSH8mdAsKg4TPw1HwYvBnavAqhVcB5ZW3Ajp/h7aJGOGoyROFIRPKc/v37X71tyI2P/v3755kabkdTbFmgKTYPdnKH0VVaPu6vMZsPVP8nNH7OuGJORHI1TbEZ4uLiSEhIuOVzQUFBhIWFeWQNmmIzQVRUFFFRUbfdjEs8QJGqxu1NAKxe4EgzNnjb/KXxKNfGuI9SubbaLkBEcrWwsLAcCUHuXsPtaIotEwYMGMDOnTtZt26d2aVIVl27IPu1M3/ddTu0srGYe/9C+PIfMKkxbPwC0pLNrVdEXEYTKLlTdv25KiBJ3nG7zd1aj4RTu6Fhf2j4rHG126ld8OtAGFfNeN/FM+bWLiLZJn1X50uXLplcibhC+p/rjbt3Z5bWIGWB1iB5qIzug3T5PGz8DNb8DxKOGc97+UHNx4zpt5AKOVq2iGS/EydOcP78ecLCwsifP/919zYTz+R0Orl06RJxcXEUKFCAYsWK3fSazHz/VkDKAgWkPMKeCjt/MbYCOLH5r/GKnYygVKa51imJeCin00lsbCznz583uxTJZgUKFKBo0aK3DL0KSC6mgJTHOJ1waCWs+gj2/AFc+SdTtAY0HgjVHgTb3bVyRcQcdrud1NRUs8uQbOLt7X3He9IpILmYAlIednofrJ4Em7+GtMvGWGBxaPg01O0N+QqaWp6IiNyeApKLKSAJl87C+imw9lO4cNIY8/aH2j2h0bNQKMLc+kRE5CYKSC6mgCRXpSXDth9gVRTE7TDGLFao3BUaD4LwBlqnJCLiJhSQXEwBSW7idMKBRUZQ2jf/r/ES9aDJQKh8L9i0L6uIiJkUkFxMAUnuKG6XEZS2fmvs0g1QoJSxx1LtnuCnvzMiImZQQHIxBSTJkAtxsG6y8bh0ZaNJ3yCo84SxKWWBcHPrExHJYxSQXEwBSTIl9TJsmWF0lc7sNcYsNqj6gLGfUok65tYnIpJHKCC5mAKSZInDAfvmGRtPHlz213ipJsY6pYqdjJ2+RUTEJRSQXEwBSe7aiS2wahJs/wEcacZYobLQ6Dmo1R18/M2tT0QkF1JAcjEFJMk2Ccdh7SewfiokxRtj+QpCvT7Q4GkILGpufSIiuYgCkospIEm2S75g7M69OgrOHTTGrN5Q/WFjnVLRaqaWJyKSGygguZgCkriMww57ZsPKj+DI6r/Gy7Yy7vtWvp02nhQRySIFJBdTQJIccXS9cYPcnb+A02GMhVY21inVeAS8/cytT0TEwygguUhUVBRRUVHY7Xaio6MVkCRnnDsEa/4HGz+HlERjLH8INOgH9Z8C/xBz6xMR8RAKSC6mDpKYIineCEmr/wsJR40xmy/UfNRYpxRaydz6RETcnAKSiykgiansabBzpjH9dnzTX+MVOhjrlCJaaJ2SiMgtKCC5mAKSuAWnEw6vNoLS7lnAlX/KRasbQanqg+DlY2qJIiLuRAHJxRSQxO2c2Q+rP4bNX0HqJWMssJixl1Ld3pC/kKnliYi4AwUkF1NAErd16SxsmAZrPoELscaYd36o1QMaPQuFy5lbn4iIiRSQXEwBSdxeWgps/9GYfju5/cqgBSp3NabfSjXSOiURyXMUkFxMAUk8htMJMUtgVRTs/fOv8eJ1jBvkVukGNi/z6hMRyUEKSC6mgCQeKW43rJ4EW2aAPdkYCw6Hhv2hzuPgF2xufSIiLqaA5GIKSOLRLpyC9VNg7adw6bQx5hMIdXtBw2egQClz6xMRcREFJBdTQJJcITUJtn5rTL+d3mOMWWwQ2c1Yp1Syrrn1iYhkMwUkF1NAklzF4YD9C4wF3QcW/zVeqrGxQ3elLmC1mVaeiEh2UUByMQUkybVit8GqSbDte3CkGmMFI4wb5NbqDr4B5tYnInIXFJBcTAFJcr2EE7DuU1g3BZLOG2N+BaDek9DgGQgqZmZ1IiJZooDkYtkdkMbNi8ZmtTC4bYWbnpuwYC92h5Pn21e8688RybSUi7D5a+Pqt7MHjDGrN1T7hzH9VqyGufWJiGRCZr5/W3OoJrkDm9XC2HnRTFiw97rxCQv2MvZKeBIxhY8/NOgHA9fDo19DqSbG1NvWGfC/5vDZvRA911jHJCKSi6iDlAWumGJLD0PtqoQx/tHaTF0ew9h50QxrX/GWnSUR0xzbYFz5tmMmOO3GWEhFY51SzUdh+XhjUXfLETe/d8locNih9cs5WbGICKApNpdz1Rqkwd9s4tctx68eP9Usglfuicy284tkq/NHYM1/YePnkJxgjOUvDGFV4OByaD3y+pC0ZDQsevvmcRGRHKIpNg91X83i1x1PW3mQAV9tZMOhcyZVJHIHBcKh49vw/A7o+A4El4JLZ4xwZLEZYWj2i8ZrFY5ExMOog5QFruogpU+zeVktpDmu/2OpXaoAfZtF0KlqUbxsyrXihuxpsPs3WPkRHFv/17jFYtwTTuFIREymDpIHSg9Hw9pXZN87XRh25aq1qsWD8LFZ2XT4PAO/3kTL9xfzydL9xF9ONblikRvYvKDqA/DUfOjzJ1S51xhP/xnsdDScjTGvPhGRTFAHKQuyu4N0bTi6dkF2+vgzLcri523jy9WHOHMxBQB/HxsP1wvnyaZlKF3Y/65rEMl26dNqFis4r1zlZvU29lJq8SIEhJlbn4jkOVqk7WJm7YOUlGrnl83HmLI8huiTFwBj9qJ9lSL0bRZBg4hCWCzaEkDcwI1rjn4fZtwgN523PzQZaNzzzU+brYpIzlBAcjGzd9J2Op0s33eaKctjWLzn1NXxaiWC6Nssgq7Vi+PjpdlTMcntFmSnjwcWg8QTxlj+wtB8ONTvC16+5tQrInmGApKLmR2QrrUvLpGpKw7y44ajJKcZ0xhFgnx5onEZejQsRYH8PqbWJ3nQonf/Zh+kNChSFRaMgjP7jPHgUtD6X1Djn7oxroi4jAKSi0RFRREVFYXdbic6OtotAlK6sxdT+HrNIT5fdYi4xGQA/Lyt/KNOSfo0i6BcqG4yKm7Gngabv4TF7/3VUQqrCu1ehwodjPljEZFspIDkYu7UQbpRSpqD37ceZ8ryGHYcT7g63qZyGH2bRdCkXGGtUxL3knIJ1v4Plo+DpHhjrFQTaPcGlGpoamkikrsoILmYOwekdE6nk9UHzjJleQwLdp+8eqV15aKB9GkWQbdaxfH10lSGuJHL54yQtOZ/kJZkjFXqAm1fM3bnFhG5SwpILuYJAelaMacvMn1FDN9vOMqlFOPeWSEBPjzeqAw9GpUiJECLY8WNxB+DJe/Bpi+N7QEsVqj5GLR62di9W0QkixSQXMzTAlK6+EupfLPuMJ+tPMiJeOMndB8vKw/UKkGfZhFUKhpocoUi1zgVDQtHwa7fjGObLzToB81fgPyFzK1NRDySApKLeWpASpdqd/DH9limLI9hy5HzV8ebVwihT7MIWlYIxWrVOiVxE0fXw/w34OAy49g3CJoOhkbPgY82SRWRjFNAcjFPD0jpnE4nGw+fY8ryGOZsjyX99m/lwwLo0zSCB+uUwM9b65TEDTidsG+BEZRObjPGAooYWwnU6QU2b1PLExHPoIDkYrklIF3ryNlLTF95kG/XHeFCchoABfN706NhaZ5oXJqwID+TKxQBHA7Y/iMsegvOHTTGCpWFNq9A5ANg1QapInJ7CkgulhsDUrrEpFS+W3+UaStiOHruMgDeNgv31ixO32YRVC0ebHKFIkBaCmyYDktHw8Uru8kXq2XsoVSujZmViYgbU0BysdwckNLZHU7+3GGsU1p/6NzV8UZlC9G3WVnaVg7TOiUxX/IFWBUFKydAinF/QiJaGnsolahjamki4n4UkFwsLwSka205cp4py2OYte0E9isLlcoUzk+fZhH8o05J/H29TK5Q8ryLp2HpB8YNce0pxljk/dDmVQgpb2ppIuI+FJBcLK8FpHTHz1/ms1UH+WbNYRKSjHVKQX5ePNawFL0al6F4gXwmVyh53rlDsOgd2Pot4ASLDeo8Aa3+DwKLml2diJhMAcnF8mpASncxOY0fNx5l6vIYDp65BIDNaqFL9WL0bRZBrfAC5hYoErsdFv4boucYx175oNGz0HQI5CtgamkiYh4FJBfL6wEpncPhZOHuOKYsj2HVgTNXx+uVLkjfZhF0qFoUm9YpiZkOrYR5r8PRtcaxXwFjo8kGT4O3rswUyWsUkFxMAelmO47HM2V5DL9tOU6q3fgrVbJgPno3KcMj9cMJ9NM+NWISpxP2/AEL3oRTu42xoBLGrUtqPgY2raETySsUkFxMAen24hKS+GL1Ib5cfYhzl1IBCPD14pH64fRuUobwQvlNrlDyLIcdtsww1iglHDXGQioZN8Ot3BUs6naK5HYKSC6mgPT3klLt/LTxGFNXxLAvzrj82mqBjlWL0rdZBHVLF8Sib0hihtQkWDcZln0Al69sYVGygbE1QJmmppYmIq6lgORiCkgZ53A4Wbr3FFOWx7Bs7+mr4zVLBtO3eVk6VyuKt027H4sJkuJhxQRYPQlSjYsNqNAB2r4ORauZW5uIuIQCkospIGXNnthEpi6P4efNx0hJcwBQLNiPXk3K8Fj9UgTn1zolMUFiLCwZbezM7bQDFqjxT2j9LyhYxuTiRCQ7KSC5mALS3Tl9IZmvVh/mi9UHOX3B2NQvv4+Nh+qW5MmmEUSE6A7tYoIz+2HhW7DjJ+PY6g31+0Lz4RAQam5tIpItFJBcTAEpeySn2fl183GmLI9hd2wiYKyTbVu5CH2bRdCobCGtU5Kcd3wTzH8TDiwyjn0CoMkgaDwAfAPNrU1E7ooCkospIGUvp9PJyv1nmLI8hoW7466ORxYLom+zCO6tWRwfL61Tkhx2YDHMf8MITAD5Q6DlCKj7JHj5mFmZiGSRApKLKSC5zv5TF5i2IoYfNhwlKdVYpxQa6EuvxqXp3rA0hfz1jUlykNMJO2fCgn/D2f3GWIHS0OYVqPYQWBXcRTyJApKLKSC53rmLKXy99jCfrzrIyYRkAHy9rDxYpyR9m5WhfJimOiQH2VNh0xew+D9wIdYYK1Id2r0O5dtpDyURD6GA5GIKSDknJc3B7G0nmLI8hm3H4q+Ot6wYSt9mETSvEKJ1SpJzUi7Bmo9h+YeQfOXvY+lmxh5K4fVNLU1E/p4CkospIOU8p9PJuoPnmLzsAPN2nST9b23FIgH0aRrB/bVL4OdtM7dIyTsunYXlY2HNJ2A3OpxUvsfYlTu0krm1ichtKSC5mAKSuQ6duci0FQf5fv0RLqbYASjs70OPRqV5vFFpQgN9Ta5Q8oz4o7D4Xdj8NTgdYLFCrR7Q6v8guKTZ1YnIDRSQXEwByT3EX07lu3VHmL7yIMfOXwbAx2blvlrF6dssgirF9GcjOSRuNyz8N+z+3Ti2+ULDp6HZMMhfyNzaROQqBSQXU0ByL2l2B3N2xDJleQybDp+/Ot60fGH6NougVcUwrFatU5IccGStsTXAoRXGsW8wNBsKDfuDj27ULGI2BSQXU0ByXxsPn2PK8hjmbI/F7jD+apcN9efJphGcOH8ZP28bg9tWuOl9Exbsxe5w8nz7ijldsuQ2TifsnQcL3oST242xwGLQ8iWo3RNsuqWOiFkUkFxMAcn9HT13ic9WHmTG2iMkJqcBxjYByWkOnm5eln91rXL1tRMW7GXsvGiGta94y/AkkiUOB2z7Hha9BecPG2OFy0ObVyGym7YGEDGBApKLKSB5jgvJaXy//gjTVhzk8NlLV8crFw1k7D9rMX/XSYUjca20ZFg/DZa+D5dOG2PFaxtbA5RtZWZlInmOApKLKSB5HrvDyfxdJ5myLIa1B89e91yfpmV47d6qJlUmeUZyIqz8CFZ9BCkXjLGyrY2gVLyWmZWJ5BkKSC6mgOTZth2Np1vUcq4sUcLLauGR+uEMbluBIkF+5hYnud+FU0Y3af1UcKQaY1UfNG5fUricubWJ5HKZ+f6tGwllQlRUFJGRkdSvrx1zPdmiPXE4nEYwAkhzOPlqzWFavr+I/8zZTfylVJMrlFwtIBS6jIaB66DGI4AFdvwEUQ3g92GQeNLsCkUEdZCyRB0kz3Xjguz04+LBfhyPTwIgyM+LZ1uVp3eTMuTz0e7c4mKx22D+m7BvnnHsnR8aPQdNB4NfsLm1ieQymmJzMQUkz3S7q9XSx7vVKs7uE4nsOZkIQFigL4PbVuCR+uF429RsFRc7uBzmvQ7H1hvH+QpB8xeg/lPgralfkeyggORiCkieady8aGxWyx33QRrctgK/bjnGmD+jOXrO2J27dOH8vNChEvdUL6YNJ8W1nE5jN+4Fo+B0tDEWHA6tXoaaj4JVHU2Ru6GA5GIKSLlfSpqDb9YeZuLCvZy+kAJAZLEgXuxUiVYVQ7FoDxtxJXsabPkGFr0DiceNsdAqxs1wK3XWHkoiWaSA5GIKSHnHxeQ0pi6P4ZOlB65uONkgohAvdapE3dK6x5a4WOplWPspLBsDSeeNsfBGxtYABxYbHaWWI25+35LR4LBD65dzsFgR96eA5GIKSHnP2YspfLx4H5+tOkRKmgOAdlWK8GLHSlQqGmhydZLrXT4PKz6E1R9DmjH1S+HycGYftB55fUhaMhoWvX3zuIi4NiD9+uuvmS6offv25MuXL9Pvc1cKSHnX8fOXmbBgL9+tP4LDacx0PFC7BM+3q0h4Id2MVFws4QQs+Q9s/Byc9r/GGw+Aju8oHIn8DZcGJKs1c1fzWCwW9u7dS9myZTP1PnemgCT74i4wdt4eZm+LBcDbZqFHw9IMaF2e0EBfk6uTXO/0Plj4b9g5868xixWcDoUjkTtw+UaRsbGxOByODD3y59dP1ZL7lA8LYFKPuvw6sCnNK4SQancyfeVBWr6/iLF/7iEhSZtNiguFlId/fgb9FkJES2PMaUz9UqC0cTWciNyVTAekXr16ZWq6rGfPnuqySK5Vo2QBvujbkK+eakjNksFcSrEzYeE+Wo5exORlB0hKtf/9SUSyqkRdKNPsysGVK9t+fhqm3wNxu00rSyQ30CLtLNAUm9yK0+lk7o5Y3p+7h/2nLgJQLNiPoe0q8I86JfHSZpOS3a5dc9R0KHz9sHF1G4DVCxoPNKbbfPzNrFLEbbhsDdLly5c5e/YsJUqUuG58x44dVK2ad+6GroAkd5Jmd/DTxmOMmx/NiSu3Lykb6s+LHSrRqVpR7aEk2eN2C7Ln/AtWR/11HFQSOv8HKnfV/kmS57kkIP3www88//zzFCpUCKfTyaeffkrDhg0BqFOnDhs3brz7yj2EApJkRFKqnS9XHyJq0T7OXbkBbs2SwYzoVJmm5UNMrk483qJ377wP0qndcHQdnD9sjFXoaASlQhE5W6eIG3FJQKpVqxbz5s0jNDSU9evX06tXL0aOHEn37t2pXbs2mzZtypbiPYECkmRGYlIqny6LYfKyA1xKMdYkNSsfwosdK1EzvIC5xUnulnLJ2GRyxYfgSAUvP+P+bk2HgJeutpS8xyUBqWrVquzYsePq8ZkzZ3jwwQdp27YtM2fOVAdJ5G+cvpDMRwv38dWaQ6TajX92nasV5YUOlSgfFmBydZKrnd4Ls16AmCXGcaFy0PUDKNfG3LpEcphLLvMPCwtj69atV48LFy7MvHnz2LVr13XjInJrIQG+vHFfVRa+0Ip/1CmJxQJ/bI+lw7glvPTDVo6fv2x2iZJbhVSAJ36Bh6ZCQFE4ux++eAC+7w0Jx82uTsQtZbiDdPToUby8vChatOhNz61YsYKmTZtme3HuSh0kyQ57YhP54M89zNt5EgAfLyu9GpfmuVblKejvY3J1kmslJcDid2HNf429k3wCoPW/oMEzYPMyuzoRl9K92FxMAUmy04ZD5/jPnN2sjTkLQKCvF/1alKVvswj8ffUNS1zkxFZj2u3oWuO4SDXoOgZKNTK3LhEXypGAFBsbe8tuUl6ggCTZzel0siT6FO/P3cOO4wkAhAT4MLB1eR5rWApfL5vJFUqu5HDA5i9h3mtw+ZwxVqsntH8T/HWlpeQ+ORKQatSokWfXHikgias4HE5mbTvBmD/3cPDMJQBKFszHsPYV6VarBDar9rERF7h4Bha8YdwEF8CvALR7A+r0gkzef1PEneVIQKpevTrbtm3LUoGeTgFJXC3V7uC79Uf4cP5e4hKTAahUJJDhHSvRrkqYNpsU1ziyFn4fBievfG0vURe6joXitUwtSyS7qIPkYgpIklMup9iZvvIgHy/eR0JSGgB1ShXgpU6VaVi2sMnVSa5kT4N1k2HhW5CSCBYr1O8HbUaCX7DZ1YncFQUkF1NAkpwWfymV/y3dz9QVMSSlGndtb1UplBc7VqJqcX3TEhdIjIW5I2H7D8axfxh0fBuqP6xblojHUkByMQUkMUtcQhITFu5lxtojpDmMf7r31SzOsPYVKROiG5KKCxxYDLOGw5m9xnGZ5sbVbqGVTC1LJCtyJCDVq1eP9evXZ6lAT6eAJGY7ePoiY+dF8+sWY5M/L6uFR+qHM7htBYoE+ZlcneQ6acmwciIs/QDSLoPVG5oMhBYvgo+CuXgO7YPkYgpI4i52HI/ng7l7WLTnFAB+3laebBpB/xblCM7vbXJ1kuucOwR/vATRfxjHweHGDXArddG0m3iEHAtI7777LkWKFKFPnz7XjU+dOpVTp07x0ksvZfXUbk0BSdzNmgNnGD13DxsOGXvZBPl58Wyr8vRuUoZ8PtpDSbLZ7tlGUIo/bBxX7GQEpYJlTC1L5O/kWEAqU6YMX3/9NU2aNLlufM2aNTz66KPExMRk9dRuTQFJ3JHT6WTBrjjen7uHPScTAQgL9GVIuwr8s1443jbtZyPZKOUSLPsAVkwARyp4+UHz4dB0MHj5ml2dyC3lWEDy8/Nj165dREREXDd+4MABIiMjSUpKyuqp3ZoCkrgzu8PJL5uPMXZeNEfPGTfALVM4P8M6VOKe6sWwarNJyU6nomH2CxCz1DguXB66vA/l2phbl8gtZOb79139SBkeHs6KFStuGl+xYgXFixe/m1OLSBbZrBYerFOSBS+05M37qhIS4MPBM5cY/M0m7pm4nMV74tDSQ8k2oRXhiV/hH1MgoAic2QdfPADfPwkJJ8yuTiTL7iogPfXUUwwdOpRp06Zx6NAhDh06xNSpU3n++efp169fdtUoIlng62WjV5MyLHmxNS+0r0igrxc7TyTQe9o6Hv1k9dX1SiJ3zWKB6g/BwHXQsL+xueSOn+Cj+rBqkrH5pIiHuaspNqfTyf/93/8xYcIEUlJSAGPa7aWXXuK1117LtiLdjabYxBOdvZjCx4v38dmqQ6SkGZtNtqtShBc7VqJS0UCTq5Nc5cQWmPUCHF1nHBepZtyypFRDc+uSPC/HL/O/cOECu3btIl++fFSoUAFf39y9QE8BSTzZ8fOXmbBgL9+tP4LDafzw/0DtEjzfriLhhfKbXZ7kFg4HbPoC5r8Ol690K2v3hHajwF+3yRFzmLIPUvpp8sJNNBWQJDfYF3eBsfP2MHtbLADeNgs9GpZmYJvyhATk7h9yJAddPGOEpE1fGMf5CkK7N6D2E2DVlZWSs3JskTbAlClTqFatGn5+fvj5+VGtWjUmT558t6cVERcrHxbApB51+WVAU5qVDyHV7mT6yoO0GL2IsfOiSUxKNbtEyQ38C0O3j6DPn8ZU2+Vz8NsQmNLemIoTcVN31UF69dVXGTduHIMGDaJx48YArFq1io8++oghQ4bw1ltvZVuh7kQdJMmNlu89zei5u9l6NB6Agvm9GdC6PD0blcbPW5tNSjawp8HaT2DR25BywVjM3eBpaP0v8NNNl8X1cmyKLSQkhIkTJ/LYY49dN/7NN98waNAgTp8+ndVTuzUFJMmtnE4nc3fEMnruHg6cughA8WA/hraryIN1SuClzSYlOyScgD9HwvYfjeOAItDhbeNKuDywTEPMk2MBqWDBgqxdu5YKFSpcNx4dHU2DBg04f/58Vk/t1hSQJLdLszv4aeMxxs2P5kS8seFruVB/XuxYiY5Vi+aJtYaSA/YvgtnDjb2TACJaQJcxxt5KIi6QYwFp0KBBeHt7M3bs2OvGhw8fzuXLl4mKisrqqd2aApLkFUmpdr5cfYioRfs4d8lYk1SzZDAjOlWmafkQk6uTXCEtGVZOgKUfQFoSWL2hySBo8SL46KpKyV45GpA+//xzwsPDadSoEQCrV6/myJEjPPHEE3h7/3U38RtDlCdTQJK8JiEplclLDzB5eQyXUuwANCsfwohOlahRsoC5xUnucO4gzB4Be+cax8GljBvgVu5ialmSu+RYQGrdunWGXmexWFi4cGFWP8btKCBJXnX6QjIfLdzHV2sOkWo3vnR0qV6UFzpUolxogMnVicdzOmHPbPjjJYg/YoxV7Ayd34OCZUwtTXIHU/ZByksUkCSvO3L2EuPmR/PzpmM4ncb93x6uW5Ih7SpQLDif2eWJp0u5CEvfh5UfgSMVvPJBi+HG1JuX9uiSrFNAcjEFJBHDnthE3p+7h/m7TgLg42WlWvEgGkYU5qXOlW96/YQFe7E7nDzfXotwJQNO7TFuWXJwmXFcuDx0+QDKZWz2QuRGLg9Iffr0ydDrpk6dmtlTewQFJJHrbTh0jv/M2c3amLNXx5qUK8ynT9TD39cLMMLR2HnRDGtfkcFtK9zuVCLXczph2w8w919wMc4Yq/YPY1uAoGLm1iYex+UByWq1Urp0aWrXrs2d3v7zzz9n9tQeQQFJ5GZOp5Ml0acYPWcPO08kAJDfx8b/da7M2QspjF+wV+FIsi4pHha+Des+BacDfAKhzUio3w9sXmZXJx7C5QHpueeeY8aMGZQqVYo+ffrQs2dPChUqlOWCPY0CksjtORxOft92gldnbif+8l+3K3m0fjjv/aOGiZVJrnB8M8waBsc2GMdFqsM9YyG8galliWfIkTVIycnJ/PTTT0ydOpWVK1fStWtX+vbtS4cOHXL9JnIKSCJ/L9XuoPKrc7A7/voS8896JXmpU2UK62a4cjccDtj4Gcx/A5LOG2N1noB2b0L+vPPDumReji/SPnToENOnT+fzzz8nNTWVnTt3EhCQey/5VUAS+Xvpa468bZarWwIABOfzZkSnSjxavxQ2a+7+YUpc7OJpmP86bPrSOM5XCNq9AbUfB6tuiyM3y8z372z5G2SxWLBYLDidThwOR3acUkQ82LULsve+3YVhV65aCwnwJf5yKiN/3s6Dk1aw9eh5cwsVz+YfAt2ioM9cCKsKl8/Cb4Nhagc4sdXs6sTDZTkgJScn880339C+fXsqVarEtm3b+Oijjzh8+HCu7h6JyJ3d6mq1wW0rMKx9RU5fSKZVpVACfb3YcjSeblEreGXmNuIvpf7NWUXuoFQjeGYpdHwHfALg6Dr4pCX88X+QlGB2deKh7nqR9pNPPknPnj0pXLiwK+pzS5piE7m9cfOisVktt7xaLX0fpB4NS/HO7F3M3HwcgEL+Pvxf58o8VKckVk27yd1IOG5sCbDjylXUAUWh49vG1gC5fH2s/L0cucy/VKlS1K5d+44Lsn/66afMntojKCCJZI9V+8/w2i/b2Rt3AYB6pQsyqls1Iovr35Xcpf0LYdZwOLvfOI5oAV3GQKg2Kc3LXB6QevfunaEr1aZNm5bZU2e7Bx54gMWLF9O2bVt++OEHABITE2nTpg2pqanY7XYGDx5Mv379MnxOBSSR7JNqdzB1eQwfLtjLpRQ7NquFJxqXZlj7igT6ef/9CURuJy0ZVkyAZR9AWhJYvaHpYGg+HHzym12dmEC3GrnGokWLuHDhAp999tnVgGS320lOTiZ//vxcunSJatWqsW7dugxPEyogiWS/E/GX+ffvO5m9LRaA0EBfXulahftqFs/1W4eIi52NgT9GwN4/jePgUtBlNFTqbG5dkuNy/Co2d9a6dWsCAwOvG7PZbOTPb/z0kJSUhN1uv+OO4CLiesWC8zGpR10+79OAiBB/TiUmM2TGZrp/uoZ9cYlmlyeerFAEdP8OHvkKgkpC/GH45lH45jE4d8js6sRNuXVAWrp0Kffeey/Fixs/Qc6cOfOm10yaNImIiAj8/PyoW7cuy5Yty9C5z58/T82aNSlZsiQjRowgJCQkm6sXkaxoUTGUOUObM7xDRXy9rKw6cIZO45fx7h+7uJicZnZ54qksFqhyDwxcC02HgtUL9syGqIawbAykpZhdobgZtw5IFy9epGbNmnz00Ue3fP7bb79l6NChjBw5kk2bNtG8eXM6d+7M4cOH//bcBQoUYMuWLcTExPD1119z8uTJ7C5fRLLI18vGwDYVmD+sJe2qFCHN4eR/Sw7QfuwS/th2Qh1fyToff2j/JvRfAWWaQ9plWDAK/tsUDiwxuzpxI24dkDp37sxbb73Fgw8+eMvnx44dS9++fXnqqaeoUqUK48ePJzw8nI8//jjDn1GkSBFq1KjB0qVLb/ua5ORkEhISrnuIiOuFF8rP5F71mNKrHiUL5uN4fBLPfrWRXtPWEXP6otnliScLqwy9foMHPgH/UDgdDZ/fBz/0hcRYs6sTN+DWAelOUlJS2LBhAx06dLhuvEOHDqxcufKO7z158uTVkJOQkMDSpUupVKnSbV//7rvvEhwcfPURHh5+978BEcmwtlWKMH9YSwa3KY+PzcrS6FN0HLeUsX/uISnVbnZ54qksFqj5CAxcDw2eBosVtv8AH9WH1f8Fu6Z08zKPDUinT5/GbrdTpEiR68aLFClCbOxf6b9jx448/PDDzJ49m5IlS7Ju3TqOHj1KixYtqFmzJs2aNWPgwIHUqHH7u4y//PLLxMfHX30cOXLEZb8vEbk1P28bwzpUYu7zLWhRMZQUu4MJC/fRftwSFuzSFLnchXwFoMv70G8hFK8DyQkw5yX4tBX8MhCWjL71+5aMhkXv5mSlkoO8XHXi/fuNzbnKlSvnqo8AuOnyX6fTed3Y3Llzb/m+zZs3Z/gzfH198fXV3cdF3EFEiD+fPVmfOdtjGfX7To6cvUzfz9bTrkoRXr83kvBC2t9Gsqh4bXhqPmz8DOa/CbHbjAdA6iXjRrjployGRW9D65GmlCqul+0dpLNnz3LPPfcwceJEJkyYQNeuXTl37lx2fwwhISHYbLbrukUAcXFxN3WVRCR3sVgsdK5ejPnDWvJMy7J4WS3M33WS9uOW8NHCvSSnadpNsshqg3p9jGm3Wj3+Gl8+Dr7pDg7H9eGo5QjzahWXyvaA9PLLLzNs2DDGjx/Phx9+yLBhwxg+fHh2fww+Pj7UrVuXefPmXTc+b948mjRpku2fJyLux9/Xi5c7V+GPIc1pVLYQSakOPvgzms7jl7Fs7ymzyxNPFhAK90+CJ/+AsEhjbM8sGFVI4SiPyPaAtGPHDtq0acOoUaM4e/Ysbdu2ZdeuXVk614ULF9i8efPV6bCYmBg2b9589TL+YcOGMXnyZKZOncquXbt4/vnnOXz4MP3798+u346IeIAKRQL5pl8jPny0FqGBvhw4fZHHp6xlwFcbORF/2ezyxJOVbgLPLIUOb10ZuLLFhNWmRdy5nMvWIF27T4nD4cjSOdavX0/r1q2vHg8bNgyAXr16MX36dB555BHOnDnDqFGjOHHiBNWqVWP27NmULl367ooXEY9jsVjoVqsErSuHMW5eNJ+tPMisbSdYtCeOIW0r0KdZBN42j70uRcxk84bUK0HbYgGn09g7adfvcP/HxpYBkutk+73YRowYQa1atejevTsA33zzDevXr2fMmDHZ+TGm0r3YRNzfzuMJvPrLdjYcMtZAViwSwKhu1WhUNmP3XBS56to1Ry1ehG8fh92/Gc/ZfKD1v6DxILC5rOcg2cTUm9VevnyZgQMHcv78eQCCg4OZNGkSfn5+2fkxplJAEvEMDoeTHzYe5b0/dnP2onEriQdql+DlLpUJC8w9X5PEhW63IPvPV2DlxL+OS9Q1ukmht99TT8xnakBKl5SUhNPpJF++fK44vSmioqKIiorCbrcTHR2tgCTiIc5fSuH9uXv4eu1hnE4I9PXihQ4V6dmoNF6adpM7WfSusd7oVguyF/8HTmyGgysgOR5svkY3qckg4z3idtwiIOVm6iCJeKYtR87z6i/b2Xo0HoDIYkH8+/5q1C1d0OTKxKPFH4PfhsC+K1dVl6h3pZtU0dy65CYKSC6mgCTiuewOJzPWHWb0nD3EX04F4J/1SvJ/natQyN/H5OrEYzmdsOlLmPsvYydumy+0GQmNB6qb5EYUkFxMAUnE8525kMx/5uzmu/VHAQjO582ITpV4rH4prFbL37xb5Dbij17pJs03jkvWh26T1E1yEy4PSH369MnQ66ZOnZrZU3sEBSSR3GPDobOM/Hk7u2MTAahZMpi37q9O9ZLBJlcmHuuW3aRXoPEAdZNM5vKAZLVaKV26NLVr1+ZOb//5558ze2qPoIAkkruk2R18vuoQY+dFcyE5DYsFejQsxYsdKhOc39vs8sRTxR+FXwfD/gXGcckGxu7cIRXMrSsPc3lAeu6555gxYwalSpWiT58+9OzZk0KFCmW5YE+jgCSSO8UlJPH27F38svk4AIX9ffi/zpX5R52SmnaTrHE6YePnMHckpCSCl5/RTWr0nLpJJsiRNUjJycn89NNPTJ06lZUrV9K1a1f69u1Lhw4dsFhy9xcSBSSR3G3l/tO89ssO9sVdAKBe6YL8+/5qVCmmf++SReePwK+D4MAi4zi8obE2KaS8uXXlMTm+SPvQoUNMnz6dzz//nNTUVHbu3ElAQMDdntZtKSCJ5H4paQ6mrYjhwwV7uZRix2a10KtxGZ5vX4FAP027SRY4nbDxM5j7yjXdpFeh0bPqJuWQzHz/zpYd0iwWCxaLBafTmeX7romIuBMfLyvPtCzH/GEt6VK9KHaHk6krYmg7Zgm/bD52x/WXIrdksUDd3vDcKijbGtKS4M+RMK0LnN5ndnVygywHpOTkZL755hvat29PpUqV2LZtGx999BGHDx/Otd2jqKgoIiMjqV+/vtmliEgOKV4gH5N61OWzPg2ICPEnLjGZITM202PyGvbFJZpdnniiAuHw+M9wz3jwCYAjq+G/TWFVFDjsZlcnV9z1Iu0nn3ySnj17Urhw3rkBpKbYRPKm5DQ7nyw5wEeL9pGc5sDbZqFvs7IMblue/D66UalkwfnDV9YmLTaOSzWGblFQuJypZeVWOXKZf6lSpahdu/YdF2T/9NNPmT21R1BAEsnbjpy9xJu/7WD+rjgAigf78dq9kXSsWjTXX6QiLuB0woZp8OerkHIBvPJBu9ehwTNg1b0Cs5PLA1Lv3r0z9EVg2rRpmT21R1BAEhGA+TtP8sZvOzh67jIALSuG8uZ9VSkT4m9yZeKRzh0yukkxS4zjUk2g20fqJmUj3WrExRSQRCTd5RQ7kxbv439LDpBid+DjZaV/y3I816ocft66MkkyyemE9VNh3mvqJrmAApKLKSCJyI1iTl/ktV+2s2zvaQDCC+Xjzfuq0qZyEZMrE4907hD8OhBilhrHpZsa3aRCZc2ty8O59DL/rVu3ZupS/h07dpCWlpbZjxER8SgRIf583qcBk3rUoWiQH0fOXqbP9PX0+3w9R85eMrs88TQFS8Pjv0DXMeDtD4dWwMdNYc3/QNvp5IhMd5BsNhuxsbGEhoZm6PVBQUFs3ryZsmVzT+pVB0lE7uRichoTFuxlyvIY0hxO/LytDGxdnn4tyuLrpWk3yaRzB+GXgXBwmXFcutmVblKEqWV5IpdOsVmtVp5++mny58+foddPmjSJnTt3KiCJSJ4TfTKRV2duZ03MWQDKhvjzZreqNK+QsR8wRa5yOGD9FJj3OqReBO/80O5NqP+U1iZlgksDUqtWrTJ9GevXX39NsWLFMvUed6aAJCIZ5XQ6+WXzcd6atYvTF5IB6FqjGK92jaRosJ/J1YnHORtjXOmW3k0q0xzum6huUgZpkbaLKSCJSGYlJKUy9s9oPl91EIcT/H1sDG1Xkd5Ny+BtUwdAMuFqN+k1SL1krFFq/ybU66tu0t9QQHKRqKgooqKisNvtREdHKyCJSKbtOB7PqzO3s/HweQAqFgng392q0bBs3rkbgWSTszHG2qRDy43jMs2NtUkFy5haljtTQHIxdZBE5G44HE5+2HiU9/7YzdmLKQA8ULsEL3epTFigpt0kExwOWDcZ5r+ublIGKCC5mAKSiGSH85dSeH/uHr5eexinEwJ9vXihQ0V6NiqNl6bdJDPOHrjSTVphHJdpbtzTrWBpc+tyMwpILqaAJCLZacuR87z6y3a2Ho0HILJYEG89UI06pQqaXJl4FIcD1n4C89+AtMtGN6nDKKjbR92kKxSQXEwBSUSym93h5Ju1h3l/7h7iL6cC8Ei9cF7qXJlC/j4mVyce5cx+o5t0eKVxHNHSWJtUoJS5dbkBl+6kDdCmTRvOnz9/2+dPnz6dq/Y9EhFxNZvVQs9GpVn4QkserlsSgG/XH6HNmMV8veYwDod+lpUMKlwOes+CTv8x7uUWswQmNTbu8aaeSIZlqYNktVqJjY0lLCwMgOTkZHx9fa8+f/LkSYoXL47dbs++St2IOkgi4mrrD57llZnb2R2bCEDN8AK81a0a1UsGm1yZeJQz++GXAXB4lXGcx7tJLu8gXWv79u2UL1+eV155Bc3WiYhkj3plCvH7oGa8dk8kAb5ebDlynvuilvPqzO3EX0o1uzzxFIXLQe/Z0PFddZMy6a46SHv27KFbt27cc889zJ07l6pVqzJjxgycTqc6SCIi2SQuIYm3Z+/il83HASjs70Od0gWpVjyIIe0q3vT6CQv2Ync4eb79zc9JHnZmP8x8Do6sNo7LtjZ24S4Qbm5dOShHOkg///wznTp14pVXXuHzzz9n48aNpKSkUKtWLRYvXpzV04qIyA3Cgvz48NHafN2vIeXDAjhzMYV5O08ybv5eXv1l+3WvnbBgL2PnRWOzZu6WUJIHFC4HT86Gju+Alx8cWGR0kzZMVzfpFrLcQfLx8WHKlCn06NHj6rjdbmfEiBGMHz/+6nFupA6SiJglJc3B1BUxfDh/L5dTja+xdUoV4LM+DZi24iBj50UzrH1FBretYHKl4tZu7CaVawP3Tsj13SSXX+b/5JNP0r17d9q3b3/L52fOnMkvv/zCtGnTMntqj6CAJCJmO37+Mv/+fSd/bI+9blzhSDLMYYfVH8PCf0NaEvgEQse3oc4TkMmb0nsKt9gHafPmzdSqVcsVpzadApKIuIsl0afoNXXt1eMHapfg9XsjKZBfeydJBp3ea3STjl75e1SuLdw3AYJLmluXC+ToVWzXio+PZ9KkSdStW5d69epl56lFROQWthw5D0D6kqOfNx2j/bilzNt50ryixLOEVIA+c6DDW8bapP0LjLVJGz/P02uTsiUgLVy4kJ49e1KsWDHefPNNypQpkysv+Y+KiiIyMpL69eubXYqIyNUF2cPaV+TAu115rIGxfuRUYjL9Pl/P0BmbOHflZrgid2S1QZNB0H85lKwPyQnw6yD46iGIP2Z2dabI8hTb0aNHmT59OtOmTePkyZN069aNHj160LFjR3bt2kXNmjW1SFtExEWuDUfXrjka++ceJizchwVwAiEBvrzzQDU6VC1qWq3iYRx2WBUFC98CezL4BhlXvtXu6fFrk1y+BqlLly4sWrSINm3a0L17d+6//378/f2vPr9jxw5q1KihgCQi4iLjrlzKf6sF2RMW7OXY+ctsOHSOfXEXAOhWqzhv3FuVgrqvm2TUqWj45Tk4us44Lt8e7v0QgkuYW9ddcHlAslqtdO/enaFDh95yrZECkoiI+ZJS7Xy4YC//W7Ifh9PoJr39QDU6qpskGeWww6qPYOHbV7pJwdDpHajVwyO7SS5fpL1ixQry5ctHmzZtqFSpEqNGjWLfvn1ZKlZERFzDz9vGS50q8/NzTakQFsDpC8k888UGBn+zibNamyQZYbVB0yHQfxmUqAfJ8ca93b7+JyQcN7s6l7qry/wvXbrEjBkzmDp1KqtWraJ+/fr06NGDqlWr0r59e3WQRETcRFKqnQkL9vLfq90kH966vxqdqhUzuzTxFPY0o5u06J1ruknvQq3uHtNNMmUfpD179jBlyhS++OILTp48icViUUASEXEzW46c58UfthB90libdG/N4rx5X1UKaW2SZFTcbmNt0rENxnGFjnDveAgqbmpZGWHqRpF2u53ffvuNqVOn8uuvv2bnqd2GApKIeLLktPRu0gHsDqe6SZJ59jRYNfFKNykF/IKh03tQ8zG37ia5xU7auZkCkojkBluPnufF77ey52QiAPfUKMab91WlcICvyZWJx4jbDTOfheMbjeOKneCe8RDknmHbtJ20RUTEc9QoWYBfBzVlYOvy2KwWft96gg7jljJ72wmzSxNPEVYZ+s6Dtq+DzQei58CkhrBlhsfvwq0OUhaogyQiuc22o/EM/37L1W5S1+rFGNVN3STJhLhdV7pJm4xjN+wmqYMkIiKZUr1kML8OasqgNkY3ada2E7Qft5RZW9VNkgwKqwJ950Pb13JFN0kdpCxQB0lEcrPtx4xu0u5Yo5vUpXpRRnWrRoi6SZJRJ3ca3aQTm43jSl3gnnEQaO4mpeogiYhIllUrEcyvA5sxuG0FvKwWZm+LpcO4pfy+NXdvDCjZqEgkPLUA2rwKVm/YMxuiGsLW7zymm6QOUhaogyQieYW6SXLXbuomdb3STSqS46Wog+QiUVFRREZGUr9+fbNLERHJEbfqJrUfu4TfthxHP19LhhSJhKfmQ+tXrnSTZhlrk7Z+79bdJHWQskAdJBHJi7Yfi+fFH7ay60QCAJ2qFuXf91cjNFDdJMmg2O1GNyl2q3Fc+R7oOjbHuknqIImISLarViKYXwY0ZWg7o5s0Z0csHcYt4Vd1kySjilaDfguh9Uijm7T7d6ObtO0Ht+smqYOUBeogiUhet+N4PMO/VzdJ7sKN3aTCFaByF2g/6ubXLhkNDju0fvmuPlIdJBERcamqxYP5deD13aT245bwy+Zj6iZJxqR3k1r9C6xecGYvrPgQvut1fTdpyWhY9DZYbTlanjpIWaAOkojIX3YeT2D491vYeaWb1CGyCG89UI2wQD+TKxOPEbvtSjdpm3EcUgl6z4IN04xw1HoktBxx1x+jm9W6mAKSiMj1Uu0OPl68n4kL95Jqd1Igvzdv3leV+2oWx+LGd3cXN2JPhWVjYMl/wOkALIAz28IRKCC5nAKSiMit7TphdJN2HDe6Se0ji/D2/dUIC1I3STLoxFb4XwvAadyy5NVT2XZqrUESERFTVCkWxMwBTXmhfUW8bRbm7TxJ+3FLmblJa5Mkg6LncDUc2VOMNUgmUEASEZFs5W2zMqhtBX4b1IxqJYKIv5zK0G830+/zDcQlJJldnriz9AXZrUcanaPWI41jE0KSApKIiLhE5aJB/PxcU4Z3MLpJ83cZ3aSfNx1VN0ludm04Sl9z1HKEaSFJAUlERFzG22ZlYJvru0nPf7uFfp+v56S6SXIth/3WC7LTQ5LDnqPlaJF2FmiRtohI5qXaHXyy9ADj50eTancS5OfF6/dW5cE6JXSlm+QILdIWERG3422zMqB1eX4f1JwaJYNJSErjhe+38NRn6iaJ+1FAEhGRHFWpaCA/PduEFztWwsdmZcHuONqPXcKPG7Q2SdyHApKIiOQ4r/Ru0uBm13WT+n62nth4dZPEfApIIiJimopFjG7SiE5GN2nh7jjaj1vC9+uPqJskplJAEhERU3nZrDzXqjyzBjejZslgEpPSePGHrfSZvk7dJDGNApKIiLiFCkUC+fHZJrzUqTI+NiuL9pyi/bglfKdukphAASkToqKiiIyMpH79+maXIiKSK3nZrDzbqpzRTQovQGJSGiN+2MqT09dxIv6y2eVJHqJ9kLJA+yCJiLhemt3B5OUxjJ0XTUqag0BfL169J5KH65XUvkmSJdoHSUREPJ6XzUr/luWYPbgZtcILkJicxogft9J72jqOn1c3SVxLAUlERNxa+TBjbdLLnSvj42VlSfQpOo5byrfrDmttkriMApKIiLg9m9XCMy3LMXtwc2qXMrpJL/24jSemruWYukniAgpIIiLiMcqHBfBD/yaM7FIFXy8ry/aepuO4pcxYq26SZC8FJBER8Sg2q4V+Lcoye0hz6pQqwIXkNP7vJ3WTJHspIImIiEcqFxrA97foJn2jbpJkAwUkERHxWNd2k+qWLsiF5DRevtJNOnruktnliQdTQBIREY9XLjSA755pzCtd/+omdRq/jK/XqJskWaOAJCIiuYLNauGp5mX5Y0hz6l3pJv3r5208PkXdJMk8BSQREclVyoYG8O0zjXn1nkj8vK0s32esTfpqzSF1kyTDFJBERCTXsVkt9G0WwR9DWlCvdEEuptgZ+fN2ek5Zw5Gz6ibJ31NAEhGRXCsixJ9vn2nMa1e6SSv2naHT+KV8ufoQDoe6SXJ7CkgiIpKr2awW+lzpJtUvY3STXpmpbpLcmQKSiIjkCREh/nz7dGNev9foJq3cf4aO45fyxaqD6ibJTRSQREQkz7BaLTzZNII5Q1rQIKIQl1LsvPrLDnpMVjdJrqeAJCIieU6ZEH9m9GvEG/dGks/bxqoD6ibJ9RSQREQkT7JaLfRuGsGcoc2v6yZ1n7yaw2fUTcrrFJBERCRPK13Y6Ca9eV9V8nnbWH3gLB3HL+Wzleom5WUWp3bNyrSEhASCg4OJj48nKCjI7HJERCSbHD5ziRd/2MKamLMAlCjgR8eqRXnt3qo3vXbCgr3YHU6eb18xp8uULMrM9291kERERK4oVTg/3/RrxKhuVcnvY+PY+SSmrjjIk9PWXtdNmrBgL2PnRWOzWkysVlxJAUlEROQaVquFJxqXYc6QFjQqWwiARXtO0Xz0Ig6duXg1HA1rX5HBbSuYXK24iqbYskBTbCIieYPD4eSrNYcY9ftOUu1/fbt8vl0FhrTT1Jqn0RSbi0RFRREZGUn9+vXNLkVERHKA1Wrh8cZlWPhCK66dTNt2LJ4zF5JNq0tcTx2kLFAHSUQkb7l2zZH9ylqk0EBfxjxckxYVQ02uTjJKHSQREZFscu2ao/3vdKFno9IAnEpM5ompaxn1206SUu0mVynZTQFJRETkNm61IPut+6sxuE35q6+ZuiKG+6NWEH0y0awyxQUUkERERG7D7nDe8mq1YR0qMax9Re6rWYzC/j7sjk3k3onL+WzlQbRyJXfQGqQs0BokERFJdyoxmRd/2MLiPacAaF0plNEP1SQ00NfkyuRGWoMkIiKSQ0IDfZnWuz5v3BuJj5eVRXtO0fnDpSzaHWd2aXIXFJBERETuksVi3Pj2t4HNqFw0kNMXUnhy+jpe/2W7FnB7KAUkERGRbFKpaCAzBzTlyaZlAPhs1SHu+2g5u04kmFuYZJoCkoiISDby87bx+r1Vmf5kfUICfIk+eYFuH61gyvKY6+7nJu5NAUlERMQFWlUKY+7Q5rStHEaK3cG/f99J7+nriEtIMrs0yQAFJBERERcpHODL5F71+Pf91fD1srI0+hSdPlzGvJ0nzS5N/oYCkoiIiAtZLBYeb1SaWYObEVksiLMXU+j3+XpG/ryNyylawO2uFJBERERyQPmwQH4e0IR+zSMA+GrNYe6ZuIztx+JNrkxuRQFJREQkh/h62RjZNZIv+zYkLNCX/acu8sCkFXyydL8WcLsZBSQREZEc1qxCCHOHtqBDZBFS7U7emb2bx6euITZeC7jdhQKSiIiICQr6+/C/x+vy7oPVyedtY8W+M3T6cClztseaXZqggCQiImIai8XCYw1K8fvgZlQvEcz5S6n0/3ID//fjVi6lpJldXp6mgCQiImKycqEB/PhsE55tVQ6LBWasO8I9E5az9eh5s0vLsxSQRERE3ICPl5WXOlXmq6caUjTIjwOnL/LgpJVMWrwPuxZw5zgFJBERETfSpFwIc4Y2p0v1oqQ5nIyes4fun67m+PnLZpeWpyggiYiIuJkC+X2I6l6H0Q/VIL+PjTUxZ+k0fimztp4wu7Q8QwFJRETEDVksFv5ZL5zZg5tTM7wACUlpDPh6I8O/38KFZC3gdjUFJBERETdWJsSfH/o3ZmDr8lgs8MOGo3SdsIxNh8+ZXVqupoAkIiLi5rxtVoZ3rMSMfo0oUSAfh85c4qH/rmLigr1awO0iCkgiIiIeomHZwswe0px7axbH7nAyZl40j36yiqPnLpldWq6jgCQiIuJBgvN5M+HRWox7pCYBvl6sO3iOzuOX8cvmY2aXlqsoIImIiHgYi8XCA7VLMntwc+qUKkBichpDZmzm+W83k5CUanZ5uYICkoiIiIcqVTg/3z3TmKHtKmC1wM+bjtHlw2VsOHTW7NI8ngKSiIiIB/OyWRnariLf929MyYL5OHruMg//dxXj5kWTZneYXZ7HUkASERHJBeqWLsQfQ5rzYO0SOJzw4YK9/PN/qzh8Rgu4s0IBSUREJJcI9PNm7CO1+PDRWgT6ebHx8Hm6TFjGTxuP4nRqO4DMUEASERHJZbrVKsEfQ5pTv0xBLiSnMey7LQyesZn4y1rAnVEKSCIiIrlQyYL5mfF0Y4Z3qIjNauG3Lcfp8uEy1sZoAXdGKCCJiIjkUjarhYFtKvBD/8aULpyfY+cv8+gnq/hg7h5StYD7jhSQMiEqKorIyEjq169vdikiIiIZVrtUQWYNbs7DdUvicMJHi/bx0H9XcfD0RbNLc1sWp1ZtZVpCQgLBwcHEx8cTFBRkdjkiIiIZNmvrCV7+aSsJSWnk97Hxxn1VebhuSSwWi9mluVxmvn+rgyQiIpKHdK1RjDlDW9AwohCXUuyM+GErA7/exPlLKWaX5lYUkERERPKY4gXy8XW/RozoVAkvq4VZ207Q+cNlrNp/xuzS3IYCkoiISB5ks1p4rlV5fnquCREh/pyIT6L75NW898duUtK0gFsBSUREJA+rUbIAswY347EG4Tid8N8l+/nHxyvZf+qC2aWZSgFJREQkj8vv48W7D9bgvz3rUiC/N9uOxXPPhOV8s/Zwnt2BWwFJREREAOhUrShzhrSgafnCXE618/JP2+j/5QbOXcx7C7gVkEREROSqosF+fNGnISO7VMHbZmHujpN0+nApy/eeNru0HKWAJCIiItexWi30a1GWn59rSrlQf04mJNNzyhrenrWT5DS72eXlCAUkERERuaVqJYL5fVBzejYqBcCny2J4IGol++ISTa7M9RSQRERE5Lby+dh46/7qTH6iHoX8fdh5IoGuE5bzxepDuXoBtwKSiIiI/K12kUWYM6Q5zSuEkJzm4NWZ2+n3+XrOXEg2uzSXUEASERGRDAkL8uOzJxvw2j2R+NiszN8VR8fxy1i8J87s0rKdApKIiIhkmNVqoU+zCH4Z2JSKRQI4fSGZ3tPW8eZvO0hKzT0LuBWQREREJNOqFAvi14HN6N2kDADTVhzk/qgV7InNHQu4FZBEREQkS/y8bbxxX1Wm9a5PSIAPu2MTufej5UxfEePxC7gVkEREROSutK4cxh9DWtC6UigpaQ7e+G0nT05fx6lEz13ArYAkIiIidy000JepveszqltVfL2sLN5zik7jl7Jw90mzS8sSBSQRERHJFhaLhScal+G3Qc2oXDSQMxdT6DN9Pa/9st3jFnArIImIiEi2qlgkkJkDmtK3WQQAn686xL0Tl7PzeILJlWWcApKIiIhkOz9vG6/eE8nnfRoQGujL3rgL3B+1gsnLDuBwuP8CbgUkERERcZkWFUOZM6Q57aoUIcXu4K1Zu+g1bS1xCUlml3ZHCkgiIiLiUoUDfPn0ibq8/UA1/LytLNt7mo7jl/LnjlizS7stBSQRERFxOYvFQo+Gpfl9UHOqFg/i3KVUnv5iA//6eRuXU9xvAbfF6ek7OZkgISGB4OBg4uPjCQoKMrscERERj5KcZmfsn9H8b+kBAMqG+tMwohDFgvMxuG2Fm14/YcFe7A4nz7eveFefm5nv31539UkiIiIimeTrZePlLlVoUTGUYd9t5sCpixw8fRGHE5xOJ0Pa/RWEJizYy9h50Qy7y3CUWZpiExEREVM0LR/CnCEt6Fi1COkXto2bv5d3Zu0Crg9Ht+osuZKm2LJAU2wiIiLZx+l08t36I7zx604uX9lQ0stqIc3hzNZwlJnv3+ogiYiIiKksFguP1C/FrMHNqFEyGIA0hxMfmzXHO0fpFJBERETELZQNDaB1pTAAvG0WUuwOJizYa0otWqQtIiIibmHCgr18uGDv1Wm19DVIQI53khSQRERExHS3WpCd/qsZIUkBSURERExnv82C7PRjew7fv01XsWWBrmITERHxPLqKTUREROQuKCCJiIiI3EABSUREROQGCkgiIiIiN1BAEhEREbmBApKIiIjIDRSQRERERG6ggCQiIiJyAwUkERERkRsoIImIiIjcQPdiy4L0u7MkJCSYXImIiIhkVPr37YzcZU0BKQsSExMBCA8PN7kSERERyazExESCg4Pv+BrdrDYLHA4Hx48fJzAwEIvFku3nr1+/PuvWrcv285rF3X8/7lBfTtfg6s9zxfmz65wJCQmEh4dz5MgR3Ww6D3GHf+eeILf9f7rx9+N0OklMTKR48eJYrXdeZaQOUhZYrVZKlizpsvPbbLZc9YXb3X8/7lBfTtfg6s9zxfmz+5xBQUGm/7lLznGHf+eeILf9f7rV7+fvOkfptEjbDQ0YMMDsErKVu/9+3KG+nK7B1Z/nivO7w5+TeC79/cmY3Pb/6W5+P5piE5E8JSEhgeDgYOLj43PVT8oikr3UQRKRPMXX15fXX38dX19fs0sRETemDpKIiIjIDdRBEhEREbmBApKIiIjIDRSQRERERG6ggCQiIiJyAwUkERERkRsoIImIXOP333+nUqVKVKhQgcmTJ5tdjoiYRJf5i4hckZaWRmRkJIsWLSIoKIg6deqwZs0aChUqZHZpIpLD1EESEbli7dq1VK1alRIlShAYGEiXLl2YO3eu2WWJiAkUkEQk11i6dCn33nsvxYsXx2KxMHPmzJteM2nSJCIiIvDz86Nu3bosW7bs6nPHjx+nRIkSV49LlizJsWPHcqJ0EXEzCkgikmtcvHiRmjVr8tFHH93y+W+//ZahQ4cycuRINm3aRPPmzencuTOHDx8G4FYrDiwWi0trFhH3pIAkIrlG586deeutt3jwwQdv+fzYsWPp27cvTz31FFWqVGH8+PGEh4fz8ccfA1CiRInrOkZHjx6lWLFiOVK7iLgXBSQRyRNSUlLYsGEDHTp0uG68Q4cOrFy5EoAGDRqwfft2jh07RmJiIrNnz6Zjx45mlCsiJvMyuwARkZxw+vRp7HY7RYoUuW68SJEixMbGAuDl5cWYMWNo3bo1DoeDESNGULhwYTPKFRGTKSCJSJ5y45oip9N53dh9993Hfffdl9NliYib0RSbiOQJISEh2Gy2q92idHFxcTd1lUREFJBEJE/w8fGhbt26zJs377rxefPm0aRJE5OqEhF3pSk2Eck1Lly4wL59+64ex8TEsHnzZgoVKkSpUqUYNmwYjz/+OPXq1aNx48Z88sknHD58mP79+5tYtYi4I91qRERyjcWLF9O6deubxnv16sX06dMBY6PI0aNHc+LECapVq8a4ceNo0aJFDlcqIu5OAUlERETkBlqDJCIiInIDBSQRERGRGyggiYiIiNxAAUlERETkBgpIIiIiIjdQQBIRERG5gQKSiIiIyA0UkERERERuoIAkIrle7969sVgsWCwWZs6caVodb7zxxtU6xo8fb1odIvL3FJBExGNcG3S8vLwoVaoUzz77LOfOnfvb93bq1IkTJ07QuXPnq2Pp51q9evV1r01OTqZw4cJYLBYWL16cbfUPHz6cEydOULJkyWw7p4i4hgKSiHiU9KBz8OBBJk+ezG+//cZzzz33t+/z9fWlaNGi+Pr6XjceHh7OtGnTrhv7+eefCQgIyNa6AQICAihatCg2my3bzy0i2UsBSUQ8SnrQKVmyJB06dOCRRx7hzz//zPL5evXqxYwZM7h8+fLVsalTp9KrV6/rXnfw4EEsFgszZsygSZMm+Pn5UbVq1Zs6TDt27KBr164EBQURGBhI8+bN2b9/f5brExFzKCCJiMc6cOAAc+bMwdvbO8vnqFu3LhEREfz4448AHDlyhKVLl/L444/f8vUvvvgiL7zwAps2baJJkybcd999nDlzBoBjx47RokUL/Pz8WLhwIRs2bKBPnz6kpaVluT4RMYeX2QWIiGTG77//TkBAAHa7naSkJADGjh17V+d88sknmTp1Kj179mTatGl06dKF0NDQW7524MCB/OMf/wDg448/Zs6cOUyZMoURI0YQFRVFcHAwM2bMuBraKlaseFe1iYg51EESEY/SunVrNm/ezJo1axg0aBAdO3Zk0KBBd3XOnj17smrVKg4cOMD06dPp06fPbV/buHHjq//t5eVFvXr12LVrFwCbN2+mefPmd9XREhH3oIAkIh7F39+f8uXLU6NGDSZMmEBycjJvvvnmXZ2zcOHC3HPPPfTt25ekpKTrrnTLCIvFAkC+fPnuqg4RcR8KSCLi0V5//XU++OADjh8/flfn6dOnD4sXL+aJJ56441Vm124JkJaWxoYNG6hcuTIANWrUYNmyZaSmpt5VLSJiPgUkEfForVq1omrVqrzzzjt3dZ5OnTpx6tQpRo0adcfXRUVF8fPPP7N7924GDBjAuXPnrk7JDRw4kISEBB599FHWr1/P3r17+eKLL9izZ89d1SYiOU8BSUQ83rBhw/j00085cuRIls9hsVgICQnBx8fnjq977733+M9//kPNmjVZtmwZv/zyCyEhIYAxVbdw4UIuXLhAy5YtqVu3Lp9++qnWJIl4IF3FJiIeY/r06bcc7969O927d8/0+ZxO522fK1CgwC2fr1Klyk07b1+rRo0azJ07N9O1iIh7UQdJRPKE9O0Bfv/9d9NqeOeddwgICODw4cOm1SAiGWNx3ulHKBGRXCAuLo6EhAQAihUrhr+/f6bPcfDgQSIiIti0aRO1atXKUh1nz57l7NmzAISGhhIcHJyl84iI6ykgiYiIiNxAU2wiIiIiN1BAEhEREbmBApKIiIjIDRSQRERERG6ggCQiIiJyAwUkERERkRsoIImIiIjcQAFJRERE5AYKSCIiIiI3+H9SlyTWp10eAgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "plt.loglog(r_arr, DS0, label=\"bin_rich_0\", marker=\"x\")\n", "plt.loglog(r_arr, DS2, label=\"bin_rich_2\", marker=\"x\")\n", @@ -610,13 +805,21 @@ "plt.ylabel(\"$\\Delta\\Sigma$ [M$_\\odot \\;$Mpc$^{-2}$]\")\n", "plt.legend()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "778073c8-cdaa-40f1-8a0e-1f2da827347f", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "wrk", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "wrk" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -627,7 +830,8 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3" + "pygments_lexer": "ipython3", + "version": "3.11.4" } }, "nbformat": 4, diff --git a/examples/cluster_number_counts/generate_SDSS_ClusterCountsMass_sacc_data.py b/examples/cluster_number_counts/generate_SDSS_ClusterCountsMass_sacc_data.py new file mode 100755 index 00000000..dcac00f4 --- /dev/null +++ b/examples/cluster_number_counts/generate_SDSS_ClusterCountsMass_sacc_data.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python + +"""Defines a function to generate a SACC file for cluster number counts.""" + +# Cluster SACC file creation for SDSS Count+Mean(logM) data vector +# in Costanzi et al. 2019, https://arxiv.org/pdf/1810.09456 + +import itertools +import numpy as np +from numcosmo_py import Ncm +from typing import Any +import sacc + + +def generate_SDSSCL_sacc_file() -> Any: + """ + Generate a SACC file with SDSS cluster number counts and Mean(logM). + + According to SDSS Costanzai et al. 2019, arxiv 1810.09456. + Output: will save a fits file to disk. + """ + # settting cluster data vector + # according to this Costanzi et al. 2019 arxiv 1810.09456 + area = 10263.037032448827 + z_edges = np.array([0.1, 0.3]) + richness_edges = np.log10(np.array([20, 27.9, 37.6, 50.3, 69.3, 140])) + cluster_counts = np.array([3711, 1788, 978, 476, 223]) + mean_logM = np.array([14.111, 14.263, 14.380, 14.609, 14.928]) + std_logM1 = np.array([0.024, 0.030, 0.033, 0.036, 0.029]) + std_logM2 = np.array([0.026, 0.024, 0.026, 0.028, 0.036]) + std_counts = np.array([100, 61, 41, 27, 18]) + var_mean_logM = std_logM1**2 + std_logM2**2 + var_counts = std_counts**2 + + # ** Correlation matrix - the "large blocks" correspond to the $N_z$ redshift bins. + # In each redshift bin are the $N_{\rm richness}$ richness bins.** + + covariance = np.diag( + np.concatenate((var_counts.flatten(), var_mean_logM.flatten())) + ) + + s_count = sacc.Sacc() + bin_z_labels = [] + bin_richness_labels = [] + + survey_name = "SDSSCluster_redshift_richness" + s_count.add_tracer("survey", survey_name, area) + + for i, z_bin in enumerate(zip(z_edges[:-1], z_edges[1:])): + lower, upper = z_bin + bin_z_label = f"bin_z_{i}" + s_count.add_tracer("bin_z", bin_z_label, lower, upper) + bin_z_labels.append(bin_z_label) + + for i, richness_bin in enumerate(zip(richness_edges[:-1], richness_edges[1:])): + lower, upper = richness_bin + bin_richness_label = f"rich_{i}" + s_count.add_tracer("bin_richness", bin_richness_label, lower, upper) + bin_richness_labels.append(bin_richness_label) + + # pylint: disable-next=no-member + cluster_count = sacc.standard_types.cluster_counts + # pylint: disable-next=no-member + cluster_mean_log_mass = sacc.standard_types.cluster_mean_log_mass + + counts_and_edges = zip( + cluster_counts.flatten(), itertools.product(bin_z_labels, bin_richness_labels) + ) + + mean_logM_and_edges = zip( + mean_logM.flatten(), itertools.product(bin_z_labels, bin_richness_labels) + ) + + for counts, (bin_z_label, bin_richness_label) in counts_and_edges: + s_count.add_data_point( + cluster_count, (survey_name, bin_z_label, bin_richness_label), int(counts) + ) + + for bin_mean_logM, (bin_z_label, bin_richness_label) in mean_logM_and_edges: + s_count.add_data_point( + cluster_mean_log_mass, + (survey_name, bin_z_label, bin_richness_label), + bin_mean_logM, + ) + + # ### Then the add the covariance and save the file + + s_count.add_covariance(covariance) + s_count.to_canonical_order() + s_count.save_fits("SDSSTEST_redshift_richness_sacc_data.fits", overwrite=True) + + +if __name__ == "__main__": + Ncm.cfg_init() + generate_SDSSCL_sacc_file()