{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Sample release for GW190425\n",
    "\n",
    "This notebook serves as a basic introduction to loading and viewing data released in associaton with the publication titled\n",
    "__GW190425: Observation of a compact binary coalescence with total mass $\\sim 3.4 M_{\\odot}$__ avaliable\n",
    "through [DCC](https://dcc.ligo.org/LIGO-P190425/public) and [arXiv](https://arxiv.org/abs/2001.01761).\n",
    "\n",
    "The data used in these tutorials will be downloaded from the public DCC page [LIGO-P2000026](https://dcc.ligo.org/P2000026/public)."
   ]
  },
  {
   "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.3.0` 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"
   ]
  },
  {
   "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 = \"GW190425_posterior_samples.h5\"\n",
    "data = h5py.File(fn,'r')\n",
    "\n",
    "# print out top-level data structures\n",
    "print('Top-level data structures:',data.keys())\n",
    "# print out parametrized waveform family names (\"approximants\" in LIGO jargon).\n",
    "# HS, LS = high-spin prior and low-spin prior, respectively\n",
    "print('approximants:',data['approximant'].keys())\n",
    "\n",
    "# extract posterior samples for one of the approximants\n",
    "posterior_samples = data['posterior_samples']['PhenomDNRT-HS']\n",
    "print('data structures in posterior_samples:',posterior_samples.keys())\n",
    "pnames = [item.decode(\"utf-8\") for item in posterior_samples['parameter_names']]\n",
    "print('parameter names:',pnames)\n",
    "\n",
    "# extract the samples data into an numpy array:\n",
    "samples = np.array(posterior_samples['samples']).T\n",
    "# get samples for one of the parameters\n",
    "ind = pnames.index('luminosity_distance')\n",
    "dL = samples[ind]\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.3.0 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 6 different approximants that were used to analyze __GW190425__ and they are all stored in the data file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "fn = \"GW190425_posterior_samples.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'll pick one approximant by random and plot its respective data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "samples_dict = data.samples_dict\n",
    "posterior_samples = samples_dict[\"PhenomPNRT-HS\"]\n",
    "prior_samples = data.priors[\"samples\"][\"PhenomPNRT-HS\"]\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 for a single waveform 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.core.plots.plot import _1d_histogram_plot\n",
    "from pesummary.gw.plots.latex_labels import GWlatex_labels\n",
    "\n",
    "parameter = \"luminosity_distance\"\n",
    "latex_label = GWlatex_labels[parameter]\n",
    "\n",
    "fig = _1d_histogram_plot(\n",
    "    parameter, posterior_samples[parameter], latex_label, prior=prior_samples[parameter]\n",
    ")\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",
    "colors = ['b', 'r', 'k', 'y', 'orange', 'g']\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": {
    "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)\n",
    "plt.show()\n",
    "#plt.savefig(fn+'_corner.png')"
   ]
  },
  {
   "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": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from pesummary.gw.plots.plot import _psd_plot\n",
    "\n",
    "psd = data.psd[\"PhenomPNRT-HS\"]\n",
    "ifos = list(psd.keys())\n",
    "frequencies, strains = [], []\n",
    "for ifo in ifos:\n",
    "    frequencies.append(np.array(psd[ifo]).T[0])\n",
    "    strains.append(np.array(psd[ifo]).T[1])\n",
    "fig = _psd_plot(frequencies, strains, labels=ifos, fmin=19.4)\n",
    "fig.set_size_inches(12, 8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "The calibration envelopes that were used in this analysis can also be extracted from this file and plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from pesummary.gw.plots.plot import _calibration_envelope_plot\n",
    "\n",
    "prior = data.priors[\"calibration\"][\"PhenomPNRT-HS\"]\n",
    "calibration = data.calibration[\"PhenomPNRT-HS\"]\n",
    "frequencies = np.arange(20., 1024., 1. / 4)\n",
    "calibration_data, prior_data = [], []\n",
    "for ifo in ifos:\n",
    "    calibration_data.append(np.array(calibration[ifo]))\n",
    "    prior_data.append(np.array(prior[ifo]))\n",
    "fig = _calibration_envelope_plot(frequencies, calibration_data, ifos, prior=prior_data)\n",
    "fig.set_size_inches(16.5, 10.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "The configuration file that were used for each analysis can also be extracted from this file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "config = data.config[\"PhenomPNRT-HS\"]\n",
    "for i in config.keys():\n",
    "    print(\"[{}]\".format(i))\n",
    "    for key, item in config[i].items():\n",
    "        print(\"{}={}\".format(key, item[0].decode(\"utf-8\")))\n",
    "    print(\"\\n\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mypy37",
   "language": "python",
   "name": "mypy37"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}