{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Sample release for GW190412\n",
    "\n",
    "This notebook serves as a basic introduction to loading and viewing data released in associaton with the publication titled\n",
    "__GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses__ avaliable\n",
    "through [DCC](https://dcc.ligo.org/LIGO-P190412/public).\n",
    "\n",
    "The data used in these tutorials can be downloaded from the same [DCC page](https://dcc.ligo.org/LIGO-P190412/public).\n",
    "\n",
    "__Note:__ This notebook has been prepared for the released samples that are available since DCC version 12. Preliminary samples were distributed with preprint versions of the paper. The samples initially released with the published paper contained a minor bug that only affected the calculation of the network matched-filter SNR. This bug has been fixed with the current release."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "The released data file can be read in using the `PESummary` or `h5py` libraries. For this notebook we'll start with simple stuff using h5py. Then we'll use `PESummary v0.9.1` to read the data  files as well as for plotting. For general instructions on how to manipulate the data file and/or read this data file with `h5py`, see the [PESummary docs](https://lscsoft.docs.ligo.org/pesummary/data/reading_the_metafile.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# import useful python packages\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import h5py\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There is a known incompatibility between some older versions of numpy/scipy and seaborn. \n",
    "# Make sure you have an up-to-date version of scipy installed.\n",
    "# This notebook has been tested with the following versions\n",
    "for name, package in zip(('numpy', 'scipy', 'seaborn', 'h5py'), (np, scipy, sns, h5py)):\n",
    "    print('{} version: {}'.format(name, package.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Some simple stuff with \"vanilla\" h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# read in the data\n",
    "fn = 'GW190412_posterior_samples_v3.h5'\n",
    "data = h5py.File(fn,'r')\n",
    "\n",
    "# print out parametrized waveform family names (\"approximants\" in LIGO jargon).\n",
    "print('approximants:',data.keys())\n",
    "\n",
    "# print out top-level data structures for one approximant. Here fore example we use the combined samples\n",
    "# between IMRPhenomPv3HM and SEOBNRv4PHM. The data structure is the same for all approximants.\n",
    "print('Top-level data structures:',data['combined'].keys())\n",
    "\n",
    "# extract posterior samples for one of the approximants\n",
    "posterior_samples = data['combined']['posterior_samples']\n",
    "print('data structures in posterior_samples:',posterior_samples.dtype)\n",
    "pnames = [item for item in posterior_samples.dtype.names]\n",
    "print('parameter names:',pnames)\n",
    "\n",
    "# get samples for one of the parameters\n",
    "dL = posterior_samples['luminosity_distance']\n",
    "print('dL shape, mean, std =',dL.shape,dL.mean(),dL.std())\n",
    "\n",
    "# smooth it\n",
    "from scipy.stats.kde import gaussian_kde\n",
    "hs = gaussian_kde(dL)\n",
    "\n",
    "# histogram, and overlay the smoothed PDF\n",
    "plt.figure()\n",
    "h, b, o = plt.hist(dL,bins=100)\n",
    "hsmoothed = hs(b)*len(dL)*(b[1]-b[0])\n",
    "plt.plot(b,hsmoothed)\n",
    "plt.xlabel('luminosity distance')\n",
    "plt.ylabel('posterior PDF')\n",
    "plt.show()\n",
    "\n",
    "# release memory for the data\n",
    "#del data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Now use PESummary v0.9.1 to read the data files as well as for plotting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# import ligo-specific python packages. \n",
    "# pesummary is a ligo-specific python package for reading and plotting the results of Bayesian parameter estimation.\n",
    "# Install with \"pip install pesummary\" , and make sure you have version >= 0.3.0.\n",
    "import pesummary\n",
    "from pesummary.gw.file.read import read\n",
    "print(pesummary.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "There are 9 different approximants that were used to analyze __GW190412__ plus the combined samples of IMRPhenomPv3HM and SEOBNRv4PHM. They are all stored in the data file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "fn = 'GW190412_posterior_samples_v3.h5'\n",
    "data = read(fn)\n",
    "labels = data.labels\n",
    "print(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "To illustrate the data structure we pick the combined posteriors and plot the respective data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "samples_dict = data.samples_dict\n",
    "posterior_samples = samples_dict[\"combined\"]\n",
    "prior_samples = data.priors[\"samples\"][\"combined\"]\n",
    "parameters = posterior_samples.keys()\n",
    "print(parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "As an example, we'll show the different posterior distributions derived from the combined (precession + higher multipoles) data and the posterior distribution derived using the different approximants for the `luminosity_distance` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from pesummary.gw.plots.latex_labels import GWlatex_labels\n",
    "\n",
    "parameter = \"luminosity_distance\"\n",
    "latex_label = GWlatex_labels[parameter]\n",
    "\n",
    "fig = posterior_samples.plot(parameter, type=\"hist\", kde=False, prior=prior_samples[parameter])\n",
    "fig.set_size_inches(12, 8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from pesummary.core.plots.plot import _1d_comparison_histogram_plot\n",
    "\n",
    "samples = []\n",
    "for label in labels:\n",
    "    samples.append(samples_dict[label][parameter])\n",
    "\n",
    "\n",
    "colors = ['b', 'r', 'k', 'y', 'orange', 'g','purple','cyan','grey','violet']\n",
    "fig = _1d_comparison_histogram_plot(parameter, samples, colors, latex_label, labels, kde=True)\n",
    "fig.set_size_inches(12, 8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is another examples for the mass ratio of the combined data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameter = \"mass_ratio\"\n",
    "latex_label = GWlatex_labels[parameter]\n",
    "\n",
    "fig = posterior_samples.plot(parameter, type=\"hist\", kde=True, prior=prior_samples[parameter])\n",
    "fig.set_size_inches(12, 8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Make a corner plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from pesummary.gw.plots.plot import _make_corner_plot\n",
    "\n",
    "fig = _make_corner_plot(posterior_samples, GWlatex_labels, corner_parameters=['mass_1_source', 'mass_2_source', 'luminosity_distance', 'theta_jn', 'chi_eff', 'chi_p'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "The psds that were used for each analysis can also be extracted from this file and plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psd = data.psd[\"combined\"]\n",
    "fig = psd.plot(fmin=20.0)\n",
    "fig.set_size_inches(12, 8)\n",
    "plt.show()"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}