From 223450839b5a4d50bcbc6ae3c247e63d4d4bd2bc Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Wed, 5 Jun 2024 13:09:10 -0400 Subject: [PATCH] checkpoint --- incubator/mask/cone_extrema.ipynb | 239 +++++++++++++++++++++++++++ incubator/mask/min_max_hist.ipynb | 263 ++++++++++++++++++++++++++++++ 2 files changed, 502 insertions(+) create mode 100644 incubator/mask/cone_extrema.ipynb create mode 100644 incubator/mask/min_max_hist.ipynb diff --git a/incubator/mask/cone_extrema.ipynb b/incubator/mask/cone_extrema.ipynb new file mode 100644 index 0000000..60cbd10 --- /dev/null +++ b/incubator/mask/cone_extrema.ipynb @@ -0,0 +1,239 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [], + "source": [ + "import hipscat\n", + "import healpy as hp\n", + "import pandas as pd\n", + "from tqdm import tqdm\n", + "from hipscat.inspection import plot_pixel_list\n", + "from hipscat.pixel_math import HealpixPixel\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_pixel_list([HealpixPixel(0,11), HealpixPixel(4, 78)])" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [], + "source": [ + "step = 2" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [], + "source": [ + "## i know this is dumb. i don't care.\n", + "\n", + "def min_max_sep(bounds_a, bounds_b):\n", + " min_sep = float(\"inf\")\n", + " max_sep = 0.0\n", + "\n", + " for i in range(4*step):\n", + " for j in range (4*step):\n", + " sep_sq = (bounds_a[0][i]-bounds_b[0][j])**2 + (bounds_a[1][i]-bounds_b[1][j])**2\n", + " min_sep = min(min_sep, sep_sq)\n", + " max_sep = max(max_sep, sep_sq)\n", + " return (min_sep, max_sep)" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3933 partitions\n", + "7732278 iterations\n" + ] + } + ], + "source": [ + "\n", + "\n", + "gaia_full_partition_frame = pd.read_csv(\"gaia_partition_info.csv\")\n", + "gaia_full_partition_list = [\n", + " HealpixPixel(order, pixel)\n", + " for order, pixel in zip(\n", + " gaia_full_partition_frame[\"Norder\"],\n", + " gaia_full_partition_frame[\"Npix\"],\n", + " )\n", + " ]\n", + "num_partitions = len(gaia_full_partition_list)\n", + "print(num_partitions, \"partitions\")\n", + "print(int(.5 * num_partitions * (num_partitions - 1)), \"iterations\")" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3933/3933 [09:51<00:00, 6.64it/s] \n" + ] + } + ], + "source": [ + "\n", + "all_seps = []\n", + "\n", + "for a in tqdm(range(0, num_partitions)):\n", + " for b in range(a, num_partitions):\n", + "\n", + " bounds_a = hp.vec2dir(hp.boundaries(2**gaia_full_partition_list[a].order, gaia_full_partition_list[a].pixel, step=step, nest=True), lonlat=True)\n", + " bounds_b = hp.vec2dir(hp.boundaries(2**gaia_full_partition_list[b].order, gaia_full_partition_list[b].pixel, step=step, nest=True), lonlat=True)\n", + " all_seps.append(min_max_sep(bounds_a, bounds_b))" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "7736211" + ] + }, + "execution_count": 42, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(all_seps)" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(0.0, 506.24999999999966),\n", + " (0.0, 1236.6206396835553),\n", + " (6299.120639683556, 15411.620639683555),\n", + " (4972.652051530026, 13072.652051530025),\n", + " (8516.402051530025, 18641.402051530025),\n", + " (11289.449049432578, 23391.34255485142),\n", + " (7239.44904943258, 16351.949049432578),\n", + " (10791.34255485142, 24620.924450944753),\n", + " (6066.342554851419, 15516.34255485142),\n", + " (9433.42445094475, 26107.389504989897)]" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "all_seps[:10]" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": {}, + "outputs": [], + "source": [ + "min_map_arrays = np.array(all_seps).T\n", + "min_map_arrays[0:10]\n", + "mins = min_map_arrays[0]\n", + "maxs = min_map_arrays[1]" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3933/3933 [00:01<00:00, 2135.54it/s]\n" + ] + } + ], + "source": [ + "pix_a_order = []\n", + "pix_a_pixel = []\n", + "pix_b_order = []\n", + "pix_b_pixel = []\n", + "\n", + "for a in tqdm(range(0, num_partitions)):\n", + " for b in range(a, num_partitions):\n", + " pix_a_order.append(gaia_full_partition_list[a].order)\n", + " pix_a_pixel.append(gaia_full_partition_list[a].pixel)\n", + " pix_b_order.append(gaia_full_partition_list[b].order)\n", + " pix_b_pixel.append(gaia_full_partition_list[b].pixel)\n", + "\n", + "\n", + "big_beautiful_frame = pd.DataFrame({\"Norder_a\": pix_a_order,\n", + " \"Npix_a\": pix_a_pixel,\n", + " \"Norder_b\": pix_b_order,\n", + " \"Npix_b\": pix_b_pixel,\n", + " \"min_sep\": mins,\n", + " \"max_sep\": maxs,\n", + " })\n", + "big_beautiful_frame.to_csv(\"bbf.csv\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "hipscatenv", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/incubator/mask/min_max_hist.ipynb b/incubator/mask/min_max_hist.ipynb new file mode 100644 index 0000000..8241cf0 --- /dev/null +++ b/incubator/mask/min_max_hist.ipynb @@ -0,0 +1,263 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import hipscat\n", + "import healpy as hp\n", + "import pandas as pd\n", + "from tqdm import tqdm\n", + "from hipscat.inspection import plot_pixel_list\n", + "from hipscat.pixel_math import HealpixPixel\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAG2CAYAAAByJ/zDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLrklEQVR4nO3de1hU1f4/8PfM6HARB0RjBhKUylSUQAURK9PkROqxsI6pkZGalmcwjdL0l6J2s8xMK44ePRXW1/vJzNQwDqak4g1F8ZqeKKkjkCGMoALOrN8fxs7xgg57AzOb9+t59vO4915rzWfv800/33XbGiGEABEREZGKaBs6ACIiIiKlMcEhIiIi1WGCQ0RERKrDBIeIiIhUhwkOERERqQ4THCIiIlIdJjhERESkOkxwiIiISHWY4BAREZHqMMEhIiIi1WnQBCczMxMDBw5EQEAANBoN1q5da3dfCIHk5GT4+/vDw8MDMTExOHHihF2Z4uJixMfHw2AwwMfHB6NGjUJZWZldmYMHD+L++++Hu7s7AgMDMXv27GtiWb16NTp06AB3d3eEhoZi48aNij8vERER1Y8GTXDKy8sRFhaGlJSU696fPXs2PvjgAyxcuBC7du1Cs2bNEBsbi4sXL0pl4uPjcfjwYaSnp2P9+vXIzMzEmDFjpPsWiwUPPfQQ2rRpg+zsbLz77ruYMWMGFi1aJJXZsWMHhg0bhlGjRmH//v2Ii4tDXFwcDh06VHcPT0RERHVHOAkA4ssvv5TObTabMJlM4t1335WulZSUCDc3N7F8+XIhhBBHjhwRAMSePXukMt98843QaDTi119/FUII8Y9//EO0aNFCVFRUSGVeeeUV0b59e+n8iSeeEAMGDLCLJyoqSjz33HOKPiMRERHVjyYNnWDdSF5eHgoKChATEyNd8/b2RlRUFLKysjB06FBkZWXBx8cHERERUpmYmBhotVrs2rULgwYNQlZWFnr16gW9Xi+ViY2NxTvvvIOzZ8+iRYsWyMrKQlJSkt3vx8bGXjNkdqWKigpUVFRI5zabDcXFxWjZsiU0Go0Cb4CIiNRKCIFz584hICAAWm3dDaZcvHgRlZWVstvR6/Vwd3dXIKL647QJTkFBAQDAaDTaXTcajdK9goIC+Pn52d1v0qQJfH197coEBwdf00b1vRYtWqCgoKDG37meWbNmYebMmbV4MiIiosvy8/PRunXrOmn74sWLuM0jGGW48b9lt8pkMiEvL8+lkhynTXCc3ZQpU+x6fUpLSxEUFIT8v4XA0FR3+WJcaANFd33fd+/c0CHc0OEWtzd0CDXaU1k3fwEpZftx/4YOoUYPjA9o6BDoD5XuoqFDIABVlyz4encbNG/evM5+o7KyEmUowIvIhxsMtW6nAha8XxCIyspKJjhKMJlMAIDCwkL4+//5l3dhYSHCw8OlMkVFRXb1Ll26hOLiYqm+yWRCYWGhXZnq85uVqb5/PW5ubnBzc7vmuqGpDgb9HwmOp/6a+w2pWXPn/T9Md4NnQ4dQo6aVXg0dQo20XnX3l6QS9Lra/+VKyhJNmOA4k/qY0qDXGqDX1P6/QSEA2JSLp7447T44wcHBMJlMyMjIkK5ZLBbs2rUL0dHRAIDo6GiUlJQgOztbKrN582bYbDZERUVJZTIzM1FVVSWVSU9PR/v27dGiRQupzJW/U12m+neIiIhclU0n/3BFDZrglJWVIScnBzk5OQAuTyzOycnBqVOnoNFoMGHCBLzxxhtYt24dcnNz8fTTTyMgIABxcXEAgI4dO+Lhhx/G6NGjsXv3bmzfvh2JiYkYOnQoAgIud4k/+eST0Ov1GDVqFA4fPoyVK1di/vz5dsNL48ePR1paGt577z0cO3YMM2bMwN69e5GYmFjfr4SIiEhRSiU4kZGRCAkJueHWLs6mQYeo9u7diz59+kjn1UlHQkICUlNTMWnSJJSXl2PMmDEoKSnBfffdh7S0NLsxwKVLlyIxMRF9+/aFVqvF448/jg8++EC67+3tjW+//RZmsxndunVDq1atkJycbLdXTs+ePbFs2TJMnToV/+///T+0a9cOa9euRefOzjtnhYiIqD7t2bMHBoPrDDdrhBAckFWAxWKBt7c3SoeF/jkH52/hDRrT1bZE39PQIdzQQd/Ahg6hRjsrgxo6hBplHnXuSbx9xzj3JPLGpMKDf+U7g6pLFqzZ0QKlpaV1ljRU/7s0wasUbjLm4FQIC+aVeddprHXBaScZExERkXw2HWCTMZfZ5qI5sdNOMiYiIiKqLfbgEBERqZhNe/modX0XXCIOMMEhIiJSNaEDhIwER7jo14c4REVERESqwx4cIiIiFbPpZA5RsQeHiIiInA03+iMiIiK6AVfb6I8JDhERkYoJmUNUrjrJmAkOERGRill1AlZt7Xfrs2pcc6c/JjhEREQqJveL4DarcrHUJ04yJiIiItVhDw4REZGKye7BUS6UesUEh4iISMUaa4LDISoiIiJSHfbgEBERqZjQXT5qXV+5UOoVExwiIiIV4xAVERERkUqwB4eIiEjF2INDREREqmPTCth0Mo4/dkHmxzaJiIhIdfixTSIiInIaQuYQlXDRZVRMcIiIiFTMqgO0MhIcKxMcIiIicjayJxm7aILDScZERESkOuzBISIiUrHG2oPDBIeIiEjFhE5A6GqfpQgXnWXMISoiIiJSHfbgEBERqZhNK3OIykW3MmaCQ0REpGKy5+C4aILDISoiIiJSHfbgEBERqVhj7cFhgkNERKRiVp2AVsYqKquLrhPnEBURERGpDntwiIiIVEz2xzY5REVERETOxlXn4LRt2xYGgwFarRYtWrTAd99951B9JjhEREQqJnsfHKtysThqx44d8PLyqlVdzsEhIiIi1WGCQ0REpGLV36KSczgqMzMTAwcOREBAADQaDdauXXtNmZSUFLRt2xbu7u6IiorC7t277e5rNBo88MADiIyMxNKlSx2OgQkOERGRilXPwZFzOKq8vBxhYWFISUm57v2VK1ciKSkJ06dPx759+xAWFobY2FgUFRVJZbZt24bs7GysW7cOb731Fg4ePOhQDExwiIiI6KYsFovdUVFRccOy/fr1wxtvvIFBgwZd9/7cuXMxevRojBgxAiEhIVi4cCE8PT3xySefSGVuv/12AIC/vz/69++Pffv2ORQvExwiIiIVU6oHJzAwEN7e3tIxa9asWsVTWVmJ7OxsxMTESNe0Wi1iYmKQlZUF4HIP0Llz5wAAZWVl2Lx5Mzp16uTQ73AVFRERkYrZdAI2GTsZV9fNz8+HwWCQrru5udWqvTNnzsBqtcJoNNpdNxqNOHbsGACgsLBQ6v2xWq0YPXo0IiMjHfodJjhERER0UwaDwS7BqUt33HEHDhw4IKsNDlERERGpmFJDVJGRkQgJCbnhxOFb1apVK+h0OhQWFtpdLywshMlkktX2ldiDQ0REpGJWLaCR0Z1h/aPunj17FOnB0ev16NatGzIyMhAXFwcAsNlsyMjIQGJiouz2qzHBISIiIkWVlZXh5MmT0nleXh5ycnLg6+uLoKAgJCUlISEhAREREejevTvmzZuH8vJyjBgxQrEYmOAQERGpmNDKm2QstI7X3bt3L/r06SOdJyUlAQASEhKQmpqKIUOG4LfffkNycjIKCgoQHh6OtLS0ayYey8EEh4iISMVsOkAj51tUV8zB0el0MJvNMJvNNdbp3bs3hKg5MUpMTFR0SOpqTHCIiIhUTKkER6k5OPWFq6iIiIhIddiDQ0REpGK1/WDmlfVdEXtwiIiIVMymlX8Ayu2DU1/Yg0NEREQ35WpzcJjgEBERqZhSk4xdDRMcIiIiFbNqBSBjHo21FvvgOAPOwSEiIiLVYYJDRESkYs72sc36wiEqIiIiFbNpBTQyhplsf9R1tUnG7MEhIiIi1WEPDhERkYoJnbyVUIKrqIiIiMjZ2HTyVlHJ+RJ5Q2KCQ0REpGJCe/mQU98VOXXYVqsV06ZNQ3BwMDw8PHDnnXfi9ddft/sEuxACycnJ8Pf3h4eHB2JiYnDixAm7doqLixEfHw+DwQAfHx+MGjUKZWVldmUOHjyI+++/H+7u7ggMDMTs2bPr5RmJiIhcAVdRKeidd97BggULsGTJEnTq1Al79+7FiBEj4O3tjRdeeAEAMHv2bHzwwQdYsmQJgoODMW3aNMTGxuLIkSNwd3cHAMTHx+P06dNIT09HVVUVRowYgTFjxmDZsmUAAIvFgoceeggxMTFYuHAhcnNzMXLkSPj4+GDMmDEN9vxERERyKTVE5WqrqJw6wdmxYwceffRRDBgwAADQtm1bLF++HLt37wZwufdm3rx5mDp1Kh599FEAwGeffQaj0Yi1a9di6NChOHr0KNLS0rBnzx5EREQAAD788EP0798fc+bMQUBAAJYuXYrKykp88skn0Ov16NSpE3JycjB37lwmOERE5NJsGsgar7FpFAulXjn1EFXPnj2RkZGBH374AQBw4MABbNu2Df369QMA5OXloaCgADExMVIdb29vREVFISsrCwCQlZUFHx8fKbkBgJiYGGi1WuzatUsq06tXL+j1eqlMbGwsjh8/jrNnz143toqKClgsFruDiIiInINT9+BMnjwZFosFHTp0gE6ng9VqxZtvvon4+HgAQEFBAQDAaDTa1TMajdK9goIC+Pn52d1v0qQJfH197coEBwdf00b1vRYtWlwT26xZszBz5kwFnpKIiKju2HQAGuHHNp26B2fVqlVYunQpli1bhn379mHJkiWYM2cOlixZ0tChYcqUKSgtLZWO/Pz8hg6JiIjoGlatkH24IqfuwZk4cSImT56MoUOHAgBCQ0Px888/Y9asWUhISIDJZAIAFBYWwt/fX6pXWFiI8PBwAIDJZEJRUZFdu5cuXUJxcbFU32QyobCw0K5M9Xl1mau5ubnBzc1N/kMSERGR4py6B+f8+fPQau1D1Ol0sNlsAIDg4GCYTCZkZGRI9y0WC3bt2oXo6GgAQHR0NEpKSpCdnS2V2bx5M2w2G6KioqQymZmZqKqqksqkp6ejffv21x2eIiIichX82KYTGjhwIN58800EBQWhU6dO2L9/P+bOnYuRI0cCADQaDSZMmIA33ngD7dq1k5aJBwQEIC4uDgDQsWNHPPzwwxg9ejQWLlyIqqoqJCYmYujQoQgICAAAPPnkk5g5cyZGjRqFV155BYcOHcL8+fPx/vvvN9SjExERKcKmFUAj/NimUyc4H374IaZNm4a///3vKCoqQkBAAJ577jkkJydLZSZNmoTy8nKMGTMGJSUluO+++5CWlibtgQMAS5cuRWJiIvr27QutVovHH38cH3zwgXTf29sb3377LcxmM7p164ZWrVohOTmZS8SJiIhclEZcuS0w1ZrFYoG3tzdKh4XCoP+jP+9v4Q0a09W2RN/T0CHc0EHfwIYOoUY7K4MaOoQaZR4NaOgQatR3zO0NHQL9ocKDf+U7g6pLFqzZ0QKlpaV11itS/e+Sz4J8aDxq/xviggUlYwPrNNa64NQ9OERERCSPkDlEJbiKioiIiJyNTSugkfGpBldNcJx6FRURERFRbbAHh4iISMVsWkAjoztDuGhXCBMcIiIiFbNqBTSNcA6Oi+ZlREREVJ+40R8RERE5DZsO0Mj4YKb4oy43+iMiIiKnYeMQFREREZE6sAeHiIhIxbiKioiIiFTn8hwcGUNUMubvNCQXzcuIiIiIbow9OErLLQR0l/PG4rHesprad/cdSkQk+U+LDoq2p6R1J9s1dAg1urTa2NAh1KjPOveGDqFGrX7WNHQI9If3Cvm/hTOwWLRYI++fiFunFfImCrvoJGMmOERERCqm0cicg6MBXDHFYYJDRESkYloFlolbFYynvnAODhEREakOe3CIiIhUTKvAMnFX7MFhgkNERKRiOp2QuUxcoErBeOoLh6iIiIjopvixTSIiInIaSkwyBvixTSIiInIiSszBcUUuGjYRERHRjbEHh4iISMWUGqJyNUxwiIiIVEwjM8Fx1U81cIiKiIiIVIc9OERERCrWWCcZM8EhIiJSMY1GQCtjmMmmcc0hKiY4REREKqbTCGhlJCkaF01wXLTjiYiIiOjG2INDRESkYlrt5aP2DSgWSr1igkNERKRiWq28OThcJk5ERETkJNiDQ0REpGKNtQeHCQ4REZGKNdYEh0NURERE5JTOnz+PNm3a4OWXX3a4LntwiIiIVMyVV1G9+eab6NGjR63qsgeHiIhIxTQaIftoCCdOnMCxY8fQr1+/WtVngkNERESKyszMxMCBAxEQEACNRoO1a9deUyYlJQVt27aFu7s7oqKisHv3brv7L7/8MmbNmlXrGJjgEBERqVj1JGM5h6PKy8sRFhaGlJSU695fuXIlkpKSMH36dOzbtw9hYWGIjY1FUVERAOCrr77C3XffjbvvvrvWz805OERERCqmkbmKSvxR12Kx2F13c3ODm5vbdev069evxqGluXPnYvTo0RgxYgQAYOHChdiwYQM++eQTTJ48GTt37sSKFSuwevVqlJWVoaqqCgaDAcnJybccN3twiIiIVEwHQKeRcfzRTmBgILy9vaWjtsNHlZWVyM7ORkxMjHRNq9UiJiYGWVlZAIBZs2YhPz8fP/30E+bMmYPRo0c7lNwA7MEhIiKiW5Cfnw+DwSCd36j35mbOnDkDq9UKo9Fod91oNOLYsWOyYrwSExwiIiIVk7vRX/UQlcFgsEtw6sszzzxTq3ocoiIiIlIxrUbIPgAgMjISISEhN5w4fKtatWoFnU6HwsJCu+uFhYUwmUyy2r4Se3CIiIjopvbs2aNID45er0e3bt2QkZGBuLg4AIDNZkNGRgYSExNlt1+NCQ4REZGKKTVE5YiysjKcPHlSOs/Ly0NOTg58fX0RFBSEpKQkJCQkICIiAt27d8e8efNQXl4urapSAhMcIiIiFdNqIA0z1YbQOF5n79696NOnj3SelJQEAEhISEBqaiqGDBmC3377DcnJySgoKEB4eDjS0tKumXgsBxMcIiIiuqnIyEjodDqYzWaYzeYay/bu3RtC1JxUJSYmKjokdTUmOERERCqm0V4+5NQHlJuDU1+Y4BAREanYlSuhakM00Mc25eIycSIiIlIdJjhEREQqptUK6GQc1SuwlNoHp744PERltVqRmpqKjIwMFBUVwWaz2d3fvHmzYsERERGRPHKHqKrrqn4Ozvjx45GamooBAwagc+fO0GhqsX6MiIiI6oVSCY6rcTjBWbFiBVatWoX+/fvXRTxEREREsjmc4Oj1etx11111EQsREREpTO5OxnLqNiSHJxm/9NJLmD9//k038CEiIqKGd3knY3kH0AgmGW/btg3fffcdvvnmG3Tq1AlNmza1u79mzRrFgiMiIiLnoPpJxj4+Phg0aFBdxEJEREQK4yTjW/Tpp5/WRRxERERUBzQQ0MhIUjRoJAlOtd9++w3Hjx8HALRv3x633XabYkERERERyeHwJOPy8nKMHDkS/v7+6NWrF3r16oWAgACMGjUK58+fr4sYiYiIqJaqh6jkHK7I4QQnKSkJW7duxddff42SkhKUlJTgq6++wtatW/HSSy/VRYxERERUS0olOK62isrhBOeLL77Axx9/jH79+sFgMMBgMKB///5YvHgx/v3vfyse4K+//oqnnnoKLVu2hIeHB0JDQ7F3717pvhACycnJ8Pf3h4eHB2JiYnDixAm7NoqLixEfHw+DwQAfHx+MGjUKZWVldmUOHjyI+++/H+7u7ggMDMTs2bMVfxYiIqL6ptS3qPbs2YMjR47AbDY38BPdGocTnPPnz8NoNF5z3c/PT/EhqrNnz+Lee+9F06ZN8c033+DIkSN477330KJFC6nM7Nmz8cEHH2DhwoXYtWsXmjVrhtjYWFy8eFEqEx8fj8OHDyM9PR3r169HZmYmxowZI923WCx46KGH0KZNG2RnZ+Pdd9/FjBkzsGjRIkWfh4iIiOqHw5OMo6OjMX36dHz22Wdwd3cHAFy4cAEzZ85EdHS0osG98847CAwMtFu5FRwcLP1ZCIF58+Zh6tSpePTRRwEAn332GYxGI9auXYuhQ4fi6NGjSEtLw549exAREQEA+PDDD9G/f3/MmTMHAQEBWLp0KSorK/HJJ59Ar9ejU6dOyMnJwdy5c+0SISIiIlfTWJeJO9yDM3/+fGzfvh2tW7dG37590bdvXwQGBmLHjh2YP3++osGtW7cOERERGDx4MPz8/NClSxcsXrxYup+Xl4eCggLExMRI17y9vREVFYWsrCwAQFZWFnx8fKTkBgBiYmKg1Wqxa9cuqUyvXr2g1+ulMrGxsTh+/DjOnj173dgqKipgsVjsDiIiImfDSca3qHPnzjhx4gRmzZqF8PBwhIeH4+2338aJEyfQqVMnRYP78ccfsWDBArRr1w6bNm3C2LFj8cILL2DJkiUAgIKCAgC4ZsjMaDRK9woKCuDn52d3v0mTJvD19bUrc702rvyNq82aNQve3t7SERgYKPNpiYiISCm12gfH09MTo0ePVjqWa9hsNkREROCtt94CAHTp0gWHDh3CwoULkZCQUOe/X5MpU6YgKSlJOrdYLExyiIjI6WghoJWxWZ+cug3plhKcdevWoV+/fmjatCnWrVtXY9lHHnlEkcAAwN/fHyEhIXbXOnbsiC+++AIAYDKZAACFhYXw9/eXyhQWFiI8PFwqU1RUZNfGpUuXUFxcLNU3mUwoLCy0K1N9Xl3mam5ubnBzc6vlkxEREdUPpebgREZGQqfTwWw2u8RKqltKcOLi4qShnri4uBuW02g0sFqtSsWGe++9V9otudoPP/yANm3aALg84dhkMiEjI0NKaCwWC3bt2oWxY8cCuDwpuqSkBNnZ2ejWrRsAYPPmzbDZbIiKipLKvPrqq6iqqpI+Hpqeno727dvbrdgiIiJqrFztY5u3NAfHZrNJ81hsNtsNDyWTGwB48cUXsXPnTrz11ls4efIkli1bhkWLFkmZo0ajwYQJE/DGG29g3bp1yM3NxdNPP42AgAApEevYsSMefvhhjB49Grt378b27duRmJiIoUOHIiAgAADw5JNPQq/XY9SoUTh8+DBWrlyJ+fPn2w1BERERuSKNRsg+XJHDk4w/++wzVFRUXHO9srISn332mSJBVYuMjMSXX36J5cuXo3Pnznj99dcxb948xMfHS2UmTZqEcePGYcyYMYiMjERZWRnS0tKkJewAsHTpUnTo0AF9+/ZF//79cd9999ntcePt7Y1vv/0WeXl56NatG1566SUkJydziTgREbk8rUb+4Yo0QgiHUjOdTofTp09fszLp999/h5+fn+K9OK7CYrHA29sbpZ39YNBdzhuL/zlEVpv77r5DidAk/2nRQdH2lLTuZLuGDqFGl1Zfu7mlM4la537zQg3I70eH/38pqiPvFd68DNU96d+M0tI6G/ap/o1HC75DU4NXrdupspThK1OfOo21Lji8ikoIAY3m2nTul19+gbe3tyJBERERkTJ0GgGdjGEmm4sOUd1ygtOlSxdoNBpoNBr07dsXTZr8WdVqtSIvLw8PP/xwnQRJREREtdNYdzK+5QSnetJuTk4OYmNj4eX1Z3eXXq9H27Zt8fjjjyseIBEREdWeFvL2snHVAeZbTnCmT58OAGjbti2GDBliN4mXiIiIyJk4nJglJCQwuSEiInIRGpnfoapeJh4ZGYmQkBCkpKQ08BPdGocnGVutVrz//vtYtWoVTp06hcrKSrv7xcXFigVHRERE8ig1B0eVG/1daebMmZg7dy6GDBmC0tJSJCUl4bHHHoNWq8WMGTPqIEQiIiIixzic4CxduhSLFy/GSy+9hCZNmmDYsGH417/+heTkZOzcubMuYiQiIqJaqv7YppzDFTmc4BQUFCA0NBQA4OXlhdLSUgDAX//6V2zYsEHZ6IiIiEgWfqrhFrVu3RqnT58GANx555349ttvAVwem+PXtYmIiMgZOJzgDBo0CBkZGQCAcePGYdq0aWjXrh2efvppjBw5UvEAiYiIqPYa6xCVw6uo3n77benPQ4YMQZs2bbBjxw60a9cOAwcOVDQ4IiIikoc7Gd+CqqoqPPfcc5g2bRqCg4MBAD169ECPHj3qJDgiIiKi2nBoiKpp06b44osv6ioWIiIiUpgOQvbhihyegxMXF4e1a9fWQShERESkNDm7GMsd3mpIDs/BadeuHV577TVs374d3bp1Q7Nmzezuv/DCC4oFR0RERPLInSjcaCYZf/zxx/Dx8UF2djays7Pt7mk0GiY4RERE1OAcTnDy8vLqIg4iIiKqAxqZw0xXfmxTp9PBbDbDbDYrFV6dcTjBISIiItehhQ1a2GTVB1zvY5u1SnB++eUXrFu37rpfE587d64igRERERHVlsMJTkZGBh555BHccccdOHbsGDp37oyffvoJQgh07dq1LmIkIiKiWtJAQCNjorCcug3J4WXiU6ZMwcsvv4zc3Fy4u7vjiy++QH5+Ph544AEMHjy4LmIkIiKiWmqsn2pwOME5evQonn76aQBAkyZNcOHCBXh5eeG1117DO++8o3iARERERI5yOMFp1qyZNO/G398f//3vf6V7Z86cUS4yIiIikk2nEbIPV+TwHJwePXpg27Zt6NixI/r374+XXnoJubm5WLNmDb9JRURE5GS40d8tmjt3LsrKygAAM2fORFlZGVauXIl27dpxBRURERE5BYcTnDvuuEP6c7NmzbBw4UJFAyIiIiLlsAfHQXv37sXRo0cBACEhIejWrZtiQREREZEytJCXpDg8WddJOJzg/PLLLxg2bBi2b98OHx8fAEBJSQl69uyJFStWoHXr1krHSERERLXUWHtwHE7Mnn32WVRVVeHo0aMoLi5GcXExjh49CpvNhmeffbYuYiQiIiJyiMM9OFu3bsWOHTvQvn176Vr79u3x4Ycf4v7771c0OJd0RwugqQ4AcN7DTVZT55vqlYhIcuK8r6LtKakos2VDh1CjmDT3hg6hRqdCrQ0dQo0+z3LVTm4i16eR2YPjqjsZO5zgBAYGoqqq6prrVqsVAQEBigRFREREyuAQ1S169913MW7cOOzdu1e6tnfvXowfPx5z5sxRNDgiIiKi2nC4B+eZZ57B+fPnERUVhSZNLle/dOkSmjRpgpEjR2LkyJFS2eLiYuUiJSIiIoe54sc2S0pKEBMTg0uXLuHSpUsYP348Ro8e7VAbDic48+bNc7QKERERNRAtbNDCJqt+fWvevDkyMzPh6emJ8vJydO7cGY899hhatrz1+ZoOJzgJCQmOViEiIiK6ZTqdDp6engCAiooKCCEghGM9SbVa2vDf//4XU6dOxbBhw1BUVAQA+Oabb3D48OHaNEdERER1pHqSsZzDUZmZmRg4cCACAgKg0Wiwdu3aa8qkpKSgbdu2cHd3R1RUFHbv3m13v6SkBGFhYWjdujUmTpyIVq1aOfjcDtq6dStCQ0Oxa9curFmzRvou1YEDBzB9+nRHmyMiIqI6pBNC9uGo8vJyhIWFISUl5br3V65ciaSkJEyfPh379u1DWFgYYmNjpU4TAPDx8cGBAweQl5eHZcuWobCw0KEYHE5wJk+ejDfeeAPp6enQ6//cp+XBBx/Ezp07HW2OiIiIXIDFYrE7Kioqbli2X79+eOONNzBo0KDr3p87dy5Gjx6NESNGICQkBAsXLoSnpyc++eSTa8oajUaEhYXh+++/dyhehxOc3Nzc6wbs5+eHM2fOONocERER1SGlhqgCAwPh7e0tHbNmzapVPJWVlcjOzkZMTMyfMWq1iImJQVZWFgCgsLAQ586dAwCUlpYiMzPTboPhW+HwJGMfHx+cPn0awcHBdtf379+P22+/3dHmiIiIqA4ptdFffn4+DAaDdN3NrXa79Z85cwZWqxVGo9HuutFoxLFjxwAAP//8M8aMGSNNLh43bhxCQ0Md+h2HE5yhQ4filVdewerVq6HRaGCz2bB9+3a8/PLLePrppx1tjoiIiOqQUgmOwWCwS3DqUvfu3ZGTkyOrDYeHqN566y106NABgYGBKCsrQ0hICHr16oWePXti6tSpsoIhIiIi5xQZGYmQkJAbThy+Va1atYJOp7tm0nBhYSFMJpOstq/kcA+OXq/H4sWLkZycjNzcXJSVlaFLly5o166dYkERERGRMjTCBq2o/WZ9mj/q7tmzR5EeHL1ej27duiEjIwNxcXEAAJvNhoyMDCQmJspuv5rDCU61wMBABAYGwmq1Ijc3F2fPnkWLFi0UC4yIiIjka4iPbZaVleHkyZPSeV5eHnJycuDr64ugoCAkJSUhISEBERER6N69O+bNm4fy8nKMGDGi1nFezeEEZ8KECQgNDcWoUaNgtVrxwAMPYMeOHfD09MT69evRu3dvxYIjIiIi17N371706dNHOk9KSgJw+WsIqampGDJkCH777TckJyejoKAA4eHhSEtLu2bisRwOJzj//ve/8dRTTwEAvv76a/z44484duwYPv/8c7z66qvYvn27YsERERGRPFohoK3FZn1X1gcuz8HR6XQwm80wm8011undu/dNP62QmJio6JDU1RxOcM6cOSNNAtq4cSOeeOIJ3H333Rg5ciTmz5+veIBERERUe0oNUSk1B6e+OLyKymg04siRI7BarUhLS8Nf/vIXAMD58+eh0+kUD5CIiIjIUQ734IwYMQJPPPEE/P39odFopJ0Id+3ahQ4dOigeIBEREdWeUkNUrsbhBGfGjBno3Lkz8vPzMXjwYGknQ51Oh8mTJyseIBEREdWeUkNUjszBcQa1Wib+t7/97ZprCQkJsoMhIiIi5+Rqc3BqvQ8OEREROb/LQ1S13+iv0QxRERERketoiI3+nAETHCIiIhVrrJOMHV4mTkRERI2PUh/brC8OJzipqanXvX7p0iVMmTJFbjxERESkJCGgkXFA/LnR35EjR1xiBRVQiwTnhRdewODBg3H27Fnp2vHjxxEVFYXly5crGhwRERHJo4VN9uGKHE5w9u/fj19++QWhoaFIT09HSkoKunbtig4dOuDAgQN1ESMRERGRQxyeZHznnXdi+/btmDBhAh5++GHodDosWbIEw4YNq4v4iIiISAZOMnbAhg0bsGLFCkRHR8PHxwcff/wx/ve//ykdGxEREcmkE0L24YocTnCee+45DB48GK+88gq+//57HDx4EHq9HqGhoVi1alVdxEhEREQNzNVWUTk8RLV9+3bs2rULYWFhAACTyYSNGzciJSUFI0eOxBNPPKF4kERERFQ7Sg1RudqnGhzuwcnOzpaSmyuZzWZkZ2crEtSNvP3229BoNJgwYYJ07eLFizCbzWjZsiW8vLzw+OOPo7Cw0K7eqVOnMGDAAHh6esLPzw8TJ07EpUuX7Mps2bIFXbt2hZubG+66664bLocnIiJyJVphk324IocTnOqvh19P+/btZQVTkz179uCf//wn7rnnHrvrL774Ir7++musXr0aW7duxf/+9z889thj0n2r1YoBAwagsrISO3bswJIlS5Camork5GSpTF5eHgYMGIA+ffogJycHEyZMwLPPPotNmzbV2fMQERFR3anVpxr+/e9/Y9WqVTh16hQqKyvt7u3bt0+RwK5UVlaG+Ph4LF68GG+88YZ0vbS0FB9//DGWLVuGBx98EADw6aefomPHjti5cyd69OiBb7/9FkeOHMF//vMfGI1GhIeH4/XXX8crr7yCGTNmQK/XY+HChQgODsZ7770HAOjYsSO2bduG999/H7GxsYo/DxERUX3hKqpb9MEHH2DEiBEwGo3Yv38/unfvjpYtW+LHH39Ev3796iJGmM1mDBgwADExMXbXs7OzUVVVZXe9Q4cOCAoKQlZWFgAgKysLoaGhMBqNUpnY2FhYLBYcPnxYKnN127GxsVIb11NRUQGLxWJ3EBEROZvqj23KOVyRwz04//jHP7Bo0SIMGzYMqampmDRpEu644w4kJyejuLhY8QBXrFiBffv2Yc+ePdfcKygogF6vh4+Pj911o9GIgoICqcyVyU31/ep7NZWxWCy4cOECPDw8rvntWbNmYebMmdcGHH474H75tf5sbHVrD3kDB71ay6p/tTUb2yranpKeSWne0CHUqMTk3GPQWxc1begQiMhJsQfnFp06dQo9e/YEAHh4eODcuXMAgOHDhyv+qYb8/HyMHz8eS5cuhbu7u6JtyzVlyhSUlpZKR35+fkOHREREVGdUv0zcZDKhuLgYbdq0QVBQEHbu3ImwsDDk5eVBKJzlZWdno6ioCF27dpWuWa1WZGZm4qOPPsKmTZtQWVmJkpISu16cwsJCmEwmKd7du3fbtVu9yurKMlevvCosLITBYLhu7w1webJ1TROuiYiInIFG5koozR91Vb9M/MEHH8S6desAACNGjMCLL76Iv/zlLxgyZAgGDRqkaHB9+/ZFbm4ucnJypCMiIgLx8fHSn5s2bYqMjAypzvHjx3Hq1ClER0cDAKKjo5Gbm4uioiKpTHp6OgwGA0JCQqQyV7ZRXaa6DSIiIleltQnZhytyuAdn0aJFsNkuZ3PV+8/s2LEDjzzyCJ577jlFg2vevDk6d+5sd61Zs2Zo2bKldH3UqFFISkqCr68vDAYDxo0bh+joaPTo0QMA8NBDDyEkJATDhw/H7NmzUVBQgKlTp8JsNks9MM8//zw++ugjTJo0CSNHjsTmzZuxatUqbNiwQdHnISIiovrhcIKj1Wqh1f7Z8TN06FAMHTpU0aAc8f7770Or1eLxxx9HRUUFYmNj8Y9//EO6r9PpsH79eowdOxbR0dFo1qwZEhIS8Nprr0llgoODsWHDBrz44ouYP38+WrdujX/9619cIk5ERC5PK+RNFNa6ZgdO7fbBuXjxIg4ePIiioiKpN6faI488okhgN7Jlyxa7c3d3d6SkpNQ46alNmzbYuHFjje327t0b+/fvVyJEIiIipyF3N2JX3cnY4QQnLS0NTz/9NM6cOXPNPY1GA6vVqkhgRERERLXl8CTjcePGYfDgwTh9+jRsNpvdweSGiIjIuWghpL1wanU0lo3+CgsLkZSUdM3GeEREROR8uNHfLfrb3/52zTwYIiIiUjfVb/T30UcfYfDgwfj+++8RGhqKpk3tt4h/4YUXFAuOiIiI5FGqB8fVNvpzOMFZvnw5vv32W7i7u2PLli3QaDTSPY1GwwSHiIjIiWhsNmhsMnYyllG3ITmc4Lz66quYOXMmJk+ebLcfDhERETkfzsG5RZWVlRgyZAiTGyIiInJaDmcpCQkJWLlyZV3EQkRERAqTtURcZu9PQ3J4iMpqtWL27NnYtGkT7rnnnmsmGc+dO1ex4IiIiEgerc0GnYx5NNrGMgcnNzcXXbp0AQAcOnTI7t6VE46JiIiIGorDCc53331XF3EQERFRHWisk4xr9bFNIiIicg1am4DWJiPBkVG3IXEpFBEREakOe3CIiIhUTCts0AoZk4xl1G1ITHCIiIhUjENURERERDeg+o9tEhERkeu4vIpKzhBVI/nYJhEREbmOxjpExQSHiIhIxbSQuQ8OXDPB4RwcIiIiUh324BAREamYxmaT9T0pTWP5FhURERG5jsY6B4dDVERERKQ67MEhIiJSMX5sk4iIiFRHK3MOjpy6DYlDVERERKQ67MEhIiJSscY6yZgJDhERkYppbDZZS71ddZk4h6iIiIhIddiDo7DKrkGobKYHAGQa75bV1le/yqt/tedebKloe0q6pG/oCGr25QZdQ4dARFQrHKIiIiIi1WGCQ0RERKqjEzboRO3n0cip25A4B4eIiIicSn5+Pnr37o2QkBDcc889WL16tcNtsAeHiIhIxTRC3hCVpgF2Mm7SpAnmzZuH8PBwFBQUoFu3bujfvz+aNWt2623UYXxERETUwC7PwZGzk3H9Jzj+/v7w9/cHAJhMJrRq1QrFxcUOJTgcoiIiIiJFZWZmYuDAgQgICIBGo8HatWuvKZOSkoK2bdvC3d0dUVFR2L1793Xbys7OhtVqRWBgoEMxMMEhIiJSsepVVHIOALBYLHZHRUXFDX+zvLwcYWFhSElJue79lStXIikpCdOnT8e+ffsQFhaG2NhYFBUV2ZUrLi7G008/jUWLFjn+3A7XICIiIpehVIITGBgIb29v6Zg1a9YNf7Nfv3544403MGjQoOvenzt3LkaPHo0RI0YgJCQECxcuhKenJz755BOpTEVFBeLi4jB58mT07NnT4efmHBwiIiK6qfz8fBgMBunczc2tVu1UVlYiOzsbU6ZMka5ptVrExMQgKysLACCEwDPPPIMHH3wQw4cPr9XvMMEhIiJSMa3NJnOS8eW6BoPBLsGprTNnzsBqtcJoNNpdNxqNOHbsGABg+/btWLlyJe655x5p/s7nn3+O0NDQW/4dJjhEREQq5oo7Gd93332wyfzIJ+fgEBER0U1FRkYiJCTkhhOHb1WrVq2g0+lQWFhod72wsBAmk0lW21diDw4REZGKaa02aK0yhqj+qLtnzx5Fhqj0ej26deuGjIwMxMXFAQBsNhsyMjKQmJgou/1qTHCIiIhUrCGGqMrKynDy5EnpPC8vDzk5OfD19UVQUBCSkpKQkJCAiIgIdO/eHfPmzUN5eTlGjBhR6zivxgSHiIhIxZSaZBwZGQmdTgez2Qyz2Vxjnb1796JPnz7SeVJSEgAgISEBqampGDJkCH777TckJyejoKAA4eHhSEtLu2bisRxMcIiIiOimHBmi6t27N8RNvmGVmJio6JDU1ZjgEBERqZhW5sc2tQ3wsU0lMMEhIiJSM5vt8iGnvgviMnEiIiJSHSY4REREamYV8g8otw9OfeEQFRERkZrZxOVDTn0otw9OfWEPDhEREakOe3CIiIjUrJFOMmaCQ0REpGZWSPNoal3fBXGIioiIiG6Kk4yJiIjIeSg0ROVqk4yZ4BAREamZQquoXA0THCIiIjWz2QBr45tkzDk4REREpDpMcIiIiNSseohKzgFOMiYiIiJnYpU5RGV1zUnG7MEhIiIi1XHqBGfWrFmIjIxE8+bN4efnh7i4OBw/ftyuzMWLF2E2m9GyZUt4eXnh8ccfR2FhoV2ZU6dOYcCAAfD09ISfnx8mTpyIS5cu2ZXZsmULunbtCjc3N9x1111ITU2t68cjIiKqewoNUbkap05wtm7dCrPZjJ07dyI9PR1VVVV46KGHUF5eLpV58cUX8fXXX2P16tXYunUr/ve//+Gxxx6T7lutVgwYMACVlZXYsWMHlixZgtTUVCQnJ0tl8vLyMGDAAPTp0wc5OTmYMGECnn32WWzatKlen5eIiEhx1UNUcg4X5NRzcNLS0uzOU1NT4efnh+zsbPTq1QulpaX4+OOPsWzZMjz44IMAgE8//RQdO3bEzp070aNHD3z77bc4cuQI/vOf/8BoNCI8PByvv/46XnnlFcyYMQN6vR4LFy5EcHAw3nvvPQBAx44dsW3bNrz//vuIjY2t9+cmIiIieZy6B+dqpaWlAABfX18AQHZ2NqqqqhATEyOV6dChA4KCgpCVlQUAyMrKQmhoKIxGo1QmNjYWFosFhw8flspc2UZ1meo2rqeiogIWi8XuICIicjqNdIjKqXtwrmSz2TBhwgTce++96Ny5MwCgoKAAer0ePj4+dmWNRiMKCgqkMlcmN9X3q+/VVMZiseDChQvw8PC4Jp5Zs2Zh5syZ11z/tnsYPA2Xy3+4L6IWT/qnvz4VIKv+1W77SaNoe0p6vaKhIyAiUimFVlFFRkZCp9PBbDbDbDYrFFzdcZkEx2w249ChQ9i2bVtDhwIAmDJlCpKSkqRzi8WCwMDABoyIiIio7rjaMnGXSHASExOxfv16ZGZmonXr1tJ1k8mEyspKlJSU2PXiFBYWwmQySWV2795t1171Kqsry1y98qqwsBAGg+G6vTcA4ObmBjc3N9nPRkREVKeEzGEm4ZpDVE49B0cIgcTERHz55ZfYvHkzgoOD7e5369YNTZs2RUZGhnTt+PHjOHXqFKKjowEA0dHRyM3NRVFRkVQmPT0dBoMBISEhUpkr26guU90GERGRy+IqKudjNpuxbNkyfPXVV2jevLk0Z8bb2xseHh7w9vbGqFGjkJSUBF9fXxgMBowbNw7R0dHo0aMHAOChhx5CSEgIhg8fjtmzZ6OgoABTp06F2WyWemCef/55fPTRR5g0aRJGjhyJzZs3Y9WqVdiwYUODPTsREZEiGunXxJ26B2fBggUoLS1F79694e/vLx0rV66Uyrz//vv461//iscffxy9evWCyWTCmjVrpPs6nQ7r16+HTqdDdHQ0nnrqKTz99NN47bXXpDLBwcHYsGED0tPTERYWhvfeew//+te/uESciIjIRTl1D464hXE/d3d3pKSk1PjxrzZt2mDjxo01ttO7d2/s37/f4RiJiIicmlXIXEXlmj04Tp3gEBERkUxWIS9JcdEEx6mHqIiIiIhqgz04REREamazXT7k1HdB7MEhIiJSs+ohKjkHLu9kHBISUuOcV2fCHhwiIiK6Ke5kTERERM7DJnOzPhcdomKCQ0REpGZcRUVERESkDuzBISIiUrNGuoqKCQ4REZGaNdIhKiY4REREaib3i+Au+jVxzsEhIiIi1WEPDhERkZpxiIqIiIhUp5Hug8MhKiIiIlId9uAQERGpmU1cPuTUd0HswSEiIlIzq/hzJVWtDn5sk4iIiFSKH9skIiIi52G1AVaNvPouiAkOERGRmjXSZeKcg0NERESqwx4cIiIiNeMQFREREalOIx2iYoJDRESkZjaZPTjcyZiIiIjIObAHh4iISM2sAtByiIqIiIjUxGoDtI1vkjGHqIiIiEh12INDRESkZhyiIiIiItVppENUTHAUtlHfCXp9MwBA/4QAWW01qVAioj+9rnB7REREzooJDhERkZo10iEqTjImIiJSM5v443MNtTxsDZPgDBo0CC1atMDf/va3WtVngkNEREROZ/z48fjss89qXZ8JDhERkZrJ6b2pPhpA79690bx581rXZ4JDRESkZtUf25RzOCgzMxMDBw5EQEAANBoN1q5de02ZlJQUtG3bFu7u7oiKisLu3bsVeNg/McEhIiJSswbowSkvL0dYWBhSUlKue3/lypVISkrC9OnTsW/fPoSFhSE2NhZFRUVyn1bCVVRERER0UxaLxe7czc0Nbm5u1y3br18/9OvX74ZtzZ07F6NHj8aIESMAAAsXLsSGDRvwySefYPLkyYrEyx4cIiIiNVNoiCowMBDe3t7SMWvWrFqFU1lZiezsbMTExEjXtFotYmJikJWVpcgjA+zBISIiUjerDZCxkXH1EFV+fj4MBoN0+Ua9Nzdz5swZWK1WGI1Gu+tGoxHHjh2TzmNiYnDgwAGUl5ejdevWWL16NaKjo2/5d5jgEBER0U0ZDAa7BKeu/ec//5FVn0NUREREamaTOTz1x0Z/kZGRCAkJueHE4VvVqlUr6HQ6FBYW2l0vLCyEyWSS1faV2INDRESkZnL3sfmj/p49exTpwdHr9ejWrRsyMjIQFxcHALDZbMjIyEBiYqLs9qsxwSEiIiJFlZWV4eTJk9J5Xl4ecnJy4Ovri6CgICQlJSEhIQERERHo3r075s2bh/LycmlVlRKY4BAREamZVQCo349t7t27F3369JHOk5KSAAAJCQlITU3FkCFD8NtvvyE5ORkFBQUIDw9HWlraNROP5WCCQ0REpGYKDVFFRkZCp9PBbDbDbDbXWKV3794QoubEKDExUdEhqasxwSEiIqKbUmoOTn1hgkNERKRmCvXguBomOAorGRaBpk0vZ7jhh+TsrARMlzFkSkREBKBB5uA4A+6DQ0REpGZWAVyy1f6wKrsPTn1hDw4RERHdFOfgEBERkfOwClkjVNU7GbsaJjhERERqZrUBQsacUBdNcDgHh4iIiFSHCQ4REZGayfnQZvUBTjImIiIiZ6LQEJWrTTJmDw4RERGpDntwiIiI1KyRTjJmgkNERKRmNpnLxG/y0UxnxSEqIiIiUh0mOERERGpmtck/wFVURERE5EwuCUDOt5+Fa66iYoJDRESkZlYboJGR4bjoHBwmOAo718qGJvrL3XnTha6BoyEiImqcmOAQERGpmVWZISpXwwSHiIhIzWwC8taJuyYmOAoRf2S4S2dfgMHQFABgsVxoyJCIiMhJWSwWAH/+21Gnv9XA9RsKExyFnDt3DgAQGBjYwJEQEZGr+P333+Ht7V0nbev1ephMJgQWFMhuy2AwoHv37tBqtTCbzTCbzQpEWLc0oj7Sx0bAZrPhf//7H5o3bw6NnNnqjZDFYkFgYCDy8/Ndagmis+D7qz2+O3n4/mqvtLQUQUFBOHv2LHx8fOrsdy5evIjKykrZ7ej1eri7uysQUf1hD45CtFotWrdu3dBhuDSDwcC/JGXg+6s9vjt5+P5qT6ut2/123d3dXS4xUQp3MiYiIiLVYYJDREREqsMEhxqcm5sbpk+fDjc3t4YOxSXx/dUe3508fH+1x3dX9zjJmIiIiFSHPThERESkOkxwiIiISHWY4BAREZHqMMEhIiIi1WGCQ4qYNWsWIiMj0bx5c/j5+SEuLg7Hjx+3K3Px4kWYzWa0bNkSXl5eePzxx1FYWGhX5tSpUxgwYAA8PT3h5+eHiRMn4tKlS3ZltmzZgq5du8LNzQ133XUXUlNT6/rx6tXbb78NjUaDCRMmSNf47mr266+/4qmnnkLLli3h4eGB0NBQ7N27V7ovhEBycjL8/f3h4eGBmJgYnDhxwq6N4uJixMfHw2AwwMfHB6NGjUJZWZldmYMHD+L++++Hu7s7AgMDMXv27Hp5vrpitVoxbdo0BAcHw8PDA3feeSdef/11u+8j8d39KTMzEwMHDkRAQAA0Gg3Wrl1rd78+39Xq1avRoUMHuLu7IzQ0FBs3blT8eV2eIFJAbGys+PTTT8WhQ4dETk6O6N+/vwgKChJlZWVSmeeff14EBgaKjIwMsXfvXtGjRw/Rs2dP6f6lS5dE586dRUxMjNi/f7/YuHGjaNWqlZgyZYpU5scffxSenp4iKSlJHDlyRHz44YdCp9OJtLS0en3eurJ7927Rtm1bcc8994jx48dL1/nubqy4uFi0adNGPPPMM2LXrl3ixx9/FJs2bRInT56Uyrz99tvC29tbrF27Vhw4cEA88sgjIjg4WFy4cEEq8/DDD4uwsDCxc+dO8f3334u77rpLDBs2TLpfWloqjEajiI+PF4cOHRLLly8XHh4e4p///Ge9Pq+S3nzzTdGyZUuxfv16kZeXJ1avXi28vLzE/PnzpTJ8d3/auHGjePXVV8WaNWsEAPHll1/a3a+vd7V9+3ah0+nE7NmzxZEjR8TUqVNF06ZNRW5ubp2/A1fCBIfqRFFRkQAgtm7dKoQQoqSkRDRt2lSsXr1aKnP06FEBQGRlZQkhLv/lodVqRUFBgVRmwYIFwmAwiIqKCiGEEJMmTRKdOnWy+60hQ4aI2NjYun6kOnfu3DnRrl07kZ6eLh544AEpweG7q9krr7wi7rvvvhvet9lswmQyiXfffVe6VlJSItzc3MTy5cuFEEIcOXJEABB79uyRynzzzTdCo9GIX3/9VQghxD/+8Q/RokUL6X1W/3b79u2VfqR6M2DAADFy5Ei7a4899piIj48XQvDd1eTqBKc+39UTTzwhBgwYYBdPVFSUeO655xR9RlfHISqqE6WlpQAAX19fAEB2djaqqqoQExMjlenQoQOCgoKQlZUFAMjKykJoaCiMRqNUJjY2FhaLBYcPH5bKXNlGdZnqNlyZ2WzGgAEDrnk+vruarVu3DhERERg8eDD8/PzQpUsXLF68WLqfl5eHgoICu2f39vZGVFSU3fvz8fFBRESEVCYmJgZarRa7du2SyvTq1Qt6vV4qExsbi+PHj+Ps2bN1/Zh1omfPnsjIyMAPP/wAADhw4AC2bduGfv36AeC7c0R9viu1/resNCY4pDibzYYJEybg3nvvRefOnQEABQUF0Ov113w112g0oqCgQCpz5T/Q1fer79VUxmKx4MKFC3XxOPVixYoV2LdvH2bNmnXNPb67mv34449YsGAB2rVrh02bNmHs2LF44YUXsGTJEgB/Pv/1nv3Kd+Pn52d3v0mTJvD19XXoHbuayZMnY+jQoejQoQOaNm2KLl26YMKECYiPjwfAd+eI+nxXNyqjlnepFH5NnBRnNptx6NAhbNu2raFDcQn5+fkYP3480tPTG+1Xf+Ww2WyIiIjAW2+9BQDo0qULDh06hIULFyIhIaGBo3Nuq1atwtKlS7Fs2TJ06tQJOTk5mDBhAgICAvjuyOWxB4cUlZiYiPXr1+O7775D69atpesmkwmVlZUoKSmxK19YWAiTySSVuXplUPX5zcoYDAZ4eHgo/Tj1Ijs7G0VFRejatSuaNGmCJk2aYOvWrfjggw/QpEkTGI1Gvrsa+Pv7IyQkxO5ax44dcerUKQB/Pv/1nv3Kd1NUVGR3/9KlSyguLnboHbuaiRMnSr04oaGhGD58OF588UWpJ5Hv7tbV57u6URm1vEulMMEhRQghkJiYiC+//BKbN29GcHCw3f1u3bqhadOmyMjIkK4dP34cp06dQnR0NAAgOjoaubm5dn8BpKenw2AwSP+ARUdH27VRXaa6DVfUt29f5ObmIicnRzoiIiIQHx8v/Znv7sbuvffea7Yk+OGHH9CmTRsAQHBwMEwmk92zWywW7Nq1y+79lZSUIDs7WyqzefNm2Gw2REVFSWUyMzNRVVUllUlPT0f79u3RokWLOnu+unT+/Hlotfb/DOh0OthsNgB8d46oz3el1v+WFdfQs5xJHcaOHSu8vb3Fli1bxOnTp6Xj/PnzUpnnn39eBAUFic2bN4u9e/eK6OhoER0dLd2vXur80EMPiZycHJGWliZuu+226y51njhxojh69KhISUlRxVLnq125ikoIvrua7N69WzRp0kS8+eab4sSJE2Lp0qXC09NT/N///Z9U5u233xY+Pj7iq6++EgcPHhSPPvrodZfvdunSRezatUts27ZNtGvXzm75bklJiTAajWL48OHi0KFDYsWKFcLT09PlljpfKSEhQdx+++3SMvE1a9aIVq1aiUmTJkll+O7+dO7cObF//36xf/9+AUDMnTtX7N+/X/z8889CiPp7V9u3bxdNmjQRc+bMEUePHhXTp0/nMvHrYIJDigBw3ePTTz+Vyly4cEH8/e9/Fy1atBCenp5i0KBB4vTp03bt/PTTT6Jfv37Cw8NDtGrVSrz00kuiqqrKrsx3330nwsPDhV6vF3fccYfdb6jF1QkO313Nvv76a9G5c2fh5uYmOnToIBYtWmR332aziWnTpgmj0Sjc3NxE3759xfHjx+3K/P7772LYsGHCy8tLGAwGMWLECHHu3Dm7MgcOHBD33XefcHNzE7fffrt4++236/zZ6pLFYhHjx48XQUFBwt3dXdxxxx3i1VdftVuizHf3p+++++66f88lJCQIIer3Xa1atUrcfffdQq/Xi06dOokNGzbU2XO7Ko0QV2xZSURERKQCnINDREREqsMEh4iIiFSHCQ4RERGpDhMcIiIiUh0mOERERKQ6THCIiIhIdZjgEBERkeowwSFSiS1btkCj0VzzzarGju+FqHHiRn9EKlFZWYni4mIYjUZoNJqGDqdB9O7dG+Hh4Zg3b550je+FqHFiDw6RSuj1ephMJlX+I37lhwcdpeb3QkQ3xgSHyAn17t0b48aNw4QJE9CiRQsYjUYsXrwY5eXlGDFiBJo3b4677roL33zzjVTn6qGY1NRU+Pj4YNOmTejYsSO8vLzw8MMP4/Tp0zf83bNnzyI+Ph633XYbPDw80K5dO3z66afS/fz8fDzxxBPw8fGBr68vHn30Ufz000/S/WeeeQZxcXGYOXMmbrvtNhgMBjz//POorKyUyqSlpeG+++6Dj48PWrZsib/+9a/473//K93/6aefoNFosHLlSjzwwANwd3fH0qVL8fvvv2PYsGG4/fbb4enpidDQUCxfvtzut7du3Yr58+dDo9FAo9Hgp59+uu4Q1RdffIFOnTrBzc0Nbdu2xXvvvWf3Htq2bYu33noLI0eORPPmzREUFIRFixbd8v9+RNTwmOAQOaklS5agVatW2L17N8aNG4exY8di8ODB6NmzJ/bt24eHHnoIw4cPx/nz52/Yxvnz5zFnzhx8/vnnyMzMxKlTp/Dyyy/fsPy0adNw5MgRfPPNNzh69CgWLFiAVq1aAbjcixIbG4vmzZvj+++/x/bt26Wk6coEJiMjA0ePHsWWLVuwfPlyrFmzBjNnzpTul5eXIykpCXv37kVGRga0Wi0GDRoEm81mF8vkyZMxfvx4HD16FLGxsbh48SK6deuGDRs24NChQxgzZgyGDx+O3bt3AwDmz5+P6OhojB49GqdPn8bp06cRGBh4zTNmZ2fjiSeewNChQ5Gbm4sZM2Zg2rRpSE1NtSv33nvvISIiAvv378ff//53jB07FsePH7/x/2BE5Fwa9lufRHQ9DzzwgLjvvvuk80uXLolmzZqJ4cOHS9dOnz4tAIisrCwhxJ9fOj579qwQQohPP/1UABAnT56U6qSkpAij0XjD3x04cKAYMWLEde99/vnnon379sJms0nXKioqhIeHh9i0aZMQQoiEhATh6+srysvLpTILFiwQXl5ewmq1Xrfd3377TQAQubm5Qggh8vLyBAAxb968G8ZZbcCAAeKll16Szq/+CrsQ176XJ598UvzlL3+xKzNx4kQREhIinbdp00Y89dRT0rnNZhN+fn5iwYIFN42JiJwDe3CInNQ999wj/Vmn06Fly5YIDQ2VrhmNRgBAUVHRDdvw9PTEnXfeKZ37+/vXWH7s2LFYsWIFwsPDMWnSJOzYsUO6d+DAAZw8eRLNmzeHl5cXvLy84Ovri4sXL9oNMYWFhcHT01M6j46ORllZGfLz8wEAJ06cwLBhw3DHHXfAYDCgbdu2AIBTp07ZxRIREWF3brVa8frrryM0NBS+vr7w8vLCpk2brql3M0ePHsW9995rd+3ee+/FiRMnYLVapWtXvn+NRgOTyVTjuyMi59KkoQMgoutr2rSp3blGo7G7Vj1p9uqhnZu1IWpYONmvXz/8/PPP2LhxI9LT09G3b1+YzWbMmTMHZWVl6NatG5YuXXpNvdtuu+2WngkABg4ciDZt2mDx4sUICAiAzWZD586d7Ya5AKBZs2Z25++++y7mz5+PefPmITQ0FM2aNcOECROuqaeU6727mt41ETkXJjhEZOe2225DQkICEhIScP/992PixImYM2cOunbtipUrV8LPzw8Gg+GG9Q8cOIALFy7Aw8MDALBz5054eXkhMDAQv//+O44fP47Fixfj/vvvBwBs27btluLavn07Hn30UTz11FMALid2P/zwA0JCQqQyer3erhfmejp27Ijt27df0/bdd98NnU53S7EQkfPjEBURSZKTk/HVV1/h5MmTOHz4MNavX4+OHTsCAOLj49GqVSs8+uij+P7775GXl4ctW7bghRdewC+//CK1UVlZiVGjRuHIkSPYuHEjpk+fjsTERGi1WrRo0QItW7bEokWLcPLkSWzevBlJSUm3FFu7du2Qnp6OHTt24OjRo3juuedQWFhoV6Zt27bYtWsXfvrpJ5w5c+a6PS4vvfQSMjIy8Prrr+OHH37AkiVL8NFHH9U4+ZqIXA8THCKS6PV6TJkyBffccw969eoFnU6HFStWALg8nyczMxNBQUF47LHH0LFjR4waNQoXL16069Hp27cv2rVrh169emHIkCF45JFHMGPGDACAVqvFihUrkJ2djc6dO+PFF1/Eu+++e0uxTZ06FV27dkVsbCx69+4Nk8mEuLg4uzIvv/wydDodQkJCcNttt113fk7Xrl2xatUqrFixAp07d0ZycjJee+01PPPMM7V6Z0TknLiTMREp5plnnkFJSQnWrl3b0KEQUSPHHhwiIiJSHSY4REREpDocoiIiIiLVYQ8OERERqQ4THCIiIlIdJjhERESkOkxwiIiISHWY4BAREZHqMMEhIiIi1WGCQ0RERKrDBIeIiIhUhwkOERERqc7/BzlShvX7TSFVAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "big_beautiful_frame = pd.read_csv(\"bbf.csv\")\n", + "\n", + "# sep_bins = [1, 10, 100, 1_000, 10_000, 50_000, 155_000]\n", + "# sep_bins = [1, 10, 100, 1_000, 10_000, 20_000, 30_000, 40_000, 50_000, 75_000]\n", + "sep_bins = [1, 10, 100, 500, 1_000, 2_000, 3_000, 4_000, 5_000, 7_500, 10_000]\n", + "\n", + "H, xedges, yedges = np.histogram2d(big_beautiful_frame[\"min_sep\"], \n", + " big_beautiful_frame[\"max_sep\"],\n", + " bins=(sep_bins, sep_bins))\n", + "H = H.T\n", + "\n", + "fig, ax = plt.subplots()\n", + "X, Y = np.meshgrid(xedges, yedges)\n", + "pc = ax.pcolormesh(X, Y, H, shading='flat', cmap=\"rainbow_r\", norm=\"log\")\n", + "ax.set_title(\"\")\n", + "ax.set_xlabel('min separation')\n", + "ax.set_ylabel('max separation')\n", + "# ax.set_yscale('log')\n", + "# ax.set_yscale('log')\n", + "fig.colorbar(pc)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "154558.46597361984" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "big_beautiful_frame[\"max_sep\"].max()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([5, 6, 3, 4, 2, 2, 3, 2, 3, 3])" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sep_bins = [1, 10, 100, 500, 1_000, 2_000, 3_000, 4_000, 5_000, 7_500, 10_000, 20_000, 30_000, 40_000, 50_000,100_000, 155_000]\n", + "# len(sep_bins)\n", + "\n", + "right = np.searchsorted(sep_bins, big_beautiful_frame[\"min_sep\"])\n", + "left = np.searchsorted(sep_bins, big_beautiful_frame[\"max_sep\"])\n", + "span = left - right + 1\n", + "span[0:10]" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "bin_hist, bin_bins = np.histogram(span, bins = np.arange(1, 19))" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABLEAAAEpCAYAAABssrDdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxlklEQVR4nO3deVRV9f7/8dcB44DKcAERURTUnIcM08ypREXyKpblsLzhkA1+ccqstFKzm6GVaaWSeRMbnC01M0WvcyUOmTezMk1RKwWHBMRE5ezfHy7OzyPjQYatPh9rnbXcm8/+7PfefNgHX+z9ORbDMAwBAAAAAAAAJuZS1gUAAAAAAAAABSHEAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeoRYAAAAAAAAMD1CLAAAAAAAAJgeIRYAAAAAAABMjxALAAAAAAAApkeIBQAAAAAAANMjxAIAAA42b94si8WiZcuWFdh2wIABCgkJKfmiCiE5OVmPPPKI/Pz8ZLFYNH369LIuySnz5s2TxWJRUlJSie0j+3u7efPmEtvH9V555RVZLBadPn26wLYhISEaMGBAyRdlYtnnq6SUxRgAAKC4lCvrAgAAAIrDM888o4SEBE2YMEGBgYFq3rx5WZeEW9isWbNUvnz52z50AwCgNBFiAQCAIpszZ45sNltZlyFJ2rhxo6KiojR69OiyLqVIHnvsMfXp00dWq7XE9tGuXTv9/fffcnNzK7F93IgDBw7IxeXmeFBg1qxZ8vf3L/YQ6+WXX9aYMWOKtU8AAG4VhFgAAKDI7rjjjrIuwS4lJUU+Pj7F1t/Fixfl5uZWaqGKq6urXF1d821jGIYuXrwoDw+PIu3DxcVF7u7uRdq2NJRkgHezKFeunMqV41d0AAByc3P8qQsAgCLKnl/m0KFDGjBggHx8fOTt7a2BAwfqwoUL9nZJSUmyWCyaN29ejj4sFoteeeWVHH3++uuv+te//iVvb29VqlRJ48aNk2EYOn78uKKiouTl5aXAwEBNnTq1ULWuX79ebdq0kY+PjypWrKi6devqxRdfdGiTmZmpCRMmqHbt2rJarQoODtbzzz+vzMzMHDUPHTpU8+fPV926deXu7q6wsDBt3bq10OcuKytLL774ogIDA1WhQgV1795dx48fd2hz/ZxY2efxrbfe0gcffKBatWrJarXqnnvu0a5duxy2PXnypAYOHKhq1arJarWqSpUqioqKcnpOqOy5pAzD0MyZM2WxWBzmFDp8+LAeffRR+fr6qnz58rr33nu1evVqhz6y5wlatGiRXn75ZVWtWlXly5dXWlparvu89jinTZumGjVqyMPDQ+3bt9ePP/7o0PaHH37QgAEDVLNmTbm7uyswMFCDBg3SmTNncj2Oa48/JCRE//znP5WQkKDmzZvLw8NDs2fP1sMPP6y7777bYftu3brJYrHoiy++sK/bsWOHLBaL1qxZ43Cc186HdPDgQfXs2VOBgYFyd3dXtWrV1KdPH6Wmpjr0/+mnnyosLEweHh7y9fVVnz59coyH/Jw7dy7fn8Hs4732zqbsc7J161Y99dRT8vPzk5eXl6Kjo/XXX385bLt7925FRETI399fHh4eCg0N1aBBgxza2Gw2vfPOO2rcuLHc3d1VqVIldenSRbt377a3iY+PV4cOHRQQECCr1aoGDRooLi4uR5379+/Xli1b7OPt/vvvz/PYnRkv18+JFR8fL4vForlz5zq0e/3112WxWPTVV1/Z1/3yyy965JFH5OvrK3d3dzVv3txhPOSlsGMAAICyxp95AAC3hV69eik0NFSxsbHas2eP/vOf/yggIEBTpkwpcp+9e/dW/fr1NXnyZK1evVqvvfaafH19NXv2bHXo0EFTpkzR/PnzNXr0aN1zzz1q165dnn3t379f//znP9WkSRO9+uqrslqtOnTokL755ht7G5vNpu7du+vrr7/Wk08+qfr162vfvn2aNm2afv31V61YscKhzy1btmjx4sUaPny4rFarZs2apS5dumjnzp1q1KhRgcc3adIkWSwWvfDCC0pJSdH06dPVsWNH7d27t8A7gRYsWKD09HQ99dRTslgseuONN/Twww/r8OHD9ru3evbsqf3792vYsGEKCQlRSkqK1q9fr2PHjjk1WXy7du30ySef6LHHHlOnTp0UHR1t/1pycrLuu+8+XbhwQcOHD5efn58++ugjde/eXcuWLdNDDz3k0Ne///1vubm5afTo0crMzCzwsbuPP/5Y6enpiomJ0cWLF/XOO++oQ4cO2rdvnypXrizpajh5+PBhDRw4UIGBgdq/f78++OAD7d+/X4mJiQVO4n3gwAH17dtXTz31lJ544gnVrVtXhmFo5cqVSktLk5eXlwzD0DfffCMXFxdt27ZN3bt3lyRt27ZNLi4uat26da59X7p0SREREcrMzNSwYcMUGBioP/74Q19++aXOnTsnb29vSVfHwrhx49SrVy8NHjxYp06d0nvvvad27drp+++/L9QdcDfyMzh06FD5+PjolVde0YEDBxQXF6ejR4/aQ7mUlBR17txZlSpV0pgxY+Tj46OkpCR9/vnnDv08/vjjmjdvniIjIzV48GBduXJF27ZtU2Jion0Otbi4ODVs2FDdu3dXuXLltGrVKv3f//2fbDabYmJiJEnTp0/XsGHDVLFiRb300kuSZP9+56cw4+V6AwcO1Oeff65Ro0apU6dOCg4O1r59+zRx4kQ9/vjjevDBByVdvYa0bt1aVatW1ZgxY1ShQgUtWbJEPXr00GeffZZjrGcr7BgAAMAUDAAAbmETJkwwJBmDBg1yWP/QQw8Zfn5+9uUjR44Ykoz4+PgcfUgyJkyYkKPPJ5980r7uypUrRrVq1QyLxWJMnjzZvv6vv/4yPDw8jP79++db57Rp0wxJxqlTp/Js88knnxguLi7Gtm3bHNa///77hiTjm2++cahZkrF79277uqNHjxru7u7GQw89lG8tmzZtMiQZVatWNdLS0uzrlyxZYkgy3nnnHfu6/v37GzVq1LAvZ59HPz8/4+zZs/b1K1euNCQZq1atMgzj6nmRZLz55pv51uIMSUZMTIzDupEjRxqSHM5Zenq6ERoaaoSEhBhZWVkOx1yzZk3jwoULBe4r+zg9PDyM33//3b5+x44dhiTjmWeesa/Lrb+FCxcakoytW7fa18XHxxuSjCNHjtjX1ahRw5BkrF271mH7Xbt2GZKMr776yjAMw/jhhx8MScajjz5qtGzZ0t6ue/fuRrNmzezL2ce5adMmwzAM4/vvvzckGUuXLs3zWJOSkgxXV1dj0qRJDuv37dtnlCtXLsf66xX2ZzD7eK/9Wck+J2FhYcalS5fs69944w1DkrFy5UrDMAxj+fLlhiRj165dedaxceNGQ5IxfPjwHF+z2Wz2f+f2/YqIiDBq1qzpsK5hw4ZG+/bt89zftZwZL9nn61onTpwwfH19jU6dOhmZmZlGs2bNjOrVqxupqan2NuHh4Ubjxo2NixcvOhzXfffdZ9x55532dUUZAwAAmAWPEwIAbgtPP/20w3Lbtm115syZPB8XK4zBgwfb/+3q6qrmzZvLMAw9/vjj9vU+Pj6qW7euDh8+nG9f2XeyrFy5Ms+J0pcuXar69eurXr16On36tP3VoUMHSdKmTZsc2rdq1UphYWH25erVqysqKkoJCQnKysoq8Piio6Pl6elpX37kkUdUpUoVh8eX8tK7d2/94x//sC+3bdtWkuznwcPDQ25ubtq8eXOOx8KK01dffaUWLVqoTZs29nUVK1bUk08+qaSkJP30008O7fv37+/UfFM9evRQ1apV7cstWrRQy5YtHc7Rtf1dvHhRp0+f1r333itJ2rNnT4H7CA0NVUREhMO6Zs2aqWLFivbHQ7dt26Zq1aopOjpae/bs0YULF2QYhr7++mv7uc9N9l02CQkJOR7ty/b555/LZrOpV69eDuMuMDBQd955Z45xl5cb+Rl88sknHeZfGzJkiMqVK2c/z9k/P19++aUuX76cax+fffaZLBaLJkyYkONr194Nd+33KzU1VadPn1b79u11+PDhG368rjDjJTeBgYGaOXOm1q9fr7Zt22rv3r2aO3euvLy8JElnz57Vxo0b1atXL6Wnp9u/R2fOnFFERIQOHjyoP/74I9e+CzMGAAAwizILsbZu3apu3bopKChIFoslxyMQhWEYht566y3VqVNHVqtVVatW1aRJk4q/WADATa969eoOy9kBy40EKNf36e3tLXd3d/n7++dYX9B+evfurdatW2vw4MGqXLmy+vTpoyVLljgEWgcPHtT+/ftVqVIlh1edOnUkXZ3Y/Fp33nlnjv3UqVNHFy5c0KlTpwo8vuu3t1gsql27dqHmrCrofFutVk2ZMkVr1qxR5cqV1a5dO73xxhs6efJkgX074+jRo6pbt26O9fXr17d//VqhoaFO9Z/XOb72HJ09e1YjRoxQ5cqV5eHhoUqVKtn3U5hQJLeaXF1d1apVK23btk3S1RCrbdu2atOmjbKyspSYmKiffvpJZ8+ezTfECg0N1ahRo/Sf//xH/v7+ioiI0MyZMx3qOnjwoAzD0J133plj7P388885xl1ebuRn8PrzXLFiRVWpUsV+ntu3b6+ePXtq4sSJ8vf3V1RUlOLj4x3mivvtt98UFBQkX1/ffPf1zTffqGPHjqpQoYJ8fHxUqVIl+9x0NxpiFWa85KVPnz7q2rWrdu7cqSeeeELh4eH2rx06dEiGYWjcuHE5vkfZoV1e36fCjAEAAMyizObEysjIUNOmTTVo0CA9/PDDRepjxIgRWrdund566y01btxYZ8+e1dmzZ4u5UgDArSCvT30zDEOS8pyXKL87lnLrs6D95MXDw0Nbt27Vpk2btHr1aq1du1aLFy9Whw4dtG7dOrm6uspms6lx48Z6++23c+0jODg4332UpsKch5EjR6pbt25asWKFEhISNG7cOMXGxmrjxo1q1qxZaZXqoKif+pefXr166dtvv9Vzzz2nu+66SxUrVpTNZlOXLl3yvOuuMDW1adNGkyZN0sWLF7Vt2za99NJL8vHxUaNGjbRt2zb7HEv5hViSNHXqVA0YMEArV67UunXrNHz4cMXGxioxMVHVqlWTzWazTw6f2/e1YsWKhTgLRf/ZKAyLxaJly5YpMTFRq1atUkJCggYNGqSpU6cqMTGx0DX+9ttvCg8PV7169fT2228rODhYbm5u+uqrrzRt2rRCfb9KypkzZ+wT0P/000+y2Wz2T87Mrmv06NE57trLVrt27Tz7LmgMAABgFmUWYkVGRioyMjLPr2dmZuqll17SwoULde7cOTVq1EhTpkyxf/LLzz//rLi4OP3444/2v7A6+9dTAACyZd8Vcu7cOYf119+pU5JcXFwUHh6u8PBwvf3223r99df10ksvadOmTerYsaNq1aql//3vfwoPDy9wMnDp6h001/v1119Vvnx5VapUyentDcPQoUOH1KRJk8IfVAFq1aqlZ599Vs8++6wOHjyou+66S1OnTtWnn35aLP3XqFFDBw4cyLH+l19+sX/9RuR1jrMnpv/rr7+0YcMGTZw4UePHj893O2e1bdtWly5d0sKFC/XHH3/Yw6p27drZQ6w6deoUasLxxo0bq3Hjxnr55Zf17bffqnXr1nr//ff12muvqVatWjIMQ6Ghofa7/krbwYMH9cADD9iXz58/rxMnTtgnNc9277336t5779WkSZO0YMEC9evXT4sWLdLgwYNVq1YtJSQk6OzZs3nejbVq1SplZmbqiy++cLhzLLdHJgvzM5jbcVzv2vGSn5iYGKWnpys2NlZjx47V9OnTNWrUKElSzZo1JUl33HGHOnbs6HRdUv5jAAAAszDtnFhDhw7V9u3btWjRIv3www969NFH1aVLF/ub/6pVq1SzZk19+eWXCg0NVUhIiAYPHsydWACAIvHy8pK/v799jqFss2bNKpX95/b+ddddd0mS/ZGoXr166Y8//tCcOXNytP3777+VkZHhsG779u0Ocy4dP35cK1euVOfOnfO8K+Za2Z+klm3ZsmU6ceJEvn+EKqwLFy7o4sWLDutq1aolT09Ph0fAbtSDDz6onTt3avv27fZ1GRkZ+uCDDxQSEqIGDRrcUP8rVqxwmGto586d2rFjh/0cZZ/n6+82mj59+g3tV5JatmypO+64Q1OmTJGvr68aNmwo6Wq4lZiYqC1bthR4F1ZaWpquXLnisK5x48ZycXGxfx8efvhhubq6auLEiTmOwzAMnTlz5oaPpSAffPCBw1xXcXFxunLliv08//XXXzlqu/7np2fPnjIMQxMnTszRf/a2uX2/UlNTFR8fn2ObChUq5Ai9C1LQeMnLsmXLtHjxYk2ePFljxoxRnz599PLLL+vXX3+VJAUEBOj+++/X7NmzdeLEiRzb5/f4cGHGAAAAZlFmd2Ll59ixY4qPj9exY8cUFBQk6ert0WvXrlV8fLxef/11HT58WEePHtXSpUv18ccfKysrS88884weeeQRbdy4sYyPAABwMxo8eLAmT56swYMHq3nz5tq6dav9P4kl7dVXX9XWrVvVtWtX1ahRQykpKZo1a5aqVatmn5T8scce05IlS/T0009r06ZNat26tbKysvTLL79oyZIlSkhIUPPmze19NmrUSBERERo+fLisVqs9kMvtP/G58fX1VZs2bTRw4EAlJydr+vTpql27tp544okbPt5ff/1V4eHh6tWrlxo0aKBy5cpp+fLlSk5OVp8+fezt5s2bp4EDByo+Pl4DBgxwej9jxozRwoULFRkZqeHDh8vX11cfffSRjhw5os8++8z+OFZR1a5dW23atNGQIUOUmZmp6dOny8/PT88//7ykq+Fo9nxfly9fVtWqVbVu3TodOXLkhvYrSeXLl1dYWJgSExPVrVs3+51B7dq1U0ZGhjIyMgoMsTZu3KihQ4fq0UcfVZ06dXTlyhV98skncnV1Vc+ePSVdDRdfe+01jR07VklJSerRo4c8PT115MgRLV++XE8++aRGjx59w8eTn0uXLtnHy4EDBzRr1iy1adNG3bt3lyR99NFHmjVrlh566CHVqlVL6enpmjNnjry8vOx3az3wwAN67LHH9O677+rgwYP2xzm3bdumBx54QEOHDlXnzp3l5uambt266amnntL58+c1Z84cBQQE5AiHwsLCFBcXp9dee021a9dWQECA/UMW8lLQeMlNSkqKhgwZYq9RkmbMmKFNmzZpwIAB+vrrr+Xi4qKZM2eqTZs2aty4sZ544gnVrFlTycnJ2r59u37//Xf973//y7X/wowBAADMwpQh1r59+5SVlZXjlvXMzEz5+flJuvrsf2Zmpj7++GN7uw8//FBhYWE6cOBArpO4AgCQn/Hjx+vUqVNatmyZlixZosjISK1Zs0YBAQElvu/u3bsrKSlJc+fO1enTp+Xv76/27dtr4sSJ9k8Pc3Fx0YoVKzRt2jR9/PHHWr58ucqXL6+aNWtqxIgROd4327dvr1atWmnixIk6duyYGjRooHnz5hX6ccAXX3xRP/zwg2JjY5Wenq7w8HDNmjVL5cuXv+HjDQ4OVt++fbVhwwZ98sknKleunOrVq6clS5Y4/Mf5/PnzkqQqVaoUaT+VK1fWt99+qxdeeEHvvfeeLl68qCZNmmjVqlXq2rXrDR9HdHS0XFxcNH36dKWkpKhFixaaMWOGQ70LFizQsGHDNHPmTBmGoc6dO2vNmjX2P9TdiOy7rq799MXAwEDVrl1bhw4dKjDEatq0qSIiIrRq1Sr98ccfKl++vJo2bao1a9bYP0FRuhoG1qlTR9OmTbOHoMHBwercubM9SCpJM2bM0Pz58zV+/HhdvnxZffv21bvvvmsP7tq3b6+dO3dq0aJFSk5Olre3t1q0aKH58+c7TDcRHx+vJk2a6MMPP9Rzzz0nb29vNW/eXPfdd58kqW7dulq2bJlefvlljR49WoGBgRoyZIgqVaqkQYMGOdQ0fvx4HT16VG+88YbS09PVvn37AkOswoyX62UHXvHx8fbj9fPz0wcffKCoqCi99dZbev7559WgQQPt3r1bEydO1Lx583TmzBkFBASoWbNmDo+yXq+wYwAAADOwGMUxm+aNFmGxaPny5erRo4ckafHixerXr5/279+f43GHihUrKjAwUBMmTNDrr7/ucGv533//rfLly2vdunXq1KlTaR4CAACmYrFYFBMToxkzZpR1KTekV69eSkpK0s6dO8u6FAdJSUkKDQ3Vm2++WeJ3Id3Osu/E27Vrl8NdhjcbxgsAAMXDlHdiNWvWTFlZWUpJScnzL4itW7fWlStX9Ntvv6lWrVqSZH/k40YnaQUAAGXPMAxt3ry52CZ5BwAAwM2tzEKs8+fP69ChQ/blI0eOaO/evfL19VWdOnXUr18/RUdHa+rUqWrWrJlOnTqlDRs2qEmTJuratas6duyou+++W4MGDdL06dNls9kUExOjTp06ldkn5wAAgOJjsViUkpJS1mUAAADAJMrs0wl3796tZs2aqVmzZpKkUaNGOTyzHx8fr+joaD377LOqW7euevTooV27dtk/7tjFxUWrVq2Sv7+/2rVrp65du6p+/fpatGhRWR0SAAAAAAAASogp5sQCAAAAAAAA8lNmd2IBAAAAAAAAhUWIBQAAAAAAANMr9YndbTab/vzzT3l6espisZT27gEAAAAAAGAShmEoPT1dQUFBcnHJ/16rUg+x/vzzTwUHB5f2bgEAAAAAAGBSx48fV7Vq1fJtU+ohlqenp6SrxXl5eZX27gEAAAAAAGASaWlpCg4OtudF+Sn1ECv7EUIvLy9CLAAAAAAAABRqyikmdgcAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTK1fWBdwKQsasLusSykzS5K5lXQIAAAAAALgNcCcWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPQIsQAAAAAAAGB6hFgAAAAAAAAwPUIsAAAAAAAAmB4hFgAAAAAAAEyPEAsAAAAAAACmR4gFAAAAAAAA0yPEAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeoRYAAAAAAAAMD1CLAAAAAAAAJgeIRYAAAAAAABMjxALAAAAAAAApkeIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTu6EQa/LkybJYLBo5cmQxlQMAAAAAAADkVOQQa9euXZo9e7aaNGlSnPUAAAAAAAAAORQpxDp//rz69eunOXPm6B//+Edx1wQAAAAAAAA4KFKIFRMTo65du6pjx44Fts3MzFRaWprDCwAAAAAAAHBGOWc3WLRokfbs2aNdu3YVqn1sbKwmTpzodGEAAAAAAABANqfuxDp+/LhGjBih+fPny93dvVDbjB07VqmpqfbX8ePHi1QoAAAAAAAAbl9O3Yn13XffKSUlRXfffbd9XVZWlrZu3aoZM2YoMzNTrq6uDttYrVZZrdbiqRYAAAAAAAC3JadCrPDwcO3bt89h3cCBA1WvXj298MILOQIsAAAAAAAAoDg4FWJ5enqqUaNGDusqVKggPz+/HOsBAAAAAACA4lKkTycEAAAAAAAASpPTn054vc2bNxdDGQAAAAAAAEDeuBMLAAAAAAAApkeIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPQIsQAAAAAAAGB6hFgAAAAAAAAwPUIsAAAAAAAAmB4hFgAAAAAAAEyPEAsAAAAAAACmR4gFAAAAAAAA0yPEAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeoRYAAAAAAAAMD1CLAAAAAAAAJgeIRYAAAAAAABMjxALAAAAAAAApkeIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPScCrHi4uLUpEkTeXl5ycvLS61atdKaNWtKqjYAAAAAAABAkpMhVrVq1TR58mR999132r17tzp06KCoqCjt37+/pOoDAAAAAAAAVM6Zxt26dXNYnjRpkuLi4pSYmKiGDRsWa2EAAAAAAABANqdCrGtlZWVp6dKlysjIUKtWrfJsl5mZqczMTPtyWlpaUXcJAAAAAACA25TTE7vv27dPFStWlNVq1dNPP63ly5erQYMGebaPjY2Vt7e3/RUcHHxDBQMAAAAAAOD243SIVbduXe3du1c7duzQkCFD1L9/f/300095th87dqxSU1Ptr+PHj99QwQAAAAAAALj9OP04oZubm2rXri1JCgsL065du/TOO+9o9uzZuba3Wq2yWq03ViUAAAAAAABua07fiXU9m83mMOcVAAAAAAAAUNycuhNr7NixioyMVPXq1ZWenq4FCxZo8+bNSkhIKKn6AAAAAAAAAOdCrJSUFEVHR+vEiRPy9vZWkyZNlJCQoE6dOpVUfQAAAAAAAIBzIdaHH35YUnUAAAAAAAAAebrhObEAAAAAAACAkkaIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPTKlXUBuLmFjFld1iWUiaTJXcu6BAAAAAAAbivciQUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPQIsQAAAAAAAGB6hFgAAAAAAAAwPUIsAAAAAAAAmB4hFgAAAAAAAEyPEAsAAAAAAACmR4gFAAAAAAAA0yPEAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeoRYAAAAAAAAMD1CLAAAAAAAAJgeIRYAAAAAAABMjxALAAAAAAAApkeIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATM+pECs2Nlb33HOPPD09FRAQoB49eujAgQMlVRsAAAAAAAAgyckQa8uWLYqJiVFiYqLWr1+vy5cvq3PnzsrIyCip+gAAAAAAAACVc6bx2rVrHZbnzZungIAAfffdd2rXrl2xFgYAAAAAAABkcyrEul5qaqokydfXN882mZmZyszMtC+npaXdyC4BAAAAAABwGyryxO42m00jR45U69at1ahRozzbxcbGytvb2/4KDg4u6i4BAAAAAABwmypyiBUTE6Mff/xRixYtyrfd2LFjlZqaan8dP368qLsEAAAAAADAbapIjxMOHTpUX375pbZu3apq1arl29ZqtcpqtRapOAAAAAAAAEByMsQyDEPDhg3T8uXLtXnzZoWGhpZUXQAAAAAAAICdUyFWTEyMFixYoJUrV8rT01MnT56UJHl7e8vDw6NECgQAAAAAAACcmhMrLi5Oqampuv/++1WlShX7a/HixSVVHwAAAAAAAOD844QAAAAAAABAaSvypxMCAAAAAAAApYUQCwAAAAAAAKZHiAUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPQIsQAAAAAAAGB6hFgAAAAAAAAwPUIsAAAAAAAAmB4hFgAAAAAAAEyPEAsAAAAAAACmR4gFAAAAAAAA0yPEAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeoRYAAAAAAAAMD1CLAAAAAAAAJgeIRYAAAAAAABMjxALAAAAAAAApkeIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPQIsQAAAAAAAGB6hFgAAAAAAAAwPUIsAAAAAAAAmB4hFgAAAAAAAEyPEAsAAAAAAACmR4gFAAAAAAAA0yPEAgAAAAAAgOk5HWJt3bpV3bp1U1BQkCwWi1asWFECZQEAAAAAAAD/n9MhVkZGhpo2baqZM2eWRD0AAAAAAABADuWc3SAyMlKRkZElUQsAAAAAAACQK6dDLGdlZmYqMzPTvpyWllbSuwQAAAAAAMAtpsQndo+NjZW3t7f9FRwcXNK7BAAAAAAAwC2mxEOssWPHKjU11f46fvx4Se8SAAAAAAAAt5gSf5zQarXKarWW9G4AAAAAAABwCyvxO7EAAAAAAACAG+X0nVjnz5/XoUOH7MtHjhzR3r175evrq+rVqxdrcQAAAAAAAIBUhBBr9+7deuCBB+zLo0aNkiT1799f8+bNK7bCAAAAAAAAgGxOh1j333+/DMMoiVoAAAAAAACAXDEnFgAAAAAAAEyPEAsAAAAAAACmR4gFAAAAAAAA0yPEAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeoRYAAAAAAAAMD1CLAAAAAAAAJgeIRYAAAAAAABMjxALAAAAAAAApkeIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADTI8QCAAAAAACA6RFiAQAAAAAAwPQIsQAAAAAAAGB6hFgAAAAAAAAwPUIsAAAAAAAAmB4hFgAAAAAAAEyPEAsAAAAAAACmR4gFAAAAAAAA0yPEAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeuXKugDgZhQyZnVZl1AmkiZ3LesSAAAAAAC3Ke7EAgAAAAAAgOkRYgEAAAAAAMD0CLEAAAAAAABgeoRYAAAAAAAAMD1CLAAAAAAAAJhekUKsmTNnKiQkRO7u7mrZsqV27txZ3HUBAAAAAAAAdk6HWIsXL9aoUaM0YcIE7dmzR02bNlVERIRSUlJKoj4AAAAAAABAFsMwDGc2aNmype655x7NmDFDkmSz2RQcHKxhw4ZpzJgxBW6flpYmb29vpaamysvLq2hVm0zImNVlXQIAFLukyV3LugQAAAAAtzhncqJyznR86dIlfffddxo7dqx9nYuLizp27Kjt27fnuk1mZqYyMzPty6mpqfYibxW2zAtlXQIAFLvqzywt6xJQyn6cGFHWJQAAAOA2k50PFeYeK6dCrNOnTysrK0uVK1d2WF+5cmX98ssvuW4TGxuriRMn5lgfHBzszK4BAEAJ855e1hUAAADgdpWeni5vb+982zgVYhXF2LFjNWrUKPuyzWbT2bNn5efnJ4vFUtK7xy0qLS1NwcHBOn78+C3zWCrMibGG0sR4Q2lhrKG0MNZQmhhvKC2MteJlGIbS09MVFBRUYFunQix/f3+5uroqOTnZYX1ycrICAwNz3cZqtcpqtTqs8/HxcWa3QJ68vLy4aKBUMNZQmhhvKC2MNZQWxhpKE+MNpYWxVnwKugMrm1OfTujm5qawsDBt2LDBvs5ms2nDhg1q1aqVcxUCAAAAAAAAheT044SjRo1S//791bx5c7Vo0ULTp09XRkaGBg4cWBL1AQAAAAAAAM6HWL1799apU6c0fvx4nTx5UnfddZfWrl2bY7J3oCRZrVZNmDAhx6OqQHFjrKE0Md5QWhhrKC2MNZQmxhtKC2Ot7FiMwnyGIQAAAAAAAFCGnJoTCwAAAAAAACgLhFgAAAAAAAAwPUIsAAAAAAAAmB4hFgAAAAAAAEyPEAumExsbq3vuuUeenp4KCAhQjx49dODAgXy3mTdvniwWi8PL3d29lCrGzeqVV17JMW7q1auX7zZLly5VvXr15O7ursaNG+urr74qpWpxswsJCckx3iwWi2JiYnJtz3UNhbV161Z169ZNQUFBslgsWrFihcPXDcPQ+PHjVaVKFXl4eKhjx446ePBggf3OnDlTISEhcnd3V8uWLbVz584SOgLcTPIbb5cvX9YLL7ygxo0bq0KFCgoKClJ0dLT+/PPPfPssyvsxbn0FXdsGDBiQY9x06dKlwH65tuF6BY213H5/s1gsevPNN/Psk+taySHEguls2bJFMTExSkxM1Pr163X58mV17txZGRkZ+W7n5eWlEydO2F9Hjx4tpYpxM2vYsKHDuPn666/zbPvtt9+qb9++evzxx/X999+rR48e6tGjh3788cdSrBg3q127djmMtfXr10uSHn300Ty34bqGwsjIyFDTpk01c+bMXL/+xhtv6N1339X777+vHTt2qEKFCoqIiNDFixfz7HPx4sUaNWqUJkyYoD179qhp06aKiIhQSkpKSR0GbhL5jbcLFy5oz549GjdunPbs2aPPP/9cBw4cUPfu3Qvs15n3Y9weCrq2SVKXLl0cxs3ChQvz7ZNrG3JT0Fi7doydOHFCc+fOlcViUc+ePfPtl+taCTEAk0tJSTEkGVu2bMmzTXx8vOHt7V16ReGWMGHCBKNp06aFbt+rVy+ja9euDutatmxpPPXUU8VcGW4HI0aMMGrVqmXYbLZcv851DUUhyVi+fLl92WazGYGBgcabb75pX3fu3DnDarUaCxcuzLOfFi1aGDExMfblrKwsIygoyIiNjS2RunFzun685Wbnzp2GJOPo0aN5tnH2/Ri3n9zGWv/+/Y2oqCin+uHahoIU5roWFRVldOjQId82XNdKDndiwfRSU1MlSb6+vvm2O3/+vGrUqKHg4GBFRUVp//79pVEebnIHDx5UUFCQatasqX79+unYsWN5tt2+fbs6duzosC4iIkLbt28v6TJxi7l06ZI+/fRTDRo0SBaLJc92XNdwo44cOaKTJ086XLu8vb3VsmXLPK9dly5d0nfffeewjYuLizp27Mj1Dk5LTU2VxWKRj49Pvu2ceT8Gsm3evFkBAQGqW7euhgwZojNnzuTZlmsbikNycrJWr16txx9/vMC2XNdKBiEWTM1ms2nkyJFq3bq1GjVqlGe7unXrau7cuVq5cqU+/fRT2Ww23Xffffr9999LsVrcbFq2bKl58+Zp7dq1iouL05EjR9S2bVulp6fn2v7kyZOqXLmyw7rKlSvr5MmTpVEubiErVqzQuXPnNGDAgDzbcF1Dcci+Pjlz7Tp9+rSysrK43uGGXbx4US+88IL69u0rLy+vPNs5+34MSFcfJfz444+1YcMGTZkyRVu2bFFkZKSysrJybc+1DcXho48+kqenpx5++OF823FdKznlyroAID8xMTH68ccfC3x+uFWrVmrVqpV9+b777lP9+vU1e/Zs/fvf/y7pMnGTioyMtP+7SZMmatmypWrUqKElS5YU6q8rQFF9+OGHioyMVFBQUJ5tuK4BuJldvnxZvXr1kmEYiouLy7ct78coij59+tj/3bhxYzVp0kS1atXS5s2bFR4eXoaV4VY2d+5c9evXr8AP2+G6VnK4EwumNXToUH355ZfatGmTqlWr5tS2d9xxh5o1a6ZDhw6VUHW4Ffn4+KhOnTp5jpvAwEAlJyc7rEtOTlZgYGBplIdbxNGjR/Xf//5XgwcPdmo7rmsoiuzrkzPXLn9/f7m6unK9Q5FlB1hHjx7V+vXr870LKzcFvR8DualZs6b8/f3zHDdc23Cjtm3bpgMHDjj9O5zEda04EWLBdAzD0NChQ7V8+XJt3LhRoaGhTveRlZWlffv2qUqVKiVQIW5V58+f12+//ZbnuGnVqpU2bNjgsG79+vUOd8sABYmPj1dAQIC6du3q1HZc11AUoaGhCgwMdLh2paWlaceOHXleu9zc3BQWFuawjc1m04YNG7jeoUDZAdbBgwf13//+V35+fk73UdD7MZCb33//XWfOnMlz3HBtw4368MMPFRYWpqZNmzq9Lde14kOIBdOJiYnRp59+qgULFsjT01MnT57UyZMn9ffff9vbREdHa+zYsfblV199VevWrdPhw4e1Z88e/etf/9LRo0eLlJLj9jF69Ght2bJFSUlJ+vbbb/XQQw/J1dVVffv2lZRznI0YMUJr167V1KlT9csvv+iVV17R7t27NXTo0LI6BNxkbDab4uPj1b9/f5Ur5/hEP9c1FNX58+e1d+9e7d27V9LVydz37t2rY8eOyWKxaOTIkXrttdf0xRdfaN++fYqOjlZQUJB69Ohh7yM8PFwzZsywL48aNUpz5szRRx99pJ9//llDhgxRRkaGBg4cWMpHB7PJb7xdvnxZjzzyiHbv3q358+crKyvL/nvcpUuX7H1cP94Kej/G7Sm/sXb+/Hk999xzSkxMVFJSkjZs2KCoqCjVrl1bERER9j64tqEw8htr2dLS0rR06dI8fw/julaKyvrjEYHrScr1FR8fb2/Tvn17o3///vblkSNHGtWrVzfc3NyMypUrGw8++KCxZ8+e0i8eN5XevXsbVapUMdzc3IyqVasavXv3Ng4dOmT/+vXjzDAMY8mSJUadOnUMNzc3o2HDhsbq1atLuWrczBISEgxJxoEDB3J8jesaimrTpk25vm9mjyebzWaMGzfOqFy5smG1Wo3w8PAcY7BGjRrGhAkTHNa999579jHYokULIzExsZSOCGaW33g7cuRInr/Hbdq0yd7H9eOtoPdj3J7yG2sXLlwwOnfubFSqVMm44447jBo1ahhPPPGEcfLkSYc+uLahMAp6HzUMw5g9e7bh4eFhnDt3Ltc+uK6VHothGEaJJ2UAAAAAAADADeBxQgAAAAAAAJgeIRYAAAAAAABMjxALAAAAAAAApkeIBQAAAAAAANMjxAIAAAAAAIDpEWIBAAAAAADA9AixAAAAAAAAYHqEWAAAAAAAADA9QiwAAAAAAACYHiEWAAAAAAAATI8QCwAAAAAAAKZHiAUAAAAAAADT+3/RNj1eu0U5hgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "width = np.diff(bin_bins)\n", + "center = (bin_bins[:-1] + bin_bins[1:]) / 2\n", + "\n", + "fig, ax = plt.subplots(figsize=(15,3))\n", + "ax.bar(center, bin_hist, align='center', width=width)\n", + "# ax.set_yscale('log')\n", + "plt.title(f\"num sep bins, for pairwise hipscat pixels\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([4680685, 2606439, 223729, 43134, 55660, 36807, 31794,\n", + " 12031, 9324, 8515, 4957, 7110, 5777, 7190,\n", + " 2282, 378, 399])" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bin_hist" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "num_pairs = np.sum(bin_hist)\n", + "proportion = bin_hist / num_pairs" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([6.05035850e-01, 3.36914156e-01, 2.89197128e-02, 5.57559767e-03,\n", + " 7.19473654e-03, 4.75775544e-03, 4.10976381e-03, 1.55515407e-03,\n", + " 1.20524117e-03, 1.10066801e-03, 6.40752947e-04, 9.19054560e-04,\n", + " 7.46747988e-04, 9.29395540e-04, 2.94976443e-04, 4.88611285e-05,\n", + " 5.15756357e-05])" + ] + }, + "execution_count": 44, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "proportion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "60% of tile pairs will have only one possible separation bin. As a result, we only need the total masked count for these tile-pairs.\n", + "\n", + "\n", + "33% will have two possible bins. Once pair separations have been determined, the binning clause can be quickly short-circuited. We will know BEFORE loading the parquet tiles into memory what the possible bins are. Can pass just the dividing line between bins to counting routine! Could just return `count sep less than X`.\n", + "\n", + "< 3% will have three possible bins.\n", + "\n", + "Remaining 3% of tile pairs have 4+ possible bins. These are the smallest separations, and the ones that we will spend the most time on (and where we SHOULD be spending the most time getting the counts right!\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.0291302809605374" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "long_tail_tile_pairs = np.sum(bin_hist[3:])\n", + "long_tail_tile_pairs / np.sum(bin_hist)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "hipscatenv", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}