diff --git a/bpd/jackknife.py b/bpd/jackknife.py new file mode 100644 index 0000000..c6637c4 --- /dev/null +++ b/bpd/jackknife.py @@ -0,0 +1,59 @@ +from math import ceil +from typing import Callable + +import jax.numpy as jnp +from jax import Array, random +from tqdm import tqdm + + +def run_jackknife_shear_pipeline( + rng_key, + init_g: Array, + post_params_pos: dict, + post_params_neg: dict, + shear_pipeline: Callable, + n_jacks: int = 10, + disable_bar: bool = True, +): + """Use Jackknife to estimate the mean and std of the shear posterior. + + Args: + rng_key: Random jax key. + init_g: Initial value for shear `g`. + post_params_pos: Interim posterior galaxy parameters estimated using positive shear. + post_params_neg: Interim posterior galaxy parameters estimated using negative shear, + and otherwise same conditions and random seed as `post_params_pos`. + shear_pipeline: Function that outputs shear posterior samples from `post_params` with all + keyword arguments pre-specified. + n_jacks: Number of jackknife batches. + + Returns: + Jackknife + + """ + N, _ = post_params_pos["e1"].shape # N = n_gals, K = n_samples_per_gal + batch_size = ceil(N / n_jacks) + + g_best_list = [] + keys = random.split(rng_key, n_jacks) + + for ii in tqdm(range(n_jacks), desc="Jackknife #", disable=disable_bar): + k_ii = keys[ii] + start, end = ii * batch_size, (ii + 1) * batch_size + + _params_jack_pos = { + k: jnp.concatenate([v[:start], v[end:]]) for k, v in post_params_pos.items() + } + _params_jack_neg = { + k: jnp.concatenate([v[:start], v[end:]]) for k, v in post_params_neg.items() + } + + g_pos_ii = shear_pipeline(k_ii, _params_jack_pos, init_g) + g_neg_ii = shear_pipeline(k_ii, _params_jack_neg, -init_g) + g_best_ii = (g_pos_ii - g_neg_ii) * 0.5 + g_best_mean_ii = g_best_ii.mean(axis=0) + + g_best_list.append(g_best_mean_ii) + + g_best_means = jnp.array(g_best_list) + return g_best_means diff --git a/notebooks/check-g-mag-random.ipynb b/notebooks/old/check-g-mag-random.ipynb similarity index 100% rename from notebooks/check-g-mag-random.ipynb rename to notebooks/old/check-g-mag-random.ipynb diff --git a/notebooks/check-likelihood-toy1.ipynb b/notebooks/old/check-likelihood-toy1.ipynb similarity index 100% rename from notebooks/check-likelihood-toy1.ipynb rename to notebooks/old/check-likelihood-toy1.ipynb diff --git a/notebooks/check-new-prior1.ipynb b/notebooks/old/check-new-prior1.ipynb similarity index 100% rename from notebooks/check-new-prior1.ipynb rename to notebooks/old/check-new-prior1.ipynb diff --git a/notebooks/normalized-prior-check1.ipynb b/notebooks/old/normalized-prior-check1.ipynb similarity index 100% rename from notebooks/normalized-prior-check1.ipynb rename to notebooks/old/normalized-prior-check1.ipynb diff --git a/notebooks/shape-noise-cancellation-draft1.ipynb b/notebooks/shape-noise-cancellation-draft1.ipynb new file mode 100644 index 0000000..f822632 --- /dev/null +++ b/notebooks/shape-noise-cancellation-draft1.ipynb @@ -0,0 +1,583 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import os\n", + "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"0\"\n", + "os.environ[\"JAX_ENABLE_X64\"] = \"True\"" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from math import ceil\n", + "from tqdm import tqdm\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import jax.numpy as jnp\n", + "from jax import random\n", + "\n", + "\n", + "from bpd.pipelines.toy_ellips import pipeline_toy_ellips_samples\n", + "from bpd.pipelines.shear_inference import pipeline_shear_inference_ellipticities" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "seed = 42\n", + "key = random.key(seed)\n", + "\n", + "g1 = 0.02\n", + "g2 = 0.0\n", + "true_g = jnp.array([g1, g2])\n", + "\n", + "sigma_e = 1e-3\n", + "sigma_e_int = 4e-2\n", + "sigma_m = 1e-5\n", + "n_gals = 1000\n", + "n_samples_per_gal = 50" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "k1, k2 = random.split(key)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# positive shear\n", + "e_post_pos, _, _ = pipeline_toy_ellips_samples(k1, g1=g1, g2=g2, sigma_e=sigma_e, sigma_e_int=sigma_e_int, sigma_m=sigma_m, n_gals=n_gals, n_samples_per_gal=n_samples_per_gal)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# negative shear (same key!)\n", + "e_post_neg, _, _ = pipeline_toy_ellips_samples(k1, g1=-g1, g2=-g2, sigma_e=sigma_e, sigma_e_int=sigma_e_int, sigma_m=sigma_m, n_gals=n_gals, n_samples_per_gal=n_samples_per_gal)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(1000, 50, 2)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "e_post_pos.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "g_pos_samples = pipeline_shear_inference_ellipticities(k2, e_post_pos, true_g, sigma_e=sigma_e, sigma_e_int=sigma_e_int, n_samples=1000, initial_step_size=1e-3)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "g_neg_samples = pipeline_shear_inference_ellipticities(k2, e_post_neg, -true_g, sigma_e=sigma_e, sigma_e_int=sigma_e_int, n_samples=1000, initial_step_size=1e-3)" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(1000, 2)" + ] + }, + "execution_count": 60, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_pos_samples.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Array(3.37190099e-05, dtype=float64)" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_pos_samples[:, 0].std()" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([ 6., 4., 9., 11., 17., 34., 42., 58., 73., 70., 72., 91., 80.,\n", + " 82., 73., 74., 61., 48., 37., 23., 16., 9., 5., 4., 1.]),\n", + " array([-0.0201118 , -0.02010414, -0.02009648, -0.02008883, -0.02008117,\n", + " -0.02007351, -0.02006586, -0.0200582 , -0.02005054, -0.02004288,\n", + " -0.02003523, -0.02002757, -0.02001991, -0.02001226, -0.0200046 ,\n", + " -0.01999694, -0.01998928, -0.01998163, -0.01997397, -0.01996631,\n", + " -0.01995866, -0.019951 , -0.01994334, -0.01993568, -0.01992803,\n", + " -0.01992037]),\n", + " )" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAJGCAYAAACUbbvDAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJ3VJREFUeJzt3Xl0lfWd+PFP2EIQEgQkQA2L1cENxaUgarVTcwTLdHBkZtRqB62DVXGlVeGMYLVWwFp1tIqtVVyqUvW4te5D1RmPCIJoXXAtKFYTtzFBUGT5/v5wuD8vBOSLIWF5vc6555jnPve53+9zn5C3T+6TW5JSSgEAwDpr0dwDAADY1AgoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACBTq+YewKpWrFgR77zzTnTo0CFKSkqaezgAwGYupRQLFy6MHj16RIsW63ZuaaMLqHfeeSeqqqqaexgAwBZmwYIFse22267TuhtdQHXo0CEivphEeXl5M48GANjc1dfXR1VVVaFB1sVGF1Arf21XXl4uoACAJpPz1iFvIgcAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATK2aewDA5qn3mPsadXvzJw5t1O0BfB3OQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZWjX3AACaWu8x9zXq9uZPHNqo2wM2fs5AAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZsgJq+fLlMW7cuOjTp0+UlZXFN7/5zfj5z38eKaXCOimlGD9+fHTv3j3Kysqiuro6XnvttUYfOABAc8kKqEmTJsXkyZPj17/+dcydOzcmTZoUF110UVxxxRWFdS666KK4/PLL4+qrr44ZM2bEVlttFYMHD47PPvus0QcPANAcsv6Q5pNPPhnDhg2LoUO/+KNxvXv3jltvvTVmzpwZEV+cfbrsssvinHPOiWHDhkVExI033hiVlZVx9913xxFHHNHIwwcAaHpZZ6D23XffmDZtWrz66qsREfHcc8/FE088EYccckhERMybNy9qamqiurq68JiKiooYOHBgTJ8+vcFtLlmyJOrr64tuAAAbs6wzUGPGjIn6+vrYcccdo2XLlrF8+fL4xS9+EUcddVRERNTU1ERERGVlZdHjKisrC/etasKECXHeeeetz9iBLUhjf/wKwNeRdQbqtttui5tvvjluueWWeOaZZ+KGG26Iiy++OG644Yb1HsDYsWOjrq6ucFuwYMF6bwsAoClknYE688wzY8yYMYX3MvXr1y/efPPNmDBhQowYMSK6desWERG1tbXRvXv3wuNqa2ujf//+DW6ztLQ0SktL13P4AABNL+sM1OLFi6NFi+KHtGzZMlasWBEREX369Ilu3brFtGnTCvfX19fHjBkzYtCgQY0wXACA5pd1Bur73/9+/OIXv4iePXvGLrvsEnPmzIlLLrkkfvSjH0VERElJSZx++ulxwQUXxA477BB9+vSJcePGRY8ePeLQQw/dEOMHAGhyWQF1xRVXxLhx4+Kkk06K9957L3r06BE//vGPY/z48YV1zjrrrFi0aFEcf/zx8fHHH8f+++8fDz74YLRt27bRBw8A0BxK0pf/jPhGoL6+PioqKqKuri7Ky8ubezjAetqSrpqbP3Focw8B+BrWpz18Fh4AQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGRq1dwDAKBY7zH3Ner25k8c2qjbA5yBAgDIJqAAADIJKACATAIKACCTN5FDE/Lm4M1TY7+uwMbPGSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMrsIDIsKVZAA5nIECAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAy+SgX2IT5+BWA5uEMFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJCpVXMPAIANq/eY+xptW/MnDm20bcGmzBkoAIBM2QH1t7/9LY4++ujo3LlzlJWVRb9+/WLWrFmF+1NKMX78+OjevXuUlZVFdXV1vPbaa406aACA5pQVUP/7v/8b++23X7Ru3ToeeOCBeOmll+JXv/pVbL311oV1Lrroorj88svj6quvjhkzZsRWW20VgwcPjs8++6zRBw8A0Byy3gM1adKkqKqqiilTphSW9enTp/DfKaW47LLL4pxzzolhw4ZFRMSNN94YlZWVcffdd8cRRxzRSMMGAGg+WWeg7r333th7773jX/7lX6Jr166xxx57xDXXXFO4f968eVFTUxPV1dWFZRUVFTFw4MCYPn16g9tcsmRJ1NfXF90AADZmWQH117/+NSZPnhw77LBDPPTQQ3HiiSfGqaeeGjfccENERNTU1ERERGVlZdHjKisrC/etasKECVFRUVG4VVVVrc88AACaTFZArVixIvbcc8+48MILY4899ojjjz8+Ro4cGVdfffV6D2Ds2LFRV1dXuC1YsGC9twUA0BSyAqp79+6x8847Fy3baaed4q233oqIiG7dukVERG1tbdE6tbW1hftWVVpaGuXl5UU3AICNWVZA7bfffvHKK68ULXv11VejV69eEfHFG8q7desW06ZNK9xfX18fM2bMiEGDBjXCcAEAml/WVXhnnHFG7LvvvnHhhRfGv/7rv8bMmTPjt7/9bfz2t7+NiIiSkpI4/fTT44ILLogddtgh+vTpE+PGjYsePXrEoYceuiHGDwDQ5LIC6lvf+lbcddddMXbs2Dj//POjT58+cdlll8VRRx1VWOess86KRYsWxfHHHx8ff/xx7L///vHggw9G27ZtG33wAADNoSSllJp7EF9WX18fFRUVUVdX5/1QbHYa8zPJoDn4LDw2R+vTHj4LDwAgk4ACAMgkoAAAMgkoAIBMWVfhwZbGm74BaIgzUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZGrV3AMAYNPRe8x9jbq9+ROHNur2oKk4AwUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQKZWzT0AaEy9x9zX3EMAYAvgDBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBk+loBNXHixCgpKYnTTz+9sOyzzz6LUaNGRefOnaN9+/YxfPjwqK2t/brjBADYaKx3QD399NPxm9/8Jnbbbbei5WeccUb88Y9/jNtvvz0ef/zxeOedd+Kwww772gMFANhYrFdAffLJJ3HUUUfFNddcE1tvvXVheV1dXVx77bVxySWXxHe/+93Ya6+9YsqUKfHkk0/GU0891WiDBgBoTusVUKNGjYqhQ4dGdXV10fLZs2fH0qVLi5bvuOOO0bNnz5g+fXqD21qyZEnU19cX3QAANmbZn4U3derUeOaZZ+Lpp59e7b6amppo06ZNdOzYsWh5ZWVl1NTUNLi9CRMmxHnnnZc7DACAZpN1BmrBggVx2mmnxc033xxt27ZtlAGMHTs26urqCrcFCxY0ynYBADaUrICaPXt2vPfee7HnnntGq1atolWrVvH444/H5ZdfHq1atYrKysr4/PPP4+OPPy56XG1tbXTr1q3BbZaWlkZ5eXnRDQBgY5b1K7yDDjoonn/++aJlxx57bOy4445x9tlnR1VVVbRu3TqmTZsWw4cPj4iIV155Jd56660YNGhQ440aAKAZZQVUhw4dYtdddy1attVWW0Xnzp0Ly4877rgYPXp0dOrUKcrLy+OUU06JQYMGxT777NN4owYAaEbZbyL/Kpdeemm0aNEihg8fHkuWLInBgwfHVVdd1dhPAwDQbEpSSqm5B/Fl9fX1UVFREXV1dd4PRbbeY+5r7iEAGeZPHNrcQ4D1ag+fhQcAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGRq9M/CA4B11dgfv+SjYWgqzkABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAECmVs09AOg95r7mHgIAZHEGCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMjUqrkHAACNpfeY+xp1e/MnDm3U7bH5cAYKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATFkBNWHChPjWt74VHTp0iK5du8ahhx4ar7zyStE6n332WYwaNSo6d+4c7du3j+HDh0dtbW2jDhoAoDllBdTjjz8eo0aNiqeeeioeeeSRWLp0aRx88MGxaNGiwjpnnHFG/PGPf4zbb789Hn/88XjnnXfisMMOa/SBAwA0l1Y5Kz/44INFX19//fXRtWvXmD17dhxwwAFRV1cX1157bdxyyy3x3e9+NyIipkyZEjvttFM89dRTsc8++6y2zSVLlsSSJUsKX9fX16/PPAAAmszXeg9UXV1dRER06tQpIiJmz54dS5cujerq6sI6O+64Y/Ts2TOmT5/e4DYmTJgQFRUVhVtVVdXXGRIAwAa33gG1YsWKOP3002O//faLXXfdNSIiampqok2bNtGxY8eidSsrK6OmpqbB7YwdOzbq6uoKtwULFqzvkAAAmkTWr/C+bNSoUfHCCy/EE0888bUGUFpaGqWlpV9rGwAATWm9zkCdfPLJ8ac//SkeffTR2HbbbQvLu3XrFp9//nl8/PHHRevX1tZGt27dvtZAAQA2FlkBlVKKk08+Oe66667485//HH369Cm6f6+99orWrVvHtGnTCsteeeWVeOutt2LQoEGNM2IAgGaW9Su8UaNGxS233BL33HNPdOjQofC+poqKiigrK4uKioo47rjjYvTo0dGpU6coLy+PU045JQYNGtTgFXgAAJuirICaPHlyRER85zvfKVo+ZcqUOOaYYyIi4tJLL40WLVrE8OHDY8mSJTF48OC46qqrGmWwAAAbg6yASil95Tpt27aNK6+8Mq688sr1HhQAwMbMZ+EBAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZMr6S+QQEdF7zH3NPQQAaFbOQAEAZBJQAACZBBQAQCYBBQCQSUABAGRyFR4ArEFjXnU8f+LQRtsWzc8ZKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyNSquQfAhtd7zH3NPQQA2Kw4AwUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmXwW3kbIZ9cBbH4a+9/2+ROHNur2yOMMFABAJgEFAJBJQAEAZBJQAACZBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJBJQAEAZGrV3ANoLr3H3Neo25s/cWijbg8A2Hg5AwUAkElAAQBkElAAAJkEFABAJgEFAJBpi70KDwA2Za4mb17OQAEAZBJQAACZBBQAQCYBBQCQSUABAGRyFV4jaeyrIQCgKTXmz7Et4Yo+Z6AAADIJKACATAIKACCTgAIAyCSgAAAyCSgAgEwCCgAgk4ACAMgkoAAAMgkoAIBMAgoAIJOAAgDIJKAAADIJKACATAIKACCTgAIAyNSquQcAAGxeeo+5r1G3N3/i0EbdXmPYYGegrrzyyujdu3e0bds2Bg4cGDNnztxQTwUA0KQ2SED94Q9/iNGjR8e5554bzzzzTOy+++4xePDgeO+99zbE0wEANKkNElCXXHJJjBw5Mo499tjYeeed4+qrr4527drFddddtyGeDgCgSTX6e6A+//zzmD17dowdO7awrEWLFlFdXR3Tp09fbf0lS5bEkiVLCl/X1dVFRER9fX1jD63IiiWLN+j2AYDGsaGbYOX2U0rr/JhGD6gPPvggli9fHpWVlUXLKysr4+WXX15t/QkTJsR555232vKqqqrGHhoAsAmquKxpnmfhwoVRUVGxTus2+1V4Y8eOjdGjRxe+XrFiRXz00UfRuXPnKCkpacaRNb76+vqoqqqKBQsWRHl5eXMPp1ls6ftgS59/hH1g/lv2/CPsg41x/imlWLhwYfTo0WOdH9PoAdWlS5do2bJl1NbWFi2vra2Nbt26rbZ+aWlplJaWFi3r2LFjYw9ro1JeXr7RHDTNZUvfB1v6/CPsA/PfsucfYR9sbPNf1zNPKzX6m8jbtGkTe+21V0ybNq2wbMWKFTFt2rQYNGhQYz8dAECT2yC/whs9enSMGDEi9t577xgwYEBcdtllsWjRojj22GM3xNMBADSpDRJQhx9+eLz//vsxfvz4qKmpif79+8eDDz642hvLtzSlpaVx7rnnrvYryy3Jlr4PtvT5R9gH5r9lzz/CPthc5l+Scq7ZAwDAhwkDAOQSUAAAmQQUAEAmAQUAkElAAQBkElD/56OPPoqjjjoqysvLo2PHjnHcccfFJ598stbHfPbZZzFq1Kjo3LlztG/fPoYPH170F9ife+65OPLII6OqqirKyspip512iv/8z/9cbTuPPfZY7LnnnlFaWhrbb799XH/99UX3//d//3d8//vfjx49ekRJSUncfffdq20jpRTjx4+P7t27R1lZWVRXV8drr72WNceNeR9ERFx55ZXRu3fvaNu2bQwcODBmzpxZuG/+/PlRUlLS4O32228vrNfQ/VOnTt3k5x8R8Z3vfGe1uZ1wwglF67z11lsxdOjQaNeuXXTt2jXOPPPMWLZs2SY//48++ihOOeWU6Nu3b5SVlUXPnj3j1FNPLXw4+Upre/039X2wLmOJaPpjICLi1FNPjb322itKS0ujf//+DW7ntttui/79+0e7du2iV69e8ctf/rLB+e+0005RVlYWffv2jRtvvLHo/oa+B0pKSmLo0KGFdY455pjV7h8yZEjh/k19H1x//fWrza9t27ZF66zt58WmPv9rrrkmvv3tb8fWW28dW2+9dVRXV6/2ffJVx8A6S6SUUhoyZEjafffd01NPPZX+53/+J22//fbpyCOPXOtjTjjhhFRVVZWmTZuWZs2alfbZZ5+07777Fu6/9tpr06mnnpoee+yx9MYbb6SbbroplZWVpSuuuKKwzl//+tfUrl27NHr06PTSSy+lK664IrVs2TI9+OCDhXXuv//+9B//8R/pzjvvTBGR7rrrrtXGMnHixFRRUZHuvvvu9Nxzz6V//Md/TH369EmffvrpOs9xY94HU6dOTW3atEnXXXddevHFF9PIkSNTx44dU21tbUoppWXLlqV333236Hbeeeel9u3bp4ULFxa2ExFpypQpReut3Eeb8vxTSunAAw9MI0eOLJpbXV1d4f5ly5alXXfdNVVXV6c5c+ak+++/P3Xp0iWNHTt2k5//888/nw477LB07733ptdffz1NmzYt7bDDDmn48OFF413b67+p74N1GUtzHAMppXTKKaekX//61+mHP/xh2n333Vfbxv33359atWqVJk+enN544430pz/9KXXv3r1oH1111VWpQ4cOaerUqemNN95It956a2rfvn269957C+t8+OGHRa/tCy+8kFq2bJmmTJlSWGfEiBFpyJAhRet99NFHhfs39X0wZcqUVF5eXjS/mpqaouda28+LTX3+P/jBD9KVV16Z5syZk+bOnZuOOeaYVFFRkd5+++3COl91DKwrAZVSeumll1JEpKeffrqw7IEHHkglJSXpb3/7W4OP+fjjj1Pr1q3T7bffXlg2d+7cFBFp+vTpa3yuk046Kf393/994euzzjor7bLLLkXrHH744Wnw4MENPr6hgFqxYkXq1q1b+uUvf1k0vtLS0nTrrbeu0xw39n0wYMCANGrUqMLXy5cvTz169EgTJkxY4/P0798//ehHPypatqYA3Rzmf+CBB6bTTjttjc97//33pxYtWhT9Yzp58uRUXl6enn322U1+/qu67bbbUps2bdLSpUsLy9b0+qe06R8D6zKW5j4Gzj333AZ/eB555JHpn//5n4uWXX755WnbbbdNK1asSCmlNGjQoPTTn/60aJ3Ro0en/fbbr8GxpZTSpZdemjp06JA++eSTwrIRI0akYcOGNbh+UxwDG3ofTJkyJVVUVDQ41pTW/vPi4osv3uTnv6ply5alDh06pBtuuKGwbG3HQA6/wouI6dOnR8eOHWPvvfcuLKuuro4WLVrEjBkzGnzM7NmzY+nSpVFdXV1YtuOOO0bPnj1j+vTpa3yuurq66NSpU9Fzf3kbERGDBw9e6zZWNW/evKipqSnaTkVFRQwcOLCwna+a48a8Dz7//POYPXt20TotWrSI6urqNT7P7Nmz49lnn43jjjtutftGjRoVXbp0iQEDBsR1110XKaXNZv4333xzdOnSJXbdddcYO3ZsLF68uOh5+vXrV/SJAIMHD476+vq48847N4v5r/o85eXl0apV8QcuNPT6rxzHprwP1mUsG8sxsKolS5as9mumsrKyePvtt+PNN99c6zozZ86MpUuXNrjda6+9No444ojYaqutipY/9thj0bVr1+jbt2+ceOKJ8eGHH0ZE0x4Dq2rMffDJJ59Er169oqqqKoYNGxYvvvhi4b61/bx44IEHNov5f9nixYtj6dKlRd9vEWs+BnIIqIioqamJrl27Fi1r1apVdOrUKWpqatb4mDZt2kTHjh2LlldWVq7xMU8++WT84Q9/iOOPP75oO6t+xE1lZWXU19fHp59+us7jX/m4NY3lq+a4Me+DDz74IJYvX77W+a3q2muvjZ122in23XffouXnn39+3HbbbfHII4/E8OHD46STToorrrhis5j/D37wg/j9738fjz76aIwdOzZuuummOProo7/yeSIi3nzzzU1+/l/2wQcfxM9//vOi54lY8+u/chyb8j5Yl7FsDMdAQwYPHhx33nlnTJs2LVasWBGvvvpq/OpXv4qIiHfffbewzu9+97uYPXt2pJRi1qxZ8bvf/S6WLl0aH3zwwWrbnDlzZrzwwgvx7//+70XLhwwZEjfeeGNMmzYtJk2aFI8//ngccsghsXz58iY7BjbkPujbt29cd911cc8998Tvf//7WLFiRey7777x9ttvF8a7cnyrjre2tnaTn/+qzj777OjRo0dR4K3tGMixQT4Lb2MxZsyYmDRp0lrXmTt3bpOM5YUXXohhw4bFueeeGwcffHCTPGdExF/+8pcoKSkpfP3l/16ppqamST6TqKn2waeffhq33HJLjBs3rsFjYPz48UVfX3jhhXHqqadusPGstKHn/+UfyP369Yvu3bvHQQcdtNpr3tAx0BSa6vWvr6+PoUOHxs477xw/+9nPVjsGVn39TzvttCb7nmyOfwcWLFgQF110UVx00UWFZc11DKzJyJEj44033oh/+Id/iKVLl0Z5eXmcdtpp8bOf/SxatPji//PHjRsXNTU1sc8++0RKKSorK2PEiBFx0UUXFdb5smuvvTb69esXAwYMWKefBau+GbmpNdY+GDRoUAwaNKiw3X333Te6du0aVVVVRc/Xo0ePoq8HDx68gWe4dhviGJg4cWJMnTo1HnvssaIzV0cccUThv/v16xe77bZbfPOb34zHHnssDjrooHUe82Z9BuonP/lJzJ07d6237bbbLrp16xbvvfde0WOXLVsWH330UXTr1q3BbXfr1i0+//zz+Pjjj4uW19bWrvaYl156KQ466KA4/vjj45xzzlltO6terVBbWxvl5eVRVla2TvNc+XwNbWf//fePuXPnxgUXXBAdOnQomvvzzz8fLVq0iF122WWj3gddunSJli1bNrhOQ2O74447YvHixfFv//ZvX3kMTJ48OWpra6Nz586bzfxXGjhwYER8cVXK3Llz48QTT4y+ffsWzf+RRx6JiIhevXptFvNfuHBhDBkyJDp06BB33XVXtG7deq3HwOTJkyMi4hvf+MYm/z2wprG0bNkyzj777GY/BtampKQkJk2aFJ988km8+eabUVNTEwMGDIiIiO222y4ivvhVzXXXXReLFy+O+fPnx1tvvRW9e/eODh06xDbbbFO0vUWLFsXUqVMLv8L/qn8HOnbsGIsXL97gx0BT7oOVWrduHd/+9rfje9/7XsydOzcefvjhiIi48847i/bBokWLorKycrOZ/8UXXxwTJ06Mhx9+OHbbbbe1Pvd2220XXbp0iddff32dxxsRrsJL6f+/cXDWrFmFZQ899NA6vXHujjvuKCx7+eWXV3vj3AsvvJC6du2azjzzzAa3c9ZZZ6Vdd921aNmRRx65Xm8iv/jiiwvL6urqGnwT+ZrmuLHvgwEDBqSTTz658PXy5cvTN77xjQbfRHzggQeudvXVmlxwwQVp66233qzmv9ITTzyRIiI999xzKaX//wbiL1+19Zvf/KboDcSb8vzr6urSPvvskw488MC0aNGiNe6XL1v5+qe08f878FX7YF3G0lzHwEpregNxQ374wx+mQYMGrXWdAw44oMErxKZMmZJKS0vTBx988JXPs2DBglRSUpLuueeeDXoMrNRU+2ClZcuWpb59+6YzzjgjpbT2nxcr30S+qc9/0qRJqby8fK0XcnzZl4+BHALq/wwZMiTtscceacaMGemJJ55IO+ywQ9GL8vbbb6e+ffumGTNmFJadcMIJqWfPnunPf/5zmjVrVho0aFDRi/3888+nbbbZJh199NFFl0u+9957hXVWXr585plnprlz56Yrr7xytcuXFy5cmObMmZPmzJmTIiJdcsklac6cOenNN98srDNx4sTUsWPHdM8996S//OUvadiwYQ3+GYO1zXFj3gdTp05NpaWl6frrr08vvfRSOv7441PHjh1Xuzz3tddeSyUlJemBBx5Y7TW+99570zXXXJOef/759Nprr6WrrroqtWvXLo0fP36Tn//rr7+ezj///DRr1qw0b968dM8996TtttsuHXDAAYVtrLyE/eCDD07PPvtsevDBB9M222xTdAn7pjr/urq6NHDgwNSvX7/0+uuvFz3XsmXL1un139T3wbqMpTmOgZS++L6cM2dO+vGPf5z+7u/+rvDv2ZIlS1JKKb3//vtp8uTJae7cuWnOnDnp1FNPTW3bti16nldeeSXddNNN6dVXX00zZsxIhx9+eOrUqVOaN29eWtX++++fDj/88NWWL1y4MP30pz9N06dPT/PmzUv/9V//lfbcc8+0ww47pM8++2yz2AfnnXdeeuihh9Ibb7yRZs+enY444ojUtm3b9OKLLxbWWdvPi019/hMnTkxt2rRJd9xxR9H328o/Z7Mux8C6ElD/58MPP0xHHnlkat++fSovL0/HHnts0d8PmjdvXoqI9OijjxaWffrpp+mkk05KW2+9dWrXrl36p3/6p/Tuu+8W7j/33HNTRKx269WrV9FzP/roo6l///6pTZs2abvttiv6myUr729oOyNGjCiss2LFijRu3LhUWVmZSktL00EHHZReeeWVrDluzPsgpZSuuOKK1LNnz9SmTZs0YMCA9NRTT622ztixY1NVVVVavnz5avc98MADqX///ql9+/Zpq622Srvvvnu6+uqrC+tuyvN/66230gEHHJA6deqUSktL0/bbb5/OPPPMor8DlVJK8+fPT4ccckgqKytLXbp0ST/5yU8Kl/lvyvNf0/dIRBT+cf2q139T3wfrMpaUmv4YSOmLs8Jre23ef//9tM8++6StttoqtWvXLh100EGrze2ll15K/fv3T2VlZam8vDwNGzYsvfzyy6vto5VnPx5++OHV7lu8eHE6+OCD0zbbbJNat26devXqlUaOHFkUoZv6Pjj99NMLx0hlZWX63ve+l5555pmiddb282JTn3+vXr0afJ5zzz03pbRux8C6Kknp/67hBQBgnWzWbyIHANgQBBQAQCYBBQCQSUABAGQSUAAAmQQUAEAmAQUAkElAAQBkElAAAJkEFABAJgEFAJDp/wFiBHz9DIjLbwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# plt.hist(g_pos_samples[:, 0], bins=25)\n", + "plt.figure(figsize=(7,7))\n", + "plt.hist(g_neg_samples[:, 0], bins=25)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "g_best_samples = (g_pos_samples - g_neg_samples) * 0.5" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([ 1., 1., 2., 3., 2., 1., 6., 3., 8., 7., 24.,\n", + " 33., 58., 129., 210., 292., 131., 45., 17., 12., 4., 4.,\n", + " 3., 2., 2.]),\n", + " array([0.01999981, 0.01999981, 0.01999981, 0.01999981, 0.01999981,\n", + " 0.01999981, 0.01999981, 0.01999981, 0.01999981, 0.01999981,\n", + " 0.01999981, 0.01999981, 0.01999981, 0.01999981, 0.01999981,\n", + " 0.01999981, 0.01999981, 0.01999981, 0.01999981, 0.01999981,\n", + " 0.01999981, 0.01999981, 0.01999981, 0.01999981, 0.01999981,\n", + " 0.01999981]),\n", + " )" + ] + }, + "execution_count": 36, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGvCAYAAABxUC54AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALRtJREFUeJzt3Xl0VGWexvEnC6kkQCUGSSqREAHFEGUTIZQLYrMEiLS2uOBERKFh8AQbiIOQEUG0FQU9tgtKd58ZcCGy9FFaUKExQBgwgkRpEJQGRFkrQZEUgSGQ5J0/+lDTJWGpkJA3xfdzzj2Huve9t36/XIo8vLduVYgxxggAAMAiofVdAAAAwC8RUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1gmv7wJqoqqqSvv371fTpk0VEhJS3+UAAIDzYIzRkSNHlJSUpNDQs8+RNMiAsn//fiUnJ9d3GQAAoAb27NmjFi1anHVMgwwoTZs2lfTPBp1OZz1XAwAAzofX61VycrLv9/jZNMiAcuqyjtPpJKAAANDAnM/bMwJ6k+ybb76pDh06+IKB2+3WJ5984tt+/PhxZWdnq1mzZmrSpIkGDRqk4uJiv2Ps3r1bmZmZio6OVnx8vMaPH6+KiopAygAAAEEuoIDSokULPf/88yoqKtKGDRv0q1/9SnfccYe2bNkiSRo3bpwWL16shQsXqqCgQPv379ddd93l27+yslKZmZk6ceKEPvvsM7311luaM2eOJk+eXLtdAQCABi3EGGMu5ABxcXGaMWOG7r77bjVv3lx5eXm6++67JUnffvut2rVrp8LCQnXv3l2ffPKJbr/9du3fv18JCQmSpFmzZmnChAk6ePCgIiIizus5vV6vYmJiVFpayiUeAAAaiEB+f9f4c1AqKys1b948HT16VG63W0VFRTp58qR69+7tG5OamqqWLVuqsLBQklRYWKj27dv7wokkZWRkyOv1+mZhqlNeXi6v1+u3AACA4BVwQNm8ebOaNGkih8OhUaNG6YMPPlBaWpo8Ho8iIiIUGxvrNz4hIUEej0eS5PF4/MLJqe2ntp3JtGnTFBMT41u4xRgAgOAWcEC55pprtHHjRq1bt06PPPKIhg4dqq1bt9ZFbT65ubkqLS31LXv27KnT5wMAAPUr4NuMIyIidNVVV0mSunTpoi+++EKvvPKK7rvvPp04cUKHDx/2m0UpLi6Wy+WSJLlcLq1fv97veKfu8jk1pjoOh0MOhyPQUgEAQAN1wd/FU1VVpfLycnXp0kWNGjVSfn6+b9u2bdu0e/duud1uSZLb7dbmzZtVUlLiG7N8+XI5nU6lpaVdaCkAACBIBDSDkpubq/79+6tly5Y6cuSI8vLytGrVKi1btkwxMTEaPny4cnJyFBcXJ6fTqUcffVRut1vdu3eXJPXt21dpaWkaMmSIpk+fLo/Ho0mTJik7O5sZEgAA4BNQQCkpKdGDDz6oAwcOKCYmRh06dNCyZcvUp08fSdLLL7+s0NBQDRo0SOXl5crIyNAbb7zh2z8sLExLlizRI488IrfbrcaNG2vo0KF6+umna7crAADQoF3w56DUBz4HBQCAhueifA4KAABAXSGgAAAA6xBQAACAdQL+HBQAwOmunPhRrRzn++cza+U4QEPHDAoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1Agoo06ZNU9euXdW0aVPFx8frzjvv1LZt2/zG9OzZUyEhIX7LqFGj/Mbs3r1bmZmZio6OVnx8vMaPH6+KiooL7wYAAASF8EAGFxQUKDs7W127dlVFRYX+8z//U3379tXWrVvVuHFj37gRI0bo6aef9j2Ojo72/bmyslKZmZlyuVz67LPPdODAAT344INq1KiRnnvuuVpoCQAANHQBBZSlS5f6PZ4zZ47i4+NVVFSkHj16+NZHR0fL5XJVe4y//e1v2rp1qz799FMlJCSoU6dOeuaZZzRhwgQ99dRTioiIqEEbAAAgmFzQe1BKS0slSXFxcX7r586dq8svv1zXXXedcnNzdezYMd+2wsJCtW/fXgkJCb51GRkZ8nq92rJlS7XPU15eLq/X67cAAIDgFdAMyr+qqqrS2LFjddNNN+m6667zrf+3f/s3paSkKCkpSZs2bdKECRO0bds2vf/++5Ikj8fjF04k+R57PJ5qn2vatGmaOnVqTUsFAAANTI0DSnZ2tr7++mutWbPGb/3IkSN9f27fvr0SExPVq1cv7dy5U23atKnRc+Xm5ionJ8f32Ov1Kjk5uWaFAwAA69XoEs/o0aO1ZMkSrVy5Ui1atDjr2PT0dEnSjh07JEkul0vFxcV+Y049PtP7VhwOh5xOp98CAACCV0ABxRij0aNH64MPPtCKFSvUqlWrc+6zceNGSVJiYqIkye12a/PmzSopKfGNWb58uZxOp9LS0gIpBwAABKmALvFkZ2crLy9Pf/3rX9W0aVPfe0ZiYmIUFRWlnTt3Ki8vTwMGDFCzZs20adMmjRs3Tj169FCHDh0kSX379lVaWpqGDBmi6dOny+PxaNKkScrOzpbD4aj9DgEAQIMT0AzKm2++qdLSUvXs2VOJiYm+Zf78+ZKkiIgIffrpp+rbt69SU1P12GOPadCgQVq8eLHvGGFhYVqyZInCwsLkdrv1wAMP6MEHH/T73BQAAHBpC2gGxRhz1u3JyckqKCg453FSUlL08ccfB/LUAADgEsJ38QAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrhNd3AQBQn66c+FF9lwCgGsygAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsE1BAmTZtmrp27aqmTZsqPj5ed955p7Zt2+Y35vjx48rOzlazZs3UpEkTDRo0SMXFxX5jdu/erczMTEVHRys+Pl7jx49XRUXFhXcDAACCQkABpaCgQNnZ2fr888+1fPlynTx5Un379tXRo0d9Y8aNG6fFixdr4cKFKigo0P79+3XXXXf5tldWViozM1MnTpzQZ599prfeektz5szR5MmTa68rAADQoIUYY0xNdz548KDi4+NVUFCgHj16qLS0VM2bN1deXp7uvvtuSdK3336rdu3aqbCwUN27d9cnn3yi22+/Xfv371dCQoIkadasWZowYYIOHjyoiIiIcz6v1+tVTEyMSktL5XQ6a1o+AFj3SbLfP59Z3yUAdSaQ398X9B6U0tJSSVJcXJwkqaioSCdPnlTv3r19Y1JTU9WyZUsVFhZKkgoLC9W+fXtfOJGkjIwMeb1ebdmypdrnKS8vl9fr9VsAAEDwqnFAqaqq0tixY3XTTTfpuuuukyR5PB5FREQoNjbWb2xCQoI8Ho9vzL+Gk1PbT22rzrRp0xQTE+NbkpOTa1o2AABoAGocULKzs/X1119r3rx5tVlPtXJzc1VaWupb9uzZU+fPCQAA6k+Nvs149OjRWrJkiVavXq0WLVr41rtcLp04cUKHDx/2m0UpLi6Wy+XyjVm/fr3f8U7d5XNqzC85HA45HI6alAoAABqggGZQjDEaPXq0PvjgA61YsUKtWrXy296lSxc1atRI+fn5vnXbtm3T7t275Xa7JUlut1ubN29WSUmJb8zy5cvldDqVlpZ2Ib0AAIAgEdAMSnZ2tvLy8vTXv/5VTZs29b1nJCYmRlFRUYqJidHw4cOVk5OjuLg4OZ1OPfroo3K73erevbskqW/fvkpLS9OQIUM0ffp0eTweTZo0SdnZ2cySAAAASQEGlDfffFOS1LNnT7/1s2fP1kMPPSRJevnllxUaGqpBgwapvLxcGRkZeuONN3xjw8LCtGTJEj3yyCNyu91q3Lixhg4dqqeffvrCOgEAAEHjgj4Hpb7wOSgAagufgwJcPBftc1AAAADqAgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwTsABZfXq1Ro4cKCSkpIUEhKiRYsW+W1/6KGHFBIS4rf069fPb8yhQ4eUlZUlp9Op2NhYDR8+XGVlZRfUCAAACB4BB5SjR4+qY8eOmjlz5hnH9OvXTwcOHPAt7733nt/2rKwsbdmyRcuXL9eSJUu0evVqjRw5MvDqAQBAUAoPdIf+/furf//+Zx3jcDjkcrmq3fbNN99o6dKl+uKLL3TDDTdIkl577TUNGDBAL774opKSkgItCQAABJk6eQ/KqlWrFB8fr2uuuUaPPPKIfvrpJ9+2wsJCxcbG+sKJJPXu3VuhoaFat25dtccrLy+X1+v1WwAAQPCq9YDSr18/vf3228rPz9cLL7yggoIC9e/fX5WVlZIkj8ej+Ph4v33Cw8MVFxcnj8dT7TGnTZummJgY35KcnFzbZQMAAIsEfInnXAYPHuz7c/v27dWhQwe1adNGq1atUq9evWp0zNzcXOXk5Pgee71eQgoAAEGszm8zbt26tS6//HLt2LFDkuRyuVRSUuI3pqKiQocOHTrj+1YcDoecTqffAgAAgledB5S9e/fqp59+UmJioiTJ7Xbr8OHDKioq8o1ZsWKFqqqqlJ6eXtflAACABiDgSzxlZWW+2RBJ2rVrlzZu3Ki4uDjFxcVp6tSpGjRokFwul3bu3KnHH39cV111lTIyMiRJ7dq1U79+/TRixAjNmjVLJ0+e1OjRozV48GDu4AEAAJJqMIOyYcMGde7cWZ07d5Yk5eTkqHPnzpo8ebLCwsK0adMm/frXv1bbtm01fPhwdenSRf/zP/8jh8PhO8bcuXOVmpqqXr16acCAAbr55pv1pz/9qfa6AgAADVrAMyg9e/aUMeaM25ctW3bOY8TFxSkvLy/QpwYAAJcIvosHAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsE54fRcAAPh/V078qFaO8/3zmbVyHKC+MIMCAACsQ0ABAADW4RIPgAapti6FALATMygAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYJOKCsXr1aAwcOVFJSkkJCQrRo0SK/7cYYTZ48WYmJiYqKilLv3r21fft2vzGHDh1SVlaWnE6nYmNjNXz4cJWVlV1QIwAAIHgEHFCOHj2qjh07aubMmdVunz59ul599VXNmjVL69atU+PGjZWRkaHjx4/7xmRlZWnLli1avny5lixZotWrV2vkyJE17wIAAASV8EB36N+/v/r371/tNmOM/vCHP2jSpEm64447JElvv/22EhIStGjRIg0ePFjffPONli5dqi+++EI33HCDJOm1117TgAED9OKLLyopKekC2gEAAMGgVt+DsmvXLnk8HvXu3du3LiYmRunp6SosLJQkFRYWKjY21hdOJKl3794KDQ3VunXrqj1ueXm5vF6v3wIAAIJXrQYUj8cjSUpISPBbn5CQ4Nvm8XgUHx/vtz08PFxxcXG+Mb80bdo0xcTE+Jbk5OTaLBsAAFimQdzFk5ubq9LSUt+yZ8+e+i4JAADUoVoNKC6XS5JUXFzst764uNi3zeVyqaSkxG97RUWFDh065BvzSw6HQ06n028BAADBq1YDSqtWreRyuZSfn+9b5/V6tW7dOrndbkmS2+3W4cOHVVRU5BuzYsUKVVVVKT09vTbLAQAADVTAd/GUlZVpx44dvse7du3Sxo0bFRcXp5YtW2rs2LH6/e9/r6uvvlqtWrXSk08+qaSkJN15552SpHbt2qlfv34aMWKEZs2apZMnT2r06NEaPHgwd/AAAABJNQgoGzZs0G233eZ7nJOTI0kaOnSo5syZo8cff1xHjx7VyJEjdfjwYd18881aunSpIiMjffvMnTtXo0ePVq9evRQaGqpBgwbp1VdfrYV2AABAMAgxxpj6LiJQXq9XMTExKi0t5f0owCXqyokf1XcJVvv++cz6LgE4TSC/vxvEXTwAAODSQkABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA64TXdwEALi1XTvyovksA0AAwgwIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDq1HlCeeuophYSE+C2pqam+7cePH1d2draaNWumJk2aaNCgQSouLq7tMgAAQANWJzMo1157rQ4cOOBb1qxZ49s2btw4LV68WAsXLlRBQYH279+vu+66qy7KAAAADVSdfA5KeHi4XC7XaetLS0v1X//1X8rLy9OvfvUrSdLs2bPVrl07ff755+revXtdlAMAABqYOplB2b59u5KSktS6dWtlZWVp9+7dkqSioiKdPHlSvXv39o1NTU1Vy5YtVVhYeMbjlZeXy+v1+i0AACB41XpASU9P15w5c7R06VK9+eab2rVrl2655RYdOXJEHo9HERERio2N9dsnISFBHo/njMecNm2aYmJifEtycnJtlw0AACxS65d4+vfv7/tzhw4dlJ6erpSUFC1YsEBRUVE1OmZubq5ycnJ8j71eLyEFAIAgVue3GcfGxqpt27basWOHXC6XTpw4ocOHD/uNKS4urvY9K6c4HA45nU6/BQAABK86DyhlZWXauXOnEhMT1aVLFzVq1Ej5+fm+7du2bdPu3bvldrvruhQAANBA1Polnv/4j//QwIEDlZKSov3792vKlCkKCwvT/fffr5iYGA0fPlw5OTmKi4uT0+nUo48+KrfbzR08AADAp9YDyt69e3X//ffrp59+UvPmzXXzzTfr888/V/PmzSVJL7/8skJDQzVo0CCVl5crIyNDb7zxRm2XAQAAGrAQY4yp7yIC5fV6FRMTo9LSUt6PAjQwV078qL5LQAC+fz6zvktAEAnk9zffxQMAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOuH1XQCAhoFvIQZwMTGDAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADW4aPugSDHR9QDaIiYQQEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArMMHtQEAzqi2Pujv++cza+U4uHQwgwIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDrcxQMAqHPcDYRAMYMCAACsQ0ABAADWIaAAAADr8B4UoJZxrR0ALhwzKAAAwDrMoACWqq2ZGABoiAgoaPC4pAIAwYdLPAAAwDrMoAAAGgzbLn0y81p36jWgzJw5UzNmzJDH41HHjh312muvqVu3bvVZUlDjUsjZ2fYPHwD78e9q3am3Szzz589XTk6OpkyZoi+//FIdO3ZURkaGSkpK6qskAABgiRBjjKmPJ05PT1fXrl31+uuvS5KqqqqUnJysRx99VBMnTjzrvl6vVzExMSotLZXT6bwY5darYP2ffW39jyFYfz4AUJ/qYlYnkN/f9XKJ58SJEyoqKlJubq5vXWhoqHr37q3CwsLTxpeXl6u8vNz3uLS0VNI/G60L101ZVifHhb+W4xbWdwkAgDOoi9+xp455PnMj9RJQfvzxR1VWViohIcFvfUJCgr799tvTxk+bNk1Tp049bX1ycnKd1QgAwKUs5g91d+wjR44oJibmrGMaxF08ubm5ysnJ8T2uqqrSoUOH1KxZM4WEhMjr9So5OVl79uy5JC75/Ct6p/dLqfdLtW+J3uk9OHo3xujIkSNKSko659h6CSiXX365wsLCVFxc7Le+uLhYLpfrtPEOh0MOh8NvXWxs7GnjnE5nUJzAmqB3er+UXKp9S/RO7w3fuWZOTqmXu3giIiLUpUsX5efn+9ZVVVUpPz9fbre7PkoCAAAWqbdLPDk5ORo6dKhuuOEGdevWTX/4wx909OhRPfzww/VVEgAAsES9BZT77rtPBw8e1OTJk+XxeNSpUyctXbr0tDfOng+Hw6EpU6acdhnoUkDv9H4puVT7luid3i+93uvtc1AAAADOhC8LBAAA1iGgAAAA6xBQAACAdQgoAADAOg0ioOzbt08PPPCAmjVrpqioKLVv314bNmw46z6rVq3S9ddfL4fDoauuukpz5sy5OMXWskB7X7VqlUJCQk5bPB7PRaz6wl155ZXV9pGdnX3GfRYuXKjU1FRFRkaqffv2+vjjjy9ixbUn0N7nzJlz2tjIyMiLXPWFq6ys1JNPPqlWrVopKipKbdq00TPPPHPO7+wIhtd6TXoPlte69M+PPR87dqxSUlIUFRWlG2+8UV988cVZ9wmG8y4F3nswnfdzMpY7dOiQSUlJMQ899JBZt26d+e6778yyZcvMjh07zrjPd999Z6Kjo01OTo7ZunWree2110xYWJhZunTpRaz8wtWk95UrVxpJZtu2bebAgQO+pbKy8iJWfuFKSkr86l++fLmRZFauXFnt+LVr15qwsDAzffp0s3XrVjNp0iTTqFEjs3nz5otbeC0ItPfZs2cbp9Ppt4/H47m4RdeCZ5991jRr1swsWbLE7Nq1yyxcuNA0adLEvPLKK2fcJ1he6zXpPVhe68YYc++995q0tDRTUFBgtm/fbqZMmWKcTqfZu3dvteOD5bwbE3jvwXTez8X6gDJhwgRz8803B7TP448/bq699lq/dffdd5/JyMiozdLqXE16P/WX9+eff66bourJmDFjTJs2bUxVVVW12++9916TmZnpty49Pd38+7//+8Uor06dq/fZs2ebmJiYi1tUHcjMzDTDhg3zW3fXXXeZrKysM+4TLK/1mvQeLK/1Y8eOmbCwMLNkyRK/9ddff7154oknqt0nWM57TXoPlvN+Pqy/xPPhhx/qhhtu0D333KP4+Hh17txZf/7zn8+6T2FhoXr37u23LiMjQ4WFhXVZaq2rSe+ndOrUSYmJierTp4/Wrl1bx5XWrRMnTujdd9/VsGHDFBISUu2YYDnnv3Q+vUtSWVmZUlJSlJycrDvuuENbtmy5iFXWjhtvvFH5+fn6xz/+IUn6+9//rjVr1qh///5n3CdYzntNej+lob/WKyoqVFlZedplyaioKK1Zs6bafYLlvNek91Ma+nk/L/WdkM7F4XAYh8NhcnNzzZdffmn++Mc/msjISDNnzpwz7nP11Veb5557zm/dRx99ZCSZY8eO1XXJtaYmvX/77bdm1qxZZsOGDWbt2rXm4YcfNuHh4aaoqOgiVl675s+fb8LCwsy+ffvOOKZRo0YmLy/Pb93MmTNNfHx8XZdXp86n988++8y89dZb5quvvjKrVq0yt99+u3E6nWbPnj0XsdILV1lZaSZMmGBCQkJMeHi4CQkJOe11/EvB8lqvSe/B9Fp3u93m1ltvNfv27TMVFRXmnXfeMaGhoaZt27bVjg+W825M4L0H03k/F+sDSqNGjYzb7fZb9+ijj5ru3bufcZ9g+ctbk96r06NHD/PAAw/UZmkXVd++fc3tt99+1jHBGlDOp/dfOnHihGnTpo2ZNGlSHVVVN9577z3TokUL895775lNmzaZt99+28TFxV0S/xmpSe/Vaaiv9R07dpgePXoYSSYsLMx07drVZGVlmdTU1GrHB8t5Nybw3qvTUM/7uVh/iScxMVFpaWl+69q1a6fdu3efcR+Xy6Xi4mK/dcXFxXI6nYqKiqqTOutCTXqvTrdu3bRjx47aLO2i+eGHH/Tpp5/qt7/97VnHnemcu1yuuiyvTp1v77/UqFEjde7cucGd8/Hjx2vixIkaPHiw2rdvryFDhmjcuHGaNm3aGfcJltd6TXqvTkN9rbdp00YFBQUqKyvTnj17tH79ep08eVKtW7eudnywnHcp8N6r01DP+7lYH1Buuukmbdu2zW/dP/7xD6WkpJxxH7fbrfz8fL91y5cvl9vtrpMa60pNeq/Oxo0blZiYWJulXTSzZ89WfHy8MjMzzzouWM75vzrf3n+psrJSmzdvbnDn/NixYwoN9f8nKSwsTFVVVWfcJ1jOe016r05Dfq1LUuPGjZWYmKiff/5Zy5Yt0x133FHtuGA57//qfHuvTkM/72dU31M457J+/XoTHh5unn32WbN9+3Yzd+5cEx0dbd59913fmIkTJ5ohQ4b4Hp+6BW38+PHmm2++MTNnzmyQt6DVpPeXX37ZLFq0yGzfvt1s3rzZjBkzxoSGhppPP/20Plq4IJWVlaZly5ZmwoQJp20bMmSImThxou/x2rVrTXh4uHnxxRfNN998Y6ZMmdJgbzM2JrDep06dapYtW2Z27txpioqKzODBg01kZKTZsmXLxSz5gg0dOtRcccUVvltt33//fXP55Zebxx9/3DcmWF/rNek9mF7rS5cuNZ988on57rvvzN/+9jfTsWNHk56ebk6cOGGMCd7zbkzgvQfTeT8X6wOKMcYsXrzYXHfddcbhcJjU1FTzpz/9yW/70KFDza233uq3buXKlaZTp04mIiLCtG7d2syePfviFVyLAu39hRdeMG3atDGRkZEmLi7O9OzZ06xYseIiV107li1b5rvf/5duvfVWM3ToUL91CxYsMG3btjURERHm2muvNR999NFFqrT2BdL72LFjTcuWLU1ERIRJSEgwAwYMMF9++eVFrLZ2eL1eM2bMGNOyZUsTGRlpWrdubZ544glTXl7uGxOsr/Wa9B5Mr/X58+eb1q1bm4iICONyuUx2drY5fPiwb3uwnndjAu89mM77uYQYc46PaQQAALjIrH8PCgAAuPQQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQDUidWrV2vgwIFKSkpSSEiIFi1aVKfPd+TIEY0dO1YpKSmKiorSjTfeqC+++OKCjvm73/1OXbp0kcPhUKdOnc5rn507d+o3v/mNmjdvLqfTqXvvvfe07w768ssv1adPH8XGxqpZs2YaOXKkysrK/Mbk5+frxhtvVNOmTeVyuTRhwgRVVFT4jVmwYIE6deqk6OhopaSkaMaMGafVM3fuXHXs2FHR0dFKTEzUsGHD9NNPPwX2gwjA+++/rz59+vj6d7vdWrZsWcDHIaAAAOrE0aNH1bFjR82cOfOiPN9vf/tbLV++XO+88442b96svn37qnfv3tq3b98Z97nyyiu1atWqsx532LBhuu+++86rhqNHj6pv374KCQnRihUrtHbtWp04cUIDBw70fbfS/v371bt3b1111VVat26dli5dqi1btuihhx7yHefvf/+7BgwYoH79+umrr77S/Pnz9eGHH2rixIm+MZ988omysrI0atQoff3113rjjTf08ssv6/XXX/eNWbt2rR588EENHz5cW7Zs0cKFC7V+/XqNGDHivPqpidWrV6tPnz76+OOPVVRUpNtuu00DBw7UV199FdiB6vujbAEAwU+S+eCDD/zWHT9+3Dz22GMmKSnJREdHm27dupmVK1fW6PjHjh0zYWFhZsmSJX7rr7/+evPEE0+ccb+UlJTzes4pU6aYjh07nnPcsmXLTGhoqCktLfWtO3z4sAkJCTHLly83xhjzxz/+0cTHx5vKykrfmE2bNhlJZvv27cYYY3Jzc80NN9zgd+wPP/zQREZGGq/Xa4wx5v777zd3332335hXX33VtGjRwlRVVRljjJkxY4Zp3br1aWOuuOIKv3V//vOfTWpqqnE4HOaaa64xM2fOPGevgUhLSzNTp04NaB9mUAAA9WL06NEqLCzUvHnztGnTJt1zzz3q16+ftm/fHvCxKioqVFlZqcjISL/1UVFRWrNmTW2VfE7l5eUKCQmRw+HwrYuMjFRoaKivjvLyckVERPh9g3VUVJQk+Y2prpfjx4+rqKjorGP27t2rH374QdI/v/l5z549+vjjj2WMUXFxsf7yl79owIABvn3mzp2ryZMn69lnn9U333yj5557Tk8++aTeeuutWvmZVFVV6ciRI4qLiwtsx1qNSAAAVEO/mEH54YcfTFhYmNm3b5/fuF69epnc3NwaPYfb7Ta33nqr2bdvn6moqDDvvPOOCQ0NNW3btj3jPrU9g1JSUmKcTqcZM2aMOXr0qCkrKzOjR482kszIkSONMcZ8/fXXJjw83EyfPt2Ul5ebQ4cOmUGDBhlJ5rnnnjPG/P9MTF5enqmoqDB79+41t9xyi5Fk8vLyjDH/nImJjo42n376qamsrDTbtm0zqampRpL57LPPfDUtWLDANGnSxISHhxtJZuDAgb5vSzbGmDZt2viOecozzzxj3G73Ofs9Hy+88IK57LLLTHFxcUD7MYMCALjoNm/erMrKSrVt21ZNmjTxLQUFBdq5c6ck6dtvv1VISMhZl399T8Y777wjY4yuuOIKORwOvfrqq7r//vv9ZipGjRrl93y7d+9W//79/dZdiObNm2vhwoVavHixmjRpopiYGB0+fFjXX3+9r45rr71Wb731ll566SVFR0fL5XKpVatWSkhI8I3p27evZsyYoVGjRsnhcKht27a+WY9TY0aMGKHRo0fr9ttvV0REhLp3767Bgwf7jdm6davGjBmjyZMnq6ioSEuXLtX333+vUaNGSfrne2Z27typ4cOH+/0Mfv/73/vOgyS5XK6znofu3btX+/PIy8vT1KlTtWDBAsXHxwf0s+TbjAEAdS4kJEQffPCB7rzzTknS/PnzlZWVpS1btigsLMxvbJMmTeRyuXTixAl99913Zz1us2bN1Lx5c791R48eldfrVWJiou677z6VlZXpo48+kiSVlJTI6/X6xvbs2VMvvPCC0tPTfeuuuuqq057nqaee0qJFi7Rx48bz7vnHH39UeHi4YmNj5XK59Nhjj2n8+PF+Y4qLi9W4cWOFhITI6XRq3rx5uueee3zbjTE6cOCALrvsMn3//fdKS0vT+vXr1bVrV9+YyspKeTweNW/eXPn5+RowYIBKSkrUvHlzDRkyRMePH9fChQt949esWaNbbrlF+/fvV2hoqFwul959912/n4EkhYWFqVWrVpKkHTt2nHYH0b+KiopSSkqK37p58+Zp2LBhWrhwoTIzM8/753ZKeMB7AABwgTp37qzKykqVlJTolltuqXZMRESEUlNTAz5248aN1bhxY/38889atmyZpk+f7tsWHx/v9z/58PBwXXHFFdWGkgt1+eWXS5JWrFihkpIS/frXvz5tTEJCgiTpv//7vxUZGak+ffr4bQ8JCVFSUpIk6b333lNycrKuv/56vzFhYWG64oorfGPcbrcvtB07dkzh4eGnjZf+GX4SEhKUlJSk7777TllZWWfsJdCfz3vvvadhw4Zp3rx5NQonEgEFAFBHysrKtGPHDt/jXbt2aePGjYqLi1Pbtm2VlZWlBx98UC+99JI6d+6sgwcPKj8/Xx06dKjRL7Vly5bJGKNrrrlGO3bs0Pjx45WamqqHH364xj3s2LFDZWVl8ng8+t///V/fDEpaWpoiIiK0b98+9erVS2+//ba6desmSZo9e7batWun5s2bq7CwUGPGjNG4ceN0zTXX+I77+uuv68Ybb1STJk20fPlyjR8/Xs8//7xiY2N9Y2bMmKF+/fopNDRU77//vp5//nktWLDAFzB+/PFH/eUvf1HPnj11/PhxzZ49WwsXLlRBQYHvGAMHDtSIESP05ptvKiMjQwcOHNDYsWPVrVs3X/CZOnWqfve73ykmJkb9+vVTeXm5NmzYoJ9//lk5OTkB/8zy8vI0dOhQvfLKK0pPT5fH45H0z1mWmJiY8z9QrbwDBgCAX1i5cqWRdNoydOhQY4wxJ06cMJMnTzZXXnmladSokUlMTDS/+c1vzKZNm2r0fPPnzzetW7c2ERERxuVymezsbHP48OGz7nOuN8neeuut1fawa9cuY4wxu3btMpL8jjFhwgSTkJBgGjVqZK6++mrz0ksv+W77PWXIkCEmLi7OREREmA4dOpi33377tOe+7bbbTExMjImMjDTp6enm448/9tt+8OBB0717d9O4cWMTHR1tevXqZT7//PPTjvPqq6+atLQ0ExUVZRITE01WVpbZu3ev35i5c+eaTp06mYiICHPZZZeZHj16mPfff/+sP7tAf2anzvv54j0oAADAOtzFAwAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1/g/G/1BVir1CUQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.hist(g_best_samples[:, 0], bins=25)" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(Array(0.01999981, dtype=float64), Array(0.01998104, dtype=float64))" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_best_samples[:, 0].mean(), g_pos_samples[:, 0].mean()" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(Array(3.50319931e-10, dtype=float64), Array(3.3719104e-05, dtype=float64))" + ] + }, + "execution_count": 44, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_best_samples[:, 0].std(), g_pos_samples[:, 0].std()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Jackknife for error on the mean" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# let's start by splitting into 10 batches\n", + "\n", + "# we only need to repeat the second step of the inference" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(1000, 50, 2)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "e_post_pos.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0\n", + "1\n", + "2\n", + "3\n", + "4\n", + "5\n", + "6\n", + "7\n", + "8\n", + "9\n" + ] + } + ], + "source": [ + "from math import ceil\n", + "n_batches = 10\n", + "batch_size = ceil(n_gals / n_batches)\n", + "\n", + "# g_pos_list = [] \n", + "# g_neg_list = []\n", + "g_best_list = [] \n", + "\n", + "keys = random.split(k2, n_batches)\n", + "\n", + "for ii in range(n_batches): \n", + " print(ii)\n", + "\n", + " k_ii = keys[ii]\n", + " start, end = ii * batch_size, (ii+1) * batch_size\n", + " e_pos1 = jnp.concatenate([e_post_pos[:start], e_post_pos[end:]], axis=0)\n", + " e_neg1 = jnp.concatenate([e_post_neg[:start], e_post_neg[end:]], axis=0)\n", + "\n", + " g_pos1 = pipeline_shear_inference_ellipticities(k_ii, e_post_pos, true_g, sigma_e=sigma_e, sigma_e_int=sigma_e_int, n_samples=1000, initial_step_size=1e-3)\n", + "\n", + " g_neg1 = pipeline_shear_inference_ellipticities(k_ii, e_post_neg, -true_g, sigma_e=sigma_e, sigma_e_int=sigma_e_int, n_samples=1000, initial_step_size=1e-3)\n", + "\n", + " g_best1 = (g_pos1 - g_neg1) * 0.5\n", + " g_best_mean1 = g_best1.mean(axis=0)\n", + "\n", + " g_best_list.append(g_best_mean1)\n", + "\n", + " # g_pos_list.append(g_pos1)\n", + " # g_neg_list.append(g_neg1)\n", + "\n", + "# g_pos_jack = jnp.stack(g_pos_list, axis=0)\n", + "# g_neg_jack = jnp.stack(g_neg_list, axis=0)\n", + "g_best_means = jnp.array(g_best_list)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(10, 2)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_best_means.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Array([0.01999975, 0.01999981, 0.01999981, 0.01999938, 0.01999988,\n", + " 0.01999977, 0.01999979, 0.01999971, 0.01999985, 0.01999981], dtype=float64)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_best_means[:, 0]" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Array(0.01999976, dtype=float64)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_best_means[:, 0].mean() #g1" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Array(3.96908294e-07, dtype=float64)" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# error on the mean (jackknife expression)\n", + "jnp.sqrt(g_best_means[:, 0].var() * (n_batches-1))" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Array(-2.44035424e-07, dtype=float64)" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_best_means[:, 0].mean() - 0.02" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "101" + ] + }, + "execution_count": 58, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g_" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(900, 50, 2)" + ] + }, + "execution_count": 52, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "jnp.concatenate([e_post_pos[:100], e_post_neg[200:]], axis=0).shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "bpd_gpu3", + "language": "python", + "name": "bpd_gpu3" + }, + "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.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}