From e45973a2b1d98cf84def7fe220ac2da06b1638e3 Mon Sep 17 00:00:00 2001 From: Wei Changfeng <2681968849@qq.com> Date: Fri, 30 Aug 2019 02:15:53 +0100 Subject: [PATCH] Delete plot.ipynb --- plot.ipynb | 366 ----------------------------------------------------- 1 file changed, 366 deletions(-) delete mode 100644 plot.ipynb diff --git a/plot.ipynb b/plot.ipynb deleted file mode 100644 index f5e6f03..0000000 --- a/plot.ipynb +++ /dev/null @@ -1,366 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The purpose of this notebook is to test the use of a numerical relativity waveform in bilby.\n", - "\n", - "As an initial approach let's just do this with pre-prepared NR injection files (basically hardware injection files), and then let's graft support for this into Minke." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "10:50 bilby INFO : Running bilby version: 0.5.5:\n" - ] - } - ], - "source": [ - "import numpy as np\n", - "import bilby" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Bilby requires a function which takes the times as the first argument, and can then pass a dictionary of other parameters into the remainder of the arguments. \n", - "For now let's keep things simple and just pass in the path to the NR data file.\n", - "We could do additional things like including distance, but `minke` does this already, and there's no need to reinvent that wheel." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "def nr_injection(time, datafile):\n", - " \"\"\"\n", - " This function produces the amplitude for a given NR-derived signal at any given time for a given data file.\n", - " \n", - " Parameters\n", - " ----------\n", - " time : array-like\n", - " A time, or an array of times, at which the amplitudes should be returned.\n", - " datafile : str\n", - " The path to the data file containing the injection.\n", - " \"\"\"\n", - " \n", - " data = np.genfromtxt(datafile)\n", - " \n", - " hp = np.interp(time, data[:,0], data[:,1])\n", - " hx = np.interp(time, data[:,0], data[:,2])\n", - " \n", - " return {\"plus\": hp, \"cross\": hx}" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "times = np.linspace(0.25, .325, 1000)" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "strain = nr_injection(times, \"./rescale-h_m16_L0.18_l2m2_r300.dat\")" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9d3ikV3m/f593qmakGfW+VVu961333hu2sY0rAdYklADfEAjwgwQCAceBkBAw2JSAwUBiwBiw48JSbMA27gWXtb19VVbSzqprep/3/P6YGa3KdLXR7rmvS9dKbzuPyn7mmec8RUgpUSgUCsXSRFtsAxQKhUJROkrEFQqFYgmjRFyhUCiWMErEFQqFYgmjRFyhUCiWMErEFQqFYgmjRFyhUCiWMMbFWlgI0SilHEp9rgE1UsrRIp/xIaAZCEkp/2v613NutEKhUJQZi+mJ3yqEuDr1+R3AuQBCiJVCiAfTFwkhmoUQ9wkhfpwS6cm8BPw7UJ3la4VCoTiqWTRPHPgH4D4hxLXAbinlg0IIB3AtUDvpug8At0kpnxNCbBdCVALHp859Dfg88K3U19FpXysUCsVRzaJ54lLKGPBH4HLgN6ljXinl7YB30qXtQE/qcw9wp5TyPVLK95AU7BjwttT56V8rFArFUc1ixsTfCpwKnATcI4S4WUp5OMOlfcAy4DBQKaX0p09IKf9q8oXTv1YoFIqjncUMpywD/lZKGU3Fuk8gKdTTuQv4qhDCC2xfSAMVCoWi3BGqi6FCoVAsXVSeuEKhUCxhlIgrFArFEkaJuEKhUCxhlIgrFArFEmbBs1NuvfXWkndS29vb6e/vn0tz5pxyt7Hc7YPyt7Hc7QNl41xQjvbdcsstYvqxRUkxvOWWW0q6z+Vy0draOsfWzC3lbmO52wflb2O52wfKxrmg3Oy79dZbMx5X4RSFQqFYwigRVygUiiWMEnGFQqFYwigRVygUiiWMEnGFQqFYwigRVygUiiWMEnGFQqFYwigRVxxz6JEww7d8jP63X4D33h8utjkKxaxQIq445nDfeRvhvzyDaeUaPD/5LqGXnl5skxSKklEirjimiA8NEPjDw1S+7Z00/vt3MTQ04//1LxfbLIWiZJSIK44p/Nt/AUDVddsQJhO2C68g/OrzJMZGFtkyhaI0lIgrjhmklAT//CjWU8/G2NAMgP3CK0DXCT3/xOIap1CUSMkiLoRYKYR4cNqxZiHEfUKIH6fmZioUZUPswB4SI4PYzrpo4phx2Sq0mjoiO3csomUKRemU1MVQCOEArgVqp536AHCblPI5IcR24M7p93Z2drJ169aJr7dt28bNN99c0Lo+nw+Xy1WKyQtGudtY7vbB/NkY7txP5KKrER3H4Zn0/NA12/CPDhMpcM1j+Wc4l5S7jeVuX5qSRFxK6QVuF0JcMu1UO9CT+twjhKiUUvonX9DR0cHdd99dyrJl1xoyE+VuY7nbB/Nn4+BtD2BH0LRq9ZTjPqsZ9/af03jjtokwy2LYN5coG2dPuduXZq5j4n3AstTnMwRcoVgs9FCQ6L5dWLacPOOcZeMWAKJ731xosxSKWTMnQyGEEJcBDcBdwFeFEF5g+1w8W6GYCyK7doCewJpBxI3LV4OmEevphHOmv7lUKMqbWYm4lPKq1L+PTjr87llZpFDMA5E3XwGDAfPGrTPOaRYrxpZ2Yj0HFsEyhWJ2qBRDxTFBdN9OTCvXoFkrMp43rVxD7KASccXSQ4m44qhHSkl0/27Ma47Leo1p5Rrih/vRw6EFtEyhmD1KxBVHPfHD/ciAD/O6HCK+Yg1ISbyvewEtUyhmjxJxxVFPbP8uAMxrs4u4sW158tpDvQtik0IxVygRVxz1RPbvQpgtmFZ0ZL3G2NIOJL12hWIpoURccdQT278b0+p1CGP2ZCzNYsVQ30TcpTxxxdJCibjiqEYmEkQ79+QMpaQxti4j7upbAKsUirlDibhiSZLQJfGEnv+6ocPIUBDT6nV5rzU2txEfLP9eGQrFZOakYlOhWCjCsQTff6abh95w4Y/EOWV5Df94yTpW1tozXh9LZZuYlq3K+2xDYwv6+CgyGkGYLXNqt0IxXyhPXLFk8IVj/N0vXuUnL/Vy+spabj51OXsHffz9L19jxB/JeE+srwcAU/vKvM83NiabX8WHB+bKZIVi3lEirlgSxHWdTz7wOnsGfXzlbZv58tWb+ej5a/jO20/EG47x2V+/iZRy5n193WjVdWhVjrxrGJqSHesSg4fn3H6FYr5QIq5YEtz1bA+v9nv4/OUbuWhd48Tx9U1VfOqidbza7+EPe4Zm3Bfr68G0bGVBaxgbWwCIDykRVywdlIgryp79Q35+/HwPb93UzJWbZvb7vmpzC2sbKvnBs91TvHEpJfH+HowFxMMBDHUNoBnU5qZiSaFEXFHWSCn52mP7qLKa+MSFazNeY9AE7zi5nZ6xIK8d8kwc191j6H5vwZ64MBgx1DeSGFIxccXSQYm4oqz5075hXulz83fnrMZZYcp63aXrm7CbDfzfa4cmjh3JTFlZ8HrGhmYSamNTsYRQIq4oW+K6zrf/fIC1DZVcuyX3mKwKs4HLj2vmsX3DBKPx5P39PQAFh1MgGVJJjA2XbLNCsdAoEVeULY/uHuSQJ8yHzlmFQRN5r790fSPRhM7z3WNAclNTVNgw1DXmufMIWkrEM2W6KBTliBJxRVmiS8n/vHCQNfV2zu2oL+iere1OnFYjTxxIetKxvm6M7SsRIv8LQBpDbQMyEkEG1HhYxdJAibiiLHli3zDdo0Hee8ZKtAJF2KhpnNNRz7Ndo+ipzJRi4uEAxroGgHkLqSgPXzHXKBFXlB1SSn78wkGW11Rw8frCQyEApyyvwROOc6B/lMTwYEGVmpPR0iI+Ovci/tOXernsO0/zpsuT/2KFokCUiCvKjkFfhD2DPradsrygWPhkTlpWDcBfdiUzU4rZ1IRkOAXm3hN/td/NHU8cwBeJc8tvdxEroHmXQlEISsQVZcfuAS92s4HLj2sq+t5WZwXNDguv9LqBwhpfTcZQm4y/J0ZHil47F4/tG8Js0PjC5RvoHQ/xSp97Tp+vOHZRIq4oK8YCUbpGg1y1uQWbubQmmye117DDI5EGw8TEnkLRrBUIe+WceuJSSp48MMJpK2q4aF0jFqPG011z+yKhOHZRIq4oKx5+04UuJTec0FbyM05cVo1bGhlavjnnNJ9sGGrnNle8dzyEyxPmnI56rCYDJ7ZX89LB8Tl7vuLYRom4omxI6JIHdrhodVpZVZe5P3ghbGpOdizsbt1U0v2G2oY5DafsG/IBsLkladcJ7U46RwJ4QrE5W0Nx7KJEXFE2PNc9issTZmNT/raxuVhZbcGciNLpWFbS/Yb6ufXE9w/7MQgx8cJ0Qnty8/XNwypLRTF7lIgryob7XjtEnd3Mylrb7B40eIhVfhf7DNUl3W6oriMxNjJnOd37h/2sqLNhNib/u61tqASgayQwJ89XHNsoEVeUBf3uEM92jXLdlla0ItMKpxPr62aNr5/9ISMJvXghNtTUQjw2Z1Wb+4f8E8IN4LCaqLWZ6RkLzsnzFcc2SsQVZcEDOw6hCcG1W3M3uiqEeH8PK/0uwgnJIU+o6Ps1Zy0ACc/YrG3xR+IM+iKsaZga419ZZ6NnVIm4YvYoEVcsOpF4goffOMx5a+ppqrLO+nmxvm5WGaMAdA4XH7Iw1CRFXHfPPoPkkDv5IrKsemqIaFWtnZ6xgCrDV8waJeKKRedPe4dxh2LceGLpaYWTifX1sKqhEgEcGCk+JDLhibtHZ21L+p1AW3XFlOMr6mx4w3HGgypDRTE7lIgrFp37XutnRa2NU5fXzPpZ6ZFsVe3LaK+uoHO4eBE3VM+dJ+7yhAFodU59h5HevO0ZU5ubitmhRFyxqOwd9PGGy8v1W9uKahmbjcToEDIUxLRsFR31djpLyADRnMmsljnxxN0hKi1GHNapU4mOiLiKiytmhxJxxaJy32uHsBg1rto8cwByKcT70o2vVrK6oZK+8RCReKKoZwiDEc3hJDFHnvh0LxygyWHFYtQ4qERcMUuUiCsWDX8kzu93D3D5xqYZnmqpxPp6ADC1r2RNvZ2ElCVlgWjVdehz4Im7PCFanRUzjmtC0OywMuANz3oNxbGNEnHFovGbnYcJx3RuPLG4JlW5iPX1IOxVaDV1ExWSpXi7BmcNCc/sPHEpJYe9mT1xgBaHlcNKxBWzRIm4YlGQUnLfq4fY3OJgQ1PVnD033teNadkqhBC0V1cggIPjpXjiteju2eWJe8NxInGdpipLxvPKE1fMBUrEFYvCy31uesaC3DiLboWZiPX3YFq2AgCryUCzw0pvKZ54dS2JWYr4sD8CQH1lZhFvcVgZD8YIx4qL2SsUk1EirlgU7n/tEE6rkUs2FDd+LRe634c+Pjplms/ymooSPfEaZMCPjEVLtmckJeKNWUS82ZEMswz4lDeuKB0l4ooFZ8gX4fF9w1x9fCsWo2HOnhtLZaZMnuazvNZG71iw6MpIQ3UdwKy88eFA8gUgqyeeipUPeJSIK0pHibhiwXlgx6FZD37IRKy/BwDjpOHIy2tsBKIJRgPFedSaM1l4pHtKH6OW9sTr7eaM55tTsfL05mbgid/Tf93ZDH3mQ+hBVQSkKIySRFwI0SyEuE8I8WMhxIcmHT9dCPGiEOJeIcQX585MxdFCLKHzwOsuzl5dR3v1zNS72RDv6wajCWPTkSZay2qSRTXpHiaFYnAkC350X+kiPuyP4LAasZoyv9toqLJgEIIBbxjd58X9/duQ8RiRN17Gv/2XJa+rOLYobYghfAC4TUr5nBBiO3Bn6vhW4BDgAZ7NdGNnZydbt26d+Hrbtm3cfPPNBS3q8/lwuVwlmrwwlLuNi21f57CfMxsFl2+wZ7WjVBsDFQ7kVe/k8ODgxLEaGeWqVRbcI4O4tMJj4wnNhP/Cq0iEopin2VKofU1akLd1VOS89qb1Nhy6j95XugmddDaVb7mO8GsvEhwewdffj9BKe7O82L/nQih3G8vdvjSling70JP63COEqJRS+oGngAeBEeAPQohHpZRTtt47Ojq4++67S1rU5XLR2jr7VqXzSbnbuNj23fLEy4wEND5/XAdaljL7Um08/NjDmFavp37SvXXxBL/5vy5aW8ycX8QzE/YKXI9vp3rdBqrOOq8k+554zIXdbMx57d4/H2avX+eyfb/BfLCLlk/eQjjgZuSLn6T2tDOpOPWcgm0uxcbFpNxtLHf70pQaE+8D0rOv0gIOcCIQlVLqgHsWz1cchewb8vHaIQ83ntCWVcBLRUYjxAcOTdnUBLAYDTRUWYruK65VVoEQ6L7SR6gN+yNZNzXTNFRaGPaFibz+FypOOQshBJYTTweTmfCrL5S8tuLYoVRP/C7gq0IIL7BdCHEZ0AD0A/cIIfqB30kpVZ9NxQS/ejXZJ+Xq41vm/NkxVx/oOsb2FTPOtTkrcBUp4sJgRLNXoXtLi4nrUjISiNJQoIjr4RDWE08HQLNYsWw+kfBrSsQV+SlJxKWUA8C7s5z+c+nmKI5WfOHYnPdJmUw83TNlmicO0FZt5cWe4kvoNUc1iRJF3B2MkdBlQSIe0SFgrKBt/eaJ49YTT8fzo2+SGBvBUFtfkg2KYwMV7lAsCA/scBGO6dw0h31SJhPr6wYhMLZl9sSH/JGiuxlqjmp0b2nhlJFA7vTCNA2pNMPx+mVoNXUTxy0bk5v/0f27SlpfceygRFwx78QSOr94pZ9Tl9ewfg77pEwm3teNobEFzTqz2VR6qs7hIotqNIez5HBKemJPTT4RT3nqvmUbp/RTN61aC0IQ7dxT0vqKYwcl4op55w97BhnyR9h26vJ5WyPZM2VlxnNtqVawRW9uzsITHw8mi4tqKnKHjhqsyf+C400rp65dYcPYtoJY596S1lccOygRV8wrUkp++lIfq+rsnLWqdn7W0HXihw5O6ZkymXQrWFeRnrjBUV1ysc94KOmJV9tye+JOzwAAbmfTjHPmjvVElYgr8qBEXDGvvNQ7zv5hPzefumxOxq9lIjF0GBmJYJpUbj+ZWrsZoyYY9EWKeq5W5URGIujh4nubuINRDELgsObOHTAM9FMZCzBmdsw4Z+pYT2J4oOTNVcWxgRJxxbzy05d6qbWZuXzj3Ixfy0S6Z0qmzBRITtFprLIwWGTv7vSszVK88fFgDGeFMW8+fPxQL3URL6PMzGIxd2wAUCEVRU6UiCvmjc4RP891j/H2k9owG+fvT23yXM1sNFVZGSyy5auW7p9Sgic8HorlDaUAxF0HqdVDDIfiM86ZVq0FIHaws+j1FccOSsQV88ZPXuzFYtS44YT5SStME+vrQXNUTzStykRTlaXocIqhygmUJuLuYDTvpiZA3NVPvUlOdDycsr6zBs3hnGixq1BkQom4Yl7od4f4/a5Brt/aRnUBYjYbkpkpmUMpaZocVoZ8EfQi+oof8cSLz1AZD8WoKcQTHzpMfYWB0UCUhD7TNmP7SuL9B4teX3HsoERcMS/8zws9GDTBu0+bv7TCNPG+7pyhFEj27o7rkrEi+oprjqQnXsrGojsYpcaW+8VLJuIkRodpqLKiyyNpiZMxta+ciPkrFJlQIq6Ycw57Qmx/c4C3bWnJW3Y+WxLuMXSvJ2tmSpqm1Ci0YkIqWlUyY6TYcEpc1/GE49RU5PbEE2MjoCdoqKkEjszknIyxfSW6e4zELBpxKY5ulIgr5pz/fbEXAfz1aTNL4OeaWG8XAKYVHTmvS0+cL2a6vDAY0SodRYdTPKlNyuo8nnhiKJkj3tCYLLcfyfAuIV3ApEIqimwoEVfMKYc9IR5+w8XVm1smBgHPJ7HeVGbK8tU5r2uqSnvixZfeJ7zFNc+aqNbMExOPDydFvLklOSw6mycOqJCKIitKxBVzyp3PdCMQvO/MlQuyXryvC2GzY6hryHmdw2qkwmRg0Ft8wY/u8xZ1jztVrZkvO2XCE1/WioCMGSrGphYwmia6NCoU01EirpgzDgz7+e3OAW46sW1BvHBIeuKm5avzVoMKIVJphkV64lXOosMpaU88XzglPnwYzeHEbLdTYzMx4p8ZThEGI6a25coTV2RFibhizvju013YzAbec8bKBVsz1tudN70wTSm54prDWfR0n7QnXp1vY3NoAENDspK1vtLCcCCzbcb2lcoTV2RFibhiTtjR7+bJAyP8zekr5j0vPE3C60Z3j+aNh6dpcpRQtVlVvIh7UiLuzNM3JTE8gDEl4g2VloyeOICxbTnxwUPIxMyqToVCibhi1kgp+eafO6mzm3nHScvy3zBHxFObmqblhXni9ZWWrEU12dAcTmQoiIwVPmnQG45jNxswGrL/95JSEh8awNCYHFVXbzdnjIkDGFvaIZGYiKErFJNRIq6YNb/bNcDrLg8fPnc1FWbDgq0b60ulFxboiTfYzVmLarJhqEo3wSrcG/eGY3lH0MmAHxkKTAmnjAWjxHV9xrXG1uQLY/xwf8E2KI4dlIgrZoU/Euebf+5kc4uDqzbP/QDkXMR6uxHWCgz1M3txZyI9eT5TPnY2Jgp+ihBxTziGM09IKZ1emA6n1E+8wMz0+I3Nyd4zSsQVmVAirpgVdz3bzVggyqcuXpe37epcE+vtwrhsFUIr7M84XT2aKR87G1qqCVYxFZPeUDxvH/HEyCAAhoamKbZlbIRVW4+wWIgd7ivYhkzEEjqv9I2zf8g/q+coyouSpt0rFADdowHufaWfa45vYVPLzKEG8028rxvLCacXfH19ZTJbJFvsORPp/inFpBl6wjEaqypzXqO7kwVEBmdNyrb0C0yUjdOuFZqGsbmduKt0EfeFY3z0vh3sPJzMeb/pxDb+8eJ18zaoQ7FwKBFXlIQuJV9+ZA82k4EPn5u75H1e1vf7SIwOF7ypCVBnM6eKaooJp5QWE3fmiYknPGPJ51cnR9ble4ExtrTPyhP/yh/3sXfQx2cvW8+B4QC/fLWfZTU23nnywm1EK+YHJeKKkvjlK/28dsjDLVdspDbPRPf5YKJnSoGbmgBGg0aNzVRcOCXtiRco4lLKZDilIvd/Ld3jRlgsCGtyiPPEC0y2XPGWZYRfeR6p6wWHj9LsHvDyyO5B3nfmSq7b2oaUEpc3xHee7OSyDU3ULcLvTzF3qJi4omj63SG+81QnZ62q462b5m/sWi4mRLzAQp80dfbsRTWZEBYrGE0Fh1MC0QQJKfNmpyQ8Y2iOmolwRvIFxsxwtlzx1nZkNEJibLhg29P87C99VFqMvPvUZFtgIQSfuHAtsYTOT15UjbWWOkrEFUWhS8m/P7IbgxD882XrFy2mGus5kMxMaWot6r6GSgujRYRThBBFVW16w8nsknwbm7p7fCKUkqa+Mk+uOMVnqPgjcR7fN8wVxzVRaTli0/IaG1cc18x9rx3CXUTKpaL8UCKuKIp7X+7jL71uPnbBmgXrj5KJWM9+TCvXFB1aaKg0FxVOgeSYtkJFfKJaM1/zK8/4xKbmEdssWdMfjS2pXPEiNzef6hwhmtAzDqq++dTlROI6v901iB4J4//Dw3juvYvo/l1FraFYXFRMXFEwuwe8fOvPnZy/pp5rtxTnAc8lUkqiXfuxnXdp0ffWpYpqErrEoBX2LkJzFN4EyxtOlsbn29jUPWMzeqDXV5rZM+jLeL2hoQmMxqI98cf2DdNQaWZz68zsoTUNlWxucfDQa31c+JPPEj+wO/k9/OR7ON71AZzbPlTUWorFQXniioIIRON87tc7qbWb+ZfLNy5qalpiZBAZ8GFaubboe0up2tSqnAXniXsmwinZRVxKie5xY6ie6onX2y2MBTJXbQqDEWNja1EiHokneK57lPPXNGTN4b/m+Ba6xsPsHI1Q99mv0HbvY9guvgrvPT/A/8iDBa+lWDyUiCvyIqXkPx7dyyFPiC++9bgFa3CVjVj3AQDMq0oQ8ZIKfqoLj4mH8sfEZSiIjEbQMoRTJJmrNiG5uRkvIs3wDZeXSFznzFV1Wa8517MXayLCU+e8G9vZF6NVOaj92OexnnQG4//9lYkNZEX5okRckZefvtTLI7sH+X/nrOakZTX5b5hnYt37ADCtXFP0vZOLagpFq3Kg+zxImb9xVjqckssT1z3pQp+ZG5tJ27KnGcZd/QXZAfDSwTEMQnDSsuqM56WuE7v3e5wR7OHpqJNYIvkOQBgM1H7qiwiLlfHvfa3g9RSLgxJxRU6e6x7l2092csn6Rt5z+vzPzCyEaM9+DE2taPbcVZGZmCiqKSLNUHM4IR5HhoJ5r/WEYlSYDJiN2f9rJVIirmUIp0D2YiRjSzsyFCh4cPNLveMc11w1JStlMpHX/0L8YBeXnbASbyTOiwfHJs4ZnDU4//rDRHa8SOi5JwpaT7E4KBFXZGXvoI9/fvhNOuor+cIix8EnE+s+UFIoBSZXbRbfP6WQkIonHCsgvTAplpmyU8hhWzEZKnFdZ/eAjxOzeOEA/kceQKtycv4VF1BpMfKHPUNTzldecT3G9pV4fvpdZIY4vaI8UCKuyIjLE+Jj9+/AbjHy9eu3LGiL2VzokTDxQwdL2tSEyVWbxbSjLbx/ijccLyi9EECbFk6ptZsQ5AintKZzxfOL+GggSlyXHN/qzGKDm9CzT2C76ErMFRVcsLaeJ/YPE40fEWthMOB81weIH+wi9NQf8q6pWByUiCtm4A5G+eivdhCN63zzxq2Lmg8+nXhvN+g6plXFx8PT1NktRTbBSvVP8Rci4vl7iesTIj7VSzZqyarNrLniTa0gREEZKukXgmyNyYKP/QbiMSrfci0Al25oIhBN8FzP6JTrKs69FOPy1Xju/aHyxssUJeKKKbiDUT78y9cY8Ib5+vVb6KgvPu48n0R79gNgWrWu5GfU2c2MFpVimOopXpAnHss/ls09hqiwoVlmvjg25KjaFCYzhobmgsIpQ74IjVWWiRDNdEIv/BnTqnUTueqnLa/BWWGaEVIRmobzHe8n3ttF6NnH8q6rWHhUsY9igrFAlL//5av0uUPcdt3xnNCePZ66WMS69iEsVozNbSU/o85upns0UPD1xfQU94QK88Snpxemaai05Ex/NLYtJ+7qzWvHsD/KpubMXrgeDBDZ/TpV12478lyDxkVrG/j97kHCsQRW05HwWcU5l2C85/t47/0hFWddVHSVrC8c4/4dh3hy/whD/gg2s5HNLQ6uXGmmpSWZ+RLr2kestwsZi2Korce8dhMGZ/n9/ZUjSsQVQNJz+8ivXsPlCfH167dw2ora/DctAtH9uzCt2YAwlB6jr7WZGQtGkVIWtFmrFRgTl1LiCcdxFFJyX53551tXaWZ3lqpNAFPbcgKP/y6n7e5gFG84xuZlmSceRd58BeJxrCdO7cV+6YZGHnjdxTNdo1y8vnHiuDAYcPzV+xi77RZCLzyJ7cwLcn5/k3mhZ4wv/GYXY8Eox7c4OGV5Db5InMf2DSEDRr71WBfv23M/q7v+MuNey9bTqLr67VhPO3dWv++jHSXiCvYN+fjE/a/jj8S548atnFwGueCZkIk4sa692C+/flbPqbWbiSUk/kicqjxeM4AwGhE2e97slGAsQUKXecMpunt8YqLPdBomVW0aM3i8xrYVyIAfPccLwc6B5OCHTRlK7QHCrzyPsFiwHLd1yvGTltVQazPz6J7BKSIOYDv/LXjv+QHen99FxRnnF/Ti9/AbLr786F5W1Nq4/YYtbJz0ziAwNMSTf3iCZ8dN/NOyG7nx5Gv5f+etpcJeQXx4gMiOlwg8+jAjX/oUxtblON75t9jOf4sS8wwoET/GeaZrhM8+vJNKq5EfvOsk1jVWLbZJWYn19SAjEcxrj5vVc2ptyVzx0WC0IBGHwqo2j1Rr5gmneMcxr50+vydJumpzLBCjsWpmPNvYmmwnGz/Um1XEdw8kPfkNTZl/l+FXX8Cy+WSEeerzDZrgkvWNPPSGi0A0jt18RB6EwUjVX72P8dv/jfBLz1Bx2jk5v8cn9g/zpd/v4bQVNfzn246fkqse7dqH59ZPsOrEs7j71K38kHZ++bqL5x8b4vOXb+SEzSdh3XwSjne8n9CzT+D9xQ8Zu+0LeH/xQxzv/AC2cy8tWMwj8ei0nikAACAASURBVAQHx4IMeCMEonHiusSkCaorTNTYzLTZNCyeYXT3ePL3KwRoGsJkJmG2kbCY0Jw1CGP5SmVJlgkhmoFvAz7geSnlnbmOK8oPXUp+8mIv//1UJ2sbKvnGDVuzboKVC+nuetkEsFDSQxDGAlFW1toLuqeQdrQTza9yhFOklCTcYzPa0E7YNqkYKZOIm9qSIh47dBDLphMyPuPAsJ8mq2mKCKeJDw0Q7++h8vLrMt576YZGfvlqP08eGOGK46Z2PrRfeCXen9+F9+c/wHrq2Vm98b2DPj7/m50c1+Lga9dtmRJfj3bvZ+gzH0KzVmC/5GqaNm/hs6l1v/TIHj7481d4x8nL+PC5q7GajNjOvYSKsy8i9OxjeO/5AWNf/Re89/4Q5zv/lopzL80Ynx/0hXmmc5Snu0Z5qXeMcCx3Vk1txE17cJjW1EdLaJSW0AiVZ5wJj/8ahEBz1mCorU99NGCorUerqUerrSfqqMMnTPjj4I2BN5rAG4jgDUXxhGJEY3Fi8TixWJz3XXYCK1bMbfO4Ul9ePgDcJqV8TgixHbgzz/EJOjs72br1yNu4bdu2cfPNNxe0qM/nw+VylWjywlDuNvp8PjoP9vHE/mH63SE+cnIt562pJ+YdxeVdbOuSZPsZhoIRYpdez5AwIWbxM7bFIly1ykJwfASXIVTQPYFTz0NGo8Rcrqz2ed0hrlploSruxeXK3P9ERiL4z7uC+LI1BDM8o1pP2uYZGcSlzxxoLHXwX3wNsbjEm+Vn0GwIsb7GkNHGaOceQhdeBeu34stwvkFI/mq9nYHDh3FVzxS/6Hs+RujFpwi/8hKmVI/zycQTOr9/3cU1q61ct6WRseHBiXO6z4v/Dw8jzr8S+6XXEJRM2Nhmgm9esZwXD46xa2CI7z4yzmnLa1hZZ0u+WKw+Dvm5r6P3deN942XcLz2Htm8PlvWb0VqWMZIw0ucO0jsWYiyVedRmMXLGpgoaEn4qvKMYx4YgHEAXGlGTlYizHn9lHT5jIx69GU8MhnTJELADWFdtpHfVxVjQMckEQk8gdR1d6sSlRjRgIhoyoruy72FomsBgNmA0CcxWIyOHXZjmuPVQqSLeDvSkPvcIISqllP4cxyfo6Ojg7rvvLmlRl8tFa+vitUAthHK38cWdB/jX3x7EE4rzyYvWct3W1rKpxEyT7Wc4+MR2hMVKY1vpmSkAZn+E7Q92s3F1BecU+LsaPdRNdN9OWt774az27fQOsr07wrbzm2htyJyaGTt0EP3x7dSefDr2DM8w+SNsf6ibzWtsnJnFNm3v65iCPuqvvHbGuVA0wT279/CF85oy2jhy97cw7XyV1k/ekvX3btof4s6X+7j69OMmes2kkQ31HP7R1zEM9NL41R/OeMZ//XEvv9rl59s3ncCalUfebSRGhxn8l09TFfLT+JUfYFq+OuPPccWydl7uHecrf9zLPbv7WFFr48rjmjl9ZS2r6+1YW1sJn3AaBx5/nJcff5Yde3ewozaIz2RHkzrHxUc5NdLPSYNv0HpoN0JPAGBobMGycQvmDcdj2bgV06q1M0IkUkrGgzH6xoP0ukNI/yi6bmQwHMcbjqFL0DTQhMBm0nBoOg49QlUiTKWI4zDoVGmSKoPEYTPjrKygosKOVmFDWCvQrBVojuo5D82U+rQ+YBlwGJgs1NmOKxYZdzDK1x/fjyHkptJi5PYbtpZ1/Hs6MhYj2r2fqqv/atbPqrYlKyPHim1Hmyc7pZBwyvQp99OpSdmWqxgpV5ph54gfyZG4/2RkIkHktReT2R45Xriv3dLKT1/q5aE3XLz/zKnj74TJTNWN78H93a8knzUpw+WZrlF+9eoh3nXyMk6fJOC6z8vw5z+C7h2n4cvfzTsX9eTlNfz8PafzyO5BHnz9EN99uovvPj29m6IJGs6n3mrgDJOPU4J72Dq8i8p4CGEwYOhox3j26ZjXbMS8bhOGuoaca0JyilOt3Uyt3czW9mpcLsk1ZeyQpSlVxO8CviqE8ALbhRCXAQ3Tj8+RjYpZkNAlv37zMN95shNfJM5nz27kn0/ZmLNBUzkS6+2EWBTzutltakKyMrLaZmIsS2VkJjSHExnwIROJrNekR7NVZWk4BclCHyBrnrhR06i1m3POATW1rcg6NHn/cNJvyjS8Ota1F93nwXrSGVmfDbCi1sZpK2p4YIeL95y+csbwjMrLrsF3/92Mf++rNH/rZwizBXcoxhd/v5uOejsfPu+ISOvhEMO3fpzYoV4abr0Dy7pNOddOY9AEV25q5spNzQx4w+wZ9NE9GiCekBgNglanlXWNVaystZXdO8mFpiQRl1IOAO/OcjrbccUCIqXkuZ4xvvnEATpHAmxpdfLPl63HFvMuOQEHiO5PTp0xr5ndpmaadK54oUzkivuzxz89oRgWozZlI286ujfliWfZ2ITUmLYcvV2MbcshFiUxPJAsxZ/E/mE/NpMh4wtJ+JXnAbCecFrWZ6e54YQ2Pv3QmzzdNcL5a6Z6scJsoeYjn2XkCx/F/eNvUfOhT/G1P+3DHYpxxw1bsRiT378eCTP65U8T3fsmdZ/5j4LWzUSzw0qzw8oFa/N708ci5Zs3oygJKSXPdo/xP8/38NohD+3VFfznNZu5aF0DQghc5bJ7WSTRfbvQKh0YZlGpOZlam7k4T3yik6EbtJleLqSaX+Wbcp8Kp6T7sWSi3p57DqhpeTLEETvYOUPEDwz7WdNQmdE7Db/6AqbV63K+gKQ5b009DZVm7n25f4aIA1ScfCaVb3sn/od+zrOOtTzSX8WHzl7F+lRaY8LjZvTL/0Rk56vUfPRz2M6+OO+aitJQIn6UEI4l+NPeIe55uY99Q36aqiz848XruG5rKybD0vO8pxPd8zrmdZvm7K1zrd3MG67CpvXAtKrN6sweoS8coypfoY9nHGGvQuRIUaivtLBrIEfV5opk869YzwEqTjt34riUkv3DAS7b0DjjHj0UJLJ7B1Vve1dO+9IYNY2bT13ONx4/wAs9Y1Ni3Gmq3/cxBgZG+HqXYK01wM2rreihIKGn/4T77u+gez3UfuqL2C+4vKA1FaWhRHyJ0zUS4IHXD/HbnQN4w3FW1tr4whUbuXxj01Eh3pDcGIv1dlFRwmDkbNQVG05xTOopnkXEPeF43kKfhGdsxmzN6TRUmhkPZq/a1OyVGBqaifUcmHJ8yB/BH4mzJkNmTOSNVKl9nnj4ZG48oZ17X+7nm38+wP8uOwXjtL+nhNC4reMGoq5xPvr81xl+bHjinGn1Ohr+9Q7MHesLXk9RGktGxH0P/ZxgXOIzaZhXr8e0eh2arbBCjaONnrEAf9o7xJ/2DrN/2I9RE1y4roHrt7Zx8rLqo26jJ7LndZASy8bMxS2lUGs3E47pBKNxbBmKYqZTSE9xXzhGW3VFzufo7uzNr9LU23NXbQKYVq2dIeKdI8mmXh31dmBq/nv41ecR5pml9rkwGzU+fsEaPv3wm3z7yU4+fuGRHu4JXfKF3+5ih8vLv1+1hVPedSfhvzyDjIQxr9uEZeupR93fYbmyZEQ84R4j7gvg/t0vkweEwNi6DFPHeswdGzB3bMDUsR5DjljjUiUUTfBK/zjP94zxQs/4RAe+La1OPnHhWi7f2JQxG+FoIbJrBxgMmNcXltlQCLWTqjYLEfFCeop7w3E2FtDBMF06n436ScOcs4r4yjWEX34WGYtNhGY6h5N/F6vrKwmOTxfxF7BsPmlGqX0+LlrfyNtPbOdnf+nDYjLw3tNX4A7F+PKje3iue4x/OL+DyzYm+8Ckq0kVC8uSEfHqv/l7gi4XTe96H9EDe4h27SF2YC/RPW8SevLI1BFDYwvmDZuxbEgm9ptXr88Zfyw3pJQM+SPsdHnZOeDlDZeXN1we4rrEYtQ4oc3JdVtbuWhdA01V5TOsYT6J7tqBuWM9mjW3l1sMR/qnxGgvoN+XsNlBM+TMFfdG8sfEE55xzHm84YZ06X2uzc0VHZBIEOvvmRhV1znip95uprrCRHD8yLXx4QHifd1UvuVtOdfNxscvXEMgGudHz/XwP8/3oEuwGDU+c+l6bjhhbjaaFaWzZEQ8jaG2norTzpnSgCfh8xDr3Eu0cw/RfbuI7n59QtiF2YJpzYaUqG/GsnErhtr6xTJ/gnTb0p7RAF2jAXpGg3SPBjgw7J+Y7GLUBOsaK3lnqnhia5szZ/ra0YiMxYju24n9yhvm9LmT+6cUghBiYup9JqJxnXBMz5mdInUd3evOWuiTJu2JZ5vwA0lPHJKbm2kR7xoJpEIpUwm/+gIAlhMLj4dPWcugccsVG7nm+BZeODiOzWTg4vWNtOcJHSkWhiUj4ofcIcYCUWJjQSxGbeLDbNQwVjkxnHDalDzU+MgQ0T2vE9nzBtHdb+B7+F74v2QxhrF1OZbjT8ay5WSsW06ZU1HXpcQdjOIOxaZ8jAWiDHjDDPjCDHgjDHrDBGNHCkesJo2VtXZOW1HLcS1VbGp2sK6xaknmdM8l0c49yGikqFhuIaQ98WJzxbPFxCcKfXJ44rrPC7qeNyZeYzOhieyzNgFM7SvBaCTWvR8uvIKELukaDXDD1pmecfjV59Fq6yem+JSCEIKTltVwUpm2KT6WWTIi/v1nutFC42zvnl5+m6zushg0LCYNs0HDYjSkRL4ac+X5mM+4ENNZIIIB8LnBMwauMcTBNzA8vAOjzYaxrh5zXQOm+iaE2UJClySkRJcy+bme/FyXEEvohGKJ5Ec09W9MJxxLcFG7ke3dezN+DzU2E81VVlbU2jhjZS3NDgsrau2sqrPR7LCiqY2gGUR2vgaAZePciniNLekxF5uhks0Tnyi5z+GJ657UlPs8edpGTaPWZmYkR9WmMBoxLVtNrHsfkBxsHYnrdDRM9cRlIkHk1RexnnaO2mg8SlkyIv7OU5YxOmjitOOqiMR1onGdSFwnEk8QTehEYjqRROrriXPJD18kTjSuo0sjurWOhLkOvU4nEYsRj0RIxOIkQnH0/hB6/0E0gwGD0YDBZMZgMmHQBJpIvlgYNIFR07CaNCpMBmptZqwmAxWpjzW2KBtWO6iuMFFdkYxPOitMVFeYjrlQyFwQef0ljO0r5jwEZjJoOKzGogt+EkMDGc/5CvDE85XcT6Y+T9UmgHn9JoJP/RGp65MyU6amF0b370b3eag4+ay8ayqWJktGxDc0VeFK+Ghtbc5/cQnIRJzo/t2EX3mO8F+eJbpvJ0iJ5qjGetKZWE85C+spZ02kmmWj3LsYLiVkLEbkzVexX3zVvDy/2NJ7g8NJrDPzuyxPyhPPlSeue9xAoSJuZsiX3RMHMG84nsDvHyB+qJfOkeSsylV1tinXhF9+BjQNy7RRbIqjhyUj4vONMBixbDgey4bjcb7rgyQ87glBD7/yHMEnfgeaAcuWk7GddREVZ15QFhukRzPRvW8iwyEsW0+dl+fX2ossva90ovvcZNqlSHvijlyeeDqcUoiI2y3sOpy9ahPAsv54AKJ73qAzvIpWp3VGumT45eeSXfyOwtRbRRIl4lkwOKuxX3gF9guvQCYSRPfvJvTCnwk98yfG//s/Gf/uVzBvOJ6Ksy7CduYFGDM0yFfMjvCOF0EIrFtOmZfn19rM7BvKLZST0RxOZCSCjMdnnCvMEx9PTolx5H43B5OqNhP6jErJNMb2FQh7FZHdO+isapxRqZkYGyG6byeOd30w73qKpYsS8QIQBgOWDZuxbNiM868/TLyvm+CzjxF69nE8P7wdzw9vx9SxHtv5l6NvLS2NSzGT8I6XMK3ZiFaVeeDvbKmzmxktoQmWjIZnnEt74pU529COo1U5EYb8/+3qU7M2RwJRmh2Z6wGEpmHZfCL+11/h4MazOX/N1HeGwWcfAymxnX1R3vUUSxcl4kUihMC0fDXO5atxvuNviQ8cIvjs44SeehTPj+7Af2EXQ2MD2C64HNvZF6PZM094UeRGD/qJ7nmDqusKG91XCrU2M4Fogkg8MdE+NRdpD1pGZsaqveE4VRbjjN7bk9E9YwXFwwFaUsJ92BvOKuIA1q2nsvfNvSR0OSNHPPTUHzEuXz2r1EJF+XNsJyHPAcbmNhzX30zTN+6m+fv3Yzn+JBJDhxm/44sc2nYZI1/+J0IvPZ1zmIBiJuGXn4dEgopTc09Vnw21RRb8THjiGUW8kA6G7rzNr9IsSxXS9I0Hc15nPeVs+uzJzf7JmSl6MEBk56vYzlEtYI92lCc+h5jaVmAVJpovu4ro3p0En/gdwScfJfTMYxgamrBf9jbsl16DsWF+MmyOJkIvPolW5cS88fh5W+NIwU+MFmf+6sN0/5RM4ZTCeomPTVRa5qPZacWgCfrGcw9yNrUtp79tEwaps6L2SGZKdP8uTIDtwisLWk+xdFEiPg8IISZi6NXv/zih5/+M/5EH8P7s+3h/fhfWk87Efvl1VJx2TkHx0WMNqeuE//IM1lPOmtefT7r0vtC4uFaZjM3L6ExPvNBe4oVkpkCy4KfNWZHXEwfob15H68gQ9HfDig70cIjo/l04zrwAU+uygtZTLF1UOGWeESYTtnMvofFL36Hlhw/huOk9RLv2MvqlT3H4/dfivf8nyXJsxQSJkSF0r4eK08+f13Vqi6zaNEzExGd64vl6ictEHN3nQSswnAKwrKaC3jyeOECn5mRVaAjPT74HgO++u5HRyLzuJyjKByXiC4ixuQ3nX3+Y1v/ZTt3nvoqxuRXPj+7A9TdXMv7fXyHW37PYJpYF8UMHwWjEevL8ZvrU2IpsgmW2ICxW9AwxcV8kljNHPF3oY3DmH42WZllNBf3uEFLKrNe4g1EGAzE2rm4l9NzjDH3mg3h//gNMq9bOeb8ZRXmi3ssvAsJgxHbWhdjOupBo5158D9+L/5EH8f/mV1hPOYuqa7dhOeG0Y7LXhZSSWH8PlZtPQrPNb2aP1WTAbjYU3T9lejhFSok3lNsTT3hSszWdhRfdLKu2EYolGA1EJzobTmfvUHK6/ZbLLqLKPEjw6T9if8u1aKedV/A6iqWNEvFFxtyxnrpP3EL1ez6C/3f34//NfQz/y99jXrcJx9vfi/X08xAZRnQdrcQ696L7PNjOnbtRbLkoZer99HBKMJYgIWUeTzwl4gUMKU6zvCa52do7Hswu4oPJYqX1zU6c7/0o1e/9KJBs/6A4Njh21KHMMdTU4XzXB2n98a+p+chnSXjdjHzpUwx85B0EHv8dMjGzSvBoJPjE70HTqFigApWiS+8ziLg3lC65z+GJu9Ml94WLeHtNMtskV4bKniEfrU4rzoqlM/hEMbcoES8zhNlC5RXX0/L9+6n91BcBGPva5zn8wRvwP/JgxpLvowWZSBB88hGMLcvyNhqbK0ryxKeFU9JtaHP2Evemm18VHk5pdlgwaoLeHBkqewf9rGusKviZiqMPJeJlijAYsV94Bc3fvpf6f/kaWpWT8W9+iYG/uynlmR99xUORN18hMTo8MalmISjaE3fM9MR9kQJ6ibvHQDNMpCkWglHTaK/OnqESiMbpHQ+yvlFVBR/LKBEvc4SmUXHmBTR943+p/8LXEZYKxr72eQY+8k6Cz/wJqeuLbeKcEXjkQYS9Mu8g4bmkzmbGE44TTxT2czSkPPHJP/d0OCV3L/FxNGd10fsbq+vtdI34M57bk4qHb2hSnvixjBLxJYIQgorTz6Ppmz+l7jP/CbrO6Jc/zeDH303oxadzpqEtBRKecYLPPIb94qsQxoXbb0+X3o+nhDgfWpUTpEQGAxPHvJECOhh6Cy/0mcyahkr6xkOEojPfee3oT04ZOr51YUJPivJEifgSQ2gatnMvofm/f0HtJ29FD/gZufXjDH/mQ0T27Vxs80om8MdfQzxG5RXXL+i6E2PaCq3aTBX8TB7TdmRjM/dUH62ITc00axsqkSQn2U/n1X43HfV2tal5jKNEfIkiDAbsF72Vljvvp+bDnyHW183QJ/6G0f/6HPHBpZVeJhNx/L+9H8umEzEtX72ga9elCn5GC9zcTDfBSqQ2KiHpiRs1QUWO8XvJkvviBzOkQyVvHp5a1ZvQJW+4PJzQroY9HOsoEV/iCKORyrfeSMtdD+D4q/cRev4JDn/wBtw/vGPJlPMHn/ojiYFDVL7tnQu+dqmdDKd74g6rMWdxVsIzXlSOeJpmh5WmKgs7Dk0d0Hxg2E8gmuDEdhVKOdZRIn6UoNkqcf71h2n+/v9hu+ByfA/8lMMfuA7fg/cgY4XFexcDqet4f/EjjCtWU3HmBQu+ftEiPhFOOfIC6Yvk6ZsSiyID/pJi4gBb25y8fsgzZd/jpd5k8ZDyxBVKxI8yjPVN1H3iFpru+CmmjvW4f/B1Bv7uJoJP/7EsNz9Dzz1BvLcLx9vftyiVqTaTAYtRKzhXPJMn7gnl7mB4pOS+VBGvZsgf4bD3SGrjY/uGWNdYSVNV9oERimMDJeJHKeaO9TR86TvU3/pNMJsZ/Y/PMPSP7yc+MrjYpk0g43E8P/sextblC1ZmPx0hBHX2wgt+NHsVINC9R0TcF8ndS7yUkvvJnLEqed/j+4YBGPCGecPl5ZL1jSU9T3F0oUT8KEYIQcUpZ9H8rXuo+Yd/IT5wiMCjDzHyn/9MfODQYpuH/ze/In6wi+r3/QPCkH882nxRYyu84EcYDAizGd03aWMzTy/xhDsp4qWGU5bX2NjYXMUju5MvwA/sSP7uLlYirkCJ+DGBMBipfMu1tPzgASybTyb84pMc/tCNuO+6fdE2PxNjI3h++j2sJ5+F9Yz57Ruej2TpfeH7BsJineKJe/P0EtdnGU4BuHxjM7sHfXzziQPc83Ifl6xvZHmNLf+NiqMeJeLHEFqFDeuWk2n+/gPYL7wc34M/S25+PvTzBd38lFIydvu/QTxO9Yc+uegtd2vtpqKm3guzBd2ffPFL6BJ/JJ6ng2G6+VXpIn7DCa1sbnHwk5d6MRs0PnZBYWPeFEc/SsSPQYz1jdR+/Baavvmz5Obn929Lbn4+86cF2fz0//oXhF9+Fuf7P4apbcW8r5ePOpsZdyhKQi/sexcWK4mUJ+4roFozMT4GRhPCXnqPE4vRwDdu2Mq/X72JX7z3dJodakNTkUSJ+DGMefW61ObnHWAyM/rlTzP0T39LZNeOeVsz9PJzuO/6BtZTz6HyrTfN2zrFUGs3o8tkbLsQhMUykZ2SvidntaZnHENN7azfcVRXmLhsQ1PW3uKKYxMl4sc4yc3Ps2n+9j3UfOSzxF19DP3j+xn63IcJv/nKnK4V2fUao//xaUzLO6j7py8tehglTXrqfaEhFWE+EhP3hQvom+IeLankXqEoBCXiCiC1+XnF9bTc9SDO932MWM8Bhj/9QYY+/UFCLzw569a3weeeYPhzf4+hpp6GW++Y99FrxTBR8FNgmqGwWJChADIWwxNOdzDMHU4xlJheqFDko+h2cUKIZuDbgA94Xkp556RzpwPfArqA/VLKz8+VoYqFQauw4bjh3VS+9SYCjzyA7/6fMPJv/x+GplYqr7wB+0VvxVBbX/Dz9IAf913fIPDoQ5jXHUf9v94xqw2++aC26IHJyXi07vPgCSU98eocTah0zzjm1QvXI11xbFFKz88PALdJKZ8TQmwH7px0bitwCPAAz86BfYpFQrNaqXrbO6l8602EnnsC/29+iefH38LzP9/GctwJWM84H8txJ2DuWI8wTRUwKSXR/bsIPf1H/L9/EBkKUHXTe3C+6wMIc/nFcydEvNCCH0vye9D9Xtyh5PeerZOglLLkDoYKRSGUIuLtQE/qc48QolJKme6T+RTwIDAC/EEI8aiUcsr78M7OTrZu3Trx9bZt27j55psLWtjn85X9ANhyt7Ek+zqOg3/4V/CME+vtItjbhd51ALoOJKfV2CvRrDYQIGMxdJ8HGY+B0DBd/x4sm08kWFtPcGR0/mycBVJKrlltxRbzFrRu0FwBF17F4bFx7DErV6+y4B8bIpghxi+jEfznvoX4inUEF/B7Kve/Qyh/G8vdvjR5RVwI8XHgkkmHXgCWAYeByQIOcCLwWymlLoRwk4y5TxHxjo4O7r777pKMdblctLa2lnTvQlHuNs7KvtZW2LgJgMToMJE9rxPdu5P40GF09yhIEFYrxpZlmNdsxXr6uSXNylyMn+HzQz1Im4Frz8y/bmJsBB7fTt3ZF7DLX8lTAzpfaGvLeG2svwf98e3Unnw69gX8nsr97xDK38Zyty9NXhGXUt4O3J7+OhUT/6oQwgtsTx27DGgA+oF7hBD9wO+klOXbPk8xKwx1DdjOvhjb2RcvtilzQm0R/VOExYIkOfzYE6nPOZRBT5fcq41NxTxRdDhFSjkAvHvasUcnffnn2RqlUCw0tcX0TzFbkyLu8+BOxHJuaibcyRBSqc2vFIp8qBRDhYLkmLZCp/tgNILRhO7z4gnlE/FUyb0SccU8oURcoQDq7GbGg9GC2g4IIdAczlSKYSx/OEUINIca3qCYH5SIKxQkwymxhJzohZIPQ5WThNeDO68nPormqF7UVruKoxsl4goFxVdtag4nIZ+faELPPRDCrao1FfOLEnGFguKrNrXqOsb9IQCqbTk88dFhDHUNszdQociCEnGFgmRMHCh4OIShth63LznzMldMPDE2glarRFwxfygRVygo3hM31NTh1ZMVmtVZwilS10mMjRTVa0ahKBYl4goFSW9aExScZmioa8BrsifvzRJO0d1joCdUOEUxrygRVygAgyaorjAxXrAnXo/PlJxxWV1hznhNYmwkeW2dGmismD+UiCsUKZIDkwv3xH0mOwJJlSVz4XNibDh5rQqnKOYRJeIKRYpae+Gl94aaetzmKpwGHYOWeUJRYjQl4iqcophHlIgrFCnq7RZGCu2fUlmFx+KgmuzZLImxYRACQ3XdXJmoUMxAibhCkaKhysKwP4JeYOm9p6Ka6kQo6zWJ0WE0Zy3CWErbfoWiMJSIKxQpGirNyzVAvQAACSVJREFUxHWJJ1RYrrjHXIkz4s16Plnoo+LhivlFibhCkaLBnhy7NuSLFHS922DDEXJnPZ8YG8agCn0U84wScYUiRUNVUsRHAvlFPBxLEBJGnL6RrNckRobUpqZi3lEirlCkaKhMiviwP//m5ngqFdERGEMP+mec10NBdK8bY1P5j/dSLG2UiCsUKepT/VOGCwinpHusVMf8xIcGZpyPDyYH7CoRV8w3SsQVihRGg0atzcSwvxART3rizqifRAYRTwwdBsCgRFwxzygRVygm0VBpKUzEA2kR9xFPCfZklCeuWCiUiCsUkyhYxNOeuIxMeN2TiQ+6EBaLGpCsmHeUiCsUkyhGxO1mA/aGBuIDh2acTwy6MDS2IkTmknyFYq5QIq5QTKKh0sJYMEY8oee8bjwYo9Zmxti2grird8b5+KALY1PLfJmpUEygRFyhmER9ZTJDJV8PldFAhBqbGWPbcuKuXqR+RPSllMQP9WJsWzGvtioUoERcoZhCY1U6Vzx3SGXIF6WxyoKpbQUyEiExMjRxLjE6hAyHMCkRVywASsQViknUF1B6L6VkyB+mscqCsW05APFDPRPn4/0HATC2KxFXzD9KxBWKSbQ4rQAMeMNZr4kmdMIxncZKy4S3HTt0JC4e6+0CwNS+cv4MVShSKBFXKCbhsJqotBg5nEPEA5EEkAy9aLX1CHsVsa79E+ejB3ajVdehqYk+igVAibhCMY0WhxWXJ4eIR+MANFZZEUJgWb+Z6L43J85H9+7EvO44lV6oWBCUiCsU02h1WnF5sg97CERTnngqk8W8YTOxg53owQB6wE+8vwfz+k0LYqtCoUaOKBTTaHFaeengOFLKjN50IBpHAPWprofmDVtA14nu3wV6UuAt6zYvpMmKYxgl4grFNFocFQRjCTzhONUVphnnA9E4tXYzJkPyjaxl/WYwGAi98CQyGEBYKzBv3LLQZiuOUZSIKxTTaE1lqBz2hDKLeCQxkU8OoFVWYTv3UvwP/RwA+6XXoFXYFsZYxTGPiokrFNNocaRFPPPmZiCaoLHSMuVY1XU3gxBgMlN59dvn3UaFIo3yxBWKaaQ9cVeGNEMpJYFInMYqx5Tj5jUbaPnRw2iVDjSbfUHsVChAibhCMYMqqwmH1UjfeHDGOU8oRjShs7Z6ZrjE2KgaXikWHhVOUSgysKrOTvdoYMbxg2NJYV9RW7HQJikUGVEirlBkoKPeTudIACnllOMHU975ihq1cakoD5SIKxQZWF1vxxuOMzqtJW3veAhNCJpTcXOFYrFRIq5QZKCjvhKAzpGpIZXesSAOqxGjpv7rKMoD9ZeoUGRgdV0yw2S6iB8cD+LMkDuuUCwWJYm4EGKlEOLBDMebhRD3CSF+LIT40OzNUygWh1q7meoKE12j/oljCV3SPx7CaVUirigfik4xFEI4gGuBTGO8PwDcJqV8TgixHbhz+gWdnZ1s3bp14utt27Zx8803F7S2z+fD5XIVa/KCUu42lrt9UD423rTORlwPTNjiDce4bLmJRoteFvblolx+hrkodxvL3b40RYu4lNIL3C6EuCTD6XagJ/W5RwhRKaX0T76go6ODu+++u2hDAVwuF62trSXdu1CUu43lbh+Uj43hfUHufaWPv72kGbNR47WdA2zvjnDdltaysC8X5fIzzEW521ju9qXJG04RQnxcCLF90sfxOS7vA5alPp8h4ArFUmJzi4NYQrJvyAfAK/3jVFmM1NrNi2yZQnGEvJ64lPJ24PZc1wghLgMagLuArwohvMD2ObFQoVgktrQ5AXi+Z+z/b+9eQuMqwzCO/7+mMF6CRqiibSlddCFYsNJpMRtRcCPNQiq4EZVKUEQExVApWIdowYVQL6lCsFCExAuCF4xYigurG5FY3SiCItp2UaTQJI0Oto2PizPR6SWZmTo9eT94fjBwMjkh/yScd2a+zMxh/cqr+fbIFBtW97HMJ3uwQC76ZfeSBpq2DzR96v7/VWQWxIreCtU1fXzy/TFuW7eCwyfqbL151VJnmZ3FTzE0W8TATTdwdKrO4FuH6K0sZ2C93x/FYvEQN1vEnTdeR3XNNdRPzzHYv9bPEbdw/C6GZouoLO/h9Xs3MF0/Td8V/oemxeN74mYtpJQ8wC0sD3Ezs4x5iJuZZcxD3MwsYx7iZmYZy2qIj42NLXVCS9Ebo/dB/MbofeDGbojeNy+rIT4+Pr7UCS1Fb4zeB/Ebo/eBG7shet+8rIa4mZmdLZ17IthLbXh4+KK/4eTkJNVqtZs5XRe9MXofxG+M3gdu7IaIfbVa7bx3Xyt9iJuZWfd4OcXMLGMe4mZmGfMQNzPLmIe4mVnGwr0VbUrpemAPcBL4StJo4/oqsB2YAX4DDgAjwC/AT5J2BuvbA+yiuKH8U9JTZfR12HgSuLXxZf3ARknHA/VNAM8Ap4DP5/crQweNh4AHgePAQUnvBmjcDAwBZ4CDwEcX2i9So6TRlNJa4GVJd0frA76h6e8u6fmyGluSFOoC7AT6G9sTTdc/DqxqbH8KPAx8AIwCdwXsGwJeBV4Dtkb8HTZdfzswGK0PqAG3AD3AexF/h8D7wJWNj98J0jhIccLyZcCHC+0XrPEq4Angi6B9Fzx2IlwiLqesBn5tbE+nlHoBJI0AsymlGrAf+BJ4BHgUGEop9QTrW0cxeB4DHkgpXV5SXyeNpJQuA7ZJ2huwbwJ4G/gM+KHEvk4anwNGUkq7KQ76CI17gQrwMfD1QvtFapQ0o+Kk7DMltnXSd96xE0XEIX6E/w6GXkmz8O/Dm13APkmvUNxDOyXpb2CK8n6Wdvt+B6Yb+/1BuUtX7TYCbKMYlGVqt287UJV0B7AppVQJ2LgS2EHxyOuSL0W12fgkcFjSFmDzQvsFa1wqbfUtcOyEEO7FPo01qhcpbpG/o1h3vBa4B5gD/qJYI30TeBo4SnFLWco9yQ76ngVeAk4AP0raXUZfJ42SHkop7Qe2SJqL1kexXHYfUKc4oGoBG0co1u3PUBzgpd1LW6TxGMVy4yzwM7CveT9Jb0RrlPRCY/8JSQPR+oBNnHPslNXYSrghbmZm7Yu4nGJmZm3yEDczy5iHuJlZxjzEzcwy5iFuZpYxD3Ezs4x5iJuZZcxD3MwsY/8AFE4Lea+ZEL4AAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "plt.plot(times, strain['plus'])\n", - "plt.plot(times, strain['cross'])" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "# define parameters to inject.\n", - "# Note that a number of these (the ones on the second line) need to be there; \n", - "# bilby always looks for them even if your model ignores them.\n", - "\n", - "injection_parameters = dict(phase=0, ra=0, dec=0, psi=0, t0=0., geocent_time=0.)" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "duration = 1.0\n", - "sampling_frequency = 1024\n", - "outdir = 'outdir'\n", - "label = 'time_domain_source_model'" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [], - "source": [ - "from functools import partial" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Annoyingly I can't see an easy way of passing an argument to the model which bilby won't try to sample over.\n", - "There may be something about this in the documentation, I'll need to take a closer look. To get around this I'll use a partial function for now." - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "nr_injection_partial = partial(nr_injection, datafile=\"./rescale-h_m16_L0.18_l2m2_r300.dat\")" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "# call the waveform_generator to create our waveform model.\n", - "waveform = bilby.gw.waveform_generator.WaveformGenerator(\n", - " duration=duration, sampling_frequency=sampling_frequency,\n", - " time_domain_source_model=nr_injection_partial,\n", - " start_time=-0.5)" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/daniel/.virtualenvs/aries/bilby/lib/python3.6/site-packages/bilby/gw/detector/psd.py:356: RuntimeWarning: invalid value encountered in multiply\n", - " frequency_domain_strain = self.__power_spectral_density_interpolated(frequencies) ** 0.5 * white_noise\n", - "10:50 bilby INFO : Injected signal in H1:\n", - "10:50 bilby INFO : optimal SNR = 2.01\n", - "10:50 bilby INFO : matched filter SNR = 3.19-0.07j\n", - "10:50 bilby INFO : phase = 0\n", - "10:50 bilby INFO : ra = 0\n", - "10:50 bilby INFO : dec = 0\n", - "10:50 bilby INFO : psi = 0\n", - "10:50 bilby INFO : t0 = 0.0\n", - "10:50 bilby INFO : geocent_time = 0.0\n", - "10:50 bilby INFO : Injected signal in L1:\n", - "10:50 bilby INFO : optimal SNR = 2.37\n", - "10:50 bilby INFO : matched filter SNR = 0.97+0.12j\n", - "10:50 bilby INFO : phase = 0\n", - "10:50 bilby INFO : ra = 0\n", - "10:50 bilby INFO : dec = 0\n", - "10:50 bilby INFO : psi = 0\n", - "10:50 bilby INFO : t0 = 0.0\n", - "10:50 bilby INFO : geocent_time = 0.0\n" - ] - } - ], - "source": [ - "# inject the signal into three interferometers\n", - "ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])\n", - "ifos.set_strain_data_from_power_spectral_densities(\n", - " sampling_frequency=sampling_frequency, duration=duration,\n", - " start_time=-0.5)\n", - "ifos.inject_signal(waveform_generator=waveform,\n", - " parameters=injection_parameters);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "So, having made our injection we want to search for it using a more conventional model.\n", - "\n", - "We'll use IMRPhenomPv2. This also means that we need a list of the parameters which we'll use to produce the prior dictionary.\n", - "For simplicity I'm copying this from the Bilby 4-parameter tutorial. This could be done slightly better!" - ] - }, - { - "cell_type": "code", - "execution_count": 20, - "metadata": {}, - "outputs": [], - "source": [ - "prior_parameters = dict(\n", - " mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,\n", - " phi_12=1.7, phi_jl=0.3, luminosity_distance=15., theta_jn=0.4, psi=2.659,\n", - " phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)\n" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": {}, - "outputs": [], - "source": [ - "# Fixed arguments passed into the source model\n", - "waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',\n", - " reference_frequency=50., minimum_frequency=20.)\n", - "\n", - "\n", - "waveform_generator_bbh = bilby.gw.WaveformGenerator(\n", - " duration=duration, sampling_frequency=sampling_frequency,\n", - " frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,\n", - " waveform_arguments=waveform_arguments)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "11:11 bilby INFO : No prior given, using default BBH priors in /home/daniel/.virtualenvs/aries/bilby/lib/python3.6/site-packages/bilby/gw/prior_files/binary_black_holes.prior.\n", - "11:11 bilby INFO : Time jittering requested with non-time-marginalised likelihood, ignoring.\n", - "11:11 bilby INFO : Running for label 'label', output will be saved to 'outdir'\n", - "11:11 bilby INFO : Search parameters:\n", - "11:11 bilby INFO : mass_1 = Uniform(minimum=5, maximum=100, name='mass_1', latex_label='$m_1$', unit='$M_{\\\\odot}$', boundary=None)\n", - "11:11 bilby INFO : luminosity_distance = UniformSourceFrame(minimum=100.0, maximum=5000.0, cosmology=FlatLambdaCDM(name=\"Planck15\", H0=67.7 km / (Mpc s), Om0=0.307, Tcmb0=2.725 K, Neff=3.05, m_nu=[0. 0. 0.06] eV, Ob0=0.0486), name='luminosity_distance', latex_label='$d_L$', unit=Unit(\"Mpc\"), boundary=None)\n", - "11:11 bilby INFO : theta_jn = Sine(name='theta_jn', latex_label='$\\\\theta_{JN}$', unit=None, minimum=0, maximum=3.141592653589793, boundary='reflective')\n", - "11:11 bilby INFO : mass_ratio = Constraint(minimum=0.125, maximum=1, name='mass_ratio', latex_label='$q$', unit=None)\n", - "11:11 bilby INFO : mass_2 = 29.0\n", - "11:11 bilby INFO : a_1 = 0.4\n", - "11:11 bilby INFO : a_2 = 0.3\n", - "11:11 bilby INFO : tilt_1 = 0.5\n", - "11:11 bilby INFO : tilt_2 = 1.0\n", - "11:11 bilby INFO : phi_12 = 1.7\n", - "11:11 bilby INFO : phi_jl = 0.3\n", - "11:11 bilby INFO : dec = -1.2108\n", - "11:11 bilby INFO : ra = 1.375\n", - "11:11 bilby INFO : psi = 2.659\n", - "11:11 bilby INFO : phase = 1.3\n", - "11:11 bilby INFO : geocent_time = 1126259642.413\n", - "11:11 bilby INFO : Single likelihood evaluation took 1.556e-03 s\n", - "11:11 bilby INFO : Using sampler Dynesty with kwargs {'bound': 'multi', 'sample': 'rwalk', 'verbose': True, 'periodic': None, 'check_point_delta_t': 600, 'nlive': 1000, 'first_update': None, 'walks': 30, 'npdim': None, 'rstate': None, 'queue_size': None, 'pool': None, 'use_pool': None, 'live_points': None, 'logl_args': None, 'logl_kwargs': None, 'ptform_args': None, 'ptform_kwargs': None, 'enlarge': None, 'bootstrap': None, 'vol_dec': 0.5, 'vol_check': 2.0, 'facc': 0.5, 'slices': 5, 'update_interval': 600, 'print_func': >, 'dlogz': 0.1, 'maxiter': None, 'maxcall': None, 'logl_max': inf, 'add_live': True, 'print_progress': True, 'save_bounds': False}\n", - "11:11 bilby INFO : Checkpoint every n_check_point = 400000\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 7042| logz ratio=-9.391 +/- 0.081 | dlogz: 0.354 > 0.100000" - ] - } - ], - "source": [ - "# create the priors\n", - "priors = bilby.gw.prior.BBHPriorDict()\n", - "priors['geocent_time'] = bilby.core.prior.Uniform(\n", - " minimum=injection_parameters['geocent_time'] - 1,\n", - " maximum=injection_parameters['geocent_time'] + 1,\n", - " name='geocent_time', latex_label='$t_c$', unit='$s$')\n", - "\n", - "for key in ['mass_2', 'a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi', 'ra',\n", - " 'dec', 'geocent_time', 'phase']:\n", - " priors[key] = prior_parameters[key]\n", - "\n", - "# define likelihood\n", - "likelihood = bilby.gw.likelihood.GravitationalWaveTransient(ifos, waveform_generator_bbh)\n", - "\n", - "# launch sampler\n", - "result = bilby.core.sampler.run_sampler(\n", - " likelihood, priors, sampler='dynesty', npoints=1000,\n", - " )\n", - "\n", - "result.plot_corner()" - ] - }, - { - "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.6.6" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} -- GitLab