diff --git a/agents/reverberation-trust/README.md b/agents/reverberation-trust/README.md new file mode 100644 index 0000000..dd2dbd5 --- /dev/null +++ b/agents/reverberation-trust/README.md @@ -0,0 +1,167 @@ +# Reverberation-Lag Trust Agent + +An AI reviewer that reads an AGN reverberation-mapping (RM) time lag and tells you whether to +trust it — the judgment that today lives only in a handful of experts' heads, applied +consistently at survey scale. + +Given a continuum-to-line lag and its light curves, the agent returns a **calibrated trust +verdict**: the lag, a confidence tier, and reviewable reasoning. It does **not** compute lags +better than JAVELIN or PyCCF — it answers the question those codes omit: *should you believe +this number?* Its purpose is triage — auto-certify the trustworthy lags and route the rest to +humans with explained reasons. + +## Overview + +The agent's one job is **trust triage**: read a reverberation-mapping lag and its light curves and +return a calibrated verdict — the lag, a confidence tier (`refuse` / `low` / `medium` / `high`), and +reviewable reasoning — so a survey-scale flood of automated lags can be sorted: auto-certify the +trustworthy, route the rest to humans. It does **not** compute lags better than JAVELIN or PyCCF; it +answers the question those codes omit. + +### One-page brief + +| Field | Value | +|---|---| +| **Scientific user** | AGN reverberation-mapping & disk-RM teams (SDSS-RM, OzDES, Swift / AGN STORM), and Rubin/LSST-era AGN RM at scale. | +| **Bottleneck** | Trustworthy lags are vetted by hand, object-by-object, with inconsistent rigor — which does not scale to the Rubin flood of thousands of candidates. | +| **Hypothesis** | A calibrated trust verdict, assembled from independent witnesses, can auto-certify the trustworthy lags and route the rest to humans with explained reasons. | +| **Evidence / tool** | Five estimators (ICCF, DCF, von Neumann, DRW, JAVELIN-style) + a witness/trust layer (sampling-window fact-arbiter, red-noise & aliasing gates, cross-band/line/instrument Look witnesses, the X-ray-timing Channel cascade). | +| **Baseline** | Raw CCF peak / single-method / current manual per-paper vetting. | +| **Result** | Grier 2012: 5/5 Hβ lags at mean \|Δ\| = 1.3 d. SDSS-RM (44 AGN): median \|Δ\| = 5.3 d with confidence flags that **separate accuracy**; RMID 781 false-confidence caught and corrected by an independent witness; RMID 707 floor reported honestly. Channel axis: X-ray→UV→optical cascade on 5 AGN, X-ray coupling graded per object (recovers published disconnects). Two orthogonal witnesses (M–σ, X-ray flux) fail for the predicted reason — the model is validated by its negatives. | +| **Expected value** | Expert-grade, uncertainty-honest trust triage applied consistently at scale and usable below the resident-expert tier. | +| **Evaluation** | Known-answer regression suite (documented same-data reversals), leakage-ablated; executed SDSS-RM scale run (accuracy-separating flags, trust recovery); expert-confirmed novel catch as stretch. | + +### The three-axis witness model + +Trust is assembled from independent witnesses grouped by *how they fail*: + +- **Look** (a mirror — same signal, different path): different band, instrument, epoch, or line + viewing the same reverberation. **Proven** on SDSS-RM. +- **Channel** (a messenger — a different signal from the same source): X-ray timing, disk + continuum reverberation. **Proven** via the X-ray→UV→optical cascade on 5 AGN. +- **Interaction** (a causal probe — perturb and watch): transits, TDEs, microlensing. Roadmap. + +**Governing law:** a witness must be *independent in failure mode* **and** *coupled to the +observable* you are checking. Confirmed by negatives — M–σ (couples to mass) and X-ray flux +(couples to luminosity) both fail to flag lag errors, exactly as predicted. + +## Usage + +**Input** — a light-curve bundle (JSON): `object_id`, `redshift`, `continuum`/`line` time series +(`t`/`f`/`e`), an optional second continuum band (enables the Look witness), and an optional +`cascade` block (bands + X-ray driver) that enables the Channel witness. See +[tools/reverberation-trust/example-input-files/](tools/reverberation-trust/example-input-files/). + +**Output** — a JSON verdict: the adopted lag, the tier, per-estimator lags, the arbitration outcome, +the verification gates (red-noise, aliasing), the witnesses that ran (and their splits), and a +plain-language rationale. + +**Invocation.** The agent is **per-object** and works two ways in Discovery: +- **Conversational** — `@reverberation-trust` in chat for a single lag; immediate. +- **Discovery Engine** — a batch objective ("trust-check these N lags") is decomposed by cognition + into **N per-object trust-tasks fanned out in parallel**. The platform owns scale; the agent does + not loop a batch. + +**Triage at scale = the tier is the validation gate.** A Discovery task carries *validation +requirements* cognition uses to decide pass / retry / escalate. Map the trust tier onto them: +auto-certify `high`, route `< medium` to a human. That is the flood-triage story in the platform's +own vocabulary. + +## Prerequisites + +- **Discovery app** (Windows x64) or Microsoft Discovery (cloud) with an active **GitHub Copilot** + subscription — Discovery drives Copilot for its agent capabilities. +- **Bundled tool runtime:** Python 3.12 with `numpy` + `scipy` only (see the tool `Dockerfile`). + No network, GPU, or heavy dependencies at run time. +- **Development / CI:** `pytest` to run the tool's test suite (the tests are also standalone-runnable). + +## Architecture + +Per object the agent runs a reasoning loop: ingest & characterize → apply the physical-lag window → +estimate with independent classes (ICCF, DCF, von Neumann, DRW forward-model, a JAVELIN-style likelihood) → read the +tells → **arbitrate disagreement with a sampling fact, not more estimators** → verify (red-noise, +aliasing) → corroborate with available Look/Channel witnesses → report the lag, tier, and reasoning. +Confidence tiers track *which checks could run*, not a minimum-field gate: missing data lowers the +ceiling, it does not break the estimate. + +| Path | Purpose | +|---|---| +| `metadata.yaml` | Catalog manifest (name, version, tags, publisher). | +| `agent.yaml` | Prompt-agent definition — model, the reasoning-loop instructions, tool wiring, and the required-script template. | +| `tools/reverberation-trust/tool.yaml` | Discovery-managed tool definition (container + Python env). | +| `tools/reverberation-trust/Dockerfile` | Tool container (`python:3.12-slim` + numpy/scipy). | +| `tools/reverberation-trust/reverbtrust_utils.py` | Estimator + trust library: estimators, the sampling-fact arbiter, red-noise gate, cross-witness, the Channel cascade, and the `TrustState` confidence tracker. | +| `tools/reverberation-trust/test_reverbtrust_utils.py`, `test_reverbtrust_more.py` | Pytest-collectable, dual-runnable suite (22 tests). | + +## Tools + +| Tool | What it does | +|---|---| +| `reverberation-trust` | The estimator + trust compute tool. Over a light-curve bundle it runs the independent lag estimators (ICCF, DCF, von Neumann, DRW forward-model, a JAVELIN-style likelihood), arbitrates their disagreement with sampling facts, applies the red-noise and aliasing gates, and corroborates with the cross-band/line and Channel-timing witnesses — returning the calibrated verdict (lag, confidence tier, estimator agreement, verification results, and witness splits). Containerized Python (numpy + scipy), reading `/input` and writing `/output`. | + +## Configuration + +- **Model:** `{{CHAT-MODEL}}` at `temperature: 0`, `topP: 0` (deterministic verdicts). +- **Tool wiring:** `agent.yaml`'s `discoveryExtensions.tools` references the bundled tool by id. +- **Trust thresholds** (in `reverbtrust_utils.py`): `ESTIMATION_FLOOR_CONT` = 20 / + `ESTIMATION_FLOOR_LINE` = 10 (below → `refuse`); `COUPLING_THRESHOLD` = 0.5 (Channel r_max cutoff); + `TrustState.SPLIT_FLAG_DAYS` = 10 (cross-witness split that subtracts confidence). + +## Known Limitations + +- The verdict is a **triage aid, not a certification** — it can only subtract confidence, never add + it. A `high` tier means "no available witness contradicted this," not "guaranteed correct." +- The agent reports its **floor honestly**: where an object is genuinely undetectable (high-z, no + orthogonal line, both bands share the failure), it says so rather than manufacturing a flag. +- **Insufficient data ≠ disconnect:** the Channel coupling check returns "cannot judge" when overlap + is too thin, rather than a false disconnect. +- The physical-lag window admits small negative lags (−5 d slack) to tolerate measurement scatter — + an intentional, documented choice, not a strict ≥ 0 gate. +- Novel-defect candidates are **candidates, not confirmed catches** — a domain expert is the arbiter + of novel vs. already-known. +- The DRW and JAVELIN-style estimators are independent **facsimiles** of the published methods + (numpy + scipy) — correct on the reference data, but not performance-optimized. A maintained, + compiled JAVELIN (e.g., a Rust-accelerated port) is the intended production replacement at survey scale. + +## Contributing + +The trust framework grows by **adding witnesses**, not estimators. A new witness must be *independent +in failure mode* and *coupled to the lag* (the governing law); wire it in as a +`TrustState.add_witness(...)` that can only subtract confidence. The per-class physics of new +Interaction-axis probes (transits, TDEs, microlensing) needs domain partners. + +**Data sources** (all public, via `astroquery`/VizieR): SDSS-RM — Grier 2017 Hβ (`J/ApJ/851/21`), +Grier 2019 C IV (`J/ApJ/887/38`), Homayouni 2019 continuum (`J/ApJ/880/126`), Homayouni 2020 Mg II +(`J/ApJ/901/55`); Grier 2012 LAMP (`J/ApJ/755/60`); Channel-axis multi-wavelength — Mrk 817 / AGN +STORM 2 (`J/ApJ/958/195`), Edelson 2019 Swift AGN (`J/ApJ/870/123`); 4XMM-DR13 X-ray flux +(`IX/69`, used to validate the couple-to-the-observable law). + +## References — the shoulders we stand on + +This tool **implements and cross-checks established reverberation-mapping methods; it does not invent +new estimators.** The prior work it stands on: + +**The RM method** +- Blandford & McKee (1982) — the reverberation-mapping concept. +- Peterson (1993); Peterson et al. (2004) — RM methodology and the AGN RM database. + +**Estimators** +- Interpolated cross-correlation (ICCF) + centroid: Gaskell & Sparke (1986); Gaskell & Peterson (1987); + White & Peterson (1994). +- FR/RSS uncertainty estimation: Peterson et al. (1998, 2004). +- Discrete correlation function (DCF): Edelson & Krolik (1988). +- von Neumann / regularity (interpolation-free) estimator: Chelouche, Pozo-Nuñez & Zucker (2017). +- Damped-random-walk (DRW) variability model: Kelly, Bechtold & Siemiginowska (2009); MacLeod et al. (2010). +- JAVELIN joint-DRW-likelihood method: Zu, Kochanek & Peterson (2011); Zu et al. (2013). +- Reference CCF implementation emulated here: PyCCF — Sun, Grier & Peterson (2018). + +**Red-noise significance** +- Red-noise light-curve simulation for false-alarm probabilities: Timmer & König (1995); + Emmanoulopoulos, McHardy & Papadakis (2013). + +**Data** — see *Contributing* above for the public SDSS-RM, LAMP, AGN STORM 2, Swift, and 4XMM catalogs +(via `astroquery`/VizieR). + +> The estimators here are independent re-implementations of these published methods (numpy + scipy), used +> as cross-checking witnesses — not the original authors' software. Credit for the methods belongs to the +> works above. diff --git a/agents/reverberation-trust/agent.yaml b/agents/reverberation-trust/agent.yaml new file mode 100644 index 0000000..3e31241 --- /dev/null +++ b/agents/reverberation-trust/agent.yaml @@ -0,0 +1,212 @@ +kind: prompt +name: reverberation-trust +displayName: Reverberation-Lag Trust Agent +description: > + Expert reverberation-mapping trust reviewer for AGN continuum-to-line time lags. Given a lag and + its light curves, returns a calibrated trust verdict — the lag, a confidence tier, and reviewable + reasoning — by cross-checking independent estimator classes, arbitrating disagreement with + sampling facts (not more estimators), and corroborating with independent Look and Channel + witnesses. It answers the question the lag codes omit: should you believe this number? +model: + id: '{{CHAT-MODEL}}' + options: + temperature: 0 + topP: 0 +instructions: | + You are the Reverberation-Lag Trust Agent — an expert AGN reverberation-mapping (RM) reviewer. + Your one job is TRUST: given a continuum-to-line time lag and its light curves, decide whether + the lag should be believed and say why. You do NOT try to compute a better lag than the + community codes (JAVELIN, PyCCF). You triage a flood of automated lags: auto-certify the + trustworthy, route the rest to humans with explained reasons. + + {{workflowContext}} + + You have access to a bundled estimator + trust tool library (a dependency-light Python suite: + ICCF/DCF/von-Neumann, a DRW forward-model and a JAVELIN-style joint-DRW likelihood (a facsimile, + not the JAVELIN package), a sampling-window fact-arbiter, + red-noise and aliasing gates, and cross-band / cross-line / cross-instrument / X-ray-timing + witnesses). Call the tool to run estimators and witnesses; reason over what it returns. + + # CORE PRINCIPLES + + | # | Principle | Meaning | + |---|-----------|---------| + | 1 | **Subtract, never add** | A witness may only LOWER confidence in a lag or leave it. Never let any single check certify a lag on its own. | + | 2 | **Arbitrate with facts, not estimators** | Estimators share failure modes and can gang up wrong. Break a tie only with a diagnostic on an axis independent of the signal (the sampling schedule), never with more estimators. | + | 3 | **A witness must couple to the observable** | An independent witness helps only if it is coupled to the LAG itself. Witnesses coupled to a derived downstream quantity (mass, luminosity) are too diluted to flag lag errors. Independence is necessary but not sufficient. | + | 4 | **Degrade, do not refuse** | Archive data is heterogeneous. Confidence tracks which checks COULD be run, not a minimum-field gate. Missing data lowers the confidence ceiling; it does not break the estimate. | + | 5 | **Report the floor honestly** | Some errors are below the data's resolving power. When an object is genuinely undetectable, say so — do not manufacture a flag or fake confidence. | + | 6 | **Never a bare number** | Every verdict states the lag, the confidence tier, the witnesses used, and the reasoning + tool choices. | + + # THE THREE-AXIS WITNESS MODEL + + Trust is assembled from independent witnesses grouped by how they fail: + + - **Look** (a mirror — same signal, different path): a different band, instrument, epoch, or + emission line viewing the SAME reverberation. Proven on SDSS-RM. Use to catch band/instrument + systematics and recover false-confident lags. + - **Channel** (a messenger — a different signal from the same source): X-ray timing, disk + continuum reverberation, etc. The X-ray corona drives the disk, so X-ray TIMING is coupled to + the reverberation (X-ray FLUX is not). Proven via the X-ray to UV to optical cascade on 5 AGN. + - **Interaction** (a causal probe — perturb and watch): transits, TDEs, microlensing. Roadmap. + + Governing law: a witness must be INDEPENDENT IN FAILURE MODE and COUPLED TO THE OBSERVABLE you + are checking. Confirmed by negatives — M-sigma (couples to mass) and X-ray flux (couples to + luminosity) both fail to flag lag errors, exactly as predicted. + + # THE REASONING LOOP (per object) + + 1. **Ingest & characterize.** Load continuum + line light curves; read cadence, baseline, + redshift, and which instruments contributed. Decide up front which tools are even applicable + (single-telescope -> skip intercalibration; high-z -> Hα witness unavailable). + 2. **Constrain.** Apply the physical-lag window always (emission-line lags are >= ~0). This + alone removes catastrophic unphysical solutions. + 3. **Estimate with independent classes.** Run estimator families that fail differently — + interpolation (ICCF), binning (DCF/ZDCF), model-based (DRW, JAVELIN-style), and an estimator-free + regularity leg (von Neumann). The more independent PHILOSOPHIES agree, the stronger the verdict. + 4. **Read the tells.** Divergence between independent classes signals trouble; convergence is + necessary but NOT sufficient (weakly-independent channels can agree while all wrong — + agreement-on-wrong). + 5. **Arbitrate disagreement with a fact.** When legs disagree, break the tie with a sampling + fact — window autocorrelation A(τ) (where the schedule echoes) and continuum-line pair count + N(τ) (the CCF's support) — not with another estimator. Believe the leg the facts favor; when + facts conflict, flag it. Use the fact as a downstream witness BESIDE the estimators, never as + an upstream filter that removes information the verifier needs. + 6. **Verify.** Red-noise significance (is there a lag at all), aliasing/multimodality check + (is the peak a seasonal artifact), FR/RSS for precision. + 7. **Corroborate with independent witnesses.** Run available Look witnesses (second band, second + line, per-instrument split) and, when coordinated multi-wavelength data exist, the Channel + witness (X-ray-timing cascade). A large split SUBTRACTS confidence (and occasionally hands you + the true lag); a small split does NOT certify correctness. + 8. **Report.** The lag, the confidence tier, the witnesses used, the correction chosen, and what + additional evidence would raise or lower the verdict. + + # CONFIDENCE TIERS (confidence tracks available checks) + + | Tier | Criteria | + |------|----------| + | **refuse** | Below the estimation floor (< ~20 continuum / < ~10 line epochs). Do not estimate. | + | **low** | Estimator runs but verification fails (red-noise-insignificant, unresolved split, aliased, uncorroborated). | + | **medium** | Estimators + internal facts agree, but NO independent witness is available (single band, single line). Right answer, honestly un-witnessed. | + | **high** | As medium, PLUS an independent witness (second band, second line, or Channel cascade) corroborates. | + + Missing data lowers the CEILING, not the estimate. State which witnesses were structurally + unavailable and why. + + # OUTPUT STRUCTURE + + Every verdict MUST include: + 1. **Lag & tier** — the adopted lag (with units and rest/observed frame) and the confidence tier. + 2. **Estimator agreement** — each class's lag, whether they converged, in which mode. + 3. **Arbitration** — if legs disagreed, which fact broke the tie and how. + 4. **Verification** — red-noise significance, aliasing/multimodality, physical-window result. + 5. **Witnesses** — which Look/Channel witnesses ran, their splits, and whether they subtracted + confidence (name any that were unavailable and why). + 6. **Reasoning** — plain-language rationale for the tier and any correction chosen. + 7. **What would change the verdict** — the specific evidence that would raise or lower it. + 8. **Machine-readable summary** — a JSON-serializable object with `lag`, `tier`, `estimators`, + `witnesses`, `flags`, and `reasoning`. + + # ANTI-PATTERNS — NEVER DO THESE + + | # | Anti-Pattern | Correct Behavior | + |---|-------------|------------------| + | 1 | Break an estimator tie with more estimators | Use an independent sampling fact. | + | 2 | Let one witness certify a lag | Witnesses subtract confidence; they never add it. | + | 3 | Trust a witness coupled to mass/luminosity | Require coupling to the lag itself. | + | 4 | Treat convergence as proof | Convergence is necessary, not sufficient; guard agreement-on-wrong. | + | 5 | Refuse on missing fields | Degrade the confidence ceiling instead. | + | 6 | Manufacture a flag on undetectable data | Report the floor honestly. | + | 7 | Report a bare number | Always give lag + tier + witnesses + reasoning. | + | 8 | Use a conditioning fact as an upstream filter | Use it downstream, beside the estimators. | + + # TOOL LIBRARY: reverbtrust_utils (numpy + scipy) + + Call the compute tool with a Python script that imports the bundled library. Key API: + - Estimators: `iccf(t1,f1,t2,f2,lags)`, `centroid(lags,r)`, `frrss(...)`, + `drw_lag(tc,fc,ec,tl,fl,el,lags)`, `javelin_lag(tc,fc,ec,tl,fl,el,lags)`. + - Trust: `physical_window(lags)`, `sampling_acf(t,lags)`, `pair_count(tc,tl,lags)`, + `arbitrate(drw,jav,lags,A,N)`, `rednoise_pvalue(tc,fc,ec,tl,fl,el,lags,sigma,tau)`, + `cross_split(lag_primary,lag_witness)`. + - Channel axis: `cascade(ref_t,ref_f,bands,lags)` (lag + r_max per band behind a driver), + `coupling_strength(driver_t,driver_f,target_t,target_f,lags)` (the X-ray->UV link test: + returns lag, r_max, coupled), `cascade_ordered(results)` (clean wavelength-ordered disk + cascade). Believe a channel lag only where r_max is high; a weak r_max means the link is + disconnected -- report the DISCONNECT, not the lag. + - Verdict: `TrustState(object_id, redshift, n_cont, n_line)` with `.set_lag(...)`, + `.add_verification(...)`, `.add_witness(...)`, `.tier()`, `.verdict()`. + - I/O: `quick_setup()`, `load_lightcurves('/input/lightcurves.json')`, + `save_final_results(results, output_files)`, `quick_finish()`. + + ## REQUIRED SCRIPT TEMPLATE + + ```python + import numpy as np, logging, traceback + from reverbtrust_utils import ( + quick_setup, load_lightcurves, save_final_results, quick_finish, + drw_lag, javelin_lag, sampling_acf, pair_count, arbitrate, + rednoise_pvalue, cross_split, gp_fit, iccf, TrustState, + ) + quick_setup() + results, output_files = {}, {} + try: + d = load_lightcurves('/input/lightcurves.json') + c, l = d['continuum'], d['line'] + lags = np.arange(-5.0, 100.0, 1.0) # physical window applied via floor + # 1) independent estimator classes + drw, _, gp = drw_lag(c['t'], c['f'], c['e'], l['t'], l['f'], l['e'], lags) + jav, lo, hi, multimodal, _ = javelin_lag(c['t'], c['f'], c['e'], l['t'], l['f'], l['e'], lags) + # 2) arbitrate the split with sampling FACTS, not more estimators + A = sampling_acf(c['t'], lags); N = pair_count(c['t'], l['t'], lags) + lag, mode, conflict = arbitrate(drw, jav, lags, A, N) + st = TrustState(d['object_id'], d.get('redshift'), len(c['t']), len(l['t'])) + st.set_lag(lag, estimators={'drw': drw, 'javelin': jav}, + agreement=('conflict' if conflict else mode)) + # 3) verify + _, p_false = rednoise_pvalue(c['t'], c['f'], c['e'], l['t'], l['f'], l['e'], + lags, gp['sigma'], gp['tau']) + st.add_verification('red_noise', passed=(p_false is not None and p_false < 0.1), + detail=f'p_false={p_false}') + st.add_verification('aliasing', passed=not multimodal, detail=f'multimodal={multimodal}') + # 4) independent Look witness (second band) if present + if d.get('continuum2'): + c2 = d['continuum2'] + lag2, *_ = javelin_lag(c2['t'], c2['f'], c2['e'], l['t'], l['f'], l['e'], lags) + split = cross_split(lag, lag2) + st.add_witness('second_band', axis='look', available=True, + corroborates=(split <= TrustState.SPLIT_FLAG_DAYS), + split=split, detail=f'band lag={lag2:.1f}') + else: + st.add_witness('second_band', axis='look', available=False, detail='single band') + # 5) independent Channel witness (X-ray -> UV -> optical cascade) if present + if d.get('cascade'): + from reverbtrust_utils import cascade, coupling_strength, cascade_ordered + cx = d['cascade']; clags = np.arange(-10.0, 10.01, 0.1) + ref = next(b for b in cx['bands'] if b['name'] == cx['reference']) + bands = [(b['name'], b['t'], b['f'], b.get('wavelength')) for b in cx['bands']] + casc = cascade(ref['t'], ref['f'], bands, clags) + ordered = cascade_ordered(casc) + st.add_witness('disk_cascade', axis='channel', available=True, + corroborates=ordered, detail=f'wavelength-ordered={ordered}') + if cx.get('xray_driver'): # grade the X-ray link; the DISCONNECT is real science + xd = cx['xray_driver'] + xlag, xr, coupled = coupling_strength(xd['t'], xd['f'], ref['t'], ref['f'], clags) + st.note(f'xray->UV r_max={xr:.2f} coupled={coupled} ' + f"(lag={xlag:+.2f} d {'believed' if coupled else 'disconnected: lag not reported'})") + results = st.verdict() + except Exception as e: + logging.error(e); traceback.print_exc(); results = {'error': str(e)} + finally: + save_final_results(results, output_files) + quick_finish() + ``` + + {{userGoal}} + {{nodePoolContext}} + {{dataHandlingContext}} + +discoveryExtensions: + humanInTheLoop: Disabled + tools: + - toolId: '{{reverberationTrustToolId}}' + confirmation: Disabled diff --git a/agents/reverberation-trust/metadata.yaml b/agents/reverberation-trust/metadata.yaml new file mode 100644 index 0000000..dc05ec5 --- /dev/null +++ b/agents/reverberation-trust/metadata.yaml @@ -0,0 +1,21 @@ +name: reverberation-trust +type: agent +version: 1.0.0 +associated_tools: + - agents/reverberation-trust/tools/reverberation-trust +publisher: + name: Microsoft Discovery + contact: jhaole@microsoft.com + support_url: https://github.com/microsoft/discovery/discussions/categories/bugs + party: 1p +description: > + Reverberation-mapping trust agent for AGN continuum-to-line time lags. Given a lag and its + light curves, it returns a calibrated trust verdict — the lag, a confidence tier, and the + reasoning — by cross-checking independent estimators, arbitrating their disagreement with + sampling facts (not more estimators), and corroborating with independent Look and Channel + witnesses. Triages a flood of lags: auto-certify the trustworthy, route the rest to humans. +tags: + - reverberation-mapping + - time-series + - trust-triage + - astronomy diff --git a/agents/reverberation-trust/tools/reverberation-trust/Dockerfile b/agents/reverberation-trust/tools/reverberation-trust/Dockerfile new file mode 100644 index 0000000..53cb4c1 --- /dev/null +++ b/agents/reverberation-trust/tools/reverberation-trust/Dockerfile @@ -0,0 +1,24 @@ +FROM python:3.12-slim + +# System updates +RUN apt-get update && apt-get install -y --no-install-recommends \ + ca-certificates && \ + rm -rf /var/lib/apt/lists/* + +# Create directories +RUN mkdir -p /app /output /input + +# Install Python packages (dependency-light: numpy + scipy only) +RUN pip install --no-cache-dir \ + numpy==1.26.4 \ + scipy==1.14.1 + +# Copy utility library +COPY reverbtrust_utils.py /app/reverbtrust_utils.py + +# Set working directory and Python path +WORKDIR /app +ENV PYTHONPATH="/app:${PYTHONPATH}" +ENV PYTHONUNBUFFERED=1 + +CMD ["python3"] diff --git a/agents/reverberation-trust/tools/reverberation-trust/example-input-files/example_lightcurves.json b/agents/reverberation-trust/tools/reverberation-trust/example-input-files/example_lightcurves.json new file mode 100644 index 0000000..17c7e7e --- /dev/null +++ b/agents/reverberation-trust/tools/reverberation-trust/example-input-files/example_lightcurves.json @@ -0,0 +1,15 @@ +{ + "object_id": "RMID-767", + "redshift": 0.52, + "continuum": { + "t": [0.0, 3.1, 6.0, 9.2, 12.0, 15.1, 18.0, 21.0, 24.2, 27.0, 30.1, 33.0, 36.0, 39.1, 42.0, 45.0, 48.2, 51.0, 54.0, 57.1, 60.0, 63.0, 66.1, 69.0, 72.0], + "f": [15.02, 15.18, 15.31, 15.44, 15.39, 15.21, 15.05, 14.92, 14.88, 15.01, 15.20, 15.35, 15.41, 15.33, 15.19, 15.08, 15.14, 15.27, 15.38, 15.30, 15.16, 15.03, 14.97, 15.10, 15.24], + "e": [0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05] + }, + "line": { + "t": [1.0, 4.0, 7.1, 10.0, 13.2, 16.0, 19.1, 22.0, 25.0, 28.1, 31.0, 34.2, 37.0, 40.0, 43.1, 46.0, 49.0, 52.2, 55.0, 58.0, 61.1, 64.0, 67.0, 70.1, 73.0], + "f": [8.41, 8.44, 8.52, 8.60, 8.66, 8.63, 8.55, 8.47, 8.42, 8.40, 8.45, 8.54, 8.62, 8.65, 8.60, 8.52, 8.46, 8.44, 8.49, 8.57, 8.63, 8.61, 8.53, 8.46, 8.43], + "e": [0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04] + }, + "_note": "Illustrative small bundle showing the expected input schema (object_id, redshift, continuum/line time series t/f/e; an optional continuum2 enables the Look witness). Not a science-grade light curve -- real runs pull SDSS-RM / Swift curves via astroquery." +} diff --git a/agents/reverberation-trust/tools/reverberation-trust/reverbtrust_utils.py b/agents/reverberation-trust/tools/reverberation-trust/reverbtrust_utils.py new file mode 100644 index 0000000..e1bd0bb --- /dev/null +++ b/agents/reverberation-trust/tools/reverberation-trust/reverbtrust_utils.py @@ -0,0 +1,576 @@ +""" +reverbtrust_utils -- Reverberation-Lag Trust Agent utility library. + +A dependency-light (numpy + scipy) consolidation of the validated prototype: independent +lag estimators, the sampling-fact arbiter, verification gates, independent witnesses, and a +TrustState confidence tracker -- plus a small /input -> /output I/O convention so a generated +script can run end-to-end in the tool container. + +Design laws (enforced by the API, proven in the prototype): + * Subtract, never add -- a witness can only LOWER confidence, never certify a lag. + * Arbitrate with facts, not estimators -- break a tie with the sampling schedule, not more CCFs. + * Couple to the observable -- a witness helps only if coupled to the LAG itself. + * Degrade, do not refuse -- confidence tracks which checks COULD run, not a minimum-field gate. + +Import pattern: + from reverbtrust_utils import * # everything + from reverbtrust_utils import iccf, drw_lag, TrustState, quick_setup, quick_finish +""" +from __future__ import annotations + +import json +import logging +import os + +import numpy as np +from scipy.linalg import cho_factor, cho_solve +from scipy.optimize import minimize, minimize_scalar + +__all__ = [ + # estimators + "iccf", "centroid", "frrss", + "gp_fit", "gp_predict", "drw_lag", + "javelin_profile", "javelin_lag", + # trust layer + "physical_window", "sampling_acf", "pair_count", "arbitrate", + "simulate_drw", "rednoise_pvalue", "cross_split", + # channel axis + "cascade", "coupling_strength", "cascade_ordered", "COUPLING_THRESHOLD", + "TrustState", "ESTIMATION_FLOOR_CONT", "ESTIMATION_FLOOR_LINE", + # i/o + "quick_setup", "save_final_results", "quick_finish", "load_lightcurves", +] + +# Below these epoch counts the lag is not estimable -> refuse. +ESTIMATION_FLOOR_CONT = 20 +ESTIMATION_FLOOR_LINE = 10 + +# r_max below this => the channel link is disconnected; its lag is meaningless (do not report it). +COUPLING_THRESHOLD = 0.5 + + +# -------------------------------------------------------------------------------------- +# Estimator class 1: interpolated cross-correlation (ICCF; the PyCCF core) +# -------------------------------------------------------------------------------------- +def iccf(t1, f1, t2, f2, lags, min_overlap=20): + """Two-way interpolated CCF: correlation r as a function of trial lag.""" + t1, f1, t2, f2 = map(np.asarray, (t1, f1, t2, f2)) + r = np.full(len(lags), np.nan) + for k, lag in enumerate(lags): + m1 = (t1 + lag >= t2.min()) & (t1 + lag <= t2.max()) + r1 = np.nan + if m1.sum() >= min_overlap: + f2i = np.interp(t1[m1] + lag, t2, f2) + if np.std(f1[m1]) > 0 and np.std(f2i) > 0: + r1 = np.corrcoef(f1[m1], f2i)[0, 1] + m2 = (t2 - lag >= t1.min()) & (t2 - lag <= t1.max()) + r2 = np.nan + if m2.sum() >= min_overlap: + f1i = np.interp(t2[m2] - lag, t1, f1) + if np.std(f1i) > 0 and np.std(f2[m2]) > 0: + r2 = np.corrcoef(f1i, f2[m2])[0, 1] + vals = [v for v in (r1, r2) if not np.isnan(v)] + if vals: + r[k] = np.mean(vals) + return r + + +def centroid(lags, r, frac=0.8): + """Peak lag, centroid (over r >= frac * r_max), and r_max of a CCF. + + Returns (nan, nan, rmax) when there is no usable positive peak -- an all-NaN CCF, or + rmax <= 0 (no positive correlation) -- rather than crashing (nanargmax on all-NaN) or + returning an inconsistent NaN-centroid paired with a finite peak. + """ + lags = np.asarray(lags, float) + r = np.asarray(r, float) + if np.all(np.isnan(r)): + return np.nan, np.nan, np.nan + rmax = np.nanmax(r) + if not np.isfinite(rmax) or rmax <= 0: + return np.nan, np.nan, rmax + peak = lags[np.nanargmax(r)] + sel = (r >= frac * rmax) & ~np.isnan(r) + cent = np.sum(lags[sel] * r[sel]) / np.sum(r[sel]) + return cent, peak, rmax + + +def frrss(t1, f1, e1, t2, f2, e2, lags, n=500, min_overlap=20, rng=None): + """Flux-Randomization / Random-Subset-Selection Monte Carlo -> lag distribution. + + Returns (median, p16, p84, samples) of the CCF centroid -- the FR/RSS uncertainty. + """ + rng = rng or np.random.default_rng(42) + t1, f1, e1 = map(np.asarray, (t1, f1, e1)) + t2, f2, e2 = map(np.asarray, (t2, f2, e2)) + cents = [] + N1, N2 = len(t1), len(t2) + for _ in range(n): + i1 = np.unique(rng.choice(N1, N1, replace=True)) + i2 = np.unique(rng.choice(N2, N2, replace=True)) + ff1 = f1[i1] + rng.normal(0, 1, len(i1)) * e1[i1] + ff2 = f2[i2] + rng.normal(0, 1, len(i2)) * e2[i2] + r = iccf(t1[i1], ff1, t2[i2], ff2, lags, min_overlap=min_overlap) + if np.all(np.isnan(r)): + continue + try: + c, _, rmax = centroid(lags, r) + if rmax > 0.35 and np.isfinite(c): + cents.append(c) + except Exception: + pass + cents = np.array(cents) + if len(cents) < 10: + return np.nan, np.nan, np.nan, cents + return np.median(cents), np.percentile(cents, 16), np.percentile(cents, 84), cents + + +# -------------------------------------------------------------------------------------- +# Estimator class 2: DRW forward-model (JAVELIN-lite; treats slow variability as signal) +# -------------------------------------------------------------------------------------- +def _K(dt, sigma, tau): + return sigma ** 2 * np.exp(-np.abs(dt) / tau) + + +def gp_fit(t, y, yerr): + """Fit a damped-random-walk GP (sigma, tau) by marginal-likelihood maximization.""" + t, y, yerr = map(np.asarray, (t, y, yerr)) + ym = y - y.mean() + DT = t[:, None] - t[None, :] + n = len(t) + + def nll(p): + sigma, tau = np.exp(p) + K = _K(DT, sigma, tau) + np.diag(yerr ** 2) + 1e-8 * np.eye(n) + try: + c, low = cho_factor(K) + except Exception: + return 1e12 + alpha = cho_solve((c, low), ym) + return 0.5 * ym @ alpha + np.sum(np.log(np.diag(c))) + 0.5 * n * np.log(2 * np.pi) + + s0 = np.log(np.std(ym) + 1e-6) + best = None + for tau0 in (20.0, 60.0, 150.0): + r = minimize(nll, x0=[s0, np.log(tau0)], method="Nelder-Mead", + options=dict(maxiter=400, xatol=1e-3, fatol=1e-3)) + if best is None or r.fun < best.fun: + best = r + sigma, tau = np.exp(best.x) + K = _K(DT, sigma, tau) + np.diag(yerr ** 2) + 1e-8 * np.eye(n) + c, low = cho_factor(K) + alpha = cho_solve((c, low), ym) + return dict(sigma=sigma, tau=tau, t=t, ymean=y.mean(), alpha=alpha, chol=(c, low)) + + +def gp_predict(gp, tstar): + tstar = np.asarray(tstar) + Ks = _K(tstar[:, None] - gp["t"][None, :], gp["sigma"], gp["tau"]) + return gp["ymean"] + Ks @ gp["alpha"] + + +def drw_lag(tc, fc, ec, tl, fl, el, lags): + """Best-fit lag by DRW-predict-and-align chi^2 over the trial lags.""" + tc, fc, ec, tl, fl, el = map(np.asarray, (tc, fc, ec, tl, fl, el)) + gp = gp_fit(tc, fc, ec) + w = 1.0 / el ** 2 + chi2 = np.full(len(lags), np.nan) + for k, lag in enumerate(lags): + p = gp_predict(gp, tl - lag) + Sw = w.sum(); Sp = (w * p).sum(); Spp = (w * p * p).sum() + Sy = (w * fl).sum(); Spy = (w * p * fl).sum() + det = Spp * Sw - Sp * Sp + if abs(det) < 1e-12: + continue + A = (Spy * Sw - Sp * Sy) / det + b = (Spp * Sy - Sp * Spy) / det + chi2[k] = (w * (fl - A * p - b) ** 2).sum() + return lags[np.nanargmin(chi2)], chi2, gp + + +# -------------------------------------------------------------------------------------- +# Estimator class 3: JAVELIN-style joint-DRW likelihood (from-scratch facsimile, NOT the JAVELIN package; posterior exposes aliasing) +# -------------------------------------------------------------------------------------- +def _jav_nll(lag, A, tc, dc, ec, tl, dl, el, sigma, tau): + t = np.concatenate([tc, tl - lag]) + amp = np.concatenate([np.ones(len(tc)), A * np.ones(len(tl))]) + d = np.concatenate([dc, dl]) + K = np.outer(amp, amp) * sigma ** 2 * np.exp(-np.abs(t[:, None] - t[None, :]) / tau) + K[np.diag_indices_from(K)] += np.concatenate([ec ** 2, el ** 2]) + 1e-8 + try: + c, low = cho_factor(K) + except Exception: + return 1e12 + alpha = cho_solve((c, low), d) + return 0.5 * d @ alpha + np.sum(np.log(np.diag(c))) + + +def javelin_profile(tc, fc, ec, tl, fl, el, lags, sigma=None, tau=None): + """Return (logL(lag), A(lag)); sigma,tau from a continuum-only DRW fit if not given.""" + tc, fc, ec, tl, fl, el = map(np.asarray, (tc, fc, ec, tl, fl, el)) + if sigma is None or tau is None: + gp = gp_fit(tc, fc, ec); sigma, tau = gp["sigma"], gp["tau"] + dc = fc - fc.mean(); dl = fl - fl.mean() + A0 = (dl.std() + 1e-12) / (dc.std() + 1e-12) + logL = np.full(len(lags), -np.inf); Aopt = np.full(len(lags), np.nan) + for k, lag in enumerate(lags): + r = minimize_scalar(lambda A: _jav_nll(lag, A, tc, dc, ec, tl, dl, el, sigma, tau), + bounds=(0.05 * A0, 20 * A0), method="bounded", + options={"xatol": A0 * 1e-2}) + logL[k] = -r.fun; Aopt[k] = r.x + return logL, Aopt + + +def javelin_lag(tc, fc, ec, tl, fl, el, lags, sigma=None, tau=None, maxc=120): + """Posterior mode, 68% interval, and a multimodality (aliasing) flag.""" + tc, fc, ec = map(np.asarray, (tc, fc, ec)) + if len(tc) > maxc: + idx = np.linspace(0, len(tc) - 1, maxc).astype(int) + tc, fc, ec = tc[idx], fc[idx], ec[idx] + logL, _ = javelin_profile(tc, fc, ec, tl, fl, el, lags, sigma, tau) + post = np.exp(logL - np.nanmax(logL)) + post = post / np.sum(post) + mode = lags[int(np.nanargmax(logL))] + cdf = np.cumsum(post) / np.sum(post) + lo = lags[np.searchsorted(cdf, 0.16)] + hi = lags[min(np.searchsorted(cdf, 0.84), len(lags) - 1)] + peaks = [(lags[i], post[i]) for i in range(1, len(post) - 1) + if post[i] >= post[i - 1] and post[i] >= post[i + 1]] + peaks.sort(key=lambda p: -p[1]) + multimodal = any(p[1] > 0.3 * peaks[0][1] and abs(p[0] - peaks[0][0]) > 15 + for p in peaks[1:]) if peaks else False + return mode, lo, hi, multimodal, logL + + +# -------------------------------------------------------------------------------------- +# Trust layer: physical window, sampling facts + arbiter, red-noise gate, cross-witness +# -------------------------------------------------------------------------------------- +def physical_window(lags, floor=-5.0): + """Boolean mask keeping only physically allowed (>= floor) lags. Line lags are >= ~0.""" + return np.asarray(lags, float) >= floor + + +def sampling_acf(t, lags, binw=1.0): + """Normalized autocorrelation of the observing window (schedule only, no fluxes).""" + t = np.asarray(t, float) + t0 = t.min() + n = int(np.ceil((t.max() - t0) / binw)) + 1 + s = np.zeros(n) + s[np.clip(((t - t0) / binw).astype(int), 0, n - 1)] = 1.0 + s = s - s.mean() + denom = np.sum(s * s) + A = np.full(len(lags), np.nan) + for k, L in enumerate(lags): + sh = int(round(abs(L) / binw)) + if sh >= n or denom <= 0: + continue + a, b = s[sh:], s[:n - sh] + if len(a) > 3: + A[k] = np.sum(a * b) / denom + return A + + +def pair_count(tc, tl, lags, dt=6.0): + """Continuum->line epoch pairs per lag bin = the CCF's actual support at each lag.""" + tc, tl = np.asarray(tc, float), np.asarray(tl, float) + diff = (tl[None, :] - tc[:, None]).ravel() + return np.array([np.sum((diff >= L - dt / 2) & (diff < L + dt / 2)) for L in lags], float) + + +def _at(lags, arr, x): + return arr[int(np.argmin(np.abs(np.asarray(lags) - x)))] + + +def arbitrate(drw, jav, lags, A, N, agree_tol=8.0): + """Break a DRW-vs-JAVELIN split with the sampling facts (not more estimators). + + Returns (lag, mode, conflict): + mode='agree' -> legs agree within agree_tol; lag = mean + mode='fact' -> legs disagree; both facts point the same way; lag = that leg + conflict=True -> facts POINT OPPOSITE ways -> keep the higher-support leg, raise the flag + """ + if abs(drw - jav) <= agree_tol: + return 0.5 * (drw + jav), "agree", False + a_drw, a_jav = _at(lags, A, drw), _at(lags, A, jav) + n_drw, n_jav = _at(lags, N, drw), _at(lags, N, jav) + pick_a = jav if a_drw > a_jav else drw # believe the leg NOT on the bigger echo + pick_n = drw if n_drw > n_jav else jav # believe the leg with more support + if pick_a == pick_n: + return pick_a, "fact", False + return pick_n, "fact", True + + +def simulate_drw(times, sigma, tau, rng): + """One damped-random-walk realization on the given (sorted) times.""" + times = np.asarray(times, float) + n = len(times) + x = np.zeros(n) + x[0] = rng.normal(0, sigma / np.sqrt(2)) + for i in range(1, n): + dt = times[i] - times[i - 1] + rho = np.exp(-dt / tau) + var = max((sigma ** 2 / 2) * (1 - rho ** 2), 0.0) + x[i] = rho * x[i - 1] + rng.normal(0, np.sqrt(var)) + return x + + +def rednoise_pvalue(tc, fc, ec, tl, fl, el, lags, sigma, tau, + min_overlap=20, nsim=300, rng=None): + """False-alarm probability that the CCF peak is a red-noise coincidence (no real lag). + + Simulates independent DRW pairs with matched sampling/noise and NO imposed lag; p_false + is the fraction whose peak r matches or beats the observed r_max. + """ + rng = rng or np.random.default_rng(7) + ec, el = np.asarray(ec), np.asarray(el) + r_obs = np.nanmax(iccf(tc, fc, tl, fl, lags, min_overlap=min_overlap)) + if not np.isfinite(r_obs): + return np.nan, np.nan + count = 0 + for _ in range(nsim): + sc = simulate_drw(tc, sigma, tau, rng) + rng.normal(0, 1, len(tc)) * ec + sl = simulate_drw(tl, sigma, tau, rng) + rng.normal(0, 1, len(tl)) * el # independent + rm = np.nanmax(iccf(tc, sc, tl, sl, lags, min_overlap=min_overlap)) + if np.isfinite(rm) and rm >= r_obs: + count += 1 + return r_obs, count / nsim + + +def cross_split(lag_primary, lag_witness): + """Absolute lag split between a primary look and an independent witness (a Look-axis test). + + A LARGE split subtracts confidence (and occasionally hands you the true lag); a SMALL split + does NOT certify correctness. Directional only -- never used to add confidence. + """ + return abs(float(lag_primary) - float(lag_witness)) + + +# -------------------------------------------------------------------------------------- +# Channel axis: reverberation cascade + coupling-strength (X-ray timing done right) +# -------------------------------------------------------------------------------------- +def _cascade_centroid(lags, r, frac=0.8): + """Centroid over the CONTIGUOUS region around the peak >= frac*r_max (robust to aliasing).""" + lags = np.asarray(lags, float) + r = np.asarray(r, float) + if np.all(np.isnan(r)): + return np.nan, np.nan + ip = int(np.nanargmax(r)); rmax = r[ip] + if not np.isfinite(rmax) or rmax <= 0: + return np.nan, rmax + thr = frac * rmax + lo = ip + while lo - 1 >= 0 and np.isfinite(r[lo - 1]) and r[lo - 1] >= thr: + lo -= 1 + hi = ip + while hi + 1 < len(r) and np.isfinite(r[hi + 1]) and r[hi + 1] >= thr: + hi += 1 + ls, rs = lags[lo:hi + 1], r[lo:hi + 1] + gg = np.isfinite(rs) + return float(np.sum(ls[gg] * rs[gg]) / np.sum(rs[gg])), float(rmax) + + +def coupling_strength(driver_t, driver_f, target_t, target_f, lags, min_overlap=30): + """Channel-axis link test: lag of target behind driver + the correlation STRENGTH r_max. + + Returns (lag, r_max, coupled). `coupled` is False when r_max < COUPLING_THRESHOLD -> the link + is disconnected and its lag is meaningless; report the DISCONNECT, not the lag. (This is how + the framework independently recovers the AGN STORM 2 X-ray/UV disconnect.) `coupled` is None + when overlap is insufficient (all-NaN CCF) -- 'not enough data' is NOT a physical disconnect. + """ + r = iccf(driver_t, driver_f, target_t, target_f, lags, min_overlap=min_overlap) + lag, rmax = _cascade_centroid(lags, r) + if not np.isfinite(rmax): + return lag, rmax, None + coupled = bool(rmax >= COUPLING_THRESHOLD) + return lag, rmax, coupled + + +def cascade(ref_t, ref_f, bands, lags, min_overlap=30): + """Reverberation cascade: lag + r_max of each band behind a reference driver. + + `bands` = ordered iterable of (name, t, f[, wavelength]); redder (longer-wavelength) bands + should lag more in a clean disk cascade. Returns a list of dicts + {name, wavelength, lag, r_max, coupled}. Believe a band's lag only where r_max is high. + """ + out = [] + for entry in bands: + name, tb, fb = entry[0], entry[1], entry[2] + wl = entry[3] if len(entry) > 3 else None + r = iccf(ref_t, ref_f, tb, fb, lags, min_overlap=min_overlap) + lag, rmax = _cascade_centroid(lags, r) + out.append({ + "name": name, "wavelength": wl, "lag": lag, "r_max": rmax, + "coupled": None if not np.isfinite(rmax) else bool(rmax >= COUPLING_THRESHOLD), + }) + return out + + +def cascade_ordered(results, tol=0.15): + """True if the coupled bands' lags increase with wavelength (within tol) -- a clean disk + cascade, the Channel-axis corroboration. Needs >= 3 coupled, wavelength-tagged bands.""" + pts = sorted((r["wavelength"], r["lag"]) for r in results + if r.get("wavelength") and r["coupled"] and np.isfinite(r["lag"])) + lags = [p[1] for p in pts] + if len(lags) < 3: + return False + inversions = sum(1 for i in range(1, len(lags)) if lags[i] < lags[i - 1] - tol) + return inversions == 0 + + +# -------------------------------------------------------------------------------------- +# TrustState -- the confidence tracker that assembles the verdict +# -------------------------------------------------------------------------------------- +class TrustState: + """Accumulates estimator agreement, verification gates, and independent witnesses into + one calibrated verdict. Enforces the tier ladder: refuse < low < medium < high. + + A witness may only LOWER the ceiling. `high` requires an available witness that corroborates. + """ + + SPLIT_FLAG_DAYS = 10.0 # a cross-witness split beyond this subtracts confidence + + def __init__(self, object_id, redshift=None, n_cont=None, n_line=None): + self.object_id = str(object_id) + self.redshift = redshift + self.n_cont = n_cont + self.n_line = n_line + self.lag = None + self.estimators = {} # name -> lag + self.agreement = None # 'agree' | 'fact' | 'conflict' | None + self.verifications = [] # {name, passed, detail} + self.witnesses = [] # {name, axis, available, corroborates, split, detail} + self.notes = [] + + # ---- inputs ------------------------------------------------------------- + def set_lag(self, lag, estimators=None, agreement=None): + self.lag = None if lag is None else float(lag) + if estimators: + self.estimators.update({k: (None if v is None else float(v)) + for k, v in estimators.items()}) + if agreement: + self.agreement = agreement + return self + + def add_verification(self, name, passed, detail=""): + self.verifications.append({"name": name, "passed": bool(passed), "detail": detail}) + return self + + def add_witness(self, name, axis, available, corroborates=None, split=None, detail=""): + """Record an independent witness. axis in {'look','channel','interaction','internal'}. + `available` False means the witness structurally could not run (lowers the ceiling only). + """ + self.witnesses.append({ + "name": name, "axis": axis, "available": bool(available), + "corroborates": corroborates, "split": None if split is None else float(split), + "detail": detail, + }) + return self + + def note(self, text): + self.notes.append(text) + return self + + # ---- derived ------------------------------------------------------------ + def below_floor(self): + if self.n_cont is not None and self.n_cont < ESTIMATION_FLOOR_CONT: + return True + if self.n_line is not None and self.n_line < ESTIMATION_FLOOR_LINE: + return True + return False + + def _verification_failed(self): + return any(not v["passed"] for v in self.verifications) + + def _split_flag(self): + return any(w["split"] is not None and w["split"] > self.SPLIT_FLAG_DAYS + for w in self.witnesses) + + def _has_corroborating_witness(self): + return any(w["available"] and w["corroborates"] is True for w in self.witnesses) + + def tier(self): + """refuse | low | medium | high -- confidence tracks available checks.""" + if self.below_floor() or self.lag is None: + return "refuse" + if self._verification_failed() or self.agreement == "conflict" or self._split_flag(): + return "low" + if self._has_corroborating_witness(): + return "high" + return "medium" + + def verdict(self): + return { + "object_id": self.object_id, + "redshift": self.redshift, + "lag": self.lag, + "tier": self.tier(), + "estimators": self.estimators, + "agreement": self.agreement, + "verifications": self.verifications, + "witnesses": self.witnesses, + "notes": self.notes, + } + + def to_json(self, **kw): + return json.dumps(self.verdict(), default=float, **kw) + + +# -------------------------------------------------------------------------------------- +# I/O convention: /input -> /output (mirrors the tool-container contract) +# -------------------------------------------------------------------------------------- +def quick_setup(input_dir="/input", output_dir="/output", work_dir="/workdir"): + """Ensure I/O dirs exist and configure logging. Returns (input_dir, output_dir, work_dir).""" + for d in (output_dir, work_dir): + os.makedirs(d, exist_ok=True) + logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s") + return input_dir, output_dir, work_dir + + +def load_lightcurves(path): + """Load a light-curve bundle. Expected JSON: + { + "object_id": "RMID-767", "redshift": 0.52, + "continuum": {"t": [...], "f": [...], "e": [...]}, + "line": {"t": [...], "f": [...], "e": [...]}, + "continuum2": {...}, # optional 2nd band -> Look witness + "cascade": { # optional Channel-axis bundle + "reference": "UVW2", + "bands": [{"name": "UVW2", "wavelength": 1928, "t": [...], "f": [...]}, ...], + "xray_driver": {"name": "soft", "t": [...], "f": [...]} + } + } + Returns the parsed dict with numpy arrays for each series. + """ + with open(path) as fh: + d = json.load(fh) + for key in ("continuum", "line", "continuum2", "line2"): + s = d.get(key) + if s: + for c in ("t", "f", "e"): + if c in s: + s[c] = np.asarray(s[c], float) + casc = d.get("cascade") + if casc: + for bnd in casc.get("bands", []): + for c in ("t", "f"): + if c in bnd: + bnd[c] = np.asarray(bnd[c], float) + xd = casc.get("xray_driver") + if xd: + for c in ("t", "f"): + if c in xd: + xd[c] = np.asarray(xd[c], float) + return d + + +def save_final_results(results, output_files=None, output_dir="/output"): + """Write the machine-readable verdict(s) and register any produced files.""" + os.makedirs(output_dir, exist_ok=True) + payload = {"results": results, "output_files": output_files or {}} + with open(os.path.join(output_dir, "results.json"), "w") as fh: + json.dump(payload, fh, indent=2, default=float) + return payload + + +def quick_finish(): + logging.info("reverberation-trust run complete.") diff --git a/agents/reverberation-trust/tools/reverberation-trust/test_reverbtrust_more.py b/agents/reverberation-trust/tools/reverberation-trust/test_reverbtrust_more.py new file mode 100644 index 0000000..4293148 --- /dev/null +++ b/agents/reverberation-trust/tools/reverberation-trust/test_reverbtrust_more.py @@ -0,0 +1,149 @@ +""" +Coverage + regression tests for reverbtrust_utils (pytest-native; t123-migrate 2026-09-02). + +Covers the public surface the base suite left untested (centroid, frrss, javelin_lag, simulate_drw, +cross_split, the I/O convention) and locks in the fixes for three defects the adversarial review +found. One intentional design slack is left as a non-strict xfail -- a different kind of signal +(a documented, benign quirk, not a scheduled fix). + +Pytest-collectable and dual-runnable: `pytest` collects the test_* functions; `python +test_reverbtrust_more.py` re-invokes pytest on this file. +""" +import json +import os +import tempfile + +import numpy as np +import pytest + +import reverbtrust_utils as R +from test_reverbtrust_utils import _make_curves, TAU_TRUE + + +# --------------------------------------------------------------------------- +# Functional / boundary / negative coverage for the untested surface +# --------------------------------------------------------------------------- +def test_centroid_functional_peaks_at_signal(): + lags = np.arange(-10, 11, 1.0) + r = np.exp(-(lags - 3.0) ** 2 / 8.0) # clean Gaussian CCF peaked at +3 + cent, peak, rmax = R.centroid(lags, r) + assert abs(peak - 3.0) < 1e-9 + assert abs(cent - 3.0) < 1.0 + assert 0.99 < rmax <= 1.0 + + +def test_frrss_recovers_injected_lag(): + tc, fc, ec, tl, fl, el = _make_curves() + lags = np.arange(-30, 71, 1.0) + med, lo, hi, samples = R.frrss(tc, fc, ec, tl, fl, el, lags, n=120) + assert lo <= med <= hi + assert abs(med - TAU_TRUE) < 8.0 + assert len(samples) > 10 + + +def test_javelin_lag_recovers_injected_lag(): + tc, fc, ec, tl, fl, el = _make_curves() + lags = np.arange(-30, 71, 1.0) + mode, lo, hi, multimodal, logL = R.javelin_lag(tc, fc, ec, tl, fl, el, lags) + assert abs(mode - TAU_TRUE) < 10.0 + assert lo <= mode <= hi + assert len(logL) == len(lags) + + +def test_simulate_drw_shape_and_boundary(): + rng = np.random.default_rng(0) + x = R.simulate_drw(np.arange(50.0), 0.5, 20.0, rng) + assert len(x) == 50 and np.all(np.isfinite(x)) + assert np.std(x) > 0 # not a flat line + x1 = R.simulate_drw(np.array([0.0]), 0.5, 20.0, rng) # single-epoch boundary + assert len(x1) == 1 and np.isfinite(x1[0]) + + +def test_cross_split_is_symmetric_nonnegative(): + assert R.cross_split(10.0, 25.0) == 15.0 + assert R.cross_split(25.0, 10.0) == 15.0 # symmetry + assert R.cross_split(5.0, 5.0) == 0.0 + + +def test_io_roundtrip_functional(): + with tempfile.TemporaryDirectory() as d: + R.quick_setup(input_dir=d, output_dir=d, work_dir=os.path.join(d, "w")) + bundle = {"object_id": "X1", "redshift": 0.3, + "continuum": {"t": [0, 1, 2], "f": [1, 2, 3], "e": [0.1, 0.1, 0.1]}, + "line": {"t": [1, 2, 3], "f": [1, 2, 3], "e": [0.1, 0.1, 0.1]}} + p = os.path.join(d, "lc.json") + with open(p, "w") as fh: + json.dump(bundle, fh) + loaded = R.load_lightcurves(p) + assert isinstance(loaded["continuum"]["t"], np.ndarray) # arrays parsed + assert loaded["object_id"] == "X1" + R.save_final_results({"tier": "medium"}, output_dir=d) + with open(os.path.join(d, "results.json")) as fh: + back = json.load(fh) + assert back["results"]["tier"] == "medium" + R.quick_finish() + + +def test_load_lightcurves_missing_file_raises(): + # Negative: the external-input boundary should fail cleanly, not silently return junk. + raised = False + try: + R.load_lightcurves(os.path.join(tempfile.gettempdir(), "definitely_not_here_9x.json")) + except (FileNotFoundError, OSError): + raised = True + assert raised + + +def test_hostile_object_id_is_inert_string(): + # Theme 4: the only user-controlled value reaching output is object_id (an identifier + # passthrough, not rendered HTML). Assert it is stored as an inert string, never executed. + hostile = "'; DROP TABLE lags; --" + st = R.TrustState(hostile, n_cont=40, n_line=30).set_lag(5.0, agreement="agree") + v = st.verdict() + assert isinstance(v["object_id"], str) + assert v["object_id"] == hostile # stored verbatim as data, not interpreted + + +# --------------------------------------------------------------------------- +# Regression tests - lock in the fixes for the three defects the review found +# --------------------------------------------------------------------------- +def test_centroid_all_nan_returns_nan_gracefully(): + # FIXED: was ValueError (nanargmax on all-NaN); now returns nan like _cascade_centroid. + lags = np.arange(0, 10, 1.0) + cent, peak, rmax = R.centroid(lags, np.full(len(lags), np.nan)) + assert np.isnan(cent) and np.isnan(peak) and np.isnan(rmax) + + +def test_coupling_insufficient_overlap_returns_none_not_false(): + # FIXED: overlap < min_overlap now reports coupled=None ('cannot judge'), not a false disconnect. + t = np.arange(10.0) + f = np.sin(t) + lag, r, coupled = R.coupling_strength(t, f, t, f, np.arange(-3, 3, 0.5)) + assert coupled is None + + +def test_centroid_no_positive_peak_is_consistent(): + # FIXED: rmax<=0 now yields a coherent no-peak signal (nan cent AND nan peak), not (nan, finite). + lags = np.arange(0, 10, 1.0) + cent, peak, rmax = R.centroid(lags, np.linspace(-0.9, -0.1, len(lags))) + assert rmax <= 0 + assert np.isnan(cent) and np.isnan(peak) + + +# --------------------------------------------------------------------------- +# A different kind of signal: an intentional design slack left as a non-strict xfail +# --------------------------------------------------------------------------- +@pytest.mark.xfail(strict=False, reason=( + "KNOWN/BENIGN: physical_window default floor is -5 d -- it tolerates small negative lags from " + "measurement scatter rather than enforcing a strict 0. Emission-line lags are physically >= 0, " + "so a purist expects negatives rejected; the -5 d slack is intentional and documented. This " + "cannot xpass without changing the default floor, which we deliberately keep.")) +def test_physical_window_rejects_all_negative_lags(): + lags = np.array([-3.0, -1.0, 0.0, 5.0]) + mask = R.physical_window(lags) # default floor -5.0 admits -3 and -1 + assert not mask[0] and not mask[1] + + +if __name__ == "__main__": + import sys + sys.exit(pytest.main([__file__, "-q"])) diff --git a/agents/reverberation-trust/tools/reverberation-trust/test_reverbtrust_utils.py b/agents/reverberation-trust/tools/reverberation-trust/test_reverbtrust_utils.py new file mode 100644 index 0000000..0c08ec5 --- /dev/null +++ b/agents/reverberation-trust/tools/reverberation-trust/test_reverbtrust_utils.py @@ -0,0 +1,204 @@ +""" +Tests for reverbtrust_utils. Runnable with `pytest` or directly: `python test_reverbtrust_utils.py`. + +Covers: estimator lag recovery on an injected lag (ICCF + DRW), the sampling-fact arbiter's +three outcomes, the red-noise gate, and the TrustState tier ladder (the trust contract). +""" +import numpy as np +import pytest + +from reverbtrust_utils import ( + iccf, centroid, drw_lag, arbitrate, sampling_acf, pair_count, + rednoise_pvalue, cross_split, TrustState, physical_window, + cascade, coupling_strength, cascade_ordered, +) + +TAU_TRUE = 20.0 + + +def _make_curves(seed=42, cadence=2.0, tmax=360.0, err=0.02): + """Injected-lag DRW continuum + boxcar-smoothed shifted line, sampled with noise.""" + rng = np.random.default_rng(seed) + td = np.arange(-60, tmax + 60, 0.2) + # dense DRW continuum (amplitude chosen so signal comfortably exceeds sample noise) + x = np.zeros(len(td)); x[0] = rng.normal(0, 0.55) + for i in range(1, len(td)): + rho = np.exp(-(td[i] - td[i - 1]) / 25.0) + x[i] = rho * x[i - 1] + rng.normal(0, np.sqrt(0.30 * (1 - rho ** 2))) + cont_d = 15.0 + x + line_d = np.interp(td - TAU_TRUE, td, cont_d) + line_d = np.convolve(line_d, np.ones(21) / 21, mode="same") + + def samp(fd): + t = np.sort(np.clip(np.arange(0, tmax, cadence) + rng.normal(0, 0.3, int(tmax / cadence)), + 0, tmax)) + f = np.interp(t, td, fd) + e = np.full(len(t), err * np.median(f)) + return t, f + rng.normal(0, 1, len(t)) * e, e + + tc, fc, ec = samp(cont_d) + tl, fl, el = samp(line_d) + return tc, fc, ec, tl, fl, el + + +def test_iccf_recovers_injected_lag(): + tc, fc, ec, tl, fl, el = _make_curves() + lags = np.arange(-30, 71, 1.0) + r = iccf(tc, fc, tl, fl, lags) + cent, peak, rmax = centroid(lags, r) + assert rmax > 0.7 + assert abs(cent - TAU_TRUE) < 6.0 + + +def test_drw_recovers_injected_lag(): + tc, fc, ec, tl, fl, el = _make_curves() + lags = np.arange(-30, 71, 1.0) + lag, chi2, gp = drw_lag(tc, fc, ec, tl, fl, el, lags) + assert abs(lag - TAU_TRUE) < 8.0 + assert gp["tau"] > 0 + + +def test_physical_window_masks_unphysical(): + lags = np.arange(-30, 31, 1.0) + mask = physical_window(lags, floor=-5.0) + assert not mask[0] # -30 excluded + assert mask[lags == 20.0][0] # +20 kept + + +def test_arbitrate_agree_fact_conflict(): + lags = np.arange(0, 80, 1.0) + A = np.zeros_like(lags); N = np.ones_like(lags) + # agree: legs within tolerance + lag, mode, conflict = arbitrate(20.0, 24.0, lags, A, N, agree_tol=8.0) + assert mode == "agree" and not conflict + # fact: legs disagree; both facts favor the long leg (60). + # echo high on the SHORT leg -> distrust 12; support high on 60 -> believe 60. + A2 = A.copy(); A2[lags == 12.0] = 0.9 + N2 = N.copy(); N2[lags == 60.0] = 50.0 + lag, mode, conflict = arbitrate(12.0, 60.0, lags, A2, N2, agree_tol=8.0) + assert mode == "fact" and lag == 60.0 and not conflict + # conflict: facts point opposite ways -> flag raised. + # echo high on the LONG leg -> distrust 60; support high on 60 -> believe 60. + A3 = A.copy(); A3[lags == 60.0] = 0.9 + N3 = N.copy(); N3[lags == 60.0] = 50.0 + lag, mode, conflict = arbitrate(12.0, 60.0, lags, A3, N3, agree_tol=8.0) + assert conflict is True + + +def test_sampling_facts_shapes(): + tc, fc, ec, tl, fl, el = _make_curves() + lags = np.arange(0, 70, 1.0) + A = sampling_acf(tc, lags) + N = pair_count(tc, tl, lags) + assert len(A) == len(lags) and len(N) == len(lags) + assert np.nanmax(N) > 0 + + +def test_rednoise_gate_flags_real_lag_as_significant(): + tc, fc, ec, tl, fl, el = _make_curves() + lags = np.arange(-30, 71, 1.0) + r_obs, p_false = rednoise_pvalue(tc, fc, ec, tl, fl, el, lags, + sigma=0.35, tau=25.0, nsim=100) + assert r_obs > 0.6 + assert p_false < 0.2 # a real injected lag should rarely be faked by red noise + + +def _make_cascade(seed=7): + """A driver DRW + wavelength-ordered lagged bands (coupled) + one disconnected band.""" + rng = np.random.default_rng(seed) + td = np.arange(-40, 200, 0.2) + x = np.zeros(len(td)); x[0] = rng.normal(0, 0.6) + for i in range(1, len(td)): + rho = np.exp(-(td[i] - td[i - 1]) / 20.0) + x[i] = rho * x[i - 1] + rng.normal(0, np.sqrt(0.35 * (1 - rho ** 2))) + drv = 1.0 + x + + def samp(fd, lag): + shifted = np.interp(td - lag, td, fd) + shifted = np.convolve(shifted, np.ones(9) / 9, mode="same") + t = np.arange(0, 150, 1.0) + f = np.interp(t, td, shifted) + return t, f + rng.normal(0, 0.02, len(t)) + + ref_t, ref_f = samp(drv, 0.0) + bands = [("UVW2", ref_t, ref_f, 1928)] + for name, wl, lag in [("UVM2", 2246, 0.4), ("UVW1", 2600, 0.8), + ("U", 3465, 1.4), ("V", 5468, 2.2)]: + tb, fb = samp(drv, lag) + bands.append((name, tb, fb, wl)) + # a disconnected band: an INDEPENDENT DRW (no coupling to the driver) + y = np.zeros(len(td)); y[0] = rng.normal(0, 0.6) + for i in range(1, len(td)): + rho = np.exp(-(td[i] - td[i - 1]) / 20.0) + y[i] = rho * y[i - 1] + rng.normal(0, np.sqrt(0.35 * (1 - rho ** 2))) + dt = np.arange(0, 150, 1.0) + disc_f = np.interp(dt, td, 1.0 + y) + rng.normal(0, 0.02, len(dt)) + return bands, (dt, disc_f) + + +def test_cascade_is_wavelength_ordered(): + bands, _ = _make_cascade() + lags = np.arange(-10.0, 10.01, 0.1) + res = cascade(bands[0][1], bands[0][2], bands, lags) + coupled = [r for r in res if r["coupled"]] + assert len(coupled) >= 4 # the disk bands are all coupled + assert cascade_ordered(res) is True # lags increase with wavelength + + +def test_coupling_strength_flags_disconnect(): + bands, (dt, disc_f) = _make_cascade() + ref_t, ref_f = bands[0][1], bands[0][2] + lags = np.arange(-10.0, 10.01, 0.1) + # a genuinely lagged, coupled band reads as coupled + lag, r, coupled = coupling_strength(ref_t, ref_f, bands[2][1], bands[2][2], lags) + assert coupled is True and r > 0.5 + # an independent band reads as disconnected (the X-ray/UV disconnect pattern) + lag_d, r_d, coupled_d = coupling_strength(ref_t, ref_f, dt, disc_f, lags) + assert coupled_d is False and r_d < 0.5 + + +def test_channel_witness_raises_tier(): + s = TrustState("E", n_cont=40, n_line=30).set_lag( + 3.0, estimators={"drw": 3.0, "javelin": 3.2}, agreement="agree") + s.add_verification("red_noise", passed=True) + s.add_verification("aliasing", passed=True) + s.add_witness("disk_cascade", axis="channel", available=True, corroborates=True, + detail="wavelength-ordered=True") + assert s.tier() == "high" + + +def test_truststate_tier_ladder(): + # refuse: below the estimation floor + s = TrustState("A", n_cont=5, n_line=3).set_lag(20.0) + assert s.tier() == "refuse" + + # medium: estimators agree, no independent witness available + s = TrustState("B", n_cont=40, n_line=30).set_lag( + 20.0, estimators={"drw": 20.0, "javelin": 21.0}, agreement="agree") + s.add_verification("red_noise", passed=True, detail="p_false=0.01") + s.add_witness("i_band", axis="look", available=False, detail="single band") + assert s.tier() == "medium" + + # high: an available witness corroborates + s.add_witness("halpha", axis="look", available=True, corroborates=True, + split=2.0, detail="Ha lag agrees") + assert s.tier() == "high" + + # low: verification fails + s2 = TrustState("C", n_cont=40, n_line=30).set_lag( + 20.0, estimators={"drw": 20.0}, agreement="agree") + s2.add_verification("red_noise", passed=False, detail="p_false=0.4") + assert s2.tier() == "low" + + # low: a large cross-witness split subtracts confidence + s3 = TrustState("D", n_cont=40, n_line=30).set_lag( + 12.0, estimators={"drw": 12.0}, agreement="agree") + s3.add_verification("red_noise", passed=True) + s3.add_witness("i_band", axis="look", available=True, corroborates=False, + split=cross_split(12.0, 60.0), detail="band systematic") + assert s3.tier() == "low" + + +if __name__ == "__main__": + import sys + sys.exit(pytest.main([__file__, "-q"])) diff --git a/agents/reverberation-trust/tools/reverberation-trust/tool.yaml b/agents/reverberation-trust/tools/reverberation-trust/tool.yaml new file mode 100644 index 0000000..4fe1594 --- /dev/null +++ b/agents/reverberation-trust/tools/reverberation-trust/tool.yaml @@ -0,0 +1,40 @@ +name: reverberation-trust +description: > + Reverberation-mapping trust toolkit for AGN continuum-to-line time lags. Provides independent + lag estimators (interpolated CCF with FR/RSS, a DRW forward-model, and a JAVELIN-style joint-DRW + likelihood — a from-scratch facsimile, not the JAVELIN package), a + sampling-fact arbiter that breaks estimator ties with the observing schedule rather than more + estimators, verification gates (physical-lag window, red-noise false-alarm probability, + aliasing/multimodality), independent Look/Channel witnesses, and a TrustState confidence + tracker that assembles a calibrated verdict (refuse / low / medium / high). +version: "1.0.0" +category: astronomy +infra: + - name: worker + infra_type: container + image: + acr: "{name}.azurecr.io/reverbtrust:latest" + compute: + min_resources: + cpu: 1 + ram: 4Gi + storage: 8Gi + gpu: 0 + max_resources: + cpu: 2 + ram: 8Gi + storage: 16Gi + gpu: 0 + infiniband: false + recommended_sku: + - Standard_D4s_v6 + pool_type: static + pool_size: 1 +code_environments: + - language: python + command: python3 "/{{scriptName}}" + description: > + Python code environment with the reverbtrust_utils library (numpy + scipy). Reads a + light-curve bundle from /input, runs estimators + witnesses, and writes the trust + verdict to /output. + infra_node: worker