diff --git a/.gitignore b/.gitignore index a33ef616..dc4b6322 100644 --- a/.gitignore +++ b/.gitignore @@ -62,6 +62,9 @@ instance/ # Scrapy stuff: .scrapy +# Intermediate documentation build files +docs/source/api_ref/ + # Sphinx documentation docs/_build/ diff --git a/docs/docgen.sh b/docs/docgen.sh index 8abb405f..2f2448dc 100755 --- a/docs/docgen.sh +++ b/docs/docgen.sh @@ -18,21 +18,24 @@ cd "$(dirname "$0")" # Remove previously created symlink file, if it exists path_symlink_file=$PWD/$SYMLINK_FILE -if [ -f "$path_symlink_file" ] ; then - rm "$path_symlink_file" +if [ -L "$path_symlink_file" ] || [ -f "$path_symlink_file" ] ; then + rm -f "$path_symlink_file" fi +rm -rf source/apidoc_modules/ +python3 generate_api_ref.py + # Install requirements (yes, this is also needed for building the docs) -pip install -r requirements.txt +# pip install -r requirements.txt # Cleanup previous material make clean # Generate apidoc automatically created documentation -sphinx-apidoc -e -P -T -f -o ./source/apidoc_modules ../holodeck +# sphinx-apidoc -e -P -T -f -o ./source/apidoc_modules ../holodeck # Run the normal sphinx build -make html +make html SPHINXOPTS="-j 1" # create convenience symlink ln -s $PWD/build/html/index.html $path_symlink_file # move back to previous directory cd - -# Done! \ No newline at end of file +# Done! diff --git a/docs/generate_api_ref.py b/docs/generate_api_ref.py new file mode 100755 index 00000000..5f0ad20a --- /dev/null +++ b/docs/generate_api_ref.py @@ -0,0 +1,189 @@ +#!/usr/bin/env python3 +from astropy.coordinates import matrix_utilities +import os +import sys + +# Define directories relative to this script +DOCS_DIR = os.path.dirname(os.path.abspath(__file__)) +SOURCE_DIR = os.path.join(DOCS_DIR, "source", "api_ref") +REPO_ROOT = os.path.dirname(DOCS_DIR) +CODE_DIR = os.path.join(REPO_ROOT, "holodeck") + +# 1. Global Ignore Configuration +IGNORE_DIRS = {"__pycache__", "tests", "data", "source", "build"} +IGNORE_FILES = { + "__init__.py", + "param_spaces_old.py", + "relations.py", +} + +# 2. Categorization Mapping for holodeck.rst +CATEGORIES = { + "Core Package Foundations": [ + "holodeck.constants", + "holodeck.cyutils", + "holodeck.discrete_cyutils", + "holodeck.extensions", + "holodeck.logger", + ], + "Astrophysics & Population Modeling": [ + "holodeck.accretion", + "holodeck.discrete", + "holodeck.galaxy_profiles", + "holodeck.host_relations", + "holodeck.pop_observational", + "holodeck.sams", + ], + "Gravitational Waves & Evolution": [ + "holodeck.anisotropy", + "holodeck.gravwaves", + "holodeck.hardening", + "holodeck.single_sources", + ], + "Analysis, Libraries & Tools": [ + "holodeck.detstats", + "holodeck.gps", + "holodeck.librarian", + "holodeck.observations", + "holodeck.param_spaces", + "holodeck.plot", + "holodeck.utils", + ], +} + +def clean_and_prepare(): + """Ensure a fresh, clean api_ref directory target.""" + print(f"Cleaning target directory: {SOURCE_DIR}") + if os.path.exists(SOURCE_DIR): + import shutil + shutil.rmtree(SOURCE_DIR) + os.makedirs(SOURCE_DIR, exist_ok=True) + +def write_rst_file(module_name, is_package=False, submodules=None, pure_structure=False): + """Generate individual .rst files for modules and packages cleanly.""" + file_path = os.path.join(SOURCE_DIR, f"{module_name}.rst") + + title = f"{module_name} package" if is_package else f"{module_name} module" + underline = "=" * len(title) + + content = [ + title, + underline, + "" + ] + + # Only try to extract docstrings if it's an importable Python file/package + if not pure_structure: + content.extend([ + f".. automodule:: {module_name}", + " :members:", + " :undoc-members:", + " :show-inheritance:", + " :member-order: bysource", + "" + ]) + + # If it's a directory package containing sub-items, link them together + if is_package and submodules: + content.extend([ + "Submodules", + "----------", + "", + ".. toctree::", + " :maxdepth: 4", + "" + ]) + for sub in sorted(submodules): + content.append(f" {sub}") + content.append("") + + with open(file_path, "w") as f: + f.write("\n".join(content)) + + +def scan_codebase(): + """Traverse holodeck codebase and construct the .rst files dynamically.""" + print(f"Scanning codebase: {CODE_DIR}") + + # Track items discovered to link in top-level maps + top_level_modules = set() + + for root, dirs, files in os.walk(CODE_DIR): + # Apply directory filters in-place + dirs[:] = [d for d in dirs if d not in IGNORE_DIRS] + + # Determine package identity path + rel_path = os.path.relpath(root, REPO_ROOT) + package_prefix = rel_path.replace(os.sep, ".") + +# If we are inside a subpackage directory (like holodeck/sams) + if root != CODE_DIR: + # ONLY treat this directory as an importable package if it contains an __init__.py + has_init = "__init__.py" in files + + submodules = [] + for f in files: + if (f.endswith(".py") or f.endswith(".pyx")) and f not in IGNORE_FILES: + mod_base = os.path.splitext(f)[0] + full_mod_name = f"{package_prefix}.{mod_base}" + write_rst_file(full_mod_name, is_package=False) + submodules.append(full_mod_name) + + if submodules and has_init: + write_rst_file(package_prefix, is_package=True, submodules=submodules) + elif submodules and not has_init: + # If there's no __init__.py, don't use .. automodule:: on the directory path itself + write_rst_file(package_prefix, is_package=True, submodules=submodules, pure_structure=True) + continue + + # Processing the top-level holodeck/ root folder items + for f in files: + if (f.endswith(".py") or f.endswith(".pyx")) and f not in IGNORE_FILES: + mod_base = os.path.splitext(f)[0] + full_mod_name = f"holodeck.{mod_base}" + write_rst_file(full_mod_name, is_package=False) + top_level_modules.add(full_mod_name) + + for d in dirs: + # Mark top level tracking directory targets (like holodeck.sams) + top_level_modules.add(f"holodeck.{d}") + +def generate_master_index(): + """Regenerate the unified holodeck.rst entry point layout.""" + master_path = os.path.join(SOURCE_DIR, "holodeck.rst") + print(f"Generating master index layout: {master_path}") + + content = [ + "======================", + "holodeck API Reference", + "======================", + "", + ".. automodule:: holodeck", + " :members:", + " :private-members:", + " :undoc-members:", + " :show-inheritance:", + " :member-order: bysource", + "" + ] + + # Iterate clean structural layout blocks + for category_name, modules in CATEGORIES.items(): + content.extend([ + ".. toctree::", + " :maxdepth: 2", + f" :caption: {category_name}", + "" + ]) + for mod in sorted(modules): + content.append(f" {mod}") + content.extend(["", ""]) + + with open(master_path, "w") as f: + f.write("\n".join(content)) + +if __name__ == "__main__": + clean_and_prepare() + scan_codebase() + generate_master_index() + print("✨ API Reference Generation Complete!") diff --git a/docs/requirements.txt b/docs/requirements.txt index 64398f04..a2456780 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,8 +1,7 @@ -r ../requirements-dev.txt ipykernel nbsphinx -numpydoc sphinx sphinx-math-dollar sphinx-rtd-theme -# sphinxcontrib-bibtex \ No newline at end of file +# sphinxcontrib-bibtex diff --git a/docs/source/api_ref/holodeck.discrete.evolution.rst b/docs/source/api_ref/holodeck.discrete.evolution.rst deleted file mode 100644 index 901e5888..00000000 --- a/docs/source/api_ref/holodeck.discrete.evolution.rst +++ /dev/null @@ -1,11 +0,0 @@ -=========================== -holodeck.discrete.evolution -=========================== - -.. automodule:: holodeck.discrete.evolution - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource - \ No newline at end of file diff --git a/docs/source/api_ref/holodeck.discrete.population.rst b/docs/source/api_ref/holodeck.discrete.population.rst deleted file mode 100644 index be5fdb5f..00000000 --- a/docs/source/api_ref/holodeck.discrete.population.rst +++ /dev/null @@ -1,11 +0,0 @@ -============================ -holodeck.discrete.population -============================ - -.. automodule:: holodeck.discrete.population - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource - \ No newline at end of file diff --git a/docs/source/api_ref/holodeck.discrete.rst b/docs/source/api_ref/holodeck.discrete.rst deleted file mode 100644 index 56ba54a6..00000000 --- a/docs/source/api_ref/holodeck.discrete.rst +++ /dev/null @@ -1,16 +0,0 @@ -======================== -holodeck.discrete module -======================== - -.. automodule:: holodeck.discrete - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource - -.. toctree:: - :maxdepth: 2 - - holodeck.discrete.evolution - holodeck.discrete.population diff --git a/docs/source/api_ref/holodeck.galaxy_profiles.rst b/docs/source/api_ref/holodeck.galaxy_profiles.rst deleted file mode 100644 index 51203ad7..00000000 --- a/docs/source/api_ref/holodeck.galaxy_profiles.rst +++ /dev/null @@ -1,10 +0,0 @@ -======================== -holodeck.galaxy_profiles -======================== - -.. automodule:: holodeck.galaxy_profiles - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource diff --git a/docs/source/api_ref/holodeck.gravwaves.rst b/docs/source/api_ref/holodeck.gravwaves.rst deleted file mode 100644 index 0a966031..00000000 --- a/docs/source/api_ref/holodeck.gravwaves.rst +++ /dev/null @@ -1,11 +0,0 @@ -================== -holodeck.gravwaves -================== - -.. automodule:: holodeck.gravwaves - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource - \ No newline at end of file diff --git a/docs/source/api_ref/holodeck.host_relations.rst b/docs/source/api_ref/holodeck.host_relations.rst deleted file mode 100644 index eee1b4a3..00000000 --- a/docs/source/api_ref/holodeck.host_relations.rst +++ /dev/null @@ -1,10 +0,0 @@ -======================= -holodeck.host_relations -======================= - -.. automodule:: holodeck.host_relations - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource diff --git a/docs/source/api_ref/holodeck.librarian.lib_tools.rst b/docs/source/api_ref/holodeck.librarian.lib_tools.rst deleted file mode 100644 index 160b12fe..00000000 --- a/docs/source/api_ref/holodeck.librarian.lib_tools.rst +++ /dev/null @@ -1,11 +0,0 @@ -============================ -holodeck.librarian.lib_tools -============================ - -.. automodule:: holodeck.librarian.lib_tools - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource - \ No newline at end of file diff --git a/docs/source/api_ref/holodeck.librarian.rst b/docs/source/api_ref/holodeck.librarian.rst deleted file mode 100644 index 9ff9eb56..00000000 --- a/docs/source/api_ref/holodeck.librarian.rst +++ /dev/null @@ -1,21 +0,0 @@ -========================= -holodeck.librarian module -========================= - -.. automodule:: holodeck.librarian - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource - -.. toctree:: - :maxdepth: 1 - - holodeck.librarian.combine - holodeck.librarian.fit_spectra - holodeck.librarian.gen_lib - holodeck.librarian.lib_tools - holodeck.librarian.param_spaces - holodeck.librarian.param_spaces_classic - holodeck.librarian.posterior_populations diff --git a/docs/source/api_ref/holodeck.rst b/docs/source/api_ref/holodeck.rst deleted file mode 100644 index cdf41c5e..00000000 --- a/docs/source/api_ref/holodeck.rst +++ /dev/null @@ -1,10 +0,0 @@ -====================== -holodeck API Reference -====================== - -.. automodule:: holodeck - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource diff --git a/docs/source/api_ref/holodeck.sams.components.rst b/docs/source/api_ref/holodeck.sams.components.rst deleted file mode 100644 index f235b8a1..00000000 --- a/docs/source/api_ref/holodeck.sams.components.rst +++ /dev/null @@ -1,10 +0,0 @@ -======================== -holodeck.sams.components -======================== - -.. automodule:: holodeck.sams.components - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource diff --git a/docs/source/api_ref/holodeck.sams.sam_cyutils.rst b/docs/source/api_ref/holodeck.sams.sam_cyutils.rst deleted file mode 100644 index 168ae0ab..00000000 --- a/docs/source/api_ref/holodeck.sams.sam_cyutils.rst +++ /dev/null @@ -1,10 +0,0 @@ -========================= -holodeck.sams.sam_cyutils -========================= - -.. automodule:: holodeck.sams.sam_cyutils - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource diff --git a/docs/source/api_ref/holodeck.utils.rst b/docs/source/api_ref/holodeck.utils.rst deleted file mode 100644 index 3e1ceeb6..00000000 --- a/docs/source/api_ref/holodeck.utils.rst +++ /dev/null @@ -1,10 +0,0 @@ -============== -holodeck.utils -============== - -.. automodule:: holodeck.utils - :members: - :private-members: - :undoc-members: - :show-inheritance: - :member-order: bysource diff --git a/docs/source/biblio.rst b/docs/source/biblio.rst index d38f8100..93658239 100644 --- a/docs/source/biblio.rst +++ b/docs/source/biblio.rst @@ -73,7 +73,36 @@ References ========== These are provided here for easy copy-and-paste usage in other files. -.. [Behroozi2013] : Behroozi, Wechsler & Conroy 2013. ApJ, 770, 1. + +.. [N15NP] Afzal et al. (2023), ApJL, 951, L11. + The NANOGrav 15 yr Data Set: Search for Signals from New Physics + https://ui.adsabs.harvard.edu/abs/2023ApJ...951L..11A/abstract + +.. [N15GWB] Agazie et al. (2023), ApJL, 951, L8. + The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background + https://ui.adsabs.harvard.edu/abs/2023ApJ...951L...8A/abstract + +.. [N15data] Agazie et al. (2023), ApJL, 951, L9. + The NANOGrav 15 yr Data Set: Observations and Timing of 68 Millisecond Pulsars + https://ui.adsabs.harvard.edu/abs/2023ApJ...951L...9A/abstract + +.. [N15detchar] Agazie et al. (2023), ApJL, 951, L10. + The NANOGrav 15 yr Data Set: Detector Characterization and Noise Budget + https://ui.adsabs.harvard.edu/abs/2023ApJ...951L..10A/abstract + +.. [N15CWs] Agazie et al. (2023), ApJL, 951, L50. + The NANOGrav 15 yr Data Set: Continuous Wave Signals from Individual Spinning Supermassive Black Holes + https://ui.adsabs.harvard.edu/abs/2023ApJ...951L..50A/abstract + +.. [N15astro] Agazie et al. (2023), ApJL, 952, L37. + The NANOGrav 15 yr Data Set: Constraints on Supermassive Black Hole Binaries from the Gravitational-wave Background + https://ui.adsabs.harvard.edu/abs/2023ApJ...952L..37A/abstract + +.. [N15anisotropy] Agazie et al. (2023), ApJL, 956, L3. + The NANOGrav 15 yr Data Set: Anisotropy and Spatial Variation in the Gravitational-wave Background + https://ui.adsabs.harvard.edu/abs/2023ApJ...956L...3A/abstract + +.. [Behroozi2013] Behroozi, Wechsler & Conroy 2013. ApJ, 770, 1. The Average Star Formation Histories of Galaxies in Dark Matter Halos from z = 0-8 https://ui.adsabs.harvard.edu/abs/2013ApJ...770...57B/abstract @@ -143,6 +172,10 @@ These are provided here for easy copy-and-paste usage in other files. Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies https://ui.adsabs.harvard.edu/abs/2013ARA%26A..51..511K/abstract +.. [Leja2020] Leja et al. 2020, ApJ, 893, 111. + A New Census of the 0.2 < z < 3.0 Universe. I. The Stellar Mass Function + https://ui.adsabs.harvard.edu/abs/2020ApJ...893..111L/abstract + .. [MM2013] McConnell & Ma 2013. ApJ, 764, 2. Revisiting the Scaling Relations of Black Hole Masses and Host Galaxy Properties https://ui.adsabs.harvard.edu/abs/2013ApJ...764..184M/abstract @@ -151,7 +184,7 @@ These are provided here for easy copy-and-paste usage in other files. A Universal Density Profile from Hierarchical Clustering https://ui.adsabs.harvard.edu/abs/1997ApJ...490..493N/abstract -.. [Nelson2015] Nelson et al. (2015), A&C, 13,. +.. [Nelson2015] Nelson et al. 2015, A&C, 13,. The illustris simulation: Public data release https://ui.adsabs.harvard.edu/abs/2015A&C....13...12N @@ -167,10 +200,14 @@ These are provided here for easy copy-and-paste usage in other files. The dynamical evolution of massive black hole binaries I. Hardening in a fixed stellar background https://ui.adsabs.harvard.edu/abs/1996NewA....1...35Q/abstract -.. [Rodriguez-Gomez2015] : Rodriguez-Gomez et al. (2015), MNRAS, 449, 1. +.. [Rodriguez-Gomez2015] Rodriguez-Gomez et al. 2015, MNRAS, 449, 1. The merger rate of galaxies in the Illustris simulation: a comparison with observations and semi-empirical models https://ui.adsabs.harvard.edu/abs/2015MNRAS.449...49R +.. [Rosado2015] Rosado, Sesana, & Gair (2015), MNRAS, 451, 2417. + Expected properties of the first gravitational wave signal detected with pulsar timing arrays + https://ui.adsabs.harvard.edu/abs/2015MNRAS.451.2417R + .. [Sesana2004] Sesana, Haardt, Madau, & Volonteri 2004. ApJ, 611, 2. astro-ph/0401543. Low-Frequency Gravitational Radiation from Coalescing Massive Black Hole Binaries in Hierarchical Cosmologies http://adsabs.harvard.edu/abs/2004ApJ...611..623S @@ -189,15 +226,19 @@ These are provided here for easy copy-and-paste usage in other files. Implications for Gravitational Wave Observations https://ui.adsabs.harvard.edu/abs/2010ApJ...719..851S/abstract -.. [Sijacki2015] Sijacki et al. (2015), MNRAS, 452, 1. +.. [Sijacki2015] Sijacki et al. 2015, MNRAS, 452, 1. The Illustris simulation: the evolving population of black holes across cosmic time https://ui.adsabs.harvard.edu/abs/2015MNRAS.452..575S -.. [Springel2010] Springel (2010), MNRAS, 401, 2. +.. [Siwek2023] Siwek et al. 2023. MNRAS, 522, 2. + Binary orbital evolution in circumbinary discs + https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.2707S/abstract + +.. [Springel2010] Springel 2010, MNRAS, 401, 2. E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh https://ui.adsabs.harvard.edu/abs/2010MNRAS.401..791S -.. [Vogelsberger2014] Vogelsberger et al. (2014), MNRAS, 444, 2. +.. [Vogelsberger2014] Vogelsberger et al. 2014, MNRAS, 444, 2. Introducing the Illustris Project: simulating the coevolution of dark and visible matter in the Universe https://ui.adsabs.harvard.edu/abs/2014MNRAS.444.1518V diff --git a/docs/source/calc_gws.rst b/docs/source/calc_gws.rst index 83bd5664..4bede80e 100644 --- a/docs/source/calc_gws.rst +++ b/docs/source/calc_gws.rst @@ -1,3 +1,5 @@ +.. _gravitational waves: + ================================================ Getting Started: Calculating Gravitational Waves ================================================ diff --git a/docs/source/conf.py b/docs/source/conf.py index 10cc427c..07f5edeb 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -26,28 +26,29 @@ extensions = [ 'sphinx.ext.autodoc', - 'sphinx.ext.autosummary', 'sphinx.ext.autosectionlabel', # {path/to/page}:{title-of-section} 'sphinx.ext.extlinks', 'sphinx.ext.intersphinx', 'sphinx.ext.mathjax', 'sphinx_math_dollar', # enable dollar-signs to write latex math - # 'sphinx.ext.napoleon', + 'sphinx.ext.napoleon', 'sphinx.ext.viewcode', - 'numpydoc', # allow numpy/google style docstrings + # 'numpydoc', # allow numpy/google style docstrings # 'sphinxcontrib.bibtex', # allow bibtex-like citation management ] - -autosummary_generate = True +## commented out the following because we use a python script +# to automatically generate the api_ref/*.rst files +# autosummary_generate = True source_suffix = ['.rst', '.md'] master_doc = 'index' -exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store', 'apidoc_modules'] # NOTE: `numpy` is actually needed, otherwise things break autodoc_mock_imports = [ + "hasasia", "healpy", "sympy", # 'pytest', 'kalepy', 'astropy', 'h5py', 'kalepy', 'matplotlib', 'scipy', 'tqdm', ] @@ -85,65 +86,83 @@ autodoc_member_order = 'bysource' autodoc_default_options = {"members": True, "undoc-members": True, "private-members": True} -numpydoc_show_class_members = False +# numpydoc_show_class_members = False # Report warnings for all validation checks # numpydoc_validation_checks = {"all"} # Report warnings for all checks *except* those listed -numpydoc_validation_checks = { - "all", - # "GL08", # The object does not have a docstring - "ES01", # No extended summary found - "PR01", # Parameters { ... } not documented - "PR02", # Unknown parameters { ... } - "PR07", # Parameter has no description - "PR10", # Parameter "___" requires a space before the colon separating the parameter name and type - # "RT01", # No Returns section found - "RT03", # Return value has no description - # "SS01", # No summary found (a short summary in a single line should be present - # at the beginning of the docstring) - - "SA01", # See Also section not found - "EX01", # No examples section found - - "GL01", # Docstring text (summary) should start in the line immediately after the opening - # quotes (not in the same line, or leaving a blank line in between) - "GL02", # Closing quotes should be placed in the line after the last text in the docstring - # (do not close the quotes in the same line as the text, or leave a blank line - # between the last text and the quotes) - "GL03", # Double line break found; please use only one blank line to separate sections or - # paragraphs, and do not leave blank lines at the end of docstrings - "PR05", # Parameter "___" type should not finish with "." - "PR08", # Parameter "weights" description should start with a capital letter - "PR09", # Parameter "___" description should finish with "." - "RT02", # The first line of the Returns section should contain only the type, - # unless multiple values are being returned - - "RT04", # Return value description should start with a capital letter - "RT05", # Return value description should finish with ". - "SS02", # Summary does not start with a capital letter - "SS03", # Summary does not end with a period - "SS05", # Summary must start with infinitive verb, not third person (e.g. use "Generate" instead of "Generates") -} - - -def run_apidoc(_): - """ - - https://github.com/readthedocs/readthedocs.org/issues/1139#issuecomment-312626491 - - """ - from sphinx.ext.apidoc import main - cur_dir = os.path.abspath(os.path.dirname(__file__)) - sys.path.append(os.path.join(cur_dir, os.path.pardir)) - output_dir = os.path.join(cur_dir, "apidoc_modules") - # docs/source ==> /docs ==> holodeck/ - input_dir = os.path.join(cur_dir, os.path.pardir, os.path.pardir, "holodeck") - main(['-e', '-o', output_dir, input_dir, '--force']) - return - - +# numpydoc_validation_checks = { +# "all", +# # "GL08", # The object does not have a docstring +# "ES01", # No extended summary found +# "PR01", # Parameters { ... } not documented +# "PR02", # Unknown parameters { ... } +# "PR07", # Parameter has no description +# "PR10", # Parameter "___" requires a space before the colon separating the parameter name and type +# # "RT01", # No Returns section found +# "RT03", # Return value has no description +# # "SS01", # No summary found (a short summary in a single line should be present +# # at the beginning of the docstring) + +# "SA01", # See Also section not found +# "EX01", # No examples section found + +# "GL01", # Docstring text (summary) should start in the line immediately after the opening +# # quotes (not in the same line, or leaving a blank line in between) +# "GL02", # Closing quotes should be placed in the line after the last text in the docstring +# # (do not close the quotes in the same line as the text, or leave a blank line +# # between the last text and the quotes) +# "GL03", # Double line break found; please use only one blank line to separate sections or +# # paragraphs, and do not leave blank lines at the end of docstrings +# "PR05", # Parameter "___" type should not finish with "." +# "PR08", # Parameter "weights" description should start with a capital letter +# "PR09", # Parameter "___" description should finish with "." +# "RT02", # The first line of the Returns section should contain only the type, +# # unless multiple values are being returned + +# "RT04", # Return value description should start with a capital letter +# "RT05", # Return value description should finish with ". +# "SS02", # Summary does not start with a capital letter +# "SS03", # Summary does not end with a period +# "SS05", # Summary must start with infinitive verb, not third person (e.g. use "Generate" instead of "Generates") +# } + + +# def run_apidoc(_): +# """ + +# https://github.com/readthedocs/readthedocs.org/issues/1139#issuecomment-312626491 + +# """ +# from sphinx.ext.apidoc import main +# cur_dir = os.path.abspath(os.path.dirname(__file__)) +# sys.path.append(os.path.join(cur_dir, os.path.pardir)) +# output_dir = os.path.join(cur_dir, "apidoc_modules") +# # docs/source ==> /docs ==> holodeck/ +# input_dir = os.path.join(cur_dir, os.path.pardir, os.path.pardir, "holodeck") +# main(['-e', '-o', output_dir, input_dir, '--force']) +# return + + +# def setup(app): +# app.connect('builder-inited', run_apidoc) +# return def setup(app): - app.connect('builder-inited', run_apidoc) - return + # Connect to the skip-member event to stop the duplicate tracking pass + app.connect('autodoc-skip-member', skip_gp_duplicate_methods) + +def skip_gp_duplicate_methods(app, what, name, obj, skip, options): + # Define the specific class methods causing the 5 duplicates + duplicate_methods = { + "GaussProc": ["lnprior", "lnlike", "lnprob"], + "SamModel": ["sam_for_params", "validate_params"] + } + + # Check if the member belongs to our target classes + if hasattr(obj, '__qualname__') and '.' in obj.__qualname__: + classname, _, methodname = obj.__qualname__.rpartition('.') + if classname in duplicate_methods and methodname in duplicate_methods[classname]: + return True # Tell Sphinx to skip documenting this duplicate entry + + return skip diff --git a/docs/source/getting_started.rst b/docs/source/getting_started.rst index 05eee5f1..2059005f 100644 --- a/docs/source/getting_started.rst +++ b/docs/source/getting_started.rst @@ -289,7 +289,7 @@ At times, instead of the finite volume, we wish to calculate the distribution of First we choose an independent variable of interest to define where we sample our population. Because GW signals are our primary interest, it is typically natural to pre-define a grid of observer-frame GW frequencies. Equivalently binary separations can be defined, or redshift intervals, etc, but we typically use observer-frame GW frequencies. Each binary evolution track is interpolated to find the point in time when the binary reaches each frequency of interest. For each frequency, and each binary, we then have the binary intrinsic properties (e.g. $M, q, z$) in additional to its hardening rate ($da/dt$ or $df/dt$). If the corresponding redshift (or scale-factor or age of the Universe) takes place after the current epoch (i.e. $z=0$), then the binary can be discarded at that frequency. -Each binary, at each frequency, is then interpreted as contributing a volume-density of $n_{ij} = 1/V_\mathrm{sim}$ to the total distribution of binaries. We can then convert from volume-density to number of binaries in the same way as for the semi-analytic-model calculation :math:numref:`eq:N_from_n`, +Each binary, at each frequency, is then interpreted as contributing a volume-density of :math:`n_{ij} = 1/V_\mathrm{sim}` to the total distribution of binaries. We can then convert from volume-density to number of binaries in the same way as for the semi-analytic-model calculation :math:numref:`eq:N_from_n`, .. math:: \frac{d N}{d\ln X} = n \frac{dV_c}{dz} \frac{dz}{dt} \frac{d t}{d\ln X}, diff --git a/docs/source/getting_started/index.rst b/docs/source/getting_started/index.rst index 95b4c466..070c5ee8 100644 --- a/docs/source/getting_started/index.rst +++ b/docs/source/getting_started/index.rst @@ -23,7 +23,7 @@ The |holodeck| framework simulates populations of MBH binaries, and calculates t Additional information can be very useful. In particular, information about the host galaxy of the MBH pair can be used in the binary evolution calculation. -(2) :ref:`Binary Evolution`: Evolve the binary population from their initial conditions (i.e. large separations) until they reach the regime of interest (i.e. small separations). In the simplest models, binaries are assumed to coalesce instantaneously, and are assumed to evolve purely due to GW emission. Note that these two assumptions are contradictory. More complex, self-consistent evolution models are recommended. These models typically involve interactions between MBH binaries and their host galaxies ('environmental' interactions). Note that the effects of binary evolution can be broken up into two distinct effects: +(2) :ref:`Binary-Evolution`: Evolve the binary population from their initial conditions (i.e. large separations) until they reach the regime of interest (i.e. small separations). In the simplest models, binaries are assumed to coalesce instantaneously, and are assumed to evolve purely due to GW emission. Note that these two assumptions are contradictory. More complex, self-consistent evolution models are recommended. These models typically involve interactions between MBH binaries and their host galaxies ('environmental' interactions). Note that the effects of binary evolution can be broken up into two distinct effects: (a) The redshift at which binaries reach the given frequencies (or separations) of interest, and similarly which binaries are able to reach those frequencies before redshift zero, and @@ -35,7 +35,7 @@ The |holodeck| framework simulates populations of MBH binaries, and calculates t (a) Discretization: whether binaries are treated as discrete objects, i.e. there can only be an integer number of binaries in a given frequency bin (this often relates to whether the number-density, or total-number of binaries is used in the calculation). One can also consider the effects of cosmic variance in this category as well. - (b) Evolution: whether self-consistent models of binary evolution are considered, or if purely GW-driven evolution is assumed (see :ref:`Binary Evolution`). + (b) Evolution: whether self-consistent models of binary evolution are considered, or if purely GW-driven evolution is assumed (see :ref:`Binary-Evolution`). (c) Eccentricity: whether binaries are restricted to circular orbits, or allowed to have eccentric evolution. Eccentricity has multiple effects on binary evolution, mostly (i) by changing the rate of binary hardening, and (ii) by changing the GW frequencies corresponding to each orbital frequency. Circular binaries emit GWs at only the :math:`n=2` harmonic of the orbital frequency, while eccentric binaries emit at all integer harmonics. diff --git a/docs/source/holodeck.rst b/docs/source/holodeck.rst deleted file mode 100644 index 7de24625..00000000 --- a/docs/source/holodeck.rst +++ /dev/null @@ -1,28 +0,0 @@ -======== -holodeck -======== - -.. automodule:: holodeck - :members: - :undoc-members: - :show-inheritance: - - -Submodules ----------- - -.. toctree:: - :maxdepth: 1 - - apidoc_modules/holodeck.constants - apidoc_modules/holodeck.cosmology - apidoc_modules/holodeck.evolution - apidoc_modules/holodeck.gravwaves - apidoc_modules/holodeck.logger - apidoc_modules/holodeck.observations - apidoc_modules/holodeck.plot - apidoc_modules/holodeck.pop_observational - apidoc_modules/holodeck.population - apidoc_modules/holodeck.relations - apidoc_modules/holodeck.sam - apidoc_modules/holodeck.utils \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst index d32629bd..9e4a90c6 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -41,7 +41,7 @@ Contents of Documentation Calculating Gravitational Waves Definitions and Abbreviations Annotated Bibliography - Full Package Documentation + Full Package Documentation .. toctree:: @@ -63,8 +63,17 @@ Contents of Documentation holodeck.plot holodeck.host_relations holodeck.utils + holodeck.ems +.. toctree:: + :maxdepth: 2 + :caption: Getting Started + + getting_started/index + getting_started/libraries + getting_started/nanograv_15yr + Getting Started =============== diff --git a/holodeck/__init__.py b/holodeck/__init__.py index ec0edf24..f758f495 100644 --- a/holodeck/__init__.py +++ b/holodeck/__init__.py @@ -33,7 +33,7 @@ References ---------- -* [WMAP9] Hinshaw, Larson, Komatsu et al. 2013 +* [WMAP9]_ Hinshaw, Larson, Komatsu et al. 2013 """ diff --git a/holodeck/accretion.py b/holodeck/accretion.py index d1f9dd0f..40acb794 100644 --- a/holodeck/accretion.py +++ b/holodeck/accretion.py @@ -180,7 +180,7 @@ def pref_acc(self, mdot, evol, step, bin=None): See Also -------- - :meth: `Evolution._take_next_step()` : Relationship + :meth:`Evolution._take_next_step()` : Relationship Notes ----- diff --git a/holodeck/anisotropy.py b/holodeck/anisotropy.py index 59cde87f..1adb2102 100644 --- a/holodeck/anisotropy.py +++ b/holodeck/anisotropy.py @@ -253,7 +253,7 @@ def sph_harm_from_hc(hc_ss, hc_bg, nside = NSIDE, lmax = LMAX): def plot_ClC0_medians(fobs, Cl_best, lmax, nshow): xx = fobs*YR - fig, ax = holo.plot.figax(figsize=(8,5), xlabel=holo.plot.LABEL_GW_FREQUENCY_YR, ylabel='$C_{\ell>0}/C_0$') + fig, ax = holo.plot.figax(figsize=(8,5), xlabel=holo.plot.LABEL_GW_FREQUENCY_YR, ylabel=r'$C_{\ell>0}/C_0$') yy = Cl_best[:,:,:,1:]/Cl_best[:,:,:,0,np.newaxis] # (B,F,R,l) yy = np.median(yy, axis=-1) # (B,F,l) median over realizations @@ -689,7 +689,7 @@ def plot_ClC0_versions(fobs_gw_cents, spk=True, bayes=True, analytic=False, Cl_analytic=None, C0_analytic=None, label_analytic='analytic', 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$') + fig, ax = plot.figax(xlabel=plot.LABEL_GW_FREQUENCY_HZ, ylabel=r'$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) diff --git a/holodeck/constants.py b/holodeck/constants.py index 8bc88c0a..6ca03c85 100644 --- a/holodeck/constants.py +++ b/holodeck/constants.py @@ -7,6 +7,7 @@ DEVELOPER NOTE: To preserve fast module import times across holodeck, do NOT import `astropy` or `numpy` in this file to calculate constants at module initialization. If you need to add a new constant: + 1. Open a temporary terminal or notebook. 2. Run your astropy conversion (e.g., `import astropy as ap; print(ap.constants.m_e.cgs.value)`). 3. Copy the raw float literal value and paste it directly below. diff --git a/holodeck/detstats.py b/holodeck/detstats.py index b361a728..348d2503 100644 --- a/holodeck/detstats.py +++ b/holodeck/detstats.py @@ -17,7 +17,7 @@ import holodeck as holo from holodeck import utils, log, plot # , cosmo, anisotropy from holodeck.constants import MPC, YR -from holodeck.sams import cyutils as sam_cyutils +from holodeck.sams import sam_cyutils try: from sympy import nsolve, Symbol @@ -137,7 +137,7 @@ def _orf_ij(i, j, theta_ij): def _orf_pta(pulsars): """ Calculate the overlap reduction function matrix Gamma for a list of hasasia.Pulsar objects - Paramters + Parameters --------- pulsars : (P,) list of hasasia.Pulsar objects. @@ -274,7 +274,7 @@ def _sigma0_Bstatistic(noise, Gamma, Sh0_bg): return sigma_0B def _sigma1_Bstatistic(noise, Gamma, Sh_bg, Sh0_bg): - """ Calculate sigma_1 for the background, by summing over all pulsars and frequencies. + r"""Calculate sigma_1 for the background, by summing over all pulsars and frequencies. Assuming the B statistic, which maximizes S/N_B = mu_1/sigma_1 Parameters @@ -287,7 +287,7 @@ def _sigma1_Bstatistic(noise, Gamma, Sh_bg, Sh0_bg): Spectral density in the background. Sh0_bg : (F,R) 1Darray of scalars Value of spectral density used to construct the statistic. -s + Returns ------- sigma_1B : (R,) 1Darray @@ -759,7 +759,7 @@ def _n_unitary_vector(theta, phi, xi): """ Calculate the unitary vector n-hat for the antenna pattern functions for each of S sky realizations. - Paramters + Parameters --------- theta : (F,S,L) NDarray Spherical coordinate position of each single source. @@ -788,7 +788,7 @@ def _Omega_unitary_vector(theta, phi): """ Calculate the unitary vector n-hat for the antenna pattern functions for each of S sky realizations. - Paramters + Parameters --------- theta : (F,S,L) NDarray Spherical coordinate position of each single source. @@ -835,23 +835,22 @@ def _pi_unitary_vector(phi_i, theta_i): ###################### Antenna Pattern Functions ###################### def dotprod(vec1, vec2): - """ Calculate the dot product for NDarrays of 3D vectors, with - vector elements specified by the first index. + r"""Calculate the dot product for NDarrays of 3D vectors, with vector elements specified by the first index. - Parameters - ---------- - vec1 : (3,N1,N2,N3,...N) NDarray - vec2 : (3,N1,N2,N3,...N) NDarray + Parameters + ---------- + vec1 : (3,N1,N2,N3,...N) NDarray + vec2 : (3,N1,N2,N3,...N) NDarray - Returns - ------- - dotted : (N1,N2,N3,...N) NDarray - The dot product of the vectors specified by the first dimension, - for every N1, N2, N3,...N. + Returns + ------- + dotted : (N1,N2,N3,...N) NDarray + The dot product of the vectors specified by the first dimension, + for every N1, N2, N3,...N. Example: find the dot product of 3D vectors for every P,F,R, using NDarrays of shape (3,P,F,R) - """ + """ dotted = vec1[0,...]*vec2[0,...] + vec1[1,...]*vec2[1,...] + vec1[2,...]*vec2[2,...] return dotted @@ -1140,11 +1139,11 @@ def _amplitude(hc_ss, fobs, dfobs): ####################### SS Signal to Noise Ratio ####################### def _snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): - """ Calculate the SNR for each pulsar wrt each single source detection, + r"""Calculate the SNR for each pulsar wrt each single source detection, for S sky realizations and R strain realizations. - Paramters - --------- + Parameters + ---------- amp : (F,R,L) NDarray Dimensionless strain amplitude for loudest source at each frequency. F_iplus : (P,F,S,L) NDarray @@ -1163,6 +1162,7 @@ def _snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): freqs : (F,) 1Darray Observed frequency bin centers. + Returns ------- snr_ss : (F,R,S,L) NDarray @@ -1176,11 +1176,11 @@ def _snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): return snr_ss def _snr_ss_5dim(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): - """ Calculate the SNR for each pulsar wrt each single source detection, + r"""Calculate the SNR for each pulsar wrt each single source detection, for S sky realizations and R strain realizations. - Paramters - --------- + Parameters + ---------- amp : (F,R,L) NDarray Dimensionless strain amplitude for loudest source at each frequency. F_iplus : (P,F,S,L) NDarray @@ -2216,7 +2216,7 @@ def plot_sample_nn(fobs, hc_ss, hc_bg, dp_ss, dp_bg, df_ss, df_bg, nn): def plot_detprob(dp_ss_all, dp_bg_all, nsamps): """ Plot detection probability for many samples. - Paramaters + Parameters ---------- dp_ss_all : (N,R, S) NDarray Single source detection probably of each strain and sky realization of each sample. @@ -2230,7 +2230,7 @@ def plot_detprob(dp_ss_all, dp_bg_all, nsamps): """ fig, ax = plt.subplots(figsize=(6.5,4)) ax.set_xlabel('Param Space Sample') - ax.set_ylabel('Detection Probability, $\gamma$') + ax.set_ylabel(r'Detection Probability, $\gamma$') ax.errorbar(np.arange(nsamps), np.mean(dp_bg_all, axis=1), yerr = np.std(dp_bg_all, axis=1), linestyle='', marker='d', capsize=5, color='cornflowerblue', alpha=0.5, @@ -2251,7 +2251,7 @@ def plot_detprob(dp_ss_all, dp_bg_all, nsamps): def plot_detfrac(df_ss, df_bg, nsamps, thresh): """ Plot detection fraction for many samples. - Paramaters + Parameters ---------- df_all : (N,) NDarray Fraction of strain and sky realization with a 'dp_ss'>'thresh'. diff --git a/holodeck/discrete/evolution.py b/holodeck/discrete/evolution.py index aa950525..33323320 100644 --- a/holodeck/discrete/evolution.py +++ b/holodeck/discrete/evolution.py @@ -279,6 +279,7 @@ def at(self, xpar, targets, params=None, coal=False, lin_interp=None): argument. The behavior of this function is broken into three sub-functions, that are only used here: + * :meth:`Evolution._at__inputs` : parse the input arguments. * :meth:`Evolution._at__index_frac` : find the indices in the evolutionary tracks bounding the target interpolation locations, and also the fractional distance to interpolate @@ -292,6 +293,7 @@ def at(self, xpar, targets, params=None, coal=False, lin_interp=None): String specifying the variable of interpolation. targets : array_like, Locations to interpolate to: + * if ``xpar == sepa`` : binary separation, units of [cm], * if ``xpar == fobs`` : binary orbital freq, observer-frame, units of [1/sec], params : None or (list of str) @@ -302,10 +304,11 @@ def at(self, xpar, targets, params=None, coal=False, lin_interp=None): Interpolated values for other binaries are set to `np.nan`. lin_interp : None or bool, Interpolate parameters in linear space. + * True : all parameters interpolated in lin-lin space. * False: all parameters interpolated in log-log space. * None : parameters are interpolated in log-log space, unless they're included in the - :attr:`Evolution._LIN_INTERP_PARS` attribute. + :attr:`Evolution._LIN_INTERP_PARS` attribute. Returns ------- diff --git a/holodeck/gps/gp_utils.py b/holodeck/gps/gp_utils.py index bf776d44..5adcbe92 100644 --- a/holodeck/gps/gp_utils.py +++ b/holodeck/gps/gp_utils.py @@ -172,7 +172,7 @@ def train_gp(spectra_file, y_is_variance=False, kernel="ExpSquaredKernel", mpi=True): - """Train gaussian processes on the first `nfreqs` of the GWB in `spectra_file`. + """Train gaussian processes on the first `nfreqs` of the GWB in ``spectra_file``. Parameters ---------- @@ -195,7 +195,7 @@ def train_gp(spectra_file, kernel : str, optional The type of kernel to use for the GP kernel_opts : dict, optional - The options to pass when constructing the kernel. Unpacks as **kwargs to george.kernels + The options to pass when constructing the kernel. Unpacks as ``**kwargs`` to george.kernels mpi : bool, optional Whether to use MPI or Python's multiprocessing module diff --git a/holodeck/hardening.py b/holodeck/hardening.py index fbf4ed19..1d25952f 100644 --- a/holodeck/hardening.py +++ b/holodeck/hardening.py @@ -2,19 +2,29 @@ To-Do (Hardening) ----------------- -*Dynamical_Friction_NFW + +* Dynamical_Friction_NFW + * Allow stellar-density profiles to also be specified (instead of using a hard-coded Dehnen profile) * Generalize calculation of stellar characteristic radius. Make self-consistent with stellar-profile, and user-specifiable. + * Evolution + * `_sample_universe()` : sample in comoving-volume instead of redshift + * Sesana_Scattering + * Allow stellar-density profile (or otherwise the binding-radius) to be user-specified and flexible. Currently hard-coded to Dehnen profile estimate. + * _SHM06 + * Interpolants of hardening parameters return 1D arrays. + * Fixed_Time_2PL + * Handle `rchar` better with respect to interpolation. Currently not an interpolation variable, which restricts it's usage. * This class should be separated into a generic `_Fixed_Time` class that can use any @@ -235,7 +245,7 @@ def consistent(self): class CBD_Torques(_Hardening): - """Binary Orbital Evolution based on Hydrodynamic Simulations by Siwek+23. + """Orbital evolution of binaries in circumbinary discs by Siwek+23. https://academic.oup.com/mnras/article/522/2/2707/7131464 This module uses data from Siwek+23, which supplies rates of change of @@ -787,21 +797,22 @@ def _attenuation_BBR1980(self, sepa, m1, m2, mstar): class Fixed_Time_2PL(_Hardening): - """Provide a binary hardening rate such that the total lifetime matches a given value. + r"""Provide a binary hardening rate such that the total lifetime matches a given value. This class uses a phenomenological functional form (defined in :meth:`Fixed_Time.function`) to model the hardening rate ($da/dt$) of binaries. The functional form is, .. math:: - \\dot{a} = - A * (1.0 + x)^{-g_2 - 1} / x^{g_1 - 1}, + \dot{a} = - A * (1.0 + x)^{-g_2 - 1} / x^{g_1 - 1}, - where :math:`x \\equiv a / r_\\mathrm{char}` is the binary separation scaled to a characteristic - transition radius (:math:`r_\\mathrm{char}`) between two power-law indices $g_1$ and $g_2$. There is + where :math:`x \equiv a / r_\mathrm{char}` is the binary separation scaled to a characteristic + transition radius (:math:`r_\mathrm{char}`) between two power-law indices $g_1$ and $g_2$. There is also an overall normalization $A$ that is calculated to yield the desired binary lifetimes. NOTE/BUG: the actual binary lifetimes tend to be 1-5% shorter than the requested value. The normalization for each binary, to produce the desired lifetime, is calculated as follows: + (1) A set of random test binary parameters are chosen. (2) The normalization constants are determined, using least-squares optimization, to yield the desired lifetime. @@ -811,6 +822,7 @@ class Fixed_Time_2PL(_Hardening): Construction/Initialization: note that in addition to the standard :meth:`Fixed_Time.__init__` constructor, there are two additional constructors are provided: + * :meth:`Fixed_Time.from_pop` - accept a :class:`holodeck.population._Discrete_Population`, * :meth:`Fixed_Time.from_sam` - accept a :class:`holodeck.sam.Semi_Analytic_Model`. @@ -828,7 +840,7 @@ class Fixed_Time_2PL(_Hardening): def __init__(self, time, mtot, mrat, redz, sepa_init, rchar=100.0*PC, gamma_inner=-1.0, gamma_outer=+1.5, progress=False, interpolate_norm=False): - """Initialize `Fixed_Time` instance for the given binary properties and function parameters. + r"""Initialize `Fixed_Time` instance for the given binary properties and function parameters. Parameters ---------- @@ -841,7 +853,7 @@ def __init__(self, time, mtot, mrat, redz, sepa_init, mtot : array_like Binary total-mass [gram]. mrat : array_like - Binary mass-ratio $q \\equiv m_2 / m_1 \\leq 1$. + Binary mass-ratio $q \equiv m_2 / m_1 \leq 1$. redz : array_like Binary Redshift. NOTE: this is only used as an argument to callable `rchar` and `time` values. @@ -950,10 +962,12 @@ def from_pop(cls, pop, time, **kwargs): Input population, from which to use masses, redshifts and separations. time : float, callable or array_like Total merger time of binaries, units of [sec], specifiable in the following ways: + * float : uniform merger time for all binaries * callable : function `time(mtot, mrat, redz)` which returns the total merger time * array_like : (N,) matching the shape of `mtot` (etc) giving the merger time for each binary + **kwargs : dict Additional keyword-argument pairs passed to the `Fixed_Time` initialization method. @@ -975,15 +989,19 @@ def from_sam(cls, sam, time, sepa_init=1e4*PC, **kwargs): Input population, from which to use masses, redshifts and separations. time : float, callable or array_like Total merger time of binaries, units of [sec], specifiable in the following ways: + * float : uniform merger time for all binaries * callable : function `time(mtot, mrat, redz)` which returns the total merger time * array_like : (N,) matching the shape of `mtot` (etc) giving the merger time for each binary + sepa_init : float or array_like Initial binary separation. Units of [cm]. + * float : initial separation applied to all binaries, * array_like : initial separations for all binaries, shaped (N,) matching the number binaries. + **kwargs : dict Additional keyword-argument pairs passed to the `Fixed_Time` initialization method. @@ -1035,7 +1053,7 @@ def dadt_dedt(self, evo, step, bin=None): @classmethod def _dadt_dedt(cls, mtot, mrat, sepa, norm, rchar, gamma_inner, gamma_outer): - """Calculate hardening rate for the given raw parameters. + r"""Calculate hardening rate for the given raw parameters. Parameters ---------- @@ -1078,17 +1096,17 @@ def _dadt_dedt(cls, mtot, mrat, sepa, norm, rchar, gamma_inner, gamma_outer): @classmethod def function(cls, norm, xx, gamma_inner, gamma_outer): - """Hardening rate given the parameters for this hardening model. + r"""Hardening rate given the parameters for this hardening model. The functional form is, .. math:: - \dot{a} = - A * (1.0 + x)^{-g_out + g_in} / x^{g_in - 1}, + \dot{a} = - A * (1.0 + x)^{-g_{\mathrm{out}} + g_{\mathrm{in}}} / x^{g_{\mathrm{in}} - 1}, - Where $A$ is an overall normalization, and x \\equiv a / r_\\mathrm{char}$ is the binary - separation scaled to a characteristic transition radius ($r_\\mathrm{char}$) between two - power-law indices $g_inner$ and $g_outer$. + Where $A$ is an overall normalization, and $x \equiv a / r_{\mathrm{char}}$ is the binary + separation scaled to a characteristic transition radius ($r_{\mathrm{char}}$) between two + power-law indices $g_{\mathrm{inner}}$ and $g_{\mathrm{outer}}$. Parameters ---------- @@ -1110,9 +1128,10 @@ def function(cls, norm, xx, gamma_inner, gamma_outer): @classmethod def _calculate_norm_interpolant(cls, rchar, gamma_inner, gamma_outer): - """Generate interpolants to map from binary parameters to hardening rate normalization. + r"""Generate interpolants to map from binary parameters to hardening rate normalization. Interpolants are calculated as follows: + (1) A set of random test binary parameters and lifetimes are chosen. (2) The normalizations to yield those binary lifetimes are calculated with least-squares optimization. @@ -1127,8 +1146,10 @@ def _calculate_norm_interpolant(cls, rchar, gamma_inner, gamma_outer): ---------- rchar : scalar or array_like #! Possible that only a scalar value is currently working! Characteristic radius separating the two power-law regimes, in units of [cm]: + * scalar : uniform radius for all binaries * array_like : characteristic radius for each binary. + gamma_inner : scalar Power-law of hardening timescale in the stellar-scattering regime, (small separations: r < rchar), at times referred to internally as `g1`. @@ -1247,7 +1268,7 @@ def convert_points_to_interp_vals(points): @classmethod def _get_norm_chunk(cls, target_time, *args, progress=True, **kwargs): - """Calculate normalizations in 'chunks' of the input arrays, to obtain the target lifetime. + r"""Calculate normalizations in 'chunks' of the input arrays, to obtain the target lifetime. Calculates normalizations for groups of parameters of size `chunk` at a time. Loops over these chunks until all inputs have been processed. Calls :meth:`Fixed_Time._get_norm` to @@ -1304,7 +1325,7 @@ def _get_norm_chunk(cls, target_time, *args, progress=True, **kwargs): @classmethod def _get_norm(cls, target_time, *args, guess=1e7, max_err=1e-6): - """Calculate normalizations of the input arrays, to obtain the target binary lifetime. + r"""Calculate normalizations of the input arrays, to obtain the target binary lifetime. Uses deterministic least-squares optimization to find the best normalization values, using `scipy.optimize.newton`. @@ -1358,7 +1379,7 @@ def time_total(self, mt, mr): @classmethod def _time_total(cls, norm, mt, mr, rchar, gamma_inner, gamma_outer, sepa_init, num=123): - """For the given parameters, integrate the binary evolution to find total lifetime. + r"""For the given parameters, integrate the binary evolution to find total lifetime. Parameters ---------- @@ -1427,7 +1448,7 @@ def _time_total(cls, norm, mt, mr, rchar, gamma_inner, gamma_outer, sepa_init, n class Fixed_Time(Fixed_Time_2PL): - """Legacy Version of `Fixed_Time` class such that outer power-law in time (a/dadt) is g1+g2+1. + r"""Legacy Version of `Fixed_Time` class such that outer power-law in time ($a/(da/dt)^{-1}$) is $g_1+g_2+1$. """ def __init__(self, *args, **kwargs): @@ -1450,16 +1471,16 @@ def __init__(self, *args, **kwargs): @classmethod def function(cls, norm, xx, g1, g2): - """Hardening rate given the parameters for this hardening model. + r"""Hardening rate given the parameters for this hardening model. The functional form is, .. math:: - \\dot{a} = - A * (1.0 + x)^{-g_2 - 1} / x^{g_1 - 1}, + \dot{a} = - A * (1.0 + x)^{-g_2 - 1} / x^{g_1 - 1}, - Where $A$ is an overall normalization, and x \\equiv a / r_\\mathrm{char}$ is the binary - separation scaled to a characteristic transition radius ($r_\\mathrm{char}$) between two + Where $A$ is an overall normalization, and $x \equiv a / r_{\mathrm{char}}$ is the binary + separation scaled to a characteristic transition radius ($r_\mathrm{char}$) between two power-law indices $g_1$ and $g_2$. Parameters @@ -1468,7 +1489,7 @@ def function(cls, norm, xx, g1, g2): Hardening rate normalization, units of [cm/s]. xx : array_like Dimensionless binary separation, the semi-major axis in units of the characteristic - (i.e. transition) radius of the model `rchar`. + (i.e. transition) radius of the model ``rchar``. g1 : scalar Power-law of hardening timescale in the stellar-scattering regime, (small separations: r < rchar). @@ -1796,7 +1817,7 @@ def _init_h(self): class _Siwek2023: - """ Hardening rates from circumbinary disk simulations as in [Siwek2023]_. + r""" Hardening rates from circumbinary disk simulations as in [Siwek2023]_. Mass ratios and eccentricities must be provided. diff --git a/holodeck/logger.py b/holodeck/logger.py index 17e1a8d6..0b521937 100644 --- a/holodeck/logger.py +++ b/holodeck/logger.py @@ -147,11 +147,11 @@ def _set_level(self, lvl): def log_to_file(logger, file_level=DEBUG, file_name=None, base_name='holodeck', path=None): """Add a `FileHandler` to the given logger, to log to an output (text) file. - If `file_name` IS given, then it is used as the output filename. - If `file_name` is NOT an absolute path, then the file is placed in the `path` directory. - If `file_name` IS an absolute path, then that's where the file is placed. - If `file_name` is NOT given, then a filename is constructed based on the `base_name`, the logging level, and - the processor rank (if this is a parallel job). + * If `file_name` IS given, then it is used as the output filename. + * If `file_name` is NOT an absolute path, then the file is placed in the `path` directory. + * If `file_name` IS an absolute path, then that's where the file is placed. + * If `file_name` is NOT given, then a filename is constructed based on the `base_name`, the logging level, and + the processor rank (if this is a parallel job). """ diff --git a/holodeck/param_spaces.py b/holodeck/param_spaces.py index 8dc4f30e..f056c730 100644 --- a/holodeck/param_spaces.py +++ b/holodeck/param_spaces.py @@ -3,7 +3,7 @@ import holodeck as holo from holodeck.constants import PC, GYR, MSOL -from holodeck.librarian import ( +from holodeck.librarian.lib_tools import ( _Param_Space, PD_Uniform, PD_Piecewise_Uniform_Density, PD_Normal, # PD_Uniform_Log, diff --git a/holodeck/pop_observational.py b/holodeck/pop_observational.py index b250540a..c685e43b 100644 --- a/holodeck/pop_observational.py +++ b/holodeck/pop_observational.py @@ -26,6 +26,7 @@ import holodeck as holo from holodeck import _PATH_DATA, cosmo, log from holodeck.constants import MSOL +from holodeck.discrete.population import _Population_Discrete # physical constants for natural units c = G = 1 c = 2.99792458*(10**8) @@ -41,7 +42,7 @@ _DEF_OBSERVATIONAL_FNAME = "observational_2mass_galaxy-catalog_extended.npz" -class BP_Observational(holo.population._Population_Discrete): +class BP_Observational(_Population_Discrete): FREQ_MIN = 1e-9 # Hz, minimum of PTA band of interest diff --git a/holodeck/sams/components.py b/holodeck/sams/components.py index d0ec276e..40f9e9d5 100644 --- a/holodeck/sams/components.py +++ b/holodeck/sams/components.py @@ -4,7 +4,7 @@ References ---------- * [Sesana2008]_ Sesana, Vecchio, Colacino 2008. -* [Rodriguez-Gomez+2015]_ Rodriguez-Gomez, Genel, Vogelsberger, et al. 2015 +* [Rodriguez-Gomez2015]_ Rodriguez-Gomez, Genel, Vogelsberger, et al. 2015 The merger rate of galaxies in the Illustris simulation: a comparison with observations and semi-empirical models https://ui.adsabs.harvard.edu/abs/2015MNRAS.449...49R/abstract * [Chen2019]_ Chen, Sesana, Conselice 2019. @@ -421,7 +421,7 @@ def __call__(self, mass, mrat, redz): class GMR_Power_Law(_Galaxy_Merger_Rate): """Galaxy Merger Rate - based on multiple power-laws. - See [Rodriguez-Gomez+2015], Table 1. + See [Rodriguez-Gomez2015], Table 1. "merger rate as a function of descendant stellar mass M_star, progenitor stellar mass ratio mu_star" """ @@ -527,7 +527,7 @@ def __call__(self, mtot, mrat, redz): class GMR_Illustris(GMR_Power_Law): """Galaxy Merger Rate - based on fits to Illustris cosmological simulations. - See [Rodriguez-Gomez+2015], Table 1. + See [Rodriguez-Gomez2015], Table 1. "merger rate as a function of descendant stellar mass M_star, progenitor stellar mass ratio mu_star" """ diff --git a/holodeck/sams/sam.py b/holodeck/sams/sam.py index a96d65f5..3d58c92c 100644 --- a/holodeck/sams/sam.py +++ b/holodeck/sams/sam.py @@ -1331,6 +1331,7 @@ def add_scatter_to_masses(mtot, mrat, dens, scatter, refine=4, log=None): """Add the given scatter to masses m1 and m2, for the given distribution of binaries. The procedure is as follows (see `dev-notebooks/sam-ndens-scatter.ipynb`): + * (1) The density is first interpolated to a uniform, regular grid in (m1, m2) space. A 2nd-order interpolant is used first. A 0th-order interpolant is used to fill-in bad values. In-between, a 1st-order interpolant is used if `linear_interp_backup` is True. @@ -1340,6 +1341,7 @@ def add_scatter_to_masses(mtot, mrat, dens, scatter, refine=4, log=None): Parameters ---------- + mtot : (M,) ndarray Total masses in grams. mrat : (Q,) ndarray diff --git a/holodeck/sams/sam_cyutils.pyx b/holodeck/sams/sam_cyutils.pyx index 9eb02715..f3b4067d 100644 --- a/holodeck/sams/sam_cyutils.pyx +++ b/holodeck/sams/sam_cyutils.pyx @@ -127,9 +127,11 @@ def integrate_differential_number_3dx1d(edges, dnum): Arguments --------- edges : (4,) array_like w/ lengths M, Q, Z, F+1 - Grid edges of `mtot`, `mrat`, `redz`, and `freq` - NOTE: `mtot` should be passed as regular `mtot`, NOT log10(mtot) - `freq` should be passed as regular `freq`, NOT ln(freq) + Grid edges of ``mtot``, ``mrat``, ``redz``, and ``freq``. + + * NOTE: ``mtot`` should be passed as regular ``mtot``, NOT ``log10(mtot)``. + * NOTE: ``freq`` should be passed as regular ``freq``, NOT ``ln(freq)``. + dnum : (M, Q, Z, F) Differential number of binaries, dN/[dlog10M dq qz dlnf] where 'N' is in units of dimensionless number. @@ -838,8 +840,8 @@ def snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): Antenna pattern function for each pulsar. iotas : (F,S,L) NDarray Inclination, used to calculate: - a_pol = 1 + np.cos(iotas) **2 - b_pol = -2 * np.cos(iotas) + ``a_pol = 1 + np.cos(iotas) **2`` + ``b_pol = -2 * np.cos(iotas)`` dur : scalar Duration of observations. Phi_0 : (F,S,L) NDarray @@ -894,8 +896,8 @@ cdef int _snr_ss( Antenna pattern function for each pulsar. iotas : (F,S,L) NDarray Inclination, used to calculate: - a_pol = 1 + np.cos(iotas) **2 - b_pol = -2 * np.cos(iotas) + ``a_pol = 1 + np.cos(iotas) **2`` + ``b_pol = -2 * np.cos(iotas)`` dur : scalar Duration of observations. Phi_0 : (F,S,L) NDarray diff --git a/holodeck/single_sources.py b/holodeck/single_sources.py index 08bef42c..88eac17f 100644 --- a/holodeck/single_sources.py +++ b/holodeck/single_sources.py @@ -29,8 +29,8 @@ log.setLevel(logging.INFO) par_names = np.array(['mtot', 'mrat', 'redz_init', 'redz_final', 'dcom_final', 'sepa_final', 'angs_final']) -par_labels = np.array(['Total Mass $M$ ($M_\odot$)', 'Mass Ratio $q$', 'Initial Redshift $z_i$', 'Final Redshift $z_f$', - 'Final Comoving Distance $d_c$ (Mpc)', 'Final Separation (pc)', 'Final Angular Separation (rad)']) +par_labels = np.array([r'Total Mass $M$ ($M_\odot$)', r'Mass Ratio $q$', r'Initial Redshift $z_i$', r'Final Redshift $z_f$', + r'Final Comoving Distance $d_c$ (Mpc)', r'Final Separation (pc)', r'Final Angular Separation (rad)']) par_units = np.array([1/MSOL, 1, 1, 1, 1/MPC, 1/PC, 1]) ################################################### @@ -39,7 +39,7 @@ def ss_gws_redz(edges, redz, number, realize, loudest = 1, params = False): - """ Calculate strain from the loudest single sources and background. + r"""Calculate strain from the loudest single sources and background. Parameters @@ -56,8 +56,10 @@ def ss_gws_redz(edges, redz, number, realize, loudest = 1, params = False): `dnum` over each bin. realize : int Specification of how many discrete realizations to construct. - loudest : int + loudest : int, optional (default: 1) Number of loudest single sources to separate from background. + params : bool, optional (default: False) + Whether to return the astrophysical parameter arrays. Returns @@ -67,12 +69,13 @@ def ss_gws_redz(edges, redz, number, realize, loudest = 1, params = False): hc_bg : (F, R) NDarray of scalars Characteristic strain of the GWB. sspar : (4, F, R, L) NDarray of scalars - Astrophysical parametes (total mass, mass ratio, initial redshift, final redshift) of each + Astrophysical parameters (total mass, mass ratio, initial redshift, final redshift) of each loud single sources, for each frequency and realization. Returned only if params = True. bgpar : (7, F, R) NDarray of scalars Average effective binary astrophysical parameters (total mass, mass ratio, initial redshift, - final redshift, final comoving distances, final separation, final angular separation) for background sources at each frequency and realization, + final redshift, final comoving distances, final separation, final angular separation) for + background sources at each frequency and realization, Returned only if params = True. """ @@ -178,7 +181,7 @@ def ss_gws_redz(edges, redz, number, realize, loudest = 1, params = False): # version for running libraries def ss_gws(edges, number, realize, loudest = 1, params = False): - """ Calculate strain from the loudest single sources and background. + r"""Calculate strain from the loudest single sources and background. Parameters @@ -190,11 +193,13 @@ def ss_gws(edges, number, realize, loudest = 1, params = False): number : (M, Q, Z, F) ndarray of scalars The number of binaries in each bin of parameter space. This is calculated by integrating `dnum` over each bin. - realize : int, + realize : int Specification of how many discrete realizations to construct. - loudest : int + loudest : int, optional (default: 1) Number of loudest single sources to separate from background. - + params : bool, optional (default: False) + Whether to return the astrophysical parameter arrays. + Returns ------- @@ -282,8 +287,7 @@ def ss_gws(edges, number, realize, loudest = 1, params = False): def loudest_by_cython(edges, number, realize, loudest, round = True, params = False): - """ More efficient way to calculate strain from numbered - grid integrated + r"""More efficient way to calculate strain from numbered grid integrated Parameters @@ -295,16 +299,17 @@ def loudest_by_cython(edges, number, realize, loudest, round = True, params = Fa number : (M, Q, Z, F) ndarray of scalars The number of binaries in each bin of parameter space. This is calculated by integrating `dnum` over each bin. - realize : int, + realize : int Specification of how many discrete realizations to construct. TODO: Set up option for `bool` value, to get multiple sources without realizing. That makes no sense though. loudest : int Number of loudest single sources to separate from background. - round : bool + round : bool, optional (default: True) Specification of whether to discretize the sample if realize is False, by rounding number of binaries in each bin to integers. Does nothing if realize is True. - + params : bool, optional (default: False) + Whether to return the astrophysical parameter arrays. Returns ------- @@ -388,8 +393,7 @@ def loudest_by_cython(edges, number, realize, loudest, round = True, params = Fa def ss_by_cdefs(edges, number, realize, round = True, params = False): - """ More efficient way to calculate strain from numbered - grid integrated + r"""More efficient way to calculate strain from numbered grid integrated Parameters @@ -401,15 +405,16 @@ def ss_by_cdefs(edges, number, realize, round = True, params = False): number : (M, Q, Z, F) ndarray of scalars The number of binaries in each bin of parameter space. This is calculated by integrating `dnum` over each bin. - realize : bool or int, + realize : bool or int Specification of how to construct one or more discrete realizations. If a `bool` value, then whether or not to construct a realization. If a `int` value, then how many discrete realizations to construct. - round : bool + round : bool, optional (default: True) Specification of whether to discretize the sample if realize is False, by rounding number of binaries in each bin to integers. Does nothing if realize is True. - + params : bool, optional (default: False) + Whether to return the astrophysical parameter arrays. Returns ------- @@ -426,11 +431,11 @@ def ss_by_cdefs(edges, number, realize, round = True, params = False): hsamp : (M, Q, Z, F) NDarray Strain amplitude of a single source in every bin (regardless of if that bin actually has any sources.) - bgpar : (3, F, R) NDarray of scalars + bgpar : (3, F, R) NDarray of scalars, optional Average effective binary astrophysical parameters for background sources at each frequency and realization, returned only if params = True. - sspar : (3, F, R) NDarray of scalars + sspar : (3, F, R) NDarray of scalars, optional Astrophysical parametes of single sources at each frequency for each realizations, returned only if params = True. @@ -520,8 +525,7 @@ def ss_by_cdefs(edges, number, realize, round = True, params = False): def ss_by_ndars(edges, number, realize, round = True): - """ More efficient way to calculate strain from numbered - grid integrated + r"""More efficient way to calculate strain from numbered grid integrated Parameters @@ -537,17 +541,10 @@ def ss_by_ndars(edges, number, realize, round = True): Specification of how to construct one or more discrete realizations. If a `bool` value, then whether or not to construct a realization. If a `int` value, then how many discrete realizations to construct. - round : bool + round : bool, optional (default: True) Specification of whether to discretize the sample if realize is False, by rounding number of binaries in each bin to integers. Does nothing if realize is True. - ss : bool - Whether or not to separate the loudest single source in each frequency bin. - sum : bool - Whether or not to sum the strain at a given frequency over all bins. - print_test : bool - Whether or not to print variable as they are calculated, for dev purposes. - Returns ------- @@ -569,7 +566,8 @@ def ss_by_ndars(edges, number, realize, round = True): The number of binaries in each bin after the loudest single source at each frequency is subtracted out. - + Notes + ----- In the unlikely scenario that there are two equal hsmaxes (at same OR dif frequencies), ssidx calculation will go wrong Could avoid this by using argwhere for each f_idx column separately. @@ -736,8 +734,7 @@ def ss_by_ndars(edges, number, realize, round = True): def h2fdf(edges): - """ More efficient way to calculate strain from numbered - grid integrated + r"""More efficient way to calculate strain from numbered grid integrated Parameters @@ -750,9 +747,9 @@ def h2fdf(edges): Returns ------- - h2fdf : (4,) ndarray of scalars + h2fdf : (M, Q, Z, F) ndarray of scalars strain amplitude squared x frequency x frequency bin width for a single - source in each bing + source in each bin """ @@ -805,20 +802,20 @@ def h2fdf(edges): def all_sspars(fobs_gw_cents, sspar): - """ Calculate all single source parameters incl. - ['mtot' 'mrat' 'redz_init' 'redz_final' 'dcom_final' 'sepa_final' 'angs_final'] + r"""Calculate all single source parameters. + Includes 'mtot', 'mrat', 'redz_init', 'redz_final', 'dcom_final', 'sepa_final', 'angs_final'. Parameters ---------- fobs_gw_cents : (F,) 1Darray of scalars Observed gw frequency bin centers. - sspar : (4, F,R,L) NDarray + sspar : (4, F, R, L) NDarray Single source parameters as calculated by ss_gws_redz(). Includes mtot, mrat, redz_init, redz_final. Returns ------- - sspar_all : (7,F,R,L) NDarray + sspar_all : (7, F, R, L) NDarray All single source parameters, corresponding to those in bgpar as calculated by ss_gws_redz(). Includes mtot, mrat, redz_init, redz_final, dcom_final, sepa_final (cm), and angs_final. """ @@ -842,22 +839,26 @@ def parameters_from_indices(edges, ssidx): frequency given their indices and edges. Parameters - ----------- + ---------- edges : (4,) list of 1darrays A list containing the edges along each dimension. The four dimensions correspond to total mass, mass ratio, redshift, and observer-frame orbital frequency. The length of each of the four arrays is M+1, Q+1, Z+1, F+1. - ssidx : (3, F, R) or (3,F) ndarray + ssidx : (3, F, R) or (3, F) ndarray The indices of loudest single sources at each frequency of each realization in the format: [[M indices], [q indices], [z indices], [f indices], [r indices]] Returns - ---------- - m_arr : (F,) or (F,R) Ndarray of scalars - q_arr : (F,) or (F,R) Ndarray of scalars - z_arr : (F,) or (F,R) Ndarray of scalars - f_arr : (F,) or (F,R) Ndarray of scalars + ------- + m_arr : (F,) or (F, R) Ndarray of scalars + Total masses. + q_arr : (F,) or (F, R) Ndarray of scalars + Mass ratios. + z_arr : (F,) or (F, R) Ndarray of scalars + Redshifts. + f_arr : (F,) or (F, R) Ndarray of scalars + Frequencies. """ # Frequency bin midpoints @@ -879,8 +880,7 @@ def parameters_from_indices(edges, ssidx): def ss_by_loops(edges, number, realize=False, round=True, print_test = False): - """ Inefficient way to calculate strain from numbered - grid integrated + r"""Inefficient way to calculate strain from numbered grid integrated Parameters ---------- @@ -891,14 +891,14 @@ def ss_by_loops(edges, number, realize=False, round=True, print_test = False): number : (M, Q, Z, F) ndarray The number of binaries in each bin of parameter space. This is calculated by integrating `dnum` over each bin. - realize : bool or int, + realize : bool or int, optional (default=False) Specification of how to construct one or more discrete realizations. If a `bool` value, then whether or not to construct a realization. If a `int` value, then how many discrete realizations to construct. - round : bool + round : bool, optional (default=True) Specification of whether to discretize the sample if realize is False, by rounding number of binaries in each bin to integers. - print_test : bool + print_test : bool, optional (default=False) Whether or not to print variable as they are calculated, for dev purposes. @@ -1076,8 +1076,7 @@ def ss_by_loops(edges, number, realize=False, round=True, print_test = False): def gws_by_ndars(edges, number, realize, round = True, sum = True, print_test = False): - """ More efficient way to calculate strain from numbered - grid integrated + r"""More efficient way to calculate strain from numbered grid integrated Parameters ---------- @@ -1092,14 +1091,14 @@ def gws_by_ndars(edges, number, realize, round = True, sum = True, print_test = Specification of how to construct one or more discrete realizations. If a `bool` value, then whether or not to construct a realization. If an `int` value, then how many discrete realizations to construct. - round : bool + round : bool, optional (default=True) Specification of whether to discretize the sample if realize is False, by rounding number of binaries in each bin to integers. This has no impact if realize is true. NOTE: should add a warning if round and realize are both True - sum : bool + sum : bool, optional (default=True) Whether or not to sum the strain at a given frequency over all bins. - print_test : bool + print_test : bool, optional (default=False) Whether or not to print variable as they are calculated, for dev purposes. @@ -1211,8 +1210,7 @@ def gws_by_ndars(edges, number, realize, round = True, sum = True, print_test = def unrealized_ss_by_ndars(edges, number, realize, round = True, print_test = False): - """ More efficient way to calculate strain from numbered - grid integrated + r"""More efficient way to calculate strain from numbered grid integrated Parameters @@ -1224,19 +1222,15 @@ def unrealized_ss_by_ndars(edges, number, realize, round = True, print_test = Fa number : (M, Q, Z, F) ndarray of scalars The number of binaries in each bin of parameter space. This is calculated by integrating `dnum` over each bin. - realize : bool or int, + realize : bool or int Specification of how to construct one or more discrete realizations. If a `bool` value, then whether or not to construct a realization. If a `int` value, then how many discrete realizations to construct. - round : bool + round : bool, optional (default=True) Specification of whether to discretize the sample if realize is False, by rounding number of binaries in each bin to integers. Does nothing if realize is True. - ss : bool - Whether or not to separate the loudest single source in each frequency bin. - sum : bool - Whether or not to sum the strain at a given frequency over all bins. - print_test : bool + print_test : bool, optional (default=False) Whether or not to print variable as they are calculated, for dev purposes. @@ -1266,7 +1260,8 @@ def unrealized_ss_by_ndars(edges, number, realize, round = True, print_test = Fa format: [[M indices], [q indices], [z indices], [f indices]] - + Notes + ----- Potential BUG: In the unlikely scenario that there are two equal hsmaxes (at same OR dif frequencies), ssidx calculation will go wrong Could avoid this by using argwhere for each f_idx column separately. @@ -1274,9 +1269,10 @@ def unrealized_ss_by_ndars(edges, number, realize, round = True, print_test = Fa values for that hsmax and raises a warning/assertion error - TODO: Calculate sspar - TODO: Implement realizations - TODO: Implement not summing, or remove option + * TODO: Calculate sspar + * TODO: Implement realizations + * TODO: Implement not summing, or remove option + """ if(print_test): @@ -1465,8 +1461,7 @@ def max_test(hsmax, hsamp): print('max_test passed') def ssidx_test(hsmax, hsamp, ssidx, print_test): - """ - Test ssidx in hsamp gives the same values as hsmax + r"""Test ssidx in hsamp gives the same values as hsmax Parameters ---------- @@ -1474,7 +1469,10 @@ def ssidx_test(hsmax, hsamp, ssidx, print_test): Maximum strain amplitude of a single source at each frequency. hsamp : (M, Q, Z, F,) ndarray of scalar Strain amplitude of a source in each bin - ssidx : (F, 4) ndarray + ssidx : (5, F, R) ndarray + Single source indices. + print_test : bool + Whether to print output variables during test. """ @@ -1490,8 +1488,7 @@ def ssidx_test(hsmax, hsamp, ssidx, print_test): def ssnew_test(hsmax, hsamp, ssnew, print_test): - """ - Test ssnew in hsamp gives the same values as hsmax + r"""Test ssnew in hsamp gives the same values as hsmax Parameters ---------- @@ -1500,8 +1497,9 @@ def ssnew_test(hsmax, hsamp, ssnew, print_test): hsamp : (M, Q, Z, F,) ndarray of scalar Strain amplitude of a source in each bin ssnew : (4, F) ndarray - - + Unraveled single source indices. + print_test : bool + Whether to print details. """ maxes = hsamp[ssnew[0], ssnew[1], ssnew[2], ssnew[3]] if(print_test): @@ -1514,8 +1512,7 @@ def ssnew_test(hsmax, hsamp, ssnew, print_test): def number_test(num, bgnum, fobs, exname='', plot_test=False): - ''' - Plots num - bgnum, where number is the ndarray of + r"""Plots num - bgnum, where number is the ndarray of integer number of sources in each bin, i.e. after rounding or Poisson sampling @@ -1527,19 +1524,13 @@ def number_test(num, bgnum, fobs, exname='', plot_test=False): bgnum : (M, Q, Z, F) array number of background sources in each bin, after single source subtraction - fobs : (F) array + fobs : (F,) array frequencies of each F, for ax titles - exname : String + exname : String, optional (Default='') name of example - plot_test : Bool - whether or not to print values a - - - Returns - ----------- - None - - ''' + plot_test : Bool, optional (Default=False) + Whether or not to plot the results + """ if np.all(num%1 == 0) != True: warnings.warn("num contains at least one non-integer value") difs = num - bgnum assert len(difs[np.where(difs>0)]) == len(difs[0,0,0,:]), "More than one bin per frequency found with a single source subtracted." @@ -1554,7 +1545,7 @@ def number_test(num, bgnum, fobs, exname='', plot_test=False): # print(num[...,0].shape) for f in range(len(fobs)): ax[f].scatter(bins, (num[...,f] - bgnum[...,f])) - ax[f].set_title('$f_\mathrm{obs}$ = %dnHz' % (fobs[f]*10**9)) + ax[f].set_title(rf'$f_{{\mathrm{{obs}}}}$ = {fobs[f]*10e9:.0f}nHz') ax[f].set_xlabel('bin') fig.tight_layout() print('number test passed') @@ -1589,11 +1580,11 @@ def quadratic_sum_test(hc_bg, hc_ss, hc_tt, print_test): -def run_ndars_tests(edges,number, fobs, exname='', print_test=False, +def run_ndars_tests(edges, number, fobs, exname='', print_test=False, loop_comparison = True): - ''' - Call tests for some edges, number - Paramaters + r"""Call tests for some edges, number. + + Parameters ---------- edges : (4,) list of 1D arrays Mass, ratio, redshift, and frequency edges of bins @@ -1601,16 +1592,24 @@ def run_ndars_tests(edges,number, fobs, exname='', print_test=False, Number of binaries in each bin fobs : (F,) array of scalars Observed frequency bin centers - exname : String + exname : String, optional (Default='') Name of example (used for number plots) Returns - ------ - hsamp - hsmax - ssidx - bgnum - ''' + ------- + hc_bg : array_like + Background characteristic strain. + hc_ss : array_like + Single-source characteristic strain. + hsamp : array_like + Strain amplitude. + ssidx : array_like + Single-source indices. + hsmax : array_like + Maximum strain. + bgnum : array_like + Background count. + """ hc_bg, hc_ss, hsamp, ssidx, hsmax, bgnum = ss_by_ndars(edges, number, realize=False, round=True) max_test(hsmax, hsamp) @@ -1640,15 +1639,16 @@ def run_ndars_tests(edges,number, fobs, exname='', print_test=False, ################################################### def example(dur, cad, mtot, mrat, redz, print_test): - ''' + r"""Example for single sources. + 1) Choose the frequency bins at which to calculate the GWB, same as in semi-analytic-models.ipynb 2) Build Semi-Analytic-Model with super simple parameters 3) Get SAM edges and numbers as in sam.gwb() Parameters - ----------oiuyuiuytd + ---------- dur : scalar - Duration of observation in secnods (multiply by YR) + Duration of observation in seconds (multiply by YR) cad : scalar Cadence of observations in seconds (multiply by YR) mtot : (3,) list of scalars @@ -1657,7 +1657,8 @@ def example(dur, cad, mtot, mrat, redz, print_test): Min, max, and steps for mass ratio. redz : (3,) list of scalars Min, max, and steps for redshift. - print_test : + print_test : bool, optional (Default=False) + Whether or not to print values. Returns ------- @@ -1666,11 +1667,11 @@ def example(dur, cad, mtot, mrat, redz, print_test): total mass, mass ratio, redshift, and observer-frame orbital frequency. The length of each of the four arrays is M+1, Q+1, Z+1, F+1. number : (M, Q, Z, F) array - The number of binaries in each bin of parameter space. This is calculated by integrating + The number of binaries in each bin of parameter space. This is calculated by integrating `dnum` over each bin. - fobs : (F) array - observed frequency bin centers - ''' + fobs : (F,) array + Observed frequency bin centers + """ # 1) Choose the frequency bins at which to calculate the GWB, same as in semi-analytic-models.ipynb fobs = utils.nyquist_freqs(dur,cad) fobs_edges = utils.nyquist_freqs_edges(dur,cad) @@ -1710,20 +1711,27 @@ def example(dur, cad, mtot, mrat, redz, print_test): def example2(print_test = True, exname='Example 2'): - ''' + r"""Run Example 2. + Parameters - --------- - print_test : Bool + ---------- + print_test : Bool, optional (Default=True) Whether to print frequencies and edges + exname : String, optional (Default='Example 2') + Name of example Returns - --------- + ------- edges : (M,Q,Z,F) array + Edges. number : (M, Q, Z, F) array + Number counts. fobs : (F) array - observed frequency bin centers - ''' + Observed frequency bin centers. + exname : str + Name of example. + """ dur = 5.0*YR/3.1557600 cad = .5*YR/3.1455145557600 @@ -1736,20 +1744,24 @@ def example2(print_test = True, exname='Example 2'): return edges, number, fobs, exname def example3(print_test = True, exname = 'Example 3'): - ''' + r"""Run Example 3. + Parameters - --------- - print_test : Bool + ---------- + print_test : Bool, optional (Default=True) Whether to print frequencies and edges - + exname : str, optional (Default='Example 3') + Name of example Returns - --------- - edges : (M,Q,Z,F) array + ------- + edges : (M, Q, Z, F) array number : (M, Q, Z, F) array fobs : (F) array observed frequency bin centers - ''' + exname : str + Name of example. + """ dur = 5.0*YR/3.1557600 cad = .5*YR/3.1557600 @@ -1763,20 +1775,26 @@ def example3(print_test = True, exname = 'Example 3'): def example4(print_test = True, exname = 'Example 4'): - ''' - Parameters - --------- - print_test : Bool - Whether to print frequencies and edges + r"""Run Example 4. + Parameters + ---------- + print_test : Bool, optional (Default=True) + Whether to print frequencies and edges. + exname : str, optional (Default='Example 4') + Name of example. Returns - --------- - edges : (M,Q,Z,F) array + ------- + edges : (M, Q, Z, F) array + Edges. number : (M, Q, Z, F) array + Number counts. fobs : (F) array - observed frequency bin centers - ''' + Observed frequency bin centers. + exname : str + Name of example. + """ dur = 5.0*YR/3.1557600 cad = .2*YR/3.1557600 @@ -1789,20 +1807,26 @@ def example4(print_test = True, exname = 'Example 4'): def example5(print_test = True, exname = 'Example 5'): - ''' - Parameters - --------- - print_test : Bool - Whether to print frequencies and edges + r"""Run Example 5. + Parameters + ---------- + print_test : Bool, optional (Default=True) + Whether to print frequencies and edges. + exname : str, optional (Default='Example 5') + Name of example. Returns - --------- - edges : (M,Q,Z,F) array + ------- + edges : (M, Q, Z, F) array + Edges. number : (M, Q, Z, F) array + Number counts. fobs : (F) array - observed frequency bin centers - ''' + Observed frequency bin centers. + exname : str + Name of example. + """ dur = 10.0*YR cad = .2*YR @@ -1821,15 +1845,28 @@ def example5(print_test = True, exname = 'Example 5'): def plot_medians(ax, xx, BG=None, SS=None, LABEL='', BG_COLOR='k', BG_ERRORS = False, SS_COLOR='k', SS_ERRORS = True): - """ - Parameters: - ax - xx - BG - SS - REALS - label - COLOR + r"""Plot median lines with optional errorbars. + + Parameters + ---------- + ax : pyplot ax object + Axes on which to plot. + xx : (F,) array_like + Frequencies. + BG : (F, R) array_like, optional + Background values across realizations. + SS : (F, R) array_like, optional + Single source values across realizations. + LABEL : str, optional + Label extension string. + BG_COLOR : str, optional + Background plot color. + BG_ERRORS : bool, optional + Whether to display background standard deviation errorbars. + SS_COLOR : str, optional + Single-source plot color. + SS_ERRORS : bool, optional + Whether to display single-source standard deviation errorbars. """ # plot median bg of samples if(BG is not None and BG_ERRORS == False): @@ -1853,22 +1890,25 @@ def plot_medians(ax, xx, BG=None, SS=None, LABEL='', def plot_BG(ax, xx, BG, LABEL, REALS=0, median=False, COLOR='b', rand=False): - """ - Plot the background median, middle 50% and 98% confidence intervals, - and optionally several realizations. All plotted in same color. + r"""Plot the background median, middle 50% and 98% confidence intervals. Parameters - ------------ + ---------- ax : pyplot ax object + Axes to draw on. xx : (F,) 1darray of scalars + Frequencies. BG : (F,R) Ndarray of scalars + Background characteristic strain arrays. LABEL : String - REALS : int - How many bg realizations to plot - median : bool - Whether or not to plot the bg median (same color) - COLOR : String - rand : bool + Label description. + REALS : int, optional (Default=0) + How many bg realizations to plot. + median : bool, optional (Default=False) + Whether or not to plot the bg median (same color). + COLOR : String, optional (Default='b') + Plotting color. + rand : bool, optional (Default=False) Whether to randomize which realizations are plotted (True) or plots the 0th to REALS-th realizations. """ @@ -1891,18 +1931,24 @@ def plot_BG(ax, xx, BG, LABEL, REALS=0, median=False, COLOR='b', rand=False): def plot_samples(ax, xx, BG=None, SS=None, REALS=1, LABEL=''): - """ - Plot the background and/or single sources for the first 'REALS' + r"""Plot the background and/or single sources for the first 'REALS' number of realizations, with each color corresponding to a difference realization. Parameters ---------- ax : pyplot ax object + Axes object. xx : (F,) array of scalars - BG : (F,R) ndarray or None - SS : (F,R) ndarray or None - REALS : int + Frequencies. + BG : (F, R) ndarray, optional + Background realizations. + SS : (F, R) ndarray, optional + Single source realizations. + REALS : int, optional + Number of realizations to plot. + LABEL : str, optional + Label extension suffix. """ colors = cm.rainbow(np.linspace(0,1,REALS)) for rr in range(REALS): @@ -1921,22 +1967,30 @@ def plot_samples(ax, xx, BG=None, SS=None, REALS=1, LABEL=''): edgecolor='k', alpha=0.5) def plot_std(ax, xx, BG, SS, COLOR='b', LABEL=''): - """ - Plot the standard deviations of the bg and ss characteristic strains. + r"""Plot standard deviations of background and single sources characteristic strain. Parameters - -------- + ---------- ax : pyplot ax object + Axes to draw on. xx : (F,) array of scalars + Frequencies. BG : (F, R) Ndarray of scalars + Background strain values. SS : (F, R) Ndarray of scalars - COLOR : string - LABEL : string + Single source strain values. + COLOR : string, optional (Default='b') + Color. + LABEL : string, optional (Default='') + Label suffix. + Returns ------- std_bg : (F,) array of scalars + Background standard deviations. std_ss : (F,) array of scalars + Single-source standard deviations. """ std_bg = np.std(BG, axis=1) # med_bg = np.median(BG, axis=1) @@ -1950,20 +2004,22 @@ def plot_std(ax, xx, BG, SS, COLOR='b', LABEL=''): return std_bg, std_ss def plot_IQR(ax, xx, BG=None, SS=None, COLOR='r', LABEL=''): - """ - Plot the IQR of the bg and ss characteristic strains. + r"""Plot the IQR of the background and single source characteristic strains. Parameters ---------- ax : pyplot ax object + Axes to draw on. xx : (F,) array of scalars - BG : (F, R) Ndarray of scalars - SS : (F, R) Ndarray of scalars - COLOR : string - LABEL : string - - Returns - ------- + Frequencies. + BG : (F, R) Ndarray of scalars, optional + Background values. + SS : (F, R) Ndarray of scalars, optional + Single source values. + COLOR : string, optional (Default='r') + Color of lines. + LABEL : string, optional (Default='') + Label extension suffix. """ if (BG is not None): @@ -1983,22 +2039,32 @@ def plot_percentiles(ax, xx, BG=None, SS=None, LABEL='', BG_COLOR='b', SS_COLOR='r', BG_MARKER=None, SS_MARKER=None, BG_LINESTYLE='solid', SS_LINESTYLE='solid'): - """ - Plots 25th and 75th percentiles, and IQR region between - (50% confidence interval). + r"""Plot 25th and 75th percentiles, and IQR region between. Parameters - --------- + ---------- ax : pyplot ax object + Axes object. xx : (F,) 1darray of scalars - BG : (F, R) Ndarray of scalars or None - SS : (F, R) Ndarray of scalars or None - BG_COLOR : string - SS_COLOR : string - BG_MARKER : string - SS_MARKER : string - BG_LINESTYLE : string - SS_LINESTYLE : string + Frequencies. + BG : (F, R) Ndarray of scalars, optional + Background realizations. + SS : (F, R) Ndarray of scalars, optional + Single-source realizations. + LABEL : str, optional + Label extension. + BG_COLOR : str, optional + Background line color. + SS_COLOR : str, optional + Single source line color. + BG_MARKER : str, optional + Background marker. + SS_MARKER : str, optional + Single source marker. + BG_LINESTYLE : str, optional + Background line style. + SS_LINESTYLE : str, optional + Single source line style. """ if (BG is not None): @@ -2027,23 +2093,46 @@ def plot_params(axs, xx, REALS=1, LABEL='', grid=None, BG_MEDIAN=True, SS_MEDIAN=True, BG_ERRORS=True, SS_ERRORS=True, BG_COLOR='k', SS_COLOR='mediumorchid', - TITLES = np.array([['Total Mass $M/M_\odot$', 'Mass Ratio $q$'], - ['Redshift $z$', 'Characteristic Strain $h_c$']]), - XLABEL = 'Frequency $f_\mathrm{obs}$ (1/yr)', + TITLES = np.array([[r'Total Mass $M/M_\odot$', r'Mass Ratio $q$'], + [r'Redshift $z$', r'Characteristic Strain $h_c$']]), + XLABEL = r'Frequency $f_{\mathrm{obs}}$ (1/yr)', SHOW_LEGEND = True): - """ - Plot mass, ratio, redshift, and strain in 4 separate subplots. + r"""Plot mass, ratio, redshift, and strain in 4 separate subplots. - Parameters: - ----------- + Parameters + ---------- axs : (2,2) array of pyplot ax object + Axes grid. xx : (F,) 1d array of scalars - params : (4,) 1Darray of (F,R,) NDarrays - titles : (4,) array of strings - xlabel : string - legend : bool - Whether or not to include a legend in each subplot - + Frequencies. + REALS : int, optional + Number of realizations to plot. + LABEL : str, optional + Label extension suffix. + grid : array_like, optional + Grid boundary values. + BG_PARAMS : array_like, optional + Background parameters. + SS_PARAMS : array_like, optional + Single source parameters. + BG_MEDIAN : bool, optional (Default=True) + Whether to plot the median background values. + SS_MEDIAN : bool, optional (Default=True) + Whether to plot the median single-source values. + BG_ERRORS : bool, optional (Default=True) + Whether to display background standard deviation limits. + SS_ERRORS : bool, optional (Default=True) + Whether to display single-source standard deviation limits. + BG_COLOR : str, optional (Default='k') + Background color line value. + SS_COLOR : str, optional (Default='mediumorchid') + Single source color value. + TITLES : array_like, optional + Label values for axes. + XLABEL : str, optional (Default=r'Frequency $f_{\mathrm{obs}}$ (1/yr)') + Independent variable label name. + SHOW_LEGEND : bool, optional (Default=True) + Whether to show the legend in each subpanel. """ colors = cm.rainbow(np.linspace(0,1,REALS)) for ii in range(len(axs)): @@ -2148,8 +2237,15 @@ def plot_params(axs, xx, REALS=1, LABEL='', grid=None, def threshold_hc(): - """ - Rosado+ 2015, SNR calculation + r"""[Rosado2015]_ SNR calculation. + + Returns + ------- + int: 0 + Hardcoded to return 0 for now. + + Notes + ----- S := cross correlation between pulsars S = \int_{-T/2}^{T/2} dt \int dt' s_i(t) s_j(t') Q(t,t') where T is the observation time, s_i(t) and s_j(t) are the data @@ -2162,25 +2258,47 @@ def threshold_hc(): def ss_occurence_rate(hc_ss, S_T): - """ - Parameters: - ------------ + r"""Calculate single-source occurence rate. + + Parameters + ---------- hc_ss : (F, R) NDarray of scalars characteristic strain of single sources S_T : float threshold strain? + Returns + ------ + int: 0 + Hardcoded to return 0 for now. """ return 0 def bg_occurence_rate(): - """ - According to Rosado+ 2015 - Universe may contain GWB if S >= S_T + r"""Calculate background occurence rate. + + Returns + ------ + int : 0 + Hardcoded to return 0 for now. + + Notes + ----- + According to [Rosado2015] + Universe may contain GWB if `S \ge S_T` + """ + return 0 def false_alarm_probability(): - """ + r"""Placeholder for false alarm probability. + + Returns + ------- + prob : float + + Notes + ----- TODO """ return 0 \ No newline at end of file diff --git a/holodeck/utils.py b/holodeck/utils.py index 4e50f340..7b086857 100644 --- a/holodeck/utils.py +++ b/holodeck/utils.py @@ -935,7 +935,7 @@ def gaussian_freqs(num, fmax, dur=16.03*YR): The number of frequency bins, `F` is the argument `num`. edges : (F+1,) ndarray Bin-edge frequencies for `F` bins, i.e. `F+1` bin edges. The frequency bin edges are at: - ``F_i = (1/dur) + (i+1) * (fmax - 1/dur)/(num-1) `` for i between 0 and `num`. + ``F_i = (1/dur) + (i+1) * (fmax - 1/dur)/(num-1)`` for i between 0 and `num`. The number of frequency bins, `F` is the argument `num`. """ @@ -2072,7 +2072,7 @@ def angs_from_sepa(sepa, dcom, redz): return angs def eddington_accretion(mass, eps=0.1): - """Eddington Accretion rate, $\\dot{M}_{Edd} = L_{Edd}/\\epsilon c^2$. + r"""Eddington Accretion rate, $\dot{M}_{Edd} = L_{Edd}/\epsilon c^2$. Arguments ---------