diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..f6ba9ce Binary files /dev/null and b/.DS_Store differ diff --git a/README.md b/README.md index 677d7f5..f1caf15 100644 --- a/README.md +++ b/README.md @@ -27,9 +27,9 @@ $ conda env create -f environment.yml $ conda activate shear ``` -## Install the `shrbk` library +## Install the `shrbk` library -`shrbk` is the library made for this book. +`shrbk` is the library made for this book. Install it ither in development mode ```bash diff --git a/notebooks/taylor_series.ipynb b/notebooks/taylor_series.ipynb new file mode 100755 index 0000000..8f9a155 --- /dev/null +++ b/notebooks/taylor_series.ipynb @@ -0,0 +1,219 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Taylor Series for Approximations\n", + "\n", + "Taylor series are commonly used in physics to approximate functions making them easier to handle specially when solving equations. In this notebook we give a visual example on how it works and the biases that it introduces.\n", + "\n", + "## Theoretical Formula\n", + "\n", + "Consider a function $f$ that is $n$ times differentiable in a point $a$. Then by Taylor's theorem, for any point $x$ in the domain of f, we have the Taylor expansion about the point $a$ is defined as:\n", + "\\begin{equation}\n", + "f(x) = f(a) + \\sum_{k=1}^n \\frac{f^{k}(a)}{k!}(x-a)^k + o\\left((x-a)^n\\right) \\quad,\n", + "\\end{equation}\n", + "where $f^{(k)}$ is the derivative of order $k$ of $f$. Usually, we consider $a=0$ which gives:\n", + "\\begin{equation}\n", + "f(x) = f(0) + \\sum_{k=1}^n \\frac{f^{k}(0)}{k!}(x)^k + o\\left((x)^n\\right) \\quad.\n", + "\\end{equation}\n", + "\n", + "For example, the exponential, $e$ is infinitely differentiable with $e^{(k)}=e$ and $e^0=1$. This gives us the following Taylor expansion:\n", + "\\begin{equation}\n", + "e(x) = 1 + \\sum_{k=1}^\\infty \\frac{x^k}{k!} \\quad.\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualising Taylor Expansion Approximation and its Bias\n", + "\n", + "Let us see visually how the Taylor expansion approximatees a given function. We start by defining our function below, for example we will consider the exponential function, $e$ again up to order 3." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "code_folding": [ + 0 + ] + }, + "outputs": [], + "source": [ + "#### FOLDED CELL\n", + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "from IPython.display import Markdown as md\n", + "from sympy import Symbol, series, lambdify, latex\n", + "from sympy.functions import *\n", + "from ipywidgets import interactive_output\n", + "import ipywidgets as widgets\n", + "from sympy.parsing.sympy_parser import parse_expr\n", + "import numpy as np\n", + "\n", + "x = Symbol('x')" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "order = 3\n", + "func = exp(x)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "code_folding": [ + 0 + ] + }, + "outputs": [], + "source": [ + "#### FOLDED CELL\n", + "taylor_exp = series(func,x,n=order+1)\n", + "approx = lambdify(x, sum(taylor_exp.args[:-1]), \"numpy\")\n", + "func_np = lambdify(x, func, \"numpy\")\n", + "latex_func = '$'+latex(func)+'$'\n", + "latex_taylor = '\\\\begin{equation} '+latex(taylor_exp)+' \\end{equation}'" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "variables": { + " latex_func ": "$e^{x}$", + "latex_taylor": "\\begin{equation} 1 + x + \\frac{x^{2}}{2} + \\frac{x^{3}}{6} + O\\left(x^{4}\\right) \\end{equation}" + } + }, + "source": [ + "The Taylor expansion of {{ latex_func }} is :\n", + "{{latex_taylor}}\n", + "\n", + "Now let's plot the function and its expansion while considering a point, noted $p$, to study the biais that we introduce when we approximate the function by its expansion:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "code_folding": [ + 0 + ] + }, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "5f5d027e2c284af4869d9e9c6ea083f9", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output(layout=Layout(height='650px'))" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d8b970e851ef4e688880bc1db4a46238", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(FloatSlider(value=3.0, description='Point Absciss', max=4.0, min=-4.0, step=0.2), IntSlider(val…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#### FOLDED CELL\n", + "order = widgets.IntSlider(min=0,max=20,step=1,value=3,description='Order')\n", + "x_min = -4\n", + "x_max = 4\n", + "x1 = widgets.FloatSlider(min=x_min,max=x_max,value=3,step=0.2,description='Point Absciss')\n", + "func = widgets.Text('exp(x)',description='Function')\n", + "text_offset = np.array([-0.15,2.])\n", + "\n", + "\n", + "ui = widgets.HBox([x1, order, func])\n", + "\n", + "def f(order=widgets.IntSlider(min=1,max=10,step=1,value=3)\n", + " ,x1=1.5\n", + " ,func='exp(x)'):\n", + " func_sp = parse_expr(func)\n", + " taylor_exp = series(func_sp,x,n=order+1)\n", + " approx = lambdify(x, sum(taylor_exp.args[:-1]), \"numpy\")\n", + " func_np = lambdify(x, func_sp, \"numpy\")\n", + " n_points = 1000\n", + " x_array = np.linspace(x_min,x_max,n_points)\n", + " approx_array = np.array([approx(z) for z in x_array])\n", + " func_array = np.array([func_np(z) for z in x_array])\n", + " func_x1 = func_np(x1)\n", + " approx_x1 = approx(x1)\n", + " plt.figure(42,figsize=(10,10))\n", + " plt.plot(x_array,approx_array,color='blue',label='Taylor Expansion')\n", + " plt.plot(x_array,func_array,color='green',label=func)\n", + " plt.plot(0,approx(0),color='black',marker='o')\n", + " plt.annotate(r'(0,0)',[0,approx(0)],xytext=text_offset)\n", + " plt.plot([x1,x1]\n", + " ,[-np.max(np.abs([np.min(func_array),np.max(func_array)])),min(approx_x1, func_x1)]\n", + " ,'--',color='black',marker='x')\n", + " plt.plot([x1,x1],[approx_x1, func_x1],'r--',marker='x')\n", + " plt.annotate(r'$p_{approx}$',[x1,approx(x1)],xytext=[x1,approx(x1)]-text_offset)\n", + " plt.annotate(r'$p$',[x1,func_np(x1)],xytext=[x1,func_np(x1)]-text_offset)\n", + " plt.xlim([x_min,x_max])\n", + " plt.ylim(-np.max(np.abs([np.min(func_array),np.max(func_array)]))\n", + " ,np.max(np.abs([np.min(func_array),np.max(func_array)])))\n", + " plt.legend()\n", + " plt.show()\n", + " print('Approximation bias : {}'.format(func_x1-approx_x1))\n", + " return None\n", + "interactive_plot = widgets.interactive_output(f, {'order': order, 'x1': x1, 'func': func})\n", + "interactive_plot.layout.height = '650px'\n", + "display(interactive_plot, ui)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that the further $p$ gets away from the point of the expansion (in that case $0$), the higher the approximation bias gets. Samely, the lower the order of approximation is, the higher the approximation bias gets." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.8.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/shearbook/_bibliography/z_moments_methods.bib b/shearbook/_bibliography/z_moments_methods.bib index 25d362d..98e9ec2 100644 --- a/shearbook/_bibliography/z_moments_methods.bib +++ b/shearbook/_bibliography/z_moments_methods.bib @@ -18,3 +18,38 @@ @ARTICLE{2011MNRAS.410.2156V adsurl = {https://ui.adsabs.harvard.edu/abs/2011MNRAS.410.2156V}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } + +@ARTICLE{1995ApJ.KSB, + author = {{Kaiser}, Nick and {Squires}, Gordon and {Broadhurst}, Tom}, + title = "{A Method for Weak Lensing Observations}", + journal = {\apj}, + keywords = {COSMOLOGY: OBSERVATIONS, COSMOLOGY: DARK MATTER, GALAXIES: FORMATION, COSMOLOGY: GRAVITATIONAL LENSING, COSMOLOGY: LARGE-SCALE STRUCTURE OF UNIVERSE, Astrophysics}, + year = 1995, + month = aug, + volume = {449}, + pages = {460}, + doi = {10.1086/176071}, +archivePrefix = {arXiv}, + eprint = {astro-ph/9411005}, + primaryClass = {astro-ph}, + adsurl = {https://ui.adsabs.harvard.edu/abs/1995ApJ...449..460K}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@ARTICLE{2003MNRAS.343..459H, + author = {{Hirata}, Christopher and {Seljak}, Uro{\v{s}}}, + title = "{Shear calibration biases in weak-lensing surveys}", + journal = {\mnras}, + keywords = {gravitational lensing, methods: data analysis, Astrophysics}, + year = 2003, + month = aug, + volume = {343}, + number = {2}, + pages = {459-480}, + doi = {10.1046/j.1365-8711.2003.06683.x}, +archivePrefix = {arXiv}, + eprint = {astro-ph/0301054}, + primaryClass = {astro-ph}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2003MNRAS.343..459H}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} \ No newline at end of file diff --git a/shearbook/_toc.yml b/shearbook/_toc.yml index aa5aba1..90d4cc1 100644 --- a/shearbook/_toc.yml +++ b/shearbook/_toc.yml @@ -5,6 +5,7 @@ chapters: - file: intro/about - file: intro/run_code + - file: intro/weak_lensing - part: Shape Measurement chapters: @@ -17,8 +18,6 @@ - file: methods/moments-methods sections: - file: methods/ksb - - file: methods/ksb-variants - - file: methods/viola2011 - file: methods/moments-methods-ref - file: model/model-intro sections: @@ -35,6 +34,10 @@ sections: - file: calibration/calibration-ref +- part: Appendix + chapters: + - file: appendix/taylor_series + - part: Contributing chapters: - file: contrib/contributors diff --git a/shearbook/appendix/.DS_Store b/shearbook/appendix/.DS_Store new file mode 100644 index 0000000..dc020f1 Binary files /dev/null and b/shearbook/appendix/.DS_Store differ diff --git a/shearbook/appendix/taylor_series.ipynb b/shearbook/appendix/taylor_series.ipynb new file mode 100644 index 0000000..751775f --- /dev/null +++ b/shearbook/appendix/taylor_series.ipynb @@ -0,0 +1,201 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Taylor Series for Approximations\n", + "\n", + "Taylor series are commonly used in physics to approximate functions making them easier to handle specially when solving equations. In this notebook we give a visual example on how it works and the biases that it introduces.\n", + "\n", + "## Theoretical Formula\n", + "\n", + "Consider a function $f$ that is $n$ times differentiable in a point $a$. Then by Taylor's theorem, for any point $x$ in the domain of f, we have the Taylor expansion about the point $a$ is defined as:\n", + "\\begin{equation}\n", + "f(x) = f(a) + \\sum_{k=1}^n \\frac{f^{k}(a)}{k!}(x-a)^k + o\\left((x-a)^n\\right) \\quad,\n", + "\\end{equation}\n", + "where $f^{(k)}$ is the derivative of order $k$ of $f$. Usually, we consider $a=0$ which gives:\n", + "\\begin{equation}\n", + "f(x) = f(0) + \\sum_{k=1}^n \\frac{f^{k}(0)}{k!}(x)^k + o\\left((x)^n\\right) \\quad.\n", + "\\end{equation}\n", + "\n", + "For example, the exponential, $e$ is infinitely differentiable with $e^{(k)}=e$ and $e^0=1$. This gives us the following Taylor expansion:\n", + "\\begin{equation}\n", + "e(x) = 1 + \\sum_{k=1}^\\infty \\frac{x^k}{k!} \\quad.\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualising Taylor Expansion Approximation and its Bias\n", + "\n", + "Let us see visually how the Taylor expansion approximatees a given function. We start by defining our function below, for example we will consider the exponential function, $e$ again up to order 3." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "code_folding": [ + 0 + ] + }, + "outputs": [], + "source": [ + "#### FOLDED CELL\n", + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "from IPython.display import Markdown as md\n", + "from sympy import Symbol, series, lambdify, latex\n", + "from sympy.functions import *\n", + "from ipywidgets import interactive_output\n", + "import ipywidgets as widgets\n", + "from sympy.parsing.sympy_parser import parse_expr\n", + "import numpy as np\n", + "\n", + "x = Symbol('x')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "order = 3\n", + "func = exp(x)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "code_folding": [ + 0 + ] + }, + "outputs": [], + "source": [ + "#### FOLDED CELL\n", + "taylor_exp = series(func,x,n=order+1)\n", + "approx = lambdify(x, sum(taylor_exp.args[:-1]), \"numpy\")\n", + "func_np = lambdify(x, func, \"numpy\")\n", + "latex_func = '$'+latex(func)+'$'\n", + "latex_taylor = '\\\\begin{equation} '+latex(taylor_exp)+' \\end{equation}'" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "variables": { + " latex_func ": "$e^{x}$", + "latex_taylor": "\\begin{equation} 1 + x + \\frac{x^{2}}{2} + \\frac{x^{3}}{6} + O\\left(x^{4}\\right) \\end{equation}" + } + }, + "source": [ + "The Taylor expansion of {{ latex_func }} is :\n", + "{{latex_taylor}}\n", + "\n", + "Now let's plot the function and its expansion while considering a point, noted $p$, to study the biais that we introduce when we approximate the function by its expansion:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "code_folding": [ + 0 + ] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAI/CAYAAABAoBw9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABQt0lEQVR4nO3dd3gU1eLG8fekQEioQuhdiiAlaEBFVKqNaqOrXFSs94rlKiBeflZQsSt2RQEvIDaaIiCogCIIRDrSewkEEkhPzu+PAYRLWALZyewm38/z7LPZ3dmZN4vK6zmzZ4y1VgAAAHBHiNcBAAAACjLKFgAAgIsoWwAAAC6ibAEAALiIsgUAAOAiyhYAAICLwrwOcKJy5crZmjVreh0DAADApx2JO7T7r93x1troM20bUGWrZs2aWrx4sdcxAAAAfGrxQQvtHrB7S262ZRoRAADgLBxMPag/dv2R6+0pWwAAAGdh7ua5yrbZud6esgUAAHAWZm6YqajwqFxvH1DnbOUkIyND27dvV2pqqtdR4GcRERGqWrWqwsPDvY4CAECuzdo0S1fVvErTNT1X2wd82dq+fbtKlCihmjVryhjjdRz4ibVW+/fv1/bt21WrVi2v4wAAkCtbD23Vuv3rdM/F9+S6bAX8NGJqaqrKli1L0SpgjDEqW7YsI5YAgKAye+NsSVL72u1z/Z6AL1uSKFoFFH+uAIBgM2vTLJWPKq9G5Rvl+j1BUba8tH//fsXExCgmJkYVK1ZUlSpVjj9OT0/P1T769eunSZMmnXOG//u//zvpuDExMTp48OA57y+vJk+erBEjRnh2fAAAvGCt1ayNs9S+dvuzGjAI+HO2vFa2bFktW7ZMklN6ihcvrkcffdTVY2ZlZSk0NPSk5x566CHXj5tbXbp0UZcuXbyOAQBAvlqxd4X2HtmrDrU7nNX7GNk6Bx988IGaN2+upk2b6qabblJycrKSkpJUq1YtZWRkSJISExNVs2bN44+PmT17tpo1a6bGjRurf//+SktLk+Ssnv/000+rVatW+uKLL3KV45VXXlH//v0lScuXL1ejRo2UnJys//u//9Ott96qtm3bqm7duvrggw8kSYcPH1a7du100UUXqXHjxvr2228lSZs3b1aDBg1011136cILL9TVV1+tlJQUSdIbb7yhhg0bqkmTJurZs6ckafTo0XrggQckSVu2bFG7du3UpEkTtWvXTlu3bpXkjOb961//UsuWLVW7du08jewBABAIZm2cJUlqV6vdWb2PsnUObrzxRi1atEhxcXFq0KCBPvroI5UoUUKtW7fWtGnTJEnjx4/XTTfddNKyBqmpqerXr58mTJig5cuXKzMzU++8887x1yMiIjRv3rzjpeZEr7766vEpxDZt2kiSBg4cqPXr1+vrr7/WP/7xD7333nuKjIyUJP3555+aNm2afv31Vz399NPauXOnIiIi9PXXX2vJkiWaM2eOHnnkEVlrJUl//fWX7r//fq1cuVKlS5fWl19+KUkaMWKEli5dqj///FPvvvvuKbkeeOAB3Xbbbfrzzz/Vp08f/etf/zr+2q5duzRv3jxNnTpVgwYNyuvHDgCAp2ZtmqX6ZeurWqlqZ/W+oJpGHDhQOjqj5zcxMdJrr53de1asWKGhQ4fq4MGDOnz4sK655hpJ0p133qkXX3xR3bp10yeffHJ8ROmYtWvXqlatWqpXr54k6fbbb9fbb7+tgQMHSpJ69Ohx2mPmNI0YEhKi0aNHq0mTJrr77rt1+eWXH3+ta9euKlasmIoVK6Y2bdro999/V8eOHTVkyBD9/PPPCgkJ0Y4dO7Rnzx5JUq1atRQTEyNJuvjii7V582ZJUpMmTdSnTx9169ZN3bp1OyXXr7/+qq+++kqSdOutt+qxxx47/lq3bt0UEhKihg0bHj8OAADBKD0rXT9t/kn9Yvqd9XsZ2ToH/fr101tvvaXly5dr2LBhx5cvuPzyy7V582b99NNPysrKUqNGJ39T4dgo0ulEReV+Ndpj/vrrLxUvXlw7d+486fn/PXHPGKNx48Zp3759+uOPP7Rs2TJVqFDhePaiRYse3zY0NFSZmZmSpGnTpun+++/XH3/8oYsvvvj486dz4nFP3OeZfncAAALZwu0LdSTjyFkt+XBMUI1sne0IlFuSkpJUqVIlZWRkaNy4capSpcrx12677Tb16tVLTz755Cnvu+CCC7R582atX79ederU0ZgxY3TVVVedc45Dhw7pwQcf1M8//6wHHnhAkyZN0s033yxJ+vbbbzV48GAdOXJEc+fO1YgRI/TFF1+ofPnyCg8P15w5c7Rli++LlWdnZ2vbtm1q06aNWrVqpc8//1yHDx8+aZuWLVtq/PjxuvXWWzVu3Di1atXqnH8fAAAC1ayNsxRiQtS6Zuuzfm9Qla1A8cwzz+iSSy5RjRo11LhxYyUlJR1/rU+fPho6dKh69ep1yvsiIiL0ySef6JZbblFmZqaaN2+ue+65J1fHfPXVVzV27Njjj7/55hs9/fTTuu+++1SvXj199NFHatOmja688kpJUosWLdSxY0dt3bpVTz75pCpXrqw+ffqoc+fOio2NVUxMjC644AKfx8zKylLfvn116NAhWWv10EMPqXTp0idt88Ybb6h///566aWXFB0drU8++SRXvw8AAMFk5saZal65uUpHlD7r95pAmt6JjY21ixcvPum51atXq0GDBh4lOnuTJk3St99+qzFjxniWIb+WqPCHYPvzBQAUPgkpCSr3UjkNaTVEz7R95vjzxpg/rLWxZ3o/I1t+9M9//lPfffedpk/P3bWSAABA4Ju9abaybbaurXPtOb2fsuVHb775ptcRJDkjWwAAwD++X/+9ShUtpUuqXnJO7+fbiAAAAKdhrdX3679Xh/M7KCzk3MaoKFsAAACnsXLfSu1I2qFrzz+3KUSJsgUAAHBa36//XpJ0TZ1rznkflC0AAIDT+H7992pUvpGqlqx6zvugbOUja63atm2rxMTE026zb98+XXvtuQ9VAgAA/zicfli/bP0lT1OIEmUrX02fPl1NmzZVyZIlT7tNdHS0KlWqpPnz5+djMgAA8L/mbp6r9Kz0c17y4RjKVi6NHTtWLVq0UExMjO6++24tXLhQTZo0UWpqqo4cOaILL7xQK1as0Ny5c3XllVfqhhtuUMOGDXXPPfcoOztbkjRu3Dh17dpVkrRo0aIc3y85F3AeN26cZ78rAABwphAjwyPVqnreLkVH2cqF1atXa8KECZo/f76WLVum0NBQrV27Vl26dNHQoUP12GOPqW/fvscvPP3777/r5Zdf1vLly7VhwwZ99dVXkqT58+fr4osvliQ1b978tO+PjY3VL7/84s0vCwAAJDllq22ttioaVjRP+wmqRU0Hfj9Qy3Yv8+s+YyrG6LVrX/O5zezZs/XHH3+oefPmkqSUlBSVL19e//nPf9S8eXNFRETojTfeOL59ixYtVLt2bUlSr169NG/ePN188806cOCASpQocXy7072/fPny2rlzpx9/SwAAcDbWH1ivDQkb9NClD+V5X0FVtrxirdXtt9+u4cOHn/T87t27dfjwYWVkZCg1NVVRUVGSJGPMSdsdexwWFqbs7GyFhDgDigcOHMjx/ampqSpWrJjbvxYAADiNY0s+5PV8LSnIytaZRqDc0q5dO3Xt2lUPPfSQypcvrwMHDigpKUn//Oc/9cwzz2jTpk16/PHH9dZbb0lyphE3bdqkGjVqaMKECRowYIAkqX79+tq4caPq1KkjSRowYECO71+3bt3xKUUAAJD/vl//veqcV0fnn3d+nvcVVGXLKw0bNtSzzz6rq6++WtnZ2QoPD1fXrl0VFham3r17KysrSy1bttSPP/6okJAQXXbZZRo0aJCWL19+/GR5SerYsaPmzp2rOnXq6LPPPsvx/W3bttWcOXPUsWNHj39rAAAKp+SMZM3eNFt3XXSXX/ZnrLV+2ZE/xMbG2sWLF5/03OrVq9WgQQOPEp29uXPnauTIkZo6deopr+3atUu33XabZs6c6XMfV155pb799luVKVPGrZgBI9j+fAEABd/UdVPV+b+d9UPfH9Th/A6n3c4Y84e1NvZM++PbiPmoUqVKuuuuu864qOnDDz9cKIoWAACBaOq6qSpepLiurHGlX/bHNKKftW7dWq1btz7t6927d/f5/ujoaHXr1s2/oQAAQK5YazV13VRdff7VeV7y4RhGtgAAAI6K2xOnHUk71KluJ7/tMyjKViCdVwb/4c8VABBopq6bKiOj6+te77d9BnzZioiI0P79+/mLuYCx1mr//v2KiIjwOgoAAMdNXTdVLaq0UIXiFfy2z4A/Z6tq1aravn279u3b53UU+FlERISqVq3qdQwAACRJew7v0e87ftfTbZ72634DvmyFh4erVq1aXscAAAAF3PS/psvKqlM9/52vJQXBNCIAAEB+mPrXVFUpUUVNKzT1634pWwAAoNBLy0zTDxt+UKd6nU65xnFeUbYAAECh9/OWn3U4/bDfpxAlyhYAAICmrJuiYmHF1K5WO7/vm7IFAAAKNWutJq+drHa126lYeDG/75+yBQAACrVlu5dpy6EtuuGCG1zZP2ULAAAUal+v+VohJkSd63V2Zf+ULQAAUKh9s+YbtareStFR0a7sn7IFAAAKrQ0HNmj53uWuTSFKlC0AAFCIfbPmG0lStwu6uXYMyhYAACi0vl7ztWIqxqhm6ZquHYOyBQAACqU9h/dowbYF6la/m6vHoWwBAIBCafLaybKyuqGBe+drSZQtAABQSH2z9hvVLlNbjcs3dvU4lC0AAFDoJKYlatbGWepWv5vfLzz9vyhbAACg0Pnur++UnpXu+hSiRNkCAACF0NdrvlZ0ZLQuq3qZ68eibAEAgEIlOSNZU9dN1Y0NblRoSKjrx6NsAQCAQuW7v77TkYwj6n5h93w5HmULAAAUKl+s+kLRkdG6ssaV+XI8yhYAACg0UjJSjk8hhoWE5csxKVsAAKDQ+G69M4V4S8Nb8u2YlC0AAFBoTFw5UeUiy+mqmlfl2zEpWwAAoFA4NoV4U4Ob8m0KUaJsAQCAQsKLKUSJsgUAAAqJL1Z9ke9TiBJlCwAAFAIpGSmasnaKbrwg/76FeAxlCwAAFHjHphDzayHTE1G2AABAgTdh5QRFR0bn+xSi5MeyZYwJNcYsNcZMPfr4PGPMTGPMX0fvy/jrWAAAALmVlJakyWsnq/uF3fN9ClHy78jWg5JWn/B4kKTZ1tq6kmYffQwAAJCvvlnzjVIzU9W7cW9Pju+XsmWMqSqpo6QPT3i6q6RPj/78qaRu/jgWAADA2fh8xeeqWbqmLqt6mSfH99fI1muSHpOUfcJzFay1uyTp6H15Px0LAAAgV/Ye2auZG2aqV6NeMsZ4kiHPZcsY00nSXmvtH+f4/gHGmMXGmMX79u3LaxwAAIDjJq6cqCybpT6N+3iWwR8jW5dL6mKM2SxpvKS2xpixkvYYYypJ0tH7vTm92Vr7vrU21lobGx0d7Yc4AAAAjs+Xf64mFZrowvIXepYhz2XLWjvYWlvVWltTUk9JP1pr+0qaLOn2o5vdLunbvB4LAAAgtzYmbNSv239V70benBh/jJvrbI2Q1MEY85ekDkcfAwAA5IvxK8ZLkno26ulpDr8uNmGtnStp7tGf90tq58/9AwAA5Ia1VuOWj1Or6q1Uo3QNT7OwgjwAAChwlu9drlX7Vnk+hShRtgAAQAE0Jm6MwkLCdHPDm72OQtkCAAAFS2Z2psYuH6vr616v6CjvVzqgbAEAgAJl5oaZ2n14t/o17ed1FEmULQAAUMCMjhutssXKqmO9jl5HkUTZAgAABUhCSoK+XfOtejXqpSKhRbyOI4myBQAACpCJKycqLStN/WL6eR3lOMoWAAAoMEbHjdaF0RfqokoXeR3lOMoWAAAoENbGr9Vv23/T7U1vlzHG6zjHUbYAAECB8FncZwoxIerbpK/XUU5C2QIAAEEvKztLn/35ma45/xpVKlHJ6zgnoWwBAICgN2fzHG1P3K7bm97udZRTULYAAEDQ+3DJhyoTUUZdL+jqdZRTULYAAEBQi0+O19drvtatTW5VRFiE13FOQdkCAABB7bO4z5Sela67Lr7L6yg5omwBAICgZa3VB0s+0KVVL1Wj8o28jpMjyhYAAAhaC7Yt0Jr4NbrrosAc1ZIoWwAAIIh9sOQDFS9SXN0v7O51lNOibAEAgKB0MPWgJq6cqN6Neqt4keJexzktyhYAAAhKny//XCmZKQF7YvwxlC0AABCUPlzyoWIqxujiShd7HcUnyhYAAAg6i3cu1tLdS3XXRXcF1EWnc0LZAgAAQWfUolGKCo9Sn8Z9vI5yRpQtAAAQVA6kHNB/V/xXfZv0VamIUl7HOSPKFgAACCqfLP1EqZmpuq/5fV5HyRXKFgAACBrZNlvvLH5Hraq3UpMKTbyOkyuULQAAEDR+2PCDNiRs0P3N7/c6Sq5RtgAAQNB4e9HbqhBVQTc2uNHrKLlG2QIAAEFhU8ImTVs3TXdddJeKhBbxOk6uUbYAAEBQeHfxuwoxIbo79m6vo5wVyhYAAAh4qZmp+mjpR+pSv4uqlqzqdZyzQtkCAAAB7/Pln2t/yv6gOjH+GMoWAAAIaNZavfbba2pSoYna1mrrdZyzFuZ1AAAAAF9+3PSjlu9dro+7fBzw10HMCSNbAAAgoL3626sqH1VevRr38jrKOaFsAQCAgLU2fq2m/TVN98beq4iwCK/jnBPKFgAACFhvLHxDRUKL6N7Ye72Ocs4oWwAAICAdSDmg0XGj1adxH1UoXsHrOOeMsgUAAALSB398oOSMZD14yYNeR8kTyhYAAAg4GVkZemvRW2pbq62aVmzqdZw8YekHAAAQcMavGK/tidv1bsd3vY6SZ4xsAQCAgGKt1YsLXlSj8o10fd3rvY6TZ4xsAQCAgDL9r+lasXeFPuv2WVAuYvq/GNkCAAABZcT8Eapeqrp6NurpdRS/oGwBAICAsWDbAs3bOk+PXPaIwkPDvY7jF5QtAAAQMF6Y/4LKFiurO5rd4XUUv6FsAQCAgLBy70pNXjtZ/2zxT0UVifI6jt9QtgAAQEB4acFLigyP1AMtHvA6il9RtgAAgOe2HNyiccvH6c5md6psZFmv4/gVZQsAAHhu+LzhCjEh+vfl//Y6it9RtgAAgKe2Hdqmj5d+rDua3aGqJat6HcfvKFsAAMBTL8x/QZI0qNUgj5O4g7IFAAA8syNxhz5Y8oH6xfRT9VLVvY7jCsoWAADwzEsLXlJWdpYGtxrsdRTXULYAAIAndh/erff+eE+3Nb1NtcrU8jqOayhbAADAEy/Nf0kZWRl64oonvI7iKsoWAADId7uSdumdxe+oT5M+Ov+8872O4yrKFgAAyHfP/fKcMrIz9J8r/+N1FNdRtgAAQL7afHCz3v/jfd3R7I4CP6olUbYAAEA+e/qnpxViQjT0yqFeR8kXlC0AAJBv1sSv0adxn+r+5vcXyNXic0LZAgAA+WbY3GGKDI8ssKvF54SyBQAA8sXSXUs1ceVEPXTpQ4qOivY6Tr6hbAEAgHwxdM5QlYkoo0cue8TrKPmKsgUAAFz346YfNf2v6RrcarBKRZTyOk6+omwBAABXZdtsPfrDo6pRqob+eck/vY6T78K8DgAAAAq2cX+O09LdSzXuxnGKCIvwOk6+Y2QLAAC4JiUjRU/8+IQurnSxejbq6XUcTzCyBQAAXPP6wte1LXGbPrvhM4WYwjnGUzh/awAA4Lp9R/bp+V+eV5f6XdS6Zmuv43iGsgUAAFzx1E9PKTkjWS+0f8HrKJ6ibAEAAL9bE79G7y5+VwMuHqALyl3gdRxPUbYAAIBfWWs18PuBiioSpWFXDfM6juc4QR4AAPjV5LWTNWPDDL16zauqULyC13E8x8gWAADwm5SMFD004yE1jG6o+5vf73WcgMDIFgAA8JuRC0Zq08FNmn3bbIWHhnsdxzWrV+d+W0a2AACAX2w5uEXD5w3XzQ1vVttabb2O45qFC6VWrXK/PWULAAD4xaMzH5UkvXz1yx4ncc+MGVLbtlLp0rl/D2ULAADk2eyNszVp1SQNbjVY1UtV9zqOK/77X6lTJ6lePWn+/Ny/j7IFAADyJC0zTQ9894Bqla6lf1/+b6/juOKNN6TevZ3pw7lzpYoVc/9eTpAHAAB58sL8F7Qmfo2m956uiLAIr+P4lbXSk09Kzz0n3XCD9PnnUsRZ/oqULQAAcM7Wxq/Vc788px4X9tB1da/zOo5fpadLAwZIn34q3XWX9M47Umjo2e+HaUQAAHBOrLW6Z9o9igyP1GvXvuZ1HL9KTJQ6dnSK1lNPSe+9d25FS2JkCwAAnKPRy0Zr7ua5eq/Te6pY/CxOYgpwO3ZI118vrVolffKJ1K9f3vZH2QIAAGdt35F9enTmo2pVvZXuvOhOr+P4zYoV0nXXSQcPStOmSVdfnfd9Mo0IAADO2sM/PKyktCS93+l9hZiCUSfmzHG+bZiVJf3yi3+KlkTZAgAAZ2nG+hka++dYDWo1SA2iG3gdxy/GjZOuuUaqUkX67TcpJsZ/+6ZsAQCAXDuYelB3TL5DDaMbasgVQ7yOk2fWSsOHS337Si1bSvPmSdX9vCYr52wBAIBce3jGw9p9eLe+7vF10K+plZYm3X23843DXr2ck+GLFvX/cRjZAgAAuTJt3TR9suwTPX7542pepbnXcfIkPl7q0MEpWsOGOdOIbhQtyQ9lyxhTzRgzxxiz2hiz0hjz4NHnzzPGzDTG/HX0vkze4wIAAC8kpCToril3qVH5RvrPVf/xOk6erF4tXXKJ9Pvvzorw//d/kjHuHc8fI1uZkh6x1jaQdKmk+40xDSUNkjTbWltX0uyjjwEAQLB48UXnK3qSHvz+Qe09sldfRj+goq+87nGwc/fDD9Jll0mHDzvXOOzVy/1j5rlsWWt3WWuXHP05SdJqSVUkdZX06dHNPpXULa/HAgAA+ah5c6l7d83/9FmN+XOMPijeW/XuHeo8H4RGjXIWK61e3RnVuvTS/DmuX0+QN8bUlNRM0kJJFay1uySnkBljyvvzWAAAwGVt2ujAp++qfvdb9O5V5dXv9++kiROlNm28TnZWMjOlhx+W3nxT6tTJmTosUSL/ju+3E+SNMcUlfSlpoLU28SzeN8AYs9gYs3jfvn3+igMAAPIo22arT8KHOhJmdff0vTL33ht0RevAAecah2++KT3yiPTNN/lbtCQ/lS1jTLicojXOWvvV0af3GGMqHX29kqS9Ob3XWvu+tTbWWhsbHR3tjzgAAMAP3lz4plJnfq+qR0Kcubd33jl+DlcwWL7cmfGcO1f68ENp5Mhzv5h0Xvjj24hG0keSVltrXznhpcmSbj/68+2Svs3rsQAAQP74c8+fmvb+o/rmqyIKadhIqlXLmULs3j0oCtekSc6J8CkpTtm64w7vsvjjnK3LJd0qabkxZtnR54ZIGiFpojHmDklbJd3ih2MBAACXpWSkqNeXvdRzT4Syx38qc/6Fzgv16zuFa9GigJ1OzMqS/vMf6fnnnbI1aZJUubK3mfJctqy18ySdbnWKdnndPwAAyF+PzXxMq/atUovXv1eZOtec/GKbNrkqWj179pS1Vps3b9bu3bs1atQodezY0aXEjoMHpT59pOnTpTvvlN56y72FSs8GK8gDAIDjvlnzjd5a9JYGXjJQ1xwrWlOmOLezEBcXp9q1a2vhwoUaN26cnnrqKRfS/m3VKqlFC2cdrXfekd5/PzCKlsS1EQEAwFEbEzaq3zf9FFs5ViPaj/j7hZdfdu47d87VflJSUhQfH69hw4ZJkho2bKiEhAR/xz3uq6+k22+XoqKc08latXLtUOeEkS0AAKDUzFTdPPFmhZgQfXHLFyoadu7DQitWrFDdunUVEeFcqHrJkiVq2rSpv6Iel5HhLOdw001Sw4bS4sWBV7QkRrYAAICkgd8P1NLdSzWl1xTVLF0zT/uKi4vT1q1blZqaqqysLA0bNkwvvviif4IetWOH1KOHNH++dP/9zuBboEwb/i/KFgAAhdzYP8fqvT/e06DLB6lTvU553l9cXJz69Omj1q1bKzExUUOGDNHll1/uh6SO2bOl3r2lI0ec1eDz4/qGeUHZAgCgEFu1b5Xunnq3rqpxlZ5p+4xf9hkXF6cPPvhAL7zwgl/2d0x2tjR8uLO0Q/36zvpZDRr49RCuoGwBAFBIHUw9qBsm3KASRUrovzf9V2Ehp6kFY8ac1X43bNigunXr+iHh3w4ckG691VnWoXdv6b33pOLF/XoI11C2AAAohLKys9Try17alLBJP97+oyqVqHT6jatVO6t979ixI4/pTrZwoXN+1q5d0qhR0j33SOZ0K3wGIL6NCABAITRk9hB9v/57vX3922pV/Qxf4Zswwbnls+xs6YUX/v6G4bx50r33BlfRkhjZAgCg0Bn35zi9uOBF3Rd7n+66+K4zv+Gdd5z7Hj3cDXaC3bul226TZs6UbrnFWaS0dOl8O7xfUbYAAChEFu9crDun3Kkra1yp1659zes4OZoxwylaSUlOybrzzuAbzToR04gAABQSO5N26oYJN6h8VHlNumWSwkPDvY50kvR06bHHpGuvlcqXdxYpveuu4C5aEiNbAAAUCofTD6vT552UkJKgef3nKToq2utIJ9m4UerZU1q0yDkB/pVXpGLFvE7lH5QtAAAKuMzsTPWc1FNxe+I0pdcUxVSM8TrScdY6K0s88IAUGipNmuRcfqcgoWwBAFCAWWv14HcPatpf0zTq+lG6vu71Z7+TSZP8H0xSfLwzivXll9IVVzilq0YNVw7lKcoWAAAF2Ku/vapRi0fp0cse1b3N7z23nZQr599Qkr77TurfX9q/X3rxRenhh52RrYKIE+QBACigvlz1pR794VHd1OAmvdAhD5fOGT3aufnBkSPOhaOvv97pcIsWSf/+d8EtWhJlCwCAAmn2xtnq/VVvXVr1Uo25YYxCTB7+yvdT2fr9d+mii5xlux55xClaTZvmebcBj7IFAEABs2jHInWb0E31ytbT1N5TVSzc26/1ZWRI//d/UsuWUkqKNHu2NHKkFBHhaax8wzlbAAAUIKv3rdZ1465TdGS0ZvSdofOKnedpnmXLpH79pLg4qU8f6a23gncl+HPFyBYAAAXE1kNbdfXYqxUWEqYfbv1BlUtU9ixLero0bJjUvLlz6Z2vv5bGji18RUtiZAsAgAJh9+Hd6jCmg5LSkvTzP35WnfPqeJZlyRJnNGv5cqlvX+n116XzvB1g8xRlCwCAILf3yF61/bStdiTu0Iy+M9SkQhP/HmD69FxtlpYmPfOMNGKEc7mdyZOlzp39GyUYUbYAAAhi+47sU9tP22rLoS2a3nu6Lq9+uf8PEhl5xk0WLZL+8Q9p5UpnVOuVV6QyZfwfJRhxzhYAAEEqPjle7T5rp40JGzW111RdVfMqdw40apRzy0FSkvTgg9Kll0oHDzqDYJ98QtE6EWULAIAgdCDlgDqM6aC/Dvylyb0mq02tNu4dbOJE5/Y/Jk+WGjaU3nzTuezOypXSdde5FyNYUbYAAAgye4/sVZtP22jVvlX6psc3al+7fb4ef8cO52LRXbs63y6cP196+22pVKl8jRE0OGcLAIAgsj1xu9p/1l7bErdpaq+p6nB+B9eONff661WqfXs1O/o4K0sadcsr2jV5lqaHT9fw4c5K8OHhrkUoEChbAAAEiY0JG9Xus3Y6kHJAM/rOUKvqrVw9Xqn27VXt0UeVVLu2TJmqeqD2K3pp66MaWnekVnwnnX++q4cvMChbAAAEgdX7Vqv9mPZKzUzV7NtmK7ZyrOvHbPbww1qQIsUOfUQ7Fa+X9LO+7T1S7459WMa4fvgCg3O2AAAIcL/v+F1Xjr5SWdlZ+qnfT/lStLKznWtP3/DGwxqhkqqpQ1p2aSvdOY6idbYoWwAABLCp66aqzadtVLJoSf3yj1/UqHwj14+5dKnUqpWzbtZ1Ea/oXiXq9ZIl1WzhPC195RXXj1/QULYAAAhQH/zxgbqO76oG5RpoQf8Fqlu2rqvHS0iQ7r9fio2V1q+X3u3xikZue1SDa9fW182aadvIkar26KMUrrPEOVsAAAQYa62e+ukpPfXTU7quznWaeMtEFS9S3LXjZWZKH30kDR0qHTjgFK6nn5aW9Z6lbSNHasRtt0mSypUrp6WSDs2aJT38sGt5ChpjrfU6w3GxsbF28eLFXscAAMAz6Vnpum/affpo6Uf6R8w/9F6n9xQe6t7aCj/84CzfsGKFdMUVzgKlTZvmfb89e/aUtVabN2/W7t27NWrUKHXs2DHvOw4gxpg/rLVnPIGOaUQAAAJEfHK8rh5ztT5a+pGGXjFUH3X5yLWitXq11LGjdM01UnKyNGmS9NNPORet0aNHa/To0We1/7i4ONWuXVsLFy7UuHHj9NRTT/kneBBiGhEAgACwcu9Kdf5vZ+1M2qkxN4xR3yZ9XTlOfLz01FPSO+9IUVHSiy9K//qXVLTo6d9zrGj169cvV8dISUlRfHy8hg0bJklq2LChEhIS8pjcmV41QfhVSEa2AADw2NR1U3XpR5cqJTNFP/X7yZWilZYmvfKKVLeuc03pAQOck+D//W/fRetcrFixQnXr1lVERIQkacmSJWp6dMjsqaee0oMPPqhhw4Zp7969uuiiizR48GDdeOONys7OPuW5nTt3qmXLlhoxYoR2796thx9+WPfff7+GDh2qffv26R//+Ie2b9+u/v37KyMjw7+/iJ9QtgAA8Ii1Vi/Of1Fd/ttF9crW06K7FumSqpf49RhZWdKYMdIFFzjnZl1yifTnn07hio7266GOi4uL09atW5WamqojR45o2LBheuihh7Rjxw5lZGSodOnS+u2337Ro0SL16tVLw4cPV/ny5bV///5Tnps9e7Z69uypQYMGacyYMerdu7fefvttrVmzRtHR0apevboeeeQRvfHGGwoP0OsGMY0IAIAHEtMS1f/b/vpy9ZfqfmF3fdL1E0WGR/pt/9ZK06dLgwdLy5dLzZpJ770nXX213w5xWnFxcerTp49at26txMREDRkyRJdffrn69++v119/Xfv27dO2bdu0aNEitWzZUpJ06NAhRUdHn/Lc9u3b1a1bN0nSypUr9eCDDyo9PV2RkZE6fPiwNm7cqLCwMBUv7t63NfOKsgUAQD5bvme5bpp4kzYmbNTIDiP18GUP+/VcpAULpEGDpF9+ca5f+N//St27SyH5NJ8VFxenDz74QC+88MJJz1944YUaOXKk9u/fr2bNmmnBggXav3+/vv76a915552SpLVr15703Lhx41S/fn1JUvfu3TVgwABFRkbq3//+t/71r3/p2Wef1cSJEzV37ly1bt06f37Bs8TSDwAA5KMxcWN099S7VSqilCbcPEFX1rjSb/teuVIaMkSaPFmqUEEaNky6804pr7NrycnJkqTIyNyNvFWpUkXbtm1TyBna3a233qoxY8ac8blAldulHyhbAADkg5SMFD004yG998d7uqrGVRp/83hVLF7RL/tev1565hlp7FipeHHpscekgQOdbxvCPbktW0wjAgDgsj/3/KneX/bWyn0r9fjlj+vZts8qLCTvfwVv2PB3ySpSRHroIeccrbJl/RD6BKNGjZIk3Xffff7dcSFB2QIAwCXWWr35+5t6bOZjKlOsjGb0naGrz8/7GeobN0rPPit99pkzRfivfzmjWRX9M1B2iokTJ0qibJ0ryhYAAC7Yc3iP/vHtP/Td+u/UqV4nfdzlY0VH5W2thU2bpOeekz79VAoNlR54QHr8calSJT+FhisoWwAA+NlXq7/SPVPvUVJ6kt6+/m3dG3tvnr5tuH699MIL0ujRTsm67z6nZFWu7L/McA9lCwAAP4lPjtcD0x/QhJUTdFGli/RZt890YfkLz3l/y5ZJw4c71y0MD5fuucdZ0qFKFf9lhvsoWwAA+MGXq77UvdPu1cHUg3qmzTN6/PLHz+ki0tY662MNHy59/71UooRzSZ2BA907JwvuomwBAJAHew7v0YPfP6gJKyeoWcVmmnXbLDWp0OSs92OtNG2aU7IWLHAupfPcc86UYenS/s99NubOnettgCBH2QIA4Bxk22y9/8f7GjRrkJIzkvV066c1qNWgsx7NSkuTxo+XXn7ZuaxO9erSm29K/ftLuVxDFAGOsgUAwFmK2x2nu6ferYU7FqpNzTYa1XGULih3wVntIz5eevdd6e23pd27pQsvdJZy6Nkz7yu++9vIkSMlSY8++qjHSYITZQsAgFw6nH5Yw+YM0+sLX9d5xc7TmBvGqE/jPmf1TcPVq6XXXnOKVWqqdO21zmKkHTpIfrw8ol9NnTpVEmXrXFG2AAA4g2ybrTFxYzR49mDtOrxLAy4aoOHth+u8Yufl6v3WSrNmSa+84pz0HhEh3Xqrc9J7w4buZof3KFsAAPgwb+s8PTTjIS3euVgtqrTQVz2+0qVVL83VexMTpTFjpFGjpFWrnItDP/20s4RDdN7WN0UQoWwBAJCDzQc36/FZj2viyomqUqKKxtwwRr0b91aICTnje5cvdwrWmDHSkSPSxRdLn3wi9eolFS2aD+ERUChbAACcID45XiPmjdBbv7+lEBOiYVcN079b/ltRRaJ8vi89XfrqK6dk/fKLM1XYs6ezdEPz5vkU3iXFihXzOkJQo2wBACApKS1Jr/72qkYuGKkjGUfUt0lfPdvmWVUrVc3n+zZtkj76SPrgA2nvXun886WRI6V+/aSyZfMnu9u+++47ryMENcoWAKBQS81M1TuL3tHz855XfHK8bmxwo55p84waRp/+zPXUVGcU66OPpB9/dL5F2KmTM4p19dVSyJlnGlGIULYAAIVSckayPlzyoV6c/6J2JO1Q+9rt9Xzb59W8yunn/JYtcwrWuHFSQoJUs6Zzwnu/flI13wNgQe2ZZ56RJD355JMeJwlOlC0AQKGSlJakUYtG6eVfX9a+5H26qsZV+uyGz9S2Vtsct09IkP77X6dkLVninOB+443SHXdIbdoUjlGs2bNnS6JsnSvKFgCgUEhISdAbC9/Q6wtfV0Jqgq45/xo9ccUTuqLGFadsm5bmXKdw7FjnPj1diolxLqPTu7d0Xu6W1wIkUbYAAAXcxoSNev231/Xxso91OP2wutbvqieueOKU6cLsbGnePKdgffGFdPCgsy7W/fdLfftKF13kTX4EP8oWAKDAsdZq/rb5euXXV/TNmm8UGhKqno166t8t/60mFZqctO2qVc45WOPGSVu2SFFRzjRh375S27ZSGH9TIo/4RwgAUGCkZ6Vr0qpJevW3V7V452KdV+w8DW41WPe3uF+VS1Q+vt3atc7o1cSJzgKkoaHOtwiff17q2tUpXPhb2YKyhoVHKFsAgKC34cAGfbDkA3289GPtS96n+mXr652O7+i2prcpMjxS0qkFS5Iuv1x6/XWpRw9nyhA5+/LLL72OENQoWwCAoJSRlaEp66bo3cXvaubGmQo1oepSv4vuvvhudTi/g0JMSI4Fq1Urp2DddJNUpYq3vwMKB8oWACCorI1fq8/iPtMnyz7RrsO7VK1kNT3d+mndcdEdqhhVWQsXSkM+lCZPllavdhYcPTaCRcE6N4MHD5YkDR8+3OMkwYmyBQAIePuT92vCygn6LO4zLdyxUCEmRNfVuU73xN6jqypfpzk/huo/D0lTpjiXzAkLk1q3lu691znZnYKVN7/++qvXEYIaZQsAEJDSMtP03frv9FncZ5q6bqoysjPUpEITvXz1y2pbvpeW/FRJHz4mdf9BSkmRSpaUrr/eOcH92mul0qW9/g0AB2ULABAw0jLTNHPjTE1cOVGT107WobRDqli8oh5o/i9dmHWr1s9rqs/elR6Jc7avXt1Zyb1rV+nKK6UiRbzND+SEsgUA8FRaZpp+2PCDvlj1hb5d+60S0xJVJqKMrq52oyru764tc9rrw6fClJTkTA+2aiWNGOGMXjVp4pyTBQQyyhYAIN8dTD2o79d/rynrpmjquqlKTEtU6aJl1DzqJhWP7661X7TVFyudYarq1Z1L5Fx7rbPIaMmSHocvhKpWrep1hKBG2QIA5Iu/9v+lKeumaMq6Kfplyy/KslkqGVZONVJvVuoft2j9zLaanVlEERHSFVdIA+5wCtYFFzB65bWxY8d6HSGoUbYAAK5IzkjWL1t+0Q8bftDUv6Zq3f51kqTytrEqbX5Mu+Z2VuLWFloTFqpLL5V6PyG1aSNdeqlUtKjH4QE/omwBAPwiKztLS3cv1cwNMzVz40zN3zpf6dnpCrVFVDy+tUL++Key13RSfGJNxcZKfXs504ItW3J5nEA3cOBASdJrr73maY5gRdkCAJwTa61Wx6/WL1t+0ayNszRzw486lH5AklTkQFOlr/mntKGDQndfocbNItWqtXT5E84UYalS3mbH2Vm2bJnXEYIaZQsAkCuZ2Zlaumupft7ys2av/0Xzt81TYuZ+SVLI4SrK/quLtLGDyiS00xUXVVCr9tLlw6SLL2ZaEIUbZQtA0EtJSdG1116rH3/8UWPHjtWzzz4rSRo6dKhuv/32U7Y/cOCAevTooc2bN6tmzZqaOHGiypQpo+XLl+vll1/W6NGj8/k3CEwJKQlatHORft74m2as/kV/JvyqdB1xXjxwvrSls7T1StUJv0KXNzxfrToZtWol1a/PCe3AiShbAILexx9/rBtvvFGHDh3SU089pcWLF8sYo4svvlhdunRRmTJlTtp+xIgRateunQYNGqQRI0ZoxIgReuGFF9S4cWNt375dW7duVfXq1T36bbyRmpmquN1xmvPX75q95nfFxf+ufdnOCe2yRtrTWNraT+WSr9AlFa/QVRdVVvMbnFGrEiW8zQ4EOsoWgKA3btw4ff7555oxY4Y6dOig8847T5LUoUMHff/99+rVq9dJ23/77beaO3euJOn2229X69at9cILL0iSOnfurPHjx+uxxx7L198hP6VkpGj5nhWavTJOP69bqj/3/65dNk42JMPZIKmitOMSFT90uy4o0UJtL4jVFR1Lq3lzqUIFb7PDG/Xq1fM6QlCjbAEIaunp6dq4caNq1qypSZMmqVq1asdfq1q1qnbs2HHKe/bs2aNKlSpJkipVqqS9e/cefy02NlYjRowoEGXLWqtdSbs1Z3Wc5qxepqU747QxOU4Hw9ZKJtvZKK24tCtWZVMeVoOSLXTF+S101aVVFBNjKFY47v333/c6QlCjbAEIavHx8Sp99IrD1tpTXjdnefJQ+fLltXPnTn9EyzfWWu0+vEe/rF6tX9as1rLtq7UxaZX2mRXKKPp3kdTB6grfH6Nq5mY1LBujK+o2VfuLa6lxoxBFRnqXHyjoKFsAglqxYsWUmpoqyRnJOjY9KEnbt29X69atT3lPhQoVtGvXLlWqVEm7du1S+fLlj7+WmpqqYsWKuR37nGRkZWrppq36eeVaLd6ySmviV2t76modDF+lrCIH/94wrbjCDjVQ2ezrVScqRrHVYtS2YRNd0rSMypfn5HWcvQEDBkhihOtcUbYABLUyZcooKytLqampuuaaazRkyBAlJCRIkn744QcNHz5ckjR48GC1aNFCN9xwg7p06aJPP/1UgwYN0qeffqquXbse39+6devUqFEjT34XSTqSlqrf1m7S7+vXa8WODVoXv17bk9crQRuUVmyzFJr598ap0Sp2uIEqmx6qXbyhGldqoCsbNNBVzaqofHkaFfxn3bp1XkcIapQtAEHv6quv1rx589S+fXs9+eSTat68uSTpP//5z/GT5ZcvX64uXbpIkgYNGqTu3bvro48+UvXq1fXFF18c39ecOXPUsWNH17IeSUvV4r+2asmGrVq5fas2xG/RjsNbtS9zs5LCNigrartkTpgOTS+pIml1VMZepCoht6h+uTpqVqOuWl/YUDH1yio83LWoAPzE5HSOg1diY2Pt4sWLvY4BIMgsXbpUr7zyisaMGXPaba655hrNmDHD537S0tJ01VVXad68eQoLO/v/F008kqY/N+7W8i07tW7XTm2M36btSVu1J22LDtqtSimyVdnF9p78JmsUcqSyojKrq2xIHVWLPF/1ytVRk2p1dEnd89W0bllFRDBKBW8dm44/cZoekjHmD2tt7Jm2Y2QLQNBr1qyZ2rRpo6ysLIWGhh5/fty4cXriiSeOr5s1btw49enT57T72bp1q0aMGHFS0bJW2nMgRSu37NbaHTu1Ye8ubU3YpV2Hd2pfyi4dyNipw9qltKI7ZSMOnLpTE6kitoaKZ1VX5fSLVCm8us4vW10Nq1TXRXWqK7ZeFZWMKuLXzwNAYKFsASgQ+vfvf9LjcePGacCAAUpOTpYkbdmyRQMGDFB2tnRlh87asHufNu3Zp20H9mnnwX3albhPew/v04HUeB369iUd0T6lhe5TVtF9UpEjpx4wO1whWRUVkVVZpUwdlc26UhVNJdUoU1nnl6+kC6pW0sV1qqt6dBmFhDAyheAWExPjdYSgxjQigKBjrXQgMVXb4hO0PT5BOxMStPtQgvYcStC+pATFH0nQvJdfVkZi0qlvLiXpodPsOCNCoWnRKpIZrShFq2RYtMoUjVb5qGhVLVVRtcpVVt1KldSgamXVrXqeioSHuPlrAghwTCMCCChZWVLS4SztO3RYexMTFZ+YpP2Hk3TgcKISkpN04HCSDiQn6lBKkhLTkpSUnqgjGUlKyUpSqk1SuklUhklSZtgh2aIJUniq7wMmnub5Q1Lnoi+pYoloVS4drWplo1UzOlp1q0SravkohdCfAPiZ62XLGHOtpNclhUr60Fo7wu1jAjg7aenZOpCYpoSkVB06nKaDh1N16EiqDiYn69DRW2LKER1OTdbh9GQdSU9WcnqykjOPKCUzWalZyUrLTlaaTVaGPaIMk6xMk6wsk6zs0GRlhx2Rwo9IRZJzmShcRiUUakoqLKSEitiSijJlVSykpqJCS6qkKaPSYWVUplgZlYsqo/Ily6hCyTKqWraMalQoo1qVSqvhJ3W0ZcuWU/Zco0YNTR70qH8/QKCA69u3ryRp7NixHicJTq6WLWNMqKS3JXWQtF3SImPMZGvtKjePCwQKa62ybJbSs9KVkZWhjOwMn/fpWeknPZeelaHk1AwdSc3QkRTn5+Q053FyWrqOpKUqJT3NuWWkKjXTuU/LSlN6VprSs1OVYdOUYVOVadOUqVRlmTRlmzRlh6TKhqZJoaknr92UW0UkhYXIZEYqJCtKodmRCrWRCreRilCkiphoFQ2JVERIlCLCIlUsLFLFw0uoRJGSKhlRQqUjSqpUsRIqE1VCZaNKqnzpEqp0XglVKVdSZUoWzfPCm88999xJ52xJUmRkpJ577rm87RgohLZv3+51hKDm9shWC0nrrbUbJckYM15SV0mUraOstbKyrv987D7bZh+/Wfv3Y3+8duLzbr6WlZ2lLJt1yn1mdmaunsvKzlKmPc3z2Zlnve/MrCxlZGUqPevvkpSZnaFM69xclx0iZUZIWUWd+8yiMtlFFZIdoVAVVaiNUJhKqIiJVlRIURUxESoSWlRFQyNU1BRVREhRFQuJULHwooosEqGookUVVTRCURFFVapYlEpFRap0ZKTKFI/SecUjdV4J51YiIlJFQouc9eVw8suxbx2e+G3E5557zue3EQHADW6XrSqStp3weLukS0638Yq9K1T/rfrHr292rCQUxJ/hnhATolATqtCQUIWFhB3/+e/7MIUoVMaGSsdu2aEyNkw2K1Q2O1TZWaGyWaHKzgxVdmYRZWUWU1ZmqLIyQpWZ4dxnZYQpK+PY+8OkrHApOzzn+6wix38ODw1XkbBwFQ0PV0R4uCKKOPfFioSrWNFwRRYposiIcEUWDVdkRLiiIsIVVSxcxYuFKyoyXCUjw1U8Mlwlo4qodPEIlSwepshIKTJSioqSIiK4HMsxffr0oVwB8JzbZSun/+Sf1DSMMQMkDZCkyCqRalaxmYwxMkffetY/yxz/P21+PvnnEBNy/Gb09+MTXzvx+UB9TTZEyYdDlZQYqsRDoUo6FKrEg2FKPBiqgwkhOnjQKCFBp9wOHpT2H5IycjnYVKKEcytd0rkvefS+eHGn1Jx4O1Z0zvQ4MlKcgA0AhYzbZWu7pGonPK4qaeeJG1hr35f0vuQs/TD+5vEuR0KgSU2V9u07+RYff/rnDh6UsrNPv7/QUKlMmb9v550nnX++83OpUn+XJl/3UVGUIgA45rLLLvM6QlBzu2wtklTXGFNL0g5JPSX1dvmYCABZWdLevdKuXafedu507vfudQrU4cM57yM0VCpXToqOdu6bNHF+Pu+8k8vU/96iophGAwB/OnZBd5wbV8uWtTbTGPOApBlyln742Fq70s1jwn3JydK2bc5t61bntmPHyYVqz56cR5/KlpUqVXJudes65enY7VixOnYrXZrRJQBA8HN9nS1r7XRJ090+DvzDWmn3bmnTppPL1LGft21zRqNOZIxUocLfJapZs79/PvFWsaJUtKg3vxcA4NzddNNNkqQvv/zS4yTBiRXkC6HUVGnzZmnDBmnjxr/vj91SUk7evlQpqXp1qVo16ZJL/v65enXnVrmyVITr6AJAgbV//36vIwQ1ylYBlZHhjE6tWSOtXevcr1/vFKsdO07eNipKql3bmda75hrn51q1pBo1nFJVsqQ3vwMAAAUBZSvIHTjglKljhepYuVq/Xso8YVHwChWcMtW+vVOmatd2vqFXu7ZUvjwnlAMA4BbKVpA4ckRatUpavvzv24oVzonox4SHS3XqSA0aSN26SRdc4Nzq13dONgcAAPmPshVgsrKkv/46uVQtX+6cS3V0AXoVKyZdeKF0/fVSw4Z/F6pataQw/kQBAH7Wrl07ryMENXPsEjKBIDY21i5evNjrGPkmM1NavVr64w/ntmSJtGyZs7SC5Cx7ULeu1Lixc2vUyLmvXdtZgwoAAHjHGPOHtTb2TNsxDpJPMjKcab8lS/4uVnFxzjcDJeck9WbNpDvvdO6bNHGmA4sV8zY3AADIG8qWS3bulH77Tfr1V+d+8eK/i1XJkk6huu8+6aKLpIsvdkawGK0CAASi6667TpL03XffeZwkOFG2/CAtzRmp+u23vwvWtm3Oa0WKOIXq3nulFi2cYnX++ayMDgAIHin/uwAjzgpl6xwkJkrz50s//+zcFi+W0tOd12rUkFq2lC69VLrsMikmhlXTAQAozChbubB/v/TLL3+Xq6VLnev+hYVJzZtLDz7oFKtLL3UuSwMAAHAMZSsHBw9Kc+ZIs2Y55WrFCuf5iAinUA0dKl11lXPpmqgoT6MCAIAAR9mSc87Vb79JM2c6BWvRImfkKipKatVK6tXLKVexsUwJAgAKn06dOnkdIagVynW2rHVGq46Vq59+cta2Cg11TmLv0MG5rM0ll3CBZQAAkDPW2fofSUlOsZo2TZo+Xdq1y3n+gguk/v2dctW6tVSqlKcxAQBAAVOgy9a6dU65mjbNOfcqI8NZ4+qaa6TrrnMKVrVqXqcEACCwtW7dWpI0d+5cT3MEqwJVtjIynFI1ebIzerV+vfN8w4bSwIFSx47Osgzh4Z7GBAAAhUjQl63kZOmHH6Svv5amTJESEpxvDbZtKz30kHOx5po1vU4JAAAKq6AsWwcPSlOnOgXr+++dwlW6tNS5s3TjjdLVV0uRkV6nBAAACKKydeCA9NVX0hdfSD/+KGVmOguI9usn3XCDszQD04MAACDQBHTZSkyUvv1WGj/emSrMzJTq1JEeftgpWC1acI1BAADc1r17d68jBLWAK1vJyc63B8ePd+7T0qTq1Z3zr3r2lJo1k4zxOiUAAIXHfffd53WEoBZQZWvTJql8eenIEaliRenuu52CdckljGABAOCV5ORkSVIkJ0Sfk4AqW4cOSXfc4RSsK690VnQHAADeuv766yWxzta5Cqiy1bSp9N57XqcAAADwn4CanONcLAAAUNAEVNkCAAAoaChbAAAALgqoc7YAAEDg6devn9cRghplCwAA+ETZyhumEQEAgE/x8fGKj4/3OkbQYmQLAAD4dPPNN0tina1zxcgWAACAiyhbAAAALqJsAQAAuIiyBQAA4CJOkAcAAD7de++9XkcIapQtAADgU48ePbyOENSYRgQAAD5t27ZN27Zt8zpG0GJkCwAA+HTrrbdKYp2tc8XIFgAAgIsoWwAAAC6ibAEAALiIsgUAAOAiTpAHAAA+PfLII15HCGqULQAA4FPnzp29jhDUmEYEAAA+rV27VmvXrvU6RtBiZAsAAPh09913S2KdrXPFyBYAAICLKFsAAAAuomwBAAC4iLIFAADgIk6QBwAAPg0dOtTrCEGNsgUAAHxq37691xGCGtOIAADAp2XLlmnZsmVexwhajGwBAACfBg4cKIl1ts4VI1sAAAAuomwBAAC4iLIFAADgIsoWAACAizhBHgAA+PT88897HSGoUbYAAIBPLVu29DpCUGMaEQAA+LRgwQItWLDA6xhBi5EtAADg05AhQySxzta5YmQLAADARZQtAAAAF1G2AAAAXETZAgAAcBEnyAMAAJ9ee+01ryMENcoWAADwKSYmxusIQY1pRAAA4NOsWbM0a9Ysr2MELUa2AACAT88++6wkqX379h4nCU6MbAEAALiIsgUAAOAiyhYAAICLKFsAAAAu4gR5AADg03vvved1hKBG2QIAAD7Vr1/f6whBjWlEAADg05QpUzRlyhSvYwQtRrYAAIBPL7/8siSpc+fOHicJToxsAQAAuIiyBQAA4CLKFgAAgIsoWwAAAC7iBHkAAODTmDFjvI4Q1ChbAADAp2rVqnkdIagxjQgAAHyaMGGCJkyY4HWMoMXIFgAA8Omdd96RJPXo0cPjJMGJkS0AAAAXUbYAAABcRNkCAABwUZ7KljHmJWPMGmPMn8aYr40xpU94bbAxZr0xZq0x5po8JwUAAAhCeT1BfqakwdbaTGPMC5IGS3rcGNNQUk9JF0qqLGmWMaaetTYrj8cDAAD5bNKkSV5HCGp5Gtmy1v5grc08+vA3SVWP/txV0nhrbZq1dpOk9ZJa5OVYAADAG+XKlVO5cuW8jhG0/HnOVn9J3x39uYqkbSe8tv3ocwAAIMiMHj1ao0eP9jpG0DrjNKIxZpakijm89IS19tuj2zwhKVPSuGNvy2F7e5r9D5A0QJKqV6+ei8gAACA/HSta/fr18zRHsDpj2bLWtvf1ujHmdkmdJLWz1h4rVNslnbi2f1VJO0+z//clvS9JsbGxORYyAACAYJXXbyNeK+lxSV2stcknvDRZUk9jTFFjTC1JdSX9npdjAQAABKO8fhvxLUlFJc00xkjSb9bae6y1K40xEyWtkjO9eD/fRAQAAIVRnsqWtbaOj9eek/RcXvYPAAAQ7LgQNQAA8Gn69OleRwhqlC0AAOBTZGSk1xGCGtdGBAAAPo0aNUqjRo3yOkbQomwBAACfJk6cqIkTJ3odI2hRtgAAAFxE2QIAAHARZQsAAMBFlC0AAAAXsfQDAADwae7cuV5HCGqMbAEAALiIsgUAAHwaOXKkRo4c6XWMoEXZAgAAPk2dOlVTp071OkbQomwBAAC4iLIFAADgIsoWAACAi1j6AQAA+FSsWDGvIwQ1yhYAAPDpu+++8zpCUGMaEQAAwEWULQAA4NMzzzyjZ555xusYQYuyBQAAfJo9e7Zmz57tdYygRdkCAABwEWULAADARZQtAAAAF7H0AwAA8Kls2bJeRwhqlC0AAODTl19+6XWEoMY0IgAAgIsoWwAAwKfBgwdr8ODBXscIWkwjAgAAn3799VevIwQ1RrYAAABcRNkCAABwEWULAADARZyzBQAAfKpatarXEYIaZQsAAPg0duxYryMENaYRAQAAXETZAgAAPg0cOFADBw70OkbQYhoRAAD4tGzZMq8jBDVGtgAAAFxE2QIAAHARZQsAAMBFnLMFAAB8qlevntcRghplCwAA+PT+++97HSGoMY0IAADgIsoWAADwacCAARowYIDXMYIW04gAAMCndevWeR0hqDGyBQAA4CLKFgAAgIsoWwAAAC7inC0AAOBTTEyM1xGCGmULAAD49Nprr3kdIagxjQgAAOAiyhYAAPCpb9++6tu3r9cxghbTiAAAwKft27d7HSGoMbIFAADgIsoWAACAiyhbAAAALuKcLQAA4NNll13mdYSgRtkCAAA+DR8+3OsIQY1pRAAAABdRtgAAgE833XSTbrrpJq9jBC2mEQEAgE/79+/3OkJQY2QLAADARZQtAAAAF1G2AAAAXMQ5WwAAwKd27dp5HSGoUbYAAIBPTz75pNcRghrTiAAAAC6ibAEAAJ+uu+46XXfddV7HCFpMIwIAAJ9SUlK8jhDUGNkCAABwEWULAADARZQtAAAAF3HOFgAA8KlTp05eRwhqlC0AAODTo48+6nWEoMY0IgAAgIsoWwAAwKfWrVurdevWXscIWpQtAAAAF1G2AAAAXETZAgAAcBFlCwAAwEUs/QAAAHzq3r271xGCGmULAAD4dN9993kdIagxjQgAAHxKTk5WcnKy1zGCFiNbAADAp+uvv16SNHfuXG+DBClGtgAAAFxE2QIAAHARZQsAAMBFlC0AAAAXcYI8AADwqV+/fl5HCGqULQAA4BNlK2+YRgQAAD7Fx8crPj7e6xhByy9lyxjzqDHGGmPKnfDcYGPMemPMWmPMNf44DgAAyH8333yzbr75Zq9jBK08TyMaY6pJ6iBp6wnPNZTUU9KFkipLmmWMqWetzcrr8QAAAIKJP0a2XpX0mCR7wnNdJY231qZZazdJWi+phR+OBQAAEFTyVLaMMV0k7bDWxv3PS1UkbTvh8fajzwEAABQqZ5xGNMbMklQxh5eekDRE0tU5vS2H52wOz8kYM0DSAEmqXr36meIAAAAElTOWLWtt+5yeN8Y0llRLUpwxRpKqSlpijGkhZySr2gmbV5W08zT7f1/S+5IUGxubYyEDAADeuffee72OENTO+QR5a+1ySeWPPTbGbJYUa62NN8ZMlvS5MeYVOSfI15X0ex6zAgAAD/To0cPrCEHNlUVNrbUrjTETJa2SlCnpfr6JCABAcNq2zTkNu1q1amfYEjnxW9my1tb8n8fPSXrOX/sHAADeuPXWWyVJc+fO9TZIkGIFeQAAABdRtgAAAFxE2QIAAHARZQsAAMBFrnwbEQAAFByPPPKI1xGCGmULAAD41LlzZ68jBDWmEQEAgE9r167V2rVrvY4RtBjZAgAAPt19992SWGfrXDGyBQAA4CLKFgAAgIsoWwAAAC6ibAEAALiIE+QBAIBPQ4cO9TpCUKNsAQAAn9q3b+91hKDGNCIAAPBp2bJlWrZsmdcxghYjWwAAwKeBAwdKYp2tc8XIFgAAgIsoWwAAAC6ibAEAALiIsgUAAOAiTpAHAAA+Pf/8815HCGqULQAA4FPLli29jhDUmEYEAAA+LViwQAsWLPA6RtBiZAsAAPg0ZMgQSayzda4Y2QIAAHARZQsAAMBFlC0AAAAXUbYAAABcxAnyAADAp9dee83rCEGNsgUAAHyKiYnxOkJQYxoRAAD4NGvWLM2aNcvrGEGLkS0AAODTs88+K0lq3769x0mCEyNbAAAALqJsAQAAuIiyBQAA4CLKFgAAgIs4QR4AAPj03nvveR0hqFG2AACAT/Xr1/c6QlBjGhEAAPg0ZcoUTZkyxesYQYuRLQAA4NPLL78sSercubPHSYITI1sAAAAuomwBAAC4iLIFAADgIsoWAACAizhBHgAA+DRmzBivIwQ1yhYAAPCpWrVqXkcIakwjAgAAnyZMmKAJEyZ4HSNoMbIFAAB8eueddyRJPXr08DhJcGJkCwAAwEWULQAAABdRtgAAAFxE2QIAAHARJ8gDAACfJk2a5HWEoEbZAgAAPpUrV87rCEGNaUQAAODT6NGjNXr0aK9jBC3KFgAA8ImylTeULQAAABdRtgAAAFxE2QIAAHARZQsAAMBFLP0AAAB8mj59utcRghplCwAA+BQZGel1hKDGNCIAAPBp1KhRGjVqlNcxghZlCwAA+DRx4kRNnDjR6xhBi7IFAADgIsoWAADI0Ysvvqg5c+ac9NycOXP04osvepQoOFG2AABAjpo3b67u3bsrISFBklO0unfvrubNm3ucLLhQtgAAQI7atGmjiRMnavXq1dq0aZO6d++uiRMnqk2bNl5HCyrGWut1huOMMUmS1nqdIwCVkxTvdYgAxOeSMz6XU/GZ5IzPJWd8LqeqLKmSpF2SdnqcJZDUt9aWONNGgbbO1lprbazXIQKNMWYxn8up+FxyxudyKj6TnPG55IzPJWd8LqcyxizOzXZMIwIAALiIsgUAAOCiQCtb73sdIEDxueSMzyVnfC6n4jPJGZ9LzvhccsbncqpcfSYBdYI8AABAQRNoI1sAAAAFSsCWLWPMo8YYa4wp53WWQGCMecYY86cxZpkx5gdjTGWvMwUCY8xLxpg1Rz+br40xpb3O5DVjzC3GmJXGmGxjTKH/5pAx5lpjzFpjzHpjzCCv8wQCY8zHxpi9xpgVXmcJFMaYasaYOcaY1Uf//XnQ60yBwBgTYYz53RgTd/RzecrrTIHEGBNqjFlqjJnqa7uALFvGmGqSOkja6nWWAPKStbaJtTZG0lRJ//E4T6CYKamRtbaJpHWSBnucJxCskHSjpJ+9DuI1Y0yopLclXSepoaRexpiG3qYKCKMlXet1iACTKekRa20DSZdKup9/ViRJaZLaWmubSoqRdK0x5lJvIwWUByWtPtNGAVm2JL0q6TFJnFB2lLU28YSHUeKzkSRZa3+w1mYeffibpKpe5gkE1trV1loWB3a0kLTeWrvRWpsuabykrh5n8py19mdJB7zOEUistbustUuO/pwk5y/QKt6m8p51HD76MPzojb9/JBljqkrqKOnDM20bcGXLGNNF0g5rbZzXWQKNMeY5Y8w2SX3EyFZO+kv6zusQCChVJG074fF28RcozsAYU1NSM0kLPY4SEI5OlS2TtFfSTGstn4vjNTkDQ9ln2tCTFeSNMbMkVczhpSckDZF0df4mCgy+Phdr7bfW2ickPWGMGSzpAUnD8jWgR870uRzd5gk50wDj8jObV3LzmUCSZHJ4jv8rx2kZY4pL+lLSwP+ZUSi0rLVZkmKOnhP7tTGmkbW2UJ/vZ4zpJGmvtfYPY0zrM23vSdmy1rbP6XljTGNJtSTFGWMkZ0poiTGmhbV2dz5G9MTpPpccfC5pmgpJ2TrT52KMuV1SJ0ntbCFZy+Qs/lkp7LZLqnbC46rium44DWNMuJyiNc5a+5XXeQKNtfagMWaunPP9CnXZknS5pC7GmOslRUgqaYwZa63tm9PGATWNaK1dbq0tb62taa2tKec/lBcVhqJ1JsaYuic87CJpjVdZAokx5lpJj0vqYq1N9joPAs4iSXWNMbWMMUUk9ZQ02eNMCEDG+T/8jySttta+4nWeQGGMiT72LW9jTDFJ7cXfP7LWDrbWVj3aVXpK+vF0RUsKsLIFn0YYY1YYY/6UM83K15Idb0kqIWnm0WUx3vU6kNeMMTcYY7ZLukzSNGPMDK8zeeXolycekDRDzgnPE621K71N5T1jzH8l/SqpvjFmuzHmDq8zBYDLJd0qqe3R/5YsOzpqUdhVkjTn6N89i+Scs+VzmQOcihXkAQAAXMTIFgAAgIsoWwAAAC6ibAEAALiIsgUAAOAiyhYAAICLKFsAAAAuomwBAAC4iLIFAADgov8HwfpYQP9txqwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Approximation order : 3\n", + "Approximation bias : 7.085536923187668\n" + ] + } + ], + "source": [ + "#### FOLDED CELL\n", + "text_offset = np.array([-0.15,2.])\n", + "\n", + "func='exp(x)'\n", + "order = 3\n", + "x1 = 3\n", + "x_min, x_max = -4,4\n", + "func_sp = parse_expr(func)\n", + "taylor_exp = series(func_sp,x,n=order+1)\n", + "approx = lambdify(x, sum(taylor_exp.args[:-1]), \"numpy\")\n", + "func_np = lambdify(x, func_sp, \"numpy\")\n", + "n_points = 1000\n", + "x_array = np.linspace(x_min,x_max,n_points)\n", + "approx_array = np.array([approx(z) for z in x_array])\n", + "func_array = np.array([func_np(z) for z in x_array])\n", + "func_x1 = func_np(x1)\n", + "approx_x1 = approx(x1)\n", + "plt.figure(42,figsize=(10,10))\n", + "plt.plot(x_array,approx_array,color='blue',label='Taylor Expansion')\n", + "plt.plot(x_array,func_array,color='green',label=func)\n", + "plt.plot(0,approx(0),color='black',marker='o')\n", + "plt.annotate(r'(0,0)',[0,approx(0)],xytext=text_offset)\n", + "plt.plot([x1,x1]\n", + " ,[-np.max(np.abs([np.min(func_array),np.max(func_array)])),min(approx_x1, func_x1)]\n", + " ,'--',color='black',marker='x')\n", + "plt.plot([x1,x1],[approx_x1, func_x1],'r--',marker='x')\n", + "plt.annotate(r'$p_{approx}$',[x1,approx(x1)],xytext=[x1,approx(x1)]-text_offset)\n", + "plt.annotate(r'$p$',[x1,func_np(x1)],xytext=[x1,func_np(x1)]-text_offset)\n", + "plt.xlim([x_min,x_max])\n", + "plt.ylim(-np.max(np.abs([np.min(func_array),np.max(func_array)]))\n", + " ,np.max(np.abs([np.min(func_array),np.max(func_array)])))\n", + "plt.legend()\n", + "plt.show()\n", + "print('Approximation order : {}'.format(order))\n", + "print('Approximation bias : {}'.format(func_x1-approx_x1))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that the further $p$ gets away from the point of the expansion (in that case $0$), the higher the approximation bias gets. Samely, the lower the order of approximation is, the higher the approximation bias gets." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.8.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/shearbook/intro/weak_lensing.ipynb b/shearbook/intro/weak_lensing.ipynb new file mode 100644 index 0000000..63ce796 --- /dev/null +++ b/shearbook/intro/weak_lensing.ipynb @@ -0,0 +1,98 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Weak Lensing\n", + "\n", + "A quick note on how the reduced shear is defined in Weak Lensing.\n", + "\n", + "![Weak Lensing Scheme](weak_lensing_scheme.jpeg)\n", + "\n", + "In the scheme above, let $\\Sigma$ be the surface mass density of the Dark Matter (L) then its lensing potential is given by:\n", + "\n", + "$$\n", + "\\Psi(\\theta) = \\frac{4G}{c^2}\n", + " \\frac{D_{l}D_{s}}{D_{ls}}\n", + " \\int \\Sigma{\\left(\\theta'\\right)} \\ln \\left\\|\n", + " \\theta - \\theta' \n", + " \\right\\| \n", + " d^2\\theta' \\quad,\n", + "$$ (lensing)\n", + "\n", + "where $G$ and $c$ are the gravitational constant and the speed of light in vacuum and $D_{l}$, $D_{s}$ and $D_{l,s}$ are the angular diameter distances between the observer and the lens (OL), the observer and the source (OS) and the lens and the source (LS), respectively.\n", + "\n", + "The deflection of the light rays is approximated well enough by the angle\n", + "\n", + "$$\n", + "\\alpha(\\theta) = \\nabla \\Psi(\\theta) \\quad.\n", + "$$\n", + "\n", + "This angle is related to the angular positions of the source $\\beta$ to the image $\\theta$ on the sky by the lens equation:\n", + "\n", + "$$\n", + "\\beta = \\theta - \\alpha(\\theta) \\quad.\n", + "$$ (lens)\n", + "\n", + "Under weak lensing assumptions (lens mapping almost constant along the solid angle of the soucre), the lens mapping is locally linear and the Jacobian matrix of the image distorsion is:\n", + "\n", + "$$\n", + "A \\equiv \\frac{\\partial \\beta}{\\partial \\theta} \n", + " = \\left( \\delta_{ij} - \\frac{\\partial^2 \\Psi(\\theta)}{\\partial \\theta_{i} \\partial \\theta_{j}} \\right)\n", + " = \\begin{pmatrix}\n", + " 1-\\kappa -\\gamma_1 & -\\gamma_2 \\\\\n", + " -\\gamma_2 & 1-\\kappa+\\gamma_! \n", + " \\end{pmatrix}\n", + "$$ (jacobian_lensing)\n", + "\n", + "with the convergence\n", + "\n", + "$$\n", + "\\kappa (\\theta) = \\frac{1}{2}\\left(\\Psi_{11} + \\Psi_{22}\\right)\n", + "$$\n", + "and the two components\n", + "\n", + "$$\n", + "\\gamma_1 = \\frac{1}{2} \\left(\\Psi_{11} - \\Psi_{22}\\right), \\quad \\gamma_2 = \\Psi_{12}\n", + "$$\n", + "\n", + "of the complex shear $\\gamma = \\gamma_1 + \\boldsymbol{i} \\gamma_2$.\n", + "\n", + "In the KSB method, the reduced shear is the measure of image distorsions and is defined as:\n", + "\n", + "$$\n", + "g = \\frac{\\gamma}{1-\\kappa}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.8.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/shearbook/intro/weak_lensing_scheme.jpeg b/shearbook/intro/weak_lensing_scheme.jpeg new file mode 100644 index 0000000..b395399 Binary files /dev/null and b/shearbook/intro/weak_lensing_scheme.jpeg differ diff --git a/shearbook/methods/ksb-variants.ipynb b/shearbook/methods/ksb-variants.ipynb deleted file mode 100644 index 7800f7e..0000000 --- a/shearbook/methods/ksb-variants.ipynb +++ /dev/null @@ -1,48 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# KSB Variants" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "```{warning}\n", - "In progress!\n", - "```" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "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.8.5" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/shearbook/methods/ksb.ipynb b/shearbook/methods/ksb.ipynb index d3e8a08..90e3dcc 100644 --- a/shearbook/methods/ksb.ipynb +++ b/shearbook/methods/ksb.ipynb @@ -4,24 +4,195 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# KSB" + "# KSB\n", + "\n", + "In this section, we introduce the Kaiser-Squire-Broadhurst method (KSB) to measure galaxies shape using the image of the galaxy surface brightness." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "KSB estimates the reduced shear $g$ using the polarisation $\\chi$. The exact relation between these quantities is :\n", + "\n", + "$$\n", + "\\chi^s = \\frac{\\chi - 2 g +g^2 \\chi^\\ast}{1 + | g |^2 - \\Re \\left( g \\chi^\\ast \\right) } \\quad,\n", + "$$(g_chi)\n", + "where $\\chi^s$ is the intrinsic (unlensed) polarisation.\n", + "\n", + "Assuming that the average of $\\chi^s$ vanishe, i.e. $\\left< \\chi^s \\right>=0$, and in a certain region $g$ is constant, $g$ is approximated by rotating the coordinate frame such that only one of its component does not vanish:\n", + "\n", + "$$\n", + "g \\approx \\frac{ \\left< \\chi \\right> }{ 2 \\left(1 - \\sigma^2_\\chi\\right) } \n", + " + \\frac{\\left< \\chi \\right>^3}{8} \\frac{1 - 5\\sigma^2_\\chi}{ \\left(1 - \\sigma^2_\\chi\\right)^4} \n", + " + \\mathcal{O} \\left( \\left< \\chi \\right>^5 \\right) \\quad,\n", + "$$\n", + "\n", + "with $\\sigma^2_\\chi$ being the standard deviation of $\\chi^s$ distribution.\n", + "\n", + "\n", + "```{note}\n", + "The polarisation $\\chi$ is used for shear measurement instead of the ellipticity $\\epsilon$ in Weak Lensing because in practice its estimator is less noisy. \n", + "```\n", + "\n", "```{warning}\n", - "In progress!\n", - "```" + "The approximation of $g$ is done by applying a Taylor Expansion of $\\chi$ around 0, i.e. for small value of $\\chi$, on the average of {eq}`g_chi`.\n", + "```\n", + "\n", + "## Standard KSB\n", + "\n", + "The core idea of KSB is to deduce the reduced shear by estimating the polarisation on raw data, i.e. images blurred by a known PSF and contaminated with additive Gaussian noise, and using a weighting window on each image then apply corrections directly in the polarisation space.\n", + "\n", + "In the following, the PSF effect and the weighting window effect are handled separately." ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "### Handling the Weighting Window\n", + "Here the polarisation $\\chi$ is estimated using the weighted moments $Q_w$ and $\\chi^s$ is the unaccessible polarisation of the source observed by a perfect telescope (no PSF effect) and without any noise.\n", + "\n", + "By introducing the shear polarisability matrix, $P^{sh}$, a first order approximation in terms of $g$ gives :\n", + "\n", + "$$\n", + "\\chi - \\chi^s = P^{sh} g\n", + "$$ (chichi)\n", + "\n", + "with\n", + "\n", + "$$\n", + "P^{sh}_{\\alpha \\beta} = - 2 \\frac{\\chi_{\\alpha}L_{\\beta}}{Q_{w,xx}+Q_{w,yy}} \n", + " - 2 \\chi_{\\alpha}\\chi_{\\beta}\n", + " + 2 \\frac{B_{\\alpha \\beta}}{Q_{w,xx}+Q_{w,yy}}\n", + " + 2 \\delta_{\\alpha \\beta}\n", + "$$\n", + "\n", + "such that $\\alpha, \\beta \\in \\{1,2\\}$ the indices of the components and $\\delta_{\\alpha \\beta}$ is the Kronecker delta and\n", + "\n", + "$$\n", + "L_\\beta = \\frac{1}{\\sigma^2} \\sum_{x,y} I(x,y) W'\\left(\\frac{(x,y)}{\\sigma}\\right) \\eta_\\beta\n", + "$$\n", + "\n", + "$$\n", + "B_{\\alpha \\beta} = \\frac{1}{\\sigma^2} \\sum_{x,y} I(x,y) W'\\left(\\frac{(x,y)}{\\sigma}\\right) \\eta_\\beta \\eta_\\beta\n", + "$$\n", + "\n", + "where $W$ is the weighting window function and $\\sigma$ its scaling factor along with\n", + "\n", + "$$\n", + "\\eta_1 = x^2 - y^2\n", + "$$\n", + "\n", + "and\n", + "\n", + "$$\n", + "\\eta_2 = 2xy.\n", + "$$\n", + "\n", + "By averaing equation {eq}`chichi`, $\\left<\\chi^s\\right> = 0$ (no intrinsic alignment assumption), $\\left = g$ (weak lensing assumption) and by multplying both sides by $\\left^{-1}$, we obtain:\n", + "\n", + "$$\n", + "g = \\left^{-1}\\left< \\chi \\right> \n", + "$$ (gpx)\n", + "\n", + "Equation {eq}`gpx` applies when the all the sources have locally constant morphology and size. In practice, this does not hold because the sources are affected by the PSF which adds instability to the morphology and size. Commonly in cosmic-shear, the averages are interchanged:\n", + "\n", + "$$\n", + "\\left< \\tilde{g} \\right> = \\left< \\left(P^{sh}\\right)^{-1} \\chi \\right> \n", + "$$ (gpx1)\n", + "\n", + "with the assumption that $\\left< \\left( P^{sh}\\right)^{-1} \\chi^s \\right> = 0$ which does not hold in reality because $P^{sh}$ depends on $\\chi$.\n", + "\n", + "Nevertheless, $\\tilde{g}$ is the solution to {eq}`chichi` with $\\chi^s=0$. This corresponds to the shear of a circular source. The tru value of $g$ is obtained by averaging $\\tilde{g}$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Handling the PSF Effect\n", + "\n", + "In KSB the PSF is supposed to be almost istropic and its effects on the source polarisation are linearly modelled by:\n", + "\n", + "$$\n", + "\\chi = \\chi^{obs} + \\chi^{g}\n", + "$$\n", + "\n", + "where $\\chi$ is the polarisation that would have been observed without the PSF effect (it is the same $\\chi$ from the section above), $\\chi^{obs}$ is the polarisation observed by the telescope and $\\chi^g$ contains the effects of the PSF which expression is:\n", + "\n", + "$$\n", + "\\chi^g = P^{sm} \\left( P^{sm,*}\\right)^{-1} \\chi^{sh,*}\n", + "$$\n", + "\n", + "$P^{sm}$ is the smear polarizability matrix and it accounts for the effects of the isotrpic part of the PSF, $P^{sm,*}$ and $\\chi^{sh,*}$ are respectively the smear polarizability and the polarization of a star and are used to account for the anisotropies of the PSF that are supposed to be small. The star is supposed to give an empirical estimation of the PSF near the source. The expression of the smear polarizabilty is:\n", + "\n", + "$$\n", + "P^{sm}_{\\alpha \\beta}=\\frac{1}{Q_{xx}+Q_{yy}}\\left[ \n", + " \\left( M + 2\\frac{Q'_{xx}+Q'_{yy}}{\\sigma^2}\n", + " \\right) \\delta_{\\alpha \\beta}\n", + " + G_{\\alpha \\beta}\n", + " -\\chi_\\alpha \\left( 2 F_\\beta + L'_\\beta \\right)\n", + " \\right]\n", + "$$\n", + "\n", + "where $L'_\\alpha$ is equivalent to $L_\\alpha$ calculated with the second order derivative of the weighting window function, $Q'_{k,k}$ (with $k \\in {x,y}$) is equivalent to the moments weighted by the first order derivative of the weighting window function and\n", + "\n", + "$$\n", + "M = \\sum_{x,y} I(x,y) W \\left( \\frac{\\left(x,y\\right)}{\\sigma} \\right) \\quad,\n", + "$$\n", + "\n", + "$$\n", + "F_\\alpha = \\frac{1}{\\sigma^2}\\sum_{x,y} I(x,y) W' \\left( \\frac{\\left(x,y\\right)}{\\sigma} \\right) \\eta_\\alpha \\quad,\n", + "$$\n", + "\n", + "$$\n", + "G_{\\alpha \\beta} = \\frac{1}{\\sigma^4} \\sum_{x,y} I(x,y) W'' \\left( \\frac{\\left(x,y\\right)}{\\sigma} \\right) \\eta_\\alpha \\eta_\\beta \\quad.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### KSB in practice\n", + "\n", + "The reduced shear is estimated by removing the PSF effect then solving equation {eq}`gpx1`.\n", + "\n", + "The original KSB {cite}`1995ApJ.KSB` had an incorrect approximation of $g$. Many ameliorations were suggested, namely from {cite}`2011MNRAS.410.2156V` which pushed the approximations up to the third order. Alternative methods to KSB exist, for example there is the re-Gaussianization approach which is used in HSC {cite}`2003MNRAS.343..459H`.\n", + "\n", + "First order approximation and limited Taylor expansion induce a large bias for the points that are distant from the expansion point, for more details see the [appendix on Taylor series](../appendix/taylor_series.ipynb). This means that for high values of $g$ the error is more important than for smaller values, as seen in the figure below.\n", + "\n", + "\n", + "```{figure} shear_estimation.png\n", + "---\n", + "width: 650px\n", + "name: ksb_shear_est\n", + "---\n", + "Shear estimation for unconvolved and almost noiseless galaxies with Sersic profile. KSB (red solidline) goes to third order approximations which are incorrect, KSB1 (violet dotted) is a first order approximation, KSBtr (green long-dashed line) is a third order approximation that is not mathematically justified and KSB3 (blue short-dashed line) is the correct third order approximation. Source {cite}`2011MNRAS.410.2156V`.\n", + "```\n", + "\n", + "Also KSB works under the hypothesis that the PSF has small anisotropies which is not always the case, specially in Radio Interferometry for example.\n", + "\n", + "Furthermore, the correction of the weighting window effect is approximative which adds another bias to the estimation of $g$ that depends on the choice of the window parameters (e.g. $\\sigma$).\n", + "\n", + "Finally, a C++ implementation of KSB by Peter Melchior can be found on this [GitHub repository](https://github.com/pmelchior/shapelens/blob/master/src/KSB.cc). And in the following code snippet, an example of how to use the GalSim KSB implementation to estimate the shear of a galaxy contained in `image` which was convolved by `psf` (assuming that both `image` and `psf` are numpy arrays):\n", + "\n", + "```python\n", + "import galsim\n", + "from galsim import Image\n", + "import galsim.hsm\n", + "\n", + "#set pixel scale\n", + "scale = 0.1\n", + "#create a galsim version of the data\n", + "image_galsim = Image(image,scale=scale)\n", + "psf_galsim = Image(psf,scale=scale)\n", + "#estimate the moments of the observation image\n", + "ell=galsim.hsm.EstimateShear(image_galsim\n", + " ,psf_galsim,shear_est='KSB')\n", + "```" + ] } ], "metadata": { @@ -40,7 +211,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.8.6" } }, "nbformat": 4, diff --git a/shearbook/methods/shear_estimation.png b/shearbook/methods/shear_estimation.png new file mode 100644 index 0000000..30bc930 Binary files /dev/null and b/shearbook/methods/shear_estimation.png differ diff --git a/shearbook/methods/viola2011.md b/shearbook/methods/viola2011.md deleted file mode 100644 index 7554207..0000000 --- a/shearbook/methods/viola2011.md +++ /dev/null @@ -1,73 +0,0 @@ -# Biases in, and corrections to, KSB shear measurements. Viola et al, 2010 - -```{warning} -In progress! -``` - -Notes and tools derived from reading groups sessions - -{cite}`2011MNRAS.410.2156V` - -## Goals of the Paper - -1. Show the theoretical and practical caveats of KSB. -2. Suggest improvements to KSB. - -## Main Results and Conclusion - -1. The assumptions in KSB introduce biaises especially for large amplitude -values. -2. Rewriting correctly the Taylor expansion of the ellipticity up to the third -order (see KSB3) gives results with biaises smaller by orders of magnitude. -3. The PSF corrections in KSB are not accurate however they partially cancel -other biaises. - -## Definitions and Equations - -*. Reduced shear, g: (eq. 7) -*. Complex ellipticity, $\chi$: (eq. 10) -*. Shear responsivity $\epsilon$: (eq. 14) -*. KSB (eq. 32) -*. KSB1 (eq. 34) -*. KSBtr (eq.36) -*. KSB3 (eq. 44) -*. Inverse Problem (eq. 50) -*. Circular PSF Formula (eq. 58) -*. Anisotropic PSf Formula (eq. 67) - -> The reduced shear, g, is preferred over the shear responsivity, -$\epsilon$, for its robustness to noise in the data. - -## Assumptions in KSB -1. The reduced shear, g, is constant and the ellipticity average is null, -$\left< \chi \right>=0$ (see eq. 12). -2. Many Taylor expansions like in the reduced shear, g (see eq. 29). -> Taylor expansions gives a polynomial approximation around a certain reference -value. The further from the reference we evaluate the expansion the bigger the -biais. In this paper, since the reference value is zero, the biais will be -greater for large amplitude values. - -3. Assume equality between (eq. 27) and (eq. 28). -4. The PSF has at most slight anisotropies, it is mostly circular. - -## Notations Clarification - -1. The angular coordinates are solid angles so they have two components. These -coordinates can be approximated by linear coordinates because the angles are -very small. -2. The star in eq. 12 stands for the conjugate. -3. Rewriting eq.48-50: - i. (eq. 48) $$\chi = f\left( \tilde{g}, \chi^{s}\right)$$ - ii. (eq. 49) $$\chi = h \left( \chi^{obs}\right)$$ - iii. (eq. 50) $$\chi^{obs} = h^{-1}\left[ f \left( \tilde{g}, h\left(\chi^{obs}\right)\right)\right]$$ - - where $\chi^{s}$ is the intrinsic ellipticity of the galaxy, $\chi^{obs}$ is - the ellipticity of the lensed galaxy and $\chi$ is the ellipticity of the - convoluted and lensed galaxy. -4. In (eq. 51): $$\chi^{sh}_{\alpha} = \chi_{\alpha} + \chi^{s}_{\alpha}$$ -5. In (eq. 52): the asterisk, *, in $P^{sm,*}$ is to indicate the PSF. -Explanation: * --> star --> dirac --> PSF --> illuminati - -## Implementation - -[Link](https://github.com/pmelchior/shapelens/blob/master/src/KSB.cc) to Peter Melchior implementation of KSB codes in C. diff --git a/shearbook/moments/moments-basic.ipynb b/shearbook/moments/moments-basic.ipynb index 28ddc82..c25594a 100644 --- a/shearbook/moments/moments-basic.ipynb +++ b/shearbook/moments/moments-basic.ipynb @@ -221,9 +221,9 @@ "Finally, we can define the second moment or quadrupole moments of an image as follows.\n", "\n", "$$\n", - " Q_{xy} = \\frac{\\sum_{x,y} I(x,y)(x - \\bar{x})(y - \\bar{y})}{\\sum_{x,y} I(x,y)}\n", + " Q_{k_1 k_2} = \\frac{\\sum_{x,y} I(x,y)(k_1 - \\bar{k_1})(k_2 - \\bar{k_2})}{\\sum_{x,y} I(x,y)}\n", "$$ (moments)\n", - "\n", + "where $k_1,k_2 \\in \\{x,y\\}$, the quadrupoles are $Q_{xx},Q_{xy}$and $Q_{yy}$ (with $Q_{xy}=Q_{yx}$).\n", "Equation {eq}`moments` is used to determine object shapes. In the following cell we a simple implementation of this equation." ] }, @@ -488,7 +488,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.8.6" }, "toc": { "base_numbering": 1, diff --git a/shearbook/moments/moments-weighted.ipynb b/shearbook/moments/moments-weighted.ipynb index d6a338e..4f6d1d2 100644 --- a/shearbook/moments/moments-weighted.ipynb +++ b/shearbook/moments/moments-weighted.ipynb @@ -5,24 +5,269 @@ "metadata": {}, "source": [ "(content:weighted)=\n", - "# Weighted Moments" + "# Weighted Moments\n", + "\n", + "Observed galaxy images are corrupted by noise. Straightforward measures of ellipticity on such images lack robustness, see figure below for example." + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ellipticity of the clean star image : -0.0040+i0.1980\n", + "Ellipticity of the noisy star image : -0.0050+i0.2190\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAF2CAYAAABzg27uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUnUlEQVR4nO3de7CtdVkH8O+DiCmCKJSGF0zNSi0zU8NSmswLXqkZp9QyTUuk0vLClBWS17Kb97zUJCNpWBqTQmVNYVmS2uQFvGSKCAnIRRBBC+TXH+978jm7vffZ57L32vvsz2fmDGut9/Zb7xme9T2/91nvqjFGAACAyQGLHgAAAGwmAjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyK6qqk6vq1EWPYyNV1ROq6t2LHgfAWlXVX1XVTy16HPtKVb2uqn590eNgexOQt7mqenxVfbCqvlxVF82F9gcWPa7lVNWTquq963mMMcafjDEesp7HAOiq6rNVdUlVHdxee2pVnbWW7ccYx44xTlm3ATZV9aaqetF6HmOMcfwY44XreQzYFQF5G6uqZyV5eZKXJLl1kjskeW2SxyxwWOumqg5c9BgAVnBgkmcuehDrTR1mqxCQt6mqukWSFyT5uTHGO8YY14wxrhtjvHOM8dwVtvm+qvqXqrqyqj5cVT/Ylj25qj5eVVdX1Weq6mlt2Q9W1YVV9eyq+sI8U/3kVcb2pHkfV1fVeXPbw3ckeV2So+fZ7ivndR9RVf9eVV+qqguq6uS2nztW1aiqp1TV55L8/RrOy06z1PP2J1TVp+bxvLCq7lxV75uP+baqOmhe95ZV9a6qurSqvjg/vl3b17dU1T/O+/m7qnpNb2FZ7fwC+73fTvKcqjpsuYVVdf+q+kBVXTX/9/5t2VlV9dT58V2q6j3zepdV1Wnz66+pqt9dss93VtUvLnOsqqrfn+v1VVX1kaq6R1X9bJInJDlxrsPvnNf/5ar69FzbPlZVP9L29aSq+ud5f1ckOXlXJ6LPUrfPjxPb58dxVfXwqvqPqrqiqp7Xtr3vXJ+vnNd99Y4aPS9/SFV9cn5fr53P1VPb8p+eP8u+WFV/U1VH7Wq87J8E5O3r6CTfkOQv1rJyVd02yRlJXpTkVkmek+TtVfWN8ypfSPLIJIcmeXKS36+q72m7uE2SWyS5bZKnJHlNVd1ymeMcnOSVSY4dYxyS5P5JPjTG+HiS45O8b4xx8zHGYfMm1yR5YpLDkjwiydOr6rgluz0myXckeeha3usyHpbk3km+L8mJSd6Q6UPi9knukeRx83oHJPnjJEdlmo3/SpJXt/28Jcn7kxye6UPiJ9v73tX5BfZvH0xyVqb/93dSVbfKVB9emal+/F6SM6rq8GX288Ik705yyyS3S/Kq+fVTkjyuqg6Y93lEkgcleesy+3hIkgcmuWum2vpjSS4fY7whyZ8kedlchx81r//pJA/IVON/I8mpVfXNbX/3S/KZJN+U5MW7OA/LuU2mz6vbJjkpyRuT/ESmuvyAJCdV1Z3mdb+W5JeSHJHpc+5BSU5o7/nPk/xKpvP4yUyfMZmXH5fkeUl+NMk3JvmnLH9+2AYE5O3r8CSXjTGuX+P6P5HkzDHGmWOMG8YYf5upoD88ScYYZ4wxPj0m78lUoB/Qtr8uyQvmWeozk3w5ybetcKwbktyjqm46xrhojHHuSoMaY5w1xvjoPKaPZCpmxyxZ7eR5hvwra3yvS/3WGONL8zjOSfLuMcZnxhhXJfmrJPeax3L5GOPtY4xrxxhXZ/ogOCZJquoOSe6T5KQxxv+MMd6b5C/bMVY9v8C2cFKSX1jmH8aPSPKpMcabxxjXjzHemuQTSR71//Yw1dqjkhw5xvjqXGsyxnh/kqsyBcYk+fEkZ40xLllhH4ck+fYkNcb4+BjjopUGPcb4szHG5+fadVqSTyW5b1vl82OMV81j35M6fF2SF48xrkvyp5nC7yvGGFfPdfncJN81j+Xfxhhnz8f6bJLX5+ufCQ9Pcu581fT6TP/guLgd52lJXjq/3+sztR9+t1nk7UlA3r4uT3JErb0f7Kgkj50vW11ZU4vDDyT55iSpqmOr6uz5cteVmQrREf14S8L4tUluvvQgY4xrMs1WHJ/koqo6o6q+faVBVdX9quof5raGq+btjliy2gVrfI8r6R8gX1nm+c3nsdysql5fVedX1ZeS/GOSw6rqRkmOTHLFGOPaFca16vkF9n9jjHOSvCvJLy9ZdGSS85e8dn6mGdWlTkxSSd5fVedW1U+3Zadk+sd45v++eYVx/H2mq1+vSXJJVb2hqg5dadxV9cSq+lCrXffIznV4b2vw5WOMr82PdwTslerwXWtqb7t4rsMvaWM5so9ljDGSXNj2c1SSV7T3cUWmc7nceWY/JyBvX+9L8tUkx61x/QuSvHmMcVj7c/AY4zer6iZJ3p7kd5Lcem5/ODNTYdltY4y/GWM8OFM4/ESmy2lJMpZZ/S2ZZmJvP8a4RaY+5aXHXW679fDsTLPi9xtjHJrpEmXm8VyU5FZVdbO2/u3b4xXP74aMHNgsnp/kZ7JzKPt8pvDW3SHJfy3deIxx8RjjZ8YYR2aaEX1tVd1lXnxqksdU1T0ztZ2dvtIgxhivHGPcO8ndM7Va7Phuyk71dJ5dfWOSn09y+Fz/z8nOdXijanCS/EGmz41vnevw89pYLsrUdpJk6rXuzzPV4actqcM3HWP8ywaNnU1EQN6m5vaAkzL1Ah83z37eeJ4Jftkym5ya5FFV9dCqulFVfcP85YnbJTkoyU2SXJrk+qo6NlMP226rqltX1aPnXuT/ztSKsWPm4JIkt+tfuMh0GfCKMcZXq+q+SR6/J8fdRw7JNJNx5dwz+PwdC8YY52dqmTi5qg6qqqOz8+XR1c4vsE2MMf4zyWlJntFePjPJXWu6LeeBVfVjSe6WabZ5J1X12FY3vpgpnH5t3veFST6Qaeb47Su1O1TVfearczfO9D2Pr2bnOnyntvrB8zEunbd9cqYZ5EU5JMmXknx5vvr49LbsjCTfOX/mHZjk5zL1N+/wuiS/UlV3T6Yvs1fVYzdo3GwyAvI2Nsb4vSTPSvJrmYrbBZlmAU5fZt0LMt3+7Xlt3ecmOWDut31GkrdlKsiPz879tbvjgEwzsZ/PdHnrmMxfsMh0F4pzk1xcVZfNr52Q5AVVdXWmwP+2PTzuvvDyJDdNclmSs5P89ZLlT8j0pZHLM30Z77RM/whY9fxuwLiBzeUFmYJnkun7DZm+BP3sTPXjxCSPHGNctsy290nyr1X15Ux1+JljjPPa8lOSfGdWaK+YHZppVviLmVo5Ls90hTBJ/ijJ3eY2hNPHGB9L8ruZrkpeMu/7n3fv7e5Tz8n0GXR1pvdw2o4F8/l6bJKXZXpPd8s0cbGjDv9Fkt9K8qdze8Y5SY7dyMGzedTUggNstJpuv/SJMcbzd7kywD5QVQ/MdMXqjmOMGxY9nkWa7+hxYZInjDH+YdHjYXMxOwUbZL5seeeqOqCqHpZpxvj0BQ8L2CbmlolnJvnD7RqO5za2w+bvzuzoTz57wcNiE/KLNrBxbpPkHZlusXdhkqePMf59sUMCtoOafmzpg0k+nOle9dvV0Zm+3H1Qko8lOW4vbgHKfkyLBQAANFosAACgEZABAKBZtQe5qvRfAKyDMcaaf0hHLQZYHyvVYjPIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQCMgAANAIyAAA0AjIAADQHLjoAcBmM8b4v8dVtcCRAGxfajGLZAYZAAAaARkAABoBGQAAGj3IsIreA5fogwNYBLWYjWYGGQAAGgEZAAAaARkAABo9yGx7S3vbdmddfXAA+4ZazGZiBhkAABoBGQAAGi0WbEu7cylvrftxiQ9g96jFbFZmkAEAoBGQAQCgEZABAKDRg8y2sK/63HbnGPrgAHamFrNVmEEGAIBGQAYAgEaLBfutjbiUt6fHd8kP2C7UYrYiM8gAANAIyAAA0AjIAADQ6EFm4Rbdn7YIG/Ge9dYBu0MtXh9q8dZkBhkAABoBGQAAGgEZAAAaARkAABoBGQAAGgEZAAAaARkAABoBGQAAGgEZAAAaARkAABoBGQAAGgEZAAAaARkAAJoDFz0AYH2MMXZ6XlULGgnA9qUWb01mkAEAoBGQAQCgEZABAKDRg8zCLe3HWtqvBcD6U4vh68wgAwBAIyADAEAjIAMAQCMgAwBAIyADAEAjIAMAQOM2byycWwltjH6e/dQpsJRavDHU4q3BDDIAADQCMgAANAIyAAA0epDZcPrcFm/p34E+ONh+1OLFU4s3LzPIAADQCMgAANBosWBDuJS3ua329+OSH+w/1OLNTS3ePMwgAwBAIyADAEAjIAMAQKMHGViV2xABLJ5avLHMIAMAQCMgAwBAo8WCdeFWQvuv/nfrEh9sbmrx/kstXl9mkAEAoBGQAQCgEZABAKDRg8w+o9dt+3HbIdh81OLtRy3e98wgAwBAIyADAEAjIAMAQKMHmX2m9zzpgdse9LnB5qMWbz9q8b5nBhkAABoBGQAAGi0WrIull3tc5gPYeGox7BkzyAAA0AjIAADQCMgAANDoQWZDrHYLGj1xW4vbCcHWpRbvP9Ti9WUGGQAAGgEZAAAaLRYsnNsQbW4u48H2oBZvbmrxxjKDDAAAjYAMAACNgAwAAI0eZDad3melB27j6XMDErV40dTixTKDDAAAjYAMAACNFgs2Nbcd2hgu5QGrUYs3hlq8eZhBBgCARkAGAIBGQAYAgEZABgCARkAGAIBGQAYAgEZABgCAxn2Q2VL89On66OfSfTiBXVGL14davHmYQQYAgEZABgCARkAGAIBGQAYAgEZABgCARkAGAIBGQAYAgEZABgCARkAGAIBGQAYAgEZABgCARkAGAIBGQAYAgEZABgCARkAGAIBGQAYAgEZABgCA5sBFDwD2VFVt+DHHGBt+zNUs4hwAdGqxWrw/MoMMAACNgAwAAI0WC9gNSy+jbfRlPpfxANRi1p8ZZAAAaARkAABoBGQAAGj0IMNe6H1o69UDp9cNYHVqMfuaGWQAAGgEZAAAaLRYwD6yr2475DIewJ5Ti9kXzCADAEAjIAMAQCMgAwBAowcZ1slG3HYIgNWpxewJM8gAANAIyAAA0GixgA3gdkEAi6cWs1ZmkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKARkAEAoBGQAQCgOXDRAwAAgL11pyTHJ/neJLdIcm2S85L8UZL37Oa+aoyx8sKqlRcCsMfGGLXWddVigJXdM8lLkxy7yjofS/KiJG9d8vpKtVhABlgAARlg7z04yTuS3HyN678oya+35wIywCYiIAPsnfskOSvJzebnNyQ5M8kpST6bqc3iR5I8MckhbbvnJvmd+bGADLCJCMgAe+ecJHefH5+f5DFJPrzMeoclOTXJI+bnNyS5S6b+5JVqsbtYAACwpTwwXw/H12ZqtVguHCfJlUl+NMkH5ucHJHnaLvYvIAMAsKWc0B6/KcmndrH+/yT5jfb8KUlussr6WiwAFkCLBcCeuzzJrebH90ryobZsabatmsrtAUk+l+S28+tHJ3mfFgsAALa6ytRXvMNH17jdDZn6lne45SrrCsgAAGwZI8nX2vODdmPb3lZx3SrrCcgAAGwpF7fHP7zGbQ7NdGu45faxlIAMAMCW8pb2+IQly6pqpz87/FSSg+fHH83O7RZLCcgAAGwpr8/UU5wkD0vyuF2sf+ckz2/PX7uL9QVkAAC2lPOSnN6evznJr2Zqo+hulOkeyO9Ncvj82qWZfjhkNW7zBrAAbvMGsHcOyxR8795euyZTcD5vXv7oJHdoy7+S5IeSnD0/91PTAJuIgAyw974pybuy85fvVnJFkuOS/FN7zU9NAwCwX/lCkmOSHJ/kIyusc2mSlya5Z3YOx6sxgwywAGaQAfa9709y70ztFdck+WySd2b6qenl7FGLBQAAbDdaLAAAoBGQAQCgEZABAKARkAEAoBGQAQCgEZABAKD5X2+LAn96yURFAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The noisy star and the clean star images are binary images.\n", + "They differ from each other by one pixel on the bottom right of the noisy star image.\n", + "This pixel induces a difference in ellipticity estimation of : 10.6% between the two images\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "def get_centroid(data):\n", + " \n", + " # Sum flux over x and y individually\n", + " sum_i = np.array([np.sum(data, axis=i) for i in (1, 0)])\n", + " \n", + " # Get range of x and y values\n", + " ranges = np.array([np.arange(i) for i in data.shape])\n", + " \n", + " # Calculate centroids\n", + " cents = np.sum(sum_i * ranges, axis=1) / np.sum(data)\n", + " \n", + " return cents.astype(int)\n", + "\n", + "def get_moments(data):\n", + " \n", + " centroid = get_centroid(data)\n", + " ranges = np.array([np.arange(i) for i in data.shape])\n", + " \n", + " x = np.outer(ranges[0] - centroid[0], np.ones(data.shape[1]))\n", + " y = np.outer(np.ones(data.shape[0]), ranges[1] - centroid[1])\n", + " \n", + " q = np.array([np.sum(data * xi * xj) for xi in (x, y) for xj in (x, y)])\n", + " q = (q / np.sum(data)).reshape(2, 2).astype('complex')\n", + " \n", + " return q\n", + "\n", + "def get_ellipticity(data, method='chi'):\n", + " \n", + " # Calculate moments\n", + " q = get_moments(data)\n", + " \n", + " # Calculate the image size.\n", + " r2 = q[0, 0] + q[1, 1]\n", + "\n", + " # Calculate the numerator\n", + " num = (q[0, 0] - q[1, 1] + 2 * np.complex(0, q[0, 1]))\n", + " \n", + " # Calculate the denominator\n", + " den = r2\n", + " \n", + " if method == 'epsilon':\n", + " den += 2 * np.sqrt(q[0, 0] * q[1, 1] - q[0, 1] ** 2)\n", + " \n", + " # Calculate the ellipticity/polarisation\n", + " ellip = num / den\n", + "\n", + " return np.around([ellip.real, ellip.imag], 3)\n", + "\n", + "#load toy image and generate noisy image\n", + "star_clean = np.load('star.npy')\n", + "n_row,n_col = star_clean.shape\n", + "star_noisy = np.copy(star_clean)\n", + "noise_coord = (n_row-3,n_col-3)\n", + "star_noisy[noise_coord] = 1\n", + "\n", + "#estimate ellipticities\n", + "e_clean = get_ellipticity(star_clean)\n", + "e_noisy = get_ellipticity(star_noisy)\n", + "print('Ellipticity of the clean star image : {:.4f}+i{:.4f}'.format(*e_clean))\n", + "print('Ellipticity of the noisy star image : {:.4f}+i{:.4f}'.format(*e_noisy))\n", + "\n", + "#compute deviation\n", + "e_clean_norm = np.linalg.norm(e_clean)\n", + "e_noisy_norm = np.linalg.norm(e_noisy)\n", + "deviation = np.abs((e_noisy_norm-e_clean_norm)/e_clean_norm)*100\n", + "\n", + "#show image\n", + "cmap = 'bone'\n", + "plt.figure(3,figsize=(10,7))\n", + "plt.subplot(121)\n", + "plt.imshow(star_clean,cmap=cmap)\n", + "plt.axis('off')\n", + "plt.title('Clean star image')\n", + "plt.subplot(122)\n", + "plt.imshow(star_noisy,cmap=cmap)\n", + "plt.axis('off')\n", + "plt.title('Noisy star image')\n", + "circle = plt.Circle(noise_coord, 2, color='r', lw = 3, fill=False)\n", + "plt.gcf().gca().add_artist(circle)\n", + "plt.tight_layout()\n", + "plt.show()\n", + "print('The noisy star and the clean star images are binary images.')\n", + "print('They differ from each other by one pixel on the bottom right of the noisy star image.')\n", + "print('This pixel induces a difference in ellipticity estimation of : {:.1f}% between the two images'.format(deviation))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "```{warning}\n", - "In progress!\n", + "One way to have more robust measurements is to reduce the weight of the pixels containing only noise in the image before computing the moments. For this purpose we use a weighting window $W$ to define the weighted quadrupole moments of an image as follows:\n", + "\n", + "$$\n", + " Q_{w,k_1 k_2} = \\frac{\\sum_{x,y} W(x,y) I(x,y) (k_1 - \\bar{k_1})(k_2 - \\bar{k_2})}{\\sum_{x,y} W(x,y) I(x,y)}\n", + "$$ (weighted_moments)\n", + "with $k_1,k_2 \\in \\{x,y\\}$.\n", + "## Weight dependent measures\n", + "\n", + "In equation {eq}`weighted_moments` the output varies with $W$ which means that the resulting shape estimation also varies depending on $W$ see example below. Usually the window is chosen to be the best fit of a 2D Gaussian on the galaxy image. The expression of 2D isotropic Gaussian window centered in 0 is :\n", + "$$\n", + "W(p) = \\frac{1}{2\\pi} \\exp\\left(-\\frac{\\|p\\|^2}{2}\\right)\n", + "$$\n", + "with $p=(x,y)$. Once the centroid, $\\bar{p}=\\left(\\bar{x},\\bar{y}\\right)$ and the scaling factor $\\sigma$ are defined for a specific image, the window is translated and rescaled as :\n", + "$$\n", + "W\\left(\\frac{p-\\bar{p}}{\\sigma}\\right) \\quad.\n", + "$$\n", + "```{note}\n", + "In the case of isotropic windows, only the norm of the point considered is required and the value of the window can be noted as $W(\\|p\\|^2)$. \n", "```" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAF2CAYAAABzg27uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsiElEQVR4nO3debBm+VkX8Od3957unumZbJAEEk1AiCRAFEF2S5ZgQHEpRUUJGNlUqFIQLBUCihGLkiigUKlSlIAxLlC4IKCSFaKICVFkSwIhMGGSzExPb7fvevzjnO73Oc/b77lL7z2fT9XUvOf+3vue857u/t3nnvN9n1/rui4AAIDe0u0+AAAAuJMokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUy19Rae25rrWutrQzbr2+tvfx2H9dhtdb+bGvtJ+7m46h/BsC9o7X2aa21X76O7+9aa8+/kcd0xP3/emvts27X/uFmUyDfZVprn9pa++nW2hOttcdaa29prX3CHXBcH9Fae21r7QOttXOttV9trX1Xa+3Zt+N4uq77wa7rPudGv25r7ftaa/8kba+21i4u+Non3azjAO4srbW/0Vr7z+Vrv7rga1/Udd2buq77Xbf2KG+N1tqzW2v/rrX2weFn1f9prb1sGLspv/jf7l8YuPcokO8irbX7I+I/RsR3RcRDEfGsiPiWiNi6zcf1/Ij4HxHxcER8fNd190fEp0TEuyLiU2/nsd0Eb4yIz0jbvzcifiMiPr18LSLi527VQQG33Rsj4lNaa8sREa21D4mI1Yh4cfna84fn3st+ICLeGxHPiYinRMSfj4hHbsQLu6PGraJAvrt8ZERE13X/quu6va7rNruu+4mu694REdFae9lwRfk7W2tnW2vvbq198vD197bW3t9a+5IrL9Zae2lr7W3DFd/3ttZecczjekVEvKXrur/add1vDsf4/q7rXtV13WuHfT3YWvuPwxXmx4fHV68u19t1rbVXtNZeMzzeaK29prX26PC+fra19oz0nt/dWjvfWvu11tqfTV9/c3q9fzS8x3OttZ9rrX1a2dfrWmv/cnidX2itXSlyqzdExEe31p46bH9aRLw2Ik6Wr/1M13U71ziOrrX2lcNVpMdba9/TWmvD2HJr7TuGqy7vjoiX5h231p7ZWvvR4c7BO1trfzGdn80r+2+t/a3W2u7wC1W01v5ua+1VB/4pAtfjZ6MviD9u2P70iPipiPjl8rV3dV33cGvtM1trv3nlm4c58Otaa+8Yrrr+69baRhr/+tba+1prD7fWvizvuLX2wDB/faC19p5hDlgaxt7TWvs9w+MvHuagFwzbL2+t/cjweKm19o2ttXcNc+3rWmsPpX38ueG1Hm2t/c0DzsUnRMT3d113seu63a7r3tZ13Y8NY1d+OTjbWrvQWvv9rbXntdb++/DaH2yt/WBr7Uw5N9/QWntHRFxsBxTJw5z+b4afG+dbfwX7I1t/lf/9w8+Cz0nP/9LW2i8Oz313a+0ryuv99XTuX97S1erW2vowb/9Ga+2R1tr3ttZOHHB+uAsokO8uvxIRe621f9Fa+7zW2oPXeM4nRsQ7ov+t/YeiL94+IfqrFl8cEd/dWjs1PPdi9L/Zn4m+GPuq1toXHuO4Pisi/t0Bz1mKiH8e/RWFD4+IzYj47kO+/pdExAMR8WHRv6+vjIjN1trJiPjHEfF5XdedjohPjoi3L3iNn43+h9RD0Z+Xf5N/+ETEH47+XJ2JiB9ddGzDLwDvib4Ijuh/4L0pIn66fG3qCtHnR/9n8rER8Scj4nOHr//FYezjo78K/SfK9/2riPjNiHjmMPb3Wmt/sOu6y8P7u3Jl+9OHY/yUtP2GieMBrlPXddvR30m7cjfpytzw5vK1qbnhT0bESyLid0TEiyLiZRERrbWXRMTXRcRnR8RHRD/nZt8V/Rz5O6OfB/58RHzpMPaGiPjMtP93x3iuuDI3fE1EfOEw9syIeDwivmfY/wsi4p9GxJ8bxp4SEVPxubdGxPe01r6otfbhZezKuTjTdd2prut+JiJaRLxyeO2Pjn6uf0X5vj8d/c+pM13X7U7s+4oviP5K9oMR8baI+PHofw49KyK+NSK+Lz33/dHPvfdHf96+s7X24uG9vyQi/mr05/z5Mb6DGBHx7dFfvPq4YfxZEfFNhzg+7nRd1/nvLvov+snj+6MvlHajL+aeMYy9LCJ+NT33hRHRXRkfvvZoRHzcgtd+VUR85/D4ucP3rgzbr4+Ily/4vt2IeEna/ssRcTYiLkTEqxd8z8dFxONp+9cj4rPS9isi4jXD4y+LvgB9UXmNk8N+/nhEnChjL4uIN0+cx8cj4mPTvv5rGntBRGxOfO/3R8R3Rj/Zvj8i7ou+aL/ytccj4jOudRzDOf3UtP26iPjG4fF/j4ivTGOfc+XPIPofGHsRcTqNvzL6qzQREX8n+l8WViLityPiayPi70fERvS/jDz1dv/d9Z//7vX/hrnkh4fHPx99MfuS8rUvGR5/ZkT8ZvreX4+IL07b/yAivnd4/M8i4u+nsY8c5obnR8Ry9DG7F6Txr4iI1w+P/0JE/Ojw+Bcj4uUR8dph+z0R8eI09gfTa3xoROwMc8o3XfmeYexkRGxHmrPLeXhwmH9+YZi33h4RnzCMPTfSz5YF3/+FEfG2cm6+7IBz30XE89Ofw0+msS+I/ufR8rB9enj+mQWv9SMR8bXp3L8yjT0/nfsW/YWm56Xx3x8Rv3a7/y767/r/cwX5LtN13S92XfeyruueHREfE/1v3K9KT8k5r83he+rXTkVEtNY+sbX2U8NtuSeiL/KeGkf3aPST6ZVj/O6u684Mx7U67Ou+1n/A7T2ttXPRX0U504Zs3gF+IPrf/l873OL6B6211a7rLkbEnxqO+32ttf/UWvuoa71Aa+2vDbfQnmitnY3+akt+r7+dHl+KiI2J23hvjP4qyAsj4t1d112K2VWiF0bEieivJC1S93Xliv4zo8/tXfGe9PiZEfFY13Xny/izhsdXrhK9OCL+T0T8ZPRXOj4pIt7Zdd0HJ44HuDHeGBGfOtzde1rXdb8a/S/3nzx87WNi+gryceaGp0bEWvlanRs+rfX55+WI+NfRZ6WfG/08+Pbhec+JiB9ufYztbPQF815EPKPuf5h7H130Jrque7zrum/suu53D9//9oj4kdb6OFnVWnt66z/k/VvDz4fXxPzPovde41un1J97H+y6bi9tR8x+Fn5ea+2trY+vnY2IP5T2X899fvy06C+Q/Fw6b/9l+Dp3OQXyXazrul+K/mrmxxzzJX4o+ivQH9Z13QMR8b3R/0Z8VP8tIv7YAc/5axHxuyLiE7v+Q3xXbrNd2d/F6CeaKz7kyoOu63a6rvuWruteEH2M4vOjv4UYXdf9eNd1nx19gf5LEfHquuPW542/Ifrblw8OxfsTcbz3GtH/gPvY6G/3vWn42i9Ef5X3pRHxs10feziq9w2vcUW+NflwRDzUWjtdxn9rePzT0Z/fPxoRb+i67v8N4y8N8Qq4VX4m+qLzyyPiLRERXdedi/7f75dHxMNd1/3aMV53am74YPRXep9Txn9r2P87oy+2vyYi3jj8kv3bw/G8ueu6/eF73ht9XO1M+m+j67rfqvtvrd0XfcziQMMv598RfaH5UPRXX6tXDl9/0fDz4Ytjfn6+1vddt9baevQRwe+I/m7rmYj4z2n/74txnCT/OXww+mL7d6dz9kDXdaeCu54C+S7SWvuo4Uros4ftD4s+l/XWY77k6eivSl5urf2+iPgzx3ydV0R/heIfttaeNRzbU6OPg+R9bUb/wYyHIuKby2u8PSK+qPUt0kb529baH2itvXC42nwu+h8Ge621Z7TW/vCQRd6K/hbaXsw7HX0M5AMRsdJa+6bos2bHMvzAeST6GMObhq910V81/to4/ifUXxcRX9P6FkkPRsQ3pn2+N/oi+JWt/1Dei6K/dfqDw/il6Ltm/KWYFcQ/Hf2tVgUy3AJd121GxP+KPrP6pjT05uFr1zM3vKy19oKhOL06fw5XRV8XEd/WWjvdWnvOsK/XpO9/Q/TRtytzwevLdkR/geTbhu+P1trTWmt/ZBj7txHx+a1vM7oWfYZ3Yf3QWvv21trHtNZWhl/qvyr6O1mPRj8P70efl77idPTz99nhZ8jXH/bE3ABrEbE+HNdua+3zoo+3XfG6iPjS1tpHD+f+ar54+OXi1dFnlp8eEdFae1Zr7XODu54C+e5yPvoP4f2P1trF6Avj/xv91dnj+OqI+NbW2vno/9G/7jgv0nXdr0R/K//ZEfHzw+u9JfqrJn97eNqroo8efHA47v9SXuZvR8Tzos/vfkv0V7ev+JDoJ+hz0d/2e0P0k/9S9O/94Yh4LPpIwVdf4xB/PCJ+LPoPOb4nIi7H0W/XVW+M/jbaW9LX3hQRT4/j/xB8dfTH+vMR8b8j4t+X8T8dfX7v4Yj44Yj45q7rfjKNvyH6SMv/TNunr+N4gKN7Q/TzwJvT165rbuj6DhCviv5zCu8c/p/9lejvwr172O8PRZ+dzceU54JrzQ3/KPo7ij8xzOFvjf7nTXRd9wvR//L9Q9FfUX08+s/BLHJf9HPU2eGYnhP9B6Gv/DL/bRHxliGW8EnRz/kvjv7O3n+K+bnvphmuqH9N9D//Ho/+QtGPpvEfi/7zHT8V/bn/mWHoSnvVbxi+/tYhHvJfo7+bx12u9Re+AACY0lr76OgvTK13h+umwV3KFWQAgAVaa3+0tbY2RN++PSL+g+L43qdABgBY7Cuizyi/K/rPuXzV7T0cbgURCwAASFxBBgCARIEMAADJopXCIiKitSZ/AXATdF136IVqzMUAN8eiudgVZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEhWbvcBwJ2n3e4DmNDd7gMAgHueK8gAAJAokAEAIBGx4Eni5scmWjv8PrruuFGJo7wPcQwAOA5XkAEAIFEgAwBAokAGAIBEBpl7yPFyxkfJDt+K1zmK6Szz1PHIJwPAIq4gAwBAokAGAIBExIK7zOFjDIeNPBwlGnEr4hjHbwE3bfy64hcAsIgryAAAkCiQAQAgUSADAEAig8wdbnFW9qA88NT49Nji3xtvRSu3qQxy1+2PtvPxXE92eTqfLJMMwJOLK8gAAJAokAEAIBGx4A5wY9qs1bG8XWMTh32dVo/tSC3hZvus0YhJJSrRpYhD1x2+PdzUPutzp6MaWsIBt97S0vLVx/v7e7fxSHgycgUZAAASBTIAACQKZAAASGSQuQ1uTHu2+czx4pzx3HPzMZSxnHubyjXPH9vhf988Sj44b8/l8NJY1+r3tbKds8z7C8cOoiUccCvkuTg/jojY3d2+1YfDk4wryAAAkCiQAQAgUSADAEAig8wtctxlnw/f23hpqWSQ0z5bya/l59Z9jDLI5bhb3ccRcsdTcia42y/54MgZ5PH+xvnk8n0lr5wzyvtzEeib0TNZHhk4vNXV9dH2VE/65fX7RttbW5du3oHxpOQKMgAAJApkAABI2lR7p9aae6Qc0/EiFXV8Kv6wVOMN5bnLy7MEUW0RlL+3xiby983tv41fJyaOdcrcv7sclej2Fj53b293PJayEvuldVttCTf63rL//L31+6bazlXT7eJMJ1k3tW54YS7mXrW+dmK2cR0tN3Pc7cKFx2/gEXKvWzQXu4IMAACJAhkAABIFMgAAJDLI3CRHWaL58LmzcXZ4nAdeXi7bSymDvLyy8LlLS3UsZ5dLK7m6nPUhM3IH5XhzS7b5LPFsu2aQ9/d309g4O7xfnrs38dyp/ddjzRnlo+WRTSeZDDJPRuulPVue75fLPD31WZGD2nxmjz32viMfJ08eMsgAAHAICmQAAEhELLiBDhejOKhdz1QrtxyVmLsdVyIXKyurs+fWGEUaW6nxi9HY6misRjXGK/LVW4BTEYvFMYoajdjd27n6eG93p4ztLhzLkYqIiN00PtUCru5/qn3cVHREy7dpIhY8WZw4cfrq4xqFOG4Ubi7+NhF3m/Lww+889HO5N4lYAADAISiQAQAgUSADAEAig8x1OHzrtsMuHx0xzh3PZ9IWZ5Bz5rgfn22vrq6NxlZXN9L31bHZ9nLJIK+s1Bzc4uWsp9QMcG7Xtrtb2rOlDPLOzvZobHd3O41dHo3V5+bX2a155ZxlrhnkiUzy1LLUR1mi+smYSZZB5k4z9TmKOr9NjdW5eTxv13l6eeHY6HMkc99X9pGPYWIZ6vr5jyjz0l6a0+pcOD2Hps+KlDaaZ88+Ety5ZJABAOAQFMgAAJAokAEAIFk5+ClwOIddTnpuidCJ5Zunsmw1c5xzxf34LEu8vjYeW03ba3VsdX3iNRfn4KYyelPLNUeMs24125azxTs7W6Ox7e3Z2PLc/sfbu7uHz0hn9diXUoSvm4vGzgbr902dD+D2O+5nR+Y+D1Lzwrm3/MRnPubH1q/5+Fr7HM/F5WdKuhbYxeIe9BHTn8fI82+di8efBxl//oO7kyvIAACQKJABACARseCIDteZav5WXWrdVm5/1ZY84+VFFy8fXeMP9RZcjk6sr983Gsvb82MnZo/XTozGVtfqbb7chmj8z2kqUjDVPmhne3zrbmt7c/Z4a3M0lvc/12bugCW9F6nHurw8vgW5m9ok1T/LfLeydiabjlXUYxPBgFttKkZRo3CjuFttz1baauboxFqdp9N8u1bm27VRFG48VuMYo5ZwU3PhRFu3iHHELccmIiK201yc4211rJ4r7k7+FAEAIFEgAwBAokAGAIBEBpljO0pLoNFY1MxxXTJ0cSu3nLmtGbTarm1j4+TVxydOnJ4YOzUeW5+Nrd83zr2tb5Tc21rKva1MZKtr7m235Hq3Z7m3rcvj3NvWpVm27fLaxfH+Vw6Xga6m2s7NZYXrktHLs+290pJu6u/AFG3f4PabasG5NLV89EpdIrrOzelzHevjOTXPxflx3a6fFamZ5NW0z6XlcQY5Z6lri839siz0Tsod51xxRMTW1qWrjy9fHs/F+XzU88jdyZ8iAAAkCmQAAEhELDjA4W+Tj76rro6Xb70vLb79FTG+lZfjFhHjVZfq6nj1FlyOVdx33+mFYydPPjAeOz27dXfiVIlY3DduUbS6niIfq4dv87a7U1Zo2koRi0vjNm+bF9JqUufH+x+1Yart84rRMdTYRJciFmVlqf3S5m0ptXnbXxqPtf3Z69a/A3kfB8t/78Qv4FaYi0Wl7Tq/TMbdSiu3HKuocbc8N8+P3Z/GxlG4Ohevpfjb8mr9GTM79rmV83bG89J2irjNzcWbF64+ru95KmLxvOd9/Gj7Xe96W3DncwUZAAASBTIAACQKZAAASGSQuWEm27zF4izb/BKms/xYbR+Ul5dePSCDPNXm7fTpM1cfnzwzzradPDP7vvtOl1zzqfE+1zZmObiVtZJBXkoZ5P2SQd4eZ5C3L8+ybpsXxkuYrp2YZd1qznlpaXFGfC5rl5a3rq2O8nKre6Xt0X5X2yItXja82z9cqz9t3eDOs7MzztzmebN+HiTP06tlaem1iVZu9fMgJ0+eufr4VHocEXHy/vQ5kgfqXLz48yH5syER03NP/vxHxDh3vHlh3OZt7Yk0F5fltetnabKu2184xp3LFWQAAEgUyAAAkIhYcCTHXSktpuIXNXKRV9Irqzfl1jp15bz5Nm+z6ERt5ZZjFacfKvGLtH3qzHhlpxMlcjG6rTcXsZi9r9o6badELPJtvY3zl0Zj+XbhVCu3/f3Fq+NFROymVe/29nYWjq0sj2+z7u3V9zV7bttf3BZqqjvbQX93RDDg9suty+ZX0pta1XRxxKLG3XKs4v6HxvP0aC5+sLSAu39x5GJ1o8YfFrd527k8ngtzrOLSufFcvLY+sXJqMhdh29td8EzuZK4gAwBAokAGAIBEgQwAAIkMMsXxlpaOGC+vWTOm4yWRx+1wlueybWmp6bk2b2vp8eLlTCMiNtZT7u30eCy3cqsZ5PufMlve9NSD4xZwp0pLuI2Ts2PYKK2OllMLtr2SD768sz3evjjL/dYWRUvLi9sH5Tzd3PLV2+P89Pb2rH3czs64ldzq6lYaG+9/ueTn8p9PzdotpaWna2uj/Pfj+MtOR1h6Gm6Omg/OLTiXyzy0kubmOhfXz4eM27zdPxrLrdzqXPzA089cfZzn5YiI02Vuzp8PWatz6Mrs2Pd3y9LSpc3bZvoMyPmynHVdwjrb201z8dxnPMbz/Yte9JlXH7/jHa9f+JrcXq4gAwBAokAGAIBEgQwAAIkMMjdFi8XLDE+NRYx7VtblPHPvzbzsdETE+lpdenS2XZclzUtI117HOXd8/1PHubcH7x/n3u7fmL3uxtr4WFdS1nq3ZHUvb48zaufWZ703l6b6a5ZloPMyqdub45zb1ua4n/H61uxYt1bHS6guL8+26znfWRq/zmH/LOsYcGepc+ZUj/qluaWmpzLIZS5OPepzf/qI8RLStddxzh2fSXnkOhYR8cCp2evctzY+npX0PnZLH+RL2+P57YmTs58ry6uLS6S9ncVZ5u3LZS7eulS2Z/PtR33UJ43GfumX3rpwn9xariADAECiQAYAgETEgmObux2XtyfG5paWbstlO93Wm1tqOi1DXVvAldtq6xuztmvrpV3PiVOz22h1+ejcyq1GKp5yqrQhOjG7lXjf+rjNW34f+6Xl2aWt8S241ZXF/xT3U/ugndKS6PLFWbu2zQuL33/E+PzUc5fPaz3n+X3027M/r/pnOVp6eurvQBmztDTcGqM4xNQcHuPWjHWZ+1E7zrl5evHS0/NzcZpDy/LRuZVbjVQ8/cx4WeqHTqZWcuvjfaymFnU7JaZ2aatE0crcmO2lVppbl8bfl5eo3jw/jrCtXRpHTvL5qefuuc994ex1Ns+Pxh555NcXHhs3nivIAACQKJABACBRIAMAQCKDzC3XSqZ1Lgc3Wpa65F+XFufelkt7spXUdq0u37y2Mcuo1UxcXj46t3GLGGeOIyIeTLm3E2t1qenZse+V1kJTObed3fHSzhdPzrJu9Vjz+6jvcWVtcYu8eu7yea3nvJWlwWOUJfY7Ntwr6r/npVGbt/E8sJy25zPI5fMhKXO7Vj4bkee02o4zfz4kt3GLGGeOIyKecnqWVz61Pm4BupIyyLslg3xhbfFcvLVbPvNxevaZjxOnxq3b8vuo73G1ZLLz+Zn7OTb6+Tc+5w899KGj7ccee9/CY+f6+ekGAACJAhkAABIRCyZNtnI70uss/l1sah/zt/wW336au1WVVqRbKSsirazNtlfXxmMbq7PbYXV1vNrKLccq7isRi6nbelW+lbexNd5nPp56rPl91Pe4vLK4LVM9d3m7nvOj/B04buRi6jW1gIPjq/NidtC/9cmxtD3XAq7OL2kuXF4dj+Vo2OpGicKlsbo6Xm3llmMVpzYOH7GoNtMqp3WfaxPHmt9HfY/5/UfUGMXiyOFBP29Pn37o6uPz5x+bfC5H5woyAAAkCmQAAEgUyAAAkMggc1NcT/uvqSWJFz3vms+deJ22tHjp6+U0tlKzdK0+N+Wc53Jmi89BHcuvW/e5PHGso/cx8f7r+HWd10Oa/zswnf0Dbowb1X7xsJ83aPVaW/m+pYnMbZvIMi+tzL5vpYytlvk2z791Lq7PXfR99bl1n/l4jvQ+plplRn2dw39eJ9vYODXavnz5wsLncjiuIAMAQKJABgCARIEMAACJDDI3RdftH/ykhd/bXfPx1POu+dyJ1+n201hZBnovje3uj3Oz+1197mx7qr9mHatLT+fXrfvcmzjW0fuYeP91/LrO6yFdz98B4PjG//aOfx1sep6Y7aOL/To42txPc9p+ncPSc+vY/u7s+3bL2E6ZU/Mce5S5uG7n1637zMdzpPdR5vTRz6aor7N43pz685A5vvFcQQYAgESBDAAAiYgFk6Zu6Ryl/ddRbhuNowCLb1XV21Z7e7vj7d0UW9gZj+1uz7Z3tsdjl3e2Z4/TsqMREZe2tkfb6yvj5Uaz3MqtRio2t8evk1+37jMfTz3W/D7qe8zvP2J8fuq5y9v1nB8lcnHcWMWNinUAY3VezC3HDvq3Pjk2ESnYq/NLii3s7ZRIw9Zsvtu5PJ77ttPYpe2t0dilrfH2hbXFc/HUUtMXti4vfN26z+2JY83vo77H/bmI3eKoxlQ0sDJP3lyuIAMAQKJABgCARIEMAACJDDK33FxOda4F2+J81v7+LE83lzneG2fCdlOWN+fDIiK2L8+yZVuXxjmzyxdn2+fWN0djqyuL/8ls7Y73kZePru3hapb5ic3Zfs5dHu8zH0891vw+6nvcLVnmfH7qucvnda590VSLIq3c4K6ys7O1cGylfKYizwVzn/lI23U+2S1z4c7ubL7bvjye+/KctnlhPPdtnr909fETJzdGY1Of/9gsc19ePrq2h6tZ5scuXpzt88Kl0Vg+nnqs+X3U95jff8T4/Mz9HJv4nE2dmy9ceDy4eVxBBgCARIEMAACJiAXHNtmaa6p129wqcotXq9ufu3WX2prV23ilJc9Wus1VowmbF2atfTbOj2+jra7Pbt0trUz/DrmTjmdja3zLbyW1U6qr49VWbjlW8fi58YpIF87OtjfLseb3Ud/jVr3Nl85PPXf5vNZzXuMh+c9rbmW/w/4d0J4I7jgXLz4x2n7ggbWrj+daue0tjrvtlkjB9vZsfpufi2djl86N57fz961ffby8Ol2u5IjbfWvro7GV1HKzro5XW7nlWMW5R8+Nj+fx2VxcjzW/j7ko3PY4jpHPz3xUcHHc7ezZR4JbxxVkAABIFMgAAJAokAEAIJFB5qboYiJ/OjEWMc5d7e4tblW2szNeInRru7bdmW1vXhhn0tZOzLJ1OXMcEbGUWgJV+2X55osnZ1mzjdW10djy0mwp7r398XvMy0dHjFu55cxxRMSFlHu7cPbiaOzSZNuhsp3OTz13+bzWcz7X9u2Qf5Z1DLi7PPHEB64+fspTnjkaG8/F05nbra00T22O57e1J2bz5lqZi5dXF8/Fezvj7O7l07M5rb7O0srsdfZ3x58H2S7tMfPnPHLmOGKcSb7w+PnR2KUnFr/H/P4jxuennrt8Xj/wgd8Ibh9XkAEAIFEgAwBAImJBUW+Lt2s+65rfmdqBdd34+/ZHqwONfy+banOzV1u5pWhCvTW1tVVWoFubxRFWz5e2P6ll0NLS4t8T98uqS3W1uvXUhmh1bfzPqaXXre3QdrbH73m0mlRp5ZZjFecfG9/Wu5jGNs+X9781jmPk81PPXT6v9ZxPtyGqKz2lFnBzbd6Ou+qeqAbcbo8++vBo++lPf87Vx/MRi3GE6/Ll2Vy0sjKOoq0sz+IQyxNtNfd2xnNNbaV24tRs3lzdKBGLNBfXyNjO5RKxmGg7l2MVc3Pxudn2pUvj9nD5/UeMz089dw8//M7gzuAKMgAAJApkAABIFMgAAJDIIHMkU0sEtzaRV55YZnhu6en9tJz0xBKmNee2vFyWel5JS0YvjdsFLS0tPtacUauZ48sXx/tc25hlkFfmMsizfXSlzdtuySBvX04Z5AvjfeRWbhdLm7eLqSVcXSZ2qtVQPXf5vNZznv88IsZ/XnN/Hw65hPTkMuXAHW/UGnJiaemIiOXllWs+jpifm0f7SG0159qxlbaWo8+DlDZv+WdTnWvqHD/6PEjZR27lljPHEREXLp6dfd/meGw+gzx73XruuHO4ggwAAIkCGQAAEgUyAAAkMsjcMN1UzrjNtmsfyv2l8fZe6j281Gof5Fl2dnkuVzzezlm3yV7HNR+cljDd3hznw+qS1TnrlnsrR0zn3nbLMqk5B1f7e+YcXO11nHPHB+XecgZ5p2SQ83md74Nceh2nfsZzy1DHxN8BOWO4Z+S+yE972oePxpbaeL5tabu1qbl4PNfkZe+3L5e5uMyFaxuz/sp1ieqpPsi1v3Lez9xcnD7XUXsd5/n30qUDMsijnvTj93XmzDOuPj579pHg9nEFGQAAEgUyAAAkIhYc4HhLT9dlhfPS093c8sTj39P281LTJRqRb0e1VtoDlTZzU23n8m22eltvZ/vk1cdbm+NbbOsbZZnUtYllUvP+S7wgty+KiNjdThGLcitx69LsdlxdPjrf8qu38WrkIkcstuYiFmmp6dLWbX+u7VtaTrqcu3HM5rhLS0dYXhruHrVV2dzcOxGryPNEXdY+v26evyIi1i6dGG2vpiWsl5YXx+/qfL9fImQ7o1ai4xhHPoY63+bt+Xjb+HW20/LS9dzl9nncXq4gAwBAokAGAIBEgQwAAIkMMsd2lLZdo2xqK3nckjsb5dd2F+eKd3eXF44ddKx5n7ulrVlehnl9q+Tc1sZt3vLy1nUJ1ak2b/U959zZzvY497yVcnA1yzbKFZeMXt3O72sq91bPR80k52PvYnErt6k2b1q+wb2jtoacmotj7vMpizPIOymrW+e+lZXyeZCVNBfX5aunPg9SW8vtHm4J7e3yOY7R2ETmuL7ufFvN8Tng9nEFGQAAEgUyAAAkIhYcUb49dfhIQzdafW383OWl8evkldvqqkv19v+U0TGU48mtfmpbnbyq3NZqva23Otoer9a3OPJRz0dtNTQV+cjHs1Nu1eXbfHV1vNrKLd/Wy69Z9zl/bOPtfC7nVtJLf85Hi1GIXMDdaveA1mRT8aq8Muf83Deb72qkoUbaplZObelaYBdlJde6sl5uMzoR+ahz8Xh+3V44FjGOVdRzJ2Jx53AFGQAAEgUyAAAkCmQAAEhkkLlhcrastvkZ585KBqwbP3cpDe/tHW5p62sZLYnclVxtGqu5t9XVWbZseblmkMf/ZJaWFmeQD3ts/XbOIE+0gJvIttVccX3uZCu3idzd3NLgo5xxzSAvzhJr7Qb3poNys7kd5PznMRZ//mJ3tOzzARnkPP+22h40ZZDLnDXV9m2qHef8HLqTHtfPcRx+vq0/G7h9XEEGAIBEgQwAAIkCGQAAEhlkrkPNlC7u+ztlLnOVf207Qk/IyaWN53pdzva5sjzuZ7mzs5rGxn2Pl5ZrBnl2sLVn81Qf5JqDy70498t7zn0y67KkuznLdsAS0dO9jmfPrfvfnzvWnO0+7nLS8shwr6jL2t8Kz3jGc0fb+TMgk0tdF1M96muP5JyXrjnjs2cfOfQ+uTu4ggwAAIkCGQAAkjZ1G7S15j4ox1Tb7Cy+5VXHpm6VLaUYQytt1ZaXy3ZuwTa3LOnsublVWz+2eMnSpRqjmDjW6YhF2c639SbiF/MtgRbf8qtRib2J507tf+oW5NHauplOsq7rDn0f2FwM0x566EMXjk3NxVWeb+faWKaxc+cfPeohcodaNBe7ggwAAIkCGQAAEgUyAAAkMsjcJIfPHE+NT+WTax64Li86zhKP88njLPP4dfL3zWegy3LSE8c6Ze7fXdreL8ti5+fOLwOd2sNNtGOb+96aK07fW7/voPz0YcdkkMdkkOHmOH36oUM/d6oFZx27dOnc9R0YdyQZZAAAOAQFMgAAJCIW3CLHi1xMtU6rK9fVlmwt7bO2hBuvgDfRZq62q6v7qDGPYxrd1iurN3UpmjAVf6irPnX1uTHx3InbiiIVN4eIBdwaGxunFo7VVm7Z7VghkFtPxAIAAA5BgQwAAIkCGQAAkpWDnwI3Qo5QjuM+Ncd6+GVBa1uz8egor7xfl3Zu+YmjsaWl/TRUMsj7U3npw/++OZV7m8oA1wxybtfWxXR2OG/PLaE6GpuOu8odA3eTy5cvXH28urp+G4+Eu4kryAAAkCiQAQAg0eaNO8DhV6C7US3hDvs6tc1bjWNMyfucilTMqdGIOFz8YWpFqDp+/NZtc88+wnPJtHkDuP20eQMAgENQIAMAQKJABgCARJs37gBT8crFLeFqjng6VztujzadQV78e+PU983bO/gp1zD9Pg7fHu5or6t1GwBc4QoyAAAkCmQAAEgUyAAAkMggc4er+dfDLkN9wKtOZpkXZ4ePlkE+nqO8r5vzXJljAJ7cXEEGAIBEgQwAAImIBXeZ47WEu2F7n16a/Ya8zlEc/3XEKABgEVeQAQAgUSADAECiQAYAgEQGmXvI4fPJo++6CVnlg9yafcoZA8BxuIIMAACJAhkAABIRC54kjhI3ON5qeWITAHBvcAUZAAASBTIAACQKZAAASGSQYY6cLwA8mbmCDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJC0rutu9zEAAMAdwxVkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAMn/B2PUnlB/GQpdAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAF2CAYAAABzg27uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsdElEQVR4nO3deZQtW10f8N/u6Q7wHiiDCCgKiCgsHFBAGTUOEDUxRpS4QEEFY6JkGWeNQzBqWDFqUOK0Eo0SIqjoksE4LOQ9HoggQhjEgUFklEHw8d4detr5o6rv+dU+XdWn+3bfe/vez2ett96p3nVO1anuW+fXVd/+7VJrDQAAoLN0uXcAAACuJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgXwNKaX8fCnlBy73fhyWK+X9XMx+lFKeWEq56bD3CbiylVIeXkr5q4t4fi2l3Psw92mf2//bUsoXXK7tw1FTIF9F+hPW2VLKLaWUD5VSXlhK+bid8Vrrv661/shFvP4XllL+uJTykVLKB0spry2lfHcp5eThvIP9udj3M6aU8lellK9Kyw/tP4zar91SSlk5qv0Ajo9SyveWUl7UfO1vRr72uFrrS2utn3xp9/LSKKXcvZTyW6WUD5RS/rGU8vpSyhP7sU/oz6crh7zNy/oLA1cfBfLV58tqrbeNiI+NiL+PiJ85jBctpTw2In4zIp4dEfeotd4hIr46Iu4eER839dxj6MaIeGRafkRE/OUuX3t5rXXzUu4YcMW6MSIeWkpZjogopdwlIlYj4jObr927X/dq9msR8Y6IuEdE3CEivja6z6OLdtiFNYxRIF+laq3noitoP3Xna6WUXyml/Ke0/F2llPeUUt5dSvnGsd/ASyklIn4yIp5Wa/2lWus/9Nv4q1rrt9Za/6Zf70GllD8ppXy4f92fLaWs9WNzVw1KKS8ppXxj//jepZQb+qsNHyilPGdn26WUnyqlvK8fe10p5f7t+ymlfFQp5QWllPf3V89fUEq5e7OtHymlvKy/Av4HpZQ7jhy+G6MrgHc8PCKevsvXbtxlPx5VSnlnKeXb+31+TynlSWk/7lBK+d1Sys2llFdGxL2aY/25pZRX9e/1VaWUz+2//nmllNen9f6of/7O8k2llC8feT/A0XtVdAXxp/fLj4iIP46Iv2q+9pZa67t3zhU7T+7vAH5Hf477x1LKc/LduVLKd6bz9dfnDZdSbldK+dX+/Pf2Usp/KKUs9WNvL6U8sH/8+P48/Kn98jeWUn6nf7xUSvmeUspbSneH8LmllI9O23hC/1ofLKV8/x7H4rMj4ldqrbfWWjdrra+ptf5eP7bzy8GHS3cX7nNKKfcqpby4f+0PlFL+dynl9s2x+e5Syusi4tayR5FcSvnhUspvlFKe1Z/vX19KuU/prvK/r5TyjlLKF6X1n1RKeVO/7ltLKd/UvN7oZ2Up5UQp5SdKKX9XSvn70kXuTu1xfDgGFMhXqVLK6eiu8L5iZPzREfHvI+ILorui8cjd1ut9cnRXin9rj81uRcS3RcQdI+JzIuKfRMS/WXCXfyQi/iAiPqrf1s6V7y+K7kPlPhFx++je0wd3ef5SRPxydFcsPj4izkbEzzbrfE1EPCki7hwRaxHxHSP7ckNE3K+U8tH9h8xnRcRzIuL26WufG+NXge4SEbeLiLtFxDdExDNLKR/Vjz0zIs5Fd4X/6/v/IiKi/zB6YUQ8I7qrLj8ZES8spdwhIv4kIu5dSrlj/+Fw/4i4eynluv5k/MCIeOnI/gBHrNa6HhF/GrNfpB8R3b/Jm5qvTV09/qqIeHREfGJEPCAinhhx4Xz9HRHxhRHxSdGdt7Ofie6cc8/ozuVfG925LqI7nz0qbf+tMTvfP6Ifj4h4akR8eT9214j4UHTnq+gL6p+LiCf0Y3eI7jw95hXRnfceV0r5+GZs51jcvtZ621rrn0REiYgf71/7U6K7K/nDzfP+VUR8Sf+8Re7cfVl0V7I/KiJeExG/H93nxN0i4mkR8Qtp3fdFxJdGxPXRHbefKqV8Zv/e9/qsfHp0n0+f3o/fLSJ+cIH94wqnQL76/E4p5cMRcXN0J9P/MrLeV0XEL9da31hrPRMR/3HiNXeutL535wullF8v3ZXiM6WUJ0RE1FpfXWt9RX/F4G+jOwFNFd7ZRnTF7V1rredqrTelr18XEfeNiFJrfVOt9T3tk2utH6y1/lat9Uyt9SMR8aO7bPuXa61/XWs9GxHPjdlVnfa1/i4i/i66q8SfFhF/0z/nZelrJ6P7MBx7L0+rtW7UWl8UEbdExCeX7jbrv4yIH+yvrLwhIv5Xet6X9Nv6tf4Y/p/ooh1f1t8R+LPoPlw+KyJeF90H70Mj4iH983b7xQG4dG6IWQH48OgK5Jc2X7thl+fteEat9d39Xbrnx+wctXO+fkOt9dZIxWN/XvnqiPjeWutH+nPvf42umN3Zp51z4cOjK0R3lh+Z9uebIuL7a63vrLWe77fxlf0v5F8ZES+otd7Yj/1ARGxPvI/H9u/7ByLibaX7e5XPHlu51vrmWusf1lrP11rfH93Fgfb8/Yxa6zv6c/EiXlpr/f2+mP6NiLhTRPznWutGRPx6RHzCzlXqWusLa61vqZ0bortY8/D+dUY/K0spJSKeHBHfVmv9h/6z58ci4nEL7iNXMAXy1efLa623j4gTEfEtEXFD6XJvrbtGlxHb8Y5d1tmxU3h97M4Xaq2P67fz5xGxk6+7T+miDe8tpdwc3YliLMbQ+q7oriK8spTyxp1biLXWF0d3JfiZEfH3pZRfLKVc3z65lHK6lPIL/S3Am6O7SnP7/sNjx3vT4zMRcduJ/dmJWexcBYqYXQl6RET8af9BsZsPNlc4drZ1p4hYieGxfnt6fNdmeWf8bv3jnStBO1d9XhLdh0j+kAMunxsj4mH9HaM79fGzl0fE5/Zfu39MX0EeO0e15+t8nrhjdHfE3t6M5/PGw/vPgeXo7oY9tJTyCdFddX5tv949IuK3+wsfH46IN0V3V/Bj2u33RfroL+S11g/VWr+n1nq//vmvje7iTdlt/VLKnfuLLu/qz9/PivnPjqnPqN3kzPPZiPhArXUrLUf0x7eU8phSyitKKf/Qv/d/mrY/9Vl5p4g4HRGvTsft//Zf55hTIF+laq1btdbnRXeCe9guq7wnhrfIpv7Q7i8j4l0R8RV7bPbn+nU/qdZ6fUR8X3RFb0TErf3/T6f1LxTutdb31lqfXGu9a3RXMv77Tsar1vqMWusDI+J+0d3K+s5dtv3t0UVBHtxve+eKza4n5AXsFMg7V4EiZleCLuSP9+n9EbEZw2Odbz++O7oPqWjG39U/bgvknStDCmS4MvxJdEXnU6K74xS11puj+7f9lIh4d631bQd43ffE+HnjAzG7A5fH39Vv/83RFdtPjYgb+6uc7+3356Za686V4HdExGNqrbdP/52stb6r3X4f4bvDIjtea/1ARPxEdIXmR0dE3WW1H++//oD+/P34mD937/a8i1ZKORFdfPAnIuJj+gs/L0rbn/qs/EB0xfb90jG7Xf+H8hxzCuSrVOn88+jyV2/aZZXnRsSTSimf0p/sRjNTtdYaXQH6Q6WUJ5fuD+JKKeWTors6sOO66KIdt5RS7hsR35xe4/3RnbAfX0pZ7q8QX/gDtVLKY8vsj+o+FN3JcKuU8tmllAeXUlajK7LPRVf0t66L7kT14T7L+0MTh2cRN0bEZ0RXfL6s/9rro8sGfl4coEDur148LyJ+uL/i/akR8XVplRdFxH1KKV9TSlkppXx1dH9k+YJ+/OXR/RLwoIh4Za31jdF9KD74IPsDHK7+9v+fRZdZzX8TcFP/tYP+O31uRDyxlPKp/fn6wvmtP688NyJ+tP+bhHv023pWev4N0d9R7Jdf0ixHRPx8/xr3iIgopdyp/wyJ6P7g+0tLKQ8r3R9ePy0m6odSytNLKffvz2PXRfdZ8OY+Bvb+6OIZ90xPuS66KNqHSyl3i90vghyVtejuuL4/IjZLKY+J7m9fdox+Vva/XPxSdJnlO0dElFLuVkr54ku29xwZBfLV5/mllFuiK1R/NCK+ri+kBvq/KH5GdH9l/ebornxEROwaG6i1Pie6LNbjo7vS8IHoThy/GF2+K6L7I5KviYiPRHfSeE7zMk+O7sT3weiuBr88jX12RPxpv++/GxH/rr/Scn3/Wh+K7rbhB6P7Tb/10xFxqt+vV0R3m+vAaq1/Hd0fbryn1vrh/mvbEfHKfp9ePv7sSd8S3W2990bEr0T3h4U72/xgdH8o8u3Rvc/viogv7a/A7NzW/POIeGP/B0ER3fft7bXW9x1wf4DDdUN0fwicJwB6af+1AxXI/fn6pyPixdGdr1/crPKt0V1AeGu/3WdHxP9s9um6tP12OSLiv0V37v2DUspHojuPPrjf/hsj4t/2r/ue6M7H74xxpyPityPiw/0+3SMi/ln/Wmei+2x6WR9LeEh0ud7PjIh/jO4PlZ83eUAOUX9F/anRfZ59KLrPsN9N43t9Vn53//VX9PGQP4ruQgbHXOkuDnKtK6V8SkS8ISJOVL19AWCOz8prhyvI17BSyr8opaz1fzzy9Ih4vn/wADDjs/LapEC+tn1TdLmrt0SX6/3m6dUB4Jrjs/IaJGIBAACJK8gAAJAokAEAIFmZGiylyF8AHIFa68KT2DgXAxyNsXOxK8gAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQLJyuXcALr9yuXfgItTLvQMAcNVxBRkAABIFMgAAJApkAABIZJC5ihznLPFBHfQ9yy4DwBhXkAEAIFEgAwBAImLBMXP4MYpSjm80o9aDRiWm3rP4BQDXNleQAQAgUSADAECiQAYAgEQGmSvQwTLBh5UlvpIzyW3m+KD7Op1d3us1ZZQBuLq5ggwAAIkCGQAAEhELLoODRxj2EylYdN2LiVSUcji/Y9a6veD2hvu6nzZved2p97z3a2oRB8DVzRVkAABIFMgAAJAokAEAIJFB5hK59NnhqfFFs8OXquVbKcujY1OZ4Lx7i+aY99rGXu958RZx8sgAHE+uIAMAQKJABgCARIEMAACJDDJH5HAyx/sZm8oV7+t19tOn+bAyylO53pwzbnK9w+zw8uhYN57HFs8rT01vvb8pq2WSATgeXEEGAIBEgQwAAImIBRfhyo5R5OW52MSCLeAuppXcfkxFFQZjTTRisPnmNWoZj2O0vxsfVru48e1FmKIagOPCFWQAAEgUyAAAkCiQAQAgkUFmny5+Guj95HrbzPFwbI/2bAd9nQWfN7ff+8hkt+3aBmNtljhlgmtdGl13LvPb5pXTNveTT44Yvs7UtNSLjs3TEg6AK4cryAAAkCiQAQAgEbFgD4ffym2v1m0LxyiasaWl4Uxy+blTY1PRjPZ5rUXf83S8YHp8e3srrzh8XoxHLAbPa8bbsen4xXhUYqol3MFn4Ov2KL3SHusCwOFyBRkAABIFMgAAJApkAABIZJA5sP1NEb1467TpdWfLS0vDsTYvvJSfO5FXbl9nkF0uTQZ5Yl8vxnA66SZLXLd2Xy8itre30+NhrnipOXbbKS/c7nd+nbY9XNvmbWi6JVx28BZwAHBpuYIMAACJAhkAABIRCxqLz3I3NTbVyq0da6MAU23Wchxifmx8eS6OkaITZSKqMRW/iGjbzu3j980mxjDdrm08RlHzWB2ODWITEVHSc+e/P7Ox7e02fjHc9dwSbnsuUpGPwcFawM2Pm2UPrkV3ucs9Lzx+73vfehn3hGuRK8gAAJAokAEAIFEgAwBAIoNMTOWOj6KVW5s5LgfMGbdjy8vtuitpbKUZW9p1vXasfd5e02SPjU1NydyOt/ncra3NC4/bXPH29vhYfl677tbWePu8nEfulsdbwi01b2uYST5YC7h2m/Mt4ExDDdeCtbWTFx5/4ic+YDD2tre97lLvDtcYV5ABACBRIAMAQKJABgCApExN8VpKEfC7Ki0+PfJ0zvhgvY6nMscRw9xvm/FdWVmdrddmh5sMcl63zRkPx9os88ro2Hx+Or/nxX/fbDPJw37Gbc54lgmezxXPxjY3N5qx4bp5fHtrmDPeSuu2r9Pu61Qmuqb9ad9HPtfMvf/mPDRcdz+noeNzyqq1LvwP0bmYa8G97vUZg+X9nIv/4k0vP7od46o2di52BRkAABIFMgAAJCIW16SDtW5rtbe8Fm3lNhWpaF+3HVtenkUjckyiHWvHV1bWRsfmt5Haw+0R48j7OnWsWvPTSadowkT8oY1Y5OU2GrG5ud4sz8a3ttp1x8emYh3zY6ldXTMt9vZEK7t2Cu1sKn6xy9oTY1cWEQuIuO99H3Lh8WGdi1/1qhcd5i5ylROxAACABSiQAQAgUSADAEAig3zNOPzppKcyyG3Ls0GWbCJzHNHmgxfPFa+2OePVtfF1B2Mnhs/L+zqXc26nnh5v8zY1XfJ8m7PZ8lzOOOWDN7fa1m3nLzze2Ggzx81yGt+Yyyevp8dtBnk8r9xmh/O+t+9jajrt9nUWbfO2dwu4K/cUJoPMteh+93vYYPlSnItf8pJnH2xnuSbIIAMAwAIUyAAAkIhYXDMWi1Hs1eYt37pq27UNxxafna5t15bjEO3Y6uqJXdeLiFhbOzm6bn7cLc/WXVtrbuvl7TcRj/Y2Xz4GS2X4viIfuzZSUJsWaKk92lbbri1FHNrYxPp6jlicG4xtbJwfXV5fH66bX7d93lT7uPnZ+xab9a+NmMzNyDcRx9jfLHtX7ilMxIJrxQMe8KgLj6+0c/HzX/DM8R3nmiBiAQAAC1AgAwBAokAGAIBEBvmqdTht3Q46nfTycptBXnyK6Nz2ZzKv1oytnRhmkNfWTu36uFs+OTo2zC4Pc87Lq+37SlNoL4/nt+famm0Nl7fS9NJbG8N8cm7f1uaD19fPpsfnRsfm1j3frLuxnyzzbH/2N2V1yiBPTEMdMTxeUy3gWsdpGmoZZK40d7rTxw+WnYsP51x8002/GVy5ZJABAGABCmQAAEgUyAAAkKzsvQpXu6k+yFPrxlxeebwPcl5eXhr+2E33QW56HadM2okTpwdjJ06capZPT6w7W1470eTnTqXc24lm39aafV/N/Z33kXvbHi5vbsx6Bm+uD/sHb5yfZXnXz45Prz3XJ7T9HuSpWCd+N56fBnu8D/Fcf+el3L94/HW2mrHSZObzulM/k3v3QQYW5Vx8NOdijidXkAEAIFEgAwBAImJx1dhHNGJyvXb66OZWVdrOVIxifqrp3IJnvAVcxPBWVTt9dG4fNHcb7+RtBssn0/LJduz07LknbzPcxonT47f15m7zpdt6yyvN75sT05tubQ4jBvm2Xr6N1y6fPzO8rXfu1tn+tMe1LE3sT6NGika07dDmbklup8fDNkjLdbYPUxGLvVotDX/uhu+jNtN0jz+v3U77/sUzIHMuPppzMceTK8gAAJAokAEAIFEgAwBAIoN8DZrLFe9j6umYaAmX24jN55NnP2rzrYTaLNn49KZ5KtK2XVCbbTt9+voLj0/ddphtO3nbU2msyc+l3NuJU830qm0Obi3lzlbG89tzmdsm97axPsu2tbm382dnU5q228/Trba5u/Z7MNDuT8254uG+tTnjrbTc5oHzc5eWmgxyzi5P/FxFDLPu7bHT9g2OhnPx0ZyLOZ58FwEAIFEgAwBAImLBwFRbt258PEaRb5MvNS3GlpfzLEftTHpro8sXc1sv38o7fX1zy+/62XNPXTe8rZdbDZ083bQ2OtXMnpRmc5pvs5Zu6223bc2G0YSNNGPT+tn1wdi5M+dm22tv663k1nrTv+/mqERtYhQ5NrG1NZw9an55fN3l5dnY9nbzvDLeIrCNdcTETHpSFHA0nIuP5lz8NU/43sHys3/tx4MrnyvIAACQKJABACBRIAMAQCKDzNz00s1gszjR5i1PNV3aNm+zbezV5m110FqoyZ2t5elNm9zb6Sa/lloG5Zxbu3ybJhN36vo09enp8alPI4a5t+XV4T+nqdZCWxvDfG7OvZ0/c34wtnpydnxWVsa30WpzvXlK1TZ3t7k1a2e0uTlsbbS11S5vpsero+vO/wyMT0M9lzMeTA0bzbqzn6WpaaeB/Xnzm189WP70T//8C4+diw/vXMzx4AoyAAAkCmQAAEhELI61fcyAdwhj7XjbAi7HKMpSO6tbbvM2vPWeW8BFRKyszG6dra2NtxZaOzEcyy2BIoazMrXtg/KtvNvcvmk7dF1qO9TM+nSquc13Ms3etNq8j6V0rLab23obTXu0c2n2prMn2tms8rEb/5127tbhZjMD3sZseXN9GJvY2Jwdn83NYWujdnl5ebY8NWNiWdpoxmb7vr118Nkc62QLuPFZ9+Zn2cvP1TsOWs7FR3Mu5nhwBRkAABIFMgAAJApkAABIZJCvEXtliy+sNzG19NzrTLSAm8+mzl6nzRzPZZDTcjv1aZ7udO3UMPfWtv3Jy3OZuNQ+KOfcIoY5uOtv04w1WbtTa7P9W2unN53Iv643bdbOrs9yvbeuDnNvSyvNlN5Jbh+0udm0Kzq/Mbq8fm6YK15P25+bbnbi+9V+7/L3uf0ZWPRnp1tOP3f7mFt6r0wyMO4BD3jUYNm5+GjOxU956o9dePyLz/i+0dfk8nIFGQAAEgUyAAAkCmQAAEhkkNlXD9qcDZ3LGaephdsekTmrOtcHeanJuKapp1eWJ6ahbnpUtssnUi6unaY0L7f9NXPW7Xanhj07b3uyyc+l3NuJNq+WDt12E4U9vzHMpJ1I05auTPTX3G76aW6maVI3zg1fc/3sMGecp02dO3bpuLbHfLmZCnx5wZ7W7c/A8Odj+Lzt7eH7KmWW55vrXjzIEwZwQPe970MGy87Fl/5c/LXf8IODsV/9H08b3SaXlivIAACQKJABACARsWBfJtvF5Wmop1rANa3jlpbHb9PP3d5fna27sjb88Z26zbd2atiiKLcdaqcsze2D2tt41ze3+fK6q/toLbTR3AJsn5ttpvZB5083t+5S+6DzZ88PxqaOR3vs8nGdO+ZtK7e0r+33skz8DLSt3caeBxyde97z0y48di6+8s7Fj/3q77zw+Jab/3Ew9nu/94uj+8bhcwUZAAASBTIAACQKZAAASGSQr0Ht9NHNYLO4nxZwKX86MWV1adt/zU1JnFvJte3iUu5ttZmiulleXVtNj9ux2fLJtWE+LLcLyo8jpqc3ze2BIqZzb+ebqUizjWbq07xuu69n0vuYf4/DdfPxaY9VPq7tMW9/XvL3q/1eDr7Pcz8DE/nkxtS01LnvW7tvtQ6PHTDkXHx8z8WPfOTjBss33PDro/vOxXMFGQAAEgUyAAAkIhbsS3vbfHzF8d+92tvik7P1lWa2vuXULm5p+LzlleZ2YVqea1+UbmutNm3M1tK67YxMbQugfCtvtbmtt5xuj22l9kC72Uy38tpt5v1p9zW/j/Y9LjXHIx+f9tgNjmsZv83aLZfRseGKi//+vfDPFbAvd7/7fUfHnIvnHadz8YMe9CUXHr/ylS8MDpcryAAAkCiQAQAgUSADAEAig8yB7Wd64Mms6tTrTrWS22Mq48G6S+NjS5PTYk/sW7O83LZMyus2Y5sH3ObUvk69x/4Lo2OLHvO9LPp9NrU0HJ073/keo2POxVfnufj+93/EYPkNb7hxdF0W4woyAAAkCmQAAEgUyAAAkMggc2DtlJ3T6073nhx93WYbeWxu+1Prbo+PbU88b3tuE+PrzvXXnOi9OfU6U9uc2tep99h/YaHtTx3HvSz6fd7PawL78773vf3C47YPsnPx1Xkuljk+fK4gAwBAokAGAIBExIJ9qbHgrfGJ23jtLb75W0yz8e26NRjb3kq3uJrbWFubw9fdTsvbW8PX2drYvPB4Y2tzMLae1j2/sTEY22imHj2/OXxultsHte+xfd7GxDbz/rT7mt9H+x63m+ORj0977AbHtTnmU9+vydu1+7mVu+jPFbAv73znXw6W73Wvz7jw2Ln46jkXc/hcQQYAgESBDAAAiQIZAAASGeRrUJtjKmU5DzbrjmecptrTtJnSvM3atNnZ3h7PWW03626lbNfmxjAD1i5vrG+kx+3YbPnc+jBndnZ9/cLjEyvDfyKry8sxZrPJnZWJ3NtGs+6t58/vuv12ud3X/D7m3+Nw3Xx82mOVj2t7zNufl/z9ar+Xg+/z3M/ARFuoxlSro6l9A6a95S2vGR174AO/eLDsXHxln4tf9aoXBUfHFWQAAEgUyAAAkIhYsC+Tt8YXnB1ou71l37b9Se1ztjaHt6a2NtJtvfY21vmN0eX1s8NbZefPpNtoJ4btgm5N7YNWlqZ/h8y35040bYeWZnf15mZkatsH5Vt3t5w7N9yffMvvzHAsv4/2PU4dj/bY5eM6d8ybdkb5+9V+L/czu9bY84DL49Wv/v3B8kMf+hUXHjsXX/5z8Y03Pje4dFxBBgCARIEMAACJAhkAABIZZIb5zzIxFm3bn2Y6zTrenibnWNvnbW03GayUu9rcarJbG7Ns11SuKyLi/NlZJuxckxdbPZmybavDfwZLKxPtg5r3lacpXWvaDk21Flpvctc595ZzbhERN996ZrbeLcP3kd9X+x7z+48YHp+5Y5eOa3vM57KH2xPfy8H3ucmaD34+Dj6FqrwyXBove9nzLjz+/M9//GDMufjoz8V/8Ie/HFw+riADAECiQAYAgETE4hqRbyvl201z60Xbmqu99b2UF0a3MRe/SLfD2rZh7fJmWt7cbNrlbORWOmuDsfNnhsurqWXQatM+aCXNyrQ00T5oe3P4Ps6fHt4OO7mWtrHc3B5Mx3l7bvam4XvOszK17YPyrbwzHzkzHLv57Ow1bh1vO9Qur7e3/NJxbY/55sT3a64F3PZ4BGdqdrypKM/cz+QE8Qs4Gi9+8bMGy49+9JMvPHYuPrxz8fOf/7PBlcEVZAAASBTIAACQKJABACCRQT7WcpZqmCtus5hTbW5yJHlqrB1vs6E5f1rbFl8LtgaLiNjcTPms9WGWbWVllvNaXRmOnbt1mG1bXp21+lleGW/708rvY3NjuG/rTUueM2uzf0LLTYuiqWO+1bzuRpputM2r5ZZBOecWEXHrzbfOxj4yHDt7y3A55+LWm/ZF6+tn0+PhWP5+REy37Mvf5/mfgfFc8UHbuu2VOZ4el1eGg8rnDOfiwzsXc+VwBRkAABIFMgAAJApkAABIZJAZ9JwtZbkdbBbrro8jhnnUPK1wN5ayZM3UxcvL41NtLi8P+0mur6dpSZeHObelZnrR5ZXZ739T/TXb97GZpizdODc+ZWpExGrKvbXbL0sp97bdHKtmetOce1s/O+xDPJjCtOmvmbNuZ24e9uU81+bezsyWz58frru+PnvdjY3hNvL3I2L4/Wu/l4M+yHM/A7Pl+b7Hbf/tqdzx9ugYcGnkvshf9IVPGow5Fx/8XPywh33lhcc33fSbweXjCjIAACQKZAAASEQsGJi79V3Gp55uW3wtlaU0NrwNnluDLTe3v9qpjfP4xsZKM5Za+Sy1t9GGv+8tNePZYOrrZgrTjdQ+qL3F1k6TupqmN11aGW5/qrXQ9ubw+Gyk6U03zo/fSmzbDuX2Qe1tvDwtakTEuXOzNkTzEYvZc/O00xG7TD2dlnNbt4i2BVwThZiYiryNTeQ2cKaPhitbez5xLj6cczGXlyvIAACQKJABACBRIAMAQCKDfA2aynROTS0dEVEmsqHbuV3cxBTEm5ttPqxpCZTyam12LbcIypnnXXd+sP1meuuUO9vaGM+9tTmzNve2kqY0XV6Z2J92etMm95anUW1zb1P7k1sN5dZBEcOcW7u8rwzyvtq8TUxDnX4+9mrrtuj00vLJcPm15xPn4sM5F3N5uYIMAACJAhkAABIRi6tGe6t5eIsr34ouE7e/5mcpG/4Oldu+tbfK8uu229jaytGI5hbbUnMbK0cwmtcp6X2VvX6/y7fim5ZjefakzfXh9tfPzSIFc7fx1ob/ZPJtvaWlZl+nWgs1sznl23qb68PWacNWR8PbeuvnU9uh5lbd1PL588PbeOvnZ7cH86x6EREbTZu3ra3Z/rSzUG2l5bmIxfb4WI3xiMV+Zs6bjlyIY8BRaM8nc5yLR5enzsVcXq4gAwBAokAGAIBEgQwAAIkMMgdu+1aizXKlNm+lzSfPfhfbaqYnjs3xvFibZS5ta7e8b9Fk21J2davJvG6mHO3G5qnB2Pp6yr2trg3GlleblnRpWuyl5X3k3raaVkMphzfX6ii1WWtbsOWWQG12uG0XNFi3ybmtp9edmlq6W54du/Z7OZUzHkw9fcC2bnutC1x6bcbWufhwzsVcXq4gAwBAokAGAIBEgQwAAIkM8lVrvC/y3PTRB+yLnHsid1+Yrbu93WbAUq/LJoLcmtqfqX1re+nmzOvW1nCjw+mShxnblZVZ1m1ledh7c3lluDyc+nqYiZua3nS7judzt5rpm3NGr93X9fVZXm1jo+lf3GTZ8nKbkcuvOzfVdLM/W1tTU03nDHLz/clTTU/0PW7XnbJ3HlleGY7aX7zp5Zd8m495zFMGy5f7XHzjjc/de6c5VlxBBgCARIEMAABJmW7x1d5D5/gajy1MRRrm26zNlpeWlkfH2nZsy8sruz7e7XVWUoxhuYk45LEchYiIWG2WV1JboLl1B2Mnhs/L+9pEKtp9z++zfc9TrYXm4iF1PA6Sb/NtzkVFcsRij/ZsuUXRRCu3qUhFO962csv73r6PQcRibrrXZurpwVTTi7eA22WNPcYvn1rrYlmicC6GvTzykY+78PhSnItf8afPP/jOckUZOxe7ggwAAIkCGQAAEgUyAAAkMsjXpINljltTGeSlNgOW1s3teCKmM8ntWM4kr8zlgxfPKw9zzu02UgZ5qdm35eF7zvu6aHu6iOkM7vZWk+vdHs/15uU2Ozw5RfRErrgda7eZ93V+LOWM26mm95FBzvY31fTxOWXJIMPReNCDvmThdfdzLn7Na//o4naMK5IMMgAALECBDAAAiYjFNWnxKMBU5GJ+bGl0LEcuylIbUxiPXLTtegbRiD3iD3ndpWbd4djweXn77dh8dGS8zduUudZCKZqw3YxNRxryDIXtrHbjswdOxTja12n3Ne/D3Gx5+fbkXPuk3LptfGx+3f2cho7PKUvEAi6N+9//EaNjU+fiyzFDIJeeiAUAACxAgQwAAIkCGQAAEhlk4iimoZ6a6nOqBVzEMJM8lwEetItrs8PtulPTWy/tul47Nj+19HjuupXH2pxba2oa5qnMb84Zt2PzeeXcLq5pwZazw217trnlg7Vym3qPB59O+vieomSQ4dK7170+Y3I8n6fe+tb/d9S7wxVABhkAABagQAYAgETEgsb0Xd9FIxcHbQHXrzAbO2D8ol1uW8ktlTQDXjs28by595WP1z7avEXbWijG25rlSEMbdxi2hxuPQrTPnYpRTD2v38G0zYO1cturddvVGqvIRCwALj8RCwAAWIACGQAAEgUyAAAkK3uvwrWljToOozk5G9rmcadzozmPutSMDLOqOdfbxGEnc6xzy+nJW23uecGccc4q94Oj616Mwb63GeSULd5PPnn+dcbXHbRum8sOL54l1soNgKuBK8gAAJAokAEAIFEgAwBAIoPMgbW50ZzHXTyPHDH9e9pcCPnCw62JbGxExPZEP+Wcld3eanLEE89rLdoX+mL6/g7ywm2ud7J/ctMzOfcvnsgr19gj230IvY73Oh4AcDm5ggwAAIkCGQAAEhEL9jDd9m2w5gFbwNVmiuThlNXjLdhKsy/bW802Sm4X17SSm5wWOz9vfN9a7f5MaWMMg7EDtk6bel7/hdHt7+d1Dhqd2F+sQgQDgMvHFWQAAEgUyAAAkCiQAQAgkUFmn3I2dLE8csR0y7PpvPJ4S7i9XmeQCW5ztZPt2ZbS43b77bpHMNX0xNhcrni44nBxHznj4XoyxwBc21xBBgCARIEMAACJiAUXYeq2+OJt3va1xdQSbj7esNSsO952Lu/6/FjdbbXdHVLEoo1HjK62xyx3i4+JUQDAGFeQAQAgUSADAECiQAYAgEQGmSNysCmqI6Zbwk1ucWLK6vZlciu3qW3s2cbtEsRqFz0Gky3g9nid/R3nRdeVOQbgeHIFGQAAEgUyAAAkCmQAAEhkkLlEFpuiOmLxTPDF9Fae7qec15t+nZxlvhh75YfHn3f42eGL61ktdwzA8ecKMgAAJApkAABIRCy4DPa6Db94S7jBs8rBpreeajO393O39l7pEB3elN1iFAAwxhVkAABIFMgAAJAokAEAIJFB5go0lXE9WD55ykGzy1eCo9nX4/P+AeAouIIMAACJAhkAABIRC46Zg8UvJl/xGEUqDs+1+J4BYDGuIAMAQKJABgCARIEMAACJDDJXkYPmag+WXb4yyBIDwGFzBRkAABIFMgAAJCIWIKYAACSuIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAECiQAYAgESBDAAAiQIZAAASBTIAACQKZAAASBTIAACQKJABACBRIAMAQKJABgCARIEMAACJAhkAABIFMgAAJApkAABIFMgAAJAokAEAIFEgAwBAokAGAIBEgQwAAIkCGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgKbXWy70PAABwxXAFGQAAEgUyAAAkCmQAAEgUyAAAkCiQAQAgUSADAEDy/wFcthAJrDLbbwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def draw_gaussian(e1, e2, sigma, x_c, y_c, shape):\n", + " \n", + " # compute centered grid\n", + " ranges = np.array([np.arange(i) for i in shape])\n", + " x = np.outer(ranges[0] - x_c, np.ones(shape[1]))\n", + " y = np.outer(np.ones(shape[0]),ranges[1] - y_c)\n", + " \n", + " # shift it to match centroid\n", + " x1 = (1-e1/2)*x - e2/2*y\n", + " y1 = (1+e1/2)*y - e2/2*x\n", + " \n", + " # compute elliptical gaussian\n", + " return np.exp(-(x1 ** 2 + y1 ** 2) / (2 * sigma))\n", + "\n", + "w50 = draw_gaussian(0,0,50,48,58, (96,96))\n", + "w80 = draw_gaussian(0,0,80,48,58, (96,96))\n", + "\n", + "plt.figure(4,figsize=(10,7))\n", + "plt.subplot(121)\n", + "plt.imshow(w50,cmap=cmap)\n", + "plt.axis('off')\n", + "plt.title('Small Gaussian Window')\n", + "plt.subplot(122)\n", + "plt.imshow(w40*star_noisy,cmap=cmap)\n", + "plt.axis('off')\n", + "plt.title('Windowed Star Image')\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "plt.figure(5,figsize=(10,7))\n", + "plt.subplot(121)\n", + "plt.imshow(w80,cmap=cmap)\n", + "plt.axis('off')\n", + "plt.title('Big Gaussian Window')\n", + "plt.subplot(122)\n", + "plt.imshow(w80*star_noisy,cmap=cmap)\n", + "plt.axis('off')\n", + "plt.title('Windowed Star Image')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ellipticity of the star image with a small window: -0.0080+i0.0810\n", + "Ellipticity of the star image with a large window: -0.0070+i0.1150\n", + "The choice of the window induces a difference in ellipticity estimation of : 29.4% between the two images\n" + ] + } + ], + "source": [ + "#estimate ellipticities\n", + "e_w0 = get_ellipticity(w50*star_noisy)\n", + "e_w80 = get_ellipticity(w80*star_noisy)\n", + "print('Ellipticity of the star image with a small window: {:.4f}+i{:.4f}'.format(*e_w50))\n", + "print('Ellipticity of the star image with a large window: {:.4f}+i{:.4f}'.format(*e_w80))\n", + "\n", + "#compute deviation\n", + "e_w50_norm = np.linalg.norm(e_w50)\n", + "e_w80_norm = np.linalg.norm(e_w80)\n", + "deviation_w = np.abs((e_w50_norm-e_w80_norm)/e_w80_norm)*100\n", + "print('The choice of the window induces a difference in ellipticity estimation of : {:.1f}% between the two images'.format(deviation_w))" + ] + }, + { + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "The bias induced by the weighting window can be corrected, for more details see [MOMENT BASED METHOD] section." + ] } ], "metadata": { @@ -41,7 +286,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.8.6" } }, "nbformat": 4, diff --git a/shearbook/moments/star.npy b/shearbook/moments/star.npy new file mode 100644 index 0000000..79b403b Binary files /dev/null and b/shearbook/moments/star.npy differ diff --git a/weak_lensing_scheme.key b/weak_lensing_scheme.key new file mode 100755 index 0000000..6dc0e61 Binary files /dev/null and b/weak_lensing_scheme.key differ