From 6b96d909ba1aa266be7cbb4b9a24791c01dc2638 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 12 Apr 2024 10:35:59 -0700 Subject: [PATCH 01/48] BUG: gw-only python version of dynamic_binary_number_at_fobs did not check for coalescence. Change default components for SAM --- CHANGELOG.md | 11 +++++++++ holodeck/sams/sam.py | 59 +++++++++++++++++++++++++++++--------------- 2 files changed, 50 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 641a5596..795a5b86 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,17 @@ ## Current +* DEFAULTS: + * Semi-Analytic Models + * New `Semi_Analytic_Models` instances now use galaxy merger-rates (instead of galaxy pair-fractions and merger-times) by default. To use GPF+GMT SAMs, the user must pass in at least a GPF instance manually. + +* BUGS: + * Semi-Analytic Models + * `Semi_Analytic_Model._dynamic_binary_number_at_fobs_inconsistent` was not checking for systems that had already coalesced. Only effected GW-only evolution using the python version of the calculation. + + +## v1.5.2 - 2024/04/12 + * DEPRECATIONS * `host_relations.py`: remove the `mamp` parameter and `MASS_AMP` attributes in the MMBulge relationships, and use `mamp_log10` and `MASS_AMP_LOG10` exclusively. diff --git a/holodeck/sams/sam.py b/holodeck/sams/sam.py index 1e00841a..0fe87b49 100644 --- a/holodeck/sams/sam.py +++ b/holodeck/sams/sam.py @@ -172,25 +172,25 @@ def __init__( gsmf = utils.get_subclass_instance(gsmf, None, _Galaxy_Stellar_Mass_Function) mmbulge = utils.get_subclass_instance(mmbulge, None, host_relations._MMBulge_Relation) - # if GMR is None, then we need both GMT and GPF - if gmr is None: - gmt = utils.get_subclass_instance(gmt, GMT_Power_Law, _Galaxy_Merger_Time) - gpf = utils.get_subclass_instance(gpf, GPF_Power_Law, _Galaxy_Pair_Fraction) - # if GMR is given, GMT can still be used - for calculating stalling - else: + if gpf is None: + log.info("No galaxy pair-fraction given, using galaxy merger-rate.") gmr = utils.get_subclass_instance(gmr, GMR_Illustris, _Galaxy_Merger_Rate) gmt = utils.get_subclass_instance(gmt, None, _Galaxy_Merger_Time, allow_none=True) - # if GMR is given, GPF is not used: make sure it is not given - if (gpf is not None): - err = f"When `GMR` ({gmr}) is provided, do not provide a GPF!" + else: + if gmr is not None: + err = "Can only use one of `gpf` and `gmr`!" log.exception(err) raise ValueError(err) + log.info("Galaxy pair-fraction provided, using galaxy pair-fraction and merger-time.") + gmt = utils.get_subclass_instance(gmt, GMT_Power_Law, _Galaxy_Merger_Time) + gpf = utils.get_subclass_instance(gpf, GPF_Power_Law, _Galaxy_Pair_Fraction) + self._gsmf = gsmf #: Galaxy Stellar-Mass Function (`_Galaxy_Stellar_Mass_Function` instance) - self._mmbulge = mmbulge #: Mbh-Mbulge relation (`host_relations._MMBulge_Relation` instance) + self._gmr = gmr #: Galaxy Merger Rate (`_Galaxy_Merger_Rate` instance) self._gpf = gpf #: Galaxy Pair Fraction (`_Galaxy_Pair_Fraction` instance) self._gmt = gmt #: Galaxy Merger Time (`_Galaxy_Merger_Time` instance) - self._gmr = gmr #: Galaxy Merger Rate (`_Galaxy_Merger_Rate` instance) + self._mmbulge = mmbulge #: Mbh-Mbulge relation (`host_relations._MMBulge_Relation` instance) log.debug(f"{gsmf=}, {gmr=}, {gpf=}, {gmt=}, {mmbulge=}") # ---- Create SAM grid edges @@ -577,7 +577,13 @@ def _dynamic_binary_number_at_fobs_inconsistent(self, hard, fobs_orb): shape = dens.shape new_shape = shape + (fobs_orb.size, ) - rz = self._redz_prime[..., np.newaxis] * np.ones(new_shape) + # If a galaxy merger-time is being used, then `_redz_prime` gives the redshifts following + # galaxy merger. NOTE: `_redz_prime` values past age of universe are set to `-1.0`. + # Use these as final redshifts if available, otherwise use the initial redshifts (`redz`). + if self._redz_prime is not None: + rz = self._redz_prime[..., np.newaxis] * np.ones(new_shape) + else: + rz = self.redz[np.newaxis, np.newaxis, :, np.newaxis] * np.ones(new_shape) coal = (rz > 0.0) dc = cosmo.comoving_distance(rz[coal]).to('Mpc').value @@ -585,27 +591,40 @@ def _dynamic_binary_number_at_fobs_inconsistent(self, hard, fobs_orb): fobs_orb[np.newaxis, np.newaxis, np.newaxis, :], rz ) + # (M,) ISCO separation + risco = utils.rad_isco(self.mtot) + # (Z,) this is `(dVc/dz) * (dz/dt)` in units of [Mpc^3/s] cosmo_fact = 4 * np.pi * (SPLC/MPC) * np.square(dc) * (1.0 + rz[coal]) - # # (M, Q) calculate chirp-mass + # broadcast to full shape, then take coalescing elements + # (M,) ==> (M, 1, 1, 1) mt = self.mtot[:, np.newaxis, np.newaxis, np.newaxis] + risco = risco[:, np.newaxis, np.newaxis, np.newaxis] + # (Q,) ==> (1, Q, 1, 1) mr = self.mrat[np.newaxis, :, np.newaxis, np.newaxis] - mt, mr = [(mm * np.ones(new_shape))[coal] for mm in [mt, mr]] + # These will now be 1D with shape (C,) for 'C' coalescing elements + mt, mr, fro, risco = [(mm * np.ones(new_shape))[coal] for mm in [mt, mr, frst_orb, risco]] # Convert from observer-frame orbital freq, to rest-frame orbital freq - sa = utils.kepler_sepa_from_freq(mt, frst_orb[coal]) - # (X, M*Q*Z), hardening rate, negative values, units of [cm/sec] + sa = utils.kepler_sepa_from_freq(mt, fro) + + # (C,), hardening rate, negative values, units of [cm/sec] args = [mt, mr, sa] dadt = hard.dadt(*args) # Calculate `tau = dt/dlnf_r = f_r / (df_r/dt)` # dfdt is positive (increasing frequency) - dfdt, _ = utils.dfdt_from_dadt(dadt, sa, frst_orb=frst_orb[coal]) - tau = frst_orb[coal] / dfdt + dfdt, _ = utils.dfdt_from_dadt(dadt, sa, frst_orb=fro) + tau = fro / dfdt - # (M, Q, Z) units: [1/s] i.e. number per second + # (M, Q, Z, F) units: [1/s] i.e. number per second dnum = np.zeros(new_shape) - dnum[coal] = (dens[..., np.newaxis] * np.ones(new_shape))[coal] * cosmo_fact * tau + + # Select only the binaries are separations larger than isco; this is a subset of `coal` + live = (sa > risco) + coal[coal] = coal[coal] & live + + dnum[coal] = (dens[..., np.newaxis] * np.ones(new_shape))[coal] * cosmo_fact[live] * tau[live] return edges, dnum, rz From f951c61d4d83c102155d6bae5cb92ec4735fc223 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 12 Apr 2024 10:36:17 -0700 Subject: [PATCH 02/48] Dont raise error for empty values passed to 'stats' --- holodeck/utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/holodeck/utils.py b/holodeck/utils.py index 790583cc..510dcb47 100644 --- a/holodeck/utils.py +++ b/holodeck/utils.py @@ -1042,7 +1042,7 @@ def rk4_step(func, x0, y0, dx, args=None, check_nan=0, check_nan_max=5): return x1, y1 -def stats(vals: npt.ArrayLike, percs: Optional[npt.ArrayLike] = None, prec: int = 2, weights=None) -> str: +def stats(vals, percs=None, prec=2, weights=None) -> str: """Return a string giving quantiles of the given input data. Parameters @@ -1065,8 +1065,8 @@ def stats(vals: npt.ArrayLike, percs: Optional[npt.ArrayLike] = None, prec: int """ try: - if len(vals) == 0: # type: ignore - raise TypeError + if len(vals) == 0: #### type: ignore + return "[]" except TypeError: raise TypeError(f"`vals` (shape={np.shape(vals)}) is not iterable!") From ad9b7e49e9a71e68fbaf60e35d8b232c0abbfd45 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 12 Apr 2024 10:36:46 -0700 Subject: [PATCH 03/48] Add unittests for comparing python and cython versions of dynamic_binary_number_at_fobs for GW-only evolution --- holodeck/sams/tests/test_sam.py | 62 +++++++++++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 2 deletions(-) diff --git a/holodeck/sams/tests/test_sam.py b/holodeck/sams/tests/test_sam.py index e6889247..36beb43c 100644 --- a/holodeck/sams/tests/test_sam.py +++ b/holodeck/sams/tests/test_sam.py @@ -4,7 +4,8 @@ import numpy as np import holodeck as holo -# from holodeck.constants import MSOL, PC, YR +from holodeck import utils +from holodeck.constants import MSOL, PC, YR def _test_sam_basics(sam): @@ -101,4 +102,61 @@ def test_sam_basics__gsmf_double_chechter(): sam = holo.sams.Semi_Analytic_Model(gsmf=gsmf, shape=SHAPE) _test_sam_basics(sam) - return \ No newline at end of file + return + + +# =========================================== +# ==== Test: dynamic_binary_number ==== +# =========================================== + + +def test_dbn_gw_only(): + """Test the dynamic_binary_number method using GW-only evolution. + + (1) runs without error + (2) no redz_final values should be <= 0.0 + (3) dnum values are consistent between cython and python + (4) redz_final values are consistent between cython and python + + """ + + shape = (10, 11, 12) + sam = holo.sams.Semi_Analytic_Model(shape=shape) + hard_gw = holo.hardening.Hard_GW() + + 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 + + # (1) + + 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) + + # (2) no redz_final values should be after redshift zero (i.e. negative, '-1.0') + + assert np.all(redz_final_py > 0.0), f"Found negative redshifts in python-version: {utils.stats(redz_final_py)=}" + assert np.all(redz_final_cy > 0.0), f"Found negative redshifts in cython-version: {utils.stats(redz_final_cy)=}" + + # (3,) dnum consistent between cython- and python- versions of calculation + + bads = ~np.isclose(dnum_py, dnum_cy) + 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), f"Found {utils.frac_str(bads)} inconsistent `dnum` b/t python and cython calcs!" + + # (4,) redz_final consistent between cython- and python- versions of calculation + bads = (~np.isclose(redz_final_py, redz_final_cy)) & (dnum_py > 0.0) + if np.any(bads): + print(f"{utils.frac_str(bads)=}") + 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!" + + + return + From fb30a3f369463419e70f450071cfa285767b47c3 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 12 Apr 2024 10:37:55 -0700 Subject: [PATCH 04/48] Fill in 'redz_final' values for gw-only evolution, even when binaries have already coalesced. This should never really make a difference, but avoids confusion. --- holodeck/sams/sam_cyutils.pyx | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/holodeck/sams/sam_cyutils.pyx b/holodeck/sams/sam_cyutils.pyx index 704e7ada..8981075d 100644 --- a/holodeck/sams/sam_cyutils.pyx +++ b/holodeck/sams/sam_cyutils.pyx @@ -452,6 +452,7 @@ def dynamic_binary_number_at_fobs(fobs_orb, sam, hard, cosmo): The differential number of binaries at each grid point. Unitless. """ + nden = sam.static_binary_density shape = sam.shape + (fobs_orb.size,) @@ -480,7 +481,8 @@ def dynamic_binary_number_at_fobs(fobs_orb, sam, hard, cosmo): elif isinstance(hard, holo.hardening.Hard_GW) or issubclass(hard, holo.hardening.Hard_GW): redz_prime = sam._redz_prime - # if `sam` is using galaxy merger rate (GMR), then `redz_prime` will be `None` + # if `sam` doesn't use a galaxy merger time (GMT), then `redz_prime` will be `None`, + # set to initial redshift values instead if redz_prime is None: sam._log.info("`redz_prime` not calculated in SAM. Setting to `redz` (initial) values.") redz_prime = sam.redz[np.newaxis, np.newaxis, :] * np.ones(sam.shape) @@ -567,6 +569,7 @@ cdef int _dynamic_binary_number_at_fobs_2pwl( The differential number of binaries at each grid point. Unitless. """ + cdef int n_mtot = mtot.size cdef int n_mrat = mrat.size cdef int n_redz = redz.size @@ -849,6 +852,10 @@ cdef int _dynamic_binary_number_at_fobs_gw( target_frst_orb = ftarget * (1.0 + rzp) # if target frequency is above ISCO freq, then all future ones will be also, so: break if target_frst_orb > frst_orb_isco: + # but still fill in the final redshifts (redz_final) + for fp in range(ff+1, n_freq): + redz_final[ii, jj, kk, fp] = rzp + break # get comoving distance From 43ceb6fb4a2c5288b32890ef9346f1dfc69baf04 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 11:05:04 -0700 Subject: [PATCH 05/48] BUG: initial redshift was being used to convert from obs-frame to rest-frame frequencies, instead of evolved redshifts --- holodeck/hardening.py | 15 ++++++++++----- holodeck/sams/sam.py | 42 +++++++++++++++++++++++++++++++++++------- 2 files changed, 45 insertions(+), 12 deletions(-) diff --git a/holodeck/hardening.py b/holodeck/hardening.py index acd54cef..4a610c69 100644 --- a/holodeck/hardening.py +++ b/holodeck/hardening.py @@ -60,11 +60,6 @@ _SCATTERING_DATA_FILENAME = "SHM06_scattering_experiments.json" -# ================================================================================================= -# ==== Hardening Classes ==== -# ================================================================================================= - - class _Hardening(abc.ABC): """Base class for binary-hardening models, providing the `dadt_dedt(evo, step, ...)` method. """ @@ -84,6 +79,11 @@ def dedt(self, *args, **kwargs): return rv_dedt +# ================================================================================================= +# ==== Physical Hardening Classes ==== +# ================================================================================================= + + class Hard_GW(_Hardening): """Gravitational-wave driven binary hardening. """ @@ -712,6 +712,11 @@ def _attenuation_BBR1980(self, sepa, m1m2, mstar): return atten +# ================================================================================================= +# ==== Phenomenological Hardening Classes ==== +# ================================================================================================= + + class Fixed_Time_2PL(_Hardening): r"""Provide a binary hardening rate such that the total lifetime matches a given value. diff --git a/holodeck/sams/sam.py b/holodeck/sams/sam.py index 0fe87b49..362c37b0 100644 --- a/holodeck/sams/sam.py +++ b/holodeck/sams/sam.py @@ -227,8 +227,9 @@ def __init__( # These values are calculated as needed by the class when the corresponding methods are called self._density = None #: Binary comoving number-density self._shape = None #: Shape of the parameter-space domain (mtot, mrat, redz) - self._gmt_time = None #: GMT timescale of galaxy mergers [sec] self._redz_prime = None #: redshift following galaxy merger process + #: GMT timescale of galaxy mergers [sec], set in `static_binary_density` + self._gmt_time = None return @@ -487,7 +488,7 @@ def _dynamic_binary_number_at_fobs_consistent(self, hard, fobs_orb, steps=200, d 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)[np.newaxis, :] # (1,X) + rads = np.linspace(0.0, 1.0, steps+1)[np.newaxis, :] # (1,X) # (M, S) = (M,1) * (1,S) rads = extr[0][:, np.newaxis] + (extr[1] - extr[0])[:, np.newaxis] * rads rads = 10.0 ** rads @@ -503,29 +504,56 @@ def _dynamic_binary_number_at_fobs_consistent(self, hard, fobs_orb, steps=200, d # (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) + + # 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, :] + self._gmt_time[:, :, :, np.newaxis] + times_tot = times_evo[:, :, np.newaxis, :] + if self._gmt_time is not None: + times_tot += 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}") # convert from separations to rest-frame orbital frequencies # (M, Q, S) frst_orb_evo = utils.kepler_freq_from_sepa(mt, rads) # (M, Q, Z, S) - fobs_orb_evo = frst_orb_evo[:, :, np.newaxis, :] / (1.0 + rz) + 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) - fobs_orb_evo, redz_evo = [mm.reshape(-1, steps-1) for mm in [fobs_orb_evo[:, :, :, 1:], redz_evo]] + fobs_orb_evo, redz_evo = [tt.reshape(-1, steps+1) for tt in [fobs_orb_evo, redz_evo]] # (M*Q*Z, X) - redz_final = utils.ndinterp(fobs_orb, fobs_orb_evo, redz_evo, xlog=True, ylog=False) + redz_final = utils.ndinterp(fobs_orb, fobs_orb_evo, redz_evo, xlog=False, ylog=False) # (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) frst_orb[frst_orb < 0.0] = 0.0 From 629d2883eb4e7a1c6dbfc09f397770b88b9a36c5 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 11:05:24 -0700 Subject: [PATCH 06/48] Add message to raised error --- holodeck/utils.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/holodeck/utils.py b/holodeck/utils.py index 510dcb47..a4e02c65 100644 --- a/holodeck/utils.py +++ b/holodeck/utils.py @@ -768,8 +768,9 @@ def ndinterp(xx, xvals, yvals, xlog=False, ylog=False): """ # assert np.ndim(xx) == 1 - assert np.ndim(xvals) == 2 - assert np.shape(xvals) == np.shape(yvals) + err = f"Bad shapes! {np.shape(xvals)=} {np.shape(yvals)=}" + assert np.ndim(xvals) == 2, err + assert np.shape(xvals) == np.shape(yvals), err xx = np.asarray(xx) xvals = np.asarray(xvals) From 5fdc5de9372d28dbc2f0779420481405fe00e72e Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 12:37:26 -0700 Subject: [PATCH 07/48] Add default value for number of integration steps when constructing cosmology instance --- holodeck/__init__.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/holodeck/__init__.py b/holodeck/__init__.py index 1a188106..4b55e1a9 100644 --- a/holodeck/__init__.py +++ b/holodeck/__init__.py @@ -79,7 +79,10 @@ class Parameters: # NOTE: Must load and initialize cosmology before importing other submodules! import cosmopy # noqa -cosmo = cosmopy.Cosmology(h=Parameters.HubbleParam, Om0=Parameters.Omega0, Ob0=Parameters.OmegaBaryon) +cosmo = cosmopy.Cosmology( + h=Parameters.HubbleParam, Om0=Parameters.Omega0, Ob0=Parameters.OmegaBaryon, + size=200, +) del cosmopy # ---- Import submodules From 20b3c6430e66d8fc57b0f583df04a9c8c7543230 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 12:37:44 -0700 Subject: [PATCH 08/48] DOCS: improve docstrings on 'redz_after' --- holodeck/utils.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/holodeck/utils.py b/holodeck/utils.py index a4e02c65..55465d94 100644 --- a/holodeck/utils.py +++ b/holodeck/utils.py @@ -1774,17 +1774,17 @@ def redz_after(time, redz=None, age=None): Parameters ---------- - time : array_like in units of [sec] - Amount of time to pass. - redz : None or array_like, - Redshift of starting point after which `time` is added. - age : None or array_like, in units of [sec] - Age of the Universe at the starting point, after which `time` is added. + time : array_like, [s] + Amount of time to pass, in units of seconds. + redz : None or array_like, [] + Redshift of starting point after which `time` is added. Unitless. + age : None or array_like, [s] + Age of the Universe at the starting point, after which `time` is added. Units of seconds. Returns ------- - new_redz : array_like - Redshift of the Universe after the given amount of time. + new_redz : array_like, [] + Redshift of the Universe after the given amount of time. Unitless """ if (redz is None) == (age is None): From ce0d88b8ccd3cf60afe2a4ab9ecde233e675dd18 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 12:41:35 -0700 Subject: [PATCH 09/48] Use lookback time instead of age of the universe when interpolating for much higher accuracy. Add tests for consistent dynamic_binary_number calculation. --- holodeck/librarian/param_spaces.py | 5 +- holodeck/sams/sam.py | 97 ++++++---- holodeck/sams/sam_cyutils.pyx | 29 ++- holodeck/sams/tests/test_sam.py | 93 +++++++++- notebooks/semi-analytic-models.ipynb | 263 ++++++++++++++++++++++----- 5 files changed, 392 insertions(+), 95 deletions(-) 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, From 21ef04c5195e3d0466c9c0fe94819378b9d1c814 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 12:43:37 -0700 Subject: [PATCH 10/48] DOCS: update changelog --- CHANGELOG.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 795a5b86..c434d11c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,12 +21,15 @@ * DEFAULTS: * Semi-Analytic Models - * New `Semi_Analytic_Models` instances now use galaxy merger-rates (instead of galaxy pair-fractions and merger-times) by default. To use GPF+GMT SAMs, the user must pass in at least a GPF instance manually. + * `Semi_Analytic_Models` instances now use galaxy merger-rates (instead of galaxy pair-fractions and merger-times) by default. To use GPF+GMT SAMs, the user must pass in at least a GPF instance manually. * BUGS: * Semi-Analytic Models * `Semi_Analytic_Model._dynamic_binary_number_at_fobs_inconsistent` was not checking for systems that had already coalesced. Only effected GW-only evolution using the python version of the calculation. +* Semi-Analytic Models + * Improve accuracy of dynamic binary number calculation for consistent evolution models. + ## v1.5.2 - 2024/04/12 From b6c34f913902f888523cca7689db89119f2cb69d Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 13:01:32 -0700 Subject: [PATCH 11/48] Decrease tolerance for error on mbulge test --- holodeck/tests/test_host_relations__bulge_frac.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/holodeck/tests/test_host_relations__bulge_frac.py b/holodeck/tests/test_host_relations__bulge_frac.py index b4b21d89..fa93eb61 100644 --- a/holodeck/tests/test_host_relations__bulge_frac.py +++ b/holodeck/tests/test_host_relations__bulge_frac.py @@ -74,7 +74,7 @@ def _check_any_bulge_frac(bf, mstar, redz): print(f"{dmstar_dmbulge_true=}") print(f"{error_deriv=}") # make sure values are consistent; calculation is numerical/fits to tolerate larger errors - assert np.allclose(dmstar_dmbulge_true, dmstar_dmbulge_test, rtol=1e-3) + assert np.allclose(dmstar_dmbulge_true, dmstar_dmbulge_test, rtol=1e-2) return mbulge From c2a20ef44583fa84ca0dd95dd4f8a87b26e3257c Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 13:50:58 -0700 Subject: [PATCH 12/48] Update PS_Astro_Strong with bulge-fraction --- holodeck/host_relations.py | 2 +- holodeck/librarian/libraries.py | 33 ++++++++++++++++++ holodeck/librarian/param_spaces.py | 55 +++++++++++++++++++++++++----- 3 files changed, 80 insertions(+), 10 deletions(-) diff --git a/holodeck/host_relations.py b/holodeck/host_relations.py index 22b9b8be..fed48c0b 100644 --- a/holodeck/host_relations.py +++ b/holodeck/host_relations.py @@ -779,7 +779,7 @@ class MMBulge_KH2013(MMBulge_Standard): """ # MASS_AMP = 0.49 * 1e9 * MSOL # 0.49 + 0.06 - 0.05 in units of [Msol] - MASS_AMP_LOG10 = 8.69 + MASS_AMP_LOG10 = 8.69 # 8.69 ± 0.05 [log10(M/Msol)] approximate uncertainties! MASS_REF = MSOL * 1e11 # 1e11 Msol MASS_PLAW = 1.17 # 1.17 ± 0.08 SCATTER_DEX = 0.28 # scatter stdev in dex diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/libraries.py index 3e3d66b1..3e3f1ff4 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/libraries.py @@ -209,11 +209,44 @@ def model_for_params(self, params, sam_shape=None): @classmethod @abc.abstractmethod def _init_sam(cls, sam_shape, params): + """Initialize a :class:`holodeck.sams.sam.Semi_Analytic_Model` instance with given params. + + Arguments + --------- + sam_shape : None or int or (3,) tuple of int + Shape of the SAM grid (M, Q, Z). If: + * `None` : use default values. + * int : apply this size to each dimension of the grid. + * (3,) of int : provide size for each dimension. + params : dict + Dictionary of parameters needed to initialize SAM model. + + Returns + ------- + sam : :class:`holodeck.sams.sam.Semi_Analytic_Model` instance, + Initialized SAM model instance. + + """ raise @classmethod @abc.abstractmethod def _init_hard(cls, sam, params): + """Initialize a :class:`holodeck.hardening._Hardening` subclass instance with given params. + + Arguments + --------- + sam : :class:`holodeck.sams.sam.Semi_Analytic_Model` instance, + SAM model instance. + params : dict + Dictionary of parameters needed to initialize hardening model. + + Returns + ------- + hard : :class:`holodeck.hardening._Hardening` subclass instance, + Initialized hardening model instance. + + """ raise def save(self, path_output): diff --git a/holodeck/librarian/param_spaces.py b/holodeck/librarian/param_spaces.py index fe34b5b0..6e14f5b3 100644 --- a/holodeck/librarian/param_spaces.py +++ b/holodeck/librarian/param_spaces.py @@ -166,7 +166,7 @@ class _PS_Astro_Strong(_Param_Space): hard_sepa_init=1e4, # [pc] hard_rchar=10.0, # [pc] hard_gamma_inner=-1.0, - hard_gamma_outer=+2.5, + hard_gamma_outer=0.0, # Galaxy stellar-mass Function (``GSMF_Double_Schechter``) # Parameters are based on `double-schechter.ipynb` conversions from [Leja2020]_ @@ -196,12 +196,14 @@ class _PS_Astro_Strong(_Param_Space): # M-MBulge Relationship (``MMBulge_KH2013``) # From [KH2013]_ - mmb_mamp=0.49e9, # 0.49e9 + 0.06 - 0.05 [Msol] + mmb_mamp_log10=8.69, # 8.69±0.05 [log10(M/Msol)] approx uncertainties! 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, + bf_frac_lo=0.4, + bf_frac_hi=0.8, + bf_mstar_crit=11.0, # [log10(M_star/M_Sol)] + bf_width_dex=1.0, # [dex] ) @classmethod @@ -242,10 +244,18 @@ def _init_sam(cls, sam_shape, params): qgammam=params['gmr_qgammam'], ) + # Mbh-MBulge relationship (and bulge-fractions) + bulge_frac = holo.host_relations.BF_Sigmoid( + bulge_frac_lo=params['bf_frac_lo'], + bulge_frac_hi=params['bf_frac_hi'], + mstar_char_log10=params['bf_mstar_crit'], + width_dex=params['bf_width_dex'], + ) mmbulge = holo.host_relations.MMBulge_KH2013( - mamp=params['mmb_mamp']*MSOL, + mamp_log10=params['mmb_mamp_log10'], mplaw=params['mmb_plaw'], scatter_dex=params['mmb_scatter_dex'], + bulge_frac=bulge_frac, ) sam = holo.sams.Semi_Analytic_Model( @@ -272,7 +282,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): parameters = [ # Hardening model (phenom 2PL) PD_Uniform("hard_time", 0.1, 11.0, default=3.0), # [Gyr] - PD_Uniform("hard_gamma_inner", -1.5, +0.0, default=-1.0), + PD_Uniform("hard_gamma_inner", -2.0, +0.0, default=-1.0), # GSMF PD_Normal('gsmf_log10_phi_one_z0', -2.383, 0.028), # - 2.383 ± 0.028 @@ -298,9 +308,15 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): PD_Normal('gmr_qgammaz', +0.0611, 0.0021), # +0.0611 ± 0.0021 beta1 PD_Normal('gmr_qgammam', -0.0477, 0.0013), # -0.0477 ± 0.0013 gamma + # MMBulge # From [KH2013]_ - PD_Normal('mmb_mamp', 0.49e9, 0.055e9), # 0.49e9 + 0.06 - 0.05 [Msol] + PD_Normal('mmb_mamp_log10', 8.69, 0.05), # 8.69 ± 0.05 [log10(M/Msol)] PD_Normal('mmb_plaw', 1.17, 0.08), # 1.17 ± 0.08 + # Extra + PD_Normal('mmb_scatter_dex', 0.28, 0.05), # no uncertainties given + PD_Uniform('bf_frac_lo', 0.0, 0.4), + PD_Uniform('bf_frac_hi', 0.6, 1.0), + PD_Uniform('bf_width_dex', 0.5, 1.5), # [dex] ] _Param_Space.__init__( self, parameters, @@ -315,7 +331,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): parameters = [ # Hardening model (phenom 2PL) PD_Uniform("hard_time", 0.1, 11.0, default=3.0), # [Gyr] - PD_Uniform("hard_gamma_inner", -1.5, +0.0, default=-1.0), + PD_Uniform("hard_gamma_inner", -2.0, +0.0, default=-1.0), ] _Param_Space.__init__( self, parameters, @@ -370,13 +386,33 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): return +class PS_Astro_Strong_MMBulge_BFrac(_PS_Astro_Strong): + + def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): + parameters = [ + # MMbulge - from [KH2013]_ + PD_Normal('mmb_mamp_log10', 8.69, 0.05), # 8.69 ± 0.05 [log10(M/Msol)] + PD_Normal('mmb_plaw', 1.17, 0.08), # 1.17 ± 0.08 + PD_Normal('mmb_scatter_dex', 0.28, 0.05), # no uncertainties given + PD_Uniform('bf_frac_lo', 0.0, 0.4), + PD_Uniform('bf_frac_hi', 0.6, 1.0), + PD_Uniform('bf_width_dex', 0.5, 1.5), # [dex] + ] + _Param_Space.__init__( + self, parameters, + log=log, nsamples=nsamples, sam_shape=sam_shape, seed=seed, + ) + return + + class PS_Astro_Strong_MMBulge(_PS_Astro_Strong): def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): parameters = [ # MMbulge - from [KH2013]_ - PD_Normal('mmb_mamp', 0.49e9, 0.055e9), # 0.49e9 + 0.06 - 0.05 [Msol] + PD_Normal('mmb_mamp_log10', 8.69, 0.05), # 8.69 ± 0.05 [log10(M/Msol)] PD_Normal('mmb_plaw', 1.17, 0.08), # 1.17 ± 0.08 + PD_Normal('mmb_scatter_dex', 0.28, 0.05), # no uncertainties given ] _Param_Space.__init__( self, parameters, @@ -390,6 +426,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): 'PS_Astro_Strong_Hard': PS_Astro_Strong_Hard, 'PS_Astro_Strong_GSMF': PS_Astro_Strong_GSMF, 'PS_Astro_Strong_GMR': PS_Astro_Strong_GMR, + 'PS_Astro_Strong_MMBulge_BFrac': PS_Astro_Strong_MMBulge_BFrac, 'PS_Astro_Strong_MMBulge': PS_Astro_Strong_MMBulge, } From 8fe4a62b182742b100232acc41125ffc863946b7 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Sat, 13 Apr 2024 14:35:10 -0700 Subject: [PATCH 13/48] Cleanup --- holodeck/librarian/combine.py | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index ded49fd6..321f0617 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -100,23 +100,6 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only # ---- load parameter space from save file - ''' - if path_pspace is None: - # look for parameter-space save files - regex = "*" + holo.librarian.PSPACE_FILE_SUFFIX # "*.pspace.npz" - files = sorted(path_output.glob(regex)) - num_files = len(files) - msg = f"found {num_files} pspace.npz files in {path_output}" - log.info(msg) - if num_files != 1: - log.exception(msg) - log.exception(f"{files=}") - log.exception(f"{regex=}") - raise RuntimeError(f"{msg}") - path_pspace = files[0] - - pspace = _Param_Space.from_save(path_pspace, log) - ''' if path_pspace is None: path_pspace = path_output pspace, pspace_fname = libraries.load_pspace_from_path(path_pspace, log=log) @@ -202,6 +185,11 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only log.warning(f"Saved to {lib_path}, size: {holo.utils.get_file_size(lib_path)}") + with h5py.File(lib_path, 'r') as h5: + assert np.all(h5['fobs_cents'][()] > 0.0) + if has_gwb: + log.info(f"Checking library file: {holo.utils.stats(gwb)=}") + return lib_path From e0a8d27f8b4e0ffc53389633181da165619b8648 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 15 Apr 2024 14:10:52 -0700 Subject: [PATCH 14/48] BUG: combine fails when first file (0th) is a failure file. Also add 'verbose' command-line option. --- holodeck/librarian/combine.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index 321f0617..63faf545 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -48,8 +48,15 @@ def main(): '--gwb', action='store_true', default=False, help='only merge the key GWB data (no single source, or binary parameter data).' ) + parser.add_argument( + '--verbose', '-v', action='store_true', default=False, + help="Verbose output in logger ('DEBUG' level)." + ) args = parser.parse_args() + if args.verbose: + log.setLevel(log.DEBUG) + log.debug(f"{args=}") path = Path(args.path) @@ -236,6 +243,12 @@ def _check_files_and_load_shapes(log, path_sims, nsamp): data_keys = list(temp.keys()) log.debug(f"{ii=} {temp_fname.name=} {data_keys=}") + if 'fail' in data_keys: + err = f"File {ii=} is a failed simulation file. {temp_fname=}" + log.error(err) + log.error(f"Error in file: {temp['fail']}") + continue + if fobs_cents is None: _fobs = temp.get('fobs', None) if _fobs is not None: From d8ace3a60b5a09d311b2f6ad5bda874bb446cc73 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 15 Apr 2024 14:11:33 -0700 Subject: [PATCH 15/48] BUG: update 'fit_spectra' to use 'fobs_cents' instead of 'fobs' in library files. Also handle NaN fit values in PSD. --- holodeck/librarian/fit_spectra.py | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/holodeck/librarian/fit_spectra.py b/holodeck/librarian/fit_spectra.py index b5138d7c..8f1a1435 100644 --- a/holodeck/librarian/fit_spectra.py +++ b/holodeck/librarian/fit_spectra.py @@ -1,7 +1,7 @@ """Script and methods to fit simulated GWB spectra with analytic functions. -Usage ------ +Usage (fit_spectra.py) +---------------------- For usage information, run the script with the ``-h`` or ``--help`` arguments, i.e.:: python -m holodeck.librarian.fit_spectra -h @@ -9,6 +9,13 @@ Typically, the only argument required is the path to the folder containing the combined library file (``sam_lib.hdf5``). +This script can be run serially (on a single processor) or in parallel. To run in parallel, MPI and +``mpi4py`` are required. For parallel runs, use:: + + mpirun -np python -m holodeck.librarian.fit_spectra + + + Notes ----- As a script, this submodule runs in parallel, with the main processor loading a holodeck library @@ -88,9 +95,9 @@ def fit_library_spectra(library_path, log, recreate=False): comm = MPI.COMM_WORLD except Exception as err: comm = None - holo.log.error(f"failed to load `mpi4py` in {__file__}: {err}") - holo.log.error("`mpi4py` may not be included in the standard `requirements.txt` file") - holo.log.error("Check if you have `mpi4py` installed, and if not, please install it") + log.error(f"failed to load `mpi4py` in {__file__}: {err}") + log.error("`mpi4py` may not be included in the standard `requirements.txt` file") + log.error("Check if you have `mpi4py` installed, and if not, please install it.") raise err # ---- setup path @@ -124,7 +131,7 @@ def fit_library_spectra(library_path, log, recreate=False): # ---- load library GWB and convert to PSD with h5py.File(library_path, 'r') as library: - fobs = library['fobs'][()] + fobs = library['fobs_cents'][()] psd = holo.utils.char_strain_to_psd(fobs[np.newaxis, :, np.newaxis], library['gwb'][()]) nsamps, nfreqs, nreals = psd.shape @@ -203,7 +210,10 @@ def fit_library_spectra(library_path, log, recreate=False): all_psd = all_psd[idx] # confirm that the resorting worked correctly - assert np.all(all_psd == psd) + matches = (all_psd == psd) + # if values are NaN, then equality check will fail... skip those + skips = ~np.isfinite(all_psd) + assert np.all(matches | skips) # reshape arrays to convert back to (Samples, Realizations, ...) len_nbins_plaw = len(nbins_plaw) @@ -219,7 +229,10 @@ def fit_library_spectra(library_path, log, recreate=False): all_psd = np.moveaxis(all_psd, 1, -1) # confirm that reshaping worked correctly - assert np.all(all_psd == psd_check) + matches = (all_psd == psd_check) + # if values are NaN, then equality check will fail... skip those + skips = ~np.isfinite(all_psd) + assert np.all(matches | skips) # Report how many fits failed fails = np.any(~np.isfinite(fits_plaw), axis=-1) From f53dcce1fac5eb543a63568e6ed935d94c1d0b74 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 15 Apr 2024 14:11:54 -0700 Subject: [PATCH 16/48] Add python console script for fit_spectra -- 'holodeck_fit_spec' --- setup.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 9c92a465..17fe2157 100644 --- a/setup.py +++ b/setup.py @@ -105,6 +105,9 @@ ext_modules=cython_modules, entry_points = { - "console_scripts": ['holodeck_lib_gen = holodeck.librarian.gen_lib:main'], + "console_scripts": [ + 'holodeck_lib_gen = holodeck.librarian.gen_lib:main', + 'holodeck_fit_spec = holodeck.librarian.fit_spectra:main', + ], }, ) From de8ea7b03555f997d1bc71ef97376d7cfa6a85f8 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 15 Apr 2024 15:13:22 -0700 Subject: [PATCH 17/48] BUG: fix issue in example slurm job script. Modify 'astro-strong' param spaces --- holodeck/librarian/param_spaces.py | 21 +++++++++++++++++++++ scripts/run_holodeck_lib_gen.sh | 9 +++++++-- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/holodeck/librarian/param_spaces.py b/holodeck/librarian/param_spaces.py index 6e14f5b3..27a897ea 100644 --- a/holodeck/librarian/param_spaces.py +++ b/holodeck/librarian/param_spaces.py @@ -283,6 +283,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): # Hardening model (phenom 2PL) PD_Uniform("hard_time", 0.1, 11.0, default=3.0), # [Gyr] PD_Uniform("hard_gamma_inner", -2.0, +0.0, default=-1.0), + PD_Uniform("hard_rchar", 2.0, 20.0, default=10.0), # [pc] # GSMF PD_Normal('gsmf_log10_phi_one_z0', -2.383, 0.028), # - 2.383 ± 0.028 @@ -332,6 +333,25 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): # Hardening model (phenom 2PL) PD_Uniform("hard_time", 0.1, 11.0, default=3.0), # [Gyr] PD_Uniform("hard_gamma_inner", -2.0, +0.0, default=-1.0), + PD_Uniform("hard_rchar", 2.0, 20.0, default=10.0), # [pc] + ] + _Param_Space.__init__( + self, parameters, + log=log, nsamples=nsamples, sam_shape=sam_shape, seed=seed, + ) + return + + +class PS_Astro_Strong_Hard_All(_PS_Astro_Strong): + + def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): + parameters = [ + # Hardening model (phenom 2PL) + PD_Uniform("hard_time", 0.1, 11.0, default=3.0), # [Gyr] + PD_Uniform("hard_gamma_inner", -2.0, +0.0, default=-1.0), + PD_Uniform("hard_rchar", 2.0, 20.0, default=10.0), # [pc] + PD_Uniform_Log("hard_sepa_init", 1e3, 1e4, default=1e4), # [pc] + PD_Uniform("hard_gamma_outer", 0.0, +2.5, default=0.0), ] _Param_Space.__init__( self, parameters, @@ -424,6 +444,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): _param_spaces_dict = { 'PS_Astro_Strong_All': PS_Astro_Strong_All, 'PS_Astro_Strong_Hard': PS_Astro_Strong_Hard, + 'PS_Astro_Strong_Hard_All': PS_Astro_Strong_Hard_All, 'PS_Astro_Strong_GSMF': PS_Astro_Strong_GSMF, 'PS_Astro_Strong_GMR': PS_Astro_Strong_GMR, 'PS_Astro_Strong_MMBulge_BFrac': PS_Astro_Strong_MMBulge_BFrac, diff --git a/scripts/run_holodeck_lib_gen.sh b/scripts/run_holodeck_lib_gen.sh index 1b37d50c..e341c632 100755 --- a/scripts/run_holodeck_lib_gen.sh +++ b/scripts/run_holodeck_lib_gen.sh @@ -9,11 +9,13 @@ # The library generation is performed using the `holodeck.librarian.lib_gen` script. # When holodeck is installed, this is aliased to the command-line function `holodeck_lib_gen`. # The following two commands should be identical: +# # python -m holodeck.librarian.lib_gen # holodeck_lib_gen # # The library generation script is parallelized using MPI (specifically mpi4py in python). # The actual run command can be very simple, for example: +# # mpirun -np 16 holodeck_lib_gen -n 64 -r 100 --gwb --ss --params PS_Classic_Phenom ./output/ # # When running on a remote server with a job scheduler (e.g. SLURM), things can be a bit more @@ -21,13 +23,16 @@ # computer you can ignore these elements. # # Remember to use `nohup` and `&` to run in the background on remote server, so that it's not killed -# when you logout. +# when you logout, e.g. +# +# nohup mpirun -np 16 holodeck_lib_gen -n 64 -r 100 --gwb --ss --params PS_Classic_Phenom ./output/ & # # ------------------ # Luke Zoltan Kelley # LZKelley@berkeley.edu # ================================================================================================== + # ==== SLURM job scheduler configuration ==== # NOTE: SLURM uses a single '#' to denote its configuration. Multiple '#' marks are ignored. @@ -136,7 +141,7 @@ echo -e "Running mpirun $(date +'%Y-%m-%d|%H:%M:%S')\n" echo "" # this is the actual call to run holodeck -mpirun -np $NTASKS $COMMAND $SPACE $OUTPUT -n $NSAMPS -r $NREALS -f $NFREQS $ARGS 1> $LOG_OUT 2> $LOG_ERR & +mpirun -np $NTASKS $COMMAND $SPACE $OUTPUT -n $NSAMPS -r $NREALS -f $NFREQS $ARGS 1> $LOG_OUT 2> $LOG_ERR echo "" echo -e "Completed python script $(date +'%Y-%m-%d|%H:%M:%S')\n" From a9e42bc397a17eba797698b46f21916ba2e22ad5 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Tue, 16 Apr 2024 13:54:30 -0700 Subject: [PATCH 18/48] BUG: missing import --- holodeck/librarian/param_spaces.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/holodeck/librarian/param_spaces.py b/holodeck/librarian/param_spaces.py index 27a897ea..53b13a3b 100644 --- a/holodeck/librarian/param_spaces.py +++ b/holodeck/librarian/param_spaces.py @@ -3,7 +3,7 @@ import holodeck as holo from holodeck.constants import GYR, PC, MSOL -from holodeck.librarian.libraries import _Param_Space, PD_Uniform, PD_Normal +from holodeck.librarian.libraries import _Param_Space, PD_Uniform, PD_Normal, PD_Uniform_Log # Define a new Parameter-Space class by subclassing the base class: @@ -447,7 +447,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): 'PS_Astro_Strong_Hard_All': PS_Astro_Strong_Hard_All, 'PS_Astro_Strong_GSMF': PS_Astro_Strong_GSMF, 'PS_Astro_Strong_GMR': PS_Astro_Strong_GMR, - 'PS_Astro_Strong_MMBulge_BFrac': PS_Astro_Strong_MMBulge_BFrac, 'PS_Astro_Strong_MMBulge': PS_Astro_Strong_MMBulge, + 'PS_Astro_Strong_MMBulge_BFrac': PS_Astro_Strong_MMBulge_BFrac, } From 30b446efd7c6fa19a80555a89b37ee250184c521 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Tue, 16 Apr 2024 13:54:49 -0700 Subject: [PATCH 19/48] BUG: outdated command; also add some more docstrings --- holodeck/librarian/gen_lib.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index 0a3b489a..ca1b958e 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -3,7 +3,8 @@ This file can be run from the command-line to generate holodeck libraries, and also provides some API methods for quick/easy generation of simulations. In general, these methods are designed to run simulations for populations constructed from parameter-spaces (i.e. -:class:`~holodeck.librarian.libraries._Param_Space` subclasses). +:class:`~holodeck.librarian.libraries._Param_Space` subclasses). This script is parallelized using +``mpi4py``, but can also be run in serial. This script can be run by executing:: @@ -11,6 +12,13 @@ Run ``python -m holodeck.librarian.gen_lib -h`` for usage information. +This script is also aliased to the console command, ``holodeck_lib_gen``. + +See https://holodeck-gw.readthedocs.io/en/main/getting_started/libraries.html for more information +about generating holodeck libraries and using holodeck parameter-spaces in general. + +For an example job-submission script using slurm, see the ``scripts/run_holodeck_lib_gen.sh`` file. + """ import argparse @@ -135,7 +143,7 @@ def main(): # noqa : ignore complexity warning if args.resume: # Load pspace object from previous save log.info(f"{args.resume=} attempting to load pspace {space_class=} from {args.output=}") - space, space_fname = holo.librarian.load_pspace_from_dir(log, args.output, space_class) + space, space_fname = holo.librarian.load_pspace_from_path(args.output, space_class=space_class, log=log) log.warning(f"resume={args.resume} :: Loaded param-space save from {space_fname}") else: space = space_class(log, args.nsamples, args.sam_shape, args.seed) From 183e1fa3141b21b8a5ca2c93a0b1a4dff37a96d9 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Tue, 16 Apr 2024 16:34:11 -0700 Subject: [PATCH 20/48] DOCS: minor --- scripts/run_holodeck_lib_gen.sh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/run_holodeck_lib_gen.sh b/scripts/run_holodeck_lib_gen.sh index e341c632..0c057818 100755 --- a/scripts/run_holodeck_lib_gen.sh +++ b/scripts/run_holodeck_lib_gen.sh @@ -2,6 +2,7 @@ # ================================================================================================== # Example job script to run holodeck library generation # ----------------------------------------------------- +# Last updated: 2024-04-15 # # Excute this script on a personal computer using: ``bash run_holodeck_lib_gen.sh``, # or on a server with slurm scheduler using: ``sbatch run_holodeck_lib_gen.sh``. @@ -18,6 +19,8 @@ # # mpirun -np 16 holodeck_lib_gen -n 64 -r 100 --gwb --ss --params PS_Classic_Phenom ./output/ # +# At the simplest level, this script runs a command like the above, and provides some conveniences. +# # When running on a remote server with a job scheduler (e.g. SLURM), things can be a bit more # complicated. Example SLURM configuration is included below, but if you're running on a personal # computer you can ignore these elements. From 0b70c2ba94122655ba64d942641292e002ed2f9d Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Tue, 16 Apr 2024 17:37:54 -0700 Subject: [PATCH 21/48] More information in error message --- holodeck/librarian/libraries.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/libraries.py index 3e3f1ff4..856ebae2 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/libraries.py @@ -336,6 +336,15 @@ def from_save(cls, fname, log=None): f"and class parameter names ({space.param_names})!" ) log.exception(err) + + for pn in param_names: + if pn not in space.param_names: + log.exception(f"parameter name `{pn}` in loaded data is not in class param names!") + + for pn in space.param_names: + if pn not in param_names: + log.exception(f"parameter name `{pn}` in class param names is not in loaded parameter names!") + raise RuntimeError(err) # Store loaded parameters into the parameter-space instance From 4ee3c1a63286c8a3d9aba26557ee557f2e6074f4 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Tue, 16 Apr 2024 17:38:36 -0700 Subject: [PATCH 22/48] Create a new save file in library output directories holding the general configuration data. Load this data when resuming runs. --- holodeck/librarian/gen_lib.py | 83 +++++++++++++++++++++++++++++++---- 1 file changed, 74 insertions(+), 9 deletions(-) diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index ca1b958e..c97f17c6 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -25,28 +25,27 @@ import sys from datetime import datetime from pathlib import Path +import json import shutil + import numpy as np -# import scipy as sp -# import scipy.stats import holodeck as holo -# from holodeck import cosmo from holodeck.constants import YR import holodeck.librarian import holodeck.librarian.combine from holodeck.librarian import ( libraries, - # DEF_NUM_FBINS, DEF_NUM_LOUDEST, DEF_NUM_REALS, DEF_PTA_DUR, ) -# from holodeck.sams import sam_cyutils MAX_FAILURES = 5 # FILES_COPY_TO_OUTPUT = [__file__, holo.librarian.__file__, holo.param_spaces.__file__] FILES_COPY_TO_OUTPUT = [] +ARGS_CONFIG_FNAME = "config.json" + def main(): # noqa : ignore complexity warning """Parent method for generating libraries from the command-line. @@ -144,7 +143,11 @@ def main(): # noqa : ignore complexity warning # Load pspace object from previous save log.info(f"{args.resume=} attempting to load pspace {space_class=} from {args.output=}") space, space_fname = holo.librarian.load_pspace_from_path(args.output, space_class=space_class, log=log) - log.warning(f"resume={args.resume} :: Loaded param-space save from {space_fname}") + # Load arguments/configuration from previous save + args, config_fname = _load_config(args.output, log) + + log.warning(f"resume={args.resume} :: Loaded param-space save from {space_fname}") + log.warning(f"resume={args.resume} :: Loaded configuration save from {config_fname}") else: space = space_class(log, args.nsamples, args.sam_shape, args.seed) else: @@ -152,6 +155,10 @@ def main(): # noqa : ignore complexity warning # share parameter space across processes space = comm.bcast(space, root=0) + # If we've loaded a new `args`, then share to all processes from rank=0 + if args.resume: + args = comm.bcast(args, root=0) + args.log = log log.info( f"param_space={args.param_space}, samples={args.nsamples}, sam_shape={args.sam_shape}, nreals={args.nreals}\n" @@ -175,11 +182,14 @@ def main(): # noqa : ignore complexity warning indices = comm.scatter(indices, root=0) iterator = holo.utils.tqdm(indices) if (comm.rank == 0) else np.atleast_1d(indices) - # Save parameter space to output directory + # Save parameter space and args/configuration to output directory if (comm.rank == 0) and (not args.resume): space_fname = space.save(args.output) log.info(f"saved parameter space {space} to {space_fname}") + config_fname = _save_config(args) + log.info(f"saved configuration to {config_fname}") + comm.barrier() beg = datetime.now() log.info(f"beginning tasks at {beg}") @@ -380,10 +390,17 @@ def _setup_argparse(*args, **kwargs): raise RuntimeError("`--params` requires the `--ss` option!") if args.resume: + # Need an existing output directory to resume from if not output.exists() or not output.is_dir(): - raise FileNotFoundError(f"`--resume` is active but output path does not exist! '{output}'") + err = f"`--resume` is active but output path does not exist! '{output}'" + raise FileNotFoundError(err) + + # Don't resume if we're recreating (i.e. erasing existing files) + if args.recreate: + err = "`resume` and `recreate` cannot both be set to True!" + raise ValueError(err) - # + # run in test mode if args.TEST: msg = "==== WARNING: running in test mode, other settings being overridden! ====" print("\n" + "=" * len(msg)) @@ -429,6 +446,54 @@ def _setup_argparse(*args, **kwargs): return args +def _save_config(args): + fname = args.output.joinpath(ARGS_CONFIG_FNAME) + + # Convert `args` parameters to serializable dictionary + config = {} + for kk, vv in args._get_kwargs(): + # convert `Path` instances to strings + if isinstance(vv, Path): + vv = str(vv) + config[kk] = vv + + # Add additional entries + config['holodeck_version'] = holo.__version__ + config['holodeck_librarian_version'] = holo.librarian.__version__ + config['holodeck_git_hash'] = holo.utils.get_git_hash() + config['created'] = str(datetime.now()) + + with open(fname, 'w') as out: + json.dump(config, out) + + print(f"Saved to {fname} - {holo.utils.get_file_size(fname)}") + + return fname + + +def _load_config(path, log): + fname = Path(path).joinpath(ARGS_CONFIG_FNAME) + + with open(fname, 'r') as inp: + config = json.load(inp) + + log.info("Loaded configuration from {fname}") + + pop_keys = [ + 'holodeck_version', 'holodeck_librarian_version', 'holodeck_git_hash', 'created' + ] + for pk in pop_keys: + val = config.pop(pk) + log.info(f"\t{pk}={val}") + + pspace = config.pop('param_space') + output = config.pop('output') + + args = _setup_argparse([pspace, output], **config) + + return args, fname + + def _setup_log(comm, args): """Setup up the logging module logger for output messaging. From 6560e588cf4b948596dd63ca76c881a8ed91c648 Mon Sep 17 00:00:00 2001 From: Cayenne Matt <59581529+CayenneMatt@users.noreply.github.com> Date: Wed, 17 Apr 2024 13:45:17 -0400 Subject: [PATCH 23/48] Add redz=redz as argument to dmstar_dmbh Add redz=redz as argument to dmstar_dmbh in line 358 for use in evolving version of mmbulge function --- holodeck/sams/sam.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/holodeck/sams/sam.py b/holodeck/sams/sam.py index 3f691459..f8761dbc 100644 --- a/holodeck/sams/sam.py +++ b/holodeck/sams/sam.py @@ -358,7 +358,7 @@ def static_binary_density(self): # ==> (dMstar-tot/dMbh-tot) = (dMstar-pri / dMbh-pri) * (dMbh-pri/dMbh-tot) / (dMstar-pri / dMstar-tot) # = (dMstar-pri / dMbh-pri) * (1 / (1+q_bh)) / (1 / (1+q_star)) # = (dMstar-pri / dMbh-pri) * ((1+q_star) / (1+q_bh)) - dmstar_dmbh_pri = self._mmbulge.dmstar_dmbh(mstar_pri) # [unitless] + dmstar_dmbh_pri = self._mmbulge.dmstar_dmbh(mstar_pri, redz=redz) # [unitless] qterm = (1.0 + mstar_rat) / (1.0 + self.mrat[np.newaxis, :, np.newaxis]) dmstar_dmbh = dmstar_dmbh_pri * qterm From dd0c7b70da5b6bc9f6cba2c3af7685bb531dcbf7 Mon Sep 17 00:00:00 2001 From: Cayenne Matt <59581529+CayenneMatt@users.noreply.github.com> Date: Wed, 17 Apr 2024 13:46:38 -0400 Subject: [PATCH 24/48] Add except statement to mbh_from_mbulge Add except ValueError statement to mbh_from_mbulge (within MMBulge_Redshift) on line 848 to catch another instance where redz doesn't need to be reshaped --- holodeck/host_relations.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/holodeck/host_relations.py b/holodeck/host_relations.py index fed48c0b..7e6688c7 100644 --- a/holodeck/host_relations.py +++ b/holodeck/host_relations.py @@ -845,6 +845,8 @@ def mbh_from_mbulge(self, mbulge, redz, scatter): redz = np.broadcast_to(redz, mbulge.T.shape).T except TypeError: redz = redz + except ValueError: + redz = redz zmamp = self._mamp * (1.0 + redz)**self._zplaw mbh = _log10_relation(mbulge, zmamp, self._mplaw, scatter_dex, x0=self._mref) return mbh From 736b9ecd30fa25291b65c831d1bdfcbd7eea3a75 Mon Sep 17 00:00:00 2001 From: Cayenne Matt <59581529+CayenneMatt@users.noreply.github.com> Date: Wed, 17 Apr 2024 13:59:56 -0400 Subject: [PATCH 25/48] Update try/except statement Changed line 846 to make the try/except statement more broad except TypeError: to except: --- holodeck/host_relations.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/holodeck/host_relations.py b/holodeck/host_relations.py index 7e6688c7..b1c107f1 100644 --- a/holodeck/host_relations.py +++ b/holodeck/host_relations.py @@ -843,9 +843,7 @@ def mbh_from_mbulge(self, mbulge, redz, scatter): # NOTE: this will work for (N,) ==> (N,) or (N,) ==> (N,X) try: redz = np.broadcast_to(redz, mbulge.T.shape).T - except TypeError: - redz = redz - except ValueError: + except: redz = redz zmamp = self._mamp * (1.0 + redz)**self._zplaw mbh = _log10_relation(mbulge, zmamp, self._mplaw, scatter_dex, x0=self._mref) From b765c76d11cba0e77af766f9f19a79d5888bc139 Mon Sep 17 00:00:00 2001 From: Cayenne Matt <59581529+CayenneMatt@users.noreply.github.com> Date: Thu, 18 Apr 2024 12:57:23 -0400 Subject: [PATCH 26/48] Update param_spaces_classic.py Update to fix some differences between PS_Classic_Phenom_Uniform and PS_Uniform_09B. Changed the range of parameter space for mmbulge amplitude from (7.5, 9.5) to (7.6, 9.0) and scatter from (0, 1.2) to (0, 0.9). --- holodeck/librarian/param_spaces_classic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/holodeck/librarian/param_spaces_classic.py b/holodeck/librarian/param_spaces_classic.py index 3509d319..8e90e705 100644 --- a/holodeck/librarian/param_spaces_classic.py +++ b/holodeck/librarian/param_spaces_classic.py @@ -100,8 +100,8 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): parameters = [ PD_Uniform("gsmf_phi0_log10", -3.5, -1.5), PD_Uniform("gsmf_mchar0_log10", 10.5, 12.5), # [log10(Msol)] - PD_Uniform("mmb_mamp_log10", +7.5, +9.5), # [log10(Msol)] - PD_Uniform("mmb_scatter_dex", +0.0, +1.2), + PD_Uniform("mmb_mamp_log10", +7.6, +9.0), # [log10(Msol)] + PD_Uniform("mmb_scatter_dex", +0.0, +0.9), PD_Uniform("hard_time", 0.1, 11.0), # [Gyr] PD_Uniform("hard_gamma_inner", -1.5, +0.0), ] From b471e5a57c1e577bc9ffbd35f10d8bd2ddc7a80f Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 22 Apr 2024 13:56:50 -0700 Subject: [PATCH 27/48] BUG: fix some issues with config files, resuming runs --- holodeck/librarian/combine.py | 2 ++ holodeck/librarian/gen_lib.py | 41 ++++++++++++++++++++++++++++++----- 2 files changed, 37 insertions(+), 6 deletions(-) diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index 63faf545..106b2016 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -226,6 +226,8 @@ def _check_files_and_load_shapes(log, path_sims, nsamp): has_gwb = False has_ss = False has_params = False + # num_fail = 0 + # num_good = 0 log.info(f"Checking {nsamp} files in {path_sims}") for ii in tqdm.trange(nsamp): diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index c97f17c6..4997df8c 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -39,13 +39,16 @@ libraries, ) -MAX_FAILURES = 5 +#: maximum number of failed simulations before task terminates with error (`None`: no limit) +MAX_FAILURES = None # FILES_COPY_TO_OUTPUT = [__file__, holo.librarian.__file__, holo.param_spaces.__file__] FILES_COPY_TO_OUTPUT = [] ARGS_CONFIG_FNAME = "config.json" +comm = None + def main(): # noqa : ignore complexity warning """Parent method for generating libraries from the command-line. @@ -143,11 +146,11 @@ def main(): # noqa : ignore complexity warning # Load pspace object from previous save log.info(f"{args.resume=} attempting to load pspace {space_class=} from {args.output=}") space, space_fname = holo.librarian.load_pspace_from_path(args.output, space_class=space_class, log=log) + log.warning(f"Loaded param-space save from {space_fname}") # Load arguments/configuration from previous save args, config_fname = _load_config(args.output, log) - - log.warning(f"resume={args.resume} :: Loaded param-space save from {space_fname}") - log.warning(f"resume={args.resume} :: Loaded configuration save from {config_fname}") + log.warning(f"Loaded configuration save from {config_fname}") + args.resume = True else: space = space_class(log, args.nsamples, args.sam_shape, args.seed) else: @@ -194,6 +197,7 @@ def main(): # noqa : ignore complexity warning beg = datetime.now() log.info(f"beginning tasks at {beg}") failures = 0 + num_done = 0 # ---- iterate over each processors' jobs @@ -208,15 +212,18 @@ def main(): # noqa : ignore complexity warning log.info(msg) rv, _sim_fname = run_sam_at_pspace_num(args, space, par_num) + if rv is False: failures += 1 - if failures > MAX_FAILURES: + if (MAX_FAILURES is not None) and (failures > MAX_FAILURES): log.error("\n\n") err = f"Failed {failures} times on rank:{comm.rank}!" log.exception(err) raise RuntimeError(err) + num_done += 1 + end = datetime.now() dur = (end - beg) log.info(f"\t{comm.rank} done at {str(end)} after {str(dur)} = {dur.total_seconds()}") @@ -276,8 +283,13 @@ def run_sam_at_pspace_num(args, space, pnum): if sim_fname.exists(): log.info(f"File {sim_fname} already exists. {args.recreate=}") + temp = np.load(sim_fname) + data_keys = list(temp.keys()) + + if 'fail' in data_keys: + log.info("Existing file was a failure, re-attempting...") # skip existing files unless we specifically want to recreate them - if not args.recreate: + elif not args.recreate: return True, sim_fname # ---- run Model @@ -447,14 +459,21 @@ def _setup_argparse(*args, **kwargs): def _save_config(args): + import logging + fname = args.output.joinpath(ARGS_CONFIG_FNAME) # Convert `args` parameters to serializable dictionary config = {} for kk, vv in args._get_kwargs(): + print(kk) + args.log.warning(f"{kk=} {vv=}") # convert `Path` instances to strings if isinstance(vv, Path): vv = str(vv) + # cannot store `logging.Logger` instances (`log`) + elif isinstance(vv, logging.Logger): + continue config[kk] = vv # Add additional entries @@ -692,3 +711,13 @@ def make_pars_plot(fobs, hc_ss, hc_bg, sspar, bgpar): if __name__ == "__main__": main() + + #! the below doesn't work for catching errors... maybe because of comm.barrier() calls? + # try: + # main() + # except Exception as err: + # print(f"Exception while running gen_lib.py::main() - '{err}'!") + # if comm is not None: + # comm.Abort() + # raise + #! ---------- From 05b8b7f1fdf88db66554d4d1887874e4fd569779 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 08:50:57 -0700 Subject: [PATCH 28/48] DOCS --- holodeck/librarian/combine.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index 106b2016..3dcec7cb 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -1,7 +1,12 @@ """Combine output files from individual simulation runs into a single library hdf5 file. This file can be executed as a script (see the :func:`main` function), and also provides an API -method (:func:`sam_lib_combine`) for programatically combining libraries. +method (:func:`sam_lib_combine`) for programatically combining libraries. When running as a script +or independent program, it must be run serially (not in parallel). + +For command-line usage, run: + + python -m holodeck.librarian.combine -h """ @@ -23,6 +28,8 @@ def main(): log = holo.log + # ---- Make sure we're NOT running in parallel (MPI) + try: from mpi4py import MPI comm = MPI.COMM_WORLD @@ -34,6 +41,8 @@ def main(): log.exception(err) raise RuntimeError(err) + # ---- Setup and parse command-line arguments + parser = argparse.ArgumentParser() parser.add_argument( @@ -60,6 +69,8 @@ def main(): log.debug(f"{args=}") path = Path(args.path) + # ---- Combine library files + sam_lib_combine(path, log, recreate=args.recreate, gwb_only=args.gwb) return From f9c31c7cabb6eaa542e026c5c52fc3b552fa30dd Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 09:18:58 -0700 Subject: [PATCH 29/48] Move constants to root of librarian module --- holodeck/librarian/__init__.py | 3 ++- holodeck/librarian/combine.py | 2 +- holodeck/librarian/gen_lib.py | 4 +--- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/holodeck/librarian/__init__.py b/holodeck/librarian/__init__.py index 0c18bfde..f75dc1e6 100644 --- a/holodeck/librarian/__init__.py +++ b/holodeck/librarian/__init__.py @@ -39,7 +39,7 @@ """ -__version__ = "1.1" +__version__ = "1.2" DEF_NUM_REALS = 100 #: Default number of realizations to construct in libraries. DEF_NUM_FBINS = 40 #: Default number of frequency bins at which to calculate GW signals. @@ -51,6 +51,7 @@ FNAME_SIM_FILE = "sam-lib__p{pnum:06d}.npz" PSPACE_FILE_SUFFIX = ".pspace.npz" +ARGS_CONFIG_FNAME = "config.json" from . libraries import ( # noqa _Param_Space, _Param_Dist, diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index 106b2016..dad03e99 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -322,7 +322,7 @@ def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log # When a processor fails for a given parameter, the output file is still created with the 'fail' key added if ('fail' in temp): msg = f"file {pnum=:06d} is a failure file, setting values to NaN ({fname})" - log.warning(msg) + log.info(msg) # set all parameters to NaN for failure files. Note that this is distinct from gwb=0.0 which can be real. if gwb is not None: gwb[pnum, :, :] = np.nan diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index 4997df8c..0c7c0413 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -36,7 +36,7 @@ import holodeck.librarian import holodeck.librarian.combine from holodeck.librarian import ( - libraries, + libraries, ARGS_CONFIG_FNAME ) #: maximum number of failed simulations before task terminates with error (`None`: no limit) @@ -45,8 +45,6 @@ # FILES_COPY_TO_OUTPUT = [__file__, holo.librarian.__file__, holo.param_spaces.__file__] FILES_COPY_TO_OUTPUT = [] -ARGS_CONFIG_FNAME = "config.json" - comm = None From a3f462d943f2db03e5cc35a655d13c16a8a8f66a Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 09:19:21 -0700 Subject: [PATCH 30/48] Add script for copying sets of libraries into a single output folder for transfer --- scripts/copy_sam_libs.py | 230 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 230 insertions(+) create mode 100644 scripts/copy_sam_libs.py diff --git a/scripts/copy_sam_libs.py b/scripts/copy_sam_libs.py new file mode 100644 index 00000000..238bfe60 --- /dev/null +++ b/scripts/copy_sam_libs.py @@ -0,0 +1,230 @@ +"""Copy holodeck SAM library files from a set of output directories into a single output folder. + +For usage information, run: + + python copy_sam_libs.py -h + +Example: + + python holodeck/scripts/copy_sam_libs.py ~/scratch/astro-strong ~/scratch/ -p "astro-strong-*" + + +Notes +----- +* The search for folders and library files are NOT recursive. +* The original copies are not deleted or modified. +* The copied files are automatically renamed. +* Configuration and pspace files are also copied and renamed. + +""" + +from pathlib import Path +import re +import shutil + +import argparse + +from holodeck.librarian import ( + ARGS_CONFIG_FNAME, PSPACE_FILE_SUFFIX +) + +ALLOW_MISSING_CONF = True +ALLOW_MISSING_PSPACE = False + + +def main(): + + # load command-line arguments + args = _setup_argparse() + + # find files + if args.debug: + print("-"*10 + " finding matching files " + "-"*10) + lib_files, pspace_files, conf_files = find_files(args, "sam_lib.hdf5") + + # copy files to output directory + if args.debug: + print("-"*10 + " copying files to new directory " + "-"*10) + copy_files(args, conf_files) + copy_files(args, pspace_files) + copy_files(args, lib_files) + + return + + +def find_files(args, file_pattern): + """ + + Find folders within the 'start_path' that match 'folder_pattern' (not recursive), then return + files matching 'file_pattern' within those folders (not recursive). + """ + lib_files = [] + pspace_files = [] + conf_files = [] + start_path = _to_abs_path(args.search_path) + # Compile the regex patterns + file_regex = re.compile(re.escape(file_pattern)) + folder_pattern = args.pattern + # folder_regex = re.compile(re.escape("*" + folder_pattern)) + if args.debug: + print(f"{start_path=}") + print(f"{file_pattern=} ==> {file_regex=}") + print(f"{folder_pattern=}") + + for path in start_path.glob(folder_pattern): + if not path.is_dir(): + continue + + # for path in start_path.rglob('*') + # if not folder_regex.match(str(path)): + # continue + + if args.debug: + print(f"Found {path=} ...") + + for file in path.glob('*'): + # Check if the file matches the file pattern + if not file.is_file() or not file_regex.search(str(file)): + continue + + # store library file + lib_files.append(file) + if args.debug: + print(f"===> found {file=}") + + # get parameter-space save file + pspace_file = list(path.glob('*' + PSPACE_FILE_SUFFIX)) + if len(pspace_file) == 1: + pspace_files.append(pspace_file[0]) + else: + err = f"Could not find unambiguous parameter-space file! matches={pspace_file}" + if ALLOW_MISSING_PSPACE: + print(err) + pspace_files.append(None) + else: + raise FileNotFoundError(err) + + # get configuration file + conf_file = path.joinpath(ARGS_CONFIG_FNAME) + if conf_file.is_file(): + conf_files.append(conf_file) + else: + err = f"Could not find configuration file! '{conf_file}'" + if ALLOW_MISSING_CONF: + print(err) + conf_files.append(None) + else: + raise FileNotFoundError(err) + + if args.debug: + print(f"Found {len(lib_files)} files.") + + return lib_files, pspace_files, conf_files + + +def copy_files(args, files): + """Copy all of the given files to the output (``args.output``) directory. + """ + + for src in files: + if src is None: + continue + src = Path(src) + assert src.is_file() + + folder = src.parts[-2] + + if args.rename: + new_name = folder + "_" + src.name + dst = args.output.joinpath(new_name) + else: + new_name = src.name + dst = args.output.joinpath(folder, new_name) + + if args.debug: + print(f"{src} ==>\n\t==> {dst}") + + if dst.exists() and not args.overwrite: + print(f"destination already exists, skipping! '{dst}'\n") + print("Use `--overwrite` to overwrite the file.") + + if not args.dry_run: + shutil.copy(src, dst) + assert dst.is_file() + if not args.debug: + print(dst) + + return + + +def _setup_argparse(*args, **kwargs): + + parser = argparse.ArgumentParser() + parser.add_argument('output', metavar='output', type=str, + help='output path [created if doesnt exist].') + + parser.add_argument('search_path', type=str, + help="where to start the search for matching folders.") + + parser.add_argument('-p', '--pattern', action='store', type=str, default='*', + help="regex for folders to match (NOTE: put this in quotations!).") + + parser.add_argument('-o', '--overwrite', action='store_true', default=False, + help="overwrite existing files [otherwise raise error].") + + # If this is disabled, it will cause problems with config and pspace files... + # so don't leave it as an option for now. See hard-coded value below. + #parser.add_argument('--rename', type=str_to_bool, nargs='?', default=True, + # help='rename the sam_lib files based on their folder name.') + + parser.add_argument('--debug', action='store_true', default=False, + help='debug.') + + parser.add_argument('--dry-run', action='store_true', default=False, dest='dry_run', + help='dry-run.') + + namespace = argparse.Namespace(**kwargs) + args = parser.parse_args(*args, namespace=namespace) + + # See note above, hard-code rename as true + args.rename = True + + # ---- check / sanitize arguments + + if args.dry_run: + print("`dry-run` is enabled. Settings `debug=True` automatically.") + args.debug = True + + args.output = _to_abs_path(args.output) + if args.debug: + print(f"absolute path: {args.output=}") + if args.output.is_file(): + raise RuntimeError(f"The output path is already a file! {output}") + args.output.mkdir(parents=True, exist_ok=True) + + return args + + +def str_to_bool(v): + """Convert string to boolean value.""" + if isinstance(v, bool): + return v + if v.lower() in ('yes', 'true', 't', 'y', '1'): + return True + elif v.lower() in ('no', 'false', 'f', 'n', '0'): + return False + else: + raise argparse.ArgumentTypeError('Boolean value expected.') + +def _to_abs_path(rpath): + apath = Path(rpath).resolve() + if not apath.is_absolute: + apath = Path('.').resolve() / apath + apath = apath.resolve() + return apath + + +if __name__ == "__main__": + main() + + From 7387365035014a29eb30a39109ac6b6c2cf80ecf Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 12:22:59 -0700 Subject: [PATCH 31/48] Cleanup; move code to setup pspace to separate function --- holodeck/librarian/gen_lib.py | 125 +++++++++--------- .../devs/libraries/library-explorer.ipynb | 42 ++++-- 2 files changed, 91 insertions(+), 76 deletions(-) diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index 4997df8c..740ea60a 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -71,14 +71,6 @@ def main(): # noqa : ignore complexity warning (6) All of the individual simulation files are combined using :func:`holodeck.librarian.combine.sam_lib_combine()`. - Arguments - --------- - None - - Returns - ------- - None - """ # ---- load mpi4py module @@ -108,8 +100,12 @@ def main(): # noqa : ignore complexity warning args.log = log if comm.rank == 0: - copy_files = FILES_COPY_TO_OUTPUT + + # get parameter-space class + space = _setup_param_space(args) + # copy certain files to output directory + copy_files = FILES_COPY_TO_OUTPUT if (not args.resume) and (copy_files is not None): for fname in copy_files: src_file = Path(fname) @@ -117,47 +113,34 @@ def main(): # noqa : ignore complexity warning shutil.copyfile(src_file, dst_file) log.info(f"Copied {fname} to {dst_file}") - # ---- get parameter-space class - - if comm.rank == 0: - - # `param_space` attribute must match the name of one of the classes in `holodeck.librarian` - try: - # if a namespace is specified for the parameter space, recursively follow it - # i.e. this will work in two cases: - # - `PS_Test` : if `PS_Test` is a class loaded in `librarian` - # - `file_name.PS_Test` : as long as `file_name` is a module within `librarian` - space_name = args.param_space.split(".") - if len(space_name) > 1: - space_class = holo.librarian - for class_name in space_name: - space_class = getattr(space_class, class_name) - else: - space_class = holo.librarian.param_spaces_dict[space_name[0]] - - except Exception as err: - log.error(f"Failed to load parameter space '{args.param_space}' !") - log.error("Make sure the class is defined, and imported into the `librarian` module.") - log.error(err) - raise err - - # instantiate the parameter space class + # Load arguments/configuration from previous save if args.resume: - # Load pspace object from previous save - log.info(f"{args.resume=} attempting to load pspace {space_class=} from {args.output=}") - space, space_fname = holo.librarian.load_pspace_from_path(args.output, space_class=space_class, log=log) - log.warning(f"Loaded param-space save from {space_fname}") - # Load arguments/configuration from previous save args, config_fname = _load_config(args.output, log) log.warning(f"Loaded configuration save from {config_fname}") args.resume = True + # Save parameter space and args/configuration to output directory else: - space = space_class(log, args.nsamples, args.sam_shape, args.seed) + space_fname = space.save(args.output) + log.info(f"saved parameter space {space} to {space_fname}") + + config_fname = _save_config(args) + log.info(f"saved configuration to {config_fname}") + + # Split and distribute index numbers to all processes + npars = args.nsamples + indices = range(npars) + indices = np.random.permutation(indices) + indices = np.array_split(indices, comm.size) + num_ind_per_proc = [len(ii) for ii in indices] + log.info(f"{npars=} cores={comm.size} || max runs per core = {np.max(num_ind_per_proc)}") + else: space = None + indices = None # share parameter space across processes space = comm.bcast(space, root=0) + # If we've loaded a new `args`, then share to all processes from rank=0 if args.resume: args = comm.bcast(args, root=0) @@ -167,40 +150,20 @@ def main(): # noqa : ignore complexity warning f"param_space={args.param_space}, samples={args.nsamples}, sam_shape={args.sam_shape}, nreals={args.nreals}\n" f"nfreqs={args.nfreqs}, pta_dur={args.pta_dur} [yr]\n" ) - comm.barrier() # ---- distribute jobs to processors - # Split and distribute index numbers to all processes - if comm.rank == 0: - npars = args.nsamples - indices = range(npars) - indices = np.random.permutation(indices) - indices = np.array_split(indices, comm.size) - num_ind_per_proc = [len(ii) for ii in indices] - log.info(f"{npars=} cores={comm.size} || max runs per core = {np.max(num_ind_per_proc)}") - else: - indices = None - indices = comm.scatter(indices, root=0) iterator = holo.utils.tqdm(indices) if (comm.rank == 0) else np.atleast_1d(indices) - # Save parameter space and args/configuration to output directory - if (comm.rank == 0) and (not args.resume): - space_fname = space.save(args.output) - log.info(f"saved parameter space {space} to {space_fname}") + comm.barrier() - config_fname = _save_config(args) - log.info(f"saved configuration to {config_fname}") + # ---- iterate over each processors' jobs - comm.barrier() beg = datetime.now() log.info(f"beginning tasks at {beg}") failures = 0 num_done = 0 - - # ---- iterate over each processors' jobs - for par_num in iterator: log.info(f"{comm.rank=} {par_num=}") @@ -458,6 +421,44 @@ def _setup_argparse(*args, **kwargs): return args +def _setup_param_space(args): + log = args.log + + # ---- Determine and load the parameter-space class + + # `param_space` attribute must match the name of one of the classes in `holodeck.librarian` + try: + # if a namespace is specified for the parameter space, recursively follow it + # i.e. this will work in two cases: + # - `PS_Test` : if `PS_Test` is a class loaded in `librarian` + # - `file_name.PS_Test` : as long as `file_name` is a module within `librarian` + space_name = args.param_space.split(".") + if len(space_name) > 1: + space_class = holo.librarian + for class_name in space_name: + space_class = getattr(space_class, class_name) + else: + space_class = holo.librarian.param_spaces_dict[space_name[0]] + + except Exception as err: + log.error(f"Failed to load parameter space '{args.param_space}' !") + log.error("Make sure the class is defined, and imported into the `librarian` module.") + log.error(err) + raise err + + # ---- instantiate the parameter space class + + if args.resume: + # Load pspace object from previous save + log.info(f"{args.resume=} attempting to load pspace {space_class=} from {args.output=}") + space, space_fname = holo.librarian.load_pspace_from_path(args.output, space_class=space_class, log=log) + log.warning(f"Loaded param-space save from {space_fname}") + else: + space = space_class(log, args.nsamples, args.sam_shape, args.seed) + + return space + + def _save_config(args): import logging diff --git a/notebooks/devs/libraries/library-explorer.ipynb b/notebooks/devs/libraries/library-explorer.ipynb index 68cc0248..fd377fa6 100644 --- a/notebooks/devs/libraries/library-explorer.ipynb +++ b/notebooks/devs/libraries/library-explorer.ipynb @@ -48,13 +48,12 @@ "outputs": [], "source": [ "lib_path = Path(\n", - " \"/Users/lzkelley/Programs/nanograv/holodeck/output\",\n", - " # \"astro-strong-hard-only\",\n", - " # \"astro-strong-all\",\n", - " \"astro-strong-all_n320_r100_f40\",\n", + " \"/Users/lzkelley/programs/nanograv/holodeck/output\",\n", + " # \"astro-strong-all_n320_r100_f40\",\n", + " \"astro-strong-mmbulge-bfrac_n512_r100_f40\",\n", " \"sam_lib.hdf5\"\n", ")\n", - "assert lib_path.exists()\n", + "assert lib_path.is_file(), f\"File {lib_path} does not exist!\"\n", "print(lib_path.parent)\n", "\n", "library = h5py.File(lib_path, 'r')\n", @@ -327,7 +326,7 @@ "spars = library['sample_params'][:].copy()\n", "idx = np.ones_like(spars[:, 0], dtype=bool)\n", "for sp in spars.T:\n", - " extr = np.percentile(sp, [1, 99])\n", + " extr = np.percentile(sp[np.isfinite(sp)], [1, 99])\n", " idx = idx & (extr[0] < sp) & (sp < extr[1])\n", "\n", "names = [dd.decode() for dd in library.attrs['param_names']]\n", @@ -356,13 +355,14 @@ "source": [ "CREATE = False\n", "RECREATE = False\n", - "fits_path = holo.librarian.get_fits_path(lib_path)\n", - "print(f\"{fits_path=} {fits_path.exists()}\")\n", + "import holodeck.librarian.fit_spectra\n", + "fits_path = holo.librarian.libraries.get_fits_path(lib_path)\n", + "print(f\"{fits_path.name=} {fits_path.exists()=}\")\n", "if CREATE or RECREATE:\n", " if (CREATE and (not fits_path.exists())) or RECREATE:\n", - " holo.librarian.fit_library_spectra(lib_path, log, recreate=RECREATE)\n", + " holo.librarian.fit_spectra.fit_library_spectra(lib_path, log, recreate=RECREATE)\n", "\n", - "assert fits_path.is_file()\n", + "assert fits_path.is_file(), f\"spectral fits file does not exist! {fits_path}\"\n", "with np.load(fits_path) as fits:\n", " print(list(fits.keys()))" ] @@ -391,7 +391,7 @@ "skip = None\n", "\n", "len_nbins = len(nbins_plaw)\n", - "fig, axes = plt.subplots(figsize=[12, 6], nrows=3, ncols=len_nbins, sharey='row')\n", + "fig, axes = plt.subplots(figsize=[12, 12], nrows=3, ncols=len_nbins, sharey='row')\n", "plt.subplots_adjust(hspace=0.4, wspace=0.4)\n", "\n", "axrow = axes[0, :]\n", @@ -401,16 +401,30 @@ " ax.axvline(-15, ls='--', color='0.5', alpha=0.5, zorder=1000)\n", " ax.axhline(-13/3, ls='--', color='0.5', alpha=0.5, zorder=1000)\n", " xx, yy = med[::skip, ii, :].T\n", + " sel = np.isfinite(yy)\n", + " xx = xx[sel]\n", + " yy = yy[sel]\n", + " zz = std[::skip][sel]\n", + "\n", " kale.dist2d([xx, yy], hist=False, ax=ax, median=False)\n", "\n", " ax = axes[1, ii]\n", " ax.set(xlabel='log10(amp)', ylabel='stdev')\n", - " ax.scatter(xx, std[::skip, ii, 0], s=4, alpha=0.5)\n", + " ax.scatter(xx, zz[:, ii, 0], s=4, alpha=0.5)\n", "\n", " ax = axes[2, ii]\n", " ax.set(xlabel='gamma', ylabel='stdev')\n", - " ax.scatter(yy, std[::skip, ii, 1], s=4, alpha=0.5)\n", + " ax.scatter(yy, zz[:, ii, 1], s=4, alpha=0.5)\n", + "\n", + "ax = axes[0, 0]\n", + "ylim = np.array(ax.get_ylim())\n", + "print(f\"{ylim=}\")\n", + "\n", + "if ylim[1] > -13.0/3.0:\n", + " ylim[1] = -4.0\n", "\n", + "ax.set_ylim(ylim)\n", + "print(f\"{ylim=}\")\n", "\n", "plt.show()\n", "skip_str = \"\" if skip in [0, 1, None] else f\"_skip{skip}\"\n", @@ -1335,7 +1349,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.11.4" }, "vscode": { "interpreter": { From 195955136e2ee918edf93909d1839e9c95d35913 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 13:24:55 -0700 Subject: [PATCH 32/48] BUG: Try to fix issue with loggers not setting correct level when using both stream and file IO --- holodeck/host_relations.py | 1 - holodeck/logger.py | 16 ++++++++++++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/holodeck/host_relations.py b/holodeck/host_relations.py index b1c107f1..0c4d97f8 100644 --- a/holodeck/host_relations.py +++ b/holodeck/host_relations.py @@ -299,7 +299,6 @@ def mstar_from_mbulge(self, mbulge, redz=None, **kwargs): mstar = np.ones_like(mbulge) * mbulge / fhi # find the systems that are below the char mass, based on this assumption sel = (mstar/self._mstar_char) < 1.0 - print(f"{utils.frac_str(sel)=}") # interpolate to numerically invert the function # mstar[sel] = np.interp(mbulge[sel], self._grid_mbulge, self._grid_mstar) mstar[sel] = self._interp_mstar_from_mbulge(mbulge[sel]) diff --git a/holodeck/logger.py b/holodeck/logger.py index ed4de429..1ea128bf 100644 --- a/holodeck/logger.py +++ b/holodeck/logger.py @@ -70,10 +70,22 @@ def get_logger(name='holodeck', level_stream=logging.WARNING, tostr=sys.stdout, handler.setLevel(level_stream) logger.addHandler(handler) + # store these values for convenience for lvl in ['DEBUG', 'INFO', 'WARNING', 'ERROR']: setattr(logger, lvl, getattr(logging, lvl)) - # Make sure that the `setLevel` command reaches the stream logger - logger.setLevel = lambda xx: logger.handlers[0].setLevel(xx) + # ---- Make sure that the `setLevel` command reaches the stream logger + + # Construct a new function to replace 'setLevel' + def _set_level(self, lvl): + for handler in self.handlers: + if not isinstance(handler, logging.StreamHandler): + continue + handler.setLevel(lvl) + + return + + # replace `setLevel` function on the logger + logger.setLevel = _set_level.__get__(logger) return logger From 577912a69ef94171400f31cec1f009361660dfa7 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 13:25:35 -0700 Subject: [PATCH 33/48] When loading saved parameter-spaces, access class-name from save file --- holodeck/librarian/libraries.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/libraries.py index 856ebae2..70d2409e 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/libraries.py @@ -267,7 +267,8 @@ def save(self, path_output): """ log = self._log - my_name = self.__class__.__name__ + class_name = self.__class__.__name__ + class_vers = self.__version__ vers = holo.librarian.__version__ # make sure `path_output` is a directory, and that it exists @@ -277,16 +278,16 @@ def save(self, path_output): log.exception(err) raise ValueError(err) - fname = f"{my_name}{holo.librarian.PSPACE_FILE_SUFFIX}" + fname = f"{class_name}{holo.librarian.PSPACE_FILE_SUFFIX}" fname = path_output.joinpath(fname) - log.debug(f"{my_name=} {vers=} {fname=}") + log.debug(f"{class_name=} {vers=} {fname=}") data = {} for key in self._SAVED_ATTRIBUTES: data[key] = getattr(self, key) np.savez( - fname, class_name=my_name, librarian_version=vers, + fname, class_name=class_name, class_vers=class_vers, librarian_version=vers, **data, ) @@ -926,7 +927,7 @@ def _calc_model_details(edges, redz_final, number): def load_pspace_from_path(path, space_class=None, log=None): - """Load a `_Param_Space` subclass instance from the saved file in the given directory. + """Load a ``_Param_Space`` subclass instance from a save file with the given path. This function tries to determine the correct class based on the save file name, then uses that class to load the save itself. @@ -970,9 +971,21 @@ class to load the save itself. else: raise - # Based on the `space_fname`, try to find a matching PS (parameter-space) in `holodeck.param_spaces_dict` if space_class is None: - space_class = _get_space_class_from_space_fname(space_fname) + try: + space_class = str(np.load(space_fname)['class_name']) + print(f"{space_class=}") + space_class = holo.librarian.param_spaces_dict[space_class] + except Exception as err: + if log is not None: + log.error(f"Could not load `class_name` from save file '{space_fname}'.") + log.error(str(err)) + else: + print(f"Could not load `class_name` from save file '{space_fname}'.") + print(str(err)) + + # Based on the `space_fname`, try to find a matching PS (parameter-space) in `holodeck.param_spaces_dict` + space_class = _get_space_class_from_space_fname(space_fname) space = space_class.from_save(space_fname, log=log) return space, space_fname From 6c7054331df3899ee1f839d7f9cb0f92f8640280 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 13:26:33 -0700 Subject: [PATCH 34/48] In gen_lib, add in 'verbose' CL argument; work on setting up 'domain' exploration runs --- holodeck/librarian/gen_lib.py | 62 ++++++++++++++++++++--------------- 1 file changed, 36 insertions(+), 26 deletions(-) diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index 8f4cd703..a7bece01 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -28,7 +28,6 @@ import json import shutil - import numpy as np import holodeck as holo @@ -79,8 +78,8 @@ def main(): # noqa : ignore complexity warning except ModuleNotFoundError as err: comm = None holo.log.error(f"failed to load `mpi4py` in {__file__}: {err}") - holo.log.error("`mpi4py` may not be included in the standard `requirements.txt` file") - holo.log.error("Check if you have `mpi4py` installed, and if not, please install it") + holo.log.error("`mpi4py` may not be included in the standard `requirements.txt` file.") + holo.log.error("Check if you have `mpi4py` installed, and if not, please install it.") raise err # ---- setup arguments / settings, loggers, and outputs @@ -124,13 +123,17 @@ def main(): # noqa : ignore complexity warning config_fname = _save_config(args) log.info(f"saved configuration to {config_fname}") - # Split and distribute index numbers to all processes - npars = args.nsamples - indices = range(npars) - indices = np.random.permutation(indices) - indices = np.array_split(indices, comm.size) - num_ind_per_proc = [len(ii) for ii in indices] - log.info(f"{npars=} cores={comm.size} || max runs per core = {np.max(num_ind_per_proc)}") + # Split simulations for all processes + + if args.domain: + pass + # + else: + indices = range(args.nsamples) + indices = np.random.permutation(indices) + indices = np.array_split(indices, comm.size) + num_per = [len(ii) for ii in indices] + log.info(f"{args.nsamples=} cores={comm.size} || max sims per core = {np.max(num_per)}") else: space = None @@ -172,7 +175,8 @@ def main(): # noqa : ignore complexity warning msg = ", ".join(msg) log.info(msg) - rv, _sim_fname = run_sam_at_pspace_num(args, space, par_num) + params = space.param_dict(par_num) + rv, _sim_fname = run_sam_at_pspace_params(args, space, par_num, params) if rv is False: failures += 1 @@ -200,17 +204,17 @@ def main(): # noqa : ignore complexity warning return -def run_sam_at_pspace_num(args, space, pnum): +def run_sam_at_pspace_params(args, space, pnum, params): """Run a given simulation (index number ``pnum``) in the ``space`` parameter-space. This function performs the following: (1) Constructs the appropriate filename for this simulation, and checks if it already exist. If - the file exists and the ``args.recreate`` option is not specified, the function returns + the file exists and the ``args.recreate`` option is False, the function returns ``True``, otherwise the function runs this simulation. - (2) Calls ``space.model_for_sample_number`` to generate the semi-analytic model and hardening + (2) Calls ``space.model_for_params`` to generate the semi-analytic model and hardening instances; see the function - :func:`holodeck.librarian.libraries._Param_Space.model_for_sample_number()`. + :func:`holodeck.librarian.libraries._Param_Space.model_for_params()`. (3) Calculates populations and GW signatures from the SAM and hardening model using :func:`holodeck.librarian.libraries.run_model()`, and saves the results to an output file. (4) Optionally: some diagnostic plots are created in the :func:`make_plots()` function. @@ -218,20 +222,28 @@ def run_sam_at_pspace_num(args, space, pnum): Arguments --------- args : ``argparse.ArgumentParser`` instance - Arguments from the `gen_lib_sams.py` script. - NOTE: this should be improved. + Arguments from the ``gen_lib_sams.py`` script. space : :class:`holodeck.librarian.libraries._Param_space` instance Parameter space from which to construct populations. pnum : int Which parameter-sample from ``space`` should be run. + params : dict + Parameters for this particular simulation. Dictionary of key-value pairs obtained from + the parameter-space ``space`` and passed back to the parameter-space to produce a model. Returns ------- rv : bool ``True`` if this simulation was successfully run, ``False`` otherwise. + NOTE: an output file is created in either case. See ``Notes`` [1] below. sim_fname : ``pathlib.Path`` instance Path of the simulation save file. + Notes + ----- + [1] When simulations are run, an output file is produced. On caught failures, an output file is + produced that contains a single key: 'fail'. This designates the file as a failure. + """ log = args.log @@ -257,7 +269,7 @@ def run_sam_at_pspace_num(args, space, pnum): try: log.debug("Selecting `sam` and `hard` instances") - sam, hard = space.model_for_sample_number(pnum) + sam, hard = space.model_for_params(params) data = libraries.run_model( sam, hard, @@ -271,6 +283,7 @@ def run_sam_at_pspace_num(args, space, pnum): log.exception(f"`run_model` FAILED on {pnum=}\n") log.exception(err) rv = False + # failed simulations get an output file with a single key: 'fail' data = dict(fail=str(err)) # ---- save data to file @@ -343,8 +356,9 @@ def _setup_argparse(*args, **kwargs): parser.add_argument('--TEST', action='store_true', default=False, help='Run in test mode (NOTE: this resets other values)') - # parser.add_argument('-v', '--verbose', action='store_true', default=False, dest='verbose', - # help='verbose output [INFO]') + # metavar='LEVEL', type=int, default=20, + parser.add_argument('-v', '--verbose', metavar='LEVEL', type=int, nargs='?', const=20, default=30, + help='verbose output level (DEBUG=10, INFO=20, WARNING=30).') namespace = argparse.Namespace(**kwargs) args = parser.parse_args(*args, namespace=namespace) @@ -465,8 +479,6 @@ def _save_config(args): # Convert `args` parameters to serializable dictionary config = {} for kk, vv in args._get_kwargs(): - print(kk) - args.log.warning(f"{kk=} {vv=}") # convert `Path` instances to strings if isinstance(vv, Path): vv = str(vv) @@ -484,7 +496,7 @@ def _save_config(args): with open(fname, 'w') as out: json.dump(config, out) - print(f"Saved to {fname} - {holo.utils.get_file_size(fname)}") + args.log.warning(f"Saved to {fname} - {holo.utils.get_file_size(fname)}") return fname @@ -549,14 +561,12 @@ def _setup_log(comm, args): # ---- setup logger - # log_lvl = holo.logger.INFO if args.verbose else holo.logger.WARNING - log_lvl = holo.logger.DEBUG + log_lvl = args.verbose if comm.rank == 0 else holo.logger.DEBUG tostr = sys.stdout if comm.rank == 0 else False log = holo.logger.get_logger(name=log_name, level_stream=log_lvl, tofile=fname, tostr=tostr) log.info(f"Output path: {output}") log.info(f" log: {fname}") log.info(args) - return log From df998b50e8dfe53983e980d4e74a58fdff4ee1a1 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 24 Apr 2024 13:57:57 -0700 Subject: [PATCH 35/48] In 'fit_spectra' only require match to library file name ending, not whole name --- holodeck/librarian/fit_spectra.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/holodeck/librarian/fit_spectra.py b/holodeck/librarian/fit_spectra.py index 8f1a1435..4da08e60 100644 --- a/holodeck/librarian/fit_spectra.py +++ b/holodeck/librarian/fit_spectra.py @@ -280,7 +280,7 @@ def _find_sam_lib_in_path_tree(path, pattern=None): if (pattern is not None) and (pattern not in str(path)): return [] # if we find the library file, return it - if path.name == "sam_lib.hdf5": + if path.name.endswith("sam_lib.hdf5"): return [path] return [] From b21dbc7bde46a5145ea8b7fc31d21ae972ee89ed Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 26 Apr 2024 09:46:31 -0700 Subject: [PATCH 36/48] BUG: issue with loading saves --- holodeck/librarian/libraries.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/libraries.py index 70d2409e..8fbabfcb 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/libraries.py @@ -973,7 +973,7 @@ class to load the save itself. if space_class is None: try: - space_class = str(np.load(space_fname)['class_name']) + space_class = str(np.load(space_fname, allow_pickle=True)['class_name']) print(f"{space_class=}") space_class = holo.librarian.param_spaces_dict[space_class] except Exception as err: From 03eb20a5613c248bbfbb566ba5a854f0a7fff771 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 26 Apr 2024 10:54:56 -0700 Subject: [PATCH 37/48] In 'normalized_params', use 'None' values to produce default values. Added tests also. --- holodeck/librarian/libraries.py | 23 ++++++--- .../tests/test_libraries__param_space.py | 51 ++++++++++++++++++- 2 files changed, 65 insertions(+), 9 deletions(-) diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/libraries.py index 8fbabfcb..9fee6f73 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/libraries.py @@ -397,22 +397,23 @@ def model_for_sample_number(self, samp_num, sam_shape=None): return self.model_for_params(params, sam_shape) def normalized_params(self, vals): - """Convert input values (uniform/linear) into parameters from the stored distributions. + """Convert input values (uniform [0.0, 1.0]) into parameters from the stored distributions. For example, if this parameter space has 2 dimensions, where the distributions are: 0. 'value_a' is a uniform parameter from [-1.0, 1.0], and - 1. 'value_b' normal with mean 10.0 and stdev 1.0 + 1. 'value_b' is normally distributed with mean 10.0 and stdev 1.0 Then input values of ``[0.75, 0.5]`` are mapped to parameters ``[0.5, 10.0]``, which will be returned as ``{value_a: 0.5, value_b: 10.0}``. Arguments --------- - vals : (P,) iterable of float, - A list/iterable of `P` float values, matching the number of parameters (i.e. dimensions) - in this parameter space. Each value is passed to the corresponding distribution for - that parameter. + vals : (P,) iterable of (float or `None`), + A list/iterable of `P` values, matching the number of parameters (i.e. dimensions) + in this parameter space. If a value is `None`, then the default value for that + parameter is obtained, otherwise the value is passed to the corresponding parameter + distribution. Returns ------- @@ -424,11 +425,17 @@ def normalized_params(self, vals): if np.ndim(vals) == 0: vals = self.nparameters * [vals] assert len(vals) == self.nparameters + all_finite = np.all([(vv is None) or np.isfinite(vv) for vv in vals]) + assert all_finite, f"Not all `vals` are finite! {vals}" params = {} for ii, pname in enumerate(self.param_names): - vv = vals[ii] # desired fractional parameter value [0.0, 1.0] - ss = self._parameters[ii](vv) # convert to actual parameter values + param = self._parameters[ii] # this parameter distribution + vv = vals[ii] # fractional parameter value [0.0, 1.0] or `None` + if vv is None: + ss = param.default + else: + ss = param(vv) # convert from fractional to actual values params[pname] = ss # store to dictionary return params diff --git a/holodeck/librarian/tests/test_libraries__param_space.py b/holodeck/librarian/tests/test_libraries__param_space.py index aa7c54bd..57ef580f 100644 --- a/holodeck/librarian/tests/test_libraries__param_space.py +++ b/holodeck/librarian/tests/test_libraries__param_space.py @@ -9,9 +9,11 @@ from holodeck import sams, host_relations, hardening, librarian from holodeck.constants import GYR +from holodeck.librarian import ( + param_spaces_dict +) from holodeck.librarian.libraries import ( _Param_Space, PD_Uniform, - ) MMB_MAMP_LOG10_EXTR = [+7.5, +9.5] @@ -243,3 +245,50 @@ def test__ps_test_with_defaults(): _check_sam_hard(sam, hard, SAM_SHAPE) return + + +def _check_ps_basics(name, pspace_class): + # all default arguments + space = pspace_class() + + param_names = space.param_names + params = space._parameters + nparams = space.nparameters + assert len(param_names) == nparams + assert len(params) == nparams + + # check defaults + def_params = space.default_params() + for ii, (name, dp) in enumerate(def_params.items()): + assert dp == params[ii].default + + np.random.seed(12345) + for ii in range(nparams): + # Choose random normalized values for each parameter + # must convert to list so that we can replace values with `None` + norm_vals = np.random.uniform(0.0, 1.0, size=nparams).tolist() + # Set parameter `ii` to default value (`None`) + norm_vals[ii] = None + dist_vals = space.normalized_params(norm_vals) + for jj, (pnam, pval) in enumerate(dist_vals.items()): + # make sure default value is set + if jj == ii: + def_val = params[jj].default + err = f"default was not set for {jj}: `{pnam}`! value={pval}, def={def_val}" + assert pval == def_val, err + # make sure other values are set + else: + truth = params[jj](norm_vals[jj]) + err = f"normalized param was not set for {jj}: `{pnam}`! value={pval} {truth=}" + assert pval == truth, err + + return + +def test_basics_all_param_spaces(): + """Run `_check_ps_basics` on all registered parameters-space classes. + """ + for name, pspace_class in param_spaces_dict.items(): + print(f"Checking '{name}' ({pspace_class})") + _check_ps_basics(name, pspace_class) + + return \ No newline at end of file From 7892313bcffb8fa8d33e828959ac92651a387c23 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 26 Apr 2024 12:16:51 -0700 Subject: [PATCH 38/48] Pass existing 'log' to sam instances. Fix issue with bf_frac_lo extending to 0.0 (cant) --- holodeck/librarian/param_spaces.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/holodeck/librarian/param_spaces.py b/holodeck/librarian/param_spaces.py index 53b13a3b..cff7f5d3 100644 --- a/holodeck/librarian/param_spaces.py +++ b/holodeck/librarian/param_spaces.py @@ -129,7 +129,7 @@ def _init_sam(self, sam_shape, params): ) sam = holo.sams.Semi_Analytic_Model( - gsmf=gsmf, gmr=gmr, mmbulge=mmbulge, shape=sam_shape, + gsmf=gsmf, gmr=gmr, mmbulge=mmbulge, shape=sam_shape, log=self._log, ) return sam @@ -206,8 +206,7 @@ class _PS_Astro_Strong(_Param_Space): bf_width_dex=1.0, # [dex] ) - @classmethod - def _init_sam(cls, sam_shape, params): + def _init_sam(self, sam_shape, params): log10_phi_one = [ params['gsmf_log10_phi_one_z0'], params['gsmf_log10_phi_one_z1'], @@ -259,12 +258,11 @@ def _init_sam(cls, sam_shape, params): ) sam = holo.sams.Semi_Analytic_Model( - gsmf=gsmf, gmr=gmr, mmbulge=mmbulge, shape=sam_shape, + gsmf=gsmf, gmr=gmr, mmbulge=mmbulge, shape=sam_shape, log=self._log, ) return sam - @classmethod - def _init_hard(cls, sam, params): + def _init_hard(self, sam, params): hard = holo.hardening.Fixed_Time_2PL_SAM( sam, params['hard_time']*GYR, @@ -315,7 +313,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): PD_Normal('mmb_plaw', 1.17, 0.08), # 1.17 ± 0.08 # Extra PD_Normal('mmb_scatter_dex', 0.28, 0.05), # no uncertainties given - PD_Uniform('bf_frac_lo', 0.0, 0.4), + PD_Uniform('bf_frac_lo', 0.1, 0.4), PD_Uniform('bf_frac_hi', 0.6, 1.0), PD_Uniform('bf_width_dex', 0.5, 1.5), # [dex] ] @@ -414,7 +412,7 @@ def __init__(self, log=None, nsamples=None, sam_shape=None, seed=None): PD_Normal('mmb_mamp_log10', 8.69, 0.05), # 8.69 ± 0.05 [log10(M/Msol)] PD_Normal('mmb_plaw', 1.17, 0.08), # 1.17 ± 0.08 PD_Normal('mmb_scatter_dex', 0.28, 0.05), # no uncertainties given - PD_Uniform('bf_frac_lo', 0.0, 0.4), + PD_Uniform('bf_frac_lo', 0.1, 0.4), PD_Uniform('bf_frac_hi', 0.6, 1.0), PD_Uniform('bf_width_dex', 0.5, 1.5), # [dex] ] From ae11d212ee625e231ee547e68d2ef144c266b70b Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Fri, 26 Apr 2024 12:17:09 -0700 Subject: [PATCH 39/48] Full (trial) implementation of domain construction with 'combine' working. --- holodeck/librarian/__init__.py | 17 +++++- holodeck/librarian/combine.py | 101 +++++++++++++++++++++++++------- holodeck/librarian/gen_lib.py | 93 +++++++++++++++++++++-------- holodeck/librarian/libraries.py | 28 ++++++--- 4 files changed, 183 insertions(+), 56 deletions(-) diff --git a/holodeck/librarian/__init__.py b/holodeck/librarian/__init__.py index f75dc1e6..56f40563 100644 --- a/holodeck/librarian/__init__.py +++ b/holodeck/librarian/__init__.py @@ -49,10 +49,22 @@ FITS_NBINS_PLAW = [3, 4, 5, 10, 15] FITS_NBINS_TURN = [5, 10, 15] -FNAME_SIM_FILE = "sam-lib__p{pnum:06d}.npz" +FNAME_LIBRARY_SIM_FILE = "library__p{pnum:06d}.npz" +FNAME_DOMAIN_SIM_FILE = "domain__p{pnum:06d}.npz" +DIRNAME_LIBRARY_SIMS = "library_sims" +DIRNAME_DOMAIN_SIMS = "domain_sims" PSPACE_FILE_SUFFIX = ".pspace.npz" ARGS_CONFIG_FNAME = "config.json" +#: When constructing parameter-space domains, go this far along each dimension when parameter is +# unbounded. +PSPACE_DOMAIN_EXTREMA = [0.001, 0.999] + +class DomainNotLibraryError(Exception): + def __init__(self, message="This looks like a 'domain' not a 'library'!"): + # Call the base class constructor with the parameters it needs + super(DomainNotLibraryError, self).__init__(message) + from . libraries import ( # noqa _Param_Space, _Param_Dist, PD_Uniform, PD_Normal, @@ -62,6 +74,7 @@ param_spaces_dict = {} #: Registry of standard parameter-spaces from . import param_spaces_classic as psc_temp # noqa param_spaces_dict.update(psc_temp._param_spaces_dict) +del psc_temp from . import param_spaces as ps_temp # noqa param_spaces_dict.update(ps_temp._param_spaces_dict) -del psc_temp, ps_temp +del ps_temp diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index a3b2f968..202fa369 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -19,7 +19,10 @@ import holodeck as holo import holodeck.librarian -from holodeck.librarian import libraries +import holodeck.librarian.gen_lib +from holodeck.librarian import ( + libraries, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS, DomainNotLibraryError +) def main(): @@ -71,12 +74,19 @@ def main(): # ---- Combine library files - sam_lib_combine(path, log, recreate=args.recreate, gwb_only=args.gwb) + for library in [True, False]: + try: + sam_lib_combine(path, log, recreate=args.recreate, gwb_only=args.gwb, library=library) + except DomainNotLibraryError as err: + pass return -def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only=False): +def sam_lib_combine( + path_output, log, + path_pspace=None, recreate=False, gwb_only=False, library=True +): """Combine individual simulation files into a single library (hdf5) file. Arguments @@ -91,6 +101,8 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only path_pspace : str or None, Path to file containing _Param_Space subclass instance. If `None` then `path_output` is searched for a `_Param_Space` save file. + library : bool, + Combine simulations from a library, as opposed to simulations for 'domain' exploration. Returns ------- @@ -103,7 +115,10 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only path_output = Path(path_output) log.info(f"Path output = {path_output}") - path_sims = path_output.joinpath('sims') + if library: + path_sims = path_output.joinpath(DIRNAME_LIBRARY_SIMS) + else: + path_sims = path_output.joinpath(DIRNAME_DOMAIN_SIMS) # ---- see if a combined library already exists @@ -121,18 +136,34 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only if path_pspace is None: path_pspace = path_output pspace, pspace_fname = libraries.load_pspace_from_path(path_pspace, log=log) + args, args_fname = holo.librarian.gen_lib.load_config_from_path(path_pspace, log=log) log.info(f"loaded param space: {pspace} from '{pspace_fname}'") param_names = pspace.param_names param_samples = pspace.param_samples - nsamp, ndim = param_samples.shape - log.debug(f"{nsamp=}, {ndim=}, {param_names=}") + + if library: + if param_samples == None: # noqa : use `== None` to match arrays + log.error(f"`library`={library} but `param_samples`={param_samples}`") + err = f"`library` is True, but {path_output} looks like it's a domain." + raise DomainNotLibraryError(err) + + nsamp_all, ndim = param_samples.shape + log.debug(f"{nsamp_all=}, {ndim=}, {param_names=}") + else: + err = f"Expected 'domain' but `param_samples`={param_samples} is not `None`!" + assert param_samples == None, err # noqa : use `== None` to match arrays + ndim = pspace.nparameters + # for 'domain' simulations, this is the number of samples in each dimension + nsamp_dim = args.nsamples + # get the total number of samples + nsamp_all = ndim * nsamp_dim # ---- make sure all files exist; get shape information from files - log.info(f"checking that all {nsamp} files exist") + log.info(f"checking that all {nsamp_all} files exist") fobs_cents, fobs_edges, nreals, nloudest, has_gwb, has_ss, has_params = _check_files_and_load_shapes( - log, path_sims, nsamp + log, path_sims, nsamp_all, library ) nfreqs = fobs_cents.size log.debug(f"{nfreqs=}, {nreals=}, {nloudest=}") @@ -150,28 +181,29 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only # ---- load results from all files - gwb = np.zeros((nsamp, nfreqs, nreals)) if has_gwb else None + gwb = np.zeros((nsamp_all, nfreqs, nreals)) if has_gwb else None if (not gwb_only) and has_ss: - hc_ss = np.zeros((nsamp, nfreqs, nreals, nloudest)) - hc_bg = np.zeros((nsamp, nfreqs, nreals)) + hc_ss = np.zeros((nsamp_all, nfreqs, nreals, nloudest)) + hc_bg = np.zeros((nsamp_all, nfreqs, nreals)) else: hc_ss = None hc_bg = None if (not gwb_only) and has_params: - sspar = np.zeros((nsamp, 4, nfreqs, nreals, nloudest)) - bgpar = np.zeros((nsamp, 7, nfreqs, nreals)) + sspar = np.zeros((nsamp_all, 4, nfreqs, nreals, nloudest)) + bgpar = np.zeros((nsamp_all, 7, nfreqs, nreals)) else: sspar = None bgpar = None gwb, hc_ss, hc_bg, sspar, bgpar, bad_files = _load_library_from_all_files( - path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log, + path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log, library ) if has_gwb: log.info(f"Loaded data from all library files | {holo.utils.stats(gwb)=}") - param_samples[bad_files] = np.nan + if library: + param_samples[bad_files] = np.nan # ---- Save to concatenated output file ---- @@ -179,14 +211,31 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only with h5py.File(lib_path, 'w') as h5: h5.create_dataset('fobs_cents', data=fobs_cents) h5.create_dataset('fobs_edges', data=fobs_edges) - h5.create_dataset('sample_params', data=param_samples) + if param_samples is not None: + h5.create_dataset('sample_params', data=param_samples) if gwb is not None: + # if 'domain', reshape to include each dimension + if not library: + new_shape = (ndim, nsamp_dim) + gwb.shape[1:] + gwb = gwb.reshape(new_shape) h5.create_dataset('gwb', data=gwb) if not gwb_only: if has_ss: + # if 'domain', reshape to include each dimension + if not library: + new_shape = (ndim, nsamp_dim) + hc_ss.shape[1:] + hc_ss = hc_ss.reshape(new_shape) + new_shape = (ndim, nsamp_dim) + hc_bg.shape[1:] + hc_bg = hc_bg.reshape(new_shape) h5.create_dataset('hc_ss', data=hc_ss) h5.create_dataset('hc_bg', data=hc_bg) if has_params: + # if 'domain', reshape to include each dimension + if not library: + new_shape = (ndim, nsamp_dim) + sspar.shape[1:] + sspar = sspar.reshape(new_shape) + new_shape = (ndim, nsamp_dim) + bgpar.shape[1:] + bgpar = bgpar.reshape(new_shape) h5.create_dataset('sspar', data=sspar) h5.create_dataset('bgpar', data=bgpar) h5.attrs['param_names'] = np.array(param_names).astype('S') @@ -211,16 +260,20 @@ def sam_lib_combine(path_output, log, path_pspace=None, recreate=False, gwb_only return lib_path -def _check_files_and_load_shapes(log, path_sims, nsamp): +def _check_files_and_load_shapes(log, path_sims, nsamp, library): """Check that all `nsamp` files exist in the given path, and load info about array shapes. Arguments --------- + log : ``logging.Logger`` instance + Holodeck logger. path_sims : str Path in which individual simulation files can be found. nsamp : int Number of simulations/files that should be found. This should typically be loaded from the parameter-space object used to generate the library. + library : bool, + Whether this is a standard library or a 'domain' exploration. Returns ------- @@ -242,7 +295,7 @@ def _check_files_and_load_shapes(log, path_sims, nsamp): log.info(f"Checking {nsamp} files in {path_sims}") for ii in tqdm.trange(nsamp): - temp_fname = libraries._get_sim_fname(path_sims, ii) + temp_fname = libraries._get_sim_fname(path_sims, ii, library=library) if not temp_fname.exists(): err = f"Missing at least file number {ii} out of {nsamp} files! {temp_fname}" log.exception(err) @@ -301,7 +354,7 @@ def _check_files_and_load_shapes(log, path_sims, nsamp): return fobs_cents, fobs_edges, nreals, nloudest, has_gwb, has_ss, has_params -def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log): +def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log, library): """Load data from all individual simulation files. Arguments @@ -311,8 +364,14 @@ def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log gwb : (S, F, R) ndarray Array in which to store GWB data from all of 'S' files. S: num-samples/simulations, F: num-frequencies, R: num-realizations. - log : `logging.Logger` + hc_ss + hc_bg + sspar + bgpar + log : ``logging.Logger`` Logging instance. + library : bool + Whether this is a standard library or a domain exploration. """ if hc_bg is not None: @@ -328,7 +387,7 @@ def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log bad_files = np.zeros(nsamp, dtype=bool) #: track which files contain UN-useable data msg = None for pnum in tqdm.trange(nsamp): - fname = libraries._get_sim_fname(path_sims, pnum) + fname = libraries._get_sim_fname(path_sims, pnum, library=library) temp = np.load(fname, allow_pickle=True) # When a processor fails for a given parameter, the output file is still created with the 'fail' key added if ('fail' in temp): diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index a7bece01..decb3f60 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -35,7 +35,7 @@ import holodeck.librarian import holodeck.librarian.combine from holodeck.librarian import ( - libraries, ARGS_CONFIG_FNAME + libraries, ARGS_CONFIG_FNAME, PSPACE_DOMAIN_EXTREMA, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS ) #: maximum number of failed simulations before task terminates with error (`None`: no limit) @@ -112,7 +112,7 @@ def main(): # noqa : ignore complexity warning # Load arguments/configuration from previous save if args.resume: - args, config_fname = _load_config(args.output, log) + args, config_fname = load_config_from_path(args.output, log) log.warning(f"Loaded configuration save from {config_fname}") args.resume = True # Save parameter space and args/configuration to output directory @@ -123,24 +123,35 @@ def main(): # noqa : ignore complexity warning config_fname = _save_config(args) log.info(f"saved configuration to {config_fname}") - # Split simulations for all processes + # ---- Split simulations for all processes + # Domain: Vary only one parameter at a time to explore the domain if args.domain: - pass - # + log.info("Constructing domain-exploration indices") + # we will have `nsamples` in each of `nparameters` + domain_shape = (space.nparameters, args.nsamples) + indices = np.arange(np.product(domain_shape)) + # indices = np.random.permutation(indices) + indices = np.array_split(indices, comm.size) + # Standard Library: vary all parameters together else: + log.info("Constructing library indices") indices = range(args.nsamples) indices = np.random.permutation(indices) indices = np.array_split(indices, comm.size) - num_per = [len(ii) for ii in indices] - log.info(f"{args.nsamples=} cores={comm.size} || max sims per core = {np.max(num_per)}") + domain_shape = None + + num_per = [len(ii) for ii in indices] + log.info(f"{args.nsamples=} cores={comm.size} || max sims per core = {np.max(num_per)}") else: space = None indices = None + domain_shape = None # share parameter space across processes space = comm.bcast(space, root=0) + domain_shape = comm.bcast(domain_shape, root=0) # If we've loaded a new `args`, then share to all processes from rank=0 if args.resume: @@ -148,7 +159,8 @@ def main(): # noqa : ignore complexity warning args.log = log log.info( - f"param_space={args.param_space}, samples={args.nsamples}, sam_shape={args.sam_shape}, nreals={args.nreals}\n" + f"param_space={args.param_space}, parameters={space.nparameters}, samples={args.nsamples}\n" + f"sam_shape={args.sam_shape}, nreals={args.nreals}\n" f"nfreqs={args.nfreqs}, pta_dur={args.pta_dur} [yr]\n" ) @@ -165,18 +177,35 @@ def main(): # noqa : ignore complexity warning log.info(f"beginning tasks at {beg}") failures = 0 num_done = 0 - for par_num in iterator: + for sim_num in iterator: + log.info(f"{comm.rank=} {sim_num=}") - log.info(f"{comm.rank=} {par_num=}") - pdict = space.param_dict(par_num) - msg = [] - for kk, vv in pdict.items(): - msg.append(f"{kk}={vv:.4e}") - msg = ", ".join(msg) - log.info(msg) + # Domain: Vary only one parameter at a time to explore the domain + if args.domain: + param_num, samp_num = np.unravel_index(sim_num, domain_shape) + # start out with all default parameters (signified by `None` values) + norm_params = [None for ii in range(space.nparameters)] + # determine the extrema for this parameter. If parameter-distribution is bounded, use + # full range. If unbounded, use hardcoded values. + yext = space._parameters[param_num].extrema #: this is the extrema of the range + xext = [0.0, 1.0] #: this is the extrema of the domain + for ii in range(2): + xext[ii] = xext[ii] if np.isfinite(yext[ii]) else PSPACE_DOMAIN_EXTREMA[ii] + assert (0.0 <= xext[ii]) and (xext[ii] <= 1.0) + # replace the parameter being varied with a fractional value based on the sample number + norm_params[param_num] = xext[0] + samp_num * np.diff(xext)[0] / (domain_shape[1] - 1) + params = space.normalized_params(norm_params) + + # Library: vary all parameters together + else: + params = space.param_dict(sim_num) + msg = [] + for kk, vv in params.items(): + msg.append(f"{kk}={vv:.4e}") + msg = ", ".join(msg) + log.info(msg) - params = space.param_dict(par_num) - rv, _sim_fname = run_sam_at_pspace_params(args, space, par_num, params) + rv, _sim_fname = run_sam_at_pspace_params(args, space, sim_num, params) if rv is False: failures += 1 @@ -198,7 +227,7 @@ def main(): # noqa : ignore complexity warning if (comm.rank == 0): log.info("Concatenating outputs into single file") - holo.librarian.combine.sam_lib_combine(args.output, log) + holo.librarian.combine.sam_lib_combine(args.output, log, library=(not args.domain)) log.info("Concatenation completed") return @@ -249,7 +278,8 @@ def run_sam_at_pspace_params(args, space, pnum, params): # ---- get output filename for this simulation, check if already exists - sim_fname = libraries._get_sim_fname(args.output_sims, pnum) + library_flag = not args.domain + sim_fname = libraries._get_sim_fname(args.output_sims, pnum, library=library_flag) beg = datetime.now() log.info(f"{pnum=} :: {sim_fname=} beginning at {beg}") @@ -343,8 +373,10 @@ def _setup_argparse(*args, **kwargs): help="calculate and store SS/CW sources and the BG separately") parser.add_argument('--params', dest="params_flag", default=True, action=argparse.BooleanOptionalAction, help="calculate and store SS/BG binary parameters [NOTE: requires `--ss`]") + parser.add_argument('--domain', default=False, action='store_true', + help="instead of generating a standard library, explore each parameter.") - # how do run + # how to run parser.add_argument('--resume', action='store_true', default=False, help='resume production of a library by loading previous parameter-space from output directory') parser.add_argument('--recreate', action='store_true', default=False, @@ -417,7 +449,12 @@ def _setup_argparse(*args, **kwargs): holo.utils.mpi_print(f"output path: {output}") args.output = output - output_sims = output.joinpath("sims") + if args.domain: + sims_dirname = DIRNAME_DOMAIN_SIMS + else: + sims_dirname = DIRNAME_LIBRARY_SIMS + + output_sims = output.joinpath(sims_dirname) output_sims.mkdir(parents=True, exist_ok=True) args.output_sims = output_sims @@ -434,6 +471,12 @@ def _setup_argparse(*args, **kwargs): def _setup_param_space(args): + """Setup the parameter-space instance. + + For normal runs, identify the parameter-space class and construct a new class instance. + For 'resume' runs, load a saved parameter-space instance. + + """ log = args.log # ---- Determine and load the parameter-space class @@ -466,7 +509,9 @@ def _setup_param_space(args): space, space_fname = holo.librarian.load_pspace_from_path(args.output, space_class=space_class, log=log) log.warning(f"Loaded param-space save from {space_fname}") else: - space = space_class(log, args.nsamples, args.sam_shape, args.seed) + # we don't use standard samples when constructing a parameter-space 'domain' + nsamples = None if args.domain else args.nsamples + space = space_class(log, nsamples, args.sam_shape, args.seed) return space @@ -501,7 +546,7 @@ def _save_config(args): return fname -def _load_config(path, log): +def load_config_from_path(path, log): fname = Path(path).joinpath(ARGS_CONFIG_FNAME) with open(fname, 'r') as inp: diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/libraries.py index 9fee6f73..2d0f2682 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/libraries.py @@ -123,7 +123,7 @@ def __init__(self, parameters, log=None, nsamples=None, sam_shape=None, seed=Non param_names.append(name) if (nsamples is None) or (nparameters == 0): - log.warning(f"{self}: {nsamples=} {nparameters=} - cannot generate parameter samples.") + log.info(f"{self}: {nsamples=} {nparameters=} - cannot generate parameter samples.") uniform_samples = None param_samples = None else: @@ -206,9 +206,9 @@ def model_for_params(self, params, sam_shape=None): return sam, hard - @classmethod + # @classmethod @abc.abstractmethod - def _init_sam(cls, sam_shape, params): + def _init_sam(self, sam_shape, params): """Initialize a :class:`holodeck.sams.sam.Semi_Analytic_Model` instance with given params. Arguments @@ -229,9 +229,9 @@ def _init_sam(cls, sam_shape, params): """ raise - @classmethod + # @classmethod @abc.abstractmethod - def _init_hard(cls, sam, params): + def _init_hard(self, sam, params): """Initialize a :class:`holodeck.hardening._Hardening` subclass instance with given params. Arguments @@ -324,8 +324,13 @@ def from_save(cls, fname, log=None): pspace_class = cls # construct instance with dummy/temporary values (which will be overwritten) - nsamples = data['param_samples'].shape[0] - nparameters = data['param_samples'].shape[1] + if data['param_samples'] == None: # noqa : use ``== None`` to match arrays + nsamples = None + nparameters = None + else: + # print(f"{data['param_samples']=}") + nsamples = data['param_samples'].shape[0] + nparameters = data['param_samples'].shape[1] param_names = data['param_names'] space = pspace_class(nsamples=nsamples, log=log) if class_name != space.name: @@ -1005,8 +1010,13 @@ def _get_space_class_from_space_fname(space_fname): return space_class -def _get_sim_fname(path, pnum): - temp = holo.librarian.FNAME_SIM_FILE.format(pnum=pnum) +def _get_sim_fname(path, pnum, library=True): + if library: + temp = holo.librarian.FNAME_LIBRARY_SIM_FILE + else: + temp = holo.librarian.FNAME_DOMAIN_SIM_FILE + + temp = temp.format(pnum=pnum) temp = path.joinpath(temp) return temp From 860ff0dd83c8903ffd18a8169486621e50cc1208 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 29 Apr 2024 10:20:49 -0700 Subject: [PATCH 40/48] Add separate names for combined simulation 'library' and 'domain' --- CHANGELOG.md | 3 +++ holodeck/librarian/__init__.py | 2 ++ holodeck/librarian/libraries.py | 11 +++++++++-- 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c434d11c..2d4387d6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,6 +29,9 @@ * Semi-Analytic Models * Improve accuracy of dynamic binary number calculation for consistent evolution models. +* `holodeck.librarian` + * Added functionality to construct 'domain' sets of simulations, to explore each parameter in a parameter-space one at a time. + * NOTE: Standard library files will now be called "sam-library.hdf5" instead of "sam_lib.hdf5" ## v1.5.2 - 2024/04/12 diff --git a/holodeck/librarian/__init__.py b/holodeck/librarian/__init__.py index 56f40563..427bfdef 100644 --- a/holodeck/librarian/__init__.py +++ b/holodeck/librarian/__init__.py @@ -53,6 +53,8 @@ FNAME_DOMAIN_SIM_FILE = "domain__p{pnum:06d}.npz" DIRNAME_LIBRARY_SIMS = "library_sims" DIRNAME_DOMAIN_SIMS = "domain_sims" +FNAME_LIBRARY_COMBINED_FILE = "sam-library" # do NOT include file suffix (i.e. 'hdf5') +FNAME_DOMAIN_COMBINED_FILE = "sam-domain" # do NOT include file suffix (i.e. 'hdf5') PSPACE_FILE_SUFFIX = ".pspace.npz" ARGS_CONFIG_FNAME = "config.json" diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/libraries.py index 2d0f2682..b9261495 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/libraries.py @@ -15,6 +15,7 @@ from holodeck.constants import YR from holodeck.librarian import ( DEF_NUM_FBINS, DEF_NUM_LOUDEST, DEF_NUM_REALS, DEF_PTA_DUR, + FNAME_LIBRARY_COMBINED_FILE, FNAME_DOMAIN_COMBINED_FILE, ) PARAM_NAMES__ERROR = [ @@ -1021,8 +1022,14 @@ def _get_sim_fname(path, pnum, library=True): return temp -def get_sam_lib_fname(path, gwb_only): - fname = 'sam_lib' +def get_sam_lib_fname(path, gwb_only, library=True): + # standard 'library' + if library: + fname = FNAME_LIBRARY_COMBINED_FILE + # 'domain' of parameter space + else: + fname = FNAME_DOMAIN_COMBINED_FILE + if gwb_only: fname += "_gwb-only" lib_path = path.joinpath(fname).with_suffix(".hdf5") From 8b6d60b381afb039bdfd644964a634bea94955f0 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 29 Apr 2024 10:58:08 -0700 Subject: [PATCH 41/48] Store parameters in each simulation file. Load param-samples when combining 'domain' runs, and save to combined domain file. --- holodeck/librarian/combine.py | 36 ++++++++++++++++++++++++----------- holodeck/librarian/gen_lib.py | 2 ++ 2 files changed, 27 insertions(+), 11 deletions(-) diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index 202fa369..ede1a380 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -142,6 +142,7 @@ def sam_lib_combine( param_names = pspace.param_names param_samples = pspace.param_samples + # Standard Library: vary all parameters together if library: if param_samples == None: # noqa : use `== None` to match arrays log.error(f"`library`={library} but `param_samples`={param_samples}`") @@ -150,6 +151,8 @@ def sam_lib_combine( nsamp_all, ndim = param_samples.shape log.debug(f"{nsamp_all=}, {ndim=}, {param_names=}") + + # Domain: Vary only one parameter at a time to explore the domain else: err = f"Expected 'domain' but `param_samples`={param_samples} is not `None`!" assert param_samples == None, err # noqa : use `== None` to match arrays @@ -158,6 +161,9 @@ def sam_lib_combine( nsamp_dim = args.nsamples # get the total number of samples nsamp_all = ndim * nsamp_dim + # for 'domain', param_samples will eventually be shaped (ndim, nsamp_dim, ndim), but load + # the data first as `(nsamp_all, ndim)`, then we will reshape. + param_samples = np.zeros((nsamp_all, ndim)) # ---- make sure all files exist; get shape information from files @@ -197,8 +203,8 @@ def sam_lib_combine( sspar = None bgpar = None - gwb, hc_ss, hc_bg, sspar, bgpar, bad_files = _load_library_from_all_files( - path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log, library + gwb, hc_ss, hc_bg, sspar, bgpar, param_samples, bad_files = _load_library_from_all_files( + path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, param_samples, log, library ) if has_gwb: log.info(f"Loaded data from all library files | {holo.utils.stats(gwb)=}") @@ -211,8 +217,7 @@ def sam_lib_combine( with h5py.File(lib_path, 'w') as h5: h5.create_dataset('fobs_cents', data=fobs_cents) h5.create_dataset('fobs_edges', data=fobs_edges) - if param_samples is not None: - h5.create_dataset('sample_params', data=param_samples) + h5.create_dataset('sample_params', data=param_samples) if gwb is not None: # if 'domain', reshape to include each dimension if not library: @@ -354,7 +359,10 @@ def _check_files_and_load_shapes(log, path_sims, nsamp, library): return fobs_cents, fobs_edges, nreals, nloudest, has_gwb, has_ss, has_params -def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log, library): +def _load_library_from_all_files( + path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, param_samples, + log, library, +): """Load data from all individual simulation files. Arguments @@ -368,6 +376,7 @@ def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log hc_bg sspar bgpar + param_samples log : ``logging.Logger`` Logging instance. library : bool @@ -375,18 +384,18 @@ def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log """ if hc_bg is not None: - nsamp = hc_bg.shape[0] + nsamp_all = hc_bg.shape[0] elif gwb is not None: - nsamp = gwb.shape[0] + nsamp_all = gwb.shape[0] else: err = "Unable to get shape from either `hc_bg` or `gwb`!" log.exception(err) raise RuntimeError(err) - log.info(f"Collecting data from {nsamp} files") - bad_files = np.zeros(nsamp, dtype=bool) #: track which files contain UN-useable data + log.info(f"Collecting data from {nsamp_all} files") + bad_files = np.zeros(nsamp_all, dtype=bool) #: track which files contain UN-useable data msg = None - for pnum in tqdm.trange(nsamp): + for pnum in tqdm.trange(nsamp_all): fname = libraries._get_sim_fname(path_sims, pnum, library=library) temp = np.load(fname, allow_pickle=True) # When a processor fails for a given parameter, the output file is still created with the 'fail' key added @@ -404,6 +413,11 @@ def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log bad_files[pnum] = True continue + # for 'domain' simulations, we need to load the parameters. For 'library' runs, we + # already have them. + if not library: + param_samples[pnum, :] = temp['params'][:] + # store the GWB from this file if gwb is not None: gwb[pnum, :, :] = temp['gwb'][...] @@ -418,7 +432,7 @@ def _load_library_from_all_files(path_sims, gwb, hc_ss, hc_bg, sspar, bgpar, log log.info(f"{holo.utils.frac_str(bad_files)} files are failures") - return gwb, hc_ss, hc_bg, sspar, bgpar, bad_files + return gwb, hc_ss, hc_bg, sspar, bgpar, param_samples, bad_files if __name__ == "__main__": diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index decb3f60..48e8814d 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -307,6 +307,8 @@ def run_sam_at_pspace_params(args, space, pnum, params): gwb_flag=args.gwb_flag, singles_flag=args.ss_flag, details_flag=False, params_flag=args.params_flag, log=log, ) + data['params'] = params + data['param_names'] = space.param_names rv = True except Exception as err: From 16a8bcf541a25733314e8e81c3c2c25c33e8be44 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Mon, 29 Apr 2024 11:22:23 -0700 Subject: [PATCH 42/48] Deprecate 'holodeck.librarian.libraries' ==> 'holodeck.librarian.lib_tools' --- CHANGELOG.md | 23 ++++++++----- ...s.rst => holodeck.librarian.lib_tools.rst} | 4 +-- docs/source/api_ref/holodeck.librarian.rst | 2 +- docs/source/getting_started/libraries.rst | 16 +++++----- docs/source/header.rst | 4 +-- holodeck/librarian/__init__.py | 31 ++++++++++++------ holodeck/librarian/combine.py | 10 +++--- holodeck/librarian/fit_spectra.py | 6 ++-- holodeck/librarian/gen_lib.py | 14 ++++---- .../librarian/{libraries.py => lib_tools.py} | 6 +++- holodeck/librarian/param_spaces.py | 2 +- holodeck/librarian/param_spaces_classic.py | 2 +- ...pace.py => test_lib_tools__param_space.py} | 8 ++--- holodeck/librarian/tests/test_param_spaces.py | 4 +-- notebooks/devs/libraries/gwb_anatomy.ipynb | 8 ++--- .../devs/libraries/library-explorer.ipynb | 4 +-- notebooks/devs/libraries/param_spaces.ipynb | 3 +- notebooks/librarian.ipynb | 32 +++++++++---------- 18 files changed, 100 insertions(+), 79 deletions(-) rename docs/source/api_ref/{holodeck.librarian.libraries.rst => holodeck.librarian.lib_tools.rst} (69%) rename holodeck/librarian/{libraries.py => lib_tools.py} (99%) rename holodeck/librarian/tests/{test_libraries__param_space.py => test_lib_tools__param_space.py} (97%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2d4387d6..d6b66d66 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,19 +19,26 @@ ## Current -* DEFAULTS: - * Semi-Analytic Models - * `Semi_Analytic_Models` instances now use galaxy merger-rates (instead of galaxy pair-fractions and merger-times) by default. To use GPF+GMT SAMs, the user must pass in at least a GPF instance manually. +* DEPRECATIONS: + * `holodeck.librarian.py` ==> `holodeck.librarian.lib_tools` + * Rename submodule. All components remain the same. All `lib_tools` elements are now also imported into the `librarian` namespace. i.e. elements like `holodeck.librarian.lib_tools._Param_Space` will now be accessible via `holodeck.librarian._Param_Space`. * BUGS: * Semi-Analytic Models * `Semi_Analytic_Model._dynamic_binary_number_at_fobs_inconsistent` was not checking for systems that had already coalesced. Only effected GW-only evolution using the python version of the calculation. -* Semi-Analytic Models - * Improve accuracy of dynamic binary number calculation for consistent evolution models. -* `holodeck.librarian` - * Added functionality to construct 'domain' sets of simulations, to explore each parameter in a parameter-space one at a time. - * NOTE: Standard library files will now be called "sam-library.hdf5" instead of "sam_lib.hdf5" +* DEFAULTS: + * Semi-Analytic Models + * `Semi_Analytic_Models` instances now use galaxy merger-rates (instead of galaxy pair-fractions and merger-times) by default. To use GPF+GMT SAMs, the user must pass in at least a GPF instance manually. + +* General Changes + + * Semi-Analytic Models (`holodeck.sams`) + * Improve accuracy of dynamic binary number calculation for consistent evolution models. + + * `holodeck.librarian` + * Added functionality to construct 'domain' sets of simulations, to explore each parameter in a parameter-space one at a time. + * NOTE: Standard library files will now be called "sam-library.hdf5" instead of "sam_lib.hdf5" ## v1.5.2 - 2024/04/12 diff --git a/docs/source/api_ref/holodeck.librarian.libraries.rst b/docs/source/api_ref/holodeck.librarian.lib_tools.rst similarity index 69% rename from docs/source/api_ref/holodeck.librarian.libraries.rst rename to docs/source/api_ref/holodeck.librarian.lib_tools.rst index 69136318..160b12fe 100644 --- a/docs/source/api_ref/holodeck.librarian.libraries.rst +++ b/docs/source/api_ref/holodeck.librarian.lib_tools.rst @@ -1,8 +1,8 @@ ============================ -holodeck.librarian.libraries +holodeck.librarian.lib_tools ============================ -.. automodule:: holodeck.librarian.libraries +.. automodule:: holodeck.librarian.lib_tools :members: :private-members: :undoc-members: diff --git a/docs/source/api_ref/holodeck.librarian.rst b/docs/source/api_ref/holodeck.librarian.rst index 0ea4e99c..9ff9eb56 100644 --- a/docs/source/api_ref/holodeck.librarian.rst +++ b/docs/source/api_ref/holodeck.librarian.rst @@ -15,7 +15,7 @@ holodeck.librarian module holodeck.librarian.combine holodeck.librarian.fit_spectra holodeck.librarian.gen_lib - holodeck.librarian.libraries + holodeck.librarian.lib_tools holodeck.librarian.param_spaces holodeck.librarian.param_spaces_classic holodeck.librarian.posterior_populations diff --git a/docs/source/getting_started/libraries.rst b/docs/source/getting_started/libraries.rst index 235e8ec2..1ea61f15 100644 --- a/docs/source/getting_started/libraries.rst +++ b/docs/source/getting_started/libraries.rst @@ -12,7 +12,7 @@ Generating and Using Holodeck Libraries Libraries Overview ================== -|holodeck| 'libraries' are collections of simulations in which a certain set of parameters are varied, producing different populations and/or GW signatures at each sampled parameter value. Libraries are run from the same parameter-space and using the same hyper parameters (for example, the functional form that is assumed for the galaxy stellar-mass function). Libraries are constructed using the :py:mod:`~holodeck.librarian` module, with a 'parameter space' class that organizes the different simulations. The base-class is called :py:class:`~holodeck.librarian.libraries._Param_Space` (defined in the :py:mod:`holodeck.librarian.libraries` file), and all parameter space classes inherit from this, and should typically be prefixed by ``PS_`` to denote that they are parameter spaces. The parameter-space subclasses implement a number of parameters that are varied. Each parameter is implemented as a subclass of :py:class:`~holodeck.librarian.libraries._Param_Dist`, for example the :py:class:`~holodeck.librarian.libraries.PD_Uniform` class that implements a uniform distribution. +|holodeck| 'libraries' are collections of simulations in which a certain set of parameters are varied, producing different populations and/or GW signatures at each sampled parameter value. Libraries are run from the same parameter-space and using the same hyper parameters (for example, the functional form that is assumed for the galaxy stellar-mass function). Libraries are constructed using the :py:mod:`~holodeck.librarian` module, with a 'parameter space' class that organizes the different simulations. The base-class is called :py:class:`~holodeck.librarian.lib_tools._Param_Space` (defined in the :py:mod:`holodeck.librarian.lib_tools` file), and all parameter space classes inherit from this, and should typically be prefixed by ``PS_`` to denote that they are parameter spaces. The parameter-space subclasses implement a number of parameters that are varied. Each parameter is implemented as a subclass of :py:class:`~holodeck.librarian.lib_tools._Param_Dist`, for example the :py:class:`~holodeck.librarian.lib_tools.PD_Uniform` class that implements a uniform distribution. As an example, the fiducial library and parameter space for :doc:`the 15yr astrophysics analysis ` was the 'phenomenological uniform' library, implemented as :py:class:`~holodeck.librarian.param_spaces_classic.PS_Classic_Phenom_Uniform` (at the time, it was internally called ``PS_Uniform_09B``). This library spanned a 6D parameter space using a 'phenomenological' binary evolution model, and assuming a uniform distribution in sampling from the parameter priors. Two parameters from the galaxy stellar-mass function were varied, along with two parameters from the M-MBulge relationship, and two parameters from the hardening model. @@ -21,14 +21,14 @@ Parameter Spaces and Distributions ================================== **NOTE: currently parameter-spaces are only designed for use with SAMs.** -Parameter spaces are implemented as subclasses of the |pspace_class| class, and are generally named with a ``PS_`` prefix. Each class generates a certain number of samples using a latin hypercube to efficiently sample the parameter space. Each parameter being varied in the parameter space corresponds to parameter distribution, implemented as a :py:class:`~holodeck.librarian.libraries._Param_Dist` subclass. Each subclass is generally named with a ``PD_`` prefix. These parameter distributions convert from uniform random variables (uniform samples in $[0.0, 1.0]$ in the latin hypercube) to the desired distributions. For example, the :py:class:`~holodeck.librarian.libraries.PD_Normal(mean, stdev)` class draws from a normal (Gaussian) distribution, and the :py:class:`~holodeck.librarian.libraries.PD_Normal(min, max)` class draws from a uniform distribution. +Parameter spaces are implemented as subclasses of the |pspace_class| class, and are generally named with a ``PS_`` prefix. Each class generates a certain number of samples using a latin hypercube to efficiently sample the parameter space. Each parameter being varied in the parameter space corresponds to parameter distribution, implemented as a :py:class:`~holodeck.librarian.lib_tools._Param_Dist` subclass. Each subclass is generally named with a ``PD_`` prefix. These parameter distributions convert from uniform random variables (uniform samples in $[0.0, 1.0]$ in the latin hypercube) to the desired distributions. For example, the :py:class:`~holodeck.librarian.lib_tools.PD_Normal(mean, stdev)` class draws from a normal (Gaussian) distribution, and the :py:class:`~holodeck.librarian.lib_tools.PD_Normal(min, max)` class draws from a uniform distribution. Parameter Distributions ----------------------- -New parameter distributions should subclass :py:class:`~holodeck.librarian.libraries._Param_Dist`, and must provide a method with signature: ``_dist_func(self, xx)`` which accepts a float value ``xx`` in $[0.0, 1.0]$ and maps it to a value in the desired distribution, and returns a float value. Typically an ``__init__`` function will also be provided to set any required parameters. See the :py:class:`~holodeck.librarian.libraries.PD_Uniform` class for a simple example that maps from $[0.0, 1.0]$ to another uniform variable with a different minimum (``lo``) and maximum (``hi``) value. +New parameter distributions should subclass :py:class:`~holodeck.librarian.lib_tools._Param_Dist`, and must provide a method with signature: ``_dist_func(self, xx)`` which accepts a float value ``xx`` in $[0.0, 1.0]$ and maps it to a value in the desired distribution, and returns a float value. Typically an ``__init__`` function will also be provided to set any required parameters. See the :py:class:`~holodeck.librarian.lib_tools.PD_Uniform` class for a simple example that maps from $[0.0, 1.0]$ to another uniform variable with a different minimum (``lo``) and maximum (``hi``) value. -How the parameter distributions are used in parameter spaces is described below, but in summary, each |pspace_class| subclass will build a list of |pdist_class| subclass instances which are used to specify the domain of the parameter space. The construct for each |pdist_class| subclass must accept first the variable name, and then any additional required arguments, for example: ``PD_Normal("gsmf_phi0", -2.56, 0.4)``. The name of the variable **must match the name used in the |pspace_class|**, i.e. for the previous example, the |pspace_class| will be expecting a variable named ``gsmf_phi0``. All |pdist_class| subclasses optionally accept a ``default=`` keyword-argument, for example, ``PD_Uniform("hard_time", 0.1, 11.0, default=3.0)``. The 'default' values are provided so that |pspace_class|'s can construct a model using default parameters (see: :py:meth:`holodeck.librarian.libraries._Param_Space.default_params`), typically as a fiducial model. In the preceding example, the default 'hard_time' parameter would be 3.0. If a ``default`` is not specified in the instance constructor, then the value produced by an input of ``0.5`` is used. In the preceding example, if no ``default`` was specified, then the middle value of $(11.0 + 0.1) / 2 = 5.55$ would be used. +How the parameter distributions are used in parameter spaces is described below, but in summary, each |pspace_class| subclass will build a list of |pdist_class| subclass instances which are used to specify the domain of the parameter space. The construct for each |pdist_class| subclass must accept first the variable name, and then any additional required arguments, for example: ``PD_Normal("gsmf_phi0", -2.56, 0.4)``. The name of the variable **must match the name used in the |pspace_class|**, i.e. for the previous example, the |pspace_class| will be expecting a variable named ``gsmf_phi0``. All |pdist_class| subclasses optionally accept a ``default=`` keyword-argument, for example, ``PD_Uniform("hard_time", 0.1, 11.0, default=3.0)``. The 'default' values are provided so that |pspace_class|'s can construct a model using default parameters (see: :py:meth:`holodeck.librarian.lib_tools._Param_Space.default_params`), typically as a fiducial model. In the preceding example, the default 'hard_time' parameter would be 3.0. If a ``default`` is not specified in the instance constructor, then the value produced by an input of ``0.5`` is used. In the preceding example, if no ``default`` was specified, then the middle value of $(11.0 + 0.1) / 2 = 5.55$ would be used. Parameter Spaces ---------------- @@ -41,9 +41,9 @@ Parameter spaces must subclass |pspace_class|, and provide 4 elements: * *|ps_test_class| Example:* while this example construct a 3-dimensional parameter space (over "hard_time", "hard_gamma_inner", "mmb_mamp"), there are ``DEFAULTS`` specified for all of the parameters used to construct the GSMF, GMR, M-MBulge, and hardening models. -(1) An ``__init__()`` method that passes all required parameter distributions (:py:class:`~holodeck.librarian.libraries._Param_Dist` subclasses) to the super-class ``__init__()`` method. The list of |pdist_class| instances is where the actual parameter-space being explored is defined. Adding or removing a new element to this list of instances is all that it takes to increase or decrease the parameter space. +(1) An ``__init__()`` method that passes all required parameter distributions (:py:class:`~holodeck.librarian.lib_tools._Param_Dist` subclasses) to the super-class ``__init__()`` method. The list of |pdist_class| instances is where the actual parameter-space being explored is defined. Adding or removing a new element to this list of instances is all that it takes to increase or decrease the parameter space. - * *|ps_test_class| Example:* in this case, a 3-dimensional parameter space is constructed, using uniform distributions (:py:class:`~holodeck.librarian.libraries.PD_Uniform`) for "hard_time" and "hard_gamma_inner", and a normal (i.e. Gaussian, :py:class:`~holodeck.librarian.libraries.PD_Normal`) distribution for "hard_time". + * *|ps_test_class| Example:* in this case, a 3-dimensional parameter space is constructed, using uniform distributions (:py:class:`~holodeck.librarian.lib_tools.PD_Uniform`) for "hard_time" and "hard_gamma_inner", and a normal (i.e. Gaussian, :py:class:`~holodeck.librarian.lib_tools.PD_Normal`) distribution for "hard_time". (2) An ``_init_sam()`` function that takes the input parameters, and then constructs and returns a |sam_class| instance. @@ -92,9 +92,9 @@ Using holodeck libraries Loading a saved parameter-space instance ---------------------------------------- -TLDR: Use the :py:func:`~holodeck.librarian.libraries.load_pspace_from_path` function, passing in the path to the directory containing the save file (a ``.pspace.npz`` file). +TLDR: Use the :py:func:`~holodeck.librarian.lib_tools.load_pspace_from_path` function, passing in the path to the directory containing the save file (a ``.pspace.npz`` file). -Typically all that is needed for using/analyzing a holodeck library is the combined library output file ``sam_lib.hdf5``. A saved instance of the parameter-space class which generated the library is also saved to the output directory (as a ``.pspace.npz`` file), and can be useful for some use cases, for example if new simulations/realizations are desired from the same parameter space. The |pspace_class| provides a method to load saved instances, see the :py:meth:`~holodeck.librarian.libraries._Param_Space.from_save` method. Typically, the best way to load a saved parameter-space instance is to use the :py:func:`~holodeck.librarian.libraries.load_pspace_from_path` function. +Typically all that is needed for using/analyzing a holodeck library is the combined library output file ``sam_lib.hdf5``. A saved instance of the parameter-space class which generated the library is also saved to the output directory (as a ``.pspace.npz`` file), and can be useful for some use cases, for example if new simulations/realizations are desired from the same parameter space. The |pspace_class| provides a method to load saved instances, see the :py:meth:`~holodeck.librarian.lib_tools._Param_Space.from_save` method. Typically, the best way to load a saved parameter-space instance is to use the :py:func:`~holodeck.librarian.lib_tools.load_pspace_from_path` function. The combined holodeck library file ``sam_lib.hdf5`` --------------------------------------------------- diff --git a/docs/source/header.rst b/docs/source/header.rst index cc4f53da..09d800fb 100644 --- a/docs/source/header.rst +++ b/docs/source/header.rst @@ -1,6 +1,6 @@ .. |holodeck| replace:: :py:mod:`holodeck` .. |sam_class| replace:: :py:class:`~holodeck.sams.sam.Semi_Analytic_Model` .. |hard_class| replace:: :py:class:`~holodeck.hardening._Hardening` -.. |pspace_class| replace:: :py:class:`~holodeck.librarian.libraries._Param_Space` -.. |pdist_class| replace:: :py:class:`~holodeck.librarian.libraries._Param_Dist` +.. |pspace_class| replace:: :py:class:`~holodeck.librarian.lib_tools._Param_Space` +.. |pdist_class| replace:: :py:class:`~holodeck.librarian.lib_tools._Param_Dist` .. |ps_test_class| replace:: :py:class:`~holodeck.librarian.param_spaces.PS_Test` diff --git a/holodeck/librarian/__init__.py b/holodeck/librarian/__init__.py index 427bfdef..a98ce30e 100644 --- a/holodeck/librarian/__init__.py +++ b/holodeck/librarian/__init__.py @@ -4,10 +4,10 @@ producing different populations and/or GW signatures at each sampled parameter value. Libraries are run from the same parameter-space and using the same hyper parameters. Libraries are constructed using a 'parameter space' class that organizes the different simulations. The base-class is -:class:`~holodeck.librarian.libraries._Param_Space` (defined in the :mod:`holodeck.librarian.libraries` +:class:`~holodeck.librarian.lib_tools._Param_Space` (defined in the :mod:`holodeck.librarian.lib_tools` file). The parameter-space subclasses are given a number of different parameters to be varied. -Each parameter is implemented as a subclass of :py:class:`~holodeck.librarian.libraries._Param_Dist`, -for example the :py:class:`~holodeck.librarian.libraries.PD_Uniform` class that implements a uniform +Each parameter is implemented as a subclass of :py:class:`~holodeck.librarian.lib_tools._Param_Dist`, +for example the :py:class:`~holodeck.librarian.lib_tools.PD_Uniform` class that implements a uniform distribution. For more information, see the :doc:`'libraries' page in the getting-started guide @@ -18,9 +18,9 @@ The ``librarian`` module is composed of the following elements: * The core components of the holodeck libraries are defined in - :py:mod:`~holodeck.librarian.libraries`. Constructing simulations from parameter spaces can be + :py:mod:`~holodeck.librarian.lib_tools`. Constructing simulations from parameter spaces can be performed using the relevant parameter spaces themselves (subclasses of - :py:class:`~holodeck.librarian.libraries._Param_Space`). + :py:class:`~holodeck.librarian.lib_tools._Param_Space`). * Parameter spaces are defined in the 'param_spaces' files, particularly: @@ -39,6 +39,10 @@ """ +# ============================================================================== +# ==== Submodule Definitions ==== +# ============================================================================== + __version__ = "1.2" DEF_NUM_REALS = 100 #: Default number of realizations to construct in libraries. @@ -67,11 +71,18 @@ def __init__(self, message="This looks like a 'domain' not a 'library'!"): # Call the base class constructor with the parameters it needs super(DomainNotLibraryError, self).__init__(message) -from . libraries import ( # noqa - _Param_Space, _Param_Dist, - PD_Uniform, PD_Normal, - run_model, load_pspace_from_path, -) +# ============================================================================== +# ==== Import Submodule Components ==== +# ============================================================================== + +from . lib_tools import * # noqa +from . import gen_lib # noqa + +# from . lib_tools import ( # noqa +# _Param_Space, _Param_Dist, +# PD_Uniform, PD_Normal, +# run_model, load_pspace_from_path, +# ) param_spaces_dict = {} #: Registry of standard parameter-spaces from . import param_spaces_classic as psc_temp # noqa diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index ede1a380..9d3c69b3 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -21,7 +21,7 @@ import holodeck.librarian import holodeck.librarian.gen_lib from holodeck.librarian import ( - libraries, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS, DomainNotLibraryError + lib_tools, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS, DomainNotLibraryError ) @@ -122,7 +122,7 @@ def sam_lib_combine( # ---- see if a combined library already exists - lib_path = libraries.get_sam_lib_fname(path_output, gwb_only) + lib_path = lib_tools.get_sam_lib_fname(path_output, gwb_only) if lib_path.exists(): lvl = log.INFO if recreate else log.WARNING log.log(lvl, f"combined library already exists: {lib_path}, run with `-r` to recreate.") @@ -135,7 +135,7 @@ def sam_lib_combine( if path_pspace is None: path_pspace = path_output - pspace, pspace_fname = libraries.load_pspace_from_path(path_pspace, log=log) + pspace, pspace_fname = lib_tools.load_pspace_from_path(path_pspace, log=log) args, args_fname = holo.librarian.gen_lib.load_config_from_path(path_pspace, log=log) log.info(f"loaded param space: {pspace} from '{pspace_fname}'") @@ -300,7 +300,7 @@ def _check_files_and_load_shapes(log, path_sims, nsamp, library): log.info(f"Checking {nsamp} files in {path_sims}") for ii in tqdm.trange(nsamp): - temp_fname = libraries._get_sim_fname(path_sims, ii, library=library) + temp_fname = lib_tools._get_sim_fname(path_sims, ii, library=library) if not temp_fname.exists(): err = f"Missing at least file number {ii} out of {nsamp} files! {temp_fname}" log.exception(err) @@ -396,7 +396,7 @@ def _load_library_from_all_files( bad_files = np.zeros(nsamp_all, dtype=bool) #: track which files contain UN-useable data msg = None for pnum in tqdm.trange(nsamp_all): - fname = libraries._get_sim_fname(path_sims, pnum, library=library) + fname = lib_tools._get_sim_fname(path_sims, pnum, library=library) temp = np.load(fname, allow_pickle=True) # When a processor fails for a given parameter, the output file is still created with the 'fail' key added if ('fail' in temp): diff --git a/holodeck/librarian/fit_spectra.py b/holodeck/librarian/fit_spectra.py index 4da08e60..e57f0432 100644 --- a/holodeck/librarian/fit_spectra.py +++ b/holodeck/librarian/fit_spectra.py @@ -43,7 +43,7 @@ import holodeck as holo import holodeck.librarian -from holodeck.librarian import libraries, FITS_NBINS_PLAW, FITS_NBINS_TURN +from holodeck.librarian import lib_tools, FITS_NBINS_PLAW, FITS_NBINS_TURN from holodeck.constants import YR @@ -108,7 +108,7 @@ def fit_library_spectra(library_path, log, recreate=False): library_path = Path(library_path) if library_path.is_dir(): - library_path = libraries.get_sam_lib_fname(library_path, gwb_only=False) + library_path = lib_tools.get_sam_lib_fname(library_path, gwb_only=False) if not library_path.exists() or not library_path.is_file(): err = f"{library_path=} must point to an existing library file!" log.exception(err) @@ -118,7 +118,7 @@ def fit_library_spectra(library_path, log, recreate=False): # ---- check for existing fits file - fits_path = libraries.get_fits_path(library_path) + fits_path = lib_tools.get_fits_path(library_path) return_flag = False if fits_path.exists(): lvl = log.INFO if recreate else log.WARNING diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index 48e8814d..57dcf84f 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -3,7 +3,7 @@ This file can be run from the command-line to generate holodeck libraries, and also provides some API methods for quick/easy generation of simulations. In general, these methods are designed to run simulations for populations constructed from parameter-spaces (i.e. -:class:`~holodeck.librarian.libraries._Param_Space` subclasses). This script is parallelized using +:class:`~holodeck.librarian.lib_tools._Param_Space` subclasses). This script is parallelized using ``mpi4py``, but can also be run in serial. This script can be run by executing:: @@ -35,7 +35,7 @@ import holodeck.librarian import holodeck.librarian.combine from holodeck.librarian import ( - libraries, ARGS_CONFIG_FNAME, PSPACE_DOMAIN_EXTREMA, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS + lib_tools, ARGS_CONFIG_FNAME, PSPACE_DOMAIN_EXTREMA, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS ) #: maximum number of failed simulations before task terminates with error (`None`: no limit) @@ -243,16 +243,16 @@ def run_sam_at_pspace_params(args, space, pnum, params): ``True``, otherwise the function runs this simulation. (2) Calls ``space.model_for_params`` to generate the semi-analytic model and hardening instances; see the function - :func:`holodeck.librarian.libraries._Param_Space.model_for_params()`. + :func:`holodeck.librarian.lib_tools._Param_Space.model_for_params()`. (3) Calculates populations and GW signatures from the SAM and hardening model using - :func:`holodeck.librarian.libraries.run_model()`, and saves the results to an output file. + :func:`holodeck.librarian.lib_tools.run_model()`, and saves the results to an output file. (4) Optionally: some diagnostic plots are created in the :func:`make_plots()` function. Arguments --------- args : ``argparse.ArgumentParser`` instance Arguments from the ``gen_lib_sams.py`` script. - space : :class:`holodeck.librarian.libraries._Param_space` instance + space : :class:`holodeck.librarian.lib_tools._Param_space` instance Parameter space from which to construct populations. pnum : int Which parameter-sample from ``space`` should be run. @@ -279,7 +279,7 @@ def run_sam_at_pspace_params(args, space, pnum, params): # ---- get output filename for this simulation, check if already exists library_flag = not args.domain - sim_fname = libraries._get_sim_fname(args.output_sims, pnum, library=library_flag) + sim_fname = lib_tools._get_sim_fname(args.output_sims, pnum, library=library_flag) beg = datetime.now() log.info(f"{pnum=} :: {sim_fname=} beginning at {beg}") @@ -301,7 +301,7 @@ def run_sam_at_pspace_params(args, space, pnum, params): log.debug("Selecting `sam` and `hard` instances") sam, hard = space.model_for_params(params) - data = libraries.run_model( + data = lib_tools.run_model( sam, hard, pta_dur=args.pta_dur, nfreqs=args.nfreqs, nreals=args.nreals, nloudest=args.nloudest, gwb_flag=args.gwb_flag, singles_flag=args.ss_flag, details_flag=False, params_flag=args.params_flag, diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/lib_tools.py similarity index 99% rename from holodeck/librarian/libraries.py rename to holodeck/librarian/lib_tools.py index b9261495..0042a6a1 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/lib_tools.py @@ -26,6 +26,10 @@ "gsmf_phi0": ["gsmf_phi0_log10", None], } +# __all__ = [ +# "_Param_space" +# ] + class _Param_Space(abc.ABC): """Base class for generating holodeck libraries. Defines the parameter space and settings. @@ -60,7 +64,7 @@ def __init__(self, parameters, log=None, nsamples=None, sam_shape=None, seed=Non param_kwargs : dict Key-value pairs specifying the parameters for this model. Each key must be the name of the parameter, and each value must be a - :class:`~holodeck.librarian.libraries._Param_Dist` + :class:`~holodeck.librarian.lib_tools._Param_Dist` subclass instance with the desired distribution. Returns diff --git a/holodeck/librarian/param_spaces.py b/holodeck/librarian/param_spaces.py index cff7f5d3..74da3bd5 100644 --- a/holodeck/librarian/param_spaces.py +++ b/holodeck/librarian/param_spaces.py @@ -3,7 +3,7 @@ import holodeck as holo from holodeck.constants import GYR, PC, MSOL -from holodeck.librarian.libraries import _Param_Space, PD_Uniform, PD_Normal, PD_Uniform_Log +from holodeck.librarian.lib_tools import _Param_Space, PD_Uniform, PD_Normal, PD_Uniform_Log # Define a new Parameter-Space class by subclassing the base class: diff --git a/holodeck/librarian/param_spaces_classic.py b/holodeck/librarian/param_spaces_classic.py index 8e90e705..6f49331a 100644 --- a/holodeck/librarian/param_spaces_classic.py +++ b/holodeck/librarian/param_spaces_classic.py @@ -2,7 +2,7 @@ """ from holodeck.constants import PC, GYR -from holodeck.librarian.libraries import _Param_Space, PD_Uniform, PD_Normal +from holodeck.librarian.lib_tools import _Param_Space, PD_Uniform, PD_Normal from holodeck import sams, hardening, host_relations diff --git a/holodeck/librarian/tests/test_libraries__param_space.py b/holodeck/librarian/tests/test_lib_tools__param_space.py similarity index 97% rename from holodeck/librarian/tests/test_libraries__param_space.py rename to holodeck/librarian/tests/test_lib_tools__param_space.py index 57ef580f..7d037ccc 100644 --- a/holodeck/librarian/tests/test_libraries__param_space.py +++ b/holodeck/librarian/tests/test_lib_tools__param_space.py @@ -1,4 +1,4 @@ -"""Test holodeck/librarian/libraries.py: the parameter space base class, and parameter distributions. +"""Test holodeck/librarian/lib_tools.py: the parameter space base class, and parameter distributions. """ import numpy as np @@ -12,7 +12,7 @@ from holodeck.librarian import ( param_spaces_dict ) -from holodeck.librarian.libraries import ( +from holodeck.librarian.lib_tools import ( _Param_Space, PD_Uniform, ) @@ -205,8 +205,8 @@ def _check_sam_hard(sam, hard, sam_shape): assert isinstance(hard, hardening.Hard_GW) # Make sure model runs - import holodeck.librarian.libraries # noqa - data = librarian.libraries.run_model(sam, hard) + import holodeck.librarian.lib_tools # noqa + data = librarian.lib_tools.run_model(sam, hard) assert data is not None return diff --git a/holodeck/librarian/tests/test_param_spaces.py b/holodeck/librarian/tests/test_param_spaces.py index 9b6a15af..7c24f71c 100644 --- a/holodeck/librarian/tests/test_param_spaces.py +++ b/holodeck/librarian/tests/test_param_spaces.py @@ -14,8 +14,8 @@ def _run_param_space(param_space_class): # sam, hard = pspace.model_for_sample_number(0) sam, hard = pspace.model_for_params(pspace.default_params()) # Make sure model runs - import holodeck.librarian.libraries # noqa - data = librarian.libraries.run_model(sam, hard, singles_flag=True, details_flag=True) + import holodeck.librarian.lib_tools # noqa + data = librarian.lib_tools.run_model(sam, hard, singles_flag=True, details_flag=True) assert data is not None, "After `run_model` returned data is None!" check_keys = ['fobs_cents', 'fobs_edges', 'hc_ss', 'hc_bg', 'gwb', 'static_binary_density', 'number'] for key in check_keys: diff --git a/notebooks/devs/libraries/gwb_anatomy.ipynb b/notebooks/devs/libraries/gwb_anatomy.ipynb index 733c7334..a26ff37d 100644 --- a/notebooks/devs/libraries/gwb_anatomy.ipynb +++ b/notebooks/devs/libraries/gwb_anatomy.ipynb @@ -132,7 +132,7 @@ "outputs": [], "source": [ "import holodeck.librarian.old_param_spaces\n", - "import holodeck.librarian.libraries\n", + "import holodeck.librarian.lib_tools\n", "\n", "# construct a param_space instance, note that `nsamples` and `seed` don't matter here for how we'll use this\n", "pspace = holo.librarian.old_param_spaces.PS_Uniform_07B_Rot(holo.log, nsamples=1, sam_shape=SHAPE, seed=None)\n", @@ -148,7 +148,7 @@ "fiducial_parameters = pspace._normalized_params(fiducial_pars)\n", "fiducial_sam, fiducial_hard = pspace.model_for_params(fiducial_parameters, sam_shape=SHAPE) # <-- slow for some reason?\n", "#fiducial_sam, fiducial_hard = pspace.model_for_normalized_params(fiducial_pars)\n", - "fiducial_model_data = holo.librarian.libraries.run_model(fiducial_sam, fiducial_hard, nreals=NREALS, gwb_flag=True, details_flag=True) # Only run fiducial model once to save time!\n", + "fiducial_model_data = holo.librarian.lib_tools.run_model(fiducial_sam, fiducial_hard, nreals=NREALS, gwb_flag=True, details_flag=True) # Only run fiducial model once to save time!\n", "\n", "alldata = []\n", "allparams_list = {}\n", @@ -156,7 +156,7 @@ "fiducial_gwonly_hard = holo.hardening.Hard_GW()\n", "fiducial_sam.ZERO_DYNAMIC_STALLED_SYSTEMS = False\n", "fiducial_sam.ZERO_GMT_STALLED_SYSTEMS = True\n", - "fiducial_gwonly_model_data = holo.librarian.libraries.run_model(\n", + "fiducial_gwonly_model_data = holo.librarian.lib_tools.run_model(\n", " fiducial_sam, fiducial_gwonly_hard, nreals=NREALS, gwb_flag=True, details_flag=True\n", ")\n", "print(fiducial_gwonly_model_data.keys())\n", @@ -199,7 +199,7 @@ " # sam, hard = pspace.model_for_normalized_params(pars)\n", " sam, hard = pspace.model_for_params(parameters, sam_shape=SHAPE) # <-- slow for some reason?\n", " # run this model, retrieving binary parameters and the GWB\n", - " _data = holo.librarian.libraries.run_model(sam, hard, nreals=NREALS, gwb_flag=True, details_flag=True)\n", + " _data = holo.librarian.lib_tools.run_model(sam, hard, nreals=NREALS, gwb_flag=True, details_flag=True)\n", " data.append(_data)\n", " mingwb = np.min([mingwb, np.min(_data['gwb'])])\n", " maxgwb = np.max([mingwb, np.max(_data['gwb'])])\n", diff --git a/notebooks/devs/libraries/library-explorer.ipynb b/notebooks/devs/libraries/library-explorer.ipynb index fd377fa6..12444622 100644 --- a/notebooks/devs/libraries/library-explorer.ipynb +++ b/notebooks/devs/libraries/library-explorer.ipynb @@ -356,7 +356,7 @@ "CREATE = False\n", "RECREATE = False\n", "import holodeck.librarian.fit_spectra\n", - "fits_path = holo.librarian.libraries.get_fits_path(lib_path)\n", + "fits_path = holo.librarian.lib_tools.get_fits_path(lib_path)\n", "print(f\"{fits_path.name=} {fits_path.exists()=}\")\n", "if CREATE or RECREATE:\n", " if (CREATE and (not fits_path.exists())) or RECREATE:\n", @@ -1349,7 +1349,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.11.7" }, "vscode": { "interpreter": { diff --git a/notebooks/devs/libraries/param_spaces.ipynb b/notebooks/devs/libraries/param_spaces.ipynb index a924be2e..fe108d2d 100644 --- a/notebooks/devs/libraries/param_spaces.ipynb +++ b/notebooks/devs/libraries/param_spaces.ipynb @@ -11,7 +11,6 @@ "\n", "import holodeck as holo\n", "from holodeck import librarian\n", - "import holodeck.librarian.libraries\n", "from holodeck.constants import MSOL, GYR, YR" ] }, @@ -35,7 +34,7 @@ "metadata": {}, "outputs": [], "source": [ - "data = librarian.libraries.run_model(sam, hard)" + "data = librarian.lib_tools.run_model(sam, hard)" ] }, { diff --git a/notebooks/librarian.ipynb b/notebooks/librarian.ipynb index c181a35e..bad0b716 100644 --- a/notebooks/librarian.ipynb +++ b/notebooks/librarian.ipynb @@ -48,7 +48,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Holodeck libraries are build around their parameter-spaces, which are implemented as subclasses of the [`holodeck.librarian.libraries._Param_Space`](https://holodeck-gw.readthedocs.io/en/main/api_ref/holodeck.librarian.libraries.html#holodeck.librarian.libraries._Param_Space) base class. These subclasses should be named with `PS_` prefixes to denote that they are parameter spaces." + "Holodeck libraries are build around their parameter-spaces, which are implemented as subclasses of the [`holodeck.librarian.lib_tools._Param_Space`](https://holodeck-gw.readthedocs.io/en/main/api_ref/holodeck.librarian.lib_tools.html#holodeck.librarian.lib_tools._Param_Space) base class. These subclasses should be named with `PS_` prefixes to denote that they are parameter spaces." ] }, { @@ -142,7 +142,7 @@ "metadata": {}, "outputs": [], "source": [ - "data = holo.librarian.libraries.run_model(sam, hard, details_flag=True)\n", + "data = holo.librarian.lib_tools.run_model(sam, hard, details_flag=True)\n", "print(data.keys())\n", "data['gwb']" ] @@ -161,13 +161,13 @@ "metadata": {}, "outputs": [], "source": [ - "class PS_Lib_Test(holo.librarian.libraries._Param_Space):\n", + "class PS_Lib_Test(holo.librarian.lib_tools._Param_Space):\n", "\n", " DEFAULTS = {'gsmf_phi0_log10': -2.77}\n", "\n", " def __init__(self, log, nsamples, sam_shape):\n", " parameters = [\n", - " holo.librarian.libraries.PD_Normal(\"gsmf_phi0_log10\", -2.77, 0.3),\n", + " holo.librarian.lib_tools.PD_Normal(\"gsmf_phi0_log10\", -2.77, 0.3),\n", " ]\n", " super().__init__(\n", " parameters,\n", @@ -215,7 +215,7 @@ "metadata": {}, "outputs": [], "source": [ - "test = holo.librarian.libraries.PD_Uniform(\"test\", -10.0, 10.0)\n", + "test = holo.librarian.lib_tools.PD_Uniform(\"test\", -10.0, 10.0)\n", "assert test(0.5) == 0.0\n", "assert test(0.0) == -10.0\n", "assert test(1.0) == +10.0\n", @@ -240,7 +240,7 @@ "metadata": {}, "outputs": [], "source": [ - "test = holo.librarian.libraries.PD_Normal(\"test\", 0.0, 1.0)\n", + "test = holo.librarian.lib_tools.PD_Normal(\"test\", 0.0, 1.0)\n", "val = test(0.5)\n", "print(val)\n", "assert test(0.5) == 0.0\n", @@ -265,7 +265,7 @@ "metadata": {}, "outputs": [], "source": [ - "test = holo.librarian.libraries.PD_Lin_Log(\"test\", 0.01, 100.0, 0.1, 0.5)\n", + "test = holo.librarian.lib_tools.PD_Lin_Log(\"test\", 0.01, 100.0, 0.1, 0.5)\n", "xx = np.linspace(0.0, 1.0, 10000)\n", "yy = test(xx)\n", "print(utils.minmax(yy))\n", @@ -300,7 +300,7 @@ "\n", "fig, ax = plot.figax(scale='log')\n", "for frac in [0.2, 0.5, 0.8]:\n", - " test = holo.librarian.libraries.PD_Lin_Log(\"test\", 0.01, 100.0, crit, frac)\n", + " test = holo.librarian.lib_tools.PD_Lin_Log(\"test\", 0.01, 100.0, crit, frac)\n", " xx = test(np.random.uniform(0.0, 1.0, size=NUM))\n", " kale.dist1d(xx, ax=ax, edges=edges, density=True, probability=False)\n", " obs_frac = np.count_nonzero(xx < crit) / xx.size\n", @@ -332,7 +332,7 @@ "\n", "fig, ax = plot.figax(scale='log')\n", "for crit in [0.1, 1.0, 10.0]:\n", - " test = holo.librarian.libraries.PD_Lin_Log(\"test\", 0.01, 100.0, crit, frac)\n", + " test = holo.librarian.lib_tools.PD_Lin_Log(\"test\", 0.01, 100.0, crit, frac)\n", " xx = test(np.random.uniform(0.0, 1.0, size=NUM))\n", " kale.dist1d(xx, ax=ax, edges=edges, density=True, probability=False)\n", " obs_frac = np.count_nonzero(xx < crit) / xx.size\n", @@ -356,7 +356,7 @@ "metadata": {}, "outputs": [], "source": [ - "test = holo.librarian.libraries.PD_Log_Lin(\"test\", 0.01, 100.0, 0.1, 0.5)\n", + "test = holo.librarian.lib_tools.PD_Log_Lin(\"test\", 0.01, 100.0, 0.1, 0.5)\n", "xx = np.linspace(0.0, 1.0, 10000)\n", "yy = test(xx)\n", "print(utils.minmax(yy))\n", @@ -389,7 +389,7 @@ "\n", "fig, ax = plot.figax(scale='log')\n", "for frac in [0.2, 0.5, 0.8]:\n", - " test = holo.librarian.libraries.PD_Log_Lin(\"test\", 0.01, 100.0, crit, frac)\n", + " test = holo.librarian.lib_tools.PD_Log_Lin(\"test\", 0.01, 100.0, crit, frac)\n", " xx = test(np.random.uniform(0.0, 1.0, size=NUM))\n", " kale.dist1d(xx, ax=ax, edges=edges, density=True, probability=False)\n", " obs_frac = np.count_nonzero(xx < crit) / xx.size\n", @@ -422,7 +422,7 @@ "\n", "fig, ax = plot.figax(scale='log')\n", "for crit in [0.1, 1.0, 10.0]:\n", - " test = holo.librarian.libraries.PD_Log_Lin(\"test\", 0.01, 100.0, crit, frac)\n", + " test = holo.librarian.lib_tools.PD_Log_Lin(\"test\", 0.01, 100.0, crit, frac)\n", " xx = test(np.random.uniform(0.0, 1.0, size=NUM))\n", " kale.dist1d(xx, ax=ax, edges=edges, density=True, probability=False)\n", " obs_frac = np.count_nonzero(xx < crit) / xx.size\n", @@ -448,7 +448,7 @@ "source": [ "edges = [-1.0, 5.0, 6.0, 7.0]\n", "amps = [1.0, 2.0, 1.0]\n", - "test = holodeck.librarian.libraries.PD_Piecewise_Uniform_Mass(\"test\", edges, amps)\n", + "test = holodeck.librarian.lib_tools.PD_Piecewise_Uniform_Mass(\"test\", edges, amps)\n", "\n", "xx = np.random.uniform(size=1000)\n", "xx = np.sort(xx)\n", @@ -466,7 +466,7 @@ "outputs": [], "source": [ "edges = [-1.0, 4.0, 6.0, 7.5]\n", - "test = holodeck.librarian.libraries.PD_Piecewise_Uniform_Density(\"test\", edges, [1.0, 2.0, 1.0])\n", + "test = holodeck.librarian.lib_tools.PD_Piecewise_Uniform_Density(\"test\", edges, [1.0, 2.0, 1.0])\n", "\n", "xx = np.random.uniform(size=1000)\n", "xx = np.sort(xx)\n", @@ -491,7 +491,7 @@ "outputs": [], "source": [ "edges = [0.1, 1.0, 9.0, 11.0]\n", - "test = holodeck.librarian.libraries.PD_Piecewise_Uniform_Density(\"test\", edges, [2.5, 0.5, 1.5])\n", + "test = holodeck.librarian.lib_tools.PD_Piecewise_Uniform_Density(\"test\", edges, [2.5, 0.5, 1.5])\n", "\n", "xx = np.random.uniform(size=2000)\n", "xx = np.sort(xx)\n", @@ -514,7 +514,7 @@ "metadata": {}, "outputs": [], "source": [ - "test = holodeck.librarian.libraries.PD_Piecewise_Uniform_Density(\n", + "test = holodeck.librarian.lib_tools.PD_Piecewise_Uniform_Density(\n", " \"test\", [7.5, 8.0, 9.0, 9.5], [1.5, 1.0, 2.0]\n", ")\n", "\n", From a19531f4678bfabd7fdaf4e16559c9a64498a792 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Tue, 30 Apr 2024 11:17:51 -0700 Subject: [PATCH 43/48] BUG: minor bug fixes --- holodeck/librarian/combine.py | 14 ++++++++++++-- holodeck/librarian/gen_lib.py | 2 +- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index ede1a380..71a3039e 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -122,7 +122,7 @@ def sam_lib_combine( # ---- see if a combined library already exists - lib_path = libraries.get_sam_lib_fname(path_output, gwb_only) + lib_path = libraries.get_sam_lib_fname(path_output, gwb_only, library=library) if lib_path.exists(): lvl = log.INFO if recreate else log.WARNING log.log(lvl, f"combined library already exists: {lib_path}, run with `-r` to recreate.") @@ -217,7 +217,12 @@ def sam_lib_combine( with h5py.File(lib_path, 'w') as h5: h5.create_dataset('fobs_cents', data=fobs_cents) h5.create_dataset('fobs_edges', data=fobs_edges) + + if not library: + new_shape = (ndim, nsamp_dim) + param_samples.shape[1:] + param_samples = param_samples.reshape(new_shape) h5.create_dataset('sample_params', data=param_samples) + if gwb is not None: # if 'domain', reshape to include each dimension if not library: @@ -416,7 +421,12 @@ def _load_library_from_all_files( # for 'domain' simulations, we need to load the parameters. For 'library' runs, we # already have them. if not library: - param_samples[pnum, :] = temp['params'][:] + #! NOTE: this is just temporary until bug is fixed + try: + param_samples[pnum, :] = temp['params'][:] + except: + pvals = temp['params'][()] + param_samples[pnum, :] = np.array([pvals[str(pn)] for pn in temp['param_names']]) # store the GWB from this file if gwb is not None: diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index 48e8814d..51b475e1 100644 --- a/holodeck/librarian/gen_lib.py +++ b/holodeck/librarian/gen_lib.py @@ -307,7 +307,7 @@ def run_sam_at_pspace_params(args, space, pnum, params): gwb_flag=args.gwb_flag, singles_flag=args.ss_flag, details_flag=False, params_flag=args.params_flag, log=log, ) - data['params'] = params + data['params'] = np.array([params[pn] for pn in space.param_names]) data['param_names'] = space.param_names rv = True From b1d7405b1e502f5312f8d120f368c7aa2542c57b Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 1 May 2024 10:21:56 -0700 Subject: [PATCH 44/48] DOCS: update CHANGELOG and librarian version --- CHANGELOG.md | 3 +++ holodeck/librarian/__init__.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d6b66d66..5dee824a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,9 @@ * DEPRECATIONS: * `holodeck.librarian.py` ==> `holodeck.librarian.lib_tools` * Rename submodule. All components remain the same. All `lib_tools` elements are now also imported into the `librarian` namespace. i.e. elements like `holodeck.librarian.lib_tools._Param_Space` will now be accessible via `holodeck.librarian._Param_Space`. + * Library filenames: + * Standard library simulation files will now be saved to the 'library_sims' subdirectory, and filenames 'library__p######.npz'. Combined library files will now be 'sam-library.hdf5'. + * 'Domain' simulation files will be saved to the 'domain_sims' subdirectory, and filenames 'domain__p######.npz'. Combined domain files will now be 'sam-domain.hdf5'. * BUGS: * Semi-Analytic Models diff --git a/holodeck/librarian/__init__.py b/holodeck/librarian/__init__.py index a98ce30e..e0380401 100644 --- a/holodeck/librarian/__init__.py +++ b/holodeck/librarian/__init__.py @@ -43,7 +43,7 @@ # ==== Submodule Definitions ==== # ============================================================================== -__version__ = "1.2" +__version__ = "1.3" DEF_NUM_REALS = 100 #: Default number of realizations to construct in libraries. DEF_NUM_FBINS = 40 #: Default number of frequency bins at which to calculate GW signals. From 646985f85fd955ab8f5d3542c82b5fc216614442 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 1 May 2024 10:38:05 -0700 Subject: [PATCH 45/48] BUG: temporarily remove some rarely-used packages from the global requirements file, as they are causing problems --- holodeck/anisotropy.py | 144 +++++++++++++++++++++------------------ holodeck/detstats.py | 24 +++++-- holodeck/gps/__init__.py | 15 ++++ requirements.txt | 12 ++-- 4 files changed, 116 insertions(+), 79 deletions(-) diff --git a/holodeck/anisotropy.py b/holodeck/anisotropy.py index b6725f99..c20ea1cd 100644 --- a/holodeck/anisotropy.py +++ b/holodeck/anisotropy.py @@ -3,26 +3,38 @@ """ import numpy as np -import matplotlib as plt +# import matplotlib as plt import matplotlib.cm as cm import kalepy as kale -import healpy as hp import h5py import holodeck as holo -from holodeck import utils, cosmo, log, detstats, plot +from holodeck import utils, log, detstats, plot from holodeck.constants import YR +try: + import healpy as hp +except ImportError as err: + SUBMOD = "anisotropy" + log.error(f"Failed to import `healpy` packaged used in `{SUBMOD}` submodule!") + log.exception(err) + log.error( + f"Some required packages for `{SUBMOD}` have been temporarily disabled in the " + "global 'requirements.txt' file, so they are not installed by default! Please install " + "the required packages manually for now, and feel free to raise a github issue." + ) + raise + NSIDE = 32 NPIX = hp.nside2npix(NSIDE) LMAX = 8 -HC_REF15_10YR = 11.2*10**-15 +HC_REF15_10YR = 11.2*10**-15 def healpix_map(hc_ss, hc_bg, nside=NSIDE, seed=None, ret_seed=False): """ Build mollview array of hc^2/dOmega for a healpix map - + Parameters ---------- hc_ss : (F,R,L) NDarray @@ -36,7 +48,7 @@ def healpix_map(hc_ss, hc_bg, nside=NSIDE, seed=None, ret_seed=False): ------- moll_hc : (NPIX,) 1Darray Array of h_c^2/dOmega at every pixel for a mollview healpix map. - + NOTE: Could speed up the for-loops, but it's ok for now. """ @@ -50,7 +62,7 @@ def healpix_map(hc_ss, hc_bg, nside=NSIDE, seed=None, ret_seed=False): if seed is None: seed = np.random.randint(99999) # get a random number print(f"random seed: {seed}") # print it out so we can reuse it if desired - np.random.seed(seed) + np.random.seed(seed) # spread background evenly across pixels in moll_hc moll_hc = np.ones((nfreqs,nreals,npix)) * hc_bg[:,:,np.newaxis]**2/(npix*area) # (frequency, realization, pixel) @@ -62,12 +74,12 @@ def healpix_map(hc_ss, hc_bg, nside=NSIDE, seed=None, ret_seed=False): for ll in range(nloudest): moll_hc[ff,rr,pix_ss[ff,rr,ll]] = (moll_hc[ff,rr,pix_ss[ff,rr,ll]] + hc_ss[ff,rr,ll]**2/area) if ret_seed: - return moll_hc, seed + return moll_hc, seed return moll_hc def healpix_map_oldhc2(hc_ss, hc_bg, nside=NSIDE): """ Build mollview array of hc^2/dOmega for a healpix map - + Parameters ---------- hc_ss : (F,R,L) NDarray @@ -81,7 +93,7 @@ def healpix_map_oldhc2(hc_ss, hc_bg, nside=NSIDE): ------- moll_hc : (NPIX,) 1Darray Array of h_c^2 at every pixel for a mollview healpix map. - + NOTE: Could speed up the for-loops, but it's ok for now. """ @@ -100,12 +112,12 @@ def healpix_map_oldhc2(hc_ss, hc_bg, nside=NSIDE): for rr in range(nreals): for ll in range(nloudest): moll_hc[ff,rr,pix_ss[ff,rr,ll]] = (moll_hc[ff,rr,pix_ss[ff,rr,ll]] + hc_ss[ff,rr,ll]**2) - + return moll_hc def healpix_map_oldhc(hc_ss, hc_bg, nside=NSIDE): """ Build mollview array of strains for a healpix map - + Parameters ---------- hc_ss : (F,R,L) NDarray @@ -119,7 +131,7 @@ def healpix_map_oldhc(hc_ss, hc_bg, nside=NSIDE): ------- moll_hc : (NPIX,) 1Darray Array of strain at every pixel for a mollview healpix map. - + NOTE: Could speed up the for-loops, but it's ok for now. """ @@ -138,13 +150,13 @@ def healpix_map_oldhc(hc_ss, hc_bg, nside=NSIDE): for ll in range(nloudest): moll_hc[ff,rr,pix_ss[ff,rr,ll]] = np.sqrt(moll_hc[ff,rr,pix_ss[ff,rr,ll]]**2 + hc_ss[ff,rr,ll]**2) - + return moll_hc def healpix_hcsq_map(hc_ss, hc_bg, nside=NSIDE): """ Build mollview array of strains for a healpix map - + Parameters ---------- hc_ss : (F,R,L) NDarray @@ -158,7 +170,7 @@ def healpix_hcsq_map(hc_ss, hc_bg, nside=NSIDE): ------- moll_hc : (NPIX,) 1Darray Array of strain at every pixel for a mollview healpix map. - + NOTE: Could speed up the for-loops, but it's ok for now. """ @@ -176,13 +188,13 @@ def healpix_hcsq_map(hc_ss, hc_bg, nside=NSIDE): for rr in range(nreals): for ll in range(nloudest): moll_hc2[ff,rr,pix_ss[ff,rr,ll]] = moll_hc2[ff,rr,pix_ss[ff,rr,ll]] + hc_ss[ff,rr,ll]**2 - + return moll_hc2 def sph_harm_from_map(moll_hc, lmax=LMAX): - """ Calculate spherical harmonics from strains at every pixel of + """ Calculate spherical harmonics from strains at every pixel of a healpix mollview map. - + Parameters ---------- moll_hc : (F,R,NPIX,) 1Darray @@ -193,8 +205,8 @@ def sph_harm_from_map(moll_hc, lmax=LMAX): Returns ------- Cl : (F,R,lmax+1) NDarray - Spherical harmonic coefficients - + Spherical harmonic coefficients + """ nfreqs = len(moll_hc) nreals = len(moll_hc[0]) @@ -224,8 +236,8 @@ def sph_harm_from_hc(hc_ss, hc_bg, nside = NSIDE, lmax = LMAX): moll_hc : (F,R,NPIX,) 2Darray Array of hc^2/dOmega at every pixel for a mollview healpix map. Cl : (F,R,lmax+1) NDarray - Spherical harmonic coefficients - + Spherical harmonic coefficients + """ moll_hc = healpix_map(hc_ss, hc_bg, nside) Cl = sph_harm_from_map(moll_hc, lmax) @@ -251,13 +263,13 @@ def plot_ClC0_medians(fobs, Cl_best, lmax, nshow): percs = pp/2 percs = [50-percs, 50+percs] ax.fill_between(xx, *np.percentile(yy[:,:,ll], percs, axis=0), alpha=0.1, color=colors[ll]) - + for bb in range(0,nshow): ax.plot(xx, yy[bb,:,ll], color=colors[ll], linestyle=':', alpha=0.1, - linewidth=1) + linewidth=1) ax.legend(ncols=2) holo.plot._twin_hz(ax, nano=False) - + # ax.set_title('50%% and 98%% confidence intervals of the %d best samples \nusing realizations medians, lmax=%d' # % (nbest, lmax)) return fig @@ -288,7 +300,7 @@ def lib_anisotropy(lib_path, hc_ref_10yr=HC_REF15_10YR, nbest=100, nreals=50, lm # ---- rank samples nsort, fidx, hc_ref = detstats.rank_samples(hc_ss, hc_bg, fobs, fidx=1, hc_ref=hc_ref_10yr, ret_all=True) - + print('Ranked samples by hc_ref = %.2e at fobs = %.2f/yr' % (hc_ref, fobs[fidx]*YR)) @@ -301,7 +313,7 @@ def lib_anisotropy(lib_path, hc_ref_10yr=HC_REF15_10YR, nbest=100, nreals=50, lm print('on nn=%d out of nbest=%d' % (nn,nbest)) moll_hc_best[nn,...], Cl_best[nn,...] = sph_harm_from_hc( hc_ss[nsort[nn]], hc_bg[nsort[nn]], nside=nside, lmax=lmax, ) - + # ---- save to npz file @@ -319,10 +331,10 @@ def lib_anisotropy(lib_path, hc_ref_10yr=HC_REF15_10YR, nbest=100, nreals=50, lm np.savez(output_name, nsort=nsort, fidx=fidx, hc_ref=hc_ref, ss_shape=shape, moll_hc_best=moll_hc_best, Cl_best=Cl_best, nside=nside, lmax=lmax, fobs=fobs) - + # ---- plot median Cl/C0 - + print('Plotting Cl/C0 for median realizations') fig = plot_ClC0_medians(fobs, Cl_best, lmax, nshow=nbest) fig_name = output_dir+'/sph_harm_lmax%d_nside%d_nbest%d.png' % (lmax, nside, nbest) @@ -349,7 +361,7 @@ def lib_anisotropy_split(lib_path, hc_ref_10yr=HC_REF15_10YR, nbest=100, nreals= # ---- rank samples nsort, fidx, hc_ref = detstats.rank_samples(hc_ss, hc_bg, fobs, fidx=1, hc_ref=hc_ref_10yr, ret_all=True) - + print('Ranked samples by hc_ref = %.2e at fobs = %.2f/yr' % (hc_ref, fobs[fidx]*YR)) @@ -366,7 +378,7 @@ def lib_anisotropy_split(lib_path, hc_ref_10yr=HC_REF15_10YR, nbest=100, nreals= print('on nn=%d out of nbest=%d' % (nn,nbest)) moll_hc_best[ii,...], Cl_best[ii,...] = sph_harm_from_hc( hc_ss[nsort[nn]], hc_bg[nsort[nn]], nside=nside, lmax=lmax, ) - + # ---- save to npz file @@ -379,19 +391,19 @@ def lib_anisotropy_split(lib_path, hc_ref_10yr=HC_REF15_10YR, nbest=100, nreals= else: print('Writing to an existing directory.') - output_name =(output_dir+'/sph_harm_hc2dOm_lmax%d_ns%02d_r%d_b%02d-%-02d.npz' + output_name =(output_dir+'/sph_harm_hc2dOm_lmax%d_ns%02d_r%d_b%02d-%-02d.npz' % (lmax, nside, nreals, bestrange[0], bestrange[1]-1)) print('Saving npz file: ', output_name) np.savez(output_name, nsort=nsort, fidx=fidx, hc_ref=hc_ref, ss_shape=shape, moll_hc_best=moll_hc_best, Cl_best=Cl_best, nside=nside, lmax=lmax, fobs=fobs, split=split) - + # # ---- plot median Cl/C0 - + # print('Plotting Cl/C0 for median realizations') # fig = plot_ClC0_medians(fobs, Cl_best, lmax, nshow=(bestrange[1]-bestrange[0])) - # fig_name = (output_dir+'/sph_harm_hc2dOm_lmax%d_ns%02d_r%d_b%02d-%-02d.png' + # fig_name = (output_dir+'/sph_harm_hc2dOm_lmax%d_ns%02d_r%d_b%02d-%-02d.png' # % (lmax, nside, nreals, bestrange[0], bestrange[1]-1)) # fig.savefig(fig_name, dpi=300) @@ -415,17 +427,17 @@ def Cl_analytic_from_num(fobs_orb_edges, number, hs, realize = False, floor = Fa How many realizations to Poisson sample. floor : boolean Whether or not to round numbers down to nearest integers, if not realizing - + Returns ------- C0 : (F,R) or (F,) NDarray - C_0 + C_0 Cl : (F,R) or (F,) NDarray C_l>0 for arbitrary l using shot noise approximation """ df = np.diff(fobs_orb_edges) #: frequency bin widths - fc = kale.utils.midpoints(fobs_orb_edges) #: frequency-bin centers + fc = kale.utils.midpoints(fobs_orb_edges) #: frequency-bin centers # df = fobs_orb_widths[np.newaxis, np.newaxis, np.newaxis, :] # (M,Q,Z,F) NDarray # fc = fobs_orb_cents[np.newaxis, np.newaxis, np.newaxis, :] # (M,Q,Z,F) NDarray @@ -433,7 +445,7 @@ def Cl_analytic_from_num(fobs_orb_edges, number, hs, realize = False, floor = Fa # Poisson sample number in each bin if utils.isinteger(realize): - number = np.random.poisson(number[...,np.newaxis], + number = np.random.poisson(number[...,np.newaxis], size = (number.shape + (realize,))) df = df[...,np.newaxis] fc = fc[...,np.newaxis] @@ -455,7 +467,7 @@ def Cl_analytic_from_num(fobs_orb_edges, number, hs, realize = False, floor = Fa def strain_amp_at_bin_edges_redz(edges, redz=None): """ Calculate strain amplitude at bin edges, with final or initial redz. - + """ assert len(edges) == 4 assert np.all([np.ndim(ee) == 1 for ee in edges]) @@ -469,7 +481,7 @@ def strain_amp_at_bin_edges_redz(edges, redz=None): dc = +np.inf * np.ones_like(redz) sel = (redz > 0.0) dc[sel] = holo.cosmo.comoving_distance(redz[sel]).cgs.value - else: + else: redz = edges[2][np.newaxis,np.newaxis,:,np.newaxis] dc = holo.cosmo.comoving_distance(redz).cgs.value @@ -478,7 +490,7 @@ def strain_amp_at_bin_edges_redz(edges, redz=None): mr = (edges[1]) mc = utils.chirp_mass_mtmr(mt[:, np.newaxis], mr[np.newaxis, :]) mc = mc[:, :, np.newaxis, np.newaxis] - + # convert from observer-frame to rest-frame; still using frequency-bin centers fr = utils.frst_from_fobs(fc[np.newaxis, np.newaxis, np.newaxis, :], redz) @@ -488,7 +500,7 @@ def strain_amp_at_bin_edges_redz(edges, redz=None): def strain_amp_at_bin_centers_redz(edges, redz=None): """ Calculate strain amplitude at bin centers, with final or initial redz. - + """ assert len(edges) == 4 assert np.all([np.ndim(ee) == 1 for ee in edges]) @@ -518,7 +530,7 @@ def strain_amp_at_bin_centers_redz(edges, redz=None): mr = kale.utils.midpoints(edges[1]) mc = utils.chirp_mass_mtmr(mt[:, np.newaxis], mr[np.newaxis, :]) mc = mc[:, :, np.newaxis, np.newaxis] - + # convert from observer-frame to rest-frame; still using frequency-bin centers fr = utils.frst_from_fobs(fc[np.newaxis, np.newaxis, np.newaxis, :], redz) @@ -536,7 +548,7 @@ def Cl_analytic_from_dnum(edges, dnum, redz=None, realize=False): dN / [ dlog10M dq dz dlnf ] hs : (M,Q,Z,F) NDarray Strain amplitude of each M,q,z bin - + """ fobs_orb_edges = edges[-1] fobs_gw_edges = fobs_orb_edges * 2.0 @@ -556,7 +568,7 @@ def Cl_analytic_from_dnum(edges, dnum, redz=None, realize=False): # integrate over redshift numh2 = utils.trapz(numh2, edges[2], axis=2) # times dln(f) - numh2 = numh2 * np.diff(np.log(fobs_gw_edges)) + numh2 = numh2 * np.diff(np.log(fobs_gw_edges)) # integrate over dlog10(M) numh4 = utils.trapz(dnum*hs_edges**4, np.log10(edges[0]), axis=0) @@ -570,18 +582,18 @@ def Cl_analytic_from_dnum(edges, dnum, redz=None, realize=False): elif utils.isinteger(realize): # add reals axis hs_cents = strain_amp_at_bin_centers_redz(edges, redz)[...,np.newaxis] - df = df[:,np.newaxis] - fc = fc[:,np.newaxis] + df = df[:,np.newaxis] + fc = fc[:,np.newaxis] + - number = holo.sam_cython.integrate_differential_number_3dx1d(edges, dnum) shape = number.shape + (realize,) number = holo.gravwaves.poisson_as_needed(number[...,np.newaxis] * np.ones(shape)) - # numh2 = number * hs_cents**2 * np.diff(np.log(fobs_gw_edges))[:,np.newaxis] - # numh4 = number * hs_cents**4 * np.diff(np.log(fobs_gw_edges))[:,np.newaxis] - numh2 = number * hs_cents**2 - numh4 = number * hs_cents**4 + # numh2 = number * hs_cents**2 * np.diff(np.log(fobs_gw_edges))[:,np.newaxis] + # numh4 = number * hs_cents**4 * np.diff(np.log(fobs_gw_edges))[:,np.newaxis] + numh2 = number * hs_cents**2 + numh4 = number * hs_cents**4 else: err = "`realize` ({}) must be one of {{False, integer}}!".format(realize) raise ValueError(err) @@ -603,7 +615,7 @@ def Cl_analytic_from_dnum(edges, dnum, redz=None, realize=False): -def draw_analytic(ax, Cl, C0, fobs_gw_cents, color='tab:orange', label='Eq. 17 analytic', +def draw_analytic(ax, Cl, C0, fobs_gw_cents, color='tab:orange', label='Eq. 17 analytic', alpha=1, lw=2, ls='-.'): xx = fobs_gw_cents yy = Cl/C0 # (F,) @@ -615,7 +627,7 @@ def draw_reals(ax, Cl_many, C0_many, fobs_gw_cents, color='tab:orange', label= xx = fobs_gw_cents yy = Cl_many/C0_many # (F,R) if show_median: - ax.plot(xx, np.median(yy[:,:], axis=-1), color=color, lw=lw_median, alpha=0.75) #, label='median of samples, $l=%d$' % ll) + ax.plot(xx, np.median(yy[:,:], axis=-1), color=color, lw=lw_median, alpha=0.75) #, label='median of samples, $l=%d$' % ll) if show_ci: for pp in [50, 98]: percs = pp/2 @@ -623,7 +635,7 @@ def draw_reals(ax, Cl_many, C0_many, fobs_gw_cents, color='tab:orange', label= ax.fill_between(xx, *np.percentile(yy[:,:], percs, axis=-1), color=color, alpha=0.15) if show_reals: rr = 0 - ax.plot(xx, yy[:,rr], color=color, alpha=0.15, linestyle=ls_reals, + ax.plot(xx, yy[:,rr], color=color, alpha=0.15, linestyle=ls_reals, label = label) for rr in range(1, np.min([nshow, len(Cl_many[0])])): ax.plot(xx, yy[:,rr], color=color, alpha=0.15, linestyle=ls_reals) @@ -635,7 +647,7 @@ def draw_spk(ax, label='SP & K Rough Estimate'): def draw_bayes(ax, lmax, colors = ['k', 'b', 'r', 'g', 'c', 'm'], ms=8): xx_nihan = np.array([2.0, 4.0, 5.9, 7.9, 9.9]) *10**-9 # Hz - + ClC0_nihan = np.array([ [0.20216773, 0.14690035, 0.09676646, 0.07453352, 0.05500382, 0.03177427], [0.21201336, 0.14884939, 0.10545698, 0.07734305, 0.05257189, 0.03090662], @@ -643,10 +655,10 @@ def draw_bayes(ax, lmax, colors = ['k', 'b', 'r', 'g', 'c', 'm'], ms=8): [0.19788951, 0.15765126, 0.09615489, 0.07475364, 0.0527356 , 0.03113331], [0.20182648, 0.14745265, 0.09681202, 0.0746824 , 0.05503161, 0.0317012 ]]) for ll in range(lmax): - ax.plot(xx_nihan, ClC0_nihan[:,ll], - label = '$l=%d$' % (ll+1), + ax.plot(xx_nihan, ClC0_nihan[:,ll], + label = '$l=%d$' % (ll+1), color=colors[ll], marker='o', ms=ms) - + def draw_sim(ax, xx, Cl_best, lmax, nshow, show_ci=True, show_reals=True): yy = Cl_best[:,:,:,1:]/Cl_best[:,:,:,0,np.newaxis] # (B,F,R,l) @@ -664,25 +676,25 @@ def draw_sim(ax, xx, Cl_best, lmax, nshow, show_ci=True, show_reals=True): for bb in range(0,nshow): # if ll==0 and bb==0: # label = "individual best samples, median of realizations" - # else: + # else: label=None ax.plot(xx, yy[bb,:,ll], color=colors[ll], linestyle=':', alpha=0.1, linewidth=1, label=label) -def plot_ClC0_versions(fobs_gw_cents, spk=True, bayes=True, +def plot_ClC0_versions(fobs_gw_cents, spk=True, bayes=True, sim=True, Cl_best_sim=None, lmax_sim=None, analytic=False, Cl_analytic=None, C0_analytic=None, label_analytic='analytic', - anreals=False, Cl_anreals=None, C0_anreals=None, label_anreals=None, + anreals=False, Cl_anreals=None, C0_anreals=None, label_anreals=None, xmax = 1/YR, leg_anchor=(0,-0.15), leg_cols=3, legend=False): fig, ax = plot.figax(xlabel=plot.LABEL_GW_FREQUENCY_HZ, ylabel='$C_{\ell>0}/C_0$') if analytic: draw_analytic(ax, Cl_analytic, C0_analytic, fobs_gw_cents, label=label_analytic) if anreals: draw_reals(ax, Cl_anreals, C0_anreals, fobs_gw_cents, label=label_anreals) - + if bayes: draw_bayes(ax, lmax=6) if spk: draw_spk(ax, label='S-P & K') - if sim and (Cl_best_sim is not None) and (lmax_sim is not None): + if sim and (Cl_best_sim is not None) and (lmax_sim is not None): draw_sim(ax, fobs_gw_cents, Cl_best_sim, lmax_sim, show_ci=True, show_reals=True, nshow=10) # ax.set_ylim(10**-6, 10**0) # plot._twin_yr(ax, nano=False) diff --git a/holodeck/detstats.py b/holodeck/detstats.py index 317dab29..9f3ff18c 100644 --- a/holodeck/detstats.py +++ b/holodeck/detstats.py @@ -7,22 +7,32 @@ import numpy as np from scipy import special, integrate -from sympy import nsolve, Symbol import h5py import matplotlib.pyplot as plt import os from datetime import datetime import warnings - import holodeck as holo -from holodeck import utils, cosmo, log, plot, anisotropy -from holodeck.constants import MPC, YR +from holodeck import utils, log, plot # , cosmo, anisotropy +from holodeck.constants import YR from holodeck import cyutils -import hasasia.sensitivity as hsen -import hasasia.sim as hsim -import hasasia as has +try: + from sympy import nsolve, Symbol + import hasasia.sensitivity as hsen + import hasasia.sim as hsim + # import hasasia as has +except ImportError as err: + SUBMOD = "detstats" + log.error(f"Failed to import some packages used in `{SUBMOD}` submodule!") + log.exception(err) + log.error( + f"Some required packages for `{SUBMOD}` have been temporarily disabled in the " + "global 'requirements.txt' file, so they are not installed by default! Please install " + "the required packages manually for now, and feel free to raise a github issue." + ) + raise GAMMA_RHO_GRID_PATH = '/Users/emigardiner/GWs/holodeck/output/rho_gamma_grids' HC_REF15_10YR = 11.2*10**-15 diff --git a/holodeck/gps/__init__.py b/holodeck/gps/__init__.py index daa31812..9cce51e9 100644 --- a/holodeck/gps/__init__.py +++ b/holodeck/gps/__init__.py @@ -1,6 +1,21 @@ """ """ +from holodeck import log + +try: + import schwimmbad + import emcee + import george +except ImportError as err: + log.error("Failed to import packages used in `gps` submodule!") + log.exception(err) + log.error( + "Some required packages for the `gps` submodule have been temporarily disabled in the " + "global 'requirements.txt' file, so they are not installed by default! Please install " + "the required packages manually for now, and feel free to raise a github issue." + ) + raise # ---- Import submodules diff --git a/requirements.txt b/requirements.txt index a33db432..0ab402e5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,10 +10,10 @@ scipy psutil tqdm cython<3.0.0 -schwimmbad -emcee -george -sympy # used in the detstats submodule -hasasia # used in the detstats submodule -healpy # used in the detstats submodule +# schwimmbad # used in the `gps` submodule +# emcee # used in the `gps` submodule +# george # used in the `gps` submodule +# sympy # used in the detstats submodule +# hasasia # used in the detstats submodule +# healpy # used in the anisotropy submodule # mpi4py From 8bfddfe3958157d8ea968cbd7c07df6409d23256 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 1 May 2024 10:52:11 -0700 Subject: [PATCH 46/48] BUG: use codecov token stored to github repository secrets --- .github/workflows/coverage.yaml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/coverage.yaml b/.github/workflows/coverage.yaml index b63a60a0..b35adaa7 100644 --- a/.github/workflows/coverage.yaml +++ b/.github/workflows/coverage.yaml @@ -64,9 +64,10 @@ jobs: pytest -v --cov=./ --cov-report term --cov-report=xml --color=yes - name: Upload coverage to Codecov - uses: codecov/codecov-action@v3 + uses: codecov/codecov-action@v4 with: name: codecov-${{ matrix.os }}-py${{ matrix.python-version }} + token: ${{ secrets.CODECOV_TOKEN }} # required env_vars: OS,PYTHON fail_ci_if_error: true verbose: true From 8a93b1c7edd4ca31e3c89cc24b2d188a353d0834 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 1 May 2024 12:34:45 -0700 Subject: [PATCH 47/48] BUG: fix issue with how 'param_samples' was being compared to None. If the array is a normal array then the '==' comparison fails. --- holodeck/librarian/combine.py | 6 +++--- holodeck/librarian/lib_tools.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/holodeck/librarian/combine.py b/holodeck/librarian/combine.py index 8a084ffe..ac0b7d8b 100644 --- a/holodeck/librarian/combine.py +++ b/holodeck/librarian/combine.py @@ -140,11 +140,11 @@ def sam_lib_combine( log.info(f"loaded param space: {pspace} from '{pspace_fname}'") param_names = pspace.param_names - param_samples = pspace.param_samples + param_samples = pspace.param_samples[()] # Standard Library: vary all parameters together if library: - if param_samples == None: # noqa : use `== None` to match arrays + if param_samples is None: log.error(f"`library`={library} but `param_samples`={param_samples}`") err = f"`library` is True, but {path_output} looks like it's a domain." raise DomainNotLibraryError(err) @@ -155,7 +155,7 @@ def sam_lib_combine( # Domain: Vary only one parameter at a time to explore the domain else: err = f"Expected 'domain' but `param_samples`={param_samples} is not `None`!" - assert param_samples == None, err # noqa : use `== None` to match arrays + assert param_samples is None, err ndim = pspace.nparameters # for 'domain' simulations, this is the number of samples in each dimension nsamp_dim = args.nsamples diff --git a/holodeck/librarian/lib_tools.py b/holodeck/librarian/lib_tools.py index 0042a6a1..616a6be6 100644 --- a/holodeck/librarian/lib_tools.py +++ b/holodeck/librarian/lib_tools.py @@ -329,7 +329,7 @@ def from_save(cls, fname, log=None): pspace_class = cls # construct instance with dummy/temporary values (which will be overwritten) - if data['param_samples'] == None: # noqa : use ``== None`` to match arrays + if data['param_samples'][()] is None: nsamples = None nparameters = None else: From 0be249e027284e1d53d9aec09e764ae5c04f6655 Mon Sep 17 00:00:00 2001 From: Luke Zoltan Kelley Date: Wed, 1 May 2024 12:35:26 -0700 Subject: [PATCH 48/48] Bump to v1.6 --- CHANGELOG.md | 2 ++ holodeck/version.txt | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5dee824a..72a5b303 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,8 @@ ## Current +## v1.6 - 2024/05/01 + * DEPRECATIONS: * `holodeck.librarian.py` ==> `holodeck.librarian.lib_tools` * Rename submodule. All components remain the same. All `lib_tools` elements are now also imported into the `librarian` namespace. i.e. elements like `holodeck.librarian.lib_tools._Param_Space` will now be accessible via `holodeck.librarian._Param_Space`. diff --git a/holodeck/version.txt b/holodeck/version.txt index 4cda8f19..810ee4e9 100644 --- a/holodeck/version.txt +++ b/holodeck/version.txt @@ -1 +1 @@ -1.5.2 +1.6