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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 641a5596..72a5b303 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,35 @@ ## 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`. + * 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 + * `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. + +* 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 + * 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/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/__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 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/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/host_relations.py b/holodeck/host_relations.py index 22b9b8be..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]) @@ -779,7 +778,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 @@ -843,7 +842,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: + except: redz = redz zmamp = self._mamp * (1.0 + redz)**self._zplaw mbh = _log10_relation(mbulge, zmamp, self._mplaw, scatter_dex, x0=self._mref) diff --git a/holodeck/librarian/__init__.py b/holodeck/librarian/__init__.py index 0c18bfde..e0380401 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,7 +39,11 @@ """ -__version__ = "1.1" +# ============================================================================== +# ==== Submodule Definitions ==== +# ============================================================================== + +__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. @@ -49,18 +53,41 @@ 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" +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" + +#: 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) + +# ============================================================================== +# ==== Import Submodule Components ==== +# ============================================================================== + +from . lib_tools import * # noqa +from . import gen_lib # noqa -from . libraries import ( # noqa - _Param_Space, _Param_Dist, - PD_Uniform, PD_Normal, - run_model, load_pspace_from_path, -) +# 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 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 ded49fd6..ac0b7d8b 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 """ @@ -14,7 +19,10 @@ import holodeck as holo import holodeck.librarian -from holodeck.librarian import libraries +import holodeck.librarian.gen_lib +from holodeck.librarian import ( + lib_tools, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS, DomainNotLibraryError +) def main(): @@ -23,6 +31,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 +44,8 @@ def main(): log.exception(err) raise RuntimeError(err) + # ---- Setup and parse command-line arguments + parser = argparse.ArgumentParser() parser.add_argument( @@ -48,17 +60,33 @@ 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) - sam_lib_combine(path, log, recreate=args.recreate, gwb_only=args.gwb) + # ---- Combine library files + + 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 @@ -73,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 ------- @@ -85,11 +115,14 @@ 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 - lib_path = libraries.get_sam_lib_fname(path_output, gwb_only) + lib_path = lib_tools.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.") @@ -100,38 +133,43 @@ 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) + 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}'") param_names = pspace.param_names - param_samples = pspace.param_samples - nsamp, ndim = param_samples.shape - log.debug(f"{nsamp=}, {ndim=}, {param_names=}") + param_samples = pspace.param_samples[()] + + # Standard Library: vary all parameters together + if library: + 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) + + 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 is None, err + 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 + # 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 - 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=}") @@ -149,28 +187,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, + 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)=}") - param_samples[bad_files] = np.nan + if library: + param_samples[bad_files] = np.nan # ---- Save to concatenated output file ---- @@ -178,14 +217,35 @@ 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) + + 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: + 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') @@ -202,19 +262,28 @@ 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 -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 ------- @@ -231,10 +300,12 @@ 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): - temp_fname = libraries._get_sim_fname(path_sims, ii) + 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) @@ -248,6 +319,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: @@ -287,7 +364,10 @@ 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, param_samples, + log, library, +): """Load data from all individual simulation files. Arguments @@ -297,29 +377,36 @@ 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 + param_samples + log : ``logging.Logger`` Logging instance. + library : bool + Whether this is a standard library or a domain exploration. """ 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): - fname = libraries._get_sim_fname(path_sims, pnum) + for pnum in tqdm.trange(nsamp_all): + 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): 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 @@ -331,6 +418,16 @@ 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: + #! 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: gwb[pnum, :, :] = temp['gwb'][...] @@ -345,7 +442,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/fit_spectra.py b/holodeck/librarian/fit_spectra.py index b5138d7c..e57f0432 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 @@ -36,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 @@ -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 @@ -101,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) @@ -111,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 @@ -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) @@ -267,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 [] diff --git a/holodeck/librarian/gen_lib.py b/holodeck/librarian/gen_lib.py index 0a3b489a..29cafe9e 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.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:: @@ -11,34 +12,40 @@ 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 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, + lib_tools, ARGS_CONFIG_FNAME, PSPACE_DOMAIN_EXTREMA, DIRNAME_LIBRARY_SIMS, DIRNAME_DOMAIN_SIMS ) -# from holodeck.sams import sam_cyutils -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 = [] +comm = None + def main(): # noqa : ignore complexity warning """Parent method for generating libraries from the command-line. @@ -61,14 +68,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 @@ -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 @@ -98,8 +97,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) @@ -107,98 +110,114 @@ 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_dir(log, args.output, space_class) - log.warning(f"resume={args.resume} :: Loaded param-space save from {space_fname}") + 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 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 simulations for all processes + + # Domain: Vary only one parameter at a time to explore the domain + if args.domain: + 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) + 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: + 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" + 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" ) - 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 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() - beg = datetime.now() - log.info(f"beginning tasks at {beg}") - failures = 0 # ---- iterate over each processors' jobs - for par_num in iterator: + beg = datetime.now() + log.info(f"beginning tasks at {beg}") + failures = 0 + num_done = 0 + for sim_num in iterator: + log.info(f"{comm.rank=} {sim_num=}") + + # 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) - 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) + rv, _sim_fname = run_sam_at_pspace_params(args, space, sim_num, params) - 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()}") @@ -208,78 +227,95 @@ 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 -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.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. - NOTE: this should be improved. - space : :class:`holodeck.librarian.libraries._Param_space` instance + Arguments from the ``gen_lib_sams.py`` script. + 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. + 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 # ---- 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 = lib_tools._get_sim_fname(args.output_sims, pnum, library=library_flag) beg = datetime.now() log.info(f"{pnum=} :: {sim_fname=} beginning at {beg}") 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 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( + 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, log=log, ) + data['params'] = np.array([params[pn] for pn in space.param_names]) + data['param_names'] = space.param_names rv = True except Exception as err: 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 @@ -339,8 +375,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, @@ -352,8 +390,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) @@ -372,10 +411,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)) @@ -405,7 +451,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 @@ -421,6 +472,105 @@ def _setup_argparse(*args, **kwargs): return args +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 + + # `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: + # 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 + + +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(): + # 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 + 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) + + args.log.warning(f"Saved to {fname} - {holo.utils.get_file_size(fname)}") + + return fname + + +def load_config_from_path(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. @@ -458,14 +608,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 @@ -619,3 +767,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 + #! ---------- diff --git a/holodeck/librarian/libraries.py b/holodeck/librarian/lib_tools.py similarity index 86% rename from holodeck/librarian/libraries.py rename to holodeck/librarian/lib_tools.py index 3e3d66b1..616a6be6 100644 --- a/holodeck/librarian/libraries.py +++ b/holodeck/librarian/lib_tools.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 = [ @@ -25,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. @@ -59,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 @@ -123,7 +128,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,14 +211,47 @@ 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 + --------- + 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 + # @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 + --------- + 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): @@ -234,7 +272,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 @@ -244,16 +283,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, ) @@ -290,8 +329,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'][()] is None: + 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: @@ -303,6 +347,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 @@ -354,22 +407,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 ------- @@ -381,11 +435,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 @@ -884,7 +944,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. @@ -928,9 +988,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, allow_pickle=True)['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 @@ -943,14 +1015,25 @@ 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 -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") diff --git a/holodeck/librarian/param_spaces.py b/holodeck/librarian/param_spaces.py index 600b7a0b..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 +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: @@ -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 @@ -158,7 +158,7 @@ class _PS_Astro_Strong(_Param_Space): """ - __version__ = "0.1" + __version__ = "0.2" DEFAULTS = dict( # Hardening model (phenom 2PL) @@ -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,13 +196,17 @@ 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_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 - 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'], @@ -239,19 +243,26 @@ 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( - 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, @@ -269,7 +280,8 @@ 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), + 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 @@ -295,9 +307,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.1, 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, @@ -312,7 +330,26 @@ 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), + 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, @@ -367,13 +404,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.1, 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, @@ -385,8 +442,10 @@ 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': PS_Astro_Strong_MMBulge, + 'PS_Astro_Strong_MMBulge_BFrac': PS_Astro_Strong_MMBulge_BFrac, } diff --git a/holodeck/librarian/param_spaces_classic.py b/holodeck/librarian/param_spaces_classic.py index 3509d319..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 @@ -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), ] diff --git a/holodeck/librarian/tests/test_libraries__param_space.py b/holodeck/librarian/tests/test_lib_tools__param_space.py similarity index 80% rename from holodeck/librarian/tests/test_libraries__param_space.py rename to holodeck/librarian/tests/test_lib_tools__param_space.py index aa7c54bd..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 @@ -9,9 +9,11 @@ from holodeck import sams, host_relations, hardening, librarian from holodeck.constants import GYR -from holodeck.librarian.libraries import ( +from holodeck.librarian import ( + param_spaces_dict +) +from holodeck.librarian.lib_tools import ( _Param_Space, PD_Uniform, - ) MMB_MAMP_LOG10_EXTR = [+7.5, +9.5] @@ -203,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 @@ -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 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/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 diff --git a/holodeck/sams/sam.py b/holodeck/sams/sam.py index 1e00841a..f8761dbc 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 @@ -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 @@ -357,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 @@ -465,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 @@ -479,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)[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], @@ -499,33 +514,65 @@ 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) - # 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] - redz_evo = utils.redz_after(times_tot, redz=rz) + # ~~~~ 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) + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + # 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) + + # ---- 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: + tlbk -= self._gmt_time[:, :, :, np.newaxis] + + # (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) 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 + # (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, 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]] # (M*Q*Z, X) 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,)) + coal = (redz_final > 0.0) frst_orb = fobs_orb * (1.0 + redz_final) frst_orb[frst_orb < 0.0] = 0.0 @@ -577,7 +624,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 +638,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 diff --git a/holodeck/sams/sam_cyutils.pyx b/holodeck/sams/sam_cyutils.pyx index 704e7ada..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 @@ -452,6 +453,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 +482,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 +570,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 @@ -640,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 @@ -654,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 @@ -681,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) @@ -721,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 @@ -849,6 +875,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 diff --git a/holodeck/sams/tests/test_sam.py b/holodeck/sams/tests/test_sam.py index e6889247..1a1874a6 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,150 @@ 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 + + # (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) + + # (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 + + +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/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 diff --git a/holodeck/utils.py b/holodeck/utils.py index 790583cc..55465d94 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) @@ -1042,7 +1043,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 +1066,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!") @@ -1773,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): 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 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 68cc0248..12444622 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.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", - " 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", 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", diff --git a/notebooks/semi-analytic-models.ipynb b/notebooks/semi-analytic-models.ipynb index a93b6acd..24fcc638 100644 --- a/notebooks/semi-analytic-models.ipynb +++ b/notebooks/semi-analytic-models.ipynb @@ -411,7 +411,7 @@ "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)" ] }, { @@ -466,9 +466,19 @@ ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "id": "32", "metadata": {}, + "outputs": [], + "source": [ + "nden.shape, diff_num_gw.shape, number_gw.shape" + ] + }, + { + "cell_type": "markdown", + "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": "33", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -507,7 +517,7 @@ }, { "cell_type": "markdown", - "id": "34", + "id": "35", "metadata": {}, "source": [ "### Self-Consistent Binary Evolution (Phenomenological Model)" @@ -516,7 +526,7 @@ { "cell_type": "code", "execution_count": null, - "id": "35", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -532,7 +542,7 @@ { "cell_type": "code", "execution_count": null, - "id": "36", + "id": "37", "metadata": {}, "outputs": [], "source": [ @@ -566,7 +576,7 @@ }, { "cell_type": "markdown", - "id": "37", + "id": "38", "metadata": {}, "source": [ "## Gravitational Waves" @@ -574,7 +584,7 @@ }, { "cell_type": "markdown", - "id": "38", + "id": "39", "metadata": {}, "source": [ "### Compare GWB and CW" @@ -583,7 +593,7 @@ { "cell_type": "code", "execution_count": null, - "id": "39", + "id": "40", "metadata": {}, "outputs": [], "source": [ @@ -604,7 +614,7 @@ { "cell_type": "code", "execution_count": null, - "id": "40", + "id": "41", "metadata": {}, "outputs": [], "source": [ @@ -628,7 +638,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "41", + "id": "42", "metadata": {}, "source": [ "Calculate the distribution of GWB Amplitudes at 1/yr" @@ -637,7 +647,7 @@ { "cell_type": "code", "execution_count": null, - "id": "42", + "id": "43", "metadata": {}, "outputs": [], "source": [ @@ -680,7 +690,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "43", + "id": "44", "metadata": {}, "source": [ "## Plot GWB Amplitude Distribution vs. M-MBulge parameters" @@ -689,7 +699,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "44", + "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": "45", + "id": "46", "metadata": {}, "outputs": [], "source": [ @@ -727,7 +737,7 @@ { "attachments": {}, "cell_type": "markdown", - "id": "46", + "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": "47", + "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": { 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 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() + + diff --git a/scripts/run_holodeck_lib_gen.sh b/scripts/run_holodeck_lib_gen.sh index 1b37d50c..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``. @@ -9,25 +10,32 @@ # 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/ # +# 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. # # 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 +144,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" 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', + ], }, )