Post-processing#

After generating observation plans, post-processing computes detection probabilities and extracts key optimization metrics for each event.

Kilonova detection rate and Statistics

Calculating poisson-lognormal rate quantiles and formatting

python

This section provides the script used for generating and calculating the detection rates.

Show Script

import warnings
import numpy as np
from astropy import units as u
from astropy.table import QTable
from scipy import stats

from earthorbitplan.probability.rate import poisson_lognormal_rate_quantiles
from earthorbitplan.utils.path import get_project_root

# ------------------------------------------------------------------------------
# Suppress known warnings for cleaner output
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
warnings.filterwarnings("ignore", ".*dubious year.*")
warnings.filterwarnings(
    "ignore", "Tried to get polar motions for times after IERS data is valid.*"
)

# ------------------------------------------------------------------------------
# Load simulated event data
root = get_project_root()
events_file = root / "data" / "events.ecsv"
main_table = QTable.read(events_file)

# Get unique run names
runs = np.unique(main_table["run"])

# Filter events by objective_value cutoff
cutoff = main_table["cutoff"][0]
main_table = main_table[main_table["objective_value"] >= cutoff]

# Group events by run
event_tables_by_run = {run: main_table[main_table["run"] == run] for run in runs}

# ------------------------------------------------------------------------------
# Set merger rate priors from O3 R&P Table II (last column)
lo, mid, hi = 100, 240, 510  # In Gpc^-3 yr^-1

# Log-normal width for 90% interval
(standard_90pct_interval,) = np.diff(stats.norm.interval(0.9))
log_target_rate_mu = np.log(mid)
log_target_rate_sigma = np.log(hi / lo) / standard_90pct_interval

# Get effective rate for each run
log_simulation_effective_rate_by_run = {
    key: np.log(value.to_value(u.Gpc**-3 * u.yr**-1))
    for key, value in main_table.meta["effective_rate"].items()
}

# ------------------------------------------------------------------------------
# Compute median and quantiles for each run
prob_quantiles = np.asarray([0.5, 0.05, 0.95])  # Median, 5%, 95%
run_duration = 1.5  # years

mu = np.asarray(
    [
        log_target_rate_mu
        + np.log(run_duration)
        - log_simulation_effective_rate_by_run[run]
        + np.log(
            [
                np.sum(_)
                for _ in [
                    np.ones_like(event_tables_by_run[run]["objective_value"]),
                    event_tables_by_run[run]["detection_probability_known_position"],
                ]
            ]
        )
        for run in runs
    ]
)

# Compute Poisson-Lognormal rate quantiles for all runs
rate_quantiles = poisson_lognormal_rate_quantiles(
    prob_quantiles[np.newaxis, np.newaxis, :],
    mu.T[:, :, np.newaxis],
    log_target_rate_sigma,
)

# ------------------------------------------------------------------------------
# Utility: Format a table as reStructuredText grid table
def make_rst_table(headers, rows):
    columns = [headers] + rows
    n_cols = len(headers)
    col_widths = [max(len(str(row[i])) for row in columns) for i in range(n_cols)]

    def sep(char="+", fill="-"):
        return char + char.join(fill * (w + 2) for w in col_widths) + char

    def fmt_row(row):
        return (
            "| "
            + " | ".join(str(cell).ljust(w) for cell, w in zip(row, col_widths))
            + " |"
        )

    lines = [
        sep(),
        fmt_row(headers),
        sep("=", "="),
    ]
    for row in rows:
        lines.append(fmt_row(row))
        lines.append(sep())
    return "\n".join(lines)

# Example: Prepare headers and format quantile results
headers = ["Run"] + list(runs)
labels = ["Number of events selected", "Number of events detected"]
rst_rows = []

for label, row in zip(labels, rate_quantiles):
    formatted = [
        ":math:`{}_{{-{}}}^{{+{}}}`".format(*np.rint([mid, mid - lo, hi - mid]).astype(int))
        for mid, lo, hi in row
    ]
    rst_rows.append([label] + formatted)

rst_table = make_rst_table(headers, rst_rows)

# Print the table for RST documentation
print(rst_table)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[1], line 22
     18 # ------------------------------------------------------------------------------
     19 # Load simulated event data
     20 root = get_project_root()
     21 events_file = root / "data" / "events.ecsv"
---> 22 main_table = QTable.read(events_file)
     23 
     24 # Get unique run names
     25 runs = np.unique(main_table["run"])

File ~/checkouts/readthedocs.org/user_builds/earthorbitplan/envs/latest/lib/python3.11/site-packages/astropy/table/connect.py:62, in TableRead.__call__(self, *args, **kwargs)
     59 units = kwargs.pop("units", None)
     60 descriptions = kwargs.pop("descriptions", None)
---> 62 out = self.registry.read(cls, *args, **kwargs)
     64 # For some readers (e.g., ascii.ecsv), the returned `out` class is not
     65 # guaranteed to be the same as the desired output `cls`.  If so,
     66 # try coercing to desired class without copying (io.registry.read
     67 # would normally do a copy).  The normal case here is swapping
     68 # Table <=> QTable.
     69 if cls is not out.__class__:

File ~/checkouts/readthedocs.org/user_builds/earthorbitplan/envs/latest/lib/python3.11/site-packages/astropy/io/registry/core.py:200, in UnifiedInputRegistry.read(self, cls, format, cache, *args, **kwargs)
    196 try:
    197     ctx = get_readable_fileobj(
    198         args[0], encoding="binary", cache=cache
    199     )
--> 200     fileobj = ctx.__enter__()
    201 except OSError:
    202     raise

File ~/.asdf/installs/python/3.11.14/lib/python3.11/contextlib.py:137, in _GeneratorContextManager.__enter__(self)
    135 del self.args, self.kwds, self.func
    136 try:
--> 137     return next(self.gen)
    138 except StopIteration:
    139     raise RuntimeError("generator didn't yield") from None

File ~/checkouts/readthedocs.org/user_builds/earthorbitplan/envs/latest/lib/python3.11/site-packages/astropy/utils/data.py:364, in get_readable_fileobj(name_or_obj, encoding, cache, show_progress, remote_timeout, sources, http_headers, use_fsspec, fsspec_kwargs, close_files)
    355 if is_url:
    356     name_or_obj = download_file(
    357         name_or_obj,
    358         cache=cache,
   (...)    362         http_headers=http_headers,
    363     )
--> 364 fileobj = io.FileIO(name_or_obj, "r")
    365 if is_url and not cache:
    366     delete_fds.append(fileobj)

FileNotFoundError: [Errno 2] No such file or directory: '/home/docs/checkouts/readthedocs.org/user_builds/earthorbitplan/checkouts/latest/data/events.ecsv'

Display Detection Rate Table

Rate of Selected and Detected Events

Run	O5	O6
Number of events selected	\(43_{-26}^{+56}\)	\(55_{-33}^{+72}\)
Number of events detected	\(19_{-12}^{+26}\)	\(25_{-16}^{+34}\)

Note

This table is just an example. If the simulation file events.ecsv is updated, the table generated by the script above may differ from this static version. This table does not update automatically; it is provided for a better overview only.

Notebook