diff --git a/holodeck/librarian/param_spaces.py b/holodeck/librarian/param_spaces.py index 600b7a0b..fe34b5b0 100644 --- a/holodeck/librarian/param_spaces.py +++ b/holodeck/librarian/param_spaces.py @@ -158,7 +158,7 @@ class _PS_Astro_Strong(_Param_Space): """ - __version__ = "0.1" + __version__ = "0.2" DEFAULTS = dict( # Hardening model (phenom 2PL) @@ -199,6 +199,9 @@ class _PS_Astro_Strong(_Param_Space): mmb_mamp=0.49e9, # 0.49e9 + 0.06 - 0.05 [Msol] mmb_plaw=1.17, # 1.17 ± 0.08 mmb_scatter_dex=0.28, # no uncertainties given + # bulge fraction + bf_sigmoid_lo=0.4, + bf_sigmoid_hi=0.8, ) @classmethod diff --git a/holodeck/sams/sam.py b/holodeck/sams/sam.py index 362c37b0..3f691459 100644 --- a/holodeck/sams/sam.py +++ b/holodeck/sams/sam.py @@ -466,12 +466,18 @@ def dynamic_binary_number_at_fobs(self, hard, fobs_orb, use_cython=True, **kwarg return grid, dnum, redz_final def _dynamic_binary_number_at_fobs_consistent(self, hard, fobs_orb, steps=200, details=False): - """Get correct redshifts for full binary-number calculation. + r"""Calculate the differential number of binaries in at each grid point, at each frequency. - Slower but more correct than old `dynamic_binary_number`. - Same as new cython implementation `sam_cyutils.dynamic_binary_number_at_fobs`, which is - more than 10x faster. - LZK 2023-05-11 + See :meth:`dynamic_binary_number_at_fobs` for general information. + + This is the python implementation for binary evolution (hardening) that is self-consistent, + i.e. evolution models that are able to evolve binaries from galaxy merger until the target + frequencies. + + This function should produce the same results as the new cython implementation in: + :func:`holodeck.sams.sam_cyutils.dynamic_binary_number_at_fobs`, which is more than 10x + faster. This python implementation is maintained for diagnostic purposes, and for + functionality when cython is not available. # BUG doesn't work for Fixed_Time_2PL @@ -480,19 +486,27 @@ def _dynamic_binary_number_at_fobs_consistent(self, hard, fobs_orb, steps=200, d edges = self.edges + [fobs_orb, ] # shape: (M, Q, Z) - dens = self.static_binary_density # d3n/[dlog10(M) dq dz] units: [Mpc^-3] + dens = self.static_binary_density # d3n/[dlog10(M) dq dz] units: [cMpc^-3] + + # ---- Choose the binary separations over which to integrate the binary evolution. + + # Start at large separations (galaxy merger) and evolve to small separations (coalescense). # start from the hardening model's initial separation rmax = hard._sepa_init - # (M,) end at the ISCO + # end at the ISCO + # (M,) rmin = utils.rad_isco(self.mtot) # Choose steps for each binary, log-spaced between rmin and rmax extr = np.log10([rmax * np.ones_like(rmin), rmin]) # (2,M,) - rads = np.linspace(0.0, 1.0, steps+1)[np.newaxis, :] # (1,X) - # (M, S) = (M,1) * (1,S) + rads = np.linspace(0.0, 1.0, steps+1)[np.newaxis, :] # (1,S) + # (M, S) <== (M,1) * (1,S) rads = extr[0][:, np.newaxis] + (extr[1] - extr[0])[:, np.newaxis] * rads rads = 10.0 ** rads + # ---- Calculate binary hardening rate (da/dt) at each separation, for each grid point + + # broadcast arrays to a consistent shape # (M, Q, S) mt, mr, rads, norm = np.broadcast_arrays( self.mtot[:, np.newaxis, np.newaxis], @@ -500,33 +514,49 @@ def _dynamic_binary_number_at_fobs_consistent(self, hard, fobs_orb, steps=200, d rads[:, np.newaxis, :], hard._norm[:, :, np.newaxis], ) + # calculate hardening rate (negative values, in units of [cm/s]) dadt_evo = hard.dadt(mt, mr, rads, norm=norm) + # ---- Integrate evolution + # to find times and redshifts at which binaries reach each separation + # (M, Q, S-1) # Integrate (inverse) hardening rates to calculate total lifetime to each separation + times_evo = -utils.trapz_loglog(-1.0 / dadt_evo, rads, axis=-1, cumsum=True) + # ~~~~ RIEMANN integration ~~~~ + # times_evo = 2.0 * np.diff(rads, axis=-1) / (dadt_evo[..., 1:] + dadt_evo[..., :-1]) + # times_evo = np.cumsum(times_evo, axis=-1) + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # times_evo = -utils.trapz_loglog(-1.0 / dadt_evo, rads, axis=-1, cumsum=True) - - times_evo = 2.0 * np.diff(rads, axis=-1) / (dadt_evo[..., 1:] + dadt_evo[..., :-1]) - # for ss in range(steps): - # print(f"py {ss:03d} : {rads[8, 0, ss]:.6e} ==> {rads[8, 0, ss+1]:.6e} == {times_evo[8, 0, ss]:.6e}") - times_evo = np.cumsum(times_evo, axis=-1) # add array of zero time-delays at starting point (i.e. before the first step) # with same shape as a slice at a single step zpad = np.zeros_like(times_evo[..., 0]) times_evo = np.concatenate([zpad[..., np.newaxis], times_evo], axis=-1) - # print(f"{times_evo[8, 0, :]=}") - # Combine the binary-evolution time, with the galaxy-merger time - # (M, Q, Z, S-1) - rz = self.redz[np.newaxis, np.newaxis, :, np.newaxis] - times_tot = times_evo[:, :, np.newaxis, :] + # ---- Convert from time to redshift + + # initial redshift (of galaxy merger) + rz = self.redz[np.newaxis, np.newaxis, :, np.newaxis] # (1, 1, Z, 1) + + tlbk_init = cosmo.z_to_tlbk(rz) + tlbk = tlbk_init - times_evo[:, :, np.newaxis, :] + # Combine the binary-evolution time, with the galaxy-merger time (if it is defined) if self._gmt_time is not None: - times_tot += self._gmt_time[:, :, :, np.newaxis] + tlbk -= self._gmt_time[:, :, :, np.newaxis] - redz_evo = utils.redz_after(times_tot, redz=rz) - # for ss in range(steps): - # print(f"py {ss:03d} : t={times_evo[8, 0, ss]:.6e} z={redz_evo[8, 0, 11, ss]:.6e}") + # (M, Q, Z, S) + redz_evo = cosmo.tlbk_to_z(tlbk) + + #! age of the universe version of calculation is MUCH less accurate !# + # Use age-of-the-universe + # times_tot = times_evo[:, :, np.newaxis, :] + # # Combine the binary-evolution time, with the galaxy-merger time (if it is defined) + # if self._gmt_time is not None: + # times_tot += self._gmt_time[:, :, :, np.newaxis] + # redz_evo = utils.redz_after(times_tot, redz=rz) + #! ---------------------------------------------------------------- !# + + # ---- interpolate to target frequencies # convert from separations to rest-frame orbital frequencies # (M, Q, S) @@ -534,25 +564,14 @@ def _dynamic_binary_number_at_fobs_consistent(self, hard, fobs_orb, steps=200, d # (M, Q, Z, S) fobs_orb_evo = frst_orb_evo[:, :, np.newaxis, :] / (1.0 + redz_evo) - # ---- interpolate to target frequencies - # `ndinterp` interpolates over 1th dimension - - # print(f"{frst_orb_evo[8, 0, :]=}") - # print(f"{fobs_orb=}") - # print(f"{fobs_orb_evo[8, 0, 11, :]=}") - # print(f"{redz_evo[8, 0, 11, :]=}") - - # (M, Q, Z, S-1) ==> (M*Q*Z, S-1) + # (M, Q, Z, S) ==> (M*Q*Z, S) fobs_orb_evo, redz_evo = [tt.reshape(-1, steps+1) for tt in [fobs_orb_evo, redz_evo]] + # `ndinterp` interpolates over 1th dimension # (M*Q*Z, X) - redz_final = utils.ndinterp(fobs_orb, fobs_orb_evo, redz_evo, xlog=False, ylog=False) + redz_final = utils.ndinterp(fobs_orb, fobs_orb_evo, redz_evo, xlog=True, ylog=False) - # (M*Q*Z, X) ===> (M, Q, Z, X) + # (M, Q, Z, X) <=== (M*Q*Z, X) redz_final = redz_final.reshape(self.shape + (fobs_orb.size,)) - # _test_times = _test_times.reshape(self.shape + (fobs_orb.size,)) - - # print(f"{redz_final[8, 0, 11, :]=}") - # print(f"{_test_times[8, 0, 11, :]=}") coal = (redz_final > 0.0) frst_orb = fobs_orb * (1.0 + redz_final) diff --git a/holodeck/sams/sam_cyutils.pyx b/holodeck/sams/sam_cyutils.pyx index 8981075d..ed411270 100644 --- a/holodeck/sams/sam_cyutils.pyx +++ b/holodeck/sams/sam_cyutils.pyx @@ -144,9 +144,10 @@ def integrate_differential_number_3dx1d(edges, dnum): # each edge should have the same length as the corresponding dimension of `dnum` shape = [len(ee) for ee in edges] + err = f"Shape of edges={shape} does not match dnum={np.shape(dnum)}" # except the last edge (freq), where `dnum` should be 1-shorter shape[-1] -= 1 - assert np.shape(dnum) == tuple(shape) + assert np.shape(dnum) == tuple(shape), err # the number will be shaped as one-less the size of each dimension of `dnum` new_shape = [sh-1 for sh in dnum.shape] # except for the last dimension (freq) which is the same shape @@ -643,7 +644,11 @@ cdef int _dynamic_binary_number_at_fobs_2pwl( ) # Find time to move from left- to right- edges: dt = da / (da/dt) + # average da/dt on the left- and right- edges of the bin (i.e. trapezoid rule) dt = 2.0 * (sepa_right - sepa_left) / (dadt_left + dadt_right) + # if ii == 8 and jj == 0: + # printf("cy %03d : %.2e ==> %.2e == %.2e\n", step, sepa_left, sepa_right, dt) + time_evo += dt # ---- Iterate over starting redshift bins @@ -657,8 +662,8 @@ cdef int _dynamic_binary_number_at_fobs_2pwl( # if we pass the age of the universe, this binary has stalled, no further redshifts will work # NOTE: if `gmt_time` decreases faster than redshift bins increase the universe age, - # then systems in later `redz` bins may no longer stall, so we still need to calculate them - # i.e. we can NOT use a `break` statement here + # then systems in later `redz` bins may no longer stall, so we still need to calculate them. + # i.e. we can NOT use a `break` statement here, must use `continue` statement. if time_left > age_universe: continue @@ -684,6 +689,9 @@ cdef int _dynamic_binary_number_at_fobs_2pwl( if redz_right < 0.0: redz_right = 0.0 + # if ii == 8 and jj == 0 and kk == 11: + # printf("cy %03d : t=%.2e z=%.2e\n", step, time_right, redz_right) + # convert to frequencies fobs_orb_left = frst_orb_left / (1.0 + redz_left) fobs_orb_right = frst_orb_right / (1.0 + redz_right) @@ -724,6 +732,21 @@ cdef int _dynamic_binary_number_at_fobs_2pwl( # get comoving distance dcom = interp_at_index(new_interp_idx, new_time, tage_interp_grid, dcom_interp_grid) + # if (ii == 0) and (jj == 0) and (kk == 0): + # printf("cy f=%03d (step=%03d)\n", ff, step) + # printf( + # "fl=%.6e, f=%.6e, fr=%.6e ==> tl=%.6e, t=%.6e, tr=%.6e\n", + # fobs_orb_left, ftarget, fobs_orb_right, + # time_left, new_time, time_right + # ) + # printf( + # "interp (%d) time: %.6e, %.6e, %.6e ==> z: %.6e, %.6e, %.6e\n", + # new_interp_idx, + # tage_interp_grid[new_interp_idx], new_time, tage_interp_grid[new_interp_idx+1], + # redz_interp_grid[new_interp_idx], new_redz, redz_interp_grid[new_interp_idx+1], + # ) + # printf("======> z=%.6e\n", new_redz) + # Store redshift redz_final[ii, jj, kk, ff] = new_redz diff --git a/holodeck/sams/tests/test_sam.py b/holodeck/sams/tests/test_sam.py index 36beb43c..1a1874a6 100644 --- a/holodeck/sams/tests/test_sam.py +++ b/holodeck/sams/tests/test_sam.py @@ -128,9 +128,8 @@ def test_dbn_gw_only(): NUM_FREQS = 9 fobs_gw_cents, fobs_gw_edges = holo.utils.pta_freqs(PTA_DUR, NUM_FREQS) fobs_orb_cents = fobs_gw_cents / 2.0 - # fobs_orb_edges = fobs_gw_edges / 2.0 - # (1) + # (1) make sure it runs grid_py, dnum_py, redz_final_py = sam.dynamic_binary_number_at_fobs(hard_gw, fobs_orb_cents, use_cython=False) grid_cy, dnum_cy, redz_final_cy = sam.dynamic_binary_number_at_fobs(hard_gw, fobs_orb_cents, use_cython=True) @@ -157,6 +156,96 @@ def test_dbn_gw_only(): print(f"{utils.stats(redz_final_cy[bads])=}") assert not np.any(bads), f"Found {utils.frac_str(bads)} inconsistent `redz_final` b/t python and cython calcs!" + return + + +def test_dbn_phenom(): + """Test the dynamic_binary_number method using Phenomenological evolution. + + (1) runs without error + (2) dnum values are consistent between cython and python + (3) redz_final values are consistent between cython and python + + """ + + shape = (10, 11, 12) + sam = holo.sams.Semi_Analytic_Model(shape=shape) + TIME = 1.0e9 * YR + hard_phenom = holo.hardening.Fixed_Time_2PL_SAM(sam, TIME, num_steps=300) + + PTA_DUR = 20.0 * YR + NUM_FREQS = 9 + fobs_gw_cents, fobs_gw_edges = holo.utils.pta_freqs(PTA_DUR, NUM_FREQS) + fobs_orb_cents = fobs_gw_cents / 2.0 + fobs_orb_edges = fobs_gw_edges / 2.0 + + # we'll allow differences at very low redshifts, where numerical differences become significant + ALLOW_BADS_BELOW_REDZ = 1.0e-2 + + # (1) make sure it runs + + grid_py, dnum_py, redz_final_py = sam.dynamic_binary_number_at_fobs(hard_phenom, fobs_orb_cents, use_cython=False) + grid_cy, dnum_cy, redz_final_cy = sam.dynamic_binary_number_at_fobs(hard_phenom, fobs_orb_cents, use_cython=True) + edges_py = grid_py[:-1] + [fobs_orb_edges,] + edges_cy = grid_cy[:-1] + [fobs_orb_edges,] + + redz_not_ignore = (redz_final_py > ALLOW_BADS_BELOW_REDZ) | (redz_final_cy > ALLOW_BADS_BELOW_REDZ) + + # (2) the same dnum values are zero + + zeros_py = (dnum_py == 0.0) + zeros_cy = (dnum_cy == 0.0) + + # ignore mismastch at low-redshifts + bads = (zeros_py != zeros_cy) & redz_not_ignore + if np.any(bads): + print(f"{utils.frac_str(bads)=}") + print(f"{utils.stats(dnum_py[bads])=}") + print(f"{utils.stats(dnum_cy[bads])=}") + assert not np.any(bads), "Zero points in `dnum` do not match between python and cython!" + + # (3) dnum consistent between cython- and python- versions of calculation + + # ignore mismastch at low-redshifts + bads = ~np.isclose(dnum_py, dnum_cy, rtol=1e-1) & redz_not_ignore + if np.any(bads): + errs = (dnum_py - dnum_cy) / dnum_cy + print(f"{utils.frac_str(bads)=}") + print(f"{utils.stats(errs)=}") + print(f"{utils.stats(errs[bads])=}") + print(f"{dnum_py[bads][:10]=}") + print(f"{dnum_cy[bads][:10]=}") + print(f"{errs[bads][:10]=}") + print(f"{utils.stats(dnum_py[bads])=}") + print(f"{utils.stats(dnum_cy[bads])=}") + assert not np.any(bads), f"Found {utils.frac_str(bads)} inconsistent `dnum` b/t python and cython calcs!" + + # (3,) redz_final consistent between cython- and python- versions of calculation + + # ignore mismastch at low-redshifts + bads = (~np.isclose(redz_final_py, redz_final_cy, rtol=1e-2)) & redz_not_ignore + if np.any(bads): + print(f"{utils.frac_str(bads)=}") + print(f"{redz_final_py[bads][:10]=}") + print(f"{redz_final_cy[bads][:10]=}") + print(f"{utils.stats(redz_final_py[bads])=}") + print(f"{utils.stats(redz_final_cy[bads])=}") + assert not np.any(bads), f"Found {utils.frac_str(bads)} inconsistent `redz_final` b/t python and cython calcs!" + + # (4,) make sure that ALL numbers of binaries are consistent + + num_py = holo.sams.sam_cyutils.integrate_differential_number_3dx1d(edges_py, dnum_py) + num_cy = holo.sams.sam_cyutils.integrate_differential_number_3dx1d(edges_cy, dnum_cy) + + # Make sure that `atol` is also set to a reasonable value + bads = ~np.isclose(num_py, num_cy, rtol=1e-2, atol=1.0e-1) + if np.any(bads): + print(f"{utils.frac_str(bads)=}") + print(f"{num_py[bads][:10]=}") + print(f"{num_cy[bads][:10]=}") + print(f"{utils.stats(num_py[bads])=}") + print(f"{utils.stats(num_cy[bads])=}") + assert not np.any(bads), f"Found {utils.frac_str(bads)} inconsistent `num` b/t python and cython calcs!" return diff --git a/notebooks/semi-analytic-models.ipynb b/notebooks/semi-analytic-models.ipynb index 95d54de1..24fcc638 100644 --- a/notebooks/semi-analytic-models.ipynb +++ b/notebooks/semi-analytic-models.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "8d47f1a3", + "id": "0", "metadata": {}, "source": [ "# `holodeck` - Semi-Analytic Models" @@ -10,7 +10,7 @@ }, { "cell_type": "markdown", - "id": "af59b140", + "id": "1", "metadata": {}, "source": [ "For more information on the holodeck SAMs, see the [holodeck getting started guide](https://holodeck-gw.readthedocs.io/en/main/getting_started/index.html)." @@ -19,7 +19,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f9a9d9bc", + "id": "2", "metadata": {}, "outputs": [], "source": [ @@ -41,7 +41,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "630dd2a1", + "id": "3", "metadata": {}, "source": [ "## Quick Start" @@ -50,7 +50,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "f691ee46", + "id": "4", "metadata": {}, "source": [ "Construct a Semi-Analytic Model (SAM) using all of the default components" @@ -59,7 +59,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4f80bc42", + "id": "5", "metadata": {}, "outputs": [], "source": [ @@ -71,7 +71,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "5bb4f96a", + "id": "6", "metadata": {}, "source": [ "Choose the edges of the frequency bins at which to calculate the GWB" @@ -80,7 +80,7 @@ { "cell_type": "code", "execution_count": null, - "id": "a4e14d8e", + "id": "7", "metadata": {}, "outputs": [], "source": [ @@ -95,7 +95,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "34f11e96", + "id": "8", "metadata": {}, "source": [ "Calculate GWB from this SAM. \n", @@ -107,7 +107,7 @@ { "cell_type": "code", "execution_count": null, - "id": "abde3ae2", + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -119,7 +119,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "ebc65eca", + "id": "10", "metadata": {}, "source": [ "Plot GWB over multiple realizations" @@ -128,7 +128,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d836901c", + "id": "11", "metadata": {}, "outputs": [], "source": [ @@ -139,7 +139,7 @@ }, { "cell_type": "markdown", - "id": "0e256c02", + "id": "12", "metadata": {}, "source": [ "Slightly fancier plot:" @@ -148,7 +148,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b69aed2a", + "id": "13", "metadata": {}, "outputs": [], "source": [ @@ -185,7 +185,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "4356677c", + "id": "14", "metadata": {}, "source": [ "## Specifics" @@ -193,7 +193,7 @@ }, { "cell_type": "markdown", - "id": "3dbf39b8", + "id": "15", "metadata": {}, "source": [ "### Constructing a SAM" @@ -201,7 +201,7 @@ }, { "cell_type": "markdown", - "id": "d736718a", + "id": "16", "metadata": {}, "source": [ "SAMs are built from simple analytic models to derive the number-density of MBH binaries.\n", @@ -222,7 +222,7 @@ }, { "cell_type": "markdown", - "id": "8e8f3189", + "id": "17", "metadata": {}, "source": [ "The SAMs are initialized over a 3-dimensional parameter space of total MBH mass ($M = m_1 + m_2$), MBH mass ratio ($q = m_2 / m_1 \\leq 1$), and redshift ($z$). The `holodeck` code typically refers to the number of bins in each of these dimensions as `M`, `Q`, and `Z`; for example, the shape of the number-density of galaxy mergers will be `(M, Q, Z)`." @@ -231,7 +231,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7dacb18d", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -261,7 +261,7 @@ }, { "cell_type": "markdown", - "id": "7a971041", + "id": "19", "metadata": {}, "source": [ "### number density and the SAM grid" @@ -269,7 +269,7 @@ }, { "cell_type": "markdown", - "id": "ddd1b513", + "id": "20", "metadata": {}, "source": [ "The formation rate of MBH-MBH 'binaries' is calculated in `Semi_Analytic_Model.static_binary_density`, evaluated at the edges of the grid so that it's shape is the number of bins in each dimension, plus one, i.e. `(M+1, Q+1, Z+1)`. `static_binary_density` is implemented as a `@property` so that the first time the value is accessed it is calculated and cached, and then returned immediately on subsequent calls. \n", @@ -280,7 +280,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3cefecee", + "id": "21", "metadata": {}, "outputs": [], "source": [ @@ -301,7 +301,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4689a2eb", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -335,7 +335,7 @@ }, { "cell_type": "markdown", - "id": "dc7b2575", + "id": "23", "metadata": {}, "source": [ "### total number of binaries in a universe" @@ -343,7 +343,7 @@ }, { "cell_type": "markdown", - "id": "6a9d679e", + "id": "24", "metadata": {}, "source": [ "Above, we calculated the volumetric number-density rate of binary mergers. Here, we calculate the total number of binaries in a simulated universe at particular GW frequencies of interest. The SAM models currently assume circular binary orbits, so that the GW frequency is exactly twice the orbital frequency.\n", @@ -366,7 +366,7 @@ { "cell_type": "code", "execution_count": null, - "id": "31caefdd", + "id": "25", "metadata": {}, "outputs": [], "source": [ @@ -381,7 +381,7 @@ }, { "cell_type": "markdown", - "id": "875f617c", + "id": "26", "metadata": {}, "source": [ "For each point in the 3-dimensional SAM grid, we will be calculating the number of binaries at each frequency. So the returned values will be 4-dimensional with an additional axis with `F` frequency bins added: `(M, Q, Z, F)`." @@ -389,7 +389,7 @@ }, { "cell_type": "markdown", - "id": "894b20da", + "id": "27", "metadata": {}, "source": [ "**GW-Only Binary Evolution**" @@ -398,7 +398,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d0cd03af", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -411,13 +411,13 @@ "fobs_orb_cents = fobs_gw_cents / 2.0\n", "# `diff_num` is a 4D array with shape (M+1, Q+1, Z+1, F)\n", "# these values are evaluated at bin edges for total-mass, mass-ratio, redshift, but freq bin-centers\n", - "_edges, diff_num_gw, redz_final = sam.dynamic_binary_number_at_fobs(hard_gw, fobs_orb_cents)" + "_edges, diff_num_gw, redz_final = sam.dynamic_binary_number_at_fobs(hard_gw, fobs_orb_cents, use_cython=False)" ] }, { "cell_type": "code", "execution_count": null, - "id": "63496676", + "id": "29", "metadata": {}, "outputs": [], "source": [ @@ -442,7 +442,7 @@ }, { "cell_type": "markdown", - "id": "38215e23", + "id": "30", "metadata": {}, "source": [ "Convert from differential number of binaries to actual number of binaries" @@ -451,7 +451,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d3c9d78f", + "id": "31", "metadata": {}, "outputs": [], "source": [ @@ -465,9 +465,19 @@ "print(f\"Total number of modeled binaries: {number_gw.sum():.1e}\")" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "nden.shape, diff_num_gw.shape, number_gw.shape" + ] + }, { "cell_type": "markdown", - "id": "79774740", + "id": "33", "metadata": {}, "source": [ "Calculate the total number of binaries in certain ranges of parameter space" @@ -476,7 +486,7 @@ { "cell_type": "code", "execution_count": null, - "id": "986c695b", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -507,7 +517,7 @@ }, { "cell_type": "markdown", - "id": "0cbf07d3", + "id": "35", "metadata": {}, "source": [ "### Self-Consistent Binary Evolution (Phenomenological Model)" @@ -516,7 +526,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d704a3e7", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -532,7 +542,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3212639a", + "id": "37", "metadata": {}, "outputs": [], "source": [ @@ -566,7 +576,7 @@ }, { "cell_type": "markdown", - "id": "1cb3479a", + "id": "38", "metadata": {}, "source": [ "## Gravitational Waves" @@ -574,7 +584,7 @@ }, { "cell_type": "markdown", - "id": "857fab5c", + "id": "39", "metadata": {}, "source": [ "### Compare GWB and CW" @@ -583,7 +593,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b8098cc1", + "id": "40", "metadata": {}, "outputs": [], "source": [ @@ -604,7 +614,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b224fa4c", + "id": "41", "metadata": {}, "outputs": [], "source": [ @@ -628,7 +638,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "526f1ed3", + "id": "42", "metadata": {}, "source": [ "Calculate the distribution of GWB Amplitudes at 1/yr" @@ -637,7 +647,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e9fc9e99", + "id": "43", "metadata": {}, "outputs": [], "source": [ @@ -680,7 +690,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "4b567f98", + "id": "44", "metadata": {}, "source": [ "## Plot GWB Amplitude Distribution vs. M-MBulge parameters" @@ -689,7 +699,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "20f4383b", + "id": "45", "metadata": {}, "source": [ "Calculate GWB amplitudes at $f = 1/yr$ over a grid of M-Mbulge parameters, specifically the amplitude and power-law." @@ -698,7 +708,7 @@ { "cell_type": "code", "execution_count": null, - "id": "55cb90ca", + "id": "46", "metadata": {}, "outputs": [], "source": [ @@ -727,7 +737,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "cc13e780", + "id": "47", "metadata": {}, "source": [ "Plot the interquartile ranges for each power-law, as a function of normalization" @@ -736,7 +746,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22a7c9c1", + "id": "48", "metadata": {}, "outputs": [], "source": [ @@ -753,6 +763,159 @@ "plt.legend(title='M-MBulge Slope')\n", "plt.show()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "sam = holo.sams.Semi_Analytic_Model(shape=(10, 11, 12))\n", + "print(f\"{sam.shape=}\")\n", + "hard_gw = holo.hardening.Hard_GW()\n", + "\n", + "nden = sam.static_binary_density\n", + "\n", + "PTA_DUR = 20.0 * YR\n", + "NUM_FREQS = 9\n", + "fobs_gw_cents, fobs_gw_edges = holo.utils.pta_freqs(PTA_DUR, NUM_FREQS)\n", + "fobs_orb_cents = fobs_gw_cents / 2.0\n", + "fobs_orb_edges = fobs_gw_edges / 2.0\n", + "\n", + "grid_py, dnum_py, redz_final_py = sam.dynamic_binary_number_at_fobs(hard_gw, fobs_orb_cents, use_cython=False)\n", + "grid_cy, dnum_cy, redz_final_cy = sam.dynamic_binary_number_at_fobs(hard_gw, fobs_orb_cents, use_cython=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"grid\")\n", + "for ii in range(4):\n", + " print(ii, np.allclose(grid_py[ii], grid_cy[ii]))\n", + "\n", + "print(\"redz\", np.allclose(redz_final_py, redz_final_cy))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53", + "metadata": {}, + "outputs": [], + "source": [ + "diff_redz = (redz_final_py - redz_final_cy)\n", + "# idx_redz = (redz_final_cy > 0.0)\n", + "idx_redz = (redz_final_py > 0.0)\n", + "diff_redz[idx_redz] = diff_redz[idx_redz] / redz_final_cy[idx_redz]\n", + "\n", + "print(f\"{utils.frac_str(idx_redz)=}\")\n", + "print(f\"{utils.stats(diff_redz[~idx_redz])=}\")\n", + "print(f\"{utils.stats(diff_redz[idx_redz])=}\")\n", + "print(f\"{utils.stats(diff_redz)=}\")\n", + "\n", + "print()\n", + "\n", + "diff_dnum = (dnum_py - dnum_cy)\n", + "idx_dnum = (dnum_cy > 0.0)\n", + "diff_dnum[idx_dnum] = diff_dnum[idx_dnum] / dnum_cy[idx_dnum]\n", + "\n", + "print(f\"{utils.frac_str(idx_dnum)=}\")\n", + "print(f\"{utils.stats(diff_dnum[~idx_dnum])=}\")\n", + "print(f\"{utils.stats(diff_dnum[idx_dnum])=}\")\n", + "print(f\"{utils.stats(diff_dnum)=}\")\n", + "\n", + "print()\n", + "\n", + "print(f\"{utils.frac_str(idx_redz == idx_dnum)=}\")" + ] + }, + { + "cell_type": "markdown", + "id": "54", + "metadata": {}, + "source": [ + "## Different Redshifts" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "55", + "metadata": {}, + "outputs": [], + "source": [ + "bads_redz = ~np.isclose(redz_final_py, redz_final_cy)\n", + "print(f\"{utils.frac_str(bads_redz)=}\")\n", + "\n", + "print(redz_final_py[bads_redz][:10])\n", + "print(redz_final_cy[bads_redz][:10])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56", + "metadata": {}, + "outputs": [], + "source": [ + "# np.where(bads_redz)[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57", + "metadata": {}, + "outputs": [], + "source": [ + "SEL = 0\n", + "vals = np.meshgrid(*grid_py, indexing='ij')[SEL][bads_redz]\n", + "np.unique(vals/MSOL)" + ] + }, + { + "cell_type": "markdown", + "id": "58", + "metadata": {}, + "source": [ + "## Different Numbers" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59", + "metadata": {}, + "outputs": [], + "source": [ + "bads_dnum = ~np.isclose(dnum_py, dnum_cy)\n", + "print(f\"{utils.frac_str(bads_dnum)=}\")\n", + "\n", + "print(dnum_py[bads_dnum][:10])\n", + "print(dnum_cy[bads_dnum][:10])" + ] } ], "metadata": { @@ -774,7 +937,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.11.7" }, "toc": { "base_numbering": 1,