{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "- title: Finding a planet by subtracting starlight\n", "- author: Joseph Long\n", "- date: 2020-09-24\n", "- tags: download-notebook\n", "- slug: klip-in-numpy\n", "- summary: Implementing the KLIP algorithm of Soummer _et al._ (2012) for starlight subtraction using basic NumPy. (Basically, the document I wish I had when I started working on high-contrast imaging.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When I started graduate school, I had to get up to speed on high contrast imaging. As is typical, I had a bunch of different papers to read that—together—represented the state of the art. Turning that into a complete conceptual picture that I could implement myself took some time. (And, of course, once you 'get it' you wonder why it took you so long...)\n", "\n", "This is not quite an [explorable explanation](https://en.wikipedia.org/wiki/Explorable_explanation), but it's something I wish I had when I started graduate school. If you'd like to, you should be able to follow along with nothing more than this notebook, an internet connection, Python, NumPy, matplotlib, and scikit-image.\n", "\n", "**Table of contents**\n", "\n", "- [Introduction](#Introduction)\n", "- [Get the data](#Get-the-data)\n", "- [The preliminaries](#The-preliminaries)\n", "- [The algorithm](#The-algorithm)\n", " - [Step 1: subtract image means](#Step-1:-subtract-image-means)\n", " - [Step 2: compute the Karhunen-Loève transform](#Step-2:-compute-the-Karhunen-Lo%C3%A8ve-transform)\n", " - [Step 3: truncate the basis set](#Step-3:-truncate-the-basis-set)\n", " - [Step 4: estimate the PSF from its projection onto eigenimages](#Step-4:-estimate-the-PSF-from-its-projection-onto-eigenimages)\n", " - [Step 5: calculate the final image](#Step-5:-calculate-the-final-image)\n", "- [Apply KLIP and ADI](#Apply-KLIP-and-ADI)\n", "- [Optimizing $K_\\text{klip}$ and measuring signal-to-noise](#Optimizing-$K_\\text{klip}$-and-measuring-signal-to-noise)\n", "- [The Case of the Missing Precision](#The-Case-of-the-Missing-Precision)\n", "- [Conclusion](#Conclusion)\n", "\n", "# Introduction\n", "\n", "To image planets around other stars, we are usually impeded by the brightness of the star relative to the planet. Therefore, clever schemes to remove starlight have been the cornerstone of exoplanet direct imaging. The state of the art in starlight subtraction advanced dramatically in 2012 with the publication in the Astrophysical Journal Letters of [\"Detection and Characterization of Exoplanets and Disks Using Projections on Karhunen-Loève Eigenimages\"](https://ui.adsabs.harvard.edu/abs/2012ApJ...755L..28S/abstract) by Rémi Soummer (my [old boss](https://www.stsci.edu/stsci-research/research-topics-and-programs/russell-b-makidon-optics-laboratory/meet-the-team)!), Laurent Pueyo (my old colleague!), and James Larkin (never had the pleasure!). \n", "\n", "In their paper, they explain how the eigenvectors of the covariance matrix of a reference set of star images may be used to estimate (and subtract) starlight in a new observation. Almost simultaneously, Adam Amara and Sascha P. Quanz published their paper [\"PYNPOINT: an image processing package for finding exoplanets\"](https://ui.adsabs.harvard.edu/abs/2012MNRAS.427..948A/abstract) in the Monthly Notices of the Royal Astronomical Society, describing a closely related technique.\n", "\n", "Though the specific notation and algorithm proposed differed, the two papers applied the method of [principal component analysis](https://en.wikipedia.org/wiki/Principal_component_analysis) to starlight subtraction. Principal component analysis goes by many names, as the following quote from the Wikipedia page demonstrates:\n", "\n", "> Depending on the field of application, it is also named the discrete Karhunen–Loève transform (KLT) in signal processing, the Hotelling transform in multivariate quality control, proper orthogonal decomposition (POD) in mechanical engineering, singular value decomposition (SVD) of X (Golub and Van Loan, 1983), eigenvalue decomposition (EVD) of XTX in linear algebra, factor analysis (for a discussion of the differences between PCA and factor analysis see Ch. 7 of Jolliffe's Principal Component Analysis), Eckart–Young theorem (Harman, 1960), or empirical orthogonal functions (EOF) in meteorological science, empirical eigenfunction decomposition (Sirovich, 1987), empirical component analysis (Lorenz, 1956), quasiharmonic modes (Brooks et al., 1988), spectral decomposition in noise and vibration, and empirical modal analysis in structural dynamics.\n", "\n", "These different names generally come with minor mathematical differences, such as whether means are subtracted along one dimension or another, or eigenimages are obtained from the eigenvectors of the image-to-image covariance matrix instead of from the singular vectors of the data matrix, but we will not attempt to cover all possible variants—just the description in Soummer *et al.*\n", "\n", "There are two key aspects to these algorithms:\n", "\n", " - identifying the eigenimages, and\n", " - truncation as noise suppression.\n", "\n", "As we will see, a set of images in an array of $K \\times n_\\text{pix} \\times n_\\text{pix}$ elements can be transformed into an array of $K \\times n_\\text{pix} \\times n_\\text{pix}$ **eigenimages**. These eigenimages are ordered such that the first image captures most of the variance in the dataset (i.e. the average image), with the successive components capturing variance in the residuals after previous components have been removed.\n", "\n", "The **truncation** step relies on the assumption (which can be checked empirically) that, in a set of $K$ observations contaminated with noise transformed into $K$ corresponding eigenimages, only the first $K_\\text{klip} < K$ of those eigenimages will capture useful information for subtracting starlight from a new image outside the reference set, while the rest capture only noise.\n", "\n", "The following code sample follows the algorithm as laid out in Soummer *et al.*, following their notation (within the constraints of Python).\n", "\n", "In the following, $I_\\psi$ means a PSF intensity pattern realized for a telescope state $\\psi$. We're trying to compute $\\hat{I}_\\psi$, an approximation of the PSF intensity pattern, in hopes that subtracting it will reveal a faint astronomical signal (which they call $A$). Our observed frames $T$ are combinations of the PSF and an astronomical signal: $T = I_\\psi + A$. If our approximation is good, we subtract it off like so: $T - \\hat{I}_\\psi = A$, recovering the astronomical signal of interest. (This is a simplification of what is stated in *Soummer et al.* since, as we will see, $A$ will not be recovered perfectly in general.)\n", "\n", "To begin, we import our libraries in the usual fashion:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import skimage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Get the data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[These data](https://github.com/carlgogo/VIP_extras/tree/master/datasets) of $\\beta$ Pictoris from [NACO](http://www.eso.org/sci/facilities/paranal/decommissioned/naco.html) were made available by Carlos Gomez Gonzalez, author of [VIP](https://github.com/vortex-exoplanet/VIP) (an image processing package implementing this and other algorithms). Unlike most exoplanet direct imaging data, they have [a planet](https://en.wikipedia.org/wiki/Beta_Pictoris_b) in them, which makes them great for a demo.\n", "\n", "Download [naco_betapic_preproc.npz](https://github.com/carlgogo/VIP_extras/raw/master/datasets/naco_betapic_preproc.npz) and save it to the current directory, or just run this cell:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " % Total % Received % Xferd Average Speed Time Time Time Current\n", " Dload Upload Total Spent Left Speed\n", "100 160 100 160 0 0 237 0 --:--:-- --:--:-- --:--:-- 0:-- 237\n", "100 2177k 100 2177k 0 0 770k 0 0:00:02 0:00:02 --:--:-- 1102k\n" ] } ], "source": [ "!curl -OL https://github.com/carlgogo/VIP_extras/raw/master/datasets/naco_betapic_preproc.npz" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load the NumPy archive and assign the data cube to frames_full." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "naco_betapic_preproc = np.load('./naco_betapic_preproc.npz')\n", "frames_full = naco_betapic_preproc['cube']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Peek at the mean frame. These data were taken with a [vortex coronagraph](https://www.eso.org/sci/publications/messenger/archive/no.152-jun13/messenger-no152-8-13.pdf), so don't be surprised that the star looks more like a donut." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAD4CAYAAACJ66HnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAviklEQVR4nO2da9Bl1Vnnf89537dpaEICaUGgieDYmhCsgFKIYUpR1DCaksyHODgVg05mmLKiohPLgF8y84EqqsbJ6FxMTU+iYhkTmRgrlBUTIyblZZSkSVLDLQQSIunQoYHQ0PTt7fc9z3xYa52zzz77svbtnH12P7+ut06ffVl7ndva//XclqgqhmEYQ2O07A4YhmF0gQ1uhmEMEhvcDMMYJDa4GYYxSGxwMwxjkKwv8mIioiCLvKRhnGYoqtroR/amN12jzz//YtSxDzzw5U+q6o1NrtcVCx3cQBB2LPaSq4QMREjrOP7Y8JqrnBMugwtjkp7fMBfZT2WzcRvPP/8i93/2f0Udu772I7sbX7AjFjy4uQ86fMjhQ4f8Dz7mi7EqX/LaJH/4fRwAQ/9C39J9LBq4agxqfaToO5j+vqePSf4OYkm20fr3X4Hx6n8uCx/cDMPoOwpbW8vuRGMWPrgl7y4xd5q2jukNdZRXpXOKju3oblzWv6L9aeWWPDatCFPPJb0/q70F0eR72vT72/r3X4EBZC6ZcjMMI4XatNQwjIFig1t1yhwKaePo4JwFRdOmVpwFWe330AmRR9b7k96Weq6a9BCuAcv/vhQ5CYq+22Xf/6x2i/bVwhwKhmEME5uW1qLMoZDetuw7cCYxsVk5RvDFkr5m+nmPvsAN3h/RxNc4LwylQTxdrT7VdDCUff+L2hWkHe2mimybt9QwjCFiyu00Jebun2cnaqzgwvlNvnw9+uK2oWhj2ljRYOE8O1qWna7Fi8LYQkEMwxgcZnMbJnUCQosCTxtR1Eaegou5bp1A36zrNenDgn48dT6HHqm8qra7dmxu2OBmGMYQUTCHQjNiEucXTp27dlFsViMF16e7Z1t9WfBrKknhmtnWA+rEuQXaTZw3m5thGIPDbG6N6Y1aK2NhMWppG1YfYuPaPifv2AUl9Rd9lguOhcvsQoM4t1axwc0wjMGhIDa4GYYxPNRKHtUhmTi/Miw1hQqywy8Wfe22r7v6ymAZFCXbt3gRK1ZpGMYQUfOW1qF3qm3hzoI6xzbro1R4jTpnTG+Q0O6vO99m3XbTKjJLVdZQhFWcDoEeJeRbEG82ptwMw5jHBrcBUFTDvxFtJLivNtmKrQsW+B6XFM5cWDcKila20bo5FAzDGB4KbG0vuxeNGfbg1llCexHV7WVlNrEsBVTFjlaFrtoNNFNzMe9tmXe3oZ0u7/0pssF1YJfr1nZtGQqGYQwRcyisAAuzh6zQAiwVCUqrSNHFHBMo9qB2Sfp6HV0/ZoGbVcBCQQzDGB66mgNyChvcDMOYxUoe9YA8Q+3CnQblxBnqs1enKj6329c6f+12nBsx56Snu3FT2bx2q/SxhUDgGPqqjsxbahjGMDFv6fLJu/N1tqp7/XPjwjliUp+WWeutjLy0qHo/lMU7H1q8TlEoyNIKMFRgANPSFXiXDcNYKIobiGP+ShCRS0Tk0yLyqIg8LCK3+e3nicinRORx/3hu4pw7ROQJEXlMRN6U2P79IvKg3/ffRKQw2G+1lVsnd8R2xvt6wbDV7XJdB92micvKqZP4n/9DaeM1NrPXZVGSxJ/1s8tL0O+d7a3VqiBbwLtU9fMi8grgARH5FPDzwH2qepeI3A7cDrxbRC4HbgZeD1wE/JWIfLeqbgPvA24F/hH4OHAj8Bd5FzblZhjGPGON+ytBVQ+q6uf9/48AjwIXAzcBd/vD7gbe4v9/E/BhVT2pqk8CTwDXiMiFwDmq+g+qqsAfJs7JZLWVW1kSc6U7/mqN852lX6XeB00pqqzrtmMTK7IvNqd9+12ePbFC+71TbJ5q3tLdIrI/8Xyfqu7LOlBELgWuAu4HLlDVg+AGQBE53x92MU6ZBQ74baf8/9Pbc1ntwc0wjPbRStPS51T16rKDRORs4E+BX1XVlwrMZVk7tGB7LssZ3OrYGorO6djulL7rd2dPy75u5r4FKc2i62hUHF4c7Re0zCYmjawaDVbzyupLX9Rci6EgIrKBG9g+qKof9ZufEZELvWq7EDjktx8ALkmcvgd42m/fk7E9l6hvj4j8mvd0PCQiHxKRnUXeDsMwVhzVuL8SvEfzA8CjqvrexK57gVv8/28BPpbYfrOInCEilwF7gc/6KewREbnWt/n2xDmZlA5uInIx8CvA1ap6BbCG82bcjvN27AXu88/jiHQjR59Tp70ql9YxqmNERjVUyoiF2PNktPD4KWG0MAVZnfj3PXy+7TCmsq0wK7RiCZ/ntD+05lAArgN+DvhREfmi//tJ4C7gx0XkceDH/XNU9WHgHuAR4BPAO72nFOAXgffjnAxfocBTCvHT0nXgTBE5BZyFk4N3ANf7/XcDnwHeHdmeYRi9RVtLv1LVvyPbXgZwQ845dwJ3ZmzfD1wRe+3SW4OqfgP4LeAp4CDwoqr+JSlvB3B+1vkicquI7HceldWPejaMwdOuclsapcrN29JuAi4DDgP/R0TeFnsB7xbe59oa5b8bZQG5naWudBUUOyptL0zp0uEW6f2F1Fm5qQoRqzy1MjWtUVg2zgkRU2OurWuXXS/nnJjPrOj30YVJpucDVwwx09IfA55U1WcBROSjwBvJ93YYhrHKVAsF6S0xg9tTwLUichZwHDdP3g8cxXk57mLW21GPpa0qNJs2k7wj56uu8oTw9LkTdZPRprT4PcpSUXnKMK7BCGW4pPCFOip7ZZLwy34PHb8OPR0GN1W9X0Q+Anwelyf2Bdw082zgHhF5B24AfGuXHTUMY4GcLkv7qep7gPekNp8kx9vRCktaCT5bDdQvgFio2FL75hRWhYDlIrtX1+Ea2sJCTEXqdS4FrMnryehr+e+4vFzVvCKsYXubuUCOzXkRKLDVk2DiBlj6lWEYs5xGNrflsKSV4ONsbrP7q9hxKqmOmNWkllwCCagV4VPFDlj6nhWtT9u1bSqq/QZ9WFY6lg1uhmEMkdPCobAwWik8uZgSR3Oe0BoKK/O8nO9TlsrJU2x17FF11d7cAi4RCmVyTjrGr4oXtqi/Jd7dTNteyg4XZz+rsnhN+vwaBSMmnVuAkrPVrwzDGCw2uBmGMThUYdu8pe0RG6RYOIUqmk5UcdXntVsehjGpcZZT0Xbm2BamGDHT0bJ6dDE1zrKOmZuSaXgoD27OnZ72nib9XA3HQlgfZtXpz+BmGEY/MJtbx7TuUIhrr1oQb8b56WNjnA05x6TVTdY5RaEgee3Op4atFfQtf1+igVlifhcR57QRHDx/3fkQkbR6LA7QzQ4lKg/qTZ4baOBYmLmYJc5n0d/BzTCMpWHT0rpUCbCsFcxbbCNbCHlqjHnFkNtEBfU6UXJFJZYmam8tc3vy/2nlUWQfnLRXGlKRPCdbLXWe2F6h/ey+lAd/59PCa1tUKMiWKTfDMIaGqgXxdkKj4ovVFdui05XaTnBPK7Ys+1metzRL7U0UlYxnzim6dtpbWqnwZPo3JPP7FuVJbaImO7e5TS60oJJTNi01DGOQrL5w6+Hg1lk58ZzLNVqLtN+UxbNleVzT28ZsVb5OWsHFFCOoYp9bFHXsf4NIpFfLLTUMY6jYtLQhHS/2EtWFwj6UKJ+G3tjYOLSYNoKtrUqcW+h/8riRbKSOde2Ox6dK+xLsdOOxV3s624bbtD3T38J4twzl1xoRC6xk29yKS2YtPJEeZmc7bQguBS0X7L3HlJthGDNY+pVhGMNEsWlpbWJuC5WmrPWnh8UOhfaCgZe5xkEdZ8nkHD/NGY3cdFW1YCXyyTS0wmstCB9pNB0tO7dlaVKtrx29rhZf0wDWhzHlZhjGPDYtbZPOK++mLhd1vRrtR667UHVf3rFpR0KWQyHLcVDGKISCBEVVoLDSgb+i/twIq3RQPBrOYeq4CO2sQjmkuCDeQNHnkJ4tRKQftj0S2bTUMIyhMi6wPqwK/RvcKgXx5t1eys9N31m7CuKtYk+LCSiupO78tYO9LB0Mm2U/mypApwjXvDLcDoqqwBYzXYfVPWT1NF1CKfQhqDxNVICVoARDSwtOx6pDe2ErOe0sIshdSUj21aV/g5thGEvFQkGaUGX1okLaDtptj5j1RMuOyT4nu1xRjM1tTWY/7m2fWpU8Z80H8Y78sROb3iQda2oTG3t1MUq9ju1gayuwz01XzErZqkbTX5X6uVGwvZG2/xUwV+hyQb/WKt+vWsn1We3PBCRHX74AQTupFLpYTLkZhjGLmnKrT+t2g/h4tHK7VoU7b85iMGX72iRqzdSUHW1ii/MqLdnGuuwEYEPOnGkj2MQ29Vhim1NxaRtY0JZBwYXj/MXxJ9UmakGbvHVRM3619Wx4xWlYMWS9jvw1UxcT3wbe5LZtys0wjKGhoOZQqMkSEuab2NrqJMy3sfJ7UeHJ+SKVWTY3d/5oYnM7A4B1/xjsa2tMk+V34BTbGeoeN9gBwCk2ATg+Ojo59iQvA7ClJwHYDkpOvBc2Q1CUqaTZwpkhhm+2GGZRcvrUlldhucGSkkrJPpUrqyylFa/y8l9bQRut29zay1AQkd8D3gwcUtUr/Lb/CPw74Fl/2G+q6sf9vjuAdwDbwK+o6if99u8H/gA4E/g4cJtqcS8XY103DGOlUJWovwj+ALgxY/t/VdUr/V8Y2C4HbgZe78/5XZnGDr0PuBXY6/+y2pzBBjfDMObQsUT9lbaj+jfAtyIvexPwYVU9qapPAk8A14jIhcA5qvoPXq39IfCWssZWPHG+TUdCus24dtugrE+Z4R0509FkkOy6uCnl2shNQ3fIWQCcIWe756kpKMCZ6hwKZxHOce0F/X9svDk59iVx09KX5TAAm+KcDWGaKiM/LU44FMK+SdpVhala3vamgbNVHArl14xfMSuzL5XaTZ/cjmNBdSGJ878kIm8H9gPvUtUXgIuBf0wcc8BvO+X/n95eiCk3wzBSCNvbo6g/YLeI7E/83RpxgfcB/wy4EjgI/JfJhefRgu2FRCk3EXkV8H7gCt/ovwEeA/4EuBT4GvAzfvSNaLCtMbV6onmVNqKdEDVSorL3ZQfojkbTj2kSXJtKYwrtrnuVBlMHQlqx7dJX+Ue3fad3GgCc5VO1dq2765y55r5XGyP3eHx76nw4vOnaf27sHl8aHQbghFd0wWER1BrMJ9Vvjb2SIz8VLKRotaLUClacDywrNGS57SeoptyeU9WrKzWv+kz4v4j8b+DP/dMDwCWJQ/cAT/vtezK2FxL7q/wd4BOq+lrgDcCjwO3Afaq6F7jPPzcMY8VRWnUozOFtaIF/CTzk/38vcLOInCEil+EcB59V1YPAERG5VkQEeDvwsbLrlCo3ETkH+CHg5wFUdRPYFJGbgOv9YXcDnwHeXfrKarP8GfSc6ipIqYoJH8lbJT6t2EaSVG6zgbchzCNs3zE6e3JsUGzBtrZT3b5z9BUAvMIru7PWpu3v8ErtrHXX7tkb7vlOLxTHOlWMu/wxG8fPcY/brp0XfN9OyNG51xfeh1Mc8/tiEvLX0hvcKVXWE02HhBSUDqqn5OorquLXUaFgaltrKEBr6Vci8iHcOLFbRA4A7wGuF5Ercb39GvDv3TX1YRG5B3gE2ALeqdPqDr/INBTkL/xfITHT0u/ExaP8voi8AXgAuA24wI+oqOpBETk/58XdinPhGoaxIoxbGtxU9WczNn+g4Pg7gTsztu/HmcWiiRnc1oHvA35ZVe8Xkd+hwhRUVfcB+wBERg3uK+2V/G5MAxtb0epU+Z7QtbljJoG53sZ2hjg1tlPOmRwbvKBnsNM/d8ee7RXb2cGutp5UVrOPgQ1/yHpix5qEJ/5rdNwpxZFP3XnJe1yDgnP/d/a4kIA/lq2Zx+SyS+n3Iy+1ahI0zHTlrbzV47OUTa6a66gUeTXb4QJtbR5VGUT6Vcyv9ABwQFXv988/ghvsnglzZ/94qJsuGoaxaMYqUX99plS5qeo3ReTrIvI9qvoYcANuTvwIcAtwl38sNfDNN564G0WroSXEoc0tajJ7py/ygMYk0GelWeWRtrFteLtaUGxne08oTJVaKDh5hv+4d4xcG8G+Vnw9R1BsOxNd9Q5UNn0w54ltv8bppk++98+P+dg5gKM+IT+ol7FPrdqWVOoWCWUWSit52+O02GZYJzXRpyDQwmeUVmoRQQVpG1vys2tSKLOtuLxFcDqVPPpl4IMisgP4KvALuO/9PSLyDuAp4K3ddNEwjEWitGdzWyZRg5uqfhHIimW5odXeGIaxfPT0Um7dsOAVr2KoUz0kOVWZcyB01N/Qz1Dh40zd5R+nU8AN//Gu+T5sTKqEuC9uWK5gLTHV3Onnn9+20z1esNPN2c7bcNPF9YRP6IVNPw31gbnHt3y7PlzkDD9HfMV4GiR8bOydHD7wd+RTtI7617PpK424dvy6CqmpYIgOmGxPvMWTKWrqt5lb542C6WJBDbj5c9qo/NEf+tuzeKyem2EYMyjC9rgHUQkNWfzgNlN7qopDofqbXazCytuTEDybd4eNqCjcpCJv1ipVwQC/TlidKqxsNVtlF6YKbT04IYJi87k1p8bucdfG9JxLz3bHXPnK4wC87vznATjn1ScA2DoxPfapb54LwIMvOGfG5tivlOWVXHA+JP0WJ72TYadXfev++bNe5R1NyMjN8cu+PR82ErEOajplK/e4AidBUX23stpvQ8GmpYZhDJJx91VBOmcFbW5NyQ+qhZTC6tgmMrG55KwrkGWbmZYKcsokrCc6fZwqFgk2t0nIimPLt7HTb7/wrOld+gfPc2rpmusOArDjh77D7Tj/PPf4rcOTY1/5ha8DcN7fu1SqjacvAOBLXvEe9ZlV6zNvsbvWGWtOceoxlxJ2Up0K3JJpkv2WeLUYkuuDra2gXFKsHauS8sqbbVSkno1tcaWOJs2ZQ8EwjKEyzgwIXC0WP7jpuKZiazf9Kp3ilH1QvrdsZn/RddJJ2MlSOw2Ua1q9BFvcKHHHDba2aZqUI9jezvNGse995bSY5FXf6yrJ7Pix73Ltv9Y/7nLeWC6ZFqscXfztbtOFjwLwwx/zVWj+6SIAnng5BN1m9N/34cgp14ezTrpg5JdkWrKpSnBzGXkpXK4v2WlRy7Gv9SPNcAHFKjvHlJthGDMowpYud3Btg+V6SwupnzCcrYhKrlvUr5KE6sxyRhErzudeLqUqR6ON+X2p9kb+y7iR+Eh3+BiyjZDq5DXUORvumNf6HPvrXjOt+7frh31xlzOdgpLP/T/3/AVfvmjPqyfHjq94rfvPT7wRgG8f/z0AV33Ylcx/eWs3AM+eTKokR/CghiKYQWWuJ1biGs0VEkgpudBYQpjOJcinyiN1H2PWNNE98vys/rda8qiddpaJKTfDMGY4rdKvlkPMuo3Zaih7MZhZW0aerS0qSbrAXldWxijL5hYUSbo4ZVBsyWKVk7VG/eO6Lys0KWuU+Eg3Jtf25xKyD9z2q17lPKPn/8hUEekVzsYmX/4aAAf/yKmwrz7rvJmvu+SpybHn/ZKzw21f+wPunCv3AnDh3z4AwLnfcue8cGr6npzwztxtDY/5EiG8L+G1TvoYbGESEuenNsM5hZyh7uaIWNt0vm+jmXYrFc4Mly08p0TBdRxpoOZQMAxjcKjFuRmGMUAUmWSZrDIrNLi16yKfS4sKoRU11HiWsyDX8J+aYsF0OppXiXeUeM3TbX566h831K9RmtGXTZ9NHqapZ/pPfc+rXnJtXXnZ5Njxd7zGXeehJwD4p+deBcDnX3DBthuj6TTpmqfcIkZy+YvucWs25Snc/U8lZlZhWnr0lNt5bMvtPI4LMTnFNIh3/jWHGnB+h29/5j3OURxhipu1ylYdyqewWSaRKiw3vcuUm2EYg8RsbnWoHcQb71CYXKogYDad+pRWcDPk9Hdy907c5fJCDYrU3eR5zusZJ15z+ohRass40ZlTvg+nUn0Jd+Wt7dHsBkB3+PJEF7kwjr0XPQ7Atpe0l16QWJp2hztGDnzDnfvIkwA8902n8o5sObW0mbj8MZ/7/pKXcy9tOcV2zK+tcIoTk2NDwnyZSuoqrKPounlrNPS5jFEszlu67F40x5SbYRhzWCjIQmnrjjhru8tVcCSVWU7ZIpm9e/snblduAUTmjp08nYQ2bM08TwbxjoPdKaweFdKvCt6fkG419mEXL3l715OHXwnAZf/38emxu12CvH77twHw6ltc6aM3Pn/Edf/ci6YNv9pFAYewkSN/cxiAJ15wi4MfPhXWVpiecnzL92HTbTzsi1MeE9f+KT0+fa1BufmE+fAaJ2snZCi7POUUbG2TtLVkMYJUu3OKLaPN+QKa6WPyP4+4dUqXywCE2yoNboZhLAJV2DLlVoM8e1tE4cd54j2o84G9xQous40cD2tWn3MLIM4U6JwcPL+PpKKYL1aZPiZLuYVVrkapFxXU0yNHXODvaz69a7Jv74797pzrXGrV+MrL3Q6fOK/Hpspq9OWvALD5OVce6ctfcWrvqWMusPhYKDue6Npx/+TlsbO1nRBXLinY2rK8menXNqe0spRVqjTU3P4ItdeMggDyCoG/UdT67ZQ0aYObYRhDQ+nL5LgZK+QtDdQ7NzZxvWlJorQHdbI+p86uhJ7Vfrj2JO6t4LWO/bFhvc9tv97ndqL/QbGFMuLpe/E3jrktf/vN3dM+fMKVFf+ulx8EYP06lzCv57uEefnmc5NjN//aKbdH7nfnP/SiW/X++U13veAlPbE9teCEuLaThJXm/WuelEnPL3MU3p+tbaceY4p5zm3vaBX5KeX2tNaVYllprhqYt9QwjAEiFudmGMbwsDi3utSu51Z9itik0u0yyOuvyPxULYRBbKlLVzohbvp4XM+aHLNLQwiJr/U2mr0bv+xDQr7y8vS6I3HTz6N/7c697EsuMHfHLpeOdfjQmZNjH3rmEgAee9k5EF7YdO37mScnQ6rV1vSXcnTbTUdPptKsqqwOViVgtkq1jpjVtZpU5e0swLeDdrfNoWAYxtBQqwpSkzyHwty2Rflr8pVhdl249pi2vzbzPJC1Bmc4ImiMk3pk5ti1RO23dZ9eNdp0ais4FsKaCuHe/EJCRD3u34fDp1yA7xf8mqTBJ3AoUVX30PGQ/O4ew/Kna/4xhJw8f3Jab+1FdaEfm35lq7Bq1yS8I+EIyFNJ6dSn2YDcnNCPFlbMWgxNK/m2g9ncDMMYJKbcOmXRd7B2SyoFYtTfRG2kvlBj5sNH8HahkH41p/ZGCftZSNXysuvk2NnGdvpV3TdGIWQjkWzv6wm9uBnu3O7xqDekfevkNLH9hO/Lhleeu9Z8BeFwjrevvaBHJ+e8NDrsruNtbkG5BdvhuEBFlac8xe3L2x+z2tXpsOK8xbkZhjFIFHMo1CPafhUTHBlf8mhy+RbsZ7lpWDCXTD+XdpWROD9Vd7PniM57SUP/Q6Bv8Jpu+1XZN33pIIAjXlGF9KVT6uxnZ265skY7vH1uI/GeBHtceAxrHBz3KuzlREmiTa+61tW1c2Jrw1/P2+J8atXLoxcT57gA3O2U8tzy7YbX4/aVp0655wk7XTqFLcajWkGn9FKxWRBvJqsVK2EYxkLQyL8yROT3ROSQiDyU2HaeiHxKRB73j+cm9t0hIk+IyGMi8qbE9u8XkQf9vv8mIqXSsj/e0miajcd1PKB558wpOJhLYs5VcG5j/r7kuYnr5sWDbXv719Z46vo8ERSTv862uGNOqSsmeYY6G9yGJlbXmqwfOnudU4SyQwn7nE/92vSq69ikj04xHhNXxvwkUzWZLkAZno9Tj0mmqVP5ZYvK6KXiKqSCzbntlDJaref2B8D/AP4wse124D5VvUtEbvfP3y0ilwM3A68HLgL+SkS+W90H/z7gVuAfgY8DNwJ/UXRhU26GYczRlnJT1b8BvpXafBNwt///3cBbEts/rKonVfVJ4AngGhG5EDhHVf9BVRU3UL6FEnpsc1txGpShqaUug11tPF/scXs0+zj2SfZbOAW3nkhWX/MqLmxbI9jt1J8zjeIPdrOQHbE9KaDpVV5Qk1mlw+dKn88rtqBGp6puK/PcLNJxbXmquLiRKt7Y+p7b3lEtiHe3iOxPPN+nqvtKzrlAVQ8CqOpBETnfb78Yp8wCB/y2U/7/6e2FRA9u4iJN9wPfUNU3i8h5wJ8AlwJfA35GVV/Ib8EwjFXAeUujD39OVa9u6dJZc2Et2F5IFVlxG/Bo4nmYN+8F7vPPDcNYeYRx5F9NnvFTTfzjIb/9AHBJ4rg9wNN++56M7YVEKTcR2QP8FHAn8B/85puA6/3/7wY+A7w7pr1mJCX+YpLpy87JmuLkGf6LnAN1nBzpayend9s+UHakmzPnhFloqKW2wRmTc0Lgb3gM66Ku+bCU4ERw//fhJ96VcMqnVqUDcZOhGtOp5ey6COnXlXwt4ZzJ2hKpmm2dTfsKQiyGtNpVFtptKMi9wC3AXf7xY4ntfywi78U5FPYCn1XVbRE5IiLXAvcDbwf+e9lFYqelvw38BvCKxLa8efMMInIrzstBYQ1vwzB6QZsZCiLyIZwI2i0iB4D34Aa1e0TkHcBTwFsBVPVhEbkHeASXPv1OnQY7/iLO83omzkta6CmFiMFNRN4MHFLVB0Tk+iovzHd4H7DPtTVaydDAmDVIJ/vSii29kn3GGqd518kipF2FFefDClmF6wj4Fec3U9cN527JVLmlw1OCcluXEKA7VWEhIDcothCGMp8mlUz8Lzbwj8dTZZh2IIxznBFZrzl3f9bPtob66k6x9UMJthXEq6o/m7Prhpzj78TNENPb9wNXVLl2jHK7DvhpEflJYCdwjoj8EX7e7FVbct5sGMaKs5IqJEXp4KaqdwB3AHjl9uuq+jYR+c9kz5uNSGIDS2fWdVCv9lIxnkGNzdj05pShU1CTcBG/OxmqkSYE5K57dZdUpiHZPSi28LwolGU+RGPWbpYsGDmxsRWsOTrXfoXk+ljiPqf8oNtq9rnY0l+J4yazgnZUn+rsimWrSpM4t8x5s2EYq88AxrZqg5uqfgbnFUVVnydn3rx44ssVlQfIVk/Qr0LW6vSV8HfpoGpkNJtcn/W60uldQR0FBbedKvkNU9tY8JpujZxyW5ON6TEpW1jaE1q0klWwrU1tbfMBumWKrVHCewX1F0dbQbyxxya/R12kX7Xa5FKwkkeGYcwxgLFtFRPnk9S/Y8WlOMX1M2pxkyxvZp3ImPCtCwpOU97TjHVX08n7E9UU1jotuNw4lE3yHtexzC+isp1SbFXKDOXFrs1sa6DYYsgsgNBXOlhdPgtTboZhDA4rVmkYxmAx5dY74h0LvSA9fUtPNYqmsqGKb5je+eDe5IwtHR6Srg8XZegOU04/hY2agacCdZPTvfS+olSqNqejMWsfxLQf6xxY5bSs2HJGfWdgg5thGI2xdUtrspB6bm2FhmS30/V6pjGu/ek6n+OwAUgoOObXYshTbDGvY6IQE0G2ReuHJp/HrBVaJTC3bcpUXfcqrEL7RWv+tthPHYB2M+VmGMYMFudWl2QoSMaqUe0Q31ac+ipWgsk7fqNwgigb2GwowNyqWhlfynSoQ1Fwbd77MVFsGTa9uevVWOMgfW7dY6PVdPK4nGvG9SX7mGQ/yttJ9rmG+upAWVYoVtlbTLkZhjFHx/XcFsJyBresO82CghMDndnLJhco93w2Iq+UUiI8KctbmTwmS8GVqYwx8za39LlFyq3M7jebflU9CT5a+TX+PIrPb6L6oujQDmgrzhuGMVjM5tZbyr2lnXs803Tg0UqS9prOXtvv0pAMvz6zvUoBzaz0rrTqykqlSnS0MrnruVawb7bncV2ypsma4bQ+K7BpqWEYA6Ti6le9pT+DWycKKnlHW94q9zPUuOMWemMjbJVplTpnD0qmEeZ8qbNU2KS0ka5lHptpd8pJWaySbJ/33DU/Kj6m4Dr58W7NPLitsqCZhtncDMMYHIqiA5iX2uBmGMYc5lAYAMWOhQWFi5RMZeoazvNCNXIdDInzywJ0Ux3M3LeM5PFl1Wbrptpu1oUW854OYGyzwc0wjFks/Wol6PbunaX6JsphElQb0YcSBRe1ov2kLxkqrCTxPysZPo+sUJBQTTeGOkn76b6l3486lZBjQkOmDpjk1uJ0q2IFV0F1LbNkksK22dwMwxgaptzaZsHpV22SpWamO+u/rjqBp5nlhVIKLibgN8boEhRfafhF5mWqrOVZTKNVsFqidftix0HfZQxAuPVocDMMozeMB+BS6N/g1qqCiw/izfaapu+a8cUvp8np1e/AReqi1L5UUEYqJgg2KwG/jFoJ7iVe2STVVmyvTq9KgvdkBmPKzTCMwWFVQbpiJYpWxl8nRoU1sRkVKrm0CshTKBne3qDgmqimqNdVQSFUSaUqbSvKqzlKPc+nNXW5ZFubuzZsD8Cj0L/BzTCMpeKUmw1uhmEMELO5tUlajjeansavZBVoVt9t3nERE6xaZ7WlWv0rm+IUOSFaXDs089IttFdpWp/5XuSdH+EY6Wr6uMTpqaKm3AzDGCam3Nqk1TvUMn092QqxszCGKgGzealDyb6m12aI6sSSAk0j1kPNJ/+4pYaG9CAsRYGtHvSjKf0Z3AzD6A1DWJS5v7lOOu7oLjam7K4d/rphTFkfikj2L7aPzoJS8Jr9/plj8t7/sD35F3n9mL8q7Uy2V3o/6r/31Rhn/K0OWb2v+4pE5Gsi8qCIfFFE9vtt54nIp0Tkcf94buL4O0TkCRF5TETeVPc19HdwMwxjKYRQkJi/CvyIql6pqlf757cD96nqXuA+/xwRuRy4GXg9cCPwuyKSv4p4AaWDm4hcIiKfFpFHReRhEbnNb88deVtBRu6vMwVXTp5Kqqbsiu51sffB/GMqKbiI15Gr4Ao+hyZqrFJbmv03S/p9Kn+P21Hqq6vS5nFlxmP+GnATcLf//93AWxLbP6yqJ1X1SeAJ4Jo6F4hRblvAu1T1dcC1wDv96Jo58hqGsfpUUG67RWR/4u/WjOYU+EsReSCx/wJVPQjgH8/32y8Gvp4494DfVplSh4K/cOjEERF51F/sJuB6f9jdwGeAd9fpRPaFF3X3K45/S5JXRjs79qxO/7tKJ6p+bF6aV6HtrkKsYPceyeIYwsy4w2V5CHvmmVRgO2vN2WyeS0w187hOVZ8WkfOBT4nIlwqOzfLT15KIlbylInIpcBVwP6mR13c865xbgazR3DCMXtJuEK+qPu0fD4nIn+Gmmc+IyIV+7LgQOOQPPwBckjh9D/B0netGOxRE5GzgT4FfVdWXYs9T1X2qerUb3asET3VNmW2mmTezbH99r2wb/qvyPs39Zdi9mrRXj3j/Xdl12+9LDZZoTy6iTYeCiOwSkVeE/wM/ATwE3Avc4g+7BfiY//+9wM0icoaIXAbsBT5b53VEKTcR2cANbB9U1Y/6zXkjr2EYK864PafIBcCfiQi48eaPVfUTIvI54B4ReQfwFPBWAFV9WETuAR7B2fvfqWH174qUDm7ievUB4FFVfW9iVxh572J25DUMY6VRVNoZ3FT1q8AbMrY/D9yQc86dwJ1Nrx2j3K4Dfg54UES+6Lf9Jm5Qmxt5OyNGvtdKeo+v2VUltaqOgT9mrdByY33yuu2FMca8Rg1rKkiTxJdmP6p2nQJZbQ0/NPS0KXmkqn9HvrEsc+Q1DGOVUbaJX66xrwwrtzR9125UXbdoX3p7O3fzqBCNiFJKUxZlrHZ9mFds/TOWZ5NW5HnP0/+vSA+dB1koMG5pWrpMhjW4GYbRCi06FJbGag9ueQX9stRMK0Uw8+jGzhVDsyKb1dvPtv/l/RCq2DPzr52m2XoFWXbTvPdu9X/g9VAb3AzDGB5K+xWXl8FqD255d+6l2jbK0rmyVF58ClgeWSqmjprLU0P1FGJ5qlYdqp2b9962/B0p6tOS1yCtjrLNqWV3ojGrPbgZhtE6ippDYaH0YT3HSlSJm+uGKvFzw2eJr3UF3+dxfOJ8b1mdwc0wjAWhZnNbKG3d/WrFwq0mXXhQY9Re+rqLVYjN7ZeV6IEqC+sdSEuFKRQY10vn7BWrM7gZhrEgTLkZhjFIzFu6mvTaMRHTp/jpVlEAbixro52uZ7oZ3cbipqF1UudaokffnzAdbWs5PotzMwxjoCg1S6j1itNvcKuThtVp6lZVqiftN1FS2+MTtc9tnyWpiR6ptCIEaa1QkaVfGYYxQBS1OLcVpkrxyyp2ujx11yv11zZthF8Uvbcdv2fpz2ZFlFqaVm1uK/oeJDl9BzfDMLJRZVvNWzo8ksoqfUePKbGUp9A6U2xtF86sk/jfNS2qiCJFsqJqJRnE25Z2M2+pYRiDw6alQyX5oZYptbA964uwtDSvKoua1CkdlLUtryhl0zLvdWhWILOUntjlkva1ttKuZls3h4JhGAPElJthGINDUbbVVr8aNm1W+s06Z+FhIWVrHWQd01Vduq6UwemRdtX+VHQWU26GYQwPtfSr04eujMixgb2NVV9Zv4ucBU1Ywt1/IAG5gXRgblGSfJtqzkJBDMMYIGrT0tOGrj7otLrIe74UOg6p6JoV/XHmVdVNK7YubW4W52YYxkBRxuYtNRbOIBPvGzIAlREoU2hZiq2thPmZNgfwntrgZhhGCmVlzREJFj64KZrp8ek6bqdTshLn22ovTV77Mee06mFt2n6NNmKKEgxAccBibWxzqCk3wzAGiK2hUJPkHWil1VqSrtZUbXpuulRTV1SJ08s7to5iW3F1kaXKYuPauv3tDCMUxKzThmGkUFRPRf3FICI3ishjIvKEiNzececn2OBmGEYG48i/YkRkDfifwL8ALgd+VkQu76bPszQa3OqMyF24rVcCGWVPq9qW/+E6WdO8Lq6XbDf9V9SnsnOSx4Rzi45dIdT/E/8vi/S+9HNN/Ouih7mfT/XP4BrgCVX9qqpuAh8Gbuqg03PUHtyWOSIbhtEtGvkvgouBryeeH/DbOqeJQ2EyIgOISBiRHyk6qaohdCmu8C5YlNKocp2mIRxl12qyRmyV66wI6e9y1ne7bPX4ot9Buyou+j3fLSL7E8/3qeq+xPOsDi9k+tZkcMsakX8gfZCI3Arc2uA6hmEslEre0udU9eqC/QeASxLP9wBP1+1ZFZoMblEjsh/F9wGIyLPKyaPAc1UvtkRL3W5q9HdJxPU1/WZ2/ebmtz+89zaCJm9/xLHfUaUvOXwStnZHHlv2nnwO2CsilwHfAG4G/nWTzsXSZHCrPCKr6reJyP6Skb5XrFJ/V6mvsFr9XaW+NkVVb2yxrS0R+SXgk8Aa8Huq+nBb7RfRZHBb2ohsGMbqoKofBz6+6OvWHtyWOSIbhmGU0Sj9quaIvK/8kF6xSv1dpb7CavV3lfpqAKJ6mgbVGoYxaCz9yjCMQWKDm2EYg2Rhg9uyKgPEIiKXiMinReRREXlYRG7z288TkU+JyOP+8dxl9zUgImsi8gUR+XP/vM99fZWIfEREvuTf4x/sa39F5Nf8d+AhEfmQiOzsa1+NfBYyuK1IHuoW8C5VfR1wLfBO38fbgftUdS9wn3/eF24DHk0873Nffwf4hKq+FngDrt+966+IXAz8CnC1ql6BiwS4mR721ShmUcptaZUBYlHVg6r6ef//I7gf38W4ft7tD7sbeMtSOphCRPYAPwW8P7G5r309B/gh4AMAqrqpqofpaX9xUQRnisg6cBYuOL2vfTVyWNTgtrTKAHUQkUuBq4D7gQtU9SC4ARA4f4ldS/LbwG8wm+Hc175+J/As8Pt+Gv1+EdlFD/urqt8Afgt4CjgIvKiqf0kP+2oUs6jBbWmVAaoiImcDfwr8qqq+tOz+ZCEibwYOqeoDy+5LJOvA9wHvU9WrgKP0dFrnbWk3AZcBFwG7RORty+2VUYdFDW5LqwxQBRHZwA1sH1TVj/rNz4jIhX7/hcChZfUvwXXAT4vI13BT/B8VkT+in30F9/kfUNX7/fOP4Aa7Pvb3x4AnVfVZdXW0Pwq8kX721ShgUYPbJA9VRHbgDLT3LujaUYiI4GxCj6rqexO77gVu8f+/BfjYovuWRlXvUNU9qnop7r38a1V9Gz3sK4CqfhP4uoh8j990A67uXx/7+xRwrYic5b8TN+Dsr33sq1HAwjIUROQncXaikId650IuHImI/HPgb4EHmdqxfhNnd7sHeA3ui/9WVf3WUjqZgYhcD/y6qr5ZRF5NT/sqIlfinB87gK8Cv4C7ufauvyLyn4B/hfOgfwH4t8DZ9LCvRj6WfmUYxiCxDAXDMAaJDW6GYQwSG9wMwxgkNrgZhjFIbHAzDGOQ2OBmGMYgscHNMIxB8v8BG0HF8tsvQR8AAAAASUVORK5CYII=\n", "text/plain": [ "