<jupyter_start><jupyter_text>Introduction to PyCBC Inference 3: Results files and plotting Collin Capano In this tutorial, we take a closer look at the output files, including how to load samples from them. We then illustrate more advanced options in `pycbc_inference_plot_posterior`, and show how to create a posterior file. We will use an already complete results file.  Prerequisites<jupyter_code>import sys
!{sys.executable} -m pip install pycbc lalsuite ligo-common --no-cache-dir

# This is needed to access the executables on sciserver. On a personal machine this should be ignore.
path = %env PATH
%env PATH=$path:/home/idies/miniconda3/envs/py37/bin 
from __future__ import print_function
import os
import matplotlib
%matplotlib inline
from matplotlib import pyplot

from IPython.display import Image<jupyter_output><empty_output><jupyter_text>Get the file We will download a result file from a fully completed analysis, from [here](https://www.atlas.aei.uni-hannover.de/~work-cdcapano/pycbc_workshop_june_2020). This is the result of running `emcee_pt` on GW150914 using the standard prior provided in the online [PyCBC Inference documentation](http://pycbc.org/pycbc/latest/html/inference/examples/gw150914.html). You can see the complete configuration file [here](https://www.atlas.aei.uni-hannover.de/~work-cdcapano/pycbc_workshop_june_2020/inference-emcee_pt.ini).Note that there are two files in this directory. This is because this run was run twice with different starting seeds, to accumulate more samples. We'll download one of the files now.<jupyter_code>if not os.path.exists('H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf'):
    !wget https://www.atlas.aei.uni-hannover.de/~work-cdcapano/pycbc_workshop_june_2020/H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf<jupyter_output><empty_output><jupyter_text>The results file The results file `H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf` is an HDF file. We could use the standard python module for reading HDF files [h5py](http://docs.h5py.org/en/stable/) to examine it. However, there are classes in [pycbc.inference.io](http://pycbc.org/pycbc/latest/html/pycbc.inference.io.html) that inherit from [h5py.File](http://docs.h5py.org/en/stable/high/file.htmlfile-objects) and add several convenience functions that make it easier to read samples from the file. There is one class for each type of sampler. To load the file, we will use [pycbc.inference.io.loadfile](http://pycbc.org/pycbc/latest/html/pycbc.inference.io.htmlpycbc.inference.io.loadfile). This function automatically determines which class to use to read the file, based on what is in the file.<jupyter_code>from pycbc.inference.io import loadfile
fp = loadfile('H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf', 'r')<jupyter_output><empty_output><jupyter_text>The top-level groups in the HDF file:<jupyter_code>list(fp.keys())<jupyter_output><empty_output><jupyter_text>HDF files also have `attrs`. These are basically dictionaries that can be used to store metadata. Let's look at what metadata `pycbc_inference` stored in the file:<jupyter_code>list(fp.attrs.items())<jupyter_output><empty_output><jupyter_text>Sampler infoThe `sampler_info` group contains information about the sampler. The data that's stored in this group is sampler-specific. This is what was stored for this run, which used `emcee_pt`:<jupyter_code>list(fp['sampler_info'].keys())<jupyter_output><empty_output><jupyter_text>We see that for `emcee_pt` (as well as any MCMC sampler), burn in information and the autocorrelation time (ACT) are stored in the sampler info. We could read these directly from the file. For example:<jupyter_code>print(fp['sampler_info/burn_in_iteration'][()])<jupyter_output><empty_output><jupyter_text>However, since the ACT and burn in itertion are such important information, they are promotted to attributes. We can get them by simply doing:<jupyter_code>print(fp.act, fp.burn_in_iteration)<jupyter_output><empty_output><jupyter_text>Every group can have their own `.attrs`; let's look at the `sampler_info` group's `attrs`:<jupyter_code>list(fp['sampler_info'].attrs.items())<jupyter_output><empty_output><jupyter_text>SamplesThe samples that the sampler produced are stored in the `samples` group. There is one data set for each variable parameter.<jupyter_code>list(fp['samples'].keys())<jupyter_output><empty_output><jupyter_text>Additional data may be stored in the `samples` group `attrs`. In this case, because we used the `marinalized_phase` model, the log likelihood of the noise $\log p(d|n) = -\left/2$ is stored.<jupyter_code>list(fp['samples'].attrs.items())<jupyter_output><empty_output><jupyter_text>Lets look at the shape of the datasets in the samples group. They all have the same shape, so we can just look at one of the parameters.<jupyter_code>print(fp['samples/mass1'].shape)<jupyter_output><empty_output><jupyter_text>The shape of the dataset is `ntemps x nwalkers x n thinned iteration`. This run used 20 temperatures and 200 walkers. Due to `max-samples-per-chain` being set to 1024, the data set has been thinned to include 704 samples from each walker and temperature.If this had been an `emcee` run (which uses no temperatures), the samples datasets would have been two dimensional: `nwalkers x niterations`. If it had been a nested sampling run (e.g., with `dynesty`), the datasets would have been 1 dimensional. **The format of the samples data is sampler dependent.** This is why we have separate classes to read the results file.  Load samplesNow lets load the samples. We can use the `read_samples` function to do this. This takes as a the first argument a list of parameters to load. Here, we'll load all of the parameters. It also takes additional keyword arguments. These arguments are sampler-specific. For `emcee_pt` if we provide no additional keyword arguments, we'll get all temperatures. We just want the coldest temperature `temp=0`, as that is the posterior.If we provide no other arguments, `read_samples` will load all of the independent samples post burn-in. That is, it will get samples from all walkers, starting from the burn in iteration, and thinned by the autocorrelation time. The samples are flattened into a 1D array.<jupyter_code>samples = fp.read_samples(fp['samples'].keys(), temps=0)
print(samples.shape)<jupyter_output><empty_output><jupyter_text>So, we have 1600 independent samples. What sort of object is `samples`?<jupyter_code>type(samples)<jupyter_output><empty_output><jupyter_text>The samples have been returned as a [FieldArray](http://pycbc.org/pycbc/latest/html/pycbc.io.htmlpycbc.io.record.FieldArray). A `FieldArray` is a wrapper around [`numpy` structured arrays](https://docs.scipy.org/doc/numpy/user/basics.rec.html) (actually, it inherits from `numpy`'s [record array](https://docs.scipy.org/doc/numpy/user/basics.rec.htmlrecord-arrays), which in turn inherits from structured arrays). Structured arrays allow multiple parameters to be stored in a single array. `FieldArray` adds several convenience features, which we will look at in a bit. Since they inherit from structured arrays, a `FieldArray` can be indexed by parameter name:<jupyter_code>samples['mass1']<jupyter_output><empty_output><jupyter_text>When you do, you get back a normal numpy array of the values stored for that parameter. We can see the full list of parameters that are stored using `.fieldnames` (This is an attribute added by `FieldArray`; it does not work for general structured arrays.)<jupyter_code>print(samples.fieldnames)<jupyter_output><empty_output><jupyter_text>Plot `mass1` and `mass2` Let's make a scatter plot of the two mass parameters.<jupyter_code>mass1 = samples['mass1']
mass2 = samples['mass2']
fig, ax = pyplot.subplots()
ax.scatter(mass1, mass2)
ax.set_xlabel('mass1')
ax.set_ylabel('mass2')
fig.show()<jupyter_output><empty_output><jupyter_text>It looks like it converged, but something is odd here... we normally define `mass1 >= mass2`. But it looks like in our posterior samples, some points have `mass2 > mass1`. Lets confirm:<jupyter_code>(mass1 >= mass2).all()<jupyter_output><empty_output><jupyter_text>**The masses in the output file do not respect the convention that `mass1 >= mass2`.** This is expected: our prior (which was uniform for both masses in range $[10, 80)\,\mathrm{M}_\odot$) was symmetric between the masses.But for plotting purposes, we want `mass1 >= mass2`. To do that, we can take advantage of one `FieldArray`'s key features. The [pycbc.conversions](http://pycbc.org/pycbc/latest/html/pycbc.htmlmodule-pycbc.conversions) provides a number of useful functions for converting parameters.<jupyter_code>from pycbc import conversions<jupyter_output><empty_output><jupyter_text>Included in the conversions module, the [primary_mass](http://pycbc.org/pycbc/latest/html/pycbc.htmlpycbc.conversions.primary_mass) function pulls out the larger of the two masses in array of masses. Likewise, the [secondary_mass](http://pycbc.org/pycbc/latest/html/pycbc.htmlpycbc.conversions.secondary_mass) function pulls out the smaller of the two. This can do what we want:<jupyter_code>mass1 = conversions.primary_mass(samples['mass1'], samples['mass2'])
mass2 = conversions.secondary_mass(samples['mass1'], samples['mass2'])
(mass1 >= mass2).all()
fig, ax = pyplot.subplots()
ax.scatter(mass1, mass2)
ax.set_xlabel('mass1')
ax.set_ylabel('mass2')
fig.show()<jupyter_output><empty_output><jupyter_text>Using functions with Field Arrays We did not need to import and call the `primary_mass` and `secondary_mass` functions to apply them to the samples. Alternatively, we could of just passed these functions to the `samples` `FieldArray` as if they were parameters:<jupyter_code>mass1 = samples['primary_mass(mass1, mass2)']
mass2 = samples['secondary_mass(mass1, mass2)']
(mass1 >= mass2).all()<jupyter_output><empty_output><jupyter_text>This highlights the main purpose of `FieldArray`: * **`FieldArrays` can be passed functions of their parameters as strings.**The available functions are any function in the [pycbc.conversions](http://pycbc.org/pycbc/latest/html/pycbc.htmlmodule-pycbc.conversions), [pycbc.coordinates](http://pycbc.org/pycbc/latest/html/pycbc.htmlmodule-pycbc.coordinates), or [pycbc.cosmology](http://pycbc.org/pycbc/latest/html/pycbc.htmlmodule-pycbc.cosmology) modules. In addition, numpy math functions are available for use. The syntax is standard python.For example, if we wanted `spin1z`:<jupyter_code>samples['spin1_a * cos(spin1_polar)']<jupyter_output><empty_output><jupyter_text>Or redshift:<jupyter_code>samples['redshift(distance)']<jupyter_output><empty_output><jupyter_text>The `--parameters` option The main advantage of `FieldArrays` is that it means functions of parameters can be passed on the command line to plotting programs using the `--parameters` option. For example:<jupyter_code>!pycbc_inference_plot_posterior \
    --input-file H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf \
    --output-file mass1_mass2.png \
    --parameters 'primary_mass(mass1, mass2):mass1' 'secondary_mass(mass1, mass2):mass2' \
    --plot-scatter \
    --plot-contours \
    --plot-marginal \
    --z-arg snr
Image(filename='mass1_mass2.png', width=640, height=480)<jupyter_output><empty_output><jupyter_text>Basically, any program with a `--parameters` option can take functions of parameters. Note that this option also allows you to pass a label for the parameters. The syntax for this is `--parameters parameter_to_plot:label`. If a label is provided that is recognized parameter in the [pycbc.waveform.waveform.parameters](https://github.com/gwastro/pycbc/blob/master/pycbc/waveform/parameters.py) module, the predefined label will be used for that parameter in the plot. This is why when we passed `mass1` and `mass2` as the labels, we got $m_1$ $m_2$ in the plot.  File helpLike all executables in `pycbc`, `pycbc_inference_plot_posterior` has a help message for its available options. The help can be seen by running the program with `--help` (or `-h`):<jupyter_code>!pycbc_inference_plot_posterior --help<jupyter_output><empty_output><jupyter_text>It can be difficult to remember all of the functions and parameters that you have available to you. For this reason, if you give the option `--file-help` (or `-H`), along with an input file, you will get a file-specific help message that tells you what parameters are available for plotting in the file, along with all of the functions you can use:<jupyter_code>!pycbc_inference_plot_posterior --input-file H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf --file-help<jupyter_output><empty_output><jupyter_text>Note that this also reveals additional options that you can provide to the program. These options are sampler-specific (and thus file-specific), which is why they do not appear in the general help message. For example, passing a `--temps` option works for a file that was created by `emcee_pt`, but not for a file that was created by `emcee`. If you're ever unsure what is available, just run the `-H` option with your input file.  Challenge:Make a posterior plot of chirp mass ("mchirp") vs eta.  Creating posterior files with `extract_samples` It can get very tedious typing long functions on the command line every time you want to create a plot. Also, you may want to release your posterior samples along with a paper. Having different file formats depending on the sampler would make such a release confusing. Plus, the full samples file contains much more information than is necessary if you just want the posterior samples.The program `pycbc_inference_extract_samples` can convert a samples file into a *posterior file*. Regardless of the sampler that was used: * **Posterior files all have a "samples" group that contains 1D datasets of the posterior samples.** In process of creating the file, `pycbc_inference_extract_samples` can also be used to write out functions of the parameters as new datasets. This is because it also takes a `--parameters` option, except that here the labels part gives the name of the parameter in the new file. Nuisance parameters can be excluded, along with the other groups, such as `sampler_info`.For example, in the following, we create a posterior file containing the component masses (which we sort such that `mass1` $\geq$ `mass2`) and the log likelihood. We also add the chirp mass and symmetric mass ratio `eta`.<jupyter_code>!pycbc_inference_extract_samples --verbose \
    --input-file H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf \
    --output-file mass_posterior.hdf \
    --parameters \
        'primary_mass(mass1, mass2):mass1' \
        'secondary_mass(mass1, mass2):mass2' \
        'mchirp_from_mass1_mass2(mass1, mass2):mchirp' \
        'eta_from_mass1_mass2(mass1, mass2):eta' \
        loglikelihood \
    --skip-groups data sampler_info<jupyter_output><empty_output><jupyter_text>Let's load the posterior file and examine it.<jupyter_code>fp = loadfile('mass_posterior.hdf', 'r')
list(fp.keys())
list(fp['samples'].keys())
fp['samples']['mchirp'].shape
list(fp.attrs.items())
posterior_samples = fp.read_samples(['mass1', 'mass2'])
mass1 = posterior_samples['mass1']
mass2 = posterior_samples['mass2']
(mass1 >= mass2).all()<jupyter_output><empty_output><jupyter_text>Note that the posterior file is much smaller than the full samples file:<jupyter_code>!ls -lh H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf mass_posterior.hdf<jupyter_output><empty_output><jupyter_text>Now lets create a posterior plot of all of the mass parameters using the posterior file. We no longer need to pass the functions on the command line:<jupyter_code>!pycbc_inference_plot_posterior \
    --input-file mass_posterior.hdf \
    --output-file posterior_mass1_mass2.png \
    --parameters mass1 mass2 mchirp eta \
    --plot-scatter \
    --plot-contours \
    --plot-marginal \
    --z-arg snr
Image(filename='posterior_mass1_mass2.png', width=640, height=480)<jupyter_output><empty_output><jupyter_text>Extracting all parameters using `'*'`When we provided specific parameters to extract to `pycbc_inference_extract_samples`, only the parameters we specified were extracted. We can get all additionaly parameters by passing `'*'` to the `--parameters` option. Let's try it:<jupyter_code>!pycbc_inference_extract_samples --verbose \
    --input-file H1L1-INFERENCE_EMCEE_PT_0-1126259200-400.hdf \
    --output-file posterior.hdf \
    --parameters \
        'primary_mass(mass1, mass2):mass1' \
        'secondary_mass(mass1, mass2):mass2' \
        'mchirp_from_mass1_mass2(mass1, mass2):mchirp' \
        'eta_from_mass1_mass2(mass1, mass2):eta' \
        '*' \
    --skip-groups data sampler_info --force<jupyter_output><empty_output><jupyter_text>Let's load the file and check that we have all the parameters in the samples group:<jupyter_code>fp = loadfile('posterior.hdf', 'r')
print(sorted(fp['samples'].keys()))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09650109525782961, 'fraction_numerical': 0.018767390918240482, 'mean_word_length': 5.5984375, 'pattern_counts': {'":': 0, '<': 129, '<?xml version=': 0, '>': 137, 'https://': 6, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '284748', 'license': 'MIT', 'n_tokens_mistral': 5429, 'n_tokens_neox': 5232, 'n_words': 1757, 'path': 'PyCBC-Tutorials-master/tutorial/inference_5_results_io/IntroToPyCBCInference.ipynb', 'repo_name': 'basuparth/ICERM_Workshop', 'size': '16891'}