diff --git a/.claude/skills/ag_build_imaging_model.md b/.claude/skills/ag_build_imaging_model.md new file mode 120000 index 0000000..cb3d83c --- /dev/null +++ b/.claude/skills/ag_build_imaging_model.md @@ -0,0 +1 @@ +../../skills/ag_build_imaging_model.md \ No newline at end of file diff --git a/.claude/skills/ag_configure_search.md b/.claude/skills/ag_configure_search.md new file mode 120000 index 0000000..6c2626e --- /dev/null +++ b/.claude/skills/ag_configure_search.md @@ -0,0 +1 @@ +../../skills/ag_configure_search.md \ No newline at end of file diff --git a/.claude/skills/ag_debug_fit_failure.md b/.claude/skills/ag_debug_fit_failure.md new file mode 120000 index 0000000..3594164 --- /dev/null +++ b/.claude/skills/ag_debug_fit_failure.md @@ -0,0 +1 @@ +../../skills/ag_debug_fit_failure.md \ No newline at end of file diff --git a/.claude/skills/ag_load_results.md b/.claude/skills/ag_load_results.md new file mode 120000 index 0000000..9968c10 --- /dev/null +++ b/.claude/skills/ag_load_results.md @@ -0,0 +1 @@ +../../skills/ag_load_results.md \ No newline at end of file diff --git a/.claude/skills/ag_plot_fit.md b/.claude/skills/ag_plot_fit.md new file mode 120000 index 0000000..05bfee4 --- /dev/null +++ b/.claude/skills/ag_plot_fit.md @@ -0,0 +1 @@ +../../skills/ag_plot_fit.md \ No newline at end of file diff --git a/.claude/skills/ag_prepare_imaging_data.md b/.claude/skills/ag_prepare_imaging_data.md new file mode 120000 index 0000000..1b9a7e8 --- /dev/null +++ b/.claude/skills/ag_prepare_imaging_data.md @@ -0,0 +1 @@ +../../skills/ag_prepare_imaging_data.md \ No newline at end of file diff --git a/.claude/skills/ag_run_search.md b/.claude/skills/ag_run_search.md new file mode 120000 index 0000000..c189e5d --- /dev/null +++ b/.claude/skills/ag_run_search.md @@ -0,0 +1 @@ +../../skills/ag_run_search.md \ No newline at end of file diff --git a/.claude/skills/ag_setup_environment.md b/.claude/skills/ag_setup_environment.md new file mode 120000 index 0000000..d8c85db --- /dev/null +++ b/.claude/skills/ag_setup_environment.md @@ -0,0 +1 @@ +../../skills/ag_setup_environment.md \ No newline at end of file diff --git a/.claude/skills/ag_simulate_dataset.md b/.claude/skills/ag_simulate_dataset.md new file mode 120000 index 0000000..a457972 --- /dev/null +++ b/.claude/skills/ag_simulate_dataset.md @@ -0,0 +1 @@ +../../skills/ag_simulate_dataset.md \ No newline at end of file diff --git a/AGENTS.md b/AGENTS.md index 106bbea..60e6115 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -21,12 +21,11 @@ ask one focused question — never default to the longest possible explanation. ``` Exit 0 = documented API matches the stack. Exit 1 = genuine drift: recommend the pinned version or an audit — [`skills/ag_audit_skill_apis.md`](./skills/ag_audit_skill_apis.md) - owns that procedure. Exit 2/3 = absent/broken stack: report the interpreter, then diagnose - and repair the environment directly - ([`wiki/core/operations/installation.md`](./wiki/core/operations/installation.md) for the - install routes, [`wiki/core/operations/sandbox.md`](./wiki/core/operations/sandbox.md) for - the cache env vars). The `ag_setup_environment` skill that will own the repair itself is not - written yet; `PENDING.md` tracks it. Skip this step by default in maintainer mode. + owns that procedure. Exit 2/3 = absent/broken stack: report the interpreter, then follow + [`skills/ag_setup_environment.md`](./skills/ag_setup_environment.md), which owns the repair + (`wiki/core/operations/installation.md` for the install routes and + `wiki/core/operations/sandbox.md` for the cache env vars are the pages it cites). + Skip this step by default in maintainer mode. ## Safety invariants — default non-negotiable @@ -44,10 +43,11 @@ deliberate refactor). Two are NEVER overridden: the real-data gate and never-rew real data. **If you can't plot it yourself — no code execution, e.g. a GitHub-connector chat — the gate is not waived: ask the user to plot and inspect the data, and to confirm both (a) contaminants and (b) the mask extent, before you compose the fit.** These are the questions - every real-data run must ask, on every harness. Ground the procedure in + every real-data run must ask, on every harness. The procedure itself is owned by + [`skills/ag_prepare_imaging_data.md`](./skills/ag_prepare_imaging_data.md) — read it before + the first real-data fit; it is grounded in `autogalaxy_workspace:scripts/imaging/data_preparation/start_here.py` and the - `imaging/start_here.py` masking section until the `ag_prepare_imaging_data` skill that will own - it is written (see `PENDING.md`). Simulated data is exempt. + `imaging/start_here.py` masking section. Simulated data is exempt. - **Code gate.** A PreToolUse hook validates PyAuto* symbols against the installed library and blocks ones written from memory. If blocked, don't guess — grep `skills/` or introspect `dir()`, then re-run. The hook fires only on harnesses with hook support (Claude Code); @@ -81,10 +81,13 @@ Map every request onto one or more layers: task. Library-API skills are `ag_.md` and produce/evolve a Python script; project-workflow skills (`start-new-project.md`, `contribute-upstream.md`) drive repo-level operations. Skills starting with `_` (`_style.md`, `_bootstrap_skill.md`) are meta-skills — - don't surface them when answering science questions. **No `ag_*` skill for doing galaxy - science is written yet**: the three that exist (`ag_audit_skill_apis`, `ag_update_wiki`, - `ag_refresh_api_docs`) are maintenance workflows for this repo's own content. - `skills/README.md` lists all seven live skills and catalogues the rest by phase with the + don't surface them when answering science questions. The **core modelling loop** is live — + `ag_setup_environment`, `ag_prepare_imaging_data`, `ag_simulate_dataset`, + `ag_build_imaging_model`, `ag_configure_search`, `ag_run_search`, `ag_plot_fit`, + `ag_load_results`, `ag_debug_fit_failure` — and is what a galaxy-science request routes to. + Three further `ag_*` skills (`ag_audit_skill_apis`, `ag_update_wiki`, `ag_refresh_api_docs`) + are maintenance workflows for this repo's own content, not science workflows. + `skills/README.md` lists all sixteen live skills and catalogues the rest by phase with the `autogalaxy_workspace` script that grounds each one. Never activate a skill name you have not confirmed is a file on disk. 3. **Wiki** (`wiki/**/*.md`) — *content*: what a Sersic profile is, which searches exist, @@ -226,8 +229,8 @@ When **not** in maintainer mode, commit at natural checkpoints (a script + its - **Working directories.** Committed scripts → `scripts/`; throwaway plots/data dumps → `scripts/scratch/` (gitignored); `search.fit(...)` output → `./output/`. - **Plot path announcement.** The plot API is **functional-only**: pass - `output_path="scripts/scratch//"`, `output_filename=...`, `output_format="png"` - straight to the `aplt.*` call (e.g. `aplt.subplot_imaging_dataset`, `aplt.subplot_fit_imaging`). + `output_path="scripts/scratch//"` and `output_format="png"` straight to the `aplt.*` + call (e.g. `aplt.subplot_imaging_dataset`, `aplt.subplot_fit_imaging`). **The object-oriented plotters (`aplt.FitImagingPlotter`, `ImagingPlotter`, `GalaxyPlotter`, `GalaxiesPlotter`, `InversionPlotter`, …) and the `aplt.MatPlot2D` / `aplt.Include2D` / `aplt.Output` objects have been removed — do not use them. @@ -235,7 +238,16 @@ When **not** in maintainer mode, commit at natural checkpoints (a script + its connector chat).** Wrong: `aplt.FitImagingPlotter(fit=fit, mat_plot_2d=aplt.MatPlot2D(...)).subplot_fit_imaging()`. Right: - `aplt.subplot_fit_imaging(fit=fit, output_path="scripts/scratch/ring/", output_filename="fit", output_format="png")`. + `aplt.subplot_fit_imaging(fit=fit, output_path="scripts/scratch/ngc1300/", output_format="png")`. + **`output_filename` is not universal** — passing it to a call that does not take it raises + `TypeError`. Only `plot_array`, `plot_grid`, `subplot_imaging_dataset`, + `subplot_imaging_dataset_list`, `subplot_interferometer_dataset` and + `subplot_interferometer_dirty_images` accept it; + `subplot_galaxies` names its file with `auto_filename`; and the remaining fit and galaxy + subplots write a **fixed stem** into `output_path` (`subplot_fit_imaging` → `fit.png`), so the + *directory* is what separates one context from another. Check the signature — + [`skills/ag_plot_fit.md`](./skills/ag_plot_fit.md) and + [`wiki/core/api/plotting.md`](./wiki/core/api/plotting.md) carry the full split. If unsure a PyAuto* symbol exists, ground it against `skills/` or `dir(aplt)` — never write it from memory. Then `print(...)` the absolute path, and after running **quote that absolute path** and offer to open it (platform opener: `open` on macOS, `xdg-open` on Linux, diff --git a/PENDING.md b/PENDING.md index 227f72a..3257b76 100644 --- a/PENDING.md +++ b/PENDING.md @@ -90,29 +90,14 @@ where it is acted on — the run-time and GPU sections of `concepts/non_linear_s ## Phase 4a — the core modelling loop (9 skills) -- [ ] `skills/ag_setup_environment.md` — grounding: `guides/modeling/bug_fix.py`, plus - `wiki/core/operations/installation.md` and `wiki/core/operations/sandbox.md` (both live - since Phase 3) and the RTD installation pages they were written from. `AGENTS.md`'s - session-start step and `skills/ag_audit_skill_apis.md` route environment failures here; - until it exists they route to `--check-install` and those two pages. -- [ ] `skills/ag_prepare_imaging_data.md` — grounding: - `imaging/data_preparation/start_here.py`, `imaging/data_preparation/examples/`, - `imaging/data_preparation/gui/`. **This skill owns the real-data inspection gate's - procedure**, which `AGENTS.md` currently points at the workspace scripts for. -- [ ] `skills/ag_simulate_dataset.md` — grounding: `imaging/simulator.py`, - `imaging/simulator_sersic.py`, `imaging/simulator_sample.py`. -- [ ] `skills/ag_build_imaging_model.md` — grounding: `imaging/start_here.py`, - `imaging/modeling.py`, `guides/modeling/cookbook.py`. -- [ ] `skills/ag_configure_search.md` — grounding: `guides/modeling/searches.py`, - `guides/modeling/customize.py`. -- [ ] `skills/ag_run_search.md` — grounding: `imaging/modeling.py`, - `guides/modeling/bug_fix.py`. -- [ ] `skills/ag_plot_fit.md` — grounding: `imaging/plot.py`, `guides/plot/start_here.py`, - `guides/plot/plotters.py`. -- [ ] `skills/ag_load_results.md` — grounding: `guides/results/start_here.py`, - `guides/results/aggregator/`, `guides/results/latent_variables.py`. -- [ ] `skills/ag_debug_fit_failure.md` — grounding: `guides/modeling/bug_fix.py`, - HowToGalaxy `chapter_2_modeling/tutorial_4_dealing_with_failure`. +**Delivered.** All nine core-loop skills are on disk with `.claude/skills/` symlinks, rows in +`wiki/core/external/skill_citation_map.md` and entries in the `skills/README.md` Index: +`ag_setup_environment`, `ag_prepare_imaging_data`, `ag_simulate_dataset`, +`ag_build_imaging_model`, `ag_configure_search`, `ag_run_search`, `ag_plot_fit`, +`ag_load_results`, `ag_debug_fit_failure`. Two `AGENTS.md` hand-offs closed with them: the +session-start environment repair now routes to `ag_setup_environment`, and the real-data +inspection gate's procedure is owned by `ag_prepare_imaging_data` rather than pointing at the +workspace scripts directly. ## Phase 4b — features beyond a single smooth profile (8 skills) diff --git a/llms.txt b/llms.txt index b3ec760..797ed0e 100644 --- a/llms.txt +++ b/llms.txt @@ -18,9 +18,9 @@ from memory. Read order: [AGENTS.md](./AGENTS.md) → [skills/README.md](./skills/README.md) → the relevant skill → [wiki/](./wiki). **In chat, lead by engaging.** Ask what the user is trying to do, ask them to describe or plot their data (you can't inspect files on their machine), plan the model, and draft the scripts — do this *first*, and never open a reply with what you can't do. The assistant can also route to examples, explain galaxy morphology and surface-brightness fitting, and review pasted scripts / errors / figures. It cannot run code or guarantee drafts against the installed PyAutoGalaxy version, so switching to a local coding agent (Claude Code, Codex, or Codex web with the repo connected) is the step for **actually running the fit or iterating live** — raise that handoff when execution becomes the blocker, not as an opening disclaimer. -**⚠️ API currency — read before writing any PyAutoGalaxy code.** Older PyAutoGalaxy releases are heavily represented in model training data and their API is **out of date**. Do **not** write PyAutoGalaxy from memory. The single most common stale error: **plotting is functional now.** Use `aplt.subplot_fit_imaging(fit=fit, output_path=..., output_filename=..., output_format="png")` and `aplt.subplot_imaging_dataset(dataset=..., ...)`. The old object-oriented plotters — `aplt.FitImagingPlotter`, `aplt.ImagingPlotter`, `aplt.GalaxyPlotter`, `aplt.GalaxiesPlotter`, `aplt.InversionPlotter` — and the `aplt.MatPlot2D` / `aplt.Include2D` / `aplt.Output` configuration objects have all been **removed**; never emit them. The whole plotting surface is module-level functions, and `dir(autogalaxy.plot)` is the authoritative list. If you can't point at a live example or a `dir()` listing for a call, say it's unverified rather than guessing. +**⚠️ API currency — read before writing any PyAutoGalaxy code.** Older PyAutoGalaxy releases are heavily represented in model training data and their API is **out of date**. Do **not** write PyAutoGalaxy from memory. The single most common stale error: **plotting is functional now.** Use `aplt.subplot_fit_imaging(fit=fit, output_path=..., output_format="png")` and `aplt.subplot_imaging_dataset(dataset=..., output_path=..., output_filename=..., output_format="png")` — note that `output_filename` is accepted by the dataset subplots, `plot_array` and `plot_grid`, but **not** by `subplot_fit_imaging`, which writes a fixed `fit.png` into `output_path` (passing it raises `TypeError`). The old object-oriented plotters — `aplt.FitImagingPlotter`, `aplt.ImagingPlotter`, `aplt.GalaxyPlotter`, `aplt.GalaxiesPlotter`, `aplt.InversionPlotter` — and the `aplt.MatPlot2D` / `aplt.Include2D` / `aplt.Output` configuration objects have all been **removed**; never emit them. The whole plotting surface is module-level functions, and `dir(autogalaxy.plot)` is the authoritative list. If you can't point at a live example or a `dir()` listing for a call, say it's unverified rather than guessing. -**⚠️ This assistant is young.** Four skills exist (two meta, two project-workflow); the `ag_*` task skills and most of `wiki/core/` are still being written, phase by phase. [`skills/README.md`](./skills/README.md) and [`PENDING.md`](./PENDING.md) say exactly what is present and what is planned, and every planned item names the `autogalaxy_workspace` script that grounds it. **Answer from those grounding scripts and say that you did** — never describe a pending skill or wiki page as though you had read it. +**⚠️ This assistant is young.** Sixteen skills exist: the nine-skill core modelling loop (environment, data preparation, simulation, model building, search configuration, running the fit, plotting, results, debugging), two meta, two project-workflow and three maintenance. The Phase-4b **feature** skills (bases and MGE, pixelisations, ellipse fitting, multi-dataset, interferometry, multi-galaxy, search chaining), the literature wiki and the benchmark suite are still being written, phase by phase. [`skills/README.md`](./skills/README.md) and [`PENDING.md`](./PENDING.md) say exactly what is present and what is planned, and every planned item names the `autogalaxy_workspace` script that grounds it. **Answer from those grounding scripts and say that you did** — never describe a pending skill or wiki page as though you had read it. ## Agent instructions (canonical) @@ -32,7 +32,7 @@ Read order: [AGENTS.md](./AGENTS.md) → [skills/README.md](./skills/README.md) ## Reference -- [wiki/core/index.md](./wiki/core/index.md): curated PyAuto\* reference. The `stack/` pages (one per library) are live; concepts, API catalogues, operations and external routing arrive in a later phase, and the index itself marks the gap. +- [wiki/core/index.md](./wiki/core/index.md): curated PyAuto\* reference. The `stack/`, `api/`, `concepts/`, `operations/` and `external/` sections are live; the index itself names the pages still missing (dataset layout, HPC). - [PENDING.md](./PENDING.md): the authoritative ledger of every unwritten skill and wiki page, with its grounding script and target phase. ## Runnable examples & tutorials (elsewhere) diff --git a/skills/README.md b/skills/README.md index 24b1694..9eb619c 100644 --- a/skills/README.md +++ b/skills/README.md @@ -27,16 +27,62 @@ configured) via symlinks; the canonical files live here. ## Index — what exists today -Seven skills are written: two meta, two project-workflow, three maintenance. **Everything -else in this file is a plan, not a file** — the "Pending" section below is a catalogue of what -has not been authored yet, and deliberately does not link to anything. Every entry here that -is a link resolves; if you find a link that doesn't, that is a bug worth fixing rather than a -file worth waiting for. +Sixteen skills are written: nine for the core modelling loop, two meta, two project-workflow, +three maintenance. **Everything else in this file is a plan, not a file** — the "Pending" +section below is a catalogue of what has not been authored yet, and deliberately does not link +to anything. Every entry here that is a link resolves; if you find a link that doesn't, that is +a bug worth fixing rather than a file worth waiting for. + +The core modelling loop is live; the feature skills (bases, pixelisations, ellipse fitting, +multi-dataset, interferometry, multi-galaxy, search chaining) arrive in Phase 4b. For a request +that falls outside the nine below, answer from the installed source and the grounding scripts +named in the Pending tables, say that is what you did, and offer to author the skill via +[`_bootstrap_skill.md`](./_bootstrap_skill.md). -No skill for *doing* galaxy science exists yet: the modelling loop arrives in Phase 4a and the -feature skills in Phase 4b. Until then, answer a science request from the installed source and -the grounding scripts named in the Pending tables, say that is what you did, and offer to -author the skill via [`_bootstrap_skill.md`](./_bootstrap_skill.md). +### Galaxy modelling — the core loop + +Read in this order for an end-to-end fit; each is usable on its own. Every one is +**python-first**: the deliverable is a runnable script plus the understanding to evolve it. + +- [`ag_setup_environment.md`](./ag_setup_environment.md) — install, diagnose and repair the + PyAutoGalaxy environment so galaxy-modelling code can actually run: a fresh pip install, the + JAX and numba extras, writable caches for a sandbox, `activate.sh` interpreter resolution, + a shared/cluster checkout and the Colab entry point, ending in a saved verification script. +- [`ag_prepare_imaging_data.md`](./ag_prepare_imaging_data.md) — load a user's own CCD imaging + from FITS into an `ag.Imaging` and get it ready to fit: pixel scale, flux units, RMS + noise-map, PSF, mask extent, contaminants, over-sampling and the `info.json` sidecars. + **Owns the real-data inspection gate.** +- [`ag_simulate_dataset.md`](./ag_simulate_dataset.md) — simulate imaging (or visibilities) of a + galaxy with known truth: grid and over-sampling, PSF, exposure time and background sky, light + profiles, FITS output plus a `galaxies.json` truth record, S/N targeting, whole samples in a + loop, and the `should_simulate` convention. +- [`ag_build_imaging_model.md`](./ag_build_imaging_model.md) — compose the model for an imaging + fit: an `af.Model` / `af.Collection` tree of light profiles on one or more galaxies plus the + `ag.AnalysisImaging` that scores it — single Sersic, bulge-plus-disk, linear profiles, MGE, + prior customisation, pairing and assertions, `ag.DatasetModel`, and checking with + `print(model.info)` before spending a search. +- [`ag_configure_search.md`](./ag_configure_search.md) — choose and configure the non-linear + search: `af.Nautilus` for a quotable posterior, `af.MultiStartProdigy` for a fast MAP check, + `af.DynestyStatic` for ellipse fitting and cross-checks, plus `n_live` / `n_batch` / + `n_starts`, the output cadence, start-point initialisers, grid searches, and the + `unique_tag` resume semantics that silently reuse a fit across datasets. +- [`ag_run_search.md`](./ag_run_search.md) — drive `search.fit(model=model, analysis=analysis)` + to completion and read what it wrote: the output-folder anatomy and on-the-fly announcement, + resume behaviour, quick updates, JAX/GPU acceleration and VRAM checks, the + `PYAUTO_TEST_MODE` smoke levels, and the `if __name__ == "__main__"` parallelisation fix. +- [`ag_plot_fit.md`](./ag_plot_fit.md) — visualise a dataset, a galaxy or a fit with the + functional `aplt` API: dataset and fit subplots, individual residual / normalised-residual / + chi-squared panels, per-galaxy breakdowns, log10 stretch and fixed colour limits, overlays, + figures and FITS to disk — plus the residual-inspection discipline for judging a fit. +- [`ag_load_results.md`](./ag_load_results.md) — get a completed fit back into Python and turn + it into science: the in-session `Result`, direct `from_json` / `SamplesNest.from_table` + loading of one output folder, the `Samples` API for medians and errors, and the aggregator + for a whole sample with its `ag.agg` generators, queries and CSV/FITS/PNG exports. +- [`ag_debug_fit_failure.md`](./ag_debug_fit_failure.md) — triage a fit that crashed, stalled or + finished with parameters you do not believe, through a failure taxonomy (environment, data, + model, priors, search settings, stale result) and the probes that separate them — including + the two silent failures: a resumed fit whose identifier ignored the data, and a cached result + mistaken for a new one. ### Meta @@ -91,20 +137,6 @@ and shrinks as each phase lands. named grounding scripts, say that is what you did, and offer to author the skill via [`_bootstrap_skill.md`](./_bootstrap_skill.md). -### Phase 4a — the core modelling loop - -| Skill | Purpose | Grounding (`autogalaxy_workspace/scripts/`) | -|-------|---------|--------------------------------------------| -| `ag_setup_environment` | detect absent or broken PyAuto\* environments, install via pip or editable clones, configure caches, verify imports | `guides/modeling/bug_fix.py` + the RTD installation pages | -| `ag_prepare_imaging_data` | load and preprocess FITS imaging, decide masking for real data, measure noise, prepare the PSF | `imaging/data_preparation/start_here.py`, `imaging/data_preparation/examples/`, `imaging/data_preparation/gui/` | -| `ag_simulate_dataset` | synthesise a galaxy dataset from a ground-truth light model, including population samples | `imaging/simulator.py`, `imaging/simulator_sersic.py`, `imaging/simulator_sample.py` | -| `ag_build_imaging_model` | compose a galaxy's light model (bulge, disk, single Sersic) and wrap it in an imaging analysis | `imaging/start_here.py`, `imaging/modeling.py`, `guides/modeling/cookbook.py` | -| `ag_configure_search` | pick and tune a non-linear search or gradient optimizer for the problem at hand | `guides/modeling/searches.py`, `guides/modeling/customize.py` | -| `ag_run_search` | execute `search.fit(model=..., analysis=...)`, monitor convergence, read the live output folder | `imaging/modeling.py`, `guides/modeling/bug_fix.py` | -| `ag_plot_fit` | plot the model image, residuals, normalised residuals and chi-squared map | `imaging/plot.py`, `guides/plot/start_here.py`, `guides/plot/plotters.py` | -| `ag_load_results` | load a completed fit's galaxies, samples, dataset and FITS products from its output folder | `guides/results/start_here.py`, `guides/results/aggregator/` | -| `ag_debug_fit_failure` | diagnose a fit that didn't converge or produced unphysical structural parameters | `guides/modeling/bug_fix.py`, HowToGalaxy `chapter_2_modeling/tutorial_4_dealing_with_failure` | - ### Phase 4b — features beyond a single smooth profile | Skill | Purpose | Grounding (`autogalaxy_workspace/scripts/`) | diff --git a/skills/_style.md b/skills/_style.md index d65909c..a2ba147 100644 --- a/skills/_style.md +++ b/skills/_style.md @@ -379,17 +379,30 @@ back up the skills. ## Plot output and path announcement Skills that produce visualisations save them through the function-style plotting API -imported as `aplt` — every entry point takes `output_path` / `output_filename` -/ `output_format` kwargs directly. Three rules: +imported as `aplt` — every entry point takes `output_path` / `output_format` kwargs +directly, and *some* also take `output_filename`. Three rules: -1. **Pass `output_path` / `output_filename` / `output_format` directly to - each plot function.** Every `plot_*` and `subplot_*` helper in `aplt` accepts - these kwargs, e.g. `aplt.subplot_imaging_dataset(dataset=…, +1. **Pass `output_path` / `output_format` directly to each plot function**, plus + `output_filename` on the functions that accept it, e.g. + `aplt.subplot_imaging_dataset(dataset=…, output_path="scripts/scratch//", output_filename=…, output_format="png")`. Never rely on interactive display — the user is often running the script from a terminal where `plt.show()` flashes and vanishes. The `` slug is usually the dataset name; for general exploration any short slug works. + + **`output_filename` is not universal, and passing it where it does not exist + raises `TypeError`.** It is accepted by `plot_array`, `plot_grid`, + `subplot_imaging_dataset`, `subplot_imaging_dataset_list`, + `subplot_interferometer_dataset` and `subplot_interferometer_dirty_images`. + `subplot_galaxies` names its file with + `auto_filename` instead. Every other subplot writes a **fixed stem** into + `output_path` — `subplot_fit_imaging` → `fit.png`, + `subplot_fit_imaging_of_galaxy` → `of_galaxy_.png`, + `subplot_galaxy_images` → `galaxy_images.png` — so give each fit or each + variant of a figure its own directory rather than trying to rename the file. + Check the signature; [`../wiki/core/api/plotting.md`](../wiki/core/api/plotting.md) + tabulates the whole split. 2. **`print(...)` each plot's path** at the end of the Python recipe so the absolute location lands in stdout. Use `print(f"Saved to: {PLOT_DIR.resolve()}")` once per branch (sufficient diff --git a/skills/ag_build_imaging_model.md b/skills/ag_build_imaging_model.md new file mode 100644 index 0000000..c53d5b5 --- /dev/null +++ b/skills/ag_build_imaging_model.md @@ -0,0 +1,474 @@ +--- +name: ag_build_imaging_model +description: Compose the galaxy model for a CCD imaging fit — an `af.Model` / `af.Collection` tree of light profiles on one or more galaxies, plus the `ag.AnalysisImaging` that scores it against a loaded dataset. Covers a single Sersic, a bulge-plus-disk decomposition, why `ag.lp_linear` profiles are the default, a Multi-Gaussian Expansion for irregular morphology, prior customisation, pairing / fixing / assertions, `ag.DatasetModel` for a residual sky or astrometric offset, modelling contaminating extra galaxies with fixed centres, and checking the model with `print(model.info)` and a single likelihood evaluation before spending a search on it. Use once a dataset is loaded and masked. Requires the real-data inspection gate to have been satisfied first (`ag_prepare_imaging_data`). Not for choosing or configuring the search (`ag_configure_search`), not for running the fit, and not for visibility-plane data. +--- + +# Composing an imaging galaxy model + +The model is your statement of what this galaxy *could* be: how many components, with what +functional forms, which parameters free and which fixed, and what prior on each. It is the +half of the fit you are responsible for — the search only explores what you wrote down. A +model too simple leaves structure in the residuals and biases every parameter that has to +absorb it; a model too complex has a parameter space the search cannot map, and returns wide +errors or a local maximum. Getting this right is most of galaxy morphology. + +A fit has exactly two halves: the **model** says what could be true, and the **analysis** +holds a dataset and knows how to score a proposal against it. A search shuttles between them: + +```python +result = search.fit(model=model, analysis=analysis) +``` + +Everything in this skill builds one of those two arguments. The canonical script is +`autogalaxy_workspace:scripts/imaging/modeling.py`; the API page is +[`../wiki/core/api/analysis_objects.md`](../wiki/core/api/analysis_objects.md), and the +systematic composition reference is +`autogalaxy_workspace:scripts/guides/modeling/cookbook.py`. + +## Before you compose — the gate + +If the dataset is **real observational data**, confirm the inspection gate has been +satisfied: that the user has seen `dataset.png`, and that both (a) contaminating extra +galaxies / foreground stars / artefacts and (b) the mask extent were settled from that look. +If either is still open, stop and do it in +[`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md) first — composing a model over data +with an unaddressed contaminant, or a mask left at a silent default, produces a confident +wrong answer rather than an obvious failure. Simulated data is exempt. + +## Ask + +- *"One galaxy, or several whose light blends?"* The only structural choice. This skill + assumes a single target with optional contaminating neighbours; two or more co-dominant + galaxies is a different skill (`ag_multi_galaxy_and_cluster`, not yet written — see + [`../PENDING.md`](../PENDING.md)). +- *"A single Sersic to start, or straight to a bulge-plus-disk decomposition?"* Single Sersic + first is almost always the right answer, even when you want the decomposition: it converges + fast, tells you the galaxy's overall size and concentration, and gives you something to + compare the two-component fit against. +- *"Is the morphology smooth, or irregular — spiral arms, asymmetry, a merger?"* Smooth means + Sersic-family profiles. Irregular means a Multi-Gaussian Expansion, and beyond that a + pixelised reconstruction (`ag_pixelization`, not yet written). +- *"What redshift?"* Even a ballpark value is enough; it is needed only for physical-unit + conversions, not for the fit itself. +- *"Has the sky background been subtracted, and do you trust that subtraction?"* If the faint + outer envelope is what you care about, the answer decides whether you need a + `ag.DatasetModel`. + +## Branch — a single Sersic + +The minimum viable model, and the right first fit for almost any galaxy. + +```python +""" +Galaxy Structure: Single Sersic +=============================== + +Compose and check a one-component galaxy model for CCD imaging: an elliptical Sersic profile +whose effective radius and Sersic index measure the galaxy's size and central concentration. +This is the first fit to run on any galaxy — it converges quickly, and its inferred size and +concentration are the reference every more complex model is judged against. + +__Contents__ + +- **Imports:** Import the required libraries. +- **Dataset:** Load the masked, over-sampled imaging dataset produced during preparation. +- **Model:** Compose the galaxy's light as a single linear Sersic profile. +- **Check:** Print the model and evaluate the likelihood once before committing to a search. +""" + +""" +__Imports__ +""" +from pathlib import Path + +import autofit as af +import autogalaxy as ag + +DATASET_PATH = Path("dataset") / "imaging" / "my_galaxy" +PIXEL_SCALES = 0.1 +MASK_RADIUS = 2.5 + +""" +__Dataset__ + +The dataset arrives already prepared: loaded from FITS, contaminant noise scaled, masked to +the radius chosen by inspecting the data, and over-sampled adaptively in the centre where the +profile's intensity gradient is steep. All four decisions belong to data preparation rather +than to the model, and all four change the answer, so they are reproduced here rather than +hidden (`PyAutoArray:autoarray/dataset/imaging/dataset.py`). +""" +dataset = ag.Imaging.from_fits( + data_path=DATASET_PATH / "data.fits", + noise_map_path=DATASET_PATH / "noise_map.fits", + psf_path=DATASET_PATH / "psf.fits", + pixel_scales=PIXEL_SCALES, +) + +mask = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=dataset.pixel_scales, + radius=MASK_RADIUS, +) +dataset = dataset.apply_mask(mask=mask) + +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=dataset.grid, + sub_size_list=[8, 4, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0)], +) +dataset = dataset.apply_over_sampling(over_sample_size_lp=over_sample_size) + +""" +__Model__ + +A single elliptical Sersic. Its `sersic_index` is free, so the fit itself decides whether this +galaxy is a de Vaucouleurs-like spheroid (n ~ 4) or an exponential disk (n ~ 1) rather than +you asserting it, and `effective_radius` is the radius containing half the profile's light. +The profile is elliptical via `ell_comps` rather than an axis-ratio and position-angle pair, +which avoids the periodic-boundary pathology a position angle creates for a non-linear search +(`PyAutoGalaxy:autogalaxy/profiles/light/standard/sersic.py`). + +We use `ag.lp_linear` rather than `ag.lp`: a linear light profile has its `intensity` solved +analytically by a linear inversion inside each likelihood evaluation instead of being sampled +as a free parameter. The profile is identical; one dimension simply leaves the search's +parameter space, which improves speed, accuracy and reliability at no cost. Every workspace +modeling example uses them by default +(`PyAutoGalaxy:autogalaxy/profiles/light/linear/sersic.py`). + +`af.Model` wraps the class so its parameters become free; `af.Collection` groups models into +the named tree the analysis expects — the `galaxies` collection is a convention the analysis +relies on, not a stylistic choice (`PyAutoFit:autofit/mapper/prior_model/collection.py`). +""" +bulge = af.Model(ag.lp_linear.Sersic) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +""" +__Check__ + +`model.info` prints every free parameter with its prior, and it is the cheapest sanity check +available: a model subtly larger, smaller or more constrained than intended shows up here +rather than after hours of searching. Evaluating the likelihood once at the prior medians is +the second cheap check — it proves the model, dataset and analysis are mutually compatible +and that a finite likelihood comes back, without paying for inference. + +`ag.AnalysisImaging` holds the dataset and defines the `log_likelihood_function` the search +calls. `use_jax=True` is the default and is what makes the gradient-based searches possible +at all; pass `use_jax=False` when debugging, because NumPy tracebacks are far easier to read +than JAX ones (`PyAutoGalaxy:autogalaxy/imaging/model/analysis.py`). +""" +print(model.info) +print(f"Total free parameters = {model.total_free_parameters}") + +analysis = ag.AnalysisImaging(dataset=dataset, use_jax=True) + +log_likelihood = analysis.log_likelihood_function( + instance=model.instance_from_prior_medians() +) + +print(f"log likelihood at prior medians: {float(log_likelihood):.2f}") +``` + +That composes to **N = 6** free parameters: centre (2), `ell_comps` (2), `effective_radius`, +`sersic_index` — with `intensity` solved rather than sampled. Six parameters is a parameter +space any search handles comfortably. + +What each profile represents and when to reach for which is +[`../wiki/core/concepts/light_profiles.md`](../wiki/core/concepts/light_profiles.md); the +`ell_comps` convention has its own section there. The catalogue of every available profile, +by module, is +[`../wiki/core/api/light_profile_catalog.md`](../wiki/core/api/light_profile_catalog.md). + +## Branch — bulge plus disk + +The decomposition galaxy morphology usually wants, because the ratio of the two components' +inferred luminosities is the bulge-to-total light ratio. Adapted from +`autogalaxy_workspace:scripts/imaging/modeling.py`. + +```python +bulge = af.Model(ag.lp_linear.Sersic) +disk = af.Model(ag.lp_linear.Exponential) + +bulge.centre = disk.centre + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, disk=disk) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(model.info) +``` + +N = 9. The physics is in the asymmetry between the two components: the bulge's +`sersic_index` is free so it can find a concentrated spheroid, while the `Exponential` disk +has n fixed at 1 by construction. That contrast is what makes the decomposition identifiable +at all — two free Sersics on the same galaxy are largely degenerate, and the fit will happily +trade light between them. + +`bulge.centre = disk.centre` is the load-bearing line. It asserts the two components are +concentric, which is physically reasonable for most galaxies and removes two parameters +rather than fitting the same centre twice. Drop it only if you have a specific reason to +believe the components are offset — a lopsided or interacting system — and expect a harder +search when you do. + +The attribute names are entirely yours. `bulge`, `disk`, `bar`, `clump`, `bulge_0` — whatever +you pass to `af.Model(ag.Galaxy, ...)` is the key you address afterwards +(`model.galaxies.galaxy.bulge.sersic_index`), because `ag.Galaxy` takes `redshift` and then +arbitrary keyword arguments (`PyAutoGalaxy:autogalaxy/galaxy/galaxy.py`). Name them for what +they measure. + +There is a concise form when you do not need to touch the components first — passing a +profile *class* rather than an `af.Model` promotes it automatically +(`autogalaxy_workspace:scripts/guides/modeling/cookbook.py`): + +```python +galaxy = af.Model( + ag.Galaxy, + redshift=0.5, + bulge=ag.lp_linear.Sersic, + disk=ag.lp_linear.Exponential, + bar=ag.lp_linear.Sersic, +) +``` + +How several profiles and galaxies compose into what is actually fitted is +[`../wiki/core/concepts/galaxies.md`](../wiki/core/concepts/galaxies.md). + +## Branch — a Multi-Gaussian Expansion, for morphology a Sersic cannot fit + +When the galaxy is irregular, asymmetric or has structure a smooth profile leaves in the +residuals, decompose its light into tens of Gaussians instead. Because they are *linear* +profiles, all of their intensities are solved by the inversion, so a 20-Gaussian basis costs +the search almost nothing in dimensionality — it is far more flexible than a Sersic at a +*lower* parameter count. + +The composition is long and technical, so the library ships a helper: + +```python +bulge = ag.model_util.mge_model_from( + mask_radius=MASK_RADIUS, + total_gaussians=20, + centre_prior_is_uniform=True, +) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(model.info) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/start_here.py`, where the MGE is the +*recommended default* model — it balances speed, flexibility and accuracy well enough to fit +the vast majority of galaxies. N = 4: a shared centre and shared `ell_comps` across the +basis, with every `sigma` fixed on a logarithmic ladder from below the pixel scale to beyond +the galaxy's extent, and every intensity solved. `mask_radius` is what sets the outer end of +that ladder, so pass the radius you actually masked at. + +Other arguments worth knowing: `centre_fixed=(y, x)` pins the centre entirely (used for +contaminating neighbours, below), `gaussian_per_basis` splits the ladder into several bases +with independent ellipticities for a twisting galaxy, and `use_spherical=True` drops the +ellipticity. The physics of why fewer free parameters wins here is +[`../wiki/core/concepts/linear_light_profiles_and_mge.md`](../wiki/core/concepts/linear_light_profiles_and_mge.md); +the orthonormal-basis alternative is +[`../wiki/core/concepts/shapelets.md`](../wiki/core/concepts/shapelets.md). Both get their own +skill in a later phase (`ag_basis_profiles`). + +## Branch — priors, pairing, fixing, assertions + +Every parameter has a default prior from the configuration YAMLs +([`../wiki/core/api/configuration.md`](../wiki/core/api/configuration.md)). Override one when +you know something the default cannot. Adapted from +`autogalaxy_workspace:scripts/guides/modeling/cookbook.py`. + +```python +bulge = af.Model(ag.lp_linear.Sersic) + +bulge.centre.centre_0 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) +bulge.centre.centre_1 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) +bulge.sersic_index = af.TruncatedGaussianPrior( + mean=4.0, sigma=1.0, lower_limit=0.8, upper_limit=5.0 +) +``` + +The centre priors say "this galaxy is within 0.1" of the image centre", which is true if the +data was cut out around it and is worth asserting because it removes a large volume of +parameter space the search would otherwise wander through. The `sersic_index` prior says +"this is an early-type" — and that is a genuine statistical commitment, not a convenience: +it changes the posterior and therefore the errors you quote. If you want to *guide* a search +without moving what you infer, the start-point API in +[`ag_configure_search`](./ag_configure_search.md) is the right tool instead. + +Three other levers reduce complexity without adding a prior: + +```python +bulge = af.Model(ag.lp_linear.Sersic) +disk = af.Model(ag.lp_linear.Exponential) + +bulge.centre = disk.centre # pair: one centre, not two (N -= 2) +bulge.sersic_index = 4.0 # fix: a de Vaucouleurs bulge (N -= 1) +bulge.effective_radius = disk.effective_radius - 0.5 # offset: a relation, not a parameter + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, disk=disk) +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +model.add_assertion( + model.galaxies.galaxy.bulge.effective_radius + < model.galaxies.galaxy.disk.effective_radius +) +``` + +Assertions are added **after** the components are collected, and they are the cleanest way to +rule out the unphysical solution a decomposition is prone to: a "bulge" that is larger than +its disk, which is the search relabelling the two components rather than fitting them. A +model can also be round-tripped through JSON with `model.dict()` and +`af.Model.from_json(file=...)`, which is how you hand a model to a collaborator or edit one +by hand. + +## Branch — nuisance parameters of the data, not the galaxy + +Some free parameters belong to the dataset: a residual background sky, a sub-pixel +astrometric offset, a small rotation. These live in `ag.DatasetModel` and sit alongside +`galaxies` rather than inside it. Adapted from +`autogalaxy_workspace:scripts/imaging/features/sky_background/modeling.py`. + +```python +dataset_model = af.Model(ag.DatasetModel) +dataset_model.background_sky_level = af.UniformPrior(lower_limit=0.0, upper_limit=5.0) + +model = af.Collection( + galaxies=af.Collection(galaxy=galaxy), dataset_model=dataset_model +) +``` + +You **must** set the sky prior by hand — unlike a light profile's priors, the right range +depends entirely on your data's units and depth, so no default can be correct. This costs one +parameter and matters for exactly the measurement galaxy structure cares most about: the +faint outer envelope that sets `effective_radius` and `sersic_index`. It is worth including +even on data you believe is sky-subtracted, precisely to check that belief. The other fields +are `grid_offset=(y, x)` and `grid_rotation_angle`, which earn their keep in multi-band work +where bands are not perfectly registered. Concept page: +[`../wiki/core/concepts/sky_background_and_operated_profiles.md`](../wiki/core/concepts/sky_background_and_operated_profiles.md). + +## Branch — modelling contaminating extra galaxies + +When a neighbour's light genuinely overlaps the target's, neither masking nor noise scaling +can separate them without removing signal you need — so model it. The convention is that its +centre is **fixed** to a value measured from the data, leaving its other parameters free. +Adapted from `autogalaxy_workspace:scripts/imaging/features/extra_galaxies/modeling.py`. + +```python +extra_galaxies_centres = ag.from_json( + file_path=DATASET_PATH / "extra_galaxies_centres.json" +) + +extra_galaxies_list = [] + +for extra_galaxy_centre in extra_galaxies_centres: + + extra_galaxy = af.Model( + ag.Galaxy, redshift=0.5, bulge=ag.lp_linear.SersicSph + ) + extra_galaxy.bulge.centre = extra_galaxy_centre + + extra_galaxies_list.append(extra_galaxy) + +model = af.Collection( + galaxies=af.Collection(galaxy=galaxy), + extra_galaxies=af.Collection(extra_galaxies_list), +) +``` + +The centres come from the JSON written during data preparation. A spherical `SersicSph` costs +about two free parameters per neighbour once its intensity is solved — cheap enough that a +handful of companions is affordable. For irregular or asymmetric companions, swap the +spherical Sersic for an MGE with a pinned centre, +`ag.model_util.mge_model_from(mask_radius=MASK_RADIUS, total_gaussians=10, +centre_fixed=tuple(extra_galaxy_centre))`, which costs the same in the linear limit while +being far more flexible. Note that `extra_galaxies` is its own top-level collection, not a +member of `galaxies`. When to model versus mask versus scale is +[`../wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md). + +## Wrap with an analysis, and check it before you search + +```python +analysis = ag.AnalysisImaging(dataset=dataset, use_jax=True) +``` + +Source: `PyAutoGalaxy:autogalaxy/imaging/model/analysis.py`. The full surface is +`AnalysisImaging(dataset, adapt_images=None, cosmology=None, settings=None, +title_prefix=None, use_jax=True)`. Pass `cosmology=ag.cosmo.Planck15()` when a derived +quantity needs angular-to-physical conversion; `settings=ag.Settings(...)` only to override +an inversion or linear-solver default; `adapt_images` only for an adaptive pixelised +reconstruction. + +Two checks before you hand this to a search, both cheap and both catching a different class +of mistake: + +```python +print(model.info) + +analysis.print_vram_use(model=model, batch_size=search.batch_size) +``` + +`model.info` catches a *model* mistake — a parameter you thought was fixed, a prior wider +than intended, a component you added twice. `print_vram_use` catches a *resource* mistake: +JAX must hold the whole batched likelihood in GPU memory, and this reports the estimate in +20-30 seconds rather than letting you discover the limit as an out-of-memory error mid-compile +(`autogalaxy_workspace:scripts/imaging/modeling.py` quotes ~0.027 GB for an MGE on a +low-resolution dataset, and more than 1 GB — occasionally more than 10 GB — for a pixelised +reconstruction at high resolution). Add the single likelihood evaluation from the first branch +when you have changed anything structural; a finite number back means the three pieces fit +together. + +## Combine — where this hands off + +- **Pick and configure the search** → [`ag_configure_search`](./ag_configure_search.md). + `MultiStartProdigy` for a fast check that this model and dataset are sensible, `Nautilus` + when you need errors you can quote. +- **Run the fit** → the run-search skill (`ag_run_search`), which owns + `search.fit(model=model, analysis=analysis)` and the output folder. +- **The search fails, stalls or returns something unphysical** → the fit-debugging skill + (`ag_debug_fit_failure`). Nine times in ten the fix is here rather than in the search: a + model too complex for the data, an unpaired centre, a missing assertion, or a contaminant + that should have been handled during preparation. +- **A smooth profile leaves structure in the residuals** → the basis-profile and pixelisation + skills (`ag_basis_profiles`, `ag_pixelization`) in a later phase; the MGE branch above is + the first step in that direction and often enough on its own. +- **Several datasets of the same galaxy** → the multi-dataset skill (`ag_multi_dataset`). + Datasets are combined through a factor graph — one `af.AnalysisFactor` per dataset, folded + into an `af.FactorGraphModel`; see + [`../wiki/core/concepts/multi_wavelength.md`](../wiki/core/concepts/multi_wavelength.md). + +Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry: the components you +chose, the parameters you fixed or paired, and *why*, are the modelling assumptions a referee +will ask about. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Realism and complexity](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_3_realism_and_complexity.ipynb): + the trade-off at the heart of this skill — what happens when a model is too simple for the + data, and what happens when it is too complex for the search. +- **General reference** — [RTD: Model cookbook](https://pyautogalaxy.readthedocs.io/en/latest/general/model_cookbook.html): + the systematic reference for `af.Model` and `af.Collection` — multiple components, multiple + galaxies, prior customisation, pairing, and many-profile bases. +- **Experienced PyAutoGalaxy user** — [workspace: imaging/modeling.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/modeling.py): + the canonical bulge-plus-disk fit end to end, including the VRAM check, run-time estimation + and the annotated output-folder layout. + +## Agent procedural checklist + +1. On real data, confirm the inspection gate was satisfied; if not, return to + `ag_prepare_imaging_data`. +2. Ask: one galaxy or several; single Sersic or a decomposition; smooth or irregular. +3. Prefer `ag.lp_linear` over `ag.lp` — the intensity is solved, not sampled. +4. Start with a single Sersic even when the goal is a decomposition; keep it as the reference. +5. Pair concentric centres, fix what physics fixes, and add assertions that rule out the + relabelling solution. +6. `print(model.info)` and check `total_free_parameters` against what you intended. +7. Build `ag.AnalysisImaging`, run one likelihood evaluation, and `print_vram_use` before a + GPU run. +8. Save the script to `scripts/`, hand off to `ag_configure_search`, and offer the + `wiki/project/` entry recording the modelling assumptions. diff --git a/skills/ag_configure_search.md b/skills/ag_configure_search.md new file mode 100644 index 0000000..edf155d --- /dev/null +++ b/skills/ag_configure_search.md @@ -0,0 +1,380 @@ +--- +name: ag_configure_search +description: Choose and configure the non-linear search that fits a galaxy model to data — `af.Nautilus` for a full posterior with errors you can quote, `af.MultiStartProdigy` for a fast maximum-a-posteriori check, `af.DynestyStatic` for isophote/ellipse fitting and as a cross-check, and the MCMC / quasi-Newton / diagnostic alternatives. Covers `n_live` and `n_batch`, `n_starts` and `n_steps`, the `iterations_per_quick_update` output cadence, the start-point initialiser API, grid searches, and — critically — the `unique_tag` / unique-identifier and resume semantics, since the identifier is a hash of the model and search but **not** of the data, so re-using a tag across datasets silently resumes the wrong fit. Use once a model and analysis exist. Not for composing the model (`ag_build_imaging_model`), not for executing the fit or reading its output folder (`ag_run_search`), and not for diagnosing a fit that already failed. +--- + +# Choosing and configuring the search + +A search is the algorithm that explores your model's parameter space, calling +`analysis.log_likelihood_function(instance)` over and over and deciding where to look next. +Every search in PyAutoFit takes the same two arguments, so swapping one for another is a +one-line change: + +```python +result = search.fit(model=model, analysis=analysis) +``` + +The choice matters for one reason above all others: **some searches return a posterior and +some return a point.** A nested sampler maps the probability density of every parameter, its +errors, and the covariances between them — so it can tell you whether an inferred effective +radius of 1.6" is 1.6 ± 0.01 or 1.6 ± 0.5, and whether it trades off against the Sersic index. +A gradient optimiser hands you a single best-fit model and *nothing else*: no errors, no +covariances. Both are useful, for different questions, and confusing them is the most +consequential mistake available here. + +The runnable tour of the whole menu is +`autogalaxy_workspace:scripts/guides/modeling/searches.py`; the catalogue with every argument +is [`../wiki/core/api/searches.md`](../wiki/core/api/searches.md), and what each family +actually does is +[`../wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md). + +## Ask + +- *"Do you need error bars, or do you just want to know whether this model fits at all?"* + This is the whole decision. Errors → `Nautilus`. A quick check → `MultiStartProdigy`. +- *"How many free parameters does the model have?"* — `print(model.info)` if you don't know. + It sets `n_live`, and it is the main driver of run time. +- *"Is this dataset one of several you will fit with the same model?"* — if yes, the + `unique_tag` section below is not optional reading; it is where results silently collide. +- *"Is JAX available, and a GPU?"* — `MultiStartProdigy` requires a JAX-traceable analysis + (`use_jax=True`) and is not available without it. + +## Branch — `af.Nautilus`, the default + +The search every workspace `modeling.py` example uses, and the one to reach for when the +answer will be quoted in a paper. Extensive testing across galaxy modelling found it the most +accurate and efficient search available. Adapted from +`autogalaxy_workspace:scripts/imaging/modeling.py`. + +```python +search = af.Nautilus( + path_prefix=Path("imaging"), + name="modeling", + unique_tag=dataset_name, + n_live=200, + n_batch=50, + iterations_per_quick_update=10000, +) +``` + +Source: `PyAutoFit:autofit/non_linear/search/nest/nautilus/`; optional dependency +`nautilus-sampler`. Reference: Lange (2023), arXiv:2306.16923. + +**`n_live` is the accuracy / run-time dial**, and essentially the only setting you need to +think about. More live points map the posterior more reliably but cost more likelihood +evaluations; fewer are faster but risk converging on a local maximum. 200 is sufficient for +the vast majority of galaxy models, and more parameters want more — a bulge-plus-disk model +at N = 9 is comfortable at 200, while a model past ~30 free parameters is worth running at +400. Dropping to 100 on a simple model to save time is a legitimate move (`modeling.py` does +exactly that) but it is a deliberate accuracy trade, so say so when you make it. + +**`n_batch` is where the GPU speed-up comes from**, and it is worth understanding because +Nautilus is *gradient-free* — it never differentiates the likelihood, so unlike the optimisers +below it does not use JAX's derivatives at all. What it does do is propose points in batches, +and when the analysis is JAX-traceable PyAutoFit evaluates the whole batch through one +`jax.vmap(jax.jit(...))` call, fitting all `n_batch` models simultaneously. It is therefore +also the main control on VRAM: a bigger batch holds more models in memory at once. Pair it +with `analysis.print_vram_use(model=model, batch_size=search.batch_size)` before a GPU run +(see [`ag_build_imaging_model`](./ag_build_imaging_model.md)). + +For parallel CPU runs, `number_of_cores` scales well: roughly linearly below 8 cores and +about `0.5 * number_of_cores` above, continuing past 50 cores — which is why supercomputing +facilities make even large models tractable. + +## Branch — `af.MultiStartProdigy`, when you want an answer now + +The search every workspace `start_here.py` example uses. Far faster than Nautilus, and +returns a single best-fit model with **no errors at all**. Adapted from +`autogalaxy_workspace:scripts/imaging/start_here.py`. + +```python +search = af.MultiStartProdigy( + path_prefix=Path("imaging"), + name="start_here", + unique_tag=dataset_name, + n_starts=48, + n_steps=300, + iterations_per_quick_update=50, + live_visual_update=False, +) +``` + +Source: `PyAutoFit:autofit/non_linear/search/mle/multi_start_gradient/`; optional dependencies +`jax` and `optax`. + +The design is what makes gradient descent usable on this problem at all. Galaxy-model +parameter spaces are multi-modal, so a *single*-start optimiser routinely descends confidently +into a local maximum — which is exactly the weakness of `LBFGS` below. +`MultiStartProdigy` launches `n_starts` independent descents from broad starting points spread +across the parameter space, runs them **all in parallel** through `jax.vmap`, and returns the +best. That wide population is what reliably finds the global maximum-likelihood basin; it is +the approach introduced for galaxy modelling by GIGA-Lens (Gu, Huang et al. 2022, +arXiv:2202.07663). Prodigy adds a *learning-rate free* update rule (Mishchenko & Defazio 2024, +arXiv:2306.06101): it estimates its own step size as it runs, so there is no `learning_rate` +for you to tune — which is precisely why it is the recommended member of the family over +`MultiStartAdam`, `MultiStartADABelief` and `MultiStartLion`, all of which want one chosen. + +Three things to hold onto: + +- **It requires `use_jax=True`.** It is gradient-based, so a NumPy analysis cannot supply what + it needs. It is therefore unavailable for `ag.AnalysisEllipse`, which is not JAX-traceable. +- **`n_steps` is a ceiling, not a target.** The search stops early once the best fit stops + improving, so expect a handful of quick updates rather than `n_steps / 50` of them. +- **There is no posterior.** If the fit reports `effective_radius = 1.6"`, this search cannot + tell you the uncertainty, and no amount of extra starts will change that. + +The workflow the workspace recommends, and the one to recommend to the user: run +`MultiStartProdigy` first to confirm cheaply that the model and data are sensible, then run +`Nautilus` on the same model when you need numbers you can quote. + +## Branch — `af.DynestyStatic`, for isophotes and as a cross-check + +The nested-sampling alternative, and the search the workspace's **ellipse** examples use, +because testing showed it the most reliable for isophote fitting specifically. + +```python +search = af.DynestyStatic( + path_prefix=Path("ellipse"), + name="fit_start", + unique_tag=dataset_name, + sample="rwalk", + nlive=50, + iterations_per_quick_update=2500, +) +``` + +Source: `PyAutoFit:autofit/non_linear/search/nest/dynesty/`. Reference: Speagle (2020), +arXiv:1904.02180. + +Note `nlive`, not `n_live` — the argument names follow each upstream sampler's own convention, +which is a real source of typos when switching. Dynesty's sampler options pass straight +through (`sample`, `walks`, `bound`, `bootstrap`, `enlarge`, `update_interval`, `facc`, +`slices`, `fmove`, `max_move`), and `af.DynestyDynamic` reallocates live points as it runs — +toward the tails for a better evidence estimate, toward the bulk for better parameters — +taking `nlive_init`, `dlogz_init` and the other `_init` variants instead. + +Two reasons to reach for it. **Isophote fitting**: use it whenever the analysis is +`ag.AnalysisEllipse`, where the gradient optimisers are unavailable anyway. **Independent +cross-check**: a posterior that reproduces under a different nested sampler is far more +convincing than one that does not, and this is the cheapest such check available. Ellipse +fitting has its own page, +[`../wiki/core/api/ellipse.md`](../wiki/core/api/ellipse.md), and gets its own skill in a +later phase (`ag_ellipse_fitting` — see [`../PENDING.md`](../PENDING.md)). + +## Branch — the rest of the catalogue + +Reach for these only when you know why. All are tabulated with full argument lists in +[`../wiki/core/api/searches.md`](../wiki/core/api/searches.md) and demonstrated in +`autogalaxy_workspace:scripts/guides/modeling/searches.py`. + +| Search | When | +|---|---| +| `af.Zeus` | Ensemble MCMC **slice** sampler; the best-performing MCMC in the workspace's tests, though still behind Nautilus. Good for characterising a posterior around a mode you already found. | +| `af.Emcee` | The familiar affine-invariant ensemble sampler. Same role as Zeus, generally a little worse on these parameter spaces. | +| `af.LBFGS` / `af.BFGS` | Single-start quasi-Newton descent via SciPy. Fast in principle, but galaxy-model parameter spaces are usually too complex to use it without careful initialisation — `MultiStartProdigy` exists to fix exactly that. | +| `af.Drawer` | Not a search: draws models from the prior and evaluates them. A diagnostic for "are my priors sane and does my likelihood run at all?", and the workspace's one-draw container for collecting many completed fits into a single output folder (`total_draws=1`). | + +MCMC and MLE searches monitor convergence differently from nested samplers — both `Emcee` and +`Zeus` accept an `af.AutoCorrelationsSettings(check_for_convergence=True, check_size=100, +required_length=50, change_threshold=0.01)` and will terminate early when the chains satisfy +it. + +## Branch — the output cadence + +Every search writes results to disk **as it runs**, using the highest-likelihood model found +so far, and `iterations_per_quick_update` sets how often. This is not free, and the two +canonical values look contradictory until you notice the unit differs: + +- `Nautilus(iterations_per_quick_update=10000)` — the unit is a likelihood evaluation, and on + a fast fit writing output every few hundred would dominate the run time. +- `MultiStartProdigy(iterations_per_quick_update=50)` — the unit is a *gradient step*, a far + coarser thing, so 50 gives a useful handful of updates. + +If the log keeps announcing that it is outputting results, raise the number. The heavier +`iterations_per_full_update` controls the full pass — all visuals plus `model.results` and +`search.summary`. Defaults for both live in `autogalaxy_assistant:config/general.yaml` under +`updates:` (set to effectively infinite there, so the per-search value governs), and the +config layer is [`../wiki/core/api/configuration.md`](../wiki/core/api/configuration.md). + +`live_visual_update=False` is the default and the right choice on a headless or cluster run. +Set it `True` to also push each quick update to a live surface — a matplotlib window that +refreshes from a script, or an in-place refresh of the Jupyter cell. The disk write happens +either way. + +## The unique identifier, `unique_tag`, and resuming — read this one + +Output lands at `output////`, where the identifier is a +32-character hash. Re-running an identical configuration **resumes** the existing fit rather +than starting over, which is a genuinely useful feature — and the sharp edge of this whole +skill, because of what goes into the hash and what does not. + +**The identifier is a hash of the model and the search. It is not a hash of the data.** + +That is not a subtlety; it is a trap with a silent failure mode. Fit galaxy A, then change +only the FITS paths and fit galaxy B with the same model and the same search: the identifier +is identical, PyAutoFit finds a completed fit at that path, and reports "Fit Already +Completed". You get galaxy A's result back, labelled as galaxy B, with no error raised. + +The fix is one argument, and it is why every workspace example passes it: + +```python +search = af.Nautilus( + path_prefix=Path("imaging"), + name="modeling", + unique_tag=dataset_name, # <-- the discriminator the data itself does not supply + n_live=200, +) +``` + +**Pass your dataset name as `unique_tag`, always.** Then the two galaxies hash differently and +land in different folders. And the corollary that catches people a second time: if the *data* +changes but its name does not — you re-reduced it, fixed the noise-map, changed the mask, +re-ran the simulator with different truth — **bump the `unique_tag`** (or delete the output +folder). Otherwise you resume a fit to data that no longer exists. + +Changing the model or the search *does* change the identifier, so you never need to manage +that case by hand: a bulge-plus-disk fit cannot collide with a single-Sersic fit, and an +`n_live=200` run cannot collide with an `n_live=400` one. It is only the data axis that is +invisible to the hash. + +Verify rather than trust when it matters. Comparing the identifier across two configurations +before launching a long run costs nothing: + +```python +search.paths.model = model +search.paths.search = search + +print(search.paths.identifier) +``` + +Two runs printing the same string will share an output folder, whatever you believe about +them. The annotated tour of what lands in that folder — `files/`, `image/`, `model.info`, +`model.results`, `search.summary` — is the `__Output Folder Layout__` section of +`autogalaxy_workspace:scripts/imaging/modeling.py`, and the run-search skill +(`ag_run_search`) owns reading it. + +## Branch — a start point, for MCMC and MLE only + +MCMC and optimiser searches have a *location* in parameter space, so you can say where they +begin. Nested samplers draw from the prior and cannot use this API. Adapted from +`autogalaxy_workspace:scripts/guides/modeling/searches.py`, where an early-type galaxy is +started near n = 4. + +```python +initializer = af.InitializerParamBounds( + { + model.galaxies.galaxy.bulge.centre_0: (-0.01, 0.01), + model.galaxies.galaxy.bulge.centre_1: (-0.01, 0.01), + model.galaxies.galaxy.bulge.effective_radius: (0.9, 1.1), + model.galaxies.galaxy.bulge.sersic_index: (3.9, 4.1), + } +) + +search = af.Emcee( + path_prefix=Path("imaging"), + name="start_point", + unique_tag=dataset_name, + nwalkers=50, + nsteps=500, + initializer=initializer, +) +``` + +Parameters you do not name are drawn from their priors. `af.InitializerBall(lower_limit=0.49, +upper_limit=0.51)` is the other common choice — a tight ball in unit-prior space around the +prior centres, and the recommended initialisation for both MCMC samplers. + +The statistical point here is worth stating explicitly, because there is a tempting shortcut +that is not equivalent. You *could* achieve a similar effect by tightening priors instead — +a narrow `af.TruncatedGaussianPrior` on `sersic_index`, say. But priors change the posterior, +and therefore change the errors you quote and the model you infer. The start-point API moves +where the search *looks* without moving what you *infer*. Prefer it whenever your knowledge is +about where the answer probably is rather than about what the answer must be. + +## Branch — mapping a parameter on a grid + +To scan a parameter on a fixed grid rather than marginalising over it — a light centre, a +fixed `sersic_index`, anything you want mapped — wrap a search: + +```python +grid_search = af.SearchGridSearch( + search=search, number_of_steps=4, number_of_cores=1, result_output_interval=100 +) +``` + +One child fit runs per grid cell and an `af.GridSearchResult` comes back. Source: +`PyAutoFit:autofit/non_linear/grid/`. This is the one search whose defaults *do* live in a +YAML — `autogalaxy_assistant:config/non_linear/GridSearch.yaml` sets `number_of_cores` and +`step_size` in unit-prior values, and it is the only file in that folder because per-search +defaults ship inside PyAutoFit itself +(`autogalaxy_assistant:config/non_linear/README.md`). + +## A different lever: chaining instead of a better search + +When a model is too complex for any single search, the answer is usually not a bigger +`n_live` — it is to fit a simpler model first and pass its posterior forward as the priors of +the complex one. Fit a single Sersic, then chain into a bulge-plus-disk or pixelised fit, +which starts from a region of parameter space that is already good. +`autogalaxy_workspace:scripts/guides/modeling/chaining.py` is the walkthrough, the +`width_modifier` entries in the prior YAMLs control how wide the passed priors become, and it +gets its own skill in a later phase (`ag_chain_searches`). + +## Picking one at a glance + +| Goal | Pick | +|---|---| +| Quick check that model + data are sensible | `MultiStartProdigy(n_starts=48, n_steps=300)` | +| Results you will quote, model under ~30 free parameters | `Nautilus(n_live=200)` | +| Production run, complex or multi-modal model | `Nautilus(n_live=400)`, higher `n_batch` if VRAM allows | +| Bayesian evidence comparison between models | `Nautilus` or `DynestyStatic` | +| Ellipse / isophote fitting | `DynestyStatic(sample="rwalk", nlive=50)` | +| Independent cross-check of a posterior | `DynestyStatic` or `Zeus` | +| Posterior refinement around a known mode | `Zeus` or `Emcee` with an `initializer` | +| Map a parameter on a fixed grid | `af.SearchGridSearch` | +| Check the priors are sane and the likelihood runs | `af.Drawer` | + +## Combine — where this hands off + +- **Run the fit** → the run-search skill (`ag_run_search`), which owns + `search.fit(model=model, analysis=analysis)`, announcing the output folder at launch, and + what to open first while it runs. +- **The model isn't right yet** → [`ag_build_imaging_model`](./ag_build_imaging_model.md). + Changing the search rarely rescues a badly-specified model, and a search that will not + converge is more often a model problem than a sampler problem. +- **It ran and failed, stalled, or returned something unphysical** → the fit-debugging skill + (`ag_debug_fit_failure`), which is also where "Fit Already Completed" on data you thought was + new gets diagnosed. +- **Read the posterior that comes back** → the results-loading and fit-plotting skills + (`ag_load_results`, `ag_plot_fit`); the `Samples` API and how errors are computed is + [`../wiki/core/concepts/samples_and_posteriors.md`](../wiki/core/concepts/samples_and_posteriors.md). + +Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry recording the search and +its settings — `n_live` in particular, since it is an accuracy claim about the posterior you +are about to publish. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: The non-linear search](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_1_non_linear_search.ipynb): + what a non-linear search actually is — parameter spaces, likelihoods, priors — and the + statistical foundation under every setting on this page. The + [optional chapter](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_optional/tutorial_searches.ipynb) + tours the sampler zoo. +- **General reference** — [RTD: Configs](https://pyautogalaxy.readthedocs.io/en/latest/general/configs.html): + how the configuration files customise searches, visualisation and output, and where each is + looked up. +- **Experienced PyAutoGalaxy user** — [workspace: guides/modeling/searches.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/guides/modeling/searches.py): + every available search instantiated with its settings spelled out, plus the start-point API. + +## Agent procedural checklist + +1. Ask the one question that decides it: errors needed, or a fast check? +2. `Nautilus(n_live=200)` by default; `MultiStartProdigy` for speed (needs `use_jax=True`); + `DynestyStatic` for isophotes or a cross-check. +3. Set `n_live` from the model's free-parameter count, and say when you trade it down. +4. **Always pass `unique_tag=dataset_name`** — and bump it when the data changes but the name + does not. +5. Print `search.paths.identifier` before a long run whenever a collision is plausible. +6. Set `iterations_per_quick_update` to the search's unit — thousands for Nautilus, tens for a + gradient optimiser. Keep `live_visual_update=False` on headless runs. +7. Run `analysis.print_vram_use(model=model, batch_size=search.batch_size)` before a GPU run. +8. Hand off to `ag_run_search`, and offer the `wiki/project/` entry recording the settings. diff --git a/skills/ag_debug_fit_failure.md b/skills/ag_debug_fit_failure.md new file mode 100644 index 0000000..847699b --- /dev/null +++ b/skills/ag_debug_fit_failure.md @@ -0,0 +1,504 @@ +--- +name: ag_debug_fit_failure +description: Triage a galaxy fit that crashed, stalled, or finished with residuals and parameters you do not believe. Works through a failure taxonomy — bad environment, bad data, bad model, bad priors, bad search settings, or a stale result that was never re-run — and gives the probes that separate them: hand-evaluating the likelihood at prior medians, checking the mask extent and over-sampling scheme, reading the normalised-residual and chi-squared maps, walking the likelihood step by step, forcing the NumPy path to get a readable stack trace, and the `PYAUTO_TEST_MODE` short-circuit loop for fast iteration. Also covers the two silent failures: a resumed fit whose identifier ignored the data, and a cached result mistaken for a new one. Use when a fit raises, hangs, returns nonsense, or converges somewhere unphysical. Not for the first composition of a model (`ag_build_imaging_model`), not for routine result inspection (`ag_load_results`), and not for an environment that will not import at all (`ag_setup_environment`). +--- + +# When the fit goes wrong + +A galaxy fit fails in a small number of characteristic ways, and the debugging cost is +dominated by mis-attribution: hours spent widening priors when the mask was truncating the +outer isophotes, or rebuilding a model when the script was quietly resuming a cached result +from last week. So the first move is never a fix — it is to establish *which kind* of failure +this is. + +The discipline that makes that possible is cheap probes before expensive ones. A single +likelihood evaluation costs a fraction of a second and tells you whether the model and data +can talk to each other at all. A `PYAUTO_TEST_MODE=2` run costs seconds and proves the whole +script executes. Only after those pass is it worth spending an hour of sampling to find out +whether the *inference* is the problem. + +## Ask + +Get the symptom precisely, because the taxonomy branches on it: + +- *"Did it raise, hang, or complete?"* A raise has a traceback to read. A hang is usually + parallelisation, a JIT compile you mistook for a stall, or an update cadence writing output + constantly. A completed fit that looks wrong is the hardest case and needs the residuals. +- *"What makes you say it is wrong?"* Coherent residuals, an unphysical parameter, an error + bar that is implausibly tight or implausibly wide, or the number simply disagreeing with the + literature — each points somewhere different. +- *"Simulated data or real observations?"* On simulated data you know the truth, which makes + everything below faster. If the user has real data and no simulated control, running one is + often the fastest path to an answer ([`ag_simulate_dataset`](./ag_simulate_dataset.md)). +- *"Has this exact script run before?"* This decides whether the stale-result trap is in play. + +## The taxonomy + +Six categories, ordered by how cheap they are to rule out: + +| Category | Tell | First probe | +|---|---|---| +| **Stale output** | The fit "finished" implausibly fast, or the result is identical after you changed something | Delete the `` folder and re-run | +| **Environment** | An import error, a numba or matplotlib cache error, a parallelisation error the moment the search starts | `--check-install`, the cache variables, the main guard | +| **Bad data** | Residuals dominated by one blob, or a suspiciously large chi-squared everywhere | Plot the dataset; check the noise-map and the mask | +| **Bad model** | Coherent, structured residuals with a physical shape | Read the residual pattern; add or change a component | +| **Bad priors** | The best fit sits on a prior boundary, or the search never moves off its start | `model.info`; hand-evaluate the likelihood | +| **Bad search settings** | Wide posteriors, a low log evidence, different runs disagreeing, no errors at all | `search.summary`; raise `n_live`, or change search | + +The rest of this skill is one branch per category. + +## Branch — stale output, the silent one + +Rule this out first, every time, because it is invisible and it invalidates everything else you +would conclude. + +A fit's output folder is `output/////`, and +`` is derived from the model, the search settings and the dataset identifier. An +identical configuration therefore **resumes** rather than restarting +(`autogalaxy_workspace:scripts/imaging/modeling.py` `__Unique Identifier__`). That is a +feature for an interrupted overnight run and a trap the rest of the time, because **the +identifier does not hash the pixel values of your data**. Swap `data.fits` for a different +galaxy while keeping the model, search and `unique_tag` the same, and the next run finds a +completed fit at the same path and hands it back — in seconds, with the previous galaxy's +parameters, and no warning that anything is wrong. + +The symptoms: a fit that completes far faster than it should; a result that does not change +after you deliberately changed the data; `model.results` timestamps older than the run you +just launched. + +```bash +# Prove it rather than reasoning about it. +rm -rf output//// +``` + +Then re-run. If the answer changes, that was the bug. Going forward, change the `unique_tag` +whenever the data changes — that is what it is for. + +The related trap: test-mode output. Any active `PYAUTO_TEST_MODE` level inserts a `test_mode` +segment straight after the output root, precisely so a smoke run cannot short-circuit a later +real fit ([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)). If you +find yourself inspecting a result whose path contains `test_mode`, stop: every parameter in it +is a wiring check. Also delete any simulated `dataset/` when you toggle +`PYAUTO_SMALL_DATASETS`, or a full-resolution dataset on disk gets reused against capped grids +and produces shape errors that look like library bugs. + +## Branch — environment and the errors that fire at launch + +If the stack itself is the problem, nothing downstream is diagnostic. Confirm the install +first: + +```bash +python autoassistant/audit_skill_apis.py --check-install +``` + +Exit `0` is ready; `2` means the packages are absent from *this* interpreter; `3` means they +were found but an import raised. Repairing it is +[`ag_setup_environment`](./ag_setup_environment.md), with the routes in +[`wiki/core/operations/installation.md`](../wiki/core/operations/installation.md). + +Three failures that are environment, not science: + +- **A numba or matplotlib cache error** buried in a long trace — the libraries cannot write + their caches. Fix with `NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib` + ([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)). +- **An error the instant the search starts**, from Python's process spawning rather than from + your model. Wrap the fit in a function under a main guard, exactly as + `autogalaxy_workspace:scripts/guides/modeling/bug_fix.py` does — it is the same + `modeling.py` fit with this one structural change: + + ```python + def fit(): + from autogalaxy import jax_wrapper # Sets the JAX environment before other imports + + import autofit as af + import autogalaxy as ag + + # ... dataset, mask, over-sampling, model, search, analysis ... + result = search.fit(model=model, analysis=analysis) + + + if __name__ == "__main__": + fit() + ``` + + It applies to every dataset type, so adopt it for any modeling script that hits the error. + That script also names the project's support channel if the fix is not enough. +- **A stale API symbol** — an `AttributeError` on a PyAuto\* name recalled from an older + release. Do not guess a replacement; check it: + + ```bash + python autoassistant/audit_skill_apis.py --file scripts/run_fit.py + ``` + + [`ag_audit_skill_apis`](./ag_audit_skill_apis.md) owns that check and its bypass. + +### Is it a JAX problem? + +The single most useful lever, because it splits the space in two: + +```python +analysis = ag.AnalysisImaging(dataset=dataset, use_jax=False) +``` + +or, without editing code, `PYAUTO_DISABLE_JAX=1` to force the NumPy path on every analysis +(`autogalaxy_workspace:scripts/guides/using_jax.py` `__Disabling JAX__`). NumPy stack traces +are far easier to read than compiled ones, and you can drop a debugger or a `print` into code +JAX would otherwise compile away. If the failure disappears on NumPy, it is a JAX problem; if +it survives, JAX was never involved and you have halved the search space. + +Two JAX-specific failure shapes worth recognising: + +- **An out-of-memory error at compile time or on the first likelihood call** on a GPU. Check + the footprint before a long run: + `analysis.print_vram_use(model=model, batch_size=search.batch_size)` + (`autogalaxy_workspace:scripts/imaging/modeling.py` `__VRAM Use__`), and lower the batch + size if it is close to the card's limit. +- **A long pause at the start that looks like a hang.** The first evaluation compiles the + likelihood and its gradient; every step after that re-uses the compiled result. The log + line announcing the compile is the tell + (`autogalaxy_workspace:scripts/imaging/start_here.py` `__JAX__`). A persistent compilation + cache means this is paid once per machine, not once per process + ([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)). + +## Branch — bad data + +Plot the dataset before you plot anything else. On real observations this is not optional — +the inspection gate in [`../AGENTS.md`](../AGENTS.md) requires it before a fit is even +composed, and a fit that skipped it is a prime suspect. + +```python +aplt.subplot_imaging_dataset( + dataset=dataset, + output_path=PLOT_DIR, + output_filename="dataset", + output_format="png", +) +``` + +Four things to check, in order: + +1. **Contaminants.** A neighbouring galaxy, a foreground star or a reduction artefact inside + the mask contributes to the likelihood and biases every parameter. The fix is to scale its + noise rather than cut the pixels out — `dataset.apply_noise_scaling(mask=mask_extra_galaxies)` + keeps the pixels in the fit but makes them contribute negligibly, which avoids the + discontinuities that removing pixels creates + (`autogalaxy_workspace:scripts/imaging/modeling.py` `__Extra Galaxies Noise Scaling__`). + Alternatively shrink the circular mask so the contaminant falls outside it. Concept: + [`wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md); + the tooling is [`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md). +2. **The mask extent.** Too small and it truncates the outer isophotes, which biases the + effective radius and the Sersic index directly — the classic cause of a "wrong" Sersic + index. Too large and it drags sky and neighbours in, slowing the fit for no information. + [`wiki/core/concepts/grids_and_masks.md`](../wiki/core/concepts/grids_and_masks.md) + "Choosing the radius is a science decision". +3. **The noise-map.** If it is underestimated, chi-squared is uniformly too large and the + posterior is spuriously tight; if overestimated, the reverse. A residual map that is + featureless but scaled wrong across the whole image points here rather than at the model. +4. **`pixel_scales`.** Wrong by a factor and every angular quantity is wrong by that factor, + while the fit looks perfectly healthy. Check it against the instrument. + +## Branch — bad model + +This is what coherent residuals mean. Rebuild the fit and read the maps: + +```python +fit = result.max_log_likelihood_fit # or ag.FitImaging(dataset=dataset, galaxies=galaxies) + +print(f"chi_squared = {fit.chi_squared}") +print(f"log_likelihood = {fit.log_likelihood}") +print(f"figure_of_merit = {fit.figure_of_merit}") + +aplt.plot_array( + array=fit.normalized_residual_map, + title="Normalized Residual Map", + symmetric=True, + output_path=PLOT_DIR, + output_filename="normalized_residual_map", + output_format="png", +) +aplt.plot_array( + array=fit.chi_squared_map, + title="Chi-Squared Map", + output_path=PLOT_DIR, + output_filename="chi_squared_map", + output_format="png", +) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/plot.py` `__Fit Figures__`. Read the +**normalised** residuals first: they are the only map calibrated against the noise, so they +answer *significance* rather than *magnitude*. Scatter within roughly ±3 is consistent with +noise; coherent structure outside that is a model failure however faint it looks. + +Two quick numerical summaries worth printing alongside the figures — chi-squared per unmasked +pixel (order unity for a good fit), and the fraction of pixels beyond 3σ (a few per cent for a +good fit, tens of per cent for a badly wrong model). + +The residual patterns and what each implies are tabulated in +[`ag_plot_fit`](./ag_plot_fit.md) "The inspection discipline". The model responses: + +- **Ring at one radius** → the radial profile is too rigid. One Sersic where a bulge plus a + disk is needed, or a Sersic index that was fixed. Free it, or add the second component. +- **Four-lobed alternating pattern** → the ellipticity or position angle is wrong, or the + isophotes twist with radius. [`wiki/core/concepts/ellipse_fitting_and_multipoles.md`](../wiki/core/concepts/ellipse_fitting_and_multipoles.md). +- **Central excess** → coarse over-sampling, or a nuclear component the model lacks. Check the + over-sampling scheme before adding physics (below). +- **Clumpy, asymmetric structure no smooth profile can absorb** → this is where a + many-component basis or a pixelised reconstruction earns its keep: + [`wiki/core/concepts/linear_light_profiles_and_mge.md`](../wiki/core/concepts/linear_light_profiles_and_mge.md), + [`wiki/core/concepts/shapelets.md`](../wiki/core/concepts/shapelets.md), + [`wiki/core/concepts/inversions_and_pixelizations.md`](../wiki/core/concepts/inversions_and_pixelizations.md). +- **A uniform offset across the whole image** → an unmodelled sky background: + [`wiki/core/concepts/sky_background_and_operated_profiles.md`](../wiki/core/concepts/sky_background_and_operated_profiles.md). + +**Check the over-sampling before blaming the physics.** A steep central profile evaluated once +per pixel is systematically wrong at the centre, and the symptom is exactly the residual a +missing nuclear component would leave: + +```python +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=dataset.grid, + sub_size_list=[8, 4, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0)], +) +dataset = dataset.apply_over_sampling(over_sample_size_lp=over_sample_size) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/modeling.py` `__Over Sampling__`. Re-plot +the dataset subplot afterwards — its lower panels show the scheme that was actually applied, +so you can confirm it rather than assume it. The full treatment, including why the outer bins +never drop to a single sub-pixel, is +[`wiki/core/concepts/grids_and_masks.md`](../wiki/core/concepts/grids_and_masks.md) +"Over-sampling", with the runnable guide at +`autogalaxy_workspace:scripts/guides/advanced/over_sampling.py`. + +Adding complexity is [`ag_build_imaging_model`](./ag_build_imaging_model.md). Do it one +component at a time and re-fit: a model that gained three components at once tells you nothing +about which one mattered. + +## Branch — bad priors, and hand-evaluating the likelihood + +The most direct probe in the whole stack. A model and an analysis are enough — no search, no +waiting: + +```python +instance = model.instance_from_prior_medians() +log_likelihood = analysis.log_likelihood_function(instance=instance) + +print(f"log likelihood at prior medians = {log_likelihood}") +``` + +From [`wiki/core/api/analysis_objects.md`](../wiki/core/api/analysis_objects.md) "Evaluating +the likelihood by hand". `log_likelihood_function(instance)` is the single contract between a +model and data, so what it returns is diagnostic on its own: + +- **It raises.** The failure is in the forward model or the data shapes, not in the inference. + Force `use_jax=False` first for a readable trace, then walk the likelihood step by step + (next branch). +- **It returns `nan` or `-inf`.** Something in the model is degenerate at those parameter + values — a zero or negative size, a profile whose centre sits outside the mask, a noise-map + with zeros. +- **It returns a finite but wildly bad value, and the fit's best likelihood is barely better.** + The search never found anything, which usually means the priors do not contain the answer. +- **It returns something close to the fit's best value.** The priors are fine and the search + is doing its job; look elsewhere. + +Then compare against a known-good instance. On simulated data you have the truth, so evaluate +the likelihood there: if the true model scores much better than anything the search found, the +priors or the search settings are at fault, not the model. This is the cleanest experiment +available, and it is the reason to keep a simulated control alongside real data +([`ag_simulate_dataset`](./ag_simulate_dataset.md)). + +Read the priors themselves rather than recalling them: + +```python +print(model.info) +``` + +Three specific things to look for. **A best-fit parameter sitting on a prior boundary** means +the prior is truncating the posterior — the answer may be outside it entirely. **A prior far +wider than physically sensible** wastes the search's effort and invites multi-modality. And +**a prior narrowed to make a search behave** is a scientific choice, not a convenience: it +changes the posterior and therefore the answer. If you want to steer where a search *starts* +without changing the posterior, that is a start point, and only the MCMC and optimisation +searches accept one — nested samplers draw from the prior by construction +([`wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md) +"Priors and start points"). Defaults come from configuration +([`wiki/core/api/configuration.md`](../wiki/core/api/configuration.md)), and overriding them +per parameter is [`ag_build_imaging_model`](./ag_build_imaging_model.md). + +Also confirm the galaxy is near the coordinate origin. The default priors assume a centre near +(0.0", 0.0"); a galaxy several arcseconds off-centre needs either re-centred data or +explicitly overridden centre priors (`autogalaxy_workspace:scripts/imaging/modeling.py` +`__Coordinates__`). + +## Branch — walking the likelihood step by step + +When the hand evaluation raises or returns something inexplicable, the step-by-step +walkthrough is the ground truth. `autogalaxy_workspace:scripts/imaging/likelihood_function.py` +computes a galaxy model's log likelihood one NumPy operation at a time — the elliptical +coordinate transform, each light profile's image, the sum into a galaxy image, the blurring +grid and PSF convolution, then: + +``` +model_data = convolved_image_2d +residual_map = data - model_data +normalized_residual_map = residual_map / noise_map +chi_squared_map = normalized_residual_map ** 2 +chi_squared = sum(chi_squared_map) +noise_normalization = sum(log(2 * pi * noise_map ** 2)) +figure_of_merit = -0.5 * (chi_squared + noise_normalization) +``` + +and finishes by showing that `ag.FitImaging(dataset=dataset, galaxies=galaxies).figure_of_merit` +reproduces the same number. Running it against your own dataset localises the failure to one +step, which is far faster than reading a compiled trace. + +There are matching walkthroughs for the features whose likelihood is genuinely different, and +the difference is where their failures live: +`autogalaxy_workspace:scripts/imaging/features/linear_light_profiles/likelihood_function.py` +(the mapping matrix, data vector, curvature matrix and the positive-only reconstruction), +`.../multi_gaussian_expansion/likelihood_function.py` (the same for a Gaussian basis), and +`.../pixelization/likelihood_function.py` (the regularisation matrix and the Bayesian +complexity terms). If a fit with linear profiles or a pixelisation fails where a standard +profile fit succeeds, the extra linear-algebra steps in those scripts are the place to look. + +A real fit does not run that NumPy code — it runs the same calculation compiled. To probe the +compiled path exactly as a search drives it, `autogalaxy_workspace:scripts/guides/using_jax.py` +`__Custom Likelihood Functions__` shows the `Fitness` route, and explains why batching a vector +of parameters is a stricter test than one concrete call: a single call can quietly succeed on +code where NumPy is leaking through, and then break as soon as the search batches. + +## Branch — bad search settings + +The model and data are fine; the sampler did not do its job. Read the search's own report +first: `search.summary` in the output folder, plus `result.samples.log_evidence` and the +posterior widths. + +- **Posteriors far wider than the data should allow, or a low log evidence** → too few live + points for the model's dimensionality. `n_live=200` covers most galaxy models; raise it for + a complex one ([`wiki/core/api/searches.md`](../wiki/core/api/searches.md)). +- **No errors at all** → the fit used a gradient optimiser, which is maximum-a-posteriori and + returns one model. Re-run with `Nautilus` if the science needs uncertainties + (`autogalaxy_workspace:scripts/imaging/start_here.py` `__Posterior__`). +- **Two runs of the same configuration disagreeing** → a multi-modal parameter space with a + search that is finding different modes. Either raise `n_live`, or start simple and chain into + the complex model with the simple fit's posterior as priors + ([`wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md) + "Search chaining"; the runnable version is + `autogalaxy_workspace:scripts/guides/modeling/chaining.py`). +- **Run time far longer than expected** → total cost is (likelihood evaluation time) × + (number of evaluations), and dimensionality dominates. A many-Gaussian expansion with six + free parameters routinely beats a multi-component decomposition with thirteen, despite the + slower per-evaluation cost. [`wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md) + "Run-time estimation". +- **The log constantly reporting that it is outputting results** → `iterations_per_quick_update` + is too low and visualisation is eating the run. Raise it + ([`ag_run_search`](./ag_run_search.md)). + +Choosing and tuning the search is [`ag_configure_search`](./ag_configure_search.md). + +## Branch — the fast iteration loop + +While you are changing a script's structure, do not pay for inference. `PYAUTO_TEST_MODE` has +four levels ([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)): + +```bash +# Level 2 — sampler bypassed, likelihood called exactly once. The fastest proof that the +# model composes, the shapes agree and the likelihood evaluates. +PYAUTO_TEST_MODE=2 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/run_fit.py + +# Level 1 — a real but minimal search. Slower, and the only level that writes model.results, +# image/fit.png and a loadable samples.csv. +PYAUTO_TEST_MODE=1 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/run_fit.py +``` + +Level 2 does **not** write the fit products, which is worth knowing before you go looking for +`image/fit.png` after a level-2 run. Level 3 skips the likelihood as well, so it only checks +the non-fitting scaffolding. `PYAUTO_SKIP_CHECKS=1` is available when a smoke run trips a +runtime validation on data too small to be sane — never set it for a real fit, since those +checks are what stop a silently degenerate inversion. And **never quote a parameter measured +in test mode**: the sampler was truncated or mocked, so the numbers are wiring evidence, not +measurements. + +The loop that works: level 2 until the script runs clean → level 1 until the products look +sane → one full run on a small mask or a simulated dataset → the production fit. + +## Branch — when nothing above explains it + +Two escalations, in order. + +**Reduce to a control.** Simulate a dataset from a model you choose, fit it with that same +model, and confirm you recover the truth +([`ag_simulate_dataset`](./ag_simulate_dataset.md)). If the control fails, the problem is in +the code or the setup. If the control succeeds, the problem is in the real data or its +preparation — and you have narrowed it to a half of the space you could not distinguish +before. This is the single highest-value move in this skill and it is under-used. + +**Check the library, not yourself.** If a symbol does not resolve or an idiom that reads +correctly does not work, run the audit — a construction can be retired while every token in it +still imports: + +```bash +python autoassistant/audit_skill_apis.py --scope scripts +python autoassistant/audit_skill_apis.py --lint-idioms +``` + +[`ag_audit_skill_apis`](./ag_audit_skill_apis.md) owns both. Only after those come back clean +is "this might be a library bug" a reasonable hypothesis — at which point +[`contribute-upstream`](./contribute-upstream.md) is the route to filing it with a minimal +reproduction. + +## Combine + +- [`ag_plot_fit`](./ag_plot_fit.md) — the residual-pattern table this skill's model branch + routes through, and the figures that make a diagnosis visible. +- [`ag_build_imaging_model`](./ag_build_imaging_model.md) — acting on a "bad model" or "bad + priors" verdict. +- [`ag_configure_search`](./ag_configure_search.md) — acting on a "bad search settings" + verdict. +- [`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md) — acting on a "bad data" verdict: + masks, noise scaling, contaminant handling. +- [`ag_simulate_dataset`](./ag_simulate_dataset.md) — the control experiment. +- [`ag_setup_environment`](./ag_setup_environment.md) — the install will not import at all. +- [`ag_run_search`](./ag_run_search.md) — re-launching once the cause is fixed. + +Record the diagnosis, not just the fix. Offer (default-yes) a dated +`wiki/project/YYYY-MM-DD-.md` entry naming the symptom, the probe that localised it, and +what changed — per [`_style.md`](./_style.md) property #5. A debugging session that is not +written down gets repeated. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Dealing with failure](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_4_dealing_with_failure.ipynb): + why searches fail on galaxy models, what a local maximum looks like in practice, and how + model complexity and prior choice interact. +- **General reference** — [RTD: Likelihood function](https://pyautogalaxy.readthedocs.io/en/latest/general/likelihood_function.html): + the likelihood the library evaluates, and pointers into the step-by-step walkthroughs. +- **Experienced PyAutoGalaxy user** — [workspace: guides/modeling/bug_fix.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/guides/modeling/bug_fix.py): + the main-guard fix for the parallelisation error, as a complete runnable fit. + +## Agent procedural checklist + +1. Get the symptom: raised, hung, or completed-but-wrong; and what makes it wrong. +2. Rule out stale output first — check timestamps, then `rm -rf` the `` folder and + re-run. Check for a `test_mode` path segment. +3. Confirm the environment with `--check-install`; apply the cache variables and the main guard + if either applies. +4. Toggle `use_jax=False` (or `PYAUTO_DISABLE_JAX=1`) to split JAX from non-JAX failures and + get a readable trace. +5. Plot the dataset: contaminants, mask extent, noise-map, `pixel_scales`. +6. Hand-evaluate `analysis.log_likelihood_function(instance=model.instance_from_prior_medians())` + before running any search. +7. Read `model.info` for boundary-hugging and over-wide priors; confirm the galaxy is near the + origin. +8. Rebuild the fit and read the **normalised residual** map, then the chi-squared map; check + the over-sampling scheme before adding physics. +9. Read `search.summary` and the log evidence for search-settings failures. +10. Iterate with `PYAUTO_TEST_MODE=2`, then `=1`; never quote a test-mode parameter. +11. If unresolved, build a simulated control and fit it; run the symbol and idiom audits before + suspecting a library bug. +12. Offer the `wiki/project/` entry recording symptom, probe and fix. diff --git a/skills/ag_load_results.md b/skills/ag_load_results.md new file mode 100644 index 0000000..fd7434b --- /dev/null +++ b/skills/ag_load_results.md @@ -0,0 +1,435 @@ +--- +name: ag_load_results +description: Get a completed fit's results back into Python and turn them into science — the in-session `Result` object (max-log-likelihood instance, galaxies and fit), direct loading of one fit's output folder via `from_json` and `SamplesNest.from_table`, the `Samples` API for medians and errors, and the aggregator for a whole sample of fits with its `ag.agg` generators, queries and CSV/FITS/PNG workflow exports. Also covers the traps: linear light-profile intensities are solved rather than sampled and are absent from the samples, derived-quantity errors need posterior draws, and test-mode results are not physically meaningful. Use when a search has finished (or is mid-run) and the user wants numbers, errors, derived quantities, a results table, or a comparison across fits. Not for running the fit (`ag_run_search`), not for rendering figures of a single fit (`ag_plot_fit`), and not for diagnosing a fit that converged somewhere wrong (`ag_debug_fit_failure`). +--- + +# Reading a finished fit + +The search is over; the question is what it measured. This skill is the bridge from an +`output/` folder to a sentence you could put in a paper: *"the bulge has an effective radius +of 1.62 ± 0.08 arcseconds and a Sersic index of 3.9 ± 0.3"*. + +Loading is never the goal on its own — the user already wants something specific: a +parameter with errors, a derived quantity the fit did not parameterise directly, the +residuals, a table across fifty galaxies. The job here is to establish which, load only that, +and be honest about what the fit can and cannot support. + +Two things about the statistics are worth stating before any API. First, the default view of a +result is the **maximum log likelihood** model — the single best-fitting parameter vector. +That is not the same as the posterior, and quoting it without errors is quoting half a +measurement. Second, a nested-sampling run gives you the full posterior, so errors, +covariances and the Bayesian evidence are all available; a gradient optimiser such as +`MultiStartProdigy` is a maximum-a-posteriori method and gives you **no errors at all** +(`autogalaxy_workspace:scripts/imaging/start_here.py` `__Posterior__`). If the fit used the +optimiser and the user wants uncertainties, the answer is to re-run with `Nautilus`, not to +manufacture them. + +The concept page behind everything below is +[`wiki/core/concepts/samples_and_posteriors.md`](../wiki/core/concepts/samples_and_posteriors.md); +the aggregator's full API surface is +[`wiki/core/api/aggregator.md`](../wiki/core/api/aggregator.md). + +## Ask + +Two questions, and they choose the branch outright: + +- **What do you want out of it?** *"A parameter and its error"* → the `Samples` branch. + *"The residuals / how good is the fit"* → the fit branch, then + [`ag_plot_fit`](./ag_plot_fit.md). *"A physical quantity in kpc or solar luminosities"* → + the derived-quantities branch. *"A table over all my galaxies"* → the aggregator branch. +- **One fit, or many?** One fit you already have a path to → simple loading. More than a + handful, or you want to iterate and filter → the aggregator, which is generator-based and + keeps memory bounded no matter how many fits the tree holds. + +If they give you a path, list it and confirm which sub-folder: the layout is +`output/////`, and a parent path can hold many +fits. If the path contains a `test_mode` segment, say so immediately — every number in it is a +wiring check, not a measurement +([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)). + +## Branch — the result you already have in the session + +If the fit just ran, `search.fit(...)` returned everything, no disk access needed: + +```python +result = search.fit(model=model, analysis=analysis) + +print(result.info) # human-readable summary, parameters + errors +print(result.max_log_likelihood_instance) # best-fit model as concrete objects + +galaxies = result.max_log_likelihood_galaxies # ag.Galaxies at the best-fit values +fit = result.max_log_likelihood_fit # ag.FitImaging: model image, residuals, chi² +samples = result.samples # the full posterior +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/modeling.py` `__Result__` and +`autogalaxy_workspace:scripts/guides/results/start_here.py`. The full member list — +`result.samples_summary`, `result.model`, `result.instance`, `result.dataset`, `result.grids`, +`result.unmasked_model_image`, the per-galaxy image dictionaries — is tabulated in +[`wiki/core/api/analysis_objects.md`](../wiki/core/api/analysis_objects.md) +"What the result gives you". Sources: `PyAutoGalaxy:autogalaxy/analysis/result.py` and +`PyAutoGalaxy:autogalaxy/imaging/model/result.py`. + +`result.samples_summary` is worth knowing about: it holds pre-computed maximum-likelihood +values and 1σ/3σ errors, so reaching for it instead of re-deriving from `samples` is +substantially faster and gives the same answer. + +## Branch — loading one fit from disk + +Everything the `Result` holds was also written to the output folder, and each file loads back +into a full Python object in one line — much cheaper than re-running the search: + +```python +from pathlib import Path + +import autofit as af +import autogalaxy as ag +from autogalaxy import from_json + +FIT_PATH = Path("output") / "imaging" / "" / "sersic" / "" +files_path = FIT_PATH / "files" +image_path = FIT_PATH / "image" + +galaxies = from_json(file_path=files_path / "galaxies.json") +model = from_json(file_path=files_path / "model.json") + +samples = af.SamplesNest.from_table(filename=files_path / "samples.csv", model=model) +print(samples.max_log_likelihood()) +``` + +Adapted from `autogalaxy_workspace:scripts/guides/results/start_here.py` +"Simple Loading". If the fit ran in the same session, `search.paths.output_path` already +points at the right folder, so you never have to know the hash. + +Two things to know about what comes back: + +- **`galaxies.json` deserialises to a plain Python list of `Galaxy` objects, not an + `ag.Galaxies`.** The list is enough for `ag.FitImaging(dataset=dataset, galaxies=galaxies)` + and for indexing (`galaxies[0].bulge`), but it has no `image_2d_from`. Wrap it when you need + the collection's own methods: `galaxies = ag.Galaxies(galaxies=galaxies)`. +- **`model.json` is the *prior* model, with free parameters** — not the best-fit instance. It + is what you pass to `SamplesNest.from_table` so the sample columns can be mapped back onto + named parameters, and it is worth printing (`model.info`) to confirm what was actually + fitted. + +The `image/` folder holds the imaging products as FITS, loadable with the standard readers — +`ag.Imaging.from_fits` for `dataset.fits` (multi-HDU, so pass the `*_hdu` arguments) and +`ag.Array2D.from_fits` for the per-galaxy image stacks. The annotated tree of everything the +fit wrote is `__Output Folder Layout__` in +[`autogalaxy_workspace/scripts/imaging/modeling.py`](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/modeling.py), +also condensed in +[`wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md) +"The output folder". + +The whole `/` folder is portable: a collaborator with a compatible environment +loads it exactly as above. `files/*.json`, `files/samples.csv`, the `image/*.fits` products +and the human-readable `model.info` / `model.results` are what matter to them. + +## Branch — parameters with errors, from `Samples` + +The `Samples` object holds every accepted sample and its likelihood, which is what makes error +estimation possible at all: + +```python +samples = result.samples + +median_instance = samples.median_pdf() +u3 = samples.values_at_upper_sigma(sigma=3.0) +l3 = samples.values_at_lower_sigma(sigma=3.0) + +print(median_instance.galaxies.galaxy.bulge) +print(u3.galaxies.galaxy.bulge, "\n", l3.galaxies.galaxy.bulge) +``` + +Adapted from `autogalaxy_workspace:scripts/guides/results/start_here.py` `__Samples__`. Each +of these returns a **model instance** — the same shape of object as the model you composed, so +`instance.galaxies.galaxy.bulge.sersic_index` addresses a parameter by the names you chose. +`samples.max_log_likelihood()` gives the best-fit instance, `samples.parameter_lists` the raw +vectors, `samples.log_likelihood_list` and `samples.log_evidence` the statistics. +`samples.errors_at_upper_sigma(...)` returns the error rather than the value. + +Two disciplines when reporting: + +- **Quote the sigma level and the summary you used.** "Effective radius 1.62, +0.08/−0.07 at + 1σ from the marginalised PDF" is a result. "1.62" is not. +- **Be explicit about units.** PyAutoGalaxy works in arcseconds and instrumental flux units. + "Effective radius = 1.6 arcsec" is fine; "effective radius = 1.6" is a bug waiting to + happen. + +### Linear intensities are not in the samples + +This one catches everyone. A linear light profile (`ag.lp_linear.*`, and every basis built +from them — a Multi-Gaussian Expansion, a shapelet expansion) has its `intensity` solved by +linear inversion at every iteration rather than sampled. So in the `Samples` object those +intensities are placeholders of 1.0, and reading them there gives you nothing +(`autogalaxy_workspace:scripts/guides/results/start_here.py` `__Linear Light Profiles__`). + +The solved values live on the fit. Either take them off the result directly: + +```python +print(result.max_log_likelihood_galaxies[0].bulge.intensity) +``` + +or, if you are working from a `Samples` instance you loaded yourself, rebuild the fit and let +it solve: + +```python +ml_instance = samples.max_log_likelihood() + +fit = ag.FitImaging(dataset=dataset, galaxies=ml_instance.galaxies) +galaxies = fit.galaxies_linear_light_profiles_to_light_profiles + +print(galaxies[0].bulge.intensity) +``` + +The same conversion is what lets you plot a linear-profile model at all +([`ag_plot_fit`](./ag_plot_fit.md)). Background: +[`wiki/core/concepts/linear_light_profiles_and_mge.md`](../wiki/core/concepts/linear_light_profiles_and_mge.md). +The same principle applies to a pixelised reconstruction's pixel fluxes — +[`wiki/core/concepts/inversions_and_pixelizations.md`](../wiki/core/concepts/inversions_and_pixelizations.md). + +## Branch — derived quantities, and their errors + +Most interesting numbers are not model parameters. A bulge-to-total light ratio, a luminosity +inside an aperture, an effective radius in kiloparsecs — each is a *function* of the +parameters, so its error has to be propagated rather than read off. + +Compute the quantity from the best-fit galaxies: + +```python +galaxies = result.max_log_likelihood_galaxies + +luminosity = galaxies[0].luminosity_within_circle_from(radius=10.0) + +cosmology = ag.cosmo.Planck15() +kpc_per_arcsec = cosmology.kpc_per_arcsec_from(redshift=galaxies[0].redshift) +effective_radius_kpc = galaxies[0].bulge.effective_radius * kpc_per_arcsec +``` + +Adapted from `autogalaxy_workspace:scripts/guides/results/aggregator/models.py` and +`autogalaxy_workspace:scripts/guides/results/start_here.py` +`__Units and Cosmological Quantities__`. The angular-to-physical machinery, fluxes, magnitudes +and luminosities are +[`wiki/core/concepts/cosmology_and_units.md`](../wiki/core/concepts/cosmology_and_units.md); +the runnable guides are `autogalaxy_workspace:scripts/guides/units/cosmology.py` and +`autogalaxy_workspace:scripts/guides/units/flux.py`. + +For the **error** on a derived quantity, recompute it for galaxies drawn from the posterior +and take the spread. The aggregator's `randomly_drawn_via_pdf_gen_from` exists for exactly +this, and `all_above_weight_gen_from` / `weights_above_gen_from` give you the weighted version +(`autogalaxy_workspace:scripts/guides/results/aggregator/models.py`). Anything else — scaling +the parameter's own error, or quoting the best-fit value with no error — is guesswork. + +Quantities the library can compute but the model never sampled can also be recorded during the +fit as **latent variables**, which puts them in the samples with errors like any other +parameter. That machinery, including subclassing `ag.LatentGalaxy` for your own, is +`autogalaxy_workspace:scripts/guides/results/latent_variables.py` and +[`wiki/core/concepts/samples_and_posteriors.md`](../wiki/core/concepts/samples_and_posteriors.md) +"Latent variables". Note that any active `PYAUTO_TEST_MODE` level skips the latent pass, so a +smoke run never has them. + +## Branch — many fits, with the aggregator + +Once the unit of work is a *sample* — a hundred galaxies with one model, or one galaxy with +ten models — the aggregator gives you one query surface over all of them and, crucially, +returns **generators** rather than lists, so memory stays bounded: + +```python +from pathlib import Path +from autofit.aggregator.aggregator import Aggregator + +agg = Aggregator.from_directory(directory=Path("output") / "imaging") + +for samples in agg.values("samples"): + print(samples.parameter_lists[0]) +``` + +Adapted from `autogalaxy_workspace:scripts/guides/results/start_here.py` +"Aggregator". Note the explicit import: `Aggregator.from_directory` lives in +`PyAutoFit:autofit/aggregator/aggregator.py` and is **not** `af.Aggregator`, which is the +database-backed class. Keys worth knowing for `agg.values(...)`: `"samples"`, +`"samples_summary"` (much faster — nothing re-derived), `"samples_info"`, `"model"`, +`"search"`, `"info"`, `"covariance"`, `"cosmology"`, and the `"dataset/..."` family. + +A generator is consumed once. Remake it rather than storing it — this is why the workspace +examples build each one at the point of use. + +**Queries** compose, because each returns a new aggregator +(`autogalaxy_workspace:scripts/guides/results/aggregator/queries.py`): + +```python +unique_tag = agg.search.unique_tag +agg_query = agg.query(unique_tag == "") + +bulge = agg.model.galaxies.galaxy.bulge +agg_query = agg.query((bulge == ag.lp_linear.Sersic) & (bulge.effective_radius > 3.0)) +``` + +`agg.search` addresses fit-level fields (`name`, `unique_tag`, `path_prefix`, `is_complete`); +`agg.model` addresses the fitted model by the names you chose, so you can select on a +component's class or on an inferred value. + +### `ag.agg` — galaxy objects back out of stored fits + +PyAutoFit's aggregator loads what the search wrote; PyAutoGalaxy's `ag.agg` rebuilds +*galaxy-domain* objects from those files, as generators. The classes are `ImagingAgg`, +`InterferometerAgg`, `GalaxiesAgg`, `FitImagingAgg`, `FitInterferometerAgg`, `EllipsesAgg`, +`FitEllipseAgg` and `MultipolesAgg`: + +```python +import autogalaxy.plot as aplt + +galaxies_agg = ag.agg.GalaxiesAgg(aggregator=agg) +galaxies_gen = galaxies_agg.max_log_likelihood_gen_from() + +dataset_agg = ag.agg.ImagingAgg(aggregator=agg) +dataset_gen = dataset_agg.dataset_gen_from() + +for dataset_list, galaxies_list in zip(dataset_gen, galaxies_gen): + dataset = dataset_list[0] # one entry per analysis in the fit + galaxies = galaxies_list[0] + + fit = ag.FitImaging(dataset=dataset, galaxies=galaxies) + galaxies = fit.galaxies_linear_light_profiles_to_light_profiles + + print(galaxies[0].bulge.intensity) +``` + +Adapted from `autogalaxy_workspace:scripts/guides/results/aggregator/models.py`. The inner +`[0]` is not noise: a fit can hold several analyses (a multi-band factor graph, for instance), +so each generator yields a *list* per fit. The dataset classes expose `dataset_gen_from()`; +the model classes share `max_log_likelihood_gen_from()`, +`randomly_drawn_via_pdf_gen_from(total_samples=...)`, +`all_above_weight_gen_from(minimum_weight=...)` and +`weights_above_gen_from(minimum_weight=...)`, all built on `object_via_gen_from(...)`. +`FitImagingAgg` also takes a `settings=` argument if you need to rebuild the fit with +different settings than it ran with +(`autogalaxy_workspace:scripts/guides/results/aggregator/data_fitting.py`). + +### The three enumerations, and the workflow exports + +`ag.agg` exposes exactly three enumerations, naming panels and HDUs the completed fits +**already have on disk** — you are extracting stored output, not re-rendering it +(`PyAutoGalaxy:autogalaxy/aggregator/subplot.py`): + +| Enumeration | Members | +|---|---| +| `ag.agg.fits_fit` | `model_data`, `residual_map`, `normalized_residual_map`, `chi_squared_map` | +| `ag.agg.subplot_dataset` | `data`, `data_log_10`, `noise_map`, `psf`, `psf_log_10`, `signal_to_noise_map`, `over_sample_size_lp`, `over_sample_size_pixelization` | +| `ag.agg.subplot_fit` | `data`, `signal_to_noise_map`, `model_data`, `normalized_residual_map`, `normalized_residual_map_one_sigma`, `chi_squared_map` | + +They feed three PyAutoFit export helpers, each taking `aggregator=agg` +(`autogalaxy_workspace:scripts/guides/results/workflow/csv_make.py`, +`fits_make.py` and `png_make.py`): + +```python +import autofit as af + +workflow_path = Path("output") / "workflow" + +agg_csv = af.AggregateCSV(aggregator=agg) +agg_csv.add_variable( + argument="galaxies.galaxy.bulge.sersic_index", + value_types=[af.ValueType.Median, af.ValueType.ValuesAt3Sigma], +) +agg_csv.add_variable(argument="galaxies.galaxy.bulge.effective_radius") +agg_csv.save(path=workflow_path / "results.csv") + +agg_fits = af.AggregateFITS(aggregator=agg) +agg_fits.output_to_folder( + folder=workflow_path, + name="unique_tag", + hdus=[ag.agg.fits_fit.model_data, ag.agg.fits_fit.residual_map], +) + +agg_image = af.AggregateImages(aggregator=agg) +agg_image.output_to_folder( + folder=workflow_path, + name="unique_tag", + subplots=[ + ag.agg.subplot_fit.data, + ag.agg.subplot_fit.model_data, + ag.agg.subplot_fit.normalized_residual_map, + ], +) +``` + +`af.ValueType` has four members — `Median`, `MaxLogLikelihood`, `ValuesAt1Sigma`, +`ValuesAt3Sigma`. `AggregateCSV` also offers `add_computed_column` (a function of the loaded +objects, for a derived quantity) and `add_label_column` (an explicit per-fit label, e.g. the +dataset names from `[search.unique_tag for search in agg.values("search")]`). `extract_fits` +and `extract_image` return one combined object across all fits; `output_to_folder` writes one +file per fit instead, which is usually easier to click through when reviewing a sample. + +Full API surface, including the database-backed `af.Aggregator.from_database` route for +hundreds of fits: [`wiki/core/api/aggregator.md`](../wiki/core/api/aggregator.md). The +database walkthrough is `autogalaxy_workspace:scripts/guides/results/database/start_here.py` +— note that `from_database` has no overwrite flag, which is why that example deletes a stale +`.sqlite` itself. Start with the directory aggregator and move to the database only when +reading from disk becomes the slow part. + +The shape that makes this worth the setup: **query** the subset you care about, **load or +derive** the quantity per fit, **export** it under a uniform naming scheme. Fitting is +expensive; changing your mind about which columns you wanted should cost seconds. + +## Branch — the honest caveats + +Say these out loud rather than letting a user discover them: + +- **A test-mode result is not a measurement.** Any `PYAUTO_TEST_MODE` level truncates or mocks + the sampler, so the parameters are whatever the short-circuited run produced, the errors are + degenerate, and latent variables were skipped. The `test_mode` path segment is the tell. +- **A gradient-optimiser fit has no errors.** `MultiStartProdigy` returns one best-fit model. + If the science needs uncertainties, re-run with `Nautilus` + ([`ag_configure_search`](./ag_configure_search.md)). +- **A resumed fit may not be the fit you think.** The unique identifier does not hash the + pixel values, so a result can be a silent resume of an earlier run on different data — see + [`ag_run_search`](./ag_run_search.md) "Resuming". +- **An unconverged search still produces a `samples.csv`.** Check `search.summary` and the + log-evidence before trusting a posterior width. + +## Combine + +- [`ag_plot_fit`](./ag_plot_fit.md) — turn the loaded fit into figures, and read the residuals + properly before quoting anything. +- [`ag_debug_fit_failure`](./ag_debug_fit_failure.md) — the numbers loaded fine but they are + not believable. +- [`ag_run_search`](./ag_run_search.md) — the run that produced this output, and what it wrote + where. +- [`ag_configure_search`](./ag_configure_search.md) — the posterior is too coarse, or you need + errors an optimiser cannot give. + +There is a second, lighter route to a result planned for this assistant: an MCP tool surface +for inspecting a fit's output folder without writing a script. It is **not built yet** — +[`../PENDING.md`](../PENDING.md) tracks it with its grounding module. Until it lands, this +skill's Python is the way in; do not tell a user to invoke a tool that does not exist. + +When a load turns into real analysis — a table, a derived quantity, a comparison across fits — +offer (default-yes) the dated `wiki/project/YYYY-MM-DD-.md` entry: the science question, +what was inferred and how, and the script produced (per [`_style.md`](./_style.md) property +#5). + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Results](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_7_results.ipynb): + the `Result` object from first principles — best-fit model, galaxies, posterior samples, and + pulling parameter values with uncertainties out of them. +- **General reference** — [RTD: Fitting API reference](https://pyautogalaxy.readthedocs.io/en/latest/api/fitting.html): + the generated reference for the fit objects a loaded result rebuilds. +- **Experienced PyAutoGalaxy user** — [workspace: guides/results/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/guides/results/start_here.py): + both routes in one runnable script — direct JSON/CSV loading of one fit, then the aggregator + over many. + +## Agent procedural checklist + +1. Establish what the user wants out of the result before loading anything. +2. Establish one fit or many; pick simple loading or the aggregator accordingly. +3. Check the path for a `test_mode` segment and say so if it is there. +4. Check which search ran — no errors are available from a gradient optimiser. +5. For any linear light profile or pixelisation, get intensities from the fit, never from the + samples. +6. Quote a sigma level, a summary statistic and a unit with every number. +7. For a derived quantity's error, draw from the posterior; never scale a parameter error. +8. Write the loading recipe to a script in `scripts/` rather than leaving it inline. +9. Route figures to the fit-plotting skill and unbelievable numbers to the debug skill. +10. Offer the `wiki/project/` entry when the load became analysis. diff --git a/skills/ag_plot_fit.md b/skills/ag_plot_fit.md new file mode 100644 index 0000000..e23bdca --- /dev/null +++ b/skills/ag_plot_fit.md @@ -0,0 +1,403 @@ +--- +name: ag_plot_fit +description: Visualise a dataset, a galaxy, or a fit with the functional `aplt` plotting API — dataset and fit subplots, individual model-image / residual / normalised-residual / chi-squared panels, per-galaxy breakdowns, log10 stretch and fixed colour limits, overlays, writing figures to disk and writing FITS. Also carries the residual-inspection discipline: what a good chi-squared map looks like, what each characteristic residual pattern is telling you about the model, and how to compare two fits honestly. Use when the user wants to look at data before fitting it, look at a finished fit, produce a figure for a paper or a talk, or decide whether a fit is acceptable. Not for composing or running the fit (`ag_build_imaging_model`, `ag_run_search`), not for pulling numbers and posteriors out of a result (`ag_load_results`), and not for the diagnostic workflow when a fit has clearly failed (`ag_debug_fit_failure`). +--- + +# Looking at the data, the galaxy, and the fit + +Plotting is not decoration here — it is the main instrument of judgement. A galaxy model is +accepted or rejected on whether its residuals look like noise, and no summary statistic +replaces that look. A chi-squared of 1.05 per pixel with a coherent ring of residuals at one +arcsecond is a *worse* result than a chi-squared of 1.3 scattered randomly, because the first +one is telling you the model is structurally wrong while the second is telling you the noise +estimate is slightly off. + +The plotting surface is **functional only**. One import, then module-level functions: + +```python +import autogalaxy.plot as aplt +``` + +Quantities are computed from PyAutoGalaxy objects by their own methods, and the resulting +array or grid is passed to a function +(`autogalaxy_workspace:scripts/guides/plot/start_here.py`). Because the layer takes plain +arrays, anything the library can compute you can plot, and there is no per-object plotting +class to learn. + +There are **no plotter classes and no figure-configuration objects** — the object-oriented +plotters and the matplotlib-wrapper objects that older releases shipped were removed. Those +older releases are heavily represented in language-model training data, which makes +reconstructing them from memory the single most common stale-API error in this library. If a +call you remember is not in `dir(aplt)`, it is not part of the current API: check before you +emit it. [`wiki/core/api/plotting.md`](../wiki/core/api/plotting.md) is the authoritative +enumeration of what does exist, with a "when to use which" note on each entry. + +## Ask + +- *"Are we looking at the data, a galaxy model, or a fit?"* Three different branches below. + Data before fitting is the inspection gate; a galaxy on its own is usually a sanity check on + a model you just composed; a fit is the judgement call. +- *"On screen, or to a file?"* With no `output_path`, a figure is displayed. From a terminal + that flashes and vanishes, so default to writing files. +- *"Quick look, or a figure for a paper?"* A quick look is one subplot. A paper figure means + fixed colour limits, a log10 stretch, a chosen colormap and probably several formats. + +## Branch — the dataset, before you fit it + +One call gives the whole dataset: data, noise-map, PSF, signal-to-noise map and the +over-sampling panels. + +```python +from pathlib import Path + +import autogalaxy as ag +import autogalaxy.plot as aplt + +DATASET_PATH = Path("dataset") / "imaging" / "" +PLOT_DIR = Path("scripts") / "scratch" / "" +PLOT_DIR.mkdir(parents=True, exist_ok=True) + +dataset = ag.Imaging.from_fits( + data_path=DATASET_PATH / "data.fits", + psf_path=DATASET_PATH / "psf.fits", + noise_map_path=DATASET_PATH / "noise_map.fits", + pixel_scales=0.1, +) + +aplt.subplot_imaging_dataset( + dataset=dataset, + output_path=PLOT_DIR, + output_filename="dataset", + output_format="png", +) + +print(f"Saved to: {PLOT_DIR.resolve()}") +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/plot.py` `__Dataset Subplot__`. Individual +components are attributes, plotted with the fundamental array function: + +```python +aplt.plot_array(array=dataset.data, title="Data") +aplt.plot_array(array=dataset.noise_map, title="Noise Map") +aplt.plot_array(array=dataset.psf.kernel, title="PSF") +aplt.plot_array(array=dataset.signal_to_noise_map, title="Signal-to-Noise Map") +``` + +What to look for in this figure is the substance of +[`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md) — neighbouring galaxies, foreground +stars, reduction artefacts, and where the galaxy's emission meets the sky, which is what sets +the mask radius. Re-plot after `apply_mask` and after `apply_over_sampling`: the subplot's +lower panels update to show the over-sampling scheme, so you can see the scheme you asked for +rather than trusting it (`autogalaxy_workspace:scripts/imaging/modeling.py` +`__Over Sampling__`). Masks and over-sampling as concepts are +[`wiki/core/concepts/grids_and_masks.md`](../wiki/core/concepts/grids_and_masks.md); the +dataset object itself is [`wiki/core/api/datasets.md`](../wiki/core/api/datasets.md). + +## Branch — a galaxy or a light profile on its own + +Useful before any data is involved: does the model you just composed actually look like the +galaxy you have in mind? Every object exposes `image_2d_from(grid=grid)`, and the array goes +straight to `plot_array` +(`autogalaxy_workspace:scripts/guides/plot/plotters.py`): + +```python +grid = ag.Grid2D.uniform(shape_native=(100, 100), pixel_scales=0.05) + +bulge = ag.lp.Sersic( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.9, angle=45.0), + intensity=4.0, + effective_radius=0.6, + sersic_index=3.0, +) +disk = ag.lp.Exponential( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.7, angle=30.0), + intensity=2.0, + effective_radius=1.6, +) + +galaxy = ag.Galaxy(redshift=0.5, bulge=bulge, disk=disk) + +aplt.plot_array(array=bulge.image_2d_from(grid=grid), title="Bulge Image") +aplt.plot_array(array=galaxy.image_2d_from(grid=grid), title="Galaxy Image", use_log10=True) +``` + +A galaxy's image is the sum of its light profiles' images, so the bulge/disk decomposition is +visible by plotting each component beside the total. `aplt.subplot_galaxy_light_profiles` +does that as one figure; `aplt.subplot_galaxies` shows the summed image plus a panel per +galaxy, and takes `auto_filename` to name its file. The physics of what the effective radius +and Sersic index actually measure is +[`wiki/core/concepts/light_profiles.md`](../wiki/core/concepts/light_profiles.md), and how +several profiles and galaxies compose is +[`wiki/core/concepts/galaxies.md`](../wiki/core/concepts/galaxies.md). + +**Linear light profiles cannot be plotted directly.** An `ag.lp_linear.*` profile has no +`intensity` until the linear inversion solves for it, so the galaxy-level subplots reject +them. Get the solved profiles from a fit — `fit.galaxies_linear_light_profiles_to_light_profiles` +— and plot those instead (`autogalaxy_workspace:scripts/guides/results/start_here.py` +`__Linear Light Profiles__`). Background: +[`wiki/core/concepts/linear_light_profiles_and_mge.md`](../wiki/core/concepts/linear_light_profiles_and_mge.md). + +## Branch — the fit, and how to read the residuals + +A `FitImaging` compares a model galaxy against the data, PSF convolution included. Build one +from a dataset and galaxies, or take the best-fit one straight off a result: + +```python +fit = ag.FitImaging(dataset=dataset, galaxies=galaxies) + +# or, from a completed search: +fit = result.max_log_likelihood_fit +``` + +The summary figure is six panels — data, signal-to-noise map, model image, residual map, +normalised residual map and chi-squared map: + +```python +aplt.subplot_fit_imaging( + fit=fit, + output_path=PLOT_DIR, + output_format="png", +) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/plot.py` `__Fit Subplot__`. + +**The fit subplots do not take `output_filename`.** They write a fixed stem into +`output_path`: `subplot_fit_imaging` writes `fit.png`, and +`subplot_fit_imaging_of_galaxy` writes `of_galaxy_.png` +(`PyAutoGalaxy:autogalaxy/imaging/plot/fit_imaging_plots.py`). So the *directory* is what +distinguishes one context from another — give each fit, or each variant of a figure, its own +folder under `scripts/scratch/`. `plot_array`, `plot_grid` and `subplot_imaging_dataset` do +take `output_filename`; `subplot_galaxies` takes `auto_filename` +(`PyAutoGalaxy:autogalaxy/galaxy/plot/galaxies_plots.py`). Check the signature rather than +assuming a uniform interface. + +Individual quantities are attributes, so they go through `plot_array` — which is also how you +get per-panel control the subplot does not offer: + +```python +aplt.plot_array(array=fit.model_data, title="Model Image") +aplt.plot_array(array=fit.residual_map, title="Residual Map") +aplt.plot_array(array=fit.normalized_residual_map, title="Normalized Residual Map", symmetric=True) +aplt.plot_array(array=fit.chi_squared_map, title="Chi-Squared Map") +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/plot.py` `__Fit Figures__`. `symmetric=True` +puts zero at the centre of a diverging colour scale, which is the right choice for a residual +map — otherwise a colormap keyed to the extremes hides the sign of the structure you are +looking for. + +### The inspection discipline + +Three maps, three different questions: + +- **`residual_map`** = data − model, in data units. Answers *how much* flux is unaccounted + for, and where. Useful for judging whether a mismatch matters physically. +- **`normalized_residual_map`** = residual / noise. Answers *how significant* the mismatch is. + This is the map to read first, because it is the only one calibrated against the noise. + Values scattered within roughly ±3 are consistent with noise; coherent structure well + outside that is a model failure, however small it looks in flux. +- **`chi_squared_map`** = normalised residual squared. Answers *where the likelihood is being + paid*, since `fit.chi_squared` is its sum + (`autogalaxy_workspace:scripts/imaging/likelihood_function.py` `__Chi Squared__`). Squaring + discards the sign, so read it alongside the normalised residuals rather than instead of + them. + +Read the pattern, not just the magnitude. The characteristic ones: + +| Pattern in the normalised residuals | Usually means | +|---|---| +| Structure concentrated at the very centre | Over-sampling too coarse, or a genuinely cored/nuclear component the model lacks | +| A symmetric ring at one radius | The radial profile is wrong — one Sersic where a bulge plus disk is needed, or a fixed Sersic index | +| Four-lobed pattern, alternating sign | The ellipticity or position angle is wrong, or the isophotes twist with radius | +| A compact blob away from the centre | A neighbouring galaxy, a foreground star, or a reduction artefact left in the fit | +| Clumpy structure with no symmetry | Real asymmetric morphology no smooth profile can absorb | +| Uniformly slightly too large everywhere | The noise-map is underestimated, or a sky background is unmodelled | + +Each row routes somewhere. Ellipticity and isophote twists are +[`wiki/core/concepts/ellipse_fitting_and_multipoles.md`](../wiki/core/concepts/ellipse_fitting_and_multipoles.md); +neighbours are +[`wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md); +irreducible clumpiness is +[`wiki/core/concepts/inversions_and_pixelizations.md`](../wiki/core/concepts/inversions_and_pixelizations.md); +an unmodelled background is +[`wiki/core/concepts/sky_background_and_operated_profiles.md`](../wiki/core/concepts/sky_background_and_operated_profiles.md). +Turning any of those diagnoses into a plan is +[`ag_debug_fit_failure`](./ag_debug_fit_failure.md). + +For a blended pair or a galaxy with a modelled neighbour, the per-galaxy view is what shows +whether the decomposition is real or just two profiles trading flux: + +```python +aplt.plot_array(array=fit.model_images_of_galaxies_list[0], title="Galaxy 0 Model Image") +aplt.subplot_fit_imaging_of_galaxy(fit=fit, galaxy_index=0, output_path=PLOT_DIR, output_format="png") +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/plot.py` `__Galaxy Images__` and +`__Galaxy Subplots__`. The per-galaxy subplot shows the data with the *other* galaxies +subtracted, alongside that galaxy's model image and residuals. + +## Branch — making the figure say what you mean + +Customisation is keyword arguments, never a wrapper object +(`autogalaxy_workspace:scripts/guides/plot/start_here.py` `__Customization__`): + +```python +aplt.plot_array(array=dataset.data, title="Jet", colormap="jet") +aplt.plot_array(array=galaxy.image_2d_from(grid=grid), title="Log10", use_log10=True) +aplt.plot_array(array=dataset.data, title="Fixed Scale", vmin=0.0, vmax=1.0) +``` + +Three of these carry scientific weight: + +- **`use_log10=True`** — a galaxy's surface brightness spans orders of magnitude, so a linear + stretch shows the bright core and nothing else. The log10 view is where faint extended + emission, tidal features and a disk under a bright bulge become visible. It is available on + the subplots too, and is often the more honest default for a galaxy. +- **`vmin` / `vmax`** — fixed limits are what make two figures comparable. Comparing two fits + on auto-scaled colour bars is not a comparison. +- **`symmetric=True`** — for anything that can be negative, which is every residual map. + +Overlays are keyword arguments as well +(`autogalaxy_workspace:scripts/guides/plot/visuals.py`): `positions=` scatters points such as +light-profile centres, `grid=` draws a coordinate grid, `lines=` / `line_colors=` draw +polylines. Two accepted forms for positions: + +```python +import numpy as np + +light_profile_centres = galaxies.extract_attribute(cls=ag.LightProfile, attr_name="centre") + +aplt.plot_array( + array=galaxies.image_2d_from(grid=grid), + positions=[np.array(light_profile_centres)], + title="Image with Light Profile Centres", +) + +aplt.plot_array( + array=galaxy.image_2d_from(grid=grid), + positions=ag.Grid2DIrregular(values=[(0.0, 0.0)]), + title="Image with Centre", +) +``` + +1D radial profiles are plain matplotlib on a radially projected grid, not an `aplt` function — +including the shaded error region drawn from a set of profiles sampled from the posterior +(`autogalaxy_workspace:scripts/guides/plot/plotters.py` `__One Dimensional Plots__` and +`__Probability Density Function (PDF) Plots__`, which uses +`ag.util.error.profile_1d_median_and_error_region_via_quantile`). + +Anything not passed explicitly comes from `config/visualize/general.yaml` — the default +colormap, tick counts, colorbar label sizes, contour levels, the unit label. Editing that file +changes the look project-wide without touching code, after a session or kernel restart. See +[`wiki/core/api/configuration.md`](../wiki/core/api/configuration.md). + +## Branch — writing to disk, and announcing it + +With `output_path` set (a directory) plus `output_format`, the figure is written rather than +displayed. Without `output_filename` the file is named from the title; `output_format` also +accepts a list, which writes several formats in one call — a PNG to look at and a PDF for the +paper (`autogalaxy_workspace:scripts/guides/plot/start_here.py` `__Output__`): + +```python +aplt.plot_array( + array=dataset.data, + title="Image", + output_path=PLOT_DIR, + output_filename="data", + output_format=["png", "pdf"], +) +``` + +Three rules from [`_style.md`](./_style.md) "Plot output and path announcement", and they are +not optional: + +1. Throwaway figures go to the gitignored `scripts/scratch//`, where `` is + usually the dataset name. Never `output/` — that belongs to the fit runtime. +2. `print(f"Saved to: {PLOT_DIR.resolve()}")` at the end of the recipe, so the absolute + location lands in stdout. +3. After running, **quote that absolute path back to the user and offer to open it** — + `xdg-open` on Linux, `open` on macOS, `explorer.exe` or `wslview` from WSL. One offer per + plot run. "Plot saved" on its own is not useful to someone who cannot see your filesystem. + +To write data rather than an image, the `fits_*` functions are separate +(`autogalaxy_workspace:scripts/imaging/plot.py` `__Outputting to FITS__`): + +```python +aplt.fits_imaging(dataset=dataset, file_path=PLOT_DIR / "dataset.fits", overwrite=True) +aplt.fits_array(array=mask, file_path=DATASET_PATH / "mask.fits", overwrite=True) +``` + +`file_path` writes one multi-HDU file with named extensions; the per-component +`data_path` / `psf_path` / `noise_map_path` arguments write separate files instead. + +## Branch — the figures a fit produces on its own + +You do not need any of the above to watch a fit converge. Each analysis carries a `Visualizer` +that PyAutoFit invokes at every update interval and again at the end, writing the standard +dataset, fit and per-galaxy figures into the fit's `image/` folder using the best model found +so far (`autogalaxy_workspace:scripts/imaging/plot.py` `__Visualizer__`; +`PyAutoGalaxy:autogalaxy/imaging/model/visualizer.py`). Which figures appear is set by +`config/visualize/plots.yaml` — `dataset -> subplot_dataset`, `fit -> subplot_fit`, +`fit -> subplot_of_galaxies` — so switching a figure on or off for every fit is a YAML edit, +not a code change. The cadence is the search's `iterations_per_quick_update` +([`ag_run_search`](./ag_run_search.md)). + +Reach for this skill when you want something the automatic figures do not give you: a +different stretch, a crop, fixed limits across two fits, a specific panel at publication +size, or a per-galaxy breakdown that `plots.yaml` has switched off. + +Posterior figures — corner plots and likelihood traces — come from PyAutoFit and take +`samples` plus `path` / `filename` / `format` rather than the `output_*` names, e.g. +`aplt.corner_cornerpy(samples=result.samples)`. They belong to +[`ag_load_results`](./ag_load_results.md); the surface is tabulated in +[`wiki/core/api/plotting.md`](../wiki/core/api/plotting.md) "Search results". + +## Combine + +- [`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md) — the inspection the dataset + subplot exists to serve, and the mask that follows from it. +- [`ag_debug_fit_failure`](./ag_debug_fit_failure.md) — you have read the residuals and they + are wrong; this is the triage. +- [`ag_load_results`](./ag_load_results.md) — the numbers behind the figures, and figures + across many fits at once. +- [`ag_run_search`](./ag_run_search.md) — the fit that produced the `image/` folder, and the + cadence at which it refreshes. +- [`ag_simulate_dataset`](./ag_simulate_dataset.md) — plot a simulated fit where you know the + truth, to calibrate what an acceptable residual map looks like for your data quality. + +When a figure is going into a paper or a talk, offer (default-yes) the dated +`wiki/project/YYYY-MM-DD-.md` entry recording what it shows and which fit it came from +— per [`_style.md`](./_style.md) property #5. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Fitting](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_1_introduction/tutorial_3_fitting.ipynb): + builds up residuals, normalised residuals and the chi-squared map from first principles, so + the six panels of a fit subplot stop being opaque. +- **General reference** — [RTD: Plot API reference](https://pyautogalaxy.readthedocs.io/en/latest/api/plot.html): + the generated reference for the visualisation library. +- **Experienced PyAutoGalaxy user** — [workspace: guides/plot/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/guides/plot/start_here.py): + the plotting API in one runnable script — customisation, output, config defaults, overlays + and subplots. + +## Agent procedural checklist + +1. Establish what is being plotted: dataset, galaxy/profile, or fit. +2. Ground every call against [`wiki/core/api/plotting.md`](../wiki/core/api/plotting.md) or + `dir(aplt)`. Never write a plotter class or a figure-configuration object. +3. Check the signature for `output_filename` before passing it — the fit and galaxy subplots + do not accept it and write a fixed stem instead. +4. Write the recipe into a script; send figures to `scripts/scratch//`. +5. `print` the resolved plot directory at the end of the script. +6. For a fit, read the **normalised residual** map first, then the chi-squared map; report the + pattern, not only the number. +7. Use `symmetric=True` on residuals, `use_log10=True` for galaxy light, and fixed + `vmin`/`vmax` whenever two figures will be compared. +8. Quote the absolute path back and offer to open it — once per plot run. +9. Route any diagnosed model failure to the debug skill rather than tweaking the figure until + it looks acceptable. +10. Offer the `wiki/project/` entry for a figure that will be published. diff --git a/skills/ag_prepare_imaging_data.md b/skills/ag_prepare_imaging_data.md new file mode 100644 index 0000000..f1d4fcc --- /dev/null +++ b/skills/ag_prepare_imaging_data.md @@ -0,0 +1,484 @@ +--- +name: ag_prepare_imaging_data +description: Load a user's own CCD imaging of a galaxy from FITS into an `ag.Imaging` dataset and get it ready to fit — pixel scale, flux units, the RMS noise-map, the PSF, the mask extent, contaminating extra galaxies and foreground stars, over-sampling, and the optional `info.json` / centre sidecars. **This skill owns the real-data inspection gate**: the plot-and-settle procedure that must run before any model is composed or fitted on real observational data. Use whenever the user brings their own data (HST, JWST, Euclid, ground-based), whenever a dataset needs checking against PyAutoGalaxy's standards, or whenever a mask or contaminant decision is open. Not for simulating data (`ag_simulate_dataset`), not for composing the model itself (`ag_build_imaging_model`), and not for visibility-plane data. +--- + +# Preparing your own imaging data + +Everything downstream inherits the decisions made here. A noise-map that is a variance +rather than an RMS makes every likelihood wrong by a constant power; a mask that truncates +the outer isophotes biases the effective radius and Sersic index directly; an unmodelled +neighbour inside the mask pulls the fit toward a galaxy that isn't there. None of those +three announce themselves — the search converges, the residuals look plausible, and the +number you publish is wrong. That is why this step is a gate rather than a formality. + +Scientifically, the job is to hand the likelihood three things it can trust: an **image** +in electrons per second, a **per-pixel RMS noise-map** in the same units, and a **PSF** to +forward-convolve the model with. Statistically, the job is to define the *support* of the +likelihood — which pixels are in the sum — and to make sure every pixel inside it is +described by the model you are about to fit. The standards themselves come from +`autogalaxy_workspace:scripts/imaging/data_preparation/start_here.py`; the API catalogue +is [`../wiki/core/api/datasets.md`](../wiki/core/api/datasets.md). + +## The real-data gate — non-negotiable + +**Before composing or running any model-fit on real observational data**, plot it, show the +user the `dataset.png` path, and settle two things from that same look: + +**(a) Extra galaxies, foreground stars and artefacts.** These are the single largest source +of fit bias. Decide, explicitly, which of the three strategies below applies to each one. + +**(b) The mask extent** — the radius and shape that captures the galaxy's emission out to +where it meets the sky, without dragging in noise or contaminants. **Never leave the mask +radius as a silent default on real data.** A mask that cuts inside the outer isophotes +biases `effective_radius` and `sersic_index`, so this is a science-critical choice, not a +tidy-up. + +**If you cannot plot it yourself** — no code execution, for example a GitHub-connector chat +— **the gate is not waived**: ask the user to plot and inspect the data, and to confirm both +(a) contaminants and (b) the mask extent, before you compose the fit. These are the +questions every real-data run must ask, on every harness. + +**Simulated data is exempt** — you know the truth, including where every component is, so +there is nothing to discover by looking. Everything else in this skill still applies to +simulated data as documentation of what a real run would need. + +The gate is satisfied by looking, not by asserting. The first branch below is how you look. + +## Ask + +- *"What instrument and filter, and what is the pixel scale?"* — you cannot proceed without + it, and it is the one number nothing in the data can tell you. +- *"Is the noise-map an RMS standard deviation per pixel, or a variance / weight map / HST + WHT map?"* — if they are not certain, treat it as unknown and convert; the likelihood + assumes RMS. +- *"Are there neighbouring galaxies, foreground stars or reduction artefacts near the + target?"* — the (a) half of the gate. +- *"How far out do you care about the light — just the bright body, or the faint outer + envelope?"* — this is the (b) half of the gate phrased as a science question, which is + how it should be decided. + +## Branch — look at it first + +Load the three ingredients and plot them. Nothing else happens until this figure has been +seen. + +```python +""" +Data Inspection: () +============================================= + +Load a galaxy's CCD imaging from FITS and inspect it before any model is composed. This +script exists to satisfy the real-data inspection gate: the image, noise-map, PSF and +signal-to-noise map are plotted so that two science-critical decisions can be made from +evidence rather than defaults — which contaminating objects are present, and how far out +the mask should extend. + +__Contents__ + +- **Imports:** Import the required libraries. +- **Dataset:** Load the image, noise-map and PSF from FITS with the correct pixel scale. +- **Standards Check:** Verify the PSF is odd-sized and normalised, and report the S/N. +- **Plot:** Write the dataset subplot to disk for inspection. +""" + +""" +__Imports__ +""" +from pathlib import Path + +import autogalaxy as ag +import autogalaxy.plot as aplt + +DATASET_PATH = Path("dataset") / "imaging" / "my_galaxy" +PLOT_DIR = Path("scripts") / "scratch" / "my_galaxy" + +""" +__Dataset__ + +Three ingredients are needed for galaxy modeling: the image in electrons per second, a +per-pixel RMS noise-map in the same units, and the PSF describing the blurring imposed by +the telescope optics. The PSF is *forward-convolved* onto the model rather than divided out +of the data, which is what lets a fit distinguish a genuinely compact bulge from a +seeing-broadened one. + +`pixel_scales` converts pixels to arcseconds and must be correct for your instrument — +HST/ACS ~0.05", JWST/NIRCam ~0.03-0.06", Euclid VIS 0.1" and NISP 0.2", ground-based +0.2-0.3". Nothing in the FITS data can correct a wrong value for you; every inferred size +scales with it. Loading is handled by `ag.Imaging.from_fits` +(`PyAutoArray:autoarray/dataset/imaging/dataset.py`). +""" +PIXEL_SCALES = 0.1 + +dataset = ag.Imaging.from_fits( + data_path=DATASET_PATH / "data.fits", + noise_map_path=DATASET_PATH / "noise_map.fits", + psf_path=DATASET_PATH / "psf.fits", + pixel_scales=PIXEL_SCALES, +) + +""" +__Standards Check__ + +Two properties of the PSF are cheap to check and expensive to get wrong. An **even-sized** +kernel shifts the convolved model by half a pixel, which propagates straight into the +inferred centre; and a kernel that does not **sum to unity** does not conserve flux, which +corrupts any magnitude or luminosity derived from the fit. The signal-to-noise peak is a +sanity check on the noise-map: a galaxy usually peaks somewhere around 10-300, and a value +orders of magnitude away from that is a sign the noise-map is in the wrong units or is a +variance rather than an RMS. +""" +psf_shape = dataset.psf.kernel.shape_native + +print(f"PSF shape: {psf_shape}, odd: {all(s % 2 == 1 for s in psf_shape)}") +print(f"PSF sum: {float(dataset.psf.kernel.sum()):.6f} (should be 1.0)") +print(f"peak signal-to-noise: {float(dataset.signal_to_noise_map.max()):.1f}") + +""" +__Plot__ + +The subplot shows the image, noise-map, PSF and signal-to-noise map together, which is +exactly the set needed to answer the two gate questions: what contaminating objects are +present, and how far out the galaxy's emission extends before it meets the sky. +""" +aplt.subplot_imaging_dataset( + dataset=dataset, + output_path=PLOT_DIR, + output_filename="dataset", + output_format="png", +) + +print(f"Saved to: {PLOT_DIR.resolve()}") +``` + +Then **quote the absolute path of `dataset.png` back to the user and offer to open it** +(`xdg-open` on Linux, `open` on macOS, `explorer.exe` or `wslview` from WSL) — one offer, +not repeated nagging. Ask the two gate questions against that specific figure, naming what +you can see in it. "There's a compact object about 2" north-east of the centre — is that a +neighbour or part of the galaxy?" is a useful question; "are there any contaminants?" in the +abstract is not. + +## Branch — the three standards, and how to fix a violation + +`autogalaxy_workspace:scripts/imaging/data_preparation/start_here.py` is the checklist; the +per-ingredient tools are in `examples/data.py`, `examples/noise_map.py` and +`examples/psf.py` beside it. + +**The image.** Flux in **electrons per second**, galaxy near the centre, cut down to a +postage stamp. The units matter because default priors on light-profile `intensity` assume +them, and because magnitudes are computed from them. Conversions both ways are in the +preprocess module (`PyAutoArray:autoarray/dataset/preprocess.py`): + +```python +exposure_time_map = ag.Array2D.full( + fill_value=1000.0, + shape_native=data.shape_native, + pixel_scales=data.pixel_scales, +) + +data_counts = ag.preprocess.array_eps_to_counts( + array_eps=data, exposure_time_map=exposure_time_map +) +data_eps = ag.preprocess.array_counts_to_eps( + array_counts=data_counts, exposure_time_map=exposure_time_map +) +``` + +ADUs need the gain as well — `ag.preprocess.array_eps_to_adus(array_eps=data, gain=4.0, +exposure_time_map=exposure_time_map)` and `array_adus_to_eps` back. If your reduction +produced a real per-pixel exposure-time map, load it with `ag.Array2D.from_fits` instead of +the flat `full` above; a flat map is a good approximation for many HST observations and not +for all. Trim an oversized stamp with +`ag.preprocess.array_with_new_shape(array=data, new_shape=(80, 80))`, which crops centred. +Background subtraction is better done by your reduction pipeline, but +`background_sky_level_via_edges_from` and `background_noise_map_via_edges_from` exist for +when it wasn't — and modelling the residual sky as a free parameter is often the better +answer anyway (see the sky-background note in +[`ag_build_imaging_model`](./ag_build_imaging_model.md)). + +**The noise-map.** RMS standard deviation per pixel, in electrons per second, **including** +the Poisson contribution from the galaxy's own counts as well as background sky. This is the +one to be most careful about: reduction pipelines frequently drop the Poisson term, and a +noise-map that is a variance, an inverse variance, or an HST WHT map will produce a fit that +converges confidently to the wrong answer. The preprocess module has a converter for each +common input form — `noise_map_via_weight_map_from`, +`noise_map_via_inverse_noise_map_from`, +`noise_map_via_data_eps_and_exposure_time_map_from`, +`noise_map_via_data_eps_exposure_time_map_and_background_noise_map_from`, +`noise_map_via_data_eps_exposure_time_map_and_background_variances_from` and +`poisson_noise_via_data_eps_from`. If you cannot determine which form you have, the +instrument handbook is the authority; guessing here is not recoverable later. + +**The PSF.** Odd dimensions, normalised to unity, centred, and roughly 11×11 to 21×21 +pixels. Larger kernels (51×51) work but slow every likelihood evaluation, since convolution +runs on every model image. Load with normalisation applied: + +```python +psf = ag.Convolver.from_fits( + file_path=DATASET_PATH / "psf.fits", + hdu=0, + pixel_scales=PIXEL_SCALES, + normalize=True, +) +``` + +`ag.Imaging.from_fits` normalises the PSF internally anyway, so this is belt-and-braces. An +even-sized kernel is a genuine problem rather than a warning: +`ag.preprocess.kernel_with_odd_dimensions_from` will interpolate one for you, but it *is* +an interpolation and the right fix is to re-derive the PSF at odd size in your reduction. +Resize with `ag.preprocess.array_with_new_shape(array=psf.kernel, new_shape=(21, 21))`. + +## Branch — (a) contaminants: three strategies, one decision per object + +For each neighbouring galaxy, foreground star or artefact you identified in the figure, +choose one: + +**1. Scale its noise (the default).** Keep the pixels in the fit but blow their noise-map +values up so they contribute negligibly to the likelihood. Preferred over deleting pixels +because it leaves the pixel grid intact — which matters a great deal later if you move to a +pixelised reconstruction, where removed pixels create discontinuities. + +```python +mask_extra_galaxies = ag.Mask2D.from_fits( + file_path=DATASET_PATH / "mask_extra_galaxies.fits", + pixel_scales=dataset.pixel_scales, + invert=True, # `True` marks the pixels whose noise is scaled. +) + +dataset = dataset.apply_noise_scaling(mask=mask_extra_galaxies) + +aplt.subplot_imaging_dataset( + dataset=dataset, output_path=PLOT_DIR, output_filename="dataset_scaled", + output_format="png", +) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/modeling.py`. Note `invert=True`: the +FITS convention here marks the *scaled* region as `True`. Build the mask as a union of +circles when you know the positions, following +`autogalaxy_workspace:scripts/imaging/data_preparation/examples/optional/mask_extra_galaxies.py`: + +```python +import numpy as np + +extra_galaxies_mask = np.zeros(data.shape_native, dtype=bool) + +for centre, radius in [((1.0, 3.5), 1.5), ((-2.0, -3.5), 2.4)]: + circle = ag.Mask2D.circular( + shape_native=data.shape_native, + pixel_scales=data.pixel_scales, + centre=centre, + radius=radius, + invert=True, + ) + extra_galaxies_mask = np.logical_or(extra_galaxies_mask, circle.native) + +mask = ag.Mask2D(mask=extra_galaxies_mask, pixel_scales=data.pixel_scales) + +aplt.fits_array( + array=mask, file_path=DATASET_PATH / "mask_extra_galaxies.fits", overwrite=True +) +``` + +For irregular shapes there is a spray-paint GUI at +`autogalaxy_workspace:scripts/imaging/data_preparation/gui/mask_extra_galaxies.py`, and the +same `ag.Scribbler`-based tool appears inline at the foot of +`autogalaxy_workspace:scripts/imaging/start_here.py`. Point the user at it rather than +hand-deriving centres when the contaminant is a messy blend. + +**2. Shrink the mask so it falls outside.** The simplest option, and the right one when the +contaminant sits comfortably beyond the galaxy's emission. It costs you nothing except the +sky annulus you gave up — but check it does not also cut into the outer isophotes, because +then you have traded a known bias for a worse one. + +**3. Model it as an extra galaxy.** Give it its own light profile with a fixed centre. This +is the choice when its light genuinely overlaps the target's, so neither masking nor scaling +can separate them without also removing signal you need. That is a *model* decision, so it +belongs to [`ag_build_imaging_model`](./ag_build_imaging_model.md); what you produce here is +the list of centres it will consume: + +```python +extra_galaxies_centres = ag.Grid2DIrregular(values=[(1.0, 3.5), (-2.0, -3.5)]) + +aplt.plot_array( + array=data, + title="Data", + positions=[np.array(extra_galaxies_centres)], + output_path=PLOT_DIR, + output_filename="data_with_extra_galaxies", + output_format="png", +) + +ag.output_to_json( + obj=extra_galaxies_centres, + file_path=DATASET_PATH / "extra_galaxies_centres.json", +) +``` + +Adapted from +`autogalaxy_workspace:scripts/imaging/data_preparation/examples/optional/extra_galaxies_centres.py`. +Always plot the centres over the image before saving them — an off-by-a-sign (y,x) ordering +is obvious in a figure and invisible in a JSON file. There is a click-to-select GUI at +`autogalaxy_workspace:scripts/imaging/data_preparation/gui/extra_galaxies_centres.py` which +snaps to the brightest pixel in a 5×5 box. + +The physics of choosing between the three, and what noise scaling does to the likelihood, +is +[`../wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md). + +## Branch — (b) the mask extent + +The mask defines which pixels enter the likelihood. Choose its radius from the figure, and +say out loud what you chose and why: + +```python +MASK_RADIUS = 2.5 + +mask = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=dataset.pixel_scales, + radius=MASK_RADIUS, +) + +dataset = dataset.apply_mask(mask=mask) + +aplt.subplot_imaging_dataset( + dataset=dataset, output_path=PLOT_DIR, output_filename="dataset_masked", + output_format="png", +) +``` + +Two failure directions, and they are not symmetric. **Too small** truncates the outer +isophotes and biases `effective_radius` and `sersic_index` — a real, systematic error in +your science result. **Too large** drags in sky-dominated pixels: it slows the fit and adds +little information, but it does not bias the answer, and it keeps honest sky pixels that +constrain a background level. When in doubt, err large and re-plot; the cost is run time, +not correctness. + +Circular is the default and usually right for a single galaxy. Non-circular shapes exist +when the geometry calls for it — `ag.Mask2D.circular_annular(inner_radius=..., outer_radius=...)` +to exclude a saturated core, `ag.Mask2D.elliptical(major_axis_radius=..., axis_ratio=..., +angle=...)` and `ag.Mask2D.elliptical_annular(...)` for a strongly inclined disk. All four +are shown in +`autogalaxy_workspace:scripts/imaging/data_preparation/examples/optional/mask.py`, and +`autogalaxy_workspace:scripts/guides/modeling/customize.py` covers applying a custom mask to +a fit. A mask drawn by hand in the GUI +(`autogalaxy_workspace:scripts/imaging/data_preparation/gui/mask.py`) is written to +`mask.fits` and loaded back with `ag.Mask2D.from_fits`; for a genuinely irregular footprint, +`autogalaxy_workspace:scripts/imaging/data_preparation/manual/mask_irregular.py` is the +manual route. + +Why the radius is a science decision rather than a default, and what the slim/native +distinction means once a mask is applied, is +[`../wiki/core/concepts/grids_and_masks.md`](../wiki/core/concepts/grids_and_masks.md). + +## Branch — over-sampling + +A Sersic profile's intensity varies steeply across a single central pixel, so evaluating it +once at the pixel centre under-counts the flux. Over-sampling evaluates the profile on a +finer sub-grid and averages, and it changes the inferred parameters — this is accuracy, not +cosmetics. Do it adaptively: high in the centre where the gradient is steep, low in the +outskirts where it is flat and the cost would be wasted. + +```python +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=dataset.grid, + sub_size_list=[8, 4, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0)], +) + +dataset = dataset.apply_over_sampling(over_sample_size_lp=over_sample_size) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/modeling.py`. That reads as 8×8 +sub-sampling inside 0.3", 4×4 between 0.3" and 0.6", and 2×2 beyond — note the list ends at +2, never 1, so even the outskirts keep a floor. `centre_list` takes **every** bright +centre, so if you are modelling an extra galaxy rather than masking it, add its centre here +too; a companion evaluated at 2×2 while the target gets 8×8 is a quiet accuracy loss. A +plain integer (`over_sample_size_lp=4`) applies uniform over-sampling when you want +simplicity over tuning. The `_lp` suffix names the grid used for **l**ight-**p**rofile +evaluation, which is a different grid from a pixelisation's. + +## Branch — the optional sidecars + +None of these are required to fit; all of them make the result easier to interpret later. + +**`info.json`** — auxiliary numbers that travel with the dataset: redshift, velocity +dispersion, a stellar mass from the literature, a previous paper's measurement. Passed to +`search.fit(model=model, analysis=analysis, info=info)`, it is stored with the fit and can +be read back by the aggregator, which is what makes a population-scale comparison against +external measurements possible. + +```python +import json + +info = {"redshift": 0.5, "velocity_dispersion": 250.0, "stellar_mass": 1e11} + +with open(DATASET_PATH / "info.json", "w+") as f: + json.dump(info, f, indent=4) +``` + +Adapted from +`autogalaxy_workspace:scripts/imaging/data_preparation/examples/optional/info.py`. + +**`light_centre.json`** — the galaxy's light centre as an `ag.Grid2DIrregular`, saved with +`ag.output_to_json`, for use as a fixed value in the model. Fixing the centre removes two +free parameters and rules out solutions where the model centre wanders somewhere +unphysical; it is the standard rescue when a search will not converge. Grounded in +`autogalaxy_workspace:scripts/imaging/data_preparation/examples/optional/light_centre.py`, +with a click-to-select GUI at +`autogalaxy_workspace:scripts/imaging/data_preparation/gui/light_centre.py`. + +The on-disk layout these files belong to is not yet a wiki page — until it is, +`autogalaxy_workspace:scripts/imaging/data_preparation/start_here.py` is the ground truth +(the workspace simulates its imaging datasets on first run rather than committing them; +[`../PENDING.md`](../PENDING.md) tracks the page). + +## Combine — where this hands off + +- **The gate is satisfied and the dataset is loaded, masked and over-sampled** → + [`ag_build_imaging_model`](./ag_build_imaging_model.md). Tell it what you decided about + contaminants and the mask; it will ask otherwise. +- **You want to rehearse the whole loop before touching real data** → + [`ag_simulate_dataset`](./ag_simulate_dataset.md), which builds a dataset in exactly this + format with known truth. +- **A contaminant needs modelling rather than masking** → the model skill's extra-galaxies + branch, fed by the `extra_galaxies_centres.json` you wrote above. +- **The data looks wrong in a way you can't place** → the environment skill + ([`ag_setup_environment`](./ag_setup_environment.md)) if it is the *plot* misbehaving, or + the fit-debugging skill (`ag_debug_fit_failure`) once a fit has actually run. + +Offer (default-yes) to record the session as a dated +`wiki/project/YYYY-MM-DD-.md` entry — the mask radius you chose and *why*, and what +you did about each contaminant, are exactly the decisions you will want justified when you +write the paper. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Data](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_1_introduction/tutorial_2_data.ipynb): + what a CCD image, a noise-map and a PSF actually are, and why the PSF is convolved onto + the model rather than divided out of the data. Its companion, + [chapter 2 tutorial 6](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_6_masking.ipynb), + is the masking lecture. +- **General reference** — [RTD: New user guide](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_2_new_user_guide.html): + routes by system scale and dataset type, and is the page to hand someone deciding how + their data should be organised before they load it. +- **Experienced PyAutoGalaxy user** — [workspace: imaging/data_preparation/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/data_preparation/start_here.py): + the standards checklist with the per-ingredient example scripts and GUIs beneath it. + +## Agent procedural checklist + +1. Establish the pixel scale and the noise-map's definition before writing code. +2. Load the three ingredients, plot the dataset subplot, **print and quote its absolute + path**, and offer to open it once. +3. Run the gate against that figure: settle (a) every contaminant and (b) the mask extent + *with the user*. On a harness that cannot plot, ask the user to do both — never waive it. +4. Check the three standards: eps units, RMS-with-Poisson noise-map, odd normalised PSF. + Convert rather than assume. +5. Apply, in order: noise scaling for contaminants → the mask → adaptive over-sampling + (with every bright centre in `centre_list`). +6. Write the sidecars the model will need (`extra_galaxies_centres.json`, + `light_centre.json`, `info.json`) and plot any centres over the image before saving. +7. Save the script to `scripts/`, hand off to `ag_build_imaging_model`, and offer the + `wiki/project/` entry recording the mask and contaminant decisions. diff --git a/skills/ag_run_search.md b/skills/ag_run_search.md new file mode 100644 index 0000000..47bcefe --- /dev/null +++ b/skills/ag_run_search.md @@ -0,0 +1,415 @@ +--- +name: ag_run_search +description: Run the fit — take a composed model, an analysis and a configured search and drive `search.fit(model=model, analysis=analysis)` to completion, then read what it wrote. Covers the output-folder anatomy and the on-the-fly announcement, the unique-identifier resume behaviour and when it silently reuses a stale fit, `iterations_per_quick_update` and `live_visual_update`, JAX/GPU acceleration and when it actually pays, VRAM checks before a long GPU run, the `PYAUTO_TEST_MODE` smoke levels, and the `if __name__ == "__main__"` parallelisation fix. Use once `ag_build_imaging_model` and `ag_configure_search` have produced `model`, `analysis` and `search` objects, or when a fit is running and the user wants to know what to watch. Not for composing the model or choosing the search (those two skills), not for interpreting the finished posterior (`ag_load_results`), and not for diagnosing a fit that ran but converged somewhere unphysical (`ag_debug_fit_failure`). +--- + +# Running the fit + +This is the moment the inference actually happens. Everything before it was preparation: +a `model` says which morphological parameters are free and what priors they carry, an +`analysis` binds that model to one dataset and knows how to turn a parameter vector into a +log likelihood, and a `search` knows how to explore the parameter space. One call joins +them: + +```python +result = search.fit(model=model, analysis=analysis) +``` + +Statistically, that call is sampling (or optimising) the posterior +$P(\theta | d) \propto \mathcal{L}(d | \theta) P(\theta)$, where the likelihood is the +Gaussian-noise chi-squared of the PSF-convolved model image against the data — walked +line by line in `autogalaxy_workspace:scripts/imaging/likelihood_function.py`. Scientifically, +it is where a Sersic index, an effective radius or a bulge-to-total ratio stops being a guess +and becomes a measurement with errors. What this skill adds is everything *around* the call: +how to watch it, how long to expect, how to make it fast, and how to prove the script works +before committing hours to it. + +Two wiki pages carry the background this skill assumes. Read +[`wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md) for +what each sampler does and how run time scales with the model, and +[`wiki/core/api/analysis_objects.md`](../wiki/core/api/analysis_objects.md) for what the +returned `Result` contains. + +## Ask + +Three questions, and the answers change the branch: + +- *"Do you already have `model`, `analysis` and `search`?"* If not, route to + [`ag_build_imaging_model`](./ag_build_imaging_model.md) and + [`ag_configure_search`](./ag_configure_search.md) first. A fit assembled from a + half-remembered model is the most expensive mistake available here. +- *"Smoke test first, or straight to the production run?"* Almost always smoke test first — + it takes seconds and answers "does this script execute end to end?", which is a different + question from "is this galaxy's Sersic index 3.2?". The smoke branch is below. +- *"How long can this run, and on what?"* A `MultiStartProdigy` fit is minutes; a `Nautilus` + fit of a smooth model is tens of minutes on a GPU and can be hours on a CPU; a pixelised + fit can be far longer. That answer decides foreground vs background, laptop vs cluster, + and whether `live_visual_update` is worth switching on. + +If the data is real observational imaging, the inspection gate in +[`../AGENTS.md`](../AGENTS.md) applies before you get here at all — contaminants and mask +extent are settled in [`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md), not +retrospectively. + +## Branch — the production run + +The deliverable is one script in `scripts/` the user can re-run and edit. Adapted from +`autogalaxy_workspace:scripts/imaging/start_here.py` (the fit and result sections) and +`autogalaxy_workspace:scripts/imaging/modeling.py` (the Nautilus configuration and the VRAM +check): + +```python +""" +Galaxy Structure: Run the Fit +============================= + +Fit a galaxy's surface-brightness profile to CCD imaging: load and mask the data, compose a +linear Sersic bulge, and sample the posterior with Nautilus so the inferred effective radius +and Sersic index come with errors. + +__Contents__ + +- **Imports:** JAX environment first, then the standard trio. +- **Dataset:** Load imaging, mask it, and apply adaptive over-sampling. +- **Model:** Compose the galaxy's light profile. +- **Search:** Configure Nautilus and its update cadence. +- **Analysis:** Bind the model to the data and choose the JAX backend. +- **Model-Fit:** Run the search and announce the output folder. +- **Result:** Read the best-fit model back. +""" +from autogalaxy import jax_wrapper # Sets the JAX environment before other imports + +from pathlib import Path + +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +""" +__Dataset__ + +Galaxy modeling needs three ingredients: the image in CCD counts, a per-pixel RMS noise-map, +and the PSF. The PSF is forward-modelled rather than divided out, which is what separates a +genuinely compact bulge from a seeing-broadened one. `pixel_scales` converts pixels to +arcseconds and must match the instrument (`ag.Imaging.from_fits`, +`PyAutoArray:autoarray/dataset/imaging/dataset.py`). + +The mask sets which pixels enter the likelihood, and its radius is a science choice, not a +default: truncate the outer isophotes and the effective radius and Sersic index are biased +directly. Over-sampling evaluates the light profile several times per pixel where the +gradient is steep, which matters most in the central few tenths of an arcsecond +(`ag.util.over_sample.over_sample_size_via_radial_bins_from`). +""" +DATASET_PATH = Path("dataset") / "imaging" / "" +MASK_RADIUS = 2.5 + +dataset = ag.Imaging.from_fits( + data_path=DATASET_PATH / "data.fits", + psf_path=DATASET_PATH / "psf.fits", + noise_map_path=DATASET_PATH / "noise_map.fits", + pixel_scales=0.1, +) + +mask = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=dataset.pixel_scales, + radius=MASK_RADIUS, +) + +dataset = dataset.apply_mask(mask=mask) + +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=dataset.grid, + sub_size_list=[8, 4, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0)], +) + +dataset = dataset.apply_over_sampling(over_sample_size_lp=over_sample_size) + +""" +__Model__ + +A linear `Sersic` bulge: the `intensity` is solved by linear inversion each iteration rather +than sampled, so the search explores one fewer dimension per component for free +(`autogalaxy_workspace:scripts/imaging/modeling.py` `__Linear Light Profiles__`). +""" +bulge = af.Model(ag.lp_linear.Sersic) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(model.info) + +""" +__Search__ + +`Nautilus` is nested sampling: it returns the full posterior — every parameter's density, +its errors, and the correlations between them. `n_live` trades accuracy against run time and +200 covers the vast majority of galaxy models. `unique_tag` is conventionally the dataset +name, which is what keeps the same model fitted to different data in different folders +(`PyAutoFit:autofit/non_linear/search/nest/nautilus/search.py`). +""" +search = af.Nautilus( + path_prefix=Path("imaging"), + name="sersic", + unique_tag="", + n_live=200, + iterations_per_quick_update=10000, + live_visual_update=False, +) + +""" +__Analysis__ + +`ag.AnalysisImaging` defines the `log_likelihood_function` the search calls, and defaults to +`use_jax=True` when JAX is installed (`PyAutoGalaxy:autogalaxy/imaging/model/analysis.py`). +""" +analysis = ag.AnalysisImaging(dataset=dataset, use_jax=True) + +""" +__Model-Fit__ + +Results are written to the output folder on the fly, from the best model found so far, so the +folder is worth opening the moment the search starts. +""" +print(f"Output folder: {search.paths.output_path.resolve()}") + +result = search.fit(model=model, analysis=analysis) + +""" +__Result__ +""" +print(result.info) +print(result.max_log_likelihood_instance) +``` + +Run it with writable caches if the default locations are not writable — a sandbox, CI, or an +install imported from a Windows mount under WSL +([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)): + +```bash +NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib python scripts/run_fit.py +``` + +`search.fit` itself is `PyAutoFit:autofit/non_linear/search/abstract_search.py`. + +## Branch — what the output folder contains, and when to look + +**Announce the folder at launch, not at the end.** Quote +`search.paths.output_path.resolve()` once the fit is running and say plainly that +`model.results` and `image/fit.png` refresh as the search goes — there is nothing to wait +for. Users new to the stack sit watching a silent terminal because nobody told them. + +The annotated tree of `files/`, `image/`, `model.info`, `model.results`, `search.summary` +and the `` folder is `__Output Folder Layout__` in +[`autogalaxy_workspace/scripts/imaging/modeling.py`](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/modeling.py); +the same section appears in the `interferometer`, `multi_dataset`, `multi_galaxy` and +`cluster` `modeling.py` scripts. Link it rather than copying the tree into the conversation, +where it would rot. The condensed version is also in +[`wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md) +"The output folder". + +Name what to open first, in this order: + +1. `model.results` — the human-readable fit summary, parameters with errors. +2. `image/fit.png` — data, model image, residuals and chi-squared map in one figure. +3. `search.log` — tail it if you want to watch progress from a terminal. + +Depth follows [`_style.md`](./_style.md) "Adaptive depth": for a newcomer, or in teacher +mode, walk all three. For a returning user, quoting the path is enough. + +## Branch — resuming, and the trap inside it + +`` is derived from the model, the search settings and the dataset identifier, so +re-running an identical configuration **resumes** the existing fit rather than restarting it +(`autogalaxy_workspace:scripts/imaging/modeling.py` `__Unique Identifier__`). That is a +feature — an interrupted overnight run picks up where it stopped. + +It is also the sharpest edge in the workflow. **The identifier does not hash the pixel +values of your data.** Swap `data.fits` for a different galaxy, keep the model, the search +and the `unique_tag` identical, and the next run finds a completed fit at the same path and +returns it — silently, in seconds, with the previous galaxy's parameters. If the data +changes, change the `unique_tag`. If you are unsure whether a result is stale, delete the +`` folder and re-run rather than reasoning about it. + +Test mode namespaces its output separately: any active `PYAUTO_TEST_MODE` level inserts a +`test_mode` segment straight after the output root, so a smoke run lands in +`output/test_mode//...` and cannot short-circuit a later real fit +([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)). + +## Branch — the update cadence and the live view + +Two settings control what you see while the fit runs, both on the search +(`PyAutoFit:autofit/non_linear/search/abstract_search.py`): + +- **`iterations_per_quick_update`** — how often the current best model is visualised and + written to disk. The unit depends on the search: likelihood evaluations for `Nautilus` + (hence the large `10000` in `modeling.py`), gradient steps for `MultiStartProdigy` (hence + the much smaller `50` in `start_here.py`). Each update costs roughly ten seconds, so too + low a value and output dominates the run time; too high and you cannot see progress. If + the log is constantly reporting that it is outputting results, raise it. +- **`live_visual_update`** — additionally push the quick-update image to a live surface: a + matplotlib window that refreshes in a plain script, or a self-updating cell in Jupyter and + Colab. The `fit.png` disk write happens either way. Leave it `False` on anything headless + — an HPC node, a background process, CI. + +If the search configuration already records the choice, don't re-ask; that decision belongs +to [`ag_configure_search`](./ag_configure_search.md). + +## Branch — JAX, GPUs, and when the acceleration pays + +Imaging fits run through JAX by default. Installing the extra (`pip install autogalaxy[jax]`) +is all that is required: `ag.AnalysisImaging` then defaults to `use_jax=True`, and the search +driver wraps the likelihood in `jax.vmap(jax.jit(...))` so a whole batch of parameter vectors +evaluates in one call. You will see a one-time log line as the JIT compile starts; every +evaluation after that re-uses the compiled trace. This is described in +`autogalaxy_workspace:scripts/guides/using_jax.py` `__Auto-Enabled Modeling__`, and if JAX is +absent the analysis warns once and falls back to NumPy. + +**When it pays.** The gain scales with how much array work each likelihood evaluation does, +so it is largest for big masks, fine over-sampling, many-component bases and pixelised +reconstructions, and smallest for a tiny mask with one smooth profile. GPU gains dwarf CPU +gains: `autogalaxy_workspace:scripts/imaging/modeling.py` `__JAX__` quotes roughly ten +minutes on a GPU against an hour on a CPU for its example, with a CPU still gaining from +multithreading. Gradient-based searches like `MultiStartProdigy` are not merely faster under +JAX — they are only possible with it, because JAX is what supplies the likelihood's +derivatives and evaluates all the parallel starts in one batched call +(`autogalaxy_workspace:scripts/imaging/start_here.py` `__Multi Start Gradient Optimization__`). + +**Check VRAM before a long GPU run.** A JAX fit must fit inside the GPU's memory, and it +fails at JIT-compile time or on the first likelihood call if it does not: + +```python +analysis.print_vram_use(model=model, batch_size=search.batch_size) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/modeling.py` `__VRAM Use__`. It takes +twenty or thirty seconds, so comment it out once you know your model's footprint. Batch size +is the lever: larger batches cut wall-clock time and raise VRAM, smaller batches the reverse. +A smooth-profile fit on a modest dataset is tens of megabytes; a pixelised reconstruction on +high-resolution data can exceed ten gigabytes. + +**Turn JAX off to debug.** `ag.AnalysisImaging(dataset=dataset, use_jax=False)`, or +`PYAUTO_DISABLE_JAX=1` to force it globally without editing code. NumPy stack traces are far +easier to read, and you can drop a debugger or a `print` into code JAX would otherwise +compile. That is the first lever to pull when deciding whether a failure is a JAX problem at +all — see [`ag_debug_fit_failure`](./ag_debug_fit_failure.md). + +For running on a cluster — SLURM array jobs, one dataset per task, the CPU and GPU batch +scripts, and `number_of_cores` read from `SLURM_CPUS_PER_TASK` — the ground truth is +`autogalaxy_workspace:scripts/guides/hpc/example_cpu_and_gpu.py`. The Python modeling code is +identical for a GPU run; only the batch script changes. This repo has no HPC operations page +yet ([`../PENDING.md`](../PENDING.md)), so cite that script rather than a page that does not +exist. + +## Branch — the smoke test, before you spend hours + +Every script gets one fast structural run before a real one. `PYAUTO_TEST_MODE` has four +levels ([`wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)), and the +choice between them matters more than it looks: + +```bash +# Level 1 — sampler runs with minimal iterations. A real (tiny) search: writes +# model.results, image/fit.png and a samples.csv you can load back. +PYAUTO_TEST_MODE=1 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/run_fit.py + +# Level 2 — sampler bypassed, likelihood called exactly once. The fastest proof that +# the model composes and the likelihood evaluates. +PYAUTO_TEST_MODE=2 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/run_fit.py +``` + +**Level 2 does not write the fit products.** A level-2 run returns a `Result`, proving the +wiring, but leaves no `model.results` and no `image/fit.png` behind — so if the next thing +you want to do is plot or load the fit, use level 1. Level 3 skips the likelihood call as +well, which is only useful for checking that a script's non-fitting scaffolding runs. + +**The parameter values from any test-mode run are meaningless.** The sampler was truncated or +mocked. Never quote a structural parameter measured in test mode; say explicitly that the +number is a wiring check. + +For the fastest possible loop, combine test mode with the dataset and output short-circuits: + +```bash +PYAUTO_TEST_MODE=2 PYAUTO_SKIP_FIT_OUTPUT=1 PYAUTO_SKIP_VISUALIZATION=1 \ + PYAUTO_SMALL_DATASETS=1 PYAUTO_FAST_PLOTS=1 python scripts/run_fit.py +``` + +`PYAUTO_SMALL_DATASETS=1` caps every array, mask and grid to 16 × 16 pixels — so delete any +previously simulated `dataset/` when you toggle it, or a full-resolution dataset on disk will +be reused and mismatch the capped grids. + +## Branch — the error that appears the moment the search starts + +On some combinations of operating system and Python version, a script raises as soon as the +search begins, from Python's process spawning rather than from anything about your model. The +fix is to wrap the whole fit in a function and call it under a main guard: + +```python +def fit(): + from autogalaxy import jax_wrapper # Sets the JAX environment before other imports + + import autofit as af + import autogalaxy as ag + + # ... dataset, model, search, analysis as above ... + result = search.fit(model=model, analysis=analysis) + + +if __name__ == "__main__": + fit() +``` + +Adapted from `autogalaxy_workspace:scripts/guides/modeling/bug_fix.py`, which is the same +`modeling.py` fit with this one structural change. It works for every dataset type, so adopt +it for any modeling script that hits the error. If parallelisation still will not work after +the fix, that script points at the project's support channel. + +## Combine + +- [`ag_plot_fit`](./ag_plot_fit.md) — render the fit yourself, at whatever scaling and + cropping the science needs, rather than reading the auto-generated `image/fit.png`. +- [`ag_load_results`](./ag_load_results.md) — pull the posterior, the best-fit galaxies and + the derived quantities back into Python, in-session or from disk later. +- [`ag_debug_fit_failure`](./ag_debug_fit_failure.md) — the fit finished but the residuals or + the parameters are wrong. +- [`ag_configure_search`](./ag_configure_search.md) — the fit is too slow, or the sampler is + the wrong one for the model. +- [`ag_simulate_dataset`](./ag_simulate_dataset.md) — build a controlled dataset with known + truth and fit that first, which is the cleanest way to separate a model problem from a + data problem. + +When the script and its fit are worth keeping, offer (default-yes) to record the run in a +dated `wiki/project/YYYY-MM-DD-.md` entry: the science question, what was inferred and +with what search, and the output path — per [`_style.md`](./_style.md) property #5. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Practicalities](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_2_practicalities.ipynb): + the practical side of running a fit — output structure, reviewing results as they appear, + and managing run times. +- **General reference** — [RTD: New user guide](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_2_new_user_guide.html): + the decision tree for which fit to run first, by system scale and dataset type. +- **Experienced PyAutoGalaxy user** — [workspace: imaging/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/start_here.py): + the minimal end-to-end imaging fit, with the JAX and quick-update sections this skill + mirrors. + +## Agent procedural checklist + +1. Confirm `model`, `analysis` and `search` exist; if not, route to the model-building and + search-configuration skills. +2. On real observational data, confirm the inspection gate has been satisfied. +3. Write the script to `scripts/` in the generated-script style; never leave it inline only. +4. Validate it with `PYAUTO_TEST_MODE=1` (level 2 if you only need the likelihood proved), + with the cache variables set if the environment needs them. +5. Launch the production run; **quote the absolute output path immediately** and say results + refresh on the fly. +6. Name `model.results` and `image/fit.png` as the first two things to open. +7. If the data changed but the model did not, confirm the `unique_tag` changed too. +8. On a GPU, run the VRAM check before a long fit. +9. On completion, hand off to the results-loading or fit-plotting skill rather than + interpreting the numbers inline. +10. Offer the `wiki/project/` entry. diff --git a/skills/ag_setup_environment.md b/skills/ag_setup_environment.md new file mode 100644 index 0000000..0bba2f2 --- /dev/null +++ b/skills/ag_setup_environment.md @@ -0,0 +1,396 @@ +--- +name: ag_setup_environment +description: Install, diagnose and repair the PyAutoGalaxy environment so galaxy-modelling code can actually run — a fresh pip install into a virtualenv, the JAX and numba extras, writable caches for a sandbox or restricted filesystem, the `activate.sh` interpreter resolution, a shared/cluster checkout, and the Google Colab entry point. Use when an import raises, when `--check-install` returns exit 2 or 3, when the user is starting from nothing, when figures fail on a matplotlib backend, or when a run dies inside a numba or cache path. Produces a saved verification script that proves the install works end to end. Not for API drift against a working install (that is `ag_audit_skill_apis` — a `--check-version` exit 1 means the stack imports fine and the *docs* moved), and not for choosing a search or composing a model. +--- + +# Getting a working PyAutoGalaxy environment + +Nothing else in this workspace matters until `import autogalaxy` succeeds in the +interpreter that will run your script. This is worth a skill of its own because the +failure mode is almost never "the package is missing" — it is that the packages are +installed in a *different* interpreter than the one on `PATH`, or that a compiled-code +cache cannot be written, or that JAX resolved a CPU wheel over a CUDA one. Those three +produce three completely different error messages and one identical user experience. + +The stack is four PyPI packages — `autonerves`, `autoarray`, `autofit`, `autogalaxy` — +and one command installs all of them, because `autogalaxy` declares the other three +transitively. The rationale, the extras table and every version floor live in +[`../wiki/core/operations/installation.md`](../wiki/core/operations/installation.md); +the cache and short-circuit environment variables live in +[`../wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md). Read whichever +of those two the branch you land in points at — this skill is the procedure, they are the +facts. + +## Ask + +One question decides the branch, so ask it before typing anything: + +- *"Are we installing from scratch, repairing an install that used to work, or setting up + a hosted/shared environment (Google Colab, an HPC checkout)?"* + +And one follow-up when the answer is "repair", because it changes everything: +*"what exactly failed — an `ImportError`, a numba or cache path, a matplotlib backend, or +a JAX/GPU message?"* If they don't know, run the diagnosis branch below and read it off +the output rather than guessing. + +## Branch — diagnose before you install anything + +Never install over a broken environment; find out which interpreter you are actually in +first. The repo ships a diagnostic that answers that in one call: + +```bash +source activate.sh +python autoassistant/audit_skill_apis.py --check-install +``` + +It prints the interpreter path, the environment prefix, every resolved version, the file +`autogalaxy` was imported from, and whether the install is a wheel or an editable source +checkout (`autogalaxy_assistant:autoassistant/audit_skill_apis.py`). Three exit codes, +three different repairs: + +| Exit | Means | Go to | +|---|---|---| +| `0` | Ready — the stack imports in *this* interpreter | "Prove it works" below | +| `2` | The packages are absent from this interpreter | "Fresh install", or just activate the right venv | +| `3` | The packages were found but an import **raised** | "Repair a broken import" | + +Exit `2` is far more often the wrong `python` on `PATH` than a genuinely missing install, +which is why the output leads with the interpreter path. Check it against where you +believe you installed before you pip-install anything a second time. + +A separate check answers a different question — whether the *documented* API still matches +the installed one: + +```bash +python autoassistant/audit_skill_apis.py --check-version +``` + +Exit `1` there is **not** an environment problem. The stack imports fine and this repo's +prose has drifted from it (or vice versa); that is +[`ag_audit_skill_apis`](./ag_audit_skill_apis.md)'s job, and it owns the baseline. Don't +reinstall anything in response to it. Exit `2`/`3` from `--check-version` are the install +codes above, forwarded, so a red result is never ambiguous between the two causes. + +## Branch — a fresh install + +Python first: `requires-python = ">=3.12"` across all four packages +(`PyAutoGalaxy:pyproject.toml`), and `autogalaxy_workspace:runtime.txt` targets +`python-3.12.1`, which is the best-tested baseline. On 3.11 and below the install cannot +succeed at all — every earlier wheel on PyPI is yanked, so pip reports "no matching +distribution" rather than silently resolving something ancient. + +```bash +python3.12 -m venv .venv +source .venv/bin/activate +pip install --upgrade pip +pip install "autogalaxy[jax]" numba +``` + +That is deliberately the exact content of `autogalaxy_workspace:requirements.txt` +(`autogalaxy[jax]` plus `numba`), so an environment built this way already satisfies the +example workspace. + +Three things about that command are worth saying out loud, because each is a support +question in disguise: + +- **`--upgrade pip` first is not superstition.** Several of the scientific wheels below + need a modern resolver. +- **The `[jax]` extra is what makes fitting fast**, and it is *not* included by default. + A plain `pip install autogalaxy` gives a fully working NumPy install; the extra adds + just-in-time compilation and the GPU path. It installs **CPU** JAX — for a GPU, follow + the official JAX guide and install the CUDA build **before** `autogalaxy[jax]`, or pip + will happily resolve the CPU wheel over it. +- **`numba` is separate and optional.** It is not in the `jax` extra, which is why the + command names it explicitly. It accelerates the compiled geometry kernels in + PyAutoArray; if it will not build on your platform, the stack runs without it. + +`pynufft` only matters if you fit visibilities — skip it otherwise; the library prints +exactly what to install if you reach interferometer code without it. The full extras +table, the conda route and the editable-clone route (for reading or modifying library +source) are in +[`../wiki/core/operations/installation.md`](../wiki/core/operations/installation.md). + +### How `activate.sh` resolves the interpreter + +This repo ships [`../activate.sh`](../activate.sh), and it resolves the environment +**relative to the script's own location**, so it behaves identically whether you +`source activate.sh` from the repo root or a batch job sources it by absolute path. It +tries two locations in order: + +1. `/.venv/bin/activate` — a virtualenv created inside this clone, which is exactly + what the `python3.12 -m venv .venv` above produces. +2. `$PYAUTO_HPC_BASE/PyAuto/bin/activate` — a shared or cluster checkout. Point + `PYAUTO_HPC_BASE` at a directory holding a `PyAuto/` virtualenv alongside source + checkouts of the four libraries; the script activates the venv and prepends those + checkouts to `PYTHONPATH` in dependency order. + +If neither exists it prints a one-line diagnostic rather than failing silently. If your +venv lives somewhere else entirely (`~/venv/PyAuto`, say), edit `VENV` at the top of the +file rather than working around it every session. + +## Branch — repair a broken import (`--check-install` exit 3) + +Exit 3 means the packages resolved and then something raised. In practice it is one of +four things, in descending order of frequency. + +**1. An unwritable cache.** The libraries cache compiled code and plot configuration under +your home directory, and in a container, in CI, on a read-only home, or on an install +imported from a `/mnt/c/...` Windows mount under WSL, those writes fail — usually with a +confusing message buried inside a numba traceback: + +```bash +export NUMBA_CACHE_DIR=/tmp/numba_cache +export MPLCONFIGDIR=/tmp/matplotlib +``` + +Set them once per shell or bake them into the venv activation. `--check-install` already +sets equivalent temporary defaults for its *own* imports, but only when they are unset, so +a deliberate choice of yours is never overridden. + +**2. The JAX compilation cache.** With JAX installed, +`PyAutoNerves:autonerves/jax_wrapper.py` enables a persistent compilation cache under +`$XDG_CACHE_HOME/pyauto_jax` (or `~/.cache/pyauto_jax`) so the minutes spent compiling a +given model/data shape are paid once per machine rather than once per process. In a +sandbox, point `JAX_COMPILATION_CACHE_DIR` somewhere writable, or set it to the **empty +string** to disable the cache outright. The same module sets `XLA_FLAGS` and +`JAX_ENABLE_X64=True` *before* JAX is imported — which is why workspace scripts that use +JAX open with `from autogalaxy import jax_wrapper` ahead of their other imports. Setting +those after JAX has loaded has no effect, so an import you reorder "for tidiness" can +silently change your numerics. + +**3. The matplotlib backend.** If a figure call hangs, raises, or kills the process with +no message, the backend is misconfigured for your system. It is read from +`config/visualize/general.yaml` → `general:` → `backend:`, whose default is `default` +(your system's own). `autogalaxy_workspace:welcome.py` exists largely to surface this one +interactively, and suggests `TKAgg`, `Qt5Agg` or `Qt4Agg` as replacements. Since every +plot this workspace produces is written to disk rather than displayed (see +[`_style.md`](./_style.md) "Plot output and path announcement"), an `Agg` backend is +usually the right answer for an agent-driven session. + +**4. A dependency pushed past its cap.** The stack pins upper bounds deliberately — +`scipy<=1.17.1`, `astropy>=5.0,<=7.2.0`, `numpy>=1.24.0,<3.0.0`, `jax>=0.7.0,<0.11.0`, +and exact pins on `dynesty` and `nautilus-sampler`. These are the ones that bite, because +forcing a newer version tends to break inside a mesh routine or a sampler's internals +rather than at import. The full list is in the installation page's "Version floors and +caps". + +One more failure that looks like a broken install and is not: a **working-directory** +error. PyAutoGalaxy resolves `config/`, `dataset/` and `output/` relative to the current +directory, so running a workspace script from the wrong place produces import-shaped and +file-not-found-shaped errors that have nothing to do with the install. +`PyAutoGalaxy:docs/installation/troubleshooting.md` lists it first for good reason. Run +workspace scripts from the workspace root; run this repo's scripts from this repo's root. +A related, harmless warning: `PyAutoNerves:autonerves/workspace.py` compares a workspace +clone's version to the installed library and warns on a mismatch — +`PYAUTO_SKIP_WORKSPACE_VERSION_CHECK=1` silences it when the two intentionally diverge. + +## Branch — Google Colab + +Don't hand-write install cells. `PyAutoNerves:autonerves/setup_colab.py` exposes a +per-project entry point, `for_autogalaxy(raise_error_if_not_gpu=...)`, which installs the +stack, clones the workspace and checks whether JAX actually found a GPU. The guard that +makes it work both in and out of Colab is the one every workspace `start_here` uses +(`autogalaxy_workspace:start_here.py`): + +```python +try: + import google.colab +except ImportError: + from autogalaxy import setup_colab as _setup_colab +else: + import importlib + import subprocess + import sys + + subprocess.check_call( + [sys.executable, "-m", "pip", "install", "autonerves", "--no-deps"] + ) + _setup_colab = importlib.import_module("autonerves.setup_colab") + +_setup_colab.for_autogalaxy( + raise_error_if_not_gpu=False # True to hard-fail when no GPU was found +) +``` + +Outside Colab the same block is a no-op verification pass, so it is safe to leave at the +top of a script you also run locally. Set `raise_error_if_not_gpu=True` when a GPU is the +whole point of the run and you would rather fail loudly than discover 40 minutes later +that you were on a CPU. + +## Prove it works + +An install that imports is not yet an install that fits galaxies. Save a verification +script — it takes seconds to run and it exercises the config layer, the array structures, +a light-profile evaluation and the plot path, which between them cover every failure above. + +```python +""" +Environment Verification: PyAutoGalaxy +====================================== + +Prove a PyAutoGalaxy install is usable end to end, not merely importable. Each section +exercises one layer that fails independently: package resolution, the YAML configuration +tree, the array/grid structures, a light-profile evaluation, and writing a figure to disk. +Run this after any install, repair or interpreter change, and read the printed paths — a +wrong interpreter or an unwritable cache shows up here rather than three hours into a fit. + +__Contents__ + +- **Imports:** Import the stack and report which interpreter and files were resolved. +- **Grid:** Build a `Grid2D` and confirm the pixel-to-arcsecond conversion. +- **Light Profile:** Evaluate a Sersic profile on the grid — the first real numerical work. +- **Plot:** Write a figure to disk and print its absolute path. +""" + +""" +__Imports__ + +Importing all four packages separately is deliberate: they form a dependency chain +(`autonerves` → `autoarray` → `autofit` → `autogalaxy`), so a failure names the layer that +broke rather than reporting a generic `autogalaxy` error. `__file__` is printed because the +single most common install problem is not a missing package but a *different interpreter* +than the one you installed into. +""" +import sys +from pathlib import Path + +import autonerves +import autoarray +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +print(f"interpreter: {sys.executable}") +print(f"autonerves={autonerves.__version__} autoarray={autoarray.__version__}") +print(f"autofit={af.__version__} autogalaxy={ag.__version__}") +print(f"autogalaxy imported from: {ag.__file__}") + +""" +__Grid__ + +The `Grid2D` is the (y,x) coordinate grid every light profile is evaluated on, and +`pixel_scales` is the arcseconds-per-pixel conversion that ties it to a real detector. +Constructing one touches the configuration tree, so a config-resolution failure surfaces +here (`PyAutoArray:autoarray/structures/grids/uniform_2d.py`). +""" +grid = ag.Grid2D.uniform(shape_native=(50, 50), pixel_scales=0.1) + +print(f"grid shape: {grid.shape_native}, pixel scale: {grid.pixel_scales}") + +""" +__Light Profile__ + +An elliptical Sersic is the workhorse profile of galaxy morphology, and evaluating its +image is the first genuinely numerical operation — it is where a numba cache that cannot be +written, or a NumPy pushed past its version cap, actually fails +(`PyAutoGalaxy:autogalaxy/profiles/light/standard/sersic.py`). `ell_comps` is the +elliptical-components parameterisation used throughout the library in place of an +axis-ratio and position-angle pair. +""" +bulge = ag.lp.Sersic( + centre=(0.0, 0.0), + ell_comps=(0.2, 0.1), + intensity=1.0, + effective_radius=0.8, + sersic_index=4.0, +) + +image = bulge.image_2d_from(grid=grid) + +print(f"image sum: {float(image.sum()):.4f}, peak: {float(image.max()):.4f}") + +""" +__Plot__ + +The final layer is visualisation, which depends on the matplotlib backend and on a writable +`MPLCONFIGDIR`. The plot API is function-based: `output_path`, `output_filename` and +`output_format` are passed straight to the call, and nothing is displayed interactively — +so this works identically in a terminal, a notebook and a headless job. +""" +PLOT_DIR = Path("scripts") / "scratch" / "environment_check" + +aplt.plot_array( + array=image, + title="Sersic Light Profile", + output_path=PLOT_DIR, + output_filename="sersic", + output_format="png", +) + +print(f"Saved to: {PLOT_DIR.resolve()}") +``` + +Run it, then **quote the printed absolute path back to the user and offer to open it** +(`xdg-open` on Linux, `open` on macOS, `explorer.exe` or `wslview` from WSL). A figure the +user can actually see is the only proof that the last layer works. + +For a guided, interactive version of the same idea, `autogalaxy_workspace:welcome.py` +walks a new user through the working-directory rule and the matplotlib backend with +`input()` prompts — worth pointing at rather than reproducing when someone is setting up +their own machine for the first time. + +## Fast iteration once it works + +While you are still shaping a script rather than measuring a galaxy, the stack ships flags +that turn a fit into a structural check. `PYAUTO_TEST_MODE=1` cuts the sampler to its +minimum iterations, `=2` bypasses the sampler and calls the likelihood **once**, and `=3` +skips the likelihood too (`PyAutoNerves:autonerves/test_mode.py`). A representative +combination, with the cache variables from above: + +```bash +PYAUTO_TEST_MODE=2 PYAUTO_SMALL_DATASETS=1 \ + NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/your_script.py +``` + +Two consequences catch everyone once. Test-mode output is **namespaced** into +`output/test_mode/...`, deliberately, so a cached structural run cannot short-circuit a +later real fit with "Fit Already Completed". And the returned parameter values are **not +physically meaningful at any level** — never quote an effective radius or Sersic index +measured in test mode. `PYAUTO_DISABLE_JAX=1` is the other lever worth knowing: it forces +every analysis onto NumPy, whose tracebacks are far easier to read than JAX's when you are +isolating whether a failure is a JAX problem at all. The full flag list is in +[`../wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md). + +## Combine — where to go next + +- **The install is clean and you have your own FITS files** → the data-preparation skill + ([`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md)), which also owns the + real-data inspection gate. +- **The install is clean and you have no data yet** → simulate one with known truth + ([`ag_simulate_dataset`](./ag_simulate_dataset.md)); it is the fastest way to confirm the + whole modelling loop runs before real data is at stake. +- **`--check-version` went red** → [`ag_audit_skill_apis`](./ag_audit_skill_apis.md), which + owns the baseline and the four other currency checks. +- **A fit ran and then failed or returned nonsense** → the fit-failure debugging skill + (`ag_debug_fit_failure`), not this one; the environment is no longer the suspect once a + likelihood has been evaluated. + +Ask if you want the editable-clone route instead — it is the setup you need if you intend +to read or modify library source rather than only call it. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Visualization and setup](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_1_introduction/tutorial_0_visualization.ipynb): + the setup tutorial — working directory, matplotlib options, and how figures are + configured. The right first stop when plots are the thing misbehaving. +- **General reference** — [RTD: Installation overview](https://pyautogalaxy.readthedocs.io/en/latest/installation/overview.html): + the upstream installation index, with the pip, conda, source, numba and troubleshooting + pages beneath it. +- **Experienced PyAutoGalaxy user** — [workspace: start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/start_here.py): + the canonical imaging script, including the Colab setup block and the JAX notes — the + fastest end-to-end proof that a fresh environment really works. + +## Agent procedural checklist + +1. Ask: fresh install, repair, or hosted/shared environment? +2. `source activate.sh`; `python autoassistant/audit_skill_apis.py --check-install`. +3. Exit 2 → check the printed interpreter *before* reinstalling; then the fresh-install + branch. Exit 3 → the repair branch, starting with the cache variables. +4. Never install over a broken environment, and never re-pin a baseline to silence a check. +5. Save and run the verification script; print and then **quote** the figure's absolute + path, and offer to open it once. +6. `--check-version` exit 1 → hand off to `ag_audit_skill_apis`; do not reinstall. +7. Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry only if the + environment work was non-trivial — a genuine repair worth remembering, not a routine + install. diff --git a/skills/ag_simulate_dataset.md b/skills/ag_simulate_dataset.md new file mode 100644 index 0000000..bb46863 --- /dev/null +++ b/skills/ag_simulate_dataset.md @@ -0,0 +1,504 @@ +--- +name: ag_simulate_dataset +description: Simulate CCD imaging (or interferometer visibilities) of a galaxy with known truth — a grid and over-sampling, a PSF, an exposure time and background sky, one or more galaxies built from light profiles, then FITS output plus a `galaxies.json` truth record. Covers targeting a signal-to-noise ratio directly with the `lp_snr` profiles, simulating a whole sample in a loop for population or machine-learning work, adding a contaminating extra galaxy with its noise-scaling mask, and the `should_simulate` auto-simulation convention every workspace script uses. Use to rehearse the modelling loop before real data, to build training sets, to test a recovery, or when a user has no data yet. Not for loading or preparing real observations (`ag_prepare_imaging_data`), and not for fitting. +--- + +# Simulating a galaxy dataset with known truth + +A simulated dataset is the only situation in galaxy modelling where you know the answer. +That makes it worth far more than a convenience: it is how you establish that your model +can recover a parameter at all, how you find out what signal-to-noise you need before you +propose for telescope time, how you separate "the fit is wrong" from "the data cannot +constrain this", and how you generate the thousands of labelled images a neural network +needs. It is also the fastest way to rehearse the whole loop — compose, search, inspect — +before real data is at stake. + +Statistically, what you are building is a draw from the likelihood's own generative model: +take a noise-free galaxy image, convolve it with the PSF, scale to counts through an +exposure time, add a background sky, and draw Poisson noise. Because that is exactly the +forward model the likelihood inverts, a fit to simulated data tests the *inference* in +isolation — any bias you find is yours, not the data's. Simulated data is also **exempt from +the real-data inspection gate** (see +[`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md)) for the same reason: you already +know where every component is. + +The canonical scripts are `autogalaxy_workspace:scripts/imaging/simulator.py` (bulge, disk +and a contaminating neighbour) and `autogalaxy_workspace:scripts/imaging/simulator_sersic.py` +(a single Sersic). What each light profile means physically is +[`../wiki/core/concepts/light_profiles.md`](../wiki/core/concepts/light_profiles.md); the +dataset objects are [`../wiki/core/api/datasets.md`](../wiki/core/api/datasets.md). + +## Ask + +- *"What are you simulating for — rehearsing the workflow, a recovery test, a training set, + or a signal-to-noise feasibility study?"* This picks the branch, and the last two change + the API you want (`lp_snr`, and the sample loop). +- *"One galaxy or several, and should a contaminating neighbour be included?"* Including one + is worth doing deliberately: it is how you rehearse the noise-scaling step you will + certainly need on real data. +- *"What instrument are you standing in for?"* — it sets `pixel_scales`, the PSF width and a + plausible exposure time, and makes the simulation answer a question about *your* telescope + rather than an abstract one. + +## Branch — a single Sersic + +The smallest useful simulation, and the one to start from. Adapted from +`autogalaxy_workspace:scripts/imaging/simulator_sersic.py`. + +```python +""" +Simulator: Single Sersic +======================== + +Simulate CCD imaging of a galaxy whose light is a single elliptical Sersic profile, with +known truth parameters recorded alongside the data. The simulation applies the same forward +model a fit inverts — evaluate the profile on an over-sampled grid, convolve with the PSF, +scale through an exposure time, add a background sky and draw Poisson noise — so a fit to +the output tests the inference rather than the data. + +__Contents__ + +- **Imports:** Import the required libraries. +- **Grid:** Build the coordinate grid and its adaptive over-sampling. +- **PSF and Simulator:** Define the optical blurring and the noise properties. +- **Galaxies:** Define the true galaxy whose parameters we intend to recover. +- **Output:** Write the dataset to FITS, plus a `galaxies.json` record of the truth. +""" + +""" +__Imports__ +""" +from pathlib import Path + +import autogalaxy as ag +import autogalaxy.plot as aplt + +DATASET_PATH = Path("dataset") / "imaging" / "simple_sersic" + +""" +__Grid__ + +The `Grid2D` is the set of (y,x) arcsecond coordinates the galaxy's light is evaluated on, +and `pixel_scales` ties it to a detector — set it to the instrument you are standing in for. +Over-sampling then evaluates the profile on a finer sub-grid where its intensity gradient is +steep: a Sersic varies enormously across the central pixel, so a single evaluation at the +pixel centre under-counts the flux and produces a subtly wrong image. Simulations use higher +sub-sampling than fits (32x32 in the centre rather than 8x8) because the cost is paid once +rather than at every likelihood call +(`PyAutoArray:autoarray/operators/over_sampling/over_sample_util.py`). +""" +grid = ag.Grid2D.uniform(shape_native=(100, 100), pixel_scales=0.1) + +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=grid, + sub_size_list=[32, 8, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0)], +) + +grid = grid.apply_over_sampling(over_sample_size=over_sample_size) + +""" +__PSF and Simulator__ + +Every real CCD image is blurred by the telescope optics, so a simulation without a PSF is not +imaging data. A Gaussian kernel is a reasonable stand-in when you do not have a real PSF +(FWHM = 2.35 * sigma). The `SimulatorImaging` object then supplies the rest of what makes +data noisy rather than ideal: `exposure_time` converts electrons per second to counts and so +sets the Poisson noise level, and `background_sky_level` adds the sky's contribution to the +noise budget (`PyAutoArray:autoarray/dataset/imaging/simulator.py`). Raising the exposure +time raises the signal-to-noise; the ratio of the two controls how much of the noise is +counting statistics versus sky. +""" +psf = ag.Convolver.from_gaussian( + shape_native=(11, 11), sigma=0.1, pixel_scales=grid.pixel_scales +) + +simulator = ag.SimulatorImaging( + exposure_time=300.0, + psf=psf, + background_sky_level=0.1, + add_poisson_noise_to_data=True, +) + +""" +__Galaxies__ + +The truth. `sersic_index` sets the profile's concentration — n ~ 4 is a de Vaucouleurs-like +spheroid, n ~ 1 an exponential disk — and `effective_radius` is the radius containing half +the profile's light. For a *simulation* the axis-ratio and position-angle parameterisation is +the intuitive one, so `ag.convert.ell_comps_from` translates it into the `ell_comps` the +profile actually takes; a fit uses `ell_comps` directly, because a position angle's periodic +boundary is pathological for a non-linear search +(`PyAutoGalaxy:autogalaxy/profiles/light/standard/sersic.py`). +""" +galaxy = ag.Galaxy( + redshift=0.5, + bulge=ag.lp.Sersic( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.9, angle=45.0), + intensity=1.0, + effective_radius=0.8, + sersic_index=4.0, + ), +) + +galaxies = ag.Galaxies(galaxies=[galaxy]) + +dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=grid) + +print(f"peak signal-to-noise: {float(dataset.signal_to_noise_map.max()):.1f}") + +""" +__Output__ + +The FITS files are written in exactly the layout `ag.Imaging.from_fits` expects, so the +simulated dataset is indistinguishable from real data to everything downstream. The +`galaxies.json` record is the part people skip and regret: it stores the true light profiles +so that months later you can still state what the fit was supposed to recover. Load it back +with `ag.from_json`. +""" +aplt.fits_imaging( + dataset=dataset, + data_path=DATASET_PATH / "data.fits", + psf_path=DATASET_PATH / "psf.fits", + noise_map_path=DATASET_PATH / "noise_map.fits", + overwrite=True, +) + +ag.output_to_json(obj=galaxies, file_path=DATASET_PATH / "galaxies.json") + +aplt.subplot_imaging_dataset( + dataset=dataset, output_path=DATASET_PATH, output_format="png" +) +aplt.subplot_galaxies( + galaxies=galaxies, grid=grid, output_path=DATASET_PATH, output_format="png" +) + +print(f"Saved to: {DATASET_PATH.resolve()}") +``` + +Quote that path back and offer to open it once. The `subplot_galaxies` figure is the useful +one for a simulation: it shows the noise-free truth beside the profile's components, which +is what you compare a fit against. + +## Branch — bulge, disk, and a contaminating neighbour + +The realistic default. Two components let you pose the question galaxy morphology usually +cares about — the bulge-to-total light ratio — and a faint offset neighbour lets you +rehearse the noise-scaling step. Adapted from +`autogalaxy_workspace:scripts/imaging/simulator.py`. + +```python +extra_galaxy_centre = (2.2, 1.6) + +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=grid, + sub_size_list=[32, 8, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0), extra_galaxy_centre], +) +grid = grid.apply_over_sampling(over_sample_size=over_sample_size) + +galaxy = ag.Galaxy( + redshift=0.5, + bulge=ag.lp.Sersic( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.9, angle=45.0), + intensity=1.0, + effective_radius=0.6, + sersic_index=3.0, + ), + disk=ag.lp.Exponential( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.7, angle=30.0), + intensity=0.5, + effective_radius=1.6, + ), +) + +extra_galaxy = ag.Galaxy( + redshift=0.5, + light=ag.lp.ExponentialSph( + centre=extra_galaxy_centre, intensity=1.0, effective_radius=0.3 + ), +) + +galaxies = ag.Galaxies(galaxies=[galaxy, extra_galaxy]) +dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=grid) +``` + +Three details carry weight. The bulge is compact with a high `sersic_index` and the disk is +extended with n fixed at 1 by construction — that contrast is what makes the decomposition +identifiable at all. The neighbour's centre appears in `centre_list` as well as in its own +profile, so it is over-sampled properly rather than being the one badly-evaluated object in +the frame. And the attribute names (`bulge`, `disk`, `light`) are arbitrary labels you choose +— name them for what they measure. + +Write the neighbour's noise-scaling mask in the same script, so the modelling examples can +load it without a separate preparation step: + +```python +mask_extra_galaxies = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=dataset.pixel_scales, + centre=extra_galaxy_centre, + radius=3.0 * 0.3, # ~3x the neighbour's effective radius + invert=True, # `True` inside the circle: the region whose noise is scaled. +) + +aplt.fits_array( + array=mask_extra_galaxies, + file_path=DATASET_PATH / "mask_extra_galaxies.fits", + overwrite=True, +) +``` + +Deriving the radius from the neighbour's own `effective_radius` keeps the two in sync if you +later change it. The strategies this mask feeds into are +[`../wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md). + +## Branch — target a signal-to-noise ratio directly + +Often the question is not "what intensity?" but "what if this galaxy were detected at S/N = +20?". Choosing an `intensity` that lands there is guesswork; the `lp_snr` profiles remove it +by solving for the intensity that achieves a requested ratio, using the simulator's +`exposure_time` and `background_sky_level` to do so. Adapted from +`autogalaxy_workspace:scripts/imaging/features/simulator_manual_signal_to_noise.py`. + +```python +galaxy_0 = ag.Galaxy( + redshift=0.5, + bulge=ag.lp_snr.Sersic( + signal_to_noise_ratio=20.0, + centre=(0.0, -1.0), + ell_comps=(0.25, 0.1), + effective_radius=0.8, + sersic_index=2.5, + ), +) + +galaxy_1 = ag.Galaxy( + redshift=0.5, + bulge=ag.lp_snr.Sersic( + signal_to_noise_ratio=10.0, + centre=(0.0, 1.0), + ell_comps=(0.0, 0.1), + effective_radius=0.6, + sersic_index=3.0, + ), +) + +dataset = simulator.via_galaxies_from( + galaxies=ag.Galaxies(galaxies=[galaxy_0, galaxy_1]), grid=grid +) +``` + +Note the profiles take **no** `intensity` — that is the parameter being solved for. The +trade-off is that `exposure_time` and `background_sky_level` no longer set the S/N, only the +*balance* between counting noise and sky noise. So they still matter for realism: doubling +the exposure time shifts the noise budget toward Poisson statistics even though the S/N is +pinned. Choose them to match your instrument if the noise character matters to the study. + +This is also the profile family to reach for in a feasibility question — simulate the same +galaxy at S/N 10, 20 and 50, fit each, and see where the parameter you care about stops +being constrained. That is a far more honest answer than a single simulation at whatever +intensity happened to be typed. + +## Branch — a sample, for population or machine-learning work + +To generate many galaxies, draw each one's parameters from a distribution and loop. Adapted +from `autogalaxy_workspace:scripts/imaging/simulator_sample.py`, which draws directly from a +NumPy generator — deliberately, because these are *truths* for synthetic data, not a model +being fitted, so there is no reason to involve the model-composition API. + +```python +import numpy as np + +rng = np.random.default_rng() + + +def _clipped_ell_comp() -> float: + return float(np.clip(rng.normal(0.0, 0.2), -1.0, 1.0)) + + +def _random_galaxy() -> ag.Galaxy: + bulge = ag.lp_snr.Sersic( + centre=(0.0, 0.0), + ell_comps=(_clipped_ell_comp(), _clipped_ell_comp()), + effective_radius=float(rng.uniform(1.0, 5.0)), + sersic_index=float(np.clip(rng.normal(4.0, 1.0), 0.8, 5.0)), + signal_to_noise_ratio=float(rng.uniform(20.0, 60.0)), + ) + return ag.Galaxy(redshift=0.5, bulge=bulge) + + +total_datasets = 3 + +for sample_index in range(total_datasets): + + sample_path = DATASET_PATH / f"dataset_{sample_index}" + + galaxies = ag.Galaxies(galaxies=[_random_galaxy()]) + dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=grid) + + aplt.fits_imaging( + dataset=dataset, + data_path=sample_path / "data.fits", + psf_path=sample_path / "psf.fits", + noise_map_path=sample_path / "noise_map.fits", + overwrite=True, + ) + ag.output_to_json(obj=galaxies, file_path=sample_path / "galaxies.json") +``` + +The `lp_snr` bulge is what makes the loop usable: with a plain `intensity` drawn from a +distribution, a fraction of your sample would be invisible and another fraction saturated. +Pinning S/N per galaxy guarantees every image is detectable. Write a `galaxies.json` per +dataset — with a sample, the truth record stops being a nicety and becomes the labels. + +For a one-line random galaxy without writing your own distributions, +`ag.model_util.random_galaxy_for_simulation_from()` returns one drawn from the library's own +choices, and `ag.model_util.SIMULATOR_RANDOM_GALAXY_SUMMARY` prints what those choices are — +read that summary before using it in anything you will publish, so the prior you inherited +is one you actually endorse. + +## Branch — interferometer visibilities + +For sub-mm and radio data (ALMA, JVLA, LOFAR) the observable is visibilities in the uv-plane, +not an image, and fitting them there avoids the correlated noise a dirty image carries. +Adapted from `autogalaxy_workspace:scripts/interferometer/simulator.py`. + +```python +uv_wavelengths = ag.ndarray_via_fits_from( + file_path=Path("dataset") / "interferometer" / "uv_wavelengths" / "sma.fits", hdu=0 +) + +simulator = ag.SimulatorInterferometer( + uv_wavelengths=uv_wavelengths, + exposure_time=300.0, + noise_sigma=1000.0, + transformer_class=ag.TransformerNUFFT, +) + +real_space_grid = ag.Grid2D.uniform(shape_native=(800, 800), pixel_scales=0.05) + +dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=real_space_grid) + +aplt.subplot_interferometer_dirty_images( + dataset=dataset, output_path=DATASET_PATH, output_format="png" +) + +aplt.fits_interferometer( + dataset=dataset, + data_path=DATASET_PATH / "data.fits", + noise_map_path=DATASET_PATH / "noise_map.fits", + uv_wavelengths_path=DATASET_PATH / "uv_wavelengths.fits", + overwrite=True, +) +``` + +Four differences from imaging are worth knowing. The baselines come from a real array's +`uv_wavelengths` — the workspace ships SMA (low resolution, very fast) and ALMA files; +swapping the file swaps the instrument. There is **no PSF and no over-sampling**: the image +is evaluated in real space and Fourier-transformed, and interferometers do not observe in a +way that makes over-sampling necessary. `noise_sigma` replaces the exposure-time-driven +Poisson model. And the thing you *look at* is the dirty image, which is a diagnostic — the +fit itself happens in the uv-plane. The physics is +[`../wiki/core/concepts/interferometer_theory.md`](../wiki/core/concepts/interferometer_theory.md). +Full interferometer modelling is a separate skill (`ag_build_interferometer_model`, not yet +written — see [`../PENDING.md`](../PENDING.md)). + +## The auto-simulation convention + +Every workspace script that fits a dataset first checks whether it exists and runs the +matching simulator if it does not, so an example can be run from a fresh clone with no +manual step: + +```python +if ag.util.dataset.should_simulate(str(dataset_path)): + import subprocess + import sys + + subprocess.run([sys.executable, "scripts/imaging/simulator.py"], check=True) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/modeling.py`. Adopt the same pattern in +scripts you write, and know the sharp edge: `should_simulate` tests **existence**, not +whether the dataset on disk matches the simulator that would produce it. Change the +simulator's parameters and re-run the fit and you will silently fit the *old* data. The same +applies to `PYAUTO_SMALL_DATASETS=1`, which caps grids and masks to a small size: a +full-resolution dataset already on disk will be reused and mismatch the capped grids. Delete +the dataset directory when you change either +([`../wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md)). + +## JAX, for many simulations + +For parameter sweeps, mock-data studies or batch figure generation, build the simulator with +`use_jax=True` — the simulation then runs on JAX arrays (and the GPU, when one is +configured): + +```python +simulator_jax = ag.SimulatorImaging( + exposure_time=300.0, + psf=psf, + background_sky_level=0.1, + add_poisson_noise_to_data=True, + use_jax=True, +) + +dataset_jax = simulator_jax.via_galaxies_from(galaxies=galaxies, grid=grid) +``` + +`dataset_jax.data.array` is a `jax.Array`, and the plot and FITS helpers call +`numpy.asarray()` internally, so saving and plotting need no manual conversion. + +Wrapping the whole simulation call in a jitted function — worth it across *many* calls — +additionally requires registering the galaxy classes as JAX pytrees before the first jitted +call, and the helper that does that registration is **not in the released stack yet** (it +lives on the libraries' development branches only). Until it ships in a release, use the +eager `use_jax=True` call above; the `__JAX Variant__` section of +`autogalaxy_workspace:scripts/imaging/simulator.py` carries the jitted recipe for a +source install. **For interferometer simulations the jitted path does not work either way** — +use the eager call; `autogalaxy_workspace:scripts/interferometer/simulator.py` documents why. + +## Combine — where this hands off + +- **Fit the dataset you just made** → [`ag_build_imaging_model`](./ag_build_imaging_model.md) + to compose the model, then [`ag_configure_search`](./ag_configure_search.md) and the + run-search skill (`ag_run_search`). Compose a model that *matches* the truth first: if it + cannot recover parameters it was handed exactly, the problem is the inference, and you have + learned that cheaply. +- **Compare a fit against the truth** → load `galaxies.json` with `ag.from_json` and put the + true and inferred values side by side; the results-loading and fit-plotting skills + (`ag_load_results`, `ag_plot_fit`) do the rest. +- **You actually have real data** → [`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md), + which also owns the inspection gate this branch is exempt from. + +Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry when the simulation is +part of a study — the truth parameters and the S/N you chose *are* the experimental design, +and a recovery test is meaningless later without them. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Grids and galaxies](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_1_introduction/tutorial_1_grids_and_galaxies.ipynb): + how a light profile becomes an image on a grid, and what each Sersic parameter changes — + the foundation under every simulation. +- **General reference** — [RTD: Start here overview](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_1_start_here.html): + the core API end to end, including the simulation section, in one page. +- **Experienced PyAutoGalaxy user** — [workspace: imaging/simulator.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/simulator.py): + the full bulge-plus-disk simulator with the contaminant mask and the executed JAX variant. + +## Agent procedural checklist + +1. Ask what the simulation is *for*; it decides the branch and whether `lp_snr` is wanted. +2. Set `pixel_scales`, PSF width and exposure time from a real instrument, not defaults. +3. Over-sample with every bright centre in `centre_list` — including any neighbour. +4. Write FITS **and** `galaxies.json`; write the neighbour's `mask_extra_galaxies.fits` in + the same script if one exists. +5. Plot the dataset and the noise-free galaxies subplot; print and **quote** the absolute + path, offering to open it once. +6. Delete the dataset directory before re-running with changed parameters or with + `PYAUTO_SMALL_DATASETS=1` — `should_simulate` only checks existence. +7. Hand off to `ag_build_imaging_model`, and offer the `wiki/project/` entry when the + simulation is part of a study. diff --git a/wiki/core/api/aggregator.md b/wiki/core/api/aggregator.md index 4d45010..b0e00ae 100644 --- a/wiki/core/api/aggregator.md +++ b/wiki/core/api/aggregator.md @@ -22,7 +22,7 @@ sources: - scripts/ellipse/database.py pinned_commit: d6db2643b9f2cd418efc9473f560dc2a2d459c73 last_updated: 2026-08-01 -content_sha256: 765f038d0ede9d9517685d0bf9d9f0da2fb2be5e78a506fcce0427c8f0593902 +content_sha256: 3d3a0d22f674bca887b20671b36eb9d2b4b6cdd2eb5477b65ca56cbcc02bf6f4 --- # Aggregator @@ -37,7 +37,7 @@ Two layers stack: 1. **PyAutoFit's aggregator** loads whatever a search wrote to disk: samples, the model, the dataset, stored FITS and PNG files. -2. **PyAutoGalaxy's `ag.agg`** rebuilds *galaxy-domain* objects from those files — `ag.Galaxies`, +2. **PyAutoGalaxy's `ag.agg`** rebuilds *galaxy-domain* objects from those files — galaxy lists, `ag.FitImaging`, `ag.Ellipse` — as generators, so you can compute a galaxy quantity per fit without holding every fit in memory. @@ -141,7 +141,7 @@ Each `*Agg` class takes an aggregator and hands back a generator of PyAutoGalaxy |---|---| | `ag.agg.ImagingAgg` | `ag.Imaging` datasets as fitted | | `ag.agg.InterferometerAgg` | `ag.Interferometer` datasets | -| `ag.agg.GalaxiesAgg` | `ag.Galaxies` | +| `ag.agg.GalaxiesAgg` | a plain `list` of `ag.Galaxy` — **not** an `ag.Galaxies` | | `ag.agg.FitImagingAgg` | `ag.FitImaging` | | `ag.agg.FitInterferometerAgg` | `ag.FitInterferometer` | | `ag.agg.EllipsesAgg` | `ag.Ellipse` lists | @@ -161,6 +161,13 @@ for galaxies, dataset in zip(galaxies_gen, dataset_gen): Adapted from `autogalaxy_workspace:scripts/guides/results/aggregator/models.py`. +**`GalaxiesAgg` yields a list, not an `ag.Galaxies`** (`PyAutoGalaxy:autogalaxy/aggregator/galaxies.py` +returns `List[Galaxy]`), exactly as `ag.from_json` on a `files/galaxies.json` does. The list is +enough for indexing (`galaxies[0].bulge`), for `ag.FitImaging(dataset=…, galaxies=galaxies)`, and +for `aplt.subplot_galaxies` (which wraps it internally) — but it has no collection methods, so +`galaxies.image_2d_from(grid=grid)` raises `AttributeError`. Re-wrap when you need them: +`galaxies = ag.Galaxies(galaxies=galaxies)`. + The dataset classes expose `dataset_gen_from()`. The model classes share four generators, which is where the aggregator earns its keep for population work: diff --git a/wiki/core/api/plotting.md b/wiki/core/api/plotting.md index 5802f3b..5140d73 100644 --- a/wiki/core/api/plotting.md +++ b/wiki/core/api/plotting.md @@ -36,7 +36,7 @@ sources: - scripts/ellipse/plot.py pinned_commit: d6db2643b9f2cd418efc9473f560dc2a2d459c73 last_updated: 2026-08-01 -content_sha256: 37dc0bdc3b53aadbdbb7bb61a84fa0c2967ebda1e7443cde326ef9dde553b9dd +content_sha256: 8448417f14908a000e38bb5285fdc5b1629745e291d7c9ce50274324ded63ba6 --- # Plotting @@ -97,8 +97,28 @@ output_format=None, lines=None, ax=None)`. ## Subplot functions -Every one takes `output_path`, `output_format` and (mostly) `output_filename` directly, plus -`colormap`, `use_log10` and `title_prefix`. +Every one takes `output_path` and `output_format` directly, plus `colormap`, `use_log10` and +`title_prefix`. **They do not share a naming kwarg**, and passing the wrong one raises +`TypeError`: + +| How the file is named | Functions | +|---|---| +| `output_filename` | `subplot_imaging_dataset`, `subplot_imaging_dataset_list`, `subplot_interferometer_dataset`, `subplot_interferometer_dirty_images` (plus `plot_array` and `plot_grid` above) | +| `auto_filename` (default `"galaxies"`) | `subplot_galaxies` | +| **fixed stem**, no naming kwarg | everything else | + +The fixed stems, read off `_save_subplot(...)` in each function: +`subplot_fit_imaging` and `subplot_fit_interferometer` → `fit`; +`subplot_fit_imaging_of_galaxy` → `of_galaxy_`; +`subplot_galaxy_images` → `galaxy_images`; `subplot_galaxy_light_profiles` → `image`; +`subplot_galaxy_mass_profiles` → the quantity it drew (`convergence`, `potential`, +`deflections_y`, `deflections_x` — one file per enabled flag); +`subplot_basis_image` → `basis_image`; `subplot_adapt_images` → `adapt_images`; +`subplot_fit_real_space` → `fit_real_space`; `subplot_fit_dirty_images` → `fit_dirty_images`; +`subplot_fit_ellipse` → `fit_ellipse`; `subplot_ellipse_errors` → `ellipse_errors`. + +So for a fixed-stem subplot the **directory** is what distinguishes one context from another — +give each fit, or each variant of a figure, its own folder. ### Datasets @@ -226,8 +246,10 @@ aplt.plot_array( With no `output_path`, a figure is displayed. Pass `output_path` (a directory) and `output_format` and it is written as `{output_path}/{title}.{output_format}`; pass -`output_filename` to name it explicitly. `output_format` also accepts a **list**, which writes the -same figure in each format at once: +`output_filename` to name it explicitly — **on the functions that accept it**, which is +`plot_array`, `plot_grid` and the five subplots tabulated under "Subplot functions" above. Every +other subplot writes a fixed stem, so name the *directory* instead. `output_format` also accepts +a **list**, which writes the same figure in each format at once: ```python aplt.plot_array( diff --git a/wiki/core/concepts/galaxies.md b/wiki/core/concepts/galaxies.md index 8ee79e7..8d28e27 100644 --- a/wiki/core/concepts/galaxies.md +++ b/wiki/core/concepts/galaxies.md @@ -15,7 +15,7 @@ sources: - scripts/guides/units/cosmology.py pinned_commit: d6db2643b9f2cd418efc9473f560dc2a2d459c73 last_updated: 2026-08-01 -content_sha256: 8d74a683b7ed668ad6d16068e0fa6405dd235d7166af5b45da00639f534ad308 +content_sha256: 0144d127e9cc512639375a9bea805e10afbb4527979a54004d25d39e35ab5ab1 --- # Galaxy and Galaxies @@ -187,7 +187,9 @@ so its profiles carry the solved intensities. [`linear_light_profiles_and_mge`](./linear_light_profiles_and_mge.md). The maximum-likelihood galaxies are also written to `files/galaxies.json` in the output -folder and reload with `ag.from_json(...)`, behaving identically to the in-memory object. +folder and reload with `ag.from_json(...)` — which returns a plain `list` of `Galaxy` objects, +**not** an `ag.Galaxies`. Wrap it (`ag.Galaxies(galaxies=galaxies)`) before calling any of the +collection methods below. ## Plotting diff --git a/wiki/core/concepts/non_linear_search.md b/wiki/core/concepts/non_linear_search.md index fb43f72..852eb11 100644 --- a/wiki/core/concepts/non_linear_search.md +++ b/wiki/core/concepts/non_linear_search.md @@ -18,7 +18,7 @@ sources: - scripts/guides/results/start_here.py pinned_commit: d6db2643b9f2cd418efc9473f560dc2a2d459c73 last_updated: 2026-08-01 -content_sha256: af178e766c1577fe0aed53c2c0ae8ea8cf664a8147d755b2f8fbde44782810fb +content_sha256: 13439cc6b908483c12c5b907e921b9491bd3ee179fc328654c130e2ac1048fce --- # The non-linear search @@ -222,7 +222,7 @@ search starts, not when it finishes: ``` files/ JSON + CSV, loadable Python objects - galaxies.json max log likelihood Galaxies + galaxies.json max log likelihood galaxies (loads as a list) model.json the fitted af.Collection samples.csv every accepted sample samples_summary.json max log likelihood values + errors diff --git a/wiki/core/concepts/samples_and_posteriors.md b/wiki/core/concepts/samples_and_posteriors.md index 0aab3c6..75ceeee 100644 --- a/wiki/core/concepts/samples_and_posteriors.md +++ b/wiki/core/concepts/samples_and_posteriors.md @@ -20,7 +20,7 @@ sources: - scripts/guides/units/flux.py pinned_commit: d6db2643b9f2cd418efc9473f560dc2a2d459c73 last_updated: 2026-08-01 -content_sha256: e9dc51e8b79c37b55ba7cdee9c2ac4fd76fa00bb1783064cb3584c76314ab1b2 +content_sha256: 927b725b1dde11c399bf9f2925683f06a21131d8496e57f6e315e074c46ae89a --- # Samples and posteriors @@ -216,8 +216,12 @@ print(samples.max_log_likelihood()) ``` `autogalaxy_workspace:scripts/guides/results/start_here.py`. `search.paths.output_path` saves -you from constructing the `` by hand. The loaded objects behave identically to the -in-memory ones — `galaxies.json` in particular already carries solved intensities. +you from constructing the `` by hand. `galaxies.json` already carries the solved +intensities of any linear profile, so the profiles themselves are as good as the in-memory ones +— but **it deserialises to a plain `list` of `Galaxy` objects, not an `ag.Galaxies`**. Indexing +(`galaxies[0].bulge`) and `ag.FitImaging(dataset=…, galaxies=galaxies)` work on the list; +collection methods like `image_2d_from` do not. Re-wrap for those: +`galaxies = ag.Galaxies(galaxies=galaxies)`. ## Many fits — the aggregator diff --git a/wiki/core/external/skill_citation_map.md b/wiki/core/external/skill_citation_map.md index 82cbc6f..6ed20eb 100644 --- a/wiki/core/external/skill_citation_map.md +++ b/wiki/core/external/skill_citation_map.md @@ -9,7 +9,7 @@ sources: - PENDING.md pinned_commit: ed72fabb33e14a9a701a4d280e8775dd3a20e98c last_updated: 2026-08-01 -content_sha256: 8dc3621b0d65e14454bf56d459923b9d7a246c7087b124c157ec5f51c8f71fe5 +content_sha256: 1d29ca6a67048d5e548dd83006597944a348f31ade73b852b625800a64d1f40b --- # Skill → external resource map @@ -53,18 +53,32 @@ Default to `.ipynb` for student-leaning users and `.py` for returning PyAutoGala | `ag_audit_skill_apis` | _ | _ | _ | | `ag_update_wiki` | _ | _ | _ | | `ag_refresh_api_docs` | _ | _ | _ | - -Every cell is `_` today, and that is the correct state rather than an omission: all seven -skills are **internal to the workspace** — two meta-skills, two repo-workflow skills, and +| `ag_setup_environment` | `chapter_1_introduction/tutorial_0_visualization` | `installation/overview` | `imaging/start_here.py` | +| `ag_prepare_imaging_data` | `chapter_1_introduction/tutorial_2_data` | `overview/overview_2_new_user_guide` | `imaging/data_preparation/start_here.py` | +| `ag_simulate_dataset` | `chapter_1_introduction/tutorial_1_grids_and_galaxies` | `overview/overview_1_start_here` | `imaging/simulator.py` | +| `ag_build_imaging_model` | `chapter_2_modeling/tutorial_3_realism_and_complexity` | `general/model_cookbook` | `imaging/modeling.py` | +| `ag_configure_search` | `chapter_2_modeling/tutorial_1_non_linear_search` | `general/configs` | `guides/modeling/searches.py` | +| `ag_run_search` | `chapter_2_modeling/tutorial_2_practicalities` | `overview/overview_2_new_user_guide` | `imaging/start_here.py` | +| `ag_plot_fit` | `chapter_1_introduction/tutorial_3_fitting` | `api/plot` | `guides/plot/start_here.py` | +| `ag_load_results` | `chapter_2_modeling/tutorial_7_results` | `api/fitting` | `guides/results/start_here.py` | +| `ag_debug_fit_failure` | `chapter_2_modeling/tutorial_4_dealing_with_failure` | `general/likelihood_function` | `guides/modeling/bug_fix.py` | + +The first seven rows are entirely `_`, and that is the correct state rather than an omission: +those skills are **internal to the workspace** — two meta-skills, two repo-workflow skills, and the three maintenance skills that audit and refresh this repo's own content. None of them teaches a galaxy-modelling task, so none has an external resource that would help a user, and per `_style.md` a skill whose row is entirely `_` omits the `## Further reading` block altogether. -The first rows with real cells arrive with the modelling skills in a later phase — a -data-preparation skill will cite `imaging/data_preparation/start_here.py`, a model-building -skill will cite `general/model_cookbook` and `imaging/start_here.py`, and so on. Author the -row and the skill in the same change, from the grounding script named in +The nine core-loop rows beneath them are the first with real cells, one per skill in the +modelling loop. Each skill's `## Further reading` block was written from its row here, and +each cell was confirmed to resolve against the target repo's catalogue (or the RTD page map) +before it was recorded. A skill may cite an *extra* tutorial inline within a bullet — the +masking lecture under `ag_prepare_imaging_data`, the optional searches chapter under +`ag_configure_search` — but the row holds the one primary cell per audience. + +Rows for the Phase-4b feature skills arrive with those skills. Author the row and the skill in +the same change, from the grounding script named in [`../../../PENDING.md`](../../../PENDING.md), never from memory. ## Template for the inserted skill block diff --git a/wiki/core/operations/sandbox.md b/wiki/core/operations/sandbox.md index 7cbe5ba..61df166 100644 --- a/wiki/core/operations/sandbox.md +++ b/wiki/core/operations/sandbox.md @@ -32,7 +32,7 @@ sources: - .claude/hooks/validate_pyauto_code.py pinned_commit: ed72fabb33e14a9a701a4d280e8775dd3a20e98c last_updated: 2026-08-01 -content_sha256: e594c604f396ad7dac62d4f787b78db19c8b20981cf379ed0a5f8b82e6043b06 +content_sha256: f31b1ecc8b237bfe50ccaa550be6e93cab490cadf7448f2c9a72cc937e0b426a --- # Sandbox / restricted-environment configuration @@ -111,8 +111,16 @@ Levels 2 and 3 are the ones to reach for while iterating on a script's *structur level** — the parameter values are whatever the truncated or mocked run produced, so never quote a structural parameter measured in test mode. -Two consequences of test mode are easy to trip over: +Three consequences of test mode are easy to trip over: +- **A level-2 or level-3 bypass writes no fit products.** The bypass path + (`PyAutoFit:autofit/non_linear/search/abstract_search.py` `_fit_bypass_test_mode`) saves + `files/samples.csv` and `files/samples_summary.json` and marks the fit complete, but it never + runs `perform_update` or `analysis.save_results`, so there is **no `model.results`, no + `model.info`, and nothing under `image/`** — no `dataset.png`, no `fit.png`, no `fit.fits`. + Level 1 goes through the normal path and *does* write all of them. So use level 1, not 2, when + the thing you are checking is the output folder or a figure; use 2 or 3 when you are checking + that the script runs at all. - **Output is namespaced.** Any active level inserts a `test_mode` segment directly after the output root, so results land in `output/test_mode///` rather than beside real runs (`PyAutoFit:autofit/non_linear/paths/abstract.py`). This is deliberate: