add prost

README.md

@@ -29,11 +29,21 @@ This uses the Babamul alert broker to fetch alerts, which requires user credenti
 To fetch alerts, crossmatch against available catalogs, and save candidates, run:
 
 ```bash
-python src/thor/crossmatch_alerts.py --start 06-28-2026 --end 06-30-2026 --additional_filtering tde_filter --scan
+python src/thor/crossmatch_alerts.py --start 06-28-2026 --end 06-30-2026
 ```
 
 To apply additional TDE-specific filtering, pass `--additional_filtering tde_filter`. Add flag `--save_raw_alerts` to save raw alerts to `data/lsst_alert_download/raw_files/`, and flag `--save_results` in order to save the crossmatch details locally to `data/lsst_alert_download/`. A summary of results will be printed in command line, but the `--scan` flag can also be included to open a temp jupyter notebook in browser and use Babamul's scanning tool.
 
+Pass `--method prost` to use probabilistic host association via [astro_prost](https://github.com/alexandergagliano/galaxy-association) instead of cone search. For efficiency, this will run a 10 arcsec consearch first to remove alerts with no nearby hosts.
+
+```bash
+python src/thor/crossmatch_alerts.py --start 06-28-2026 --end 06-30-2026 --method prost --additional_filtering tde_filter --scan
+```
+
+Results are saved as a `.csv` (one row per transient) rather than `.json`. Transients with a host posterior > 0.3 are printed to the command line.
+
+The script currently uses galaxy offset as the default scoring property, with a `uniform(0, 10)` prior and `gamma(a=0.75)` likelihood. Redshift information is loaded from catalogs where available and can be included in scoring via custom priors. For more flexibility, see the Prost section of the [crossmatching notebook](docs/notebooks/Example_LSST_Catalog_Crossmatching.ipynb) or the [astro_prost documentation](https://github.com/alexandergagliano/galaxy-association) for the full list of supported properties and association options.
+
 
 ### Data
 

docs/notebooks/Example_LSST_Catalog_Crossmatching.ipynb

@@ -591,6 +591,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# test Prost out of the box / basic usage\n",
+    "\n",
     "import pandas as pd\n",
     "from astro_prost import associate_sample\n",
     "from scipy.stats import gamma, halfnorm, uniform\n",
@@ -603,7 +605,7 @@
     "})\n",
     "\n",
     "# define a set of catalogs to search -- options are glade, decals, panstarrs, and skymapper\n",
-    "catalogs = [\"decals\"]\n",
+    "catalogs = [\"decals\", \"glade\"]\n",
     "\n",
     "# define priors and likelihoods\n",
     "priorfunc_offset = uniform(loc=0, scale=10)\n",
@@ -632,14 +634,37 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from pathlib import Path                                                                                                                                                    \n",
+    "# test prost with custom patch to handle unsuported catalogs\n",
+    "\n",
+    "from pathlib import Path    \n",
+    "import pandas as pd                                                                                                                                                \n",
     "from astropy.table import Table                                                                                                                                             \n",
-    "from thor.utils.prost_catalogs import register_local_catalog                                                                                                                \n",
-    "                                                                                                                                                                            \n",
-    "CATALOG_DIR = Path.cwd().parents[1] / \"data\" / \"catalogs\"                                                                                                                   \n",
-    "                                                                                                                                                                            \n",
+    "from thor.utils.prost_catalogs import register_local_catalog\n",
+    "from astro_prost import associate_sample\n",
+    "from scipy.stats import gamma, halfnorm, uniform\n",
+    "\n",
+    "\n",
+    "\n",
+    "# custom local catalogs                                                                                                              \n",
+    "CATALOG_DIR = Path.cwd().parents[1] / \"data\" / \"catalogs\"\n",
     "df = Table.read(CATALOG_DIR / \"COSMOS2025_cut.fits\").to_pandas()                                                                                                            \n",
-    "register_local_catalog(\"cosmos2025\", df)  \n",
+    "register_local_catalog(\"cosmos2025\", df)\n",
+    "\n",
+    "# define a transient catalog \n",
+    "transient_catalog = pd.DataFrame({\n",
+    "    'name': ['MyTransient'],\n",
+    "    'ra': [150.1],\n",
+    "    'dec': [2.2]\n",
+    "})\n",
+    "\n",
+    "# define a set of catalogs to search -- options are glade, decals, panstarrs, and skymapper\n",
+    "catalogs = [\"cosmos2025\"]\n",
+    "\n",
+    "# define priors and likelihoods\n",
+    "priorfunc_offset = uniform(loc=0, scale=10)\n",
+    "likefunc_offset = gamma(a=0.75)\n",
+    "priors = {\"offset\": priorfunc_offset}\n",
+    "likes = {\"offset\": likefunc_offset}\n",
     "\n",
     "# associate\n",
     "hosts = \\\n",

src/thor/crossmatch_alerts.py

@@ -1,5 +1,6 @@
 import argparse
 import json
+import math
 import os
 import subprocess
 import tempfile
@@ -109,6 +110,72 @@ def _print_match_report(crossmatched_objects):
     print(divider)
 
 
+def _print_prost_report(results_df):
+    matched = results_df[results_df['best_cat'].notna()]
+    n = len(matched)
+    print(f"\nMatched Object IDs: {n}")
+
+    if n > 100:
+        print(f"More than 100 matches ({n} total) — skipping per-object summary.")
+        return
+
+    if n == 0:
+        return
+
+    def _isnan(v):
+        try:
+            return math.isnan(v)
+        except (TypeError, ValueError):
+            return v is None
+
+    id_w   = max(len("LSST Object ID"), matched['name'].astype(str).str.len().max())
+    cat_w  = 30
+    z_w    = 8
+    sep_w  = 10
+    post_w = 10
+
+    header = (
+        f"{'LSST Object ID':<{id_w}}  "
+        f"{'Catalog':<{cat_w}}  "
+        f"{'z':>{z_w}}  "
+        f"{'Sep (\")'  :>{sep_w}}  "
+        f"{'Posterior':>{post_w}}"
+    )
+    divider = "-" * len(header)
+    print(divider)
+    print(header)
+    print(divider)
+
+    for _, row in matched.iterrows():
+        first = True
+        for prefix in ("host", "host_2"):
+            objid = row.get(f"{prefix}_objID")
+            post_val = row.get(f"{prefix}_total_posterior")
+            if objid is None or _isnan(objid):
+                continue
+            if _isnan(post_val) or post_val < 0.3:
+                continue
+            z_val   = row.get(f"{prefix}_redshift_mean")
+            sep_val = row.get(f"{prefix}_offset_mean")
+            cat     = row.get("best_cat", "—")
+
+            z_str    = f"{z_val:.3f}"    if not _isnan(z_val)    else "—"
+            sep_str  = f"{sep_val:.2f}"  if not _isnan(sep_val)  else "—"
+            post_str = f"{post_val:.3f}" if not _isnan(post_val) else "—"
+            id_str   = str(row['name']) if first else ""
+
+            print(
+                f"{id_str:<{id_w}}  "
+                f"{cat:<{cat_w}}  "
+                f"{z_str:>{z_w}}  "
+                f"{sep_str:>{sep_w}}  "
+                f"{post_str:>{post_w}}"
+            )
+            first = False
+
+    print(divider)
+
+
 def main():
     dotenv.load_dotenv()
     parser = argparse.ArgumentParser(description="Fetch LSST alerts and crossmatch against catalogs.")
@@ -135,6 +202,12 @@ def main():
         action="store_true",
         help="Open a temporary Jupyter notebook to scan candidates with scan_alerts (default: off).",
     )
+    parser.add_argument(
+        "--method",
+        default="conesearch",
+        choices=["conesearch", "prost"],
+        help="Crossmatch method: conesearch (default) or prost (probabilistic host association).",
+    )
     args = parser.parse_args()
 
     # ── Fetch alerts ──────────────────────────────────────────────────────────
@@ -160,17 +233,51 @@ def main():
     # ── Deduplicate to unique objects ─────────────────────────────────────────
     filtered_objects = filter_functions.deduplicate_alerts(filtered_alerts)
 
+    # ── Optional additional filtering (alert-based, method-independent) ───────
+    if args.additional_filtering == "tde_filter":
+        filtered_objects = filter_functions.filter_alerts(
+            filtered_objects,
+            filter_functions.tde_filter,
+        )
+
     # ── Crossmatch against all available catalogs ─────────────────────────────
+    if args.method == "prost":
+        # Pre-filter with a fast cone search to avoid running prost on alerts
+        # with no catalog coverage at all.
+        cone_matches = filter_functions.catalog_crossmatch(
+            alerts=filtered_objects,
+            method="conesearch",
+            radius_arcsec=10.0,
+        )
+        n_before = len(filtered_objects)
+        filtered_objects = [a for a in filtered_objects if a.objectId in cone_matches]
+        n_after = len(filtered_objects)
+        print(f"10\" pre-filter: cut {n_before - n_after} candidates, running prost on {n_after}.")
+
+    import time
+    _t0 = time.monotonic()
     crossmatched_objects = filter_functions.catalog_crossmatch(
         alerts=filtered_objects,
+        method=args.method,
     )
+    _elapsed_min = (time.monotonic() - _t0) / 60
+    print(f"\nRan crossmatch in {_elapsed_min:.1f} minutes on {len(filtered_objects):,} alerts from {args.start} to {args.end}.")
 
-    # ── Optional additional filtering ─────────────────────────────────────────
-    if args.additional_filtering == "tde_filter":
-        crossmatched_objects = filter_functions.filter_alerts(
-            crossmatched_objects,
-            filter_functions.tde_filter,
-        )
+    timestamp = datetime.now(tz=timezone.utc).strftime("%Y%m%dT%H%M%SZ")
+    repo_root = Path(__file__).resolve().parents[2]
+    out_dir = repo_root / "data" / "lsst_alert_download"
+
+    # ── prost returns a DataFrame; handle separately ──────────────────────────
+    if args.method == "prost":
+        _print_prost_report(crossmatched_objects)
+        if args.save_result:
+            out_dir.mkdir(parents=True, exist_ok=True)
+            out_file = out_dir / f"crossmatch_candidates_{timestamp}.csv"
+            crossmatched_objects.to_csv(out_file, index=False)
+            print(f"\nSaved {len(crossmatched_objects):,} prost results to {out_file}")
+        if args.scan:
+            _launch_scan_notebook(crossmatched_objects['name'].tolist())
+        return
 
     # ── Report ────────────────────────────────────────────────────────────────
     if not crossmatched_objects:
@@ -192,9 +299,6 @@ def main():
         for obj_id, obj in crossmatched_objects.items()
     }
 
-    timestamp = datetime.now(tz=timezone.utc).strftime("%Y%m%dT%H%M%SZ")
-    repo_root = Path(__file__).resolve().parents[2]
-    out_dir = repo_root / "data" / "lsst_alert_download"
     out_dir.mkdir(parents=True, exist_ok=True)
     out_file = out_dir / f"crossmatch_candidates_{timestamp}.json"
 

src/thor/summarize_rubin_alerts.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 """
 Fetch and visualise LSST alerts for a given date range.