example.py

import numpy as np
import matpotlib.pyplot as plt
import corner   # From: https://github.com/dfm/corner.py


SKIP = 100   # Only plot 1/SKIP data points
fname = "./lzk-gw-data_cornish-becsy.npz"


data = np.load(fname)

# Load the number of binaries in each realization
# Note that each binary-data array is shaped as (10, 6905804), for 10 realizations, and a
#      maximum number of binaries 6905804.  In realizations with fewer binaries, the remaining
#      entires in each array are zero.
num = data['num']
reals = len(num)
print("{} Realizations with {:.2e} median binaries".format(reals, np.median(num)))

keys = ["mch", "dlum", "fobs", "frst"]
names = ["Chirp-Mass [Msol]", "Luminosity Distance [Gpc]", "f_obs [1/yr]", "f_rst [1/yr]"]
units = [1.99e+33, 3.09e+27, 3.17e-08, 3.17e-08]

# Select a random realization to plot
real = np.random.choice(reals)

vals = [data[vv][real, :num[real]] for vv in keys]
# Downsample data
if skip is not None:
    vals = [vv[::skip] for vv in vals]
# Convert to log-space
vals = [np.log10(vv/uu) for vv, uu in zip(vals, units)]

# Create figure
fig, axes = plt.subplots(figsize=[12, 12], ncols=len(keys), nrows=len(keys))
plt.subplots_adjust(hspace=0.2, wspace=0.2)
# Create corner plot using `https://github.com/dfm/corner.py`
corner.corner(np.asarray(vals).T, labels=names, fig=fig)

plt.show()