From 4704d4e051dd703ecc97bf41fbc997cb420746f7 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Thu, 1 Feb 2024 10:28:35 -0800 Subject: [PATCH] Reorganize and remove some 'analysis' stuff. --- analysis/limit-curve_witt.ipynb | 789 ------------------ analysis/pop-data.ipynb | 173 ---- .../librarian/posterior_populations.py | 0 .../nanograv-15yr-populations.ipynb | 0 4 files changed, 962 deletions(-) delete mode 100644 analysis/limit-curve_witt.ipynb delete mode 100644 analysis/pop-data.ipynb rename analysis/gen-15yr-pops/gen_holodeck_pops.py => holodeck/librarian/posterior_populations.py (100%) rename analysis/gen-15yr-pops/holodeck-pops.ipynb => notebooks/nanograv-15yr-populations.ipynb (100%) diff --git a/analysis/limit-curve_witt.ipynb b/analysis/limit-curve_witt.ipynb deleted file mode 100644 index 87e59f28..00000000 --- a/analysis/limit-curve_witt.ipynb +++ /dev/null @@ -1,789 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# %load ../notebooks/init.ipy\n", - "%reload_ext autoreload\n", - "%autoreload 2\n", - "\n", - "# Builtin packages\n", - "from importlib import reload\n", - "import logging\n", - "import os\n", - "from pathlib import Path\n", - "import sys\n", - "import warnings\n", - "\n", - "# standard secondary packages\n", - "import astropy as ap\n", - "import h5py\n", - "import matplotlib as mpl\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import scipy as sp\n", - "import scipy.stats\n", - "import tqdm.notebook as tqdm\n", - "\n", - "# development packages\n", - "import kalepy as kale\n", - "import kalepy.utils\n", - "import kalepy.plot\n", - "\n", - "# --- Holodeck ----\n", - "import holodeck as holo\n", - "import holodeck.sam\n", - "from holodeck import cosmo, utils, plot\n", - "from holodeck.constants import MSOL, PC, YR, MPC, GYR\n", - "import holodeck.gravwaves\n", - "import holodeck.evolution\n", - "import holodeck.population\n", - "\n", - "# Silence annoying numpy errors\n", - "np.seterr(divide='ignore', invalid='ignore', over='ignore')\n", - "warnings.filterwarnings(\"ignore\", category=UserWarning)\n", - "\n", - "# Plotting settings\n", - "mpl.rc('font', **{'family': 'serif', 'sans-serif': ['Times'], 'size': 15})\n", - "mpl.rc('lines', solid_capstyle='round')\n", - "mpl.rc('mathtext', fontset='cm')\n", - "mpl.style.use('default') # avoid dark backgrounds from dark theme vscode\n", - "plt.rcParams.update({'grid.alpha': 0.5})\n", - "\n", - "# Load log and set logging level\n", - "log = holo.log\n", - "log.setLevel(logging.INFO)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from holodeck.constants import NWTG, SPLC, GPC\n", - "\n", - "LABEL_FREQ_YRS = 'Frequency $[\\mathrm{yr}^{-1}]$'\n", - "LABEL_STRAIN_AMP = 'GW Strain Amplitude'\n", - "LABEL_DIST_LUM_MPC = 'Luminosity Distance $[\\mathrm{Mpc}]$'\n", - "LABEL_NDENS_MPC3 = 'Number Density $[\\mathrm{cMpc}^{-3}]$'\n", - "LABEL_MASS_TOTAL = 'Total Mass $[M_\\odot]$'\n", - "\n", - "def _twin_redz_from_dlum_mpc(ax, zcutoff=0.01):\n", - " yticks = np.array(ax.get_yticks())\n", - " ytick_labels = ax.get_yticklabels()\n", - "\n", - " zticks = cosmo.dlum_to_z(yticks * MPC)\n", - " ztick_labels = [f\"{zz:.2f}\" if zz > zcutoff else \"\" for zz in zticks]\n", - "\n", - " tw = ax.twinx()\n", - " tw.set_yticks(yticks)\n", - " tw.set(yscale='log', ylim=ax.get_ylim(), ylabel='Redshift')\n", - " tw.set_yticklabels(ztick_labels)\n", - " return tw\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Load Data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "path = Path(\"/Users/lzkelley/programs/nanograv/holodeck/data/caitlin-witt_2022-08-05_12p5yr_final-curves\")\n", - "fname = \"CW_fixCRN_UL_final.txt\"\n", - "path = path.joinpath(fname)\n", - "data = np.loadtxt(path)\n", - "\n", - "ul_freq = data[:, 0] # hz\n", - "ul_strain = data[:, 1] # hs? hc?\n", - "ul_strain_err = data[:, 2]\n", - "fig, ax = plot.figax()\n", - "\n", - "# freq = ul_freq\n", - "# freq = ul_freq * YR\n", - "NUM = 21\n", - "UNITS = 1/YR\n", - "\n", - "num = (NUM + 1) // 2\n", - "assert num >= 2\n", - "\n", - "ax.plot(ul_freq/UNITS, ul_strain, lw=1.5, alpha=0.5, color='k')\n", - "ax.errorbar(ul_freq/UNITS, ul_strain, yerr=ul_strain_err, ls='none', alpha=0.5, lw=1.0, color='r')\n", - "\n", - "imin = np.argmin(ul_strain)\n", - "xmin = ul_freq[imin]\n", - "ymin = ul_strain[imin]\n", - "print(f\"{ymin:.2e} at {UNITS/xmin:.2f} yr ({xmin/UNITS:.2e} 1/yr = {xmin:.2e} Hz)\")\n", - "ax.axvline(xmin/UNITS, color='r', alpha=0.25)\n", - "ax.axhline(ymin, color='r', alpha=0.25)\n", - "\n", - "xr = kale.utils.spacing([xmin, 100*xmin], scale='log', num=num)\n", - "yr = ymin * np.power(xr/xmin, 5.0/6.0)\n", - "ax.plot(xr/UNITS, yr, 'b--', alpha=0.25)\n", - "\n", - "# powers = [1.0, 5.0/6.0, 2.0/3.0, 0.5]\n", - "# for pp in powers:\n", - "# yr = ymin * np.power(xr/xmin, pp)\n", - "# ax.plot(xr, yr, ls='--', alpha=0.5, label=f'{pp:.4f}')\n", - "\n", - "xl = kale.utils.spacing([xmin, xmin/10.0], scale='log', num=num)\n", - "yl = ymin * np.power(xl/xmin, -1.0)\n", - "ax.plot(xl/UNITS, yl, 'b--', alpha=0.25)\n", - "# powers = [3.0/2.0, 1.0, 5.0/6.0]\n", - "# for pp in powers:\n", - "# yl = ymin * np.power(xl/xmin, -pp)\n", - "# ax.plot(xl, yl, ls='--', alpha=0.5, label=f'{-pp:.4f}')\n", - "\n", - "\n", - "plot._twin_hz(ax, fs=10)\n", - "ax.set(xlabel=LABEL_FREQ_YRS, ylabel=LABEL_STRAIN_AMP)\n", - "\n", - "# ax.legend()\n", - "plt.show()\n", - "\n", - "\n", - "# freqs = np.concatenate([xl, xr])\n", - "# strains = np.concatenate([yl, yr])\n", - "# _, iq = np.unique(freqs, return_index=True)\n", - "# freqs = freqs[iq]\n", - "# strains = strains[iq]\n", - "\n", - "freqs = ul_freq\n", - "strains = ul_strain" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Calculate Luminosity-Distance" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For a circular binary, sky and polarization averaged, the GW strain amplitude is:\n", - "$$h_s = \\frac{8}{10^{1/2}} \\frac{(G \\mathcal{M})^{5/3}}{c^4 d_c} (2\\pi f_r)^{2/3},$$\n", - "where $\\mathcal{M}$ is the rest-frame chirp-mass, $d_c$ is the comoving distance, and $f_r$ is the rest-frame *orbital* frequency. This can also be expressed in terms of observer-frame quantities as,\n", - "$$h_s = \\frac{8}{10^{1/2}} \\frac{(G \\mathcal{M_o})^{5/3}}{c^4 d_L} (2\\pi f_o)^{2/3},$$\n", - "where $d_L$ is the luminosity-distance, and $\\mathcal{M}_o = \\mathcal{M} (1+z)$ and $f_o = f_r/(1+z)$ are the observer-frame chirp-mass and *orbital* frequency respectively." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def _dist_lum_from_strain_amp_and_mchirp_obs(hs, mchirp_obs, fobs_orb):\n", - " fterm = np.power(2*np.pi*fobs_orb, 2.0/3.0)\n", - " mterm = np.power(NWTG*mchirp_obs, 5.0/3.0) / SPLC**4\n", - " dist_lum = (8 / np.sqrt(10)) * mterm * fterm / hs\n", - " return dist_lum\n", - "\n", - "def _mchirp_obs_to_mchirp_rst(dist_lum, mchirp_obs, fobs_orb):\n", - " redz = cosmo.dlum_to_z(dist_lum)\n", - " zp1 = redz + 1.0\n", - " mchirp_rst = mchirp_obs / zp1\n", - " frst_orb = fobs_orb * zp1\n", - " return mchirp_rst, frst_orb\n", - "\n", - "def _strain_from_obs(dist_lum, mchirp_obs, fobs_orb):\n", - " fterm = np.power(2*np.pi*fobs_orb, 2.0/3.0)\n", - " mterm = np.power(NWTG*mchirp_obs, 5.0/3.0) / SPLC**4\n", - " hs = (8 / np.sqrt(10)) * mterm * fterm / dist_lum\n", - " return hs\n", - "\n", - "def dist_from_strain_amp(hs, mchirp_rst, fobs_gw):\n", - " \"\"\"Calculate luminosity distance producing the given GW Strain-Amplitude, from the given binary properties.\n", - " \n", - " Parameters\n", - " ----------\n", - " hs : GW strain-amplitude\n", - " mchirp_rst : rest-frame chirp-mass [gram]\n", - " fobs_gw : observer-frame GW-frequency [1/sec]\n", - " \n", - " Returns\n", - " -------\n", - " dist_lum : luminosity-distance [cm]\n", - " dist_com : comoving-distance [cm]\n", - " redz : redshift\n", - " \n", - " \"\"\"\n", - " # convert from GW to orbital frequency (assuming circular orbit)\n", - " fobs_orb = fobs_gw / 2.0\n", - " # approximate the distance assuming z~0 and thus mc_obs ~ mc_rst\n", - " dist = _dist_lum_from_strain_amp_and_mchirp_obs(hs, mchirp_rst, fobs_orb)\n", - "\n", - " # Function that returns the difference between actual and guessed strain, given a luminosity-distance\n", - " def func(dl, verbose=False):\n", - " # get redshift from luminosity-distance\n", - " zz = cosmo.dlum_to_z(dl)\n", - " # if zz > 10.0:\n", - " # raise ValueError\n", - " # convert from rest-frame chirp-mass to observer frame\n", - " mc_obs = mchirp_rst * (1.0 + zz)\n", - " # calculate GW strain-amp for these values\n", - " hs_test = _strain_from_obs(dl, mc_obs, fobs_orb)\n", - " if verbose:\n", - " print(f\"{dl=} {zz=} {mc_obs=} {hs_test=}\")\n", - "\n", - " # return fractional error\n", - " return (hs - hs_test) / hs\n", - " \n", - " dist_lum = sp.optimize.newton(func, dist, rtol=1e-2, maxiter=int(50))\n", - " # try:\n", - " # dist_lum = sp.optimize.newton(func, dist, rtol=1e-2, maxiter=int(50))\n", - " # except:\n", - " # print(f\"{hs=:.4e}, {mchirp_rst/MSOL=:.4e}, {fobs_gw*YR=:.4e}\")\n", - " # print(f\"{dist/MPC=:.4e} {func(dist, verbose=True)=:.4e}\")\n", - " # dist_lum = sp.optimize.newton(lambda xx: func(xx, True), dist, rtol=1e-2, maxiter=int(1e3))\n", - " # raise\n", - " \n", - " redz = cosmo.dlum_to_z(dist_lum)\n", - " dist_com = cosmo.comoving_distance(redz).cgs.value\n", - " # err = func(dist_lum)\n", - " return dist_lum, dist_com, redz\n", - "\n", - "\n", - "warnings.filterwarnings(\"ignore\", category=RuntimeWarning, message=\"RMS of *\")\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plot.figax()\n", - "mchirp_list = [8.0, 8.5, 9.0, 9.5]\n", - "for mc in mchirp_list[::-1]:\n", - " mc = MSOL * (10.0 ** mc)\n", - " dl_limit, *_ = dist_from_strain_amp(strains, mc, freqs)\n", - " ax.plot(freqs*YR, dl_limit/MPC, label=f'{np.log10(mc/MSOL):.1f}')\n", - "\n", - "plot._twin_hz(ax, fs=10)\n", - "ax.set(xlabel=LABEL_FREQ_YRS, ylabel=LABEL_DIST_LUM_MPC, ylim=[1, 1e3])\n", - "ax.legend(title='Chirp Mass $[\\log_{10}(\\mathcal{M}/\\mathrm{M}_\\odot)]$')\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plot.figax(top=0.78, bottom=0.15)\n", - "mrat_list = [0.1, 0.3, 0.5]\n", - "cmap_list = ['Reds', 'Blues', 'Purples']\n", - "mtot_list = [8.5, 9.0, 9.5, 10.0]\n", - "# lines_list = ['--', '-', '--', '-']\n", - "lines_list = None\n", - "handles_mrat = []\n", - "labels_mrat = []\n", - "handles_mtot = []\n", - "labels_mtot = []\n", - "last_mtot = len(mtot_list) - 1\n", - "for ii, (mrat, cmap) in enumerate(zip(mrat_list, cmap_list)):\n", - " smap = plot.smap(mtot_list, cmap=cmap, left=0.4, right=0.8)\n", - " colors = smap.to_rgba(mtot_list)\n", - " for jj, (_mt, col) in enumerate(zip(mtot_list, colors)):\n", - " mtot = MSOL * (10.0 ** _mt)\n", - " mc = utils.chirp_mass(utils.m1m2_from_mtmr(mtot, mrat))\n", - " \n", - " dl_limit, *_ = dist_from_strain_amp(strains, mc, freqs)\n", - " label = f\"{_mt:.2f}, {mrat:.2f}\"\n", - " ls = '-' if lines_list is None else lines_list[jj]\n", - " hh, = ax.plot(freqs*YR, dl_limit/MPC, color=col, ls=ls)\n", - "\n", - " if ii == 2:\n", - " labels_mtot.append(f\"{_mt:.1f}\")\n", - " handles_mtot.append(hh)\n", - "\n", - " if jj == 1:\n", - " labels_mrat.append(f\"{mrat:.1f}\")\n", - " handles_mrat.append(hh)\n", - "\n", - "plot._twin_hz(ax, fs=10)\n", - "\n", - "ax.set(xlabel=LABEL_FREQ_YRS, ylabel=LABEL_DIST_LUM_MPC, ylim=[1.0, 1e3])\n", - "leg = ax.legend(\n", - " handles_mtot, labels_mtot, bbox_transform=fig.transFigure, bbox_to_anchor=(1.0, 1.0),\n", - " title='Total Mass $[\\log_{10}(M/\\mathrm{M}_\\odot)]$', loc='upper right', ncol=len(labels_mtot)\n", - ")\n", - "ax.legend(\n", - " handles_mrat, labels_mrat, bbox_transform=fig.transFigure, bbox_to_anchor=(0.01, 1.0),\n", - " title='Mass Ratio', loc='upper left', ncol=len(labels_mrat)\n", - ")\n", - "ax.add_artist(leg)\n", - "\n", - "_twin_redz_from_dlum_mpc(ax, zcutoff=0.01)\n", - "plt.show()\n", - "# fig.savefig('test.png')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Calculate Number-Density *of a given binary configuration*" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If we place a limit, such that there are no sources (of given binary parameters) within a comoving distance $d_c$, we can say the local density is less than $n_c = 1/V_c = [(4/3) \\pi d_c^3]^{-1}$. To consider this as a limit on the average density in some volume, that is relatively-local but larger than the explicitly measured volume, there should be some additional pre-factor to account for the confidence of having a source within this volume based on Poisson (or similar) distributions of sources.\n", - "\n", - "To be more precise, define an expected number of events $\\Lambda = n_c V_c$, then the likelihood for $N$ detections in the volume $V_c$ is, $P_N(\\Lambda) = \\frac{\\Lambda^N e^{-\\Lambda}}{N!}$. For no detections, the likelihood is $P_0(\\Lambda) = e^{- \\Lambda}$. To find an upper-limit on the occurrence rate, $\\Lambda_{ul}$, we must integrate from that limit to infinity, such that the result matches our desired confidence level $p_0$ (e.g.~if $p_0=0.95$, then there is a $95\\%$ chance that the rate is below this upper-limit): $F_{ul}(\\Lambda_{ul}) = \\int_{\\Lambda_{ul}}^\\infty e^{-\\Lambda} d\\Lambda = 1 - p_0$. Or, converting this back to the dimensional volume density, we find,\n", - "$$n_{ul} = \\frac{-\\ln (1-p_0)}{V_c}.$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Plot $95\\%$ upper-limits on number density of sources at some particular chirp-masses." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def ndens_from_strain_amp(hs, mchirp_rst, fobs_gw, conf=0.95):\n", - " assert 0.0 < conf < 1.0\n", - " try:\n", - " dist_lum, dist_com, redz = dist_from_strain_amp(hs, mchirp_rst, fobs_gw)\n", - " except Exception as err:\n", - " # print(dist_lum)\n", - " log.error(err)\n", - " dist_lum = 1e5 * MPC\n", - " redz = cosmo.dlum_to_z(dist_lum)\n", - " dist_com = cosmo.comoving_distance(redz).cgs.value\n", - " # print(dist_lum, redz, dist_com)\n", - " \n", - " vol_com = (4.0/3.0) * np.pi * (dist_com/MPC) ** 3\n", - " conf_factor = - np.log(1.0 - conf)\n", - " ndens_com_mpc3 = conf_factor / vol_com\n", - " return ndens_com_mpc3\n", - " \n", - "fig, ax = plot.figax()\n", - "mchirp_list = [8.0, 8.5, 9.0, 9.5]\n", - "ndens_limits = []\n", - "for mc in mchirp_list:\n", - " mc = MSOL * (10.0 ** mc)\n", - " ndlim = ndens_from_strain_amp(strains, mc, freqs)\n", - " ndens_limits.append(ndlim)\n", - " ax.plot(freqs*YR, ndlim, label=f'{np.log10(mc/MSOL):.1f}')\n", - "\n", - "plot._twin_hz(ax, fs=10)\n", - "ax.set(xlabel=LABEL_FREQ_YRS, ylabel=LABEL_NDENS_MPC3, ylim=[1e-9, 1e3])\n", - "ax.legend(title='Chirp Mass $[\\log_{10}(\\mathcal{M}/\\mathrm{M}_\\odot)]$')\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Plot $95\\%$ upper-limits on number density of sources at some particular total-masses and mass-ratios." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plot.figax(top=0.78, bottom=0.15, left=0.08, right=0.98)\n", - "mrat_list = [0.1, 0.3, 0.5]\n", - "cmap_list = ['Reds', 'Blues', 'Purples']\n", - "mtot_list = [8.5, 9.0, 9.5, 10.0]\n", - "# lines_list = ['--', '-', '--', '-']\n", - "lines_list = None\n", - "handles_mrat = []\n", - "labels_mrat = []\n", - "handles_mtot = []\n", - "labels_mtot = []\n", - "last_mtot = len(mtot_list) - 1\n", - "for ii, (mrat, cmap) in enumerate(zip(mrat_list, cmap_list)):\n", - " smap = plot.smap(mtot_list, cmap=cmap, left=0.4, right=0.8)\n", - " colors = smap.to_rgba(mtot_list)\n", - " for jj, (_mt, col) in enumerate(zip(mtot_list, colors)):\n", - " mtot = MSOL * (10.0 ** _mt)\n", - " mc = utils.chirp_mass(utils.m1m2_from_mtmr(mtot, mrat))\n", - " \n", - " ndens_limit = ndens_from_strain_amp(strains, mc, freqs)\n", - " label = f\"{_mt:.2f}, {mrat:.2f}\"\n", - " ls = '-' if lines_list is None else lines_list[jj]\n", - " hh, = ax.plot(freqs*YR, ndens_limit, color=col, ls=ls)\n", - "\n", - " if ii == 2:\n", - " labels_mtot.append(f\"{_mt:.1f}\")\n", - " handles_mtot.append(hh)\n", - "\n", - " if jj == 1:\n", - " labels_mrat.append(f\"{mrat:.1f}\")\n", - " handles_mrat.append(hh)\n", - "\n", - "plot._twin_hz(ax, fs=10)\n", - "\n", - "ax.set(xlabel=LABEL_FREQ_YRS, ylabel=LABEL_NDENS_MPC3, ylim=[1.0e-9, 1e2])\n", - "leg = ax.legend(\n", - " handles_mtot, labels_mtot, bbox_transform=fig.transFigure, bbox_to_anchor=(1.0, 1.0),\n", - " title='Total Mass $[\\log_{10}(M/\\mathrm{M}_\\odot)]$', loc='upper right', ncol=len(labels_mtot)\n", - ")\n", - "ax.legend(\n", - " handles_mrat, labels_mrat, bbox_transform=fig.transFigure, bbox_to_anchor=(0.01, 1.0),\n", - " title='Mass Ratio', loc='upper left', ncol=len(labels_mrat)\n", - ")\n", - "ax.add_artist(leg)\n", - "\n", - "plt.show()\n", - "fig.savefig('ndens-constraints_witt.png')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Calculate Number-Density of *all binaries, based on some population model*" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_fobs_orb = freqs / 2.0 # convert from GW to orbital frequency (assume circular orbits)\n", - "df = np.diff(freqs)\n", - "freqs_edges = freqs[1] - freqs[0]\n", - "freqs_edges = [freqs[0] - df[0]/2.0,] + utils.midpoints(freqs).tolist() + [freqs[-1] + df[-1]/2.0,]\n", - "freqs_edges = np.array(freqs_edges)\n", - "fobs_orb_edges = freqs_edges / 2.0\n", - "\n", - "pop_sam = holo.sam.Semi_Analytic_Model()\n", - "# time = 5 * GYR\n", - "# hard = holo.evolution.Fixed_Time.from_sam(pop_sam, time)\n", - "hard = holo.hardening.Hard_GW()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# `edges` is (4,) of array giving the bin-edges of total-mass, mass-ratio, redshift, freq-obs-orbital\n", - "# `d4_num_sam` is shaped (M, Q, Z, F) where each dimension length is the number of bin-edges minus 1\n", - "# this is `d^4 N / [dlog10(M) dq dz dln(f_r)]`\n", - "edges, d4_num_sam = pop_sam.dynamic_binary_number(hard, fobs_orb=fobs_orb_edges, zero_stalled=False)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Integrate over each bin to find the number of binaries in each bin. Units of number (dimensionless).\n", - "sam_number = holo.utils._integrate_grid_differential_number(edges, d4_num_sam, freq=True)\n", - "# Find the cumulative sum over redshifts, meaning the number of binaries out to that redshift\n", - "sam_number = np.cumsum(sam_number, axis=2)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Convert back to volume-density of binaries in each bin, i.e. [Mpc^-3]\n", - "dvol = np.diff(cosmo.comoving_volume(edges[2])).to('Mpc3').value\n", - "sam_ndens = sam_number / dvol[np.newaxis, np.newaxis, :, np.newaxis]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cents = [utils.midpoints(xx) for xx in [pop_sam.mtot, pop_sam.mrat]]\n", - "mc = utils.chirp_mass_mtmr(cents[0][:, np.newaxis], cents[1][np.newaxis, :])\n", - "mc = mc.flatten()\n", - "sam_ndens_chirp = np.copy(sam_ndens).reshape(-1, *sam_ndens.shape[2:])\n", - "idx = np.argsort(mc)\n", - "mc = mc[idx]\n", - "sam_ndens_chirp = sam_ndens_chirp[idx][::-1]\n", - "print(sam_ndens_chirp.shape)\n", - "sam_ndens_chirp = np.cumsum(sam_ndens_chirp, axis=0)[::-1]\n", - "print(sam_ndens_chirp.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# (Mc, z, f)\n", - "# sam_ndens_chirp\n", - "fig, ax = plot.figax()\n", - "cmap = mpl.cm.get_cmap('viridis')\n", - "\n", - "# find the most-sensitive 12.5yr frequency bin\n", - "fbin = np.argmin(ul_strain)\n", - "vals = sam_ndens_chirp[:, :, fbin]\n", - "num_mchirp = vals.shape[0]\n", - "colors = cmap(np.linspace(0.0, 1.0, num_mchirp))\n", - "log10_mc = np.log10(mc/MSOL)\n", - "\n", - "dlum = cosmo.z_to_dlum(pop_sam.redz[1:])\n", - "xx = dlum / MPC\n", - "\n", - "# Go through target chirp-masses\n", - "for mchirp in mchirp_list:\n", - " # Find the population index for this chirp-mass\n", - " idx = np.argmax(log10_mc > mchirp)\n", - " cc = colors[idx]\n", - " # plot the average number-density out to each redshift, for binaries at this chirp-mass and above\n", - " ax.plot(xx, vals[idx], color=cc, label=mchirp)\n", - "\n", - " # plot the number-density limit at this chirp-mass\n", - " _mc = MSOL * (10.0**mchirp)\n", - " ndlim = ndens_from_strain_amp(strains[fbin], _mc, freqs[fbin], conf=0.95)\n", - " ax.axhline(ndlim, color=cc, ls='--', alpha=0.5)\n", - " \n", - " dl_limit, *_ = dist_from_strain_amp(strains[fbin], _mc, freqs[fbin])\n", - " # zlim = cosmo.dlum_to_z(dl_limit)\n", - " # print(_mc/MSOL, dl_limit/MPC, zlim)\n", - " ax.axvline(dl_limit/MPC, color=cc, ls='--', alpha=0.5)\n", - " \n", - " \n", - "ax.legend()\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plot.figax()\n", - "# tw = ax.twinx()\n", - "# tw.set(yscale='log')\n", - "cmap = mpl.cm.get_cmap('viridis')\n", - "vals = sam_ndens_chirp\n", - "num = vals.shape[0]\n", - "colors = cmap(np.linspace(0.05, 0.85, len(mchirp_list)))\n", - "xx = np.log10(mc/MSOL)\n", - "for ii, mchirp in enumerate(mchirp_list):\n", - " cc = colors[ii]\n", - " idx = np.argmax(xx > mchirp)\n", - " # Select this chirp-mass index, (Z, F)\n", - " vals = sam_ndens_chirp[idx, :, :]\n", - " # Get the maximum density over redshifts\n", - " vals = np.max(vals, axis=0)\n", - "\n", - " ax.plot(freqs*YR, vals, color='0.75', ls='-', lw=3.0, alpha=0.25, zorder=3)\n", - " ax.plot(freqs*YR, vals, color=cc, ls='--')\n", - "\n", - " lim = ndens_limits[ii]\n", - " ax.plot(freqs*YR, lim, color=cc, ls='-', label=f'$10^{{{mchirp:.1f}}} \\, M_\\odot$')\n", - " \n", - " # tw.plot(freqs*YR, vals/lim, color=cc, ls=':')\n", - "\n", - "ax.set(xlabel='frequency [yr$^{-1}$]', ylabel='number density [Mpc$^{-3}$]', ylim=[1e-15, 1e5])\n", - "ax.legend(loc='lower right', title='chirp mass')\n", - "fname = \"12p5yr_cw__ndens_limit_vs_predicted.png\"\n", - "fig.savefig(fname, dpi=300)\n", - "print(f\"saved to {fname}\")\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mrat_list = [0.05, 0.1, 0.3, 0.5]\n", - "\n", - "fig, ax = plot.figax(\n", - " top=0.97, bottom=0.16, left=0.08, right=0.92,\n", - " xlabel=LABEL_MASS_TOTAL, ylabel=LABEL_NDENS_MPC3\n", - ")\n", - "tw = ax.twinx()\n", - "tw.set(yscale='log', ylabel='Ratio (upper-limit / model | dotted)')\n", - "ax.set(ylim=[1e-12, 1e7])\n", - "tw.set(ylim=[1e6, 1e10])\n", - "\n", - "mtot_centers = kale.utils.midpoints(edges[0], log=True)\n", - "xx = mtot_centers / MSOL\n", - "\n", - "for mrat in mrat_list:\n", - " qidx = np.searchsorted(edges[1], mrat)\n", - " mrat = edges[1][qidx]\n", - " \n", - " # sum over all larger mass-ratios\n", - " ndens_array = np.sum(ndens[:, qidx:], axis=1)\n", - " # cumulative sum over all larger masses, for each total-mass bin\n", - " ndens_array = np.cumsum(ndens_array[::-1])[::-1]\n", - "\n", - " ndens_limit = np.zeros_like(ndens_array)\n", - " for ii, mt in enumerate(mtot_centers):\n", - " mc = holo.utils.chirp_mass(holo.utils.m1m2_from_mtmr(mt, mrat))\n", - " # print(ii, f\"{mc/MSOL:.2e}\", fbin, freqs[fbin], strains[fbin])\n", - " ndens_limit[ii] = ndens_from_strain_amp(strains[fbin], mc, freqs[fbin], conf=0.95)\n", - "\n", - " hh, = ax.plot(xx, ndens_array, label=f'{mrat:.2f}')\n", - " cc = hh.get_color()\n", - " ax.plot(xx, ndens_limit, color=cc, ls='--')\n", - "\n", - " tw.plot(xx, ndens_limit/ndens_array, color=cc, ls=':')\n", - "\n", - "ax.legend(title='Mass Ratio')\n", - "plt.show()\n", - "fig.savefig('population-limit_witt.png')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sam_ndens.shape" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fbin = np.argmin(ul_strain)\n", - "print(f\"{fbin=}\")\n", - "ndens = sam_ndens[..., fbin]\n", - "imax = np.argmax(ndens, axis=-1)\n", - "assert np.all(imax == imax.flatten()[0])\n", - "imax = imax.flatten()[0]\n", - "ndens = ndens[..., imax]\n", - "\n", - "mrat_list = [0.05, 0.1, 0.3, 0.5]\n", - "\n", - "fig, ax = plot.figax(\n", - " top=0.97, bottom=0.16, left=0.08, right=0.92,\n", - " xlabel=LABEL_MASS_TOTAL, ylabel=LABEL_NDENS_MPC3\n", - ")\n", - "tw = ax.twinx()\n", - "tw.set(yscale='log', ylabel='Ratio (upper-limit / model | dotted)')\n", - "ax.set(ylim=[1e-12, 1e7])\n", - "tw.set(ylim=[1e6, 1e10])\n", - "\n", - "mtot_centers = kale.utils.midpoints(edges[0], log=True)\n", - "xx = mtot_centers / MSOL\n", - "\n", - "for mrat in mrat_list:\n", - " qidx = np.searchsorted(edges[1], mrat)\n", - " mrat = edges[1][qidx]\n", - " \n", - " # sum over all larger mass-ratios\n", - " ndens_array = np.sum(ndens[:, qidx:], axis=1)\n", - " # cumulative sum over all larger masses, for each total-mass bin\n", - " ndens_array = np.cumsum(ndens_array[::-1])[::-1]\n", - "\n", - " ndens_limit = np.zeros_like(ndens_array)\n", - " for ii, mt in enumerate(mtot_centers):\n", - " mc = holo.utils.chirp_mass(holo.utils.m1m2_from_mtmr(mt, mrat))\n", - " # print(ii, f\"{mc/MSOL:.2e}\", fbin, freqs[fbin], strains[fbin])\n", - " ndens_limit[ii] = ndens_from_strain_amp(strains[fbin], mc, freqs[fbin], conf=0.95)\n", - "\n", - " hh, = ax.plot(xx, ndens_array, label=f'{mrat:.2f}')\n", - " cc = hh.get_color()\n", - " ax.plot(xx, ndens_limit, color=cc, ls='--')\n", - "\n", - " tw.plot(xx, ndens_limit/ndens_array, color=cc, ls=':')\n", - "\n", - "ax.legend(title='Mass Ratio')\n", - "plt.show()\n", - "fig.savefig('population-limit_witt.png')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - " " - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3.9.13 ('py39')", - "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.10.8" - }, - "vscode": { - "interpreter": { - "hash": "61b4062b24dfb1010f420dad5aa3bd73a4d2af47d0ec44eafec465a35a9d7239" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/analysis/pop-data.ipynb b/analysis/pop-data.ipynb deleted file mode 100644 index cb9b98be..00000000 --- a/analysis/pop-data.ipynb +++ /dev/null @@ -1,173 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "import holodeck as holo\n", - "import holodeck.plot\n", - "import holodeck.utils\n", - "from holodeck.constants import YR" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# SAM Populations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fname_real = \"./pop-sam_v0.2_b54efb6__d20.0_c0.100_sh0050_st+2.000/0000.npz\"\n", - "fname_sam = \"./pop-sam_v0.2_b54efb6__d20.0_c0.100_sh0050_st+2.000/sam.npz\"\n", - "\n", - "data_sam = np.load(fname_sam)\n", - "fobs = data_sam['fobs']\n", - "print(f\"SAM frequencies: [{fobs[0]*YR:.2f}, {fobs[-1]*YR:.2f}] 1/yr, {fobs.size-1} bins\")\n", - "\n", - "data_real = np.load(fname_real)\n", - "weights = data_real['weights']\n", - "print(f\"Sample points: {weights.size:.4e}, effective number of binaries: {weights.sum():.4e}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nloud = 5\n", - "\n", - "hs, fo = holo.sam._strains_from_samples([data_real[kk] for kk in ['mtot', 'mrat', 'redz', 'fobs']])\n", - "colors = holo.plot._get_cmap('plasma')(np.linspace(0.1, 0.9, nloud))\n", - "\n", - "fig, ax = holo.plot.figax(figsize=[12, 8], xlabel='Frequency [yr$^{-1}$]', ylabel='c-Strain')\n", - "for ii in holo.utils.tqdm(range(fobs.size-1)):\n", - " fextr = [fobs[ii+jj] for jj in range(2)]\n", - " fextr = np.asarray(fextr)\n", - " cycles = 1.0 / np.diff(np.log(fextr))[0]\n", - "\n", - " idx = (fextr[0] <= fo) & (fo < fextr[1])\n", - " hs_bin = hs[idx]\n", - " fo_bin = fo[idx] \n", - " ww_bin = weights[idx]\n", - " ww = ww_bin * cycles\n", - "\n", - " tot = np.sqrt(np.sum(ww * hs_bin**2))\n", - " ax.plot(fextr*YR, tot * np.ones_like(fextr), 'k--')\n", - "\n", - " idx = np.argsort(hs_bin)[::-1]\n", - " if any(ww_bin[idx[:nloud]] > 1):\n", - " raise\n", - " \n", - " for jj, cc in enumerate(colors):\n", - " if jj > len(idx):\n", - " break\n", - " hi = idx[jj]\n", - " lo = idx[jj+1:]\n", - " gw_hi = np.sqrt(np.sum(ww[hi] * hs_bin[hi]**2))\n", - " gw_lo = np.sqrt(np.sum(ww[lo] * hs_bin[lo]**2))\n", - "\n", - " fave = np.average(fo_bin[hi], weights=hs_bin[hi])\n", - " ax.plot(fextr*YR, gw_lo * np.ones_like(fextr), color=cc, lw=0.5)\n", - " ax.scatter(fave*YR, gw_hi, marker='.', color=cc, alpha=0.5)\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Illustris Population" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fname = \"./pop-ill_v0.2_e6bc2cb__d20.0_c0.100_mamp+8.845098_mamp1.00/debug_0000.npz\"\n", - "data = np.load(fname)\n", - "print(list(data.keys()))\n", - "print(f\"Loaded {data['redz'].size:.4e} samples\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nloud = 5\n", - "fobs = data['freq_bins']\n", - "hs, fo = holo.sam._strains_from_samples([data[kk] for kk in ['mtot', 'mrat', 'redz', 'fobs']])\n", - "colors = holo.plot._get_cmap('plasma')(np.linspace(0.1, 0.9, nloud))\n", - "\n", - "fig, ax = holo.plot.figax(figsize=[12, 8], xlabel='Frequency [yr$^{-1}$]', ylabel='c-Strain')\n", - "for ii in holo.utils.tqdm(range(fobs.size-1)):\n", - " fextr = [fobs[ii+jj] for jj in range(2)]\n", - " fextr = np.asarray(fextr)\n", - " cycles = 1.0 / np.diff(np.log(fextr))[0]\n", - "\n", - " idx = (fextr[0] <= fo) & (fo < fextr[1])\n", - " hs_bin = hs[idx]\n", - " fo_bin = fo[idx] \n", - "\n", - " tot = np.sqrt(np.sum(cycles * hs_bin**2))\n", - " ax.plot(fextr*YR, tot * np.ones_like(fextr), 'k--')\n", - "\n", - " idx = np.argsort(hs_bin)[::-1]\n", - " \n", - " for jj, cc in enumerate(colors):\n", - " if jj >= len(idx):\n", - " break\n", - " hi = idx[jj]\n", - " lo = idx[jj+1:]\n", - " gw_hi = np.sqrt(np.sum(cycles * hs_bin[hi]**2))\n", - " gw_lo = np.sqrt(np.sum(cycles * hs_bin[lo]**2))\n", - "\n", - " fave = np.average(fo_bin[hi], weights=hs_bin[hi])\n", - " ax.plot(fextr*YR, gw_lo * np.ones_like(fextr), color=cc, lw=0.5)\n", - " ax.scatter(fave*YR, gw_hi, marker='.', color=cc, alpha=0.5)\n", - "\n", - "plt.show()" - ] - } - ], - "metadata": { - "interpreter": { - "hash": "61b4062b24dfb1010f420dad5aa3bd73a4d2af47d0ec44eafec465a35a9d7239" - }, - "kernelspec": { - "display_name": "Python 3.9.12 ('py39')", - "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.9.12" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/analysis/gen-15yr-pops/gen_holodeck_pops.py b/holodeck/librarian/posterior_populations.py similarity index 100% rename from analysis/gen-15yr-pops/gen_holodeck_pops.py rename to holodeck/librarian/posterior_populations.py diff --git a/analysis/gen-15yr-pops/holodeck-pops.ipynb b/notebooks/nanograv-15yr-populations.ipynb similarity index 100% rename from analysis/gen-15yr-pops/holodeck-pops.ipynb rename to notebooks/nanograv-15yr-populations.ipynb