diff --git a/.claude/skills/ag_basis_profiles.md b/.claude/skills/ag_basis_profiles.md new file mode 120000 index 0000000..f0dfeaf --- /dev/null +++ b/.claude/skills/ag_basis_profiles.md @@ -0,0 +1 @@ +../../skills/ag_basis_profiles.md \ No newline at end of file diff --git a/.claude/skills/ag_build_interferometer_model.md b/.claude/skills/ag_build_interferometer_model.md new file mode 120000 index 0000000..8913818 --- /dev/null +++ b/.claude/skills/ag_build_interferometer_model.md @@ -0,0 +1 @@ +../../skills/ag_build_interferometer_model.md \ No newline at end of file diff --git a/.claude/skills/ag_chain_searches.md b/.claude/skills/ag_chain_searches.md new file mode 120000 index 0000000..a3997ce --- /dev/null +++ b/.claude/skills/ag_chain_searches.md @@ -0,0 +1 @@ +../../skills/ag_chain_searches.md \ No newline at end of file diff --git a/.claude/skills/ag_ellipse_fitting.md b/.claude/skills/ag_ellipse_fitting.md new file mode 120000 index 0000000..751dbb9 --- /dev/null +++ b/.claude/skills/ag_ellipse_fitting.md @@ -0,0 +1 @@ +../../skills/ag_ellipse_fitting.md \ No newline at end of file diff --git a/.claude/skills/ag_light_model_extras.md b/.claude/skills/ag_light_model_extras.md new file mode 120000 index 0000000..fb0b965 --- /dev/null +++ b/.claude/skills/ag_light_model_extras.md @@ -0,0 +1 @@ +../../skills/ag_light_model_extras.md \ No newline at end of file diff --git a/.claude/skills/ag_multi_dataset.md b/.claude/skills/ag_multi_dataset.md new file mode 120000 index 0000000..d2c95ce --- /dev/null +++ b/.claude/skills/ag_multi_dataset.md @@ -0,0 +1 @@ +../../skills/ag_multi_dataset.md \ No newline at end of file diff --git a/.claude/skills/ag_multi_galaxy_and_cluster.md b/.claude/skills/ag_multi_galaxy_and_cluster.md new file mode 120000 index 0000000..258d3dd --- /dev/null +++ b/.claude/skills/ag_multi_galaxy_and_cluster.md @@ -0,0 +1 @@ +../../skills/ag_multi_galaxy_and_cluster.md \ No newline at end of file diff --git a/.claude/skills/ag_pixelization.md b/.claude/skills/ag_pixelization.md new file mode 120000 index 0000000..494fd57 --- /dev/null +++ b/.claude/skills/ag_pixelization.md @@ -0,0 +1 @@ +../../skills/ag_pixelization.md \ No newline at end of file diff --git a/AGENTS.md b/AGENTS.md index 6ec92f1..3833f5a 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -85,11 +85,16 @@ Map every request onto one or more layers: `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. + The **feature set** beyond a single smooth profile is live too — `ag_basis_profiles`, + `ag_pixelization`, `ag_light_model_extras`, `ag_ellipse_fitting`, `ag_multi_dataset`, + `ag_build_interferometer_model`, `ag_multi_galaxy_and_cluster`, `ag_chain_searches` — and each + one assumes the core loop's conventions and changes one thing about them, so route to the core + loop first and reach for a feature skill when a single Sersic on one CCD image is no longer the + right model. 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 twenty-four 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, how a pixelised reconstruction is regularised. diff --git a/PENDING.md b/PENDING.md index fba407d..e985ce3 100644 --- a/PENDING.md +++ b/PENDING.md @@ -110,25 +110,39 @@ workspace scripts directly. ## Phase 4b — features beyond a single smooth profile (8 skills) -- [ ] `skills/ag_basis_profiles.md` — grounding: `imaging/features/linear_light_profiles/`, - `imaging/features/multi_gaussian_expansion/`, `imaging/features/shapelets/`. -- [ ] `skills/ag_pixelization.md` — grounding: `imaging/features/pixelization/` (including - `galaxy_reconstruction.py` and `likelihood_function.py`). -- [ ] `skills/ag_light_model_extras.md` — grounding: `imaging/features/extra_galaxies/`, - `imaging/features/sky_background/`, `imaging/features/operated_light_profile/`. -- [ ] `skills/ag_ellipse_fitting.md` — grounding: `ellipse/modeling.py`, - `ellipse/multipoles.py`, `ellipse/database.py`. **Route to `modeling.py`: `ellipse/` has - no `start_here.py`.** -- [ ] `skills/ag_multi_dataset.md` — grounding: `multi_dataset/start_here.py`, - `multi_dataset/features/`. The idiom deny-list genuinely bites here: datasets combine - via the factor graph, never by summing analyses. -- [ ] `skills/ag_build_interferometer_model.md` — grounding: `interferometer/start_here.py`, - `interferometer/modeling.py`, `interferometer/features/`. -- [ ] `skills/ag_multi_galaxy_and_cluster.md` — grounding: `multi_galaxy/start_here.py`, - `cluster/start_here.py`. The subject is member **light**, not lensing — the phase where - lensing language most wants to creep back in. -- [ ] `skills/ag_chain_searches.md` — grounding: `guides/modeling/chaining.py`, HowToGalaxy - `chapter_3_search_chaining`. +**Delivered.** All eight feature 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_basis_profiles`, `ag_pixelization`, `ag_light_model_extras`, `ag_ellipse_fitting`, +`ag_multi_dataset`, `ag_build_interferometer_model`, `ag_multi_galaxy_and_cluster`, +`ag_chain_searches`. Each was written from the grounding scripts named above rather than from +memory, and every symbol resolves against released `2026.7.29.2`. That brings the repo to +twenty-four live skills, seventeen of them galaxy-modelling. + +Three grounding notes, recorded so nobody re-derives them: + +- `ellipse/` really has no `start_here.py`, so `ag_ellipse_fitting` routes to `modeling.py` and + says so in its own `## Further reading` block. +- The lecture series has no chapter on ellipse fitting, interferometry or multi-wavelength + fitting. Rather than invent a citation, `ag_ellipse_fitting` omits its student bullet and + points at `wiki/core/concepts/ellipse_fitting_and_multipoles.md` instead, while + `ag_multi_dataset` and `ag_build_interferometer_model` cite the HowToGalaxy tutorial that + teaches the idea each fit leans on hardest (linear profiles; the likelihood) and say plainly + that it is not a chapter on their own subject. +- Grounding `ag_chain_searches` against the wheel exposed a **stale claim** that had spread from + the workspace script into two merged wiki pages: `result.model` returns the fitted model with + its **original priors unchanged** (`samples_summary.model.mapper_via_defaults_from`, which maps + every prior to itself), *not* narrowed `TruncatedGaussianPrior`s. The narrowing lives on + `result.model_centred` and its `model_centred_absolute(a=)` / `model_centred_relative(r=)` / + `model_centred_max_lh_bounded(b=)` variants. `wiki/core/concepts/non_linear_search.md` and + `wiki/core/api/configuration.md` were corrected in this phase and re-stamped. + +### Still open + +- [ ] Upstream fix for the same claim in + `autogalaxy_workspace:scripts/guides/modeling/chaining.py`, whose prose still describes + `result.model` as producing narrowed Gaussians. A `contribute-upstream` candidate; until it + lands, `skills/ag_chain_searches.md` warns the reader that the script's description is out + of date. ## Phase 5 — `wiki/literature` corpus diff --git a/README.md b/README.md index 5fe5c26..4c36a18 100644 --- a/README.md +++ b/README.md @@ -158,8 +158,12 @@ No file in this repository links to something that does not exist yet. - **The core modelling loop — nine skills.** Environment setup, imaging data preparation, dataset simulation, model building, search configuration, running the fit, plotting the fit, loading results, and debugging a failed fit. A galaxy-science request routes into - these. Alongside them sit two meta-skills, two project-workflow skills and three - repository-maintenance skills — sixteen in total, catalogued in + these. +- **The features beyond a single smooth profile — eight skills.** Basis profiles and MGE, + pixelised reconstruction, extra galaxies / sky / operated profiles, ellipse fitting, + multi-dataset fits, interferometer modelling, multi-galaxy and cluster fields, and search + chaining. Alongside both sets sit two meta-skills, two project-workflow skills and three + repository-maintenance skills — twenty-four in total, catalogued in [`skills/README.md`](skills/README.md). - **The curated reference wiki**, [`wiki/core/`](wiki/core/index.md) — 37 pages across `stack/`, `api/`, `concepts/`, `operations/` and `external/`. Every page pins the source @@ -172,10 +176,6 @@ No file in this repository links to something that does not exist yet. **Still being written:** -- **Feature skills beyond a single smooth profile** (Phase 4b): basis profiles and MGE, - pixelisations, ellipse fitting, multi-dataset and interferometer modelling, multi-galaxy - and cluster fields, and search chaining. The underlying PyAutoGalaxy features all exist and - the reference wiki covers them; what is pending is the dedicated procedural skill for each. - **The literature wiki** (Phase 5) — a galaxy-structure science reference with cross-linked concept, entity and source pages. Until it lands, the assistant will not cite it and must not invent its contents. diff --git a/benchmarks/prompts/easy_cosj100020_imaging.md b/benchmarks/prompts/easy_cosj100020_imaging.md index b3cfdc9..bf836b3 100644 --- a/benchmarks/prompts/easy_cosj100020_imaging.md +++ b/benchmarks/prompts/easy_cosj100020_imaging.md @@ -19,14 +19,15 @@ running the fit, plotting it and reading the result are all covered by existing measures whether it *finds and follows* the built-in workflow rather than writing PyAutoGalaxy from memory. -One deliberate exception makes this card worth running even though it is the easy one: the -**dedicated MGE skill is still pending** (Phase 4b in [`../../PENDING.md`](../../PENDING.md)). -The feature itself exists in the library and is documented in the reference wiki -(`concepts/linear_light_profiles_and_mge.md`, `api/light_profile_catalog.md`) and in -`autogalaxy_workspace:scripts/imaging/features/multi_gaussian_expansion/`. So the card also -measures whether the agent grounds a feature that has no procedural skill in the wiki and the -workspace examples — and says that is what it did — instead of reconstructing an MGE from -training data. +The MGE half of the prompt is now covered too: `skills/ag_basis_profiles.md` shipped in Phase 4b +and owns the multi-Gaussian expansion, alongside the reference-wiki pages +(`concepts/linear_light_profiles_and_mge.md`, `api/light_profile_catalog.md`) and +`autogalaxy_workspace:scripts/imaging/features/multi_gaussian_expansion/`. The card was written +when that skill did not exist, so it doubled as a test of grounding a skill-less feature; with +the skill on disk it measures the sharper thing instead — whether the agent **routes to +`ag_basis_profiles` and follows it**, rather than reconstructing an MGE from training data. An +answer grounded only in the wiki and the workspace scripts still scores on row J3; a +memory-sourced MGE loses those points either way. ## Prompt @@ -78,7 +79,7 @@ stay identical (a divergence is a bug — fix the README or bump this card's `ve |---|-----------|-----| | J1 | Real-data gate honoured: dataset plotted and inspected, and **both** questions settled before any fit — contaminants (the faint 2.6" neighbour) and the mask extent, with the chosen radius justified rather than left as a default | 15 | | J2 | The sky pedestal is handled as the prompt asks and the agent explains the consequence of not doing so (it inflates the effective radius and Sersic index); the recovered level is sane against the dataset's measured value | 10 | -| J3 | Sensible model and priors for this galaxy; the MGE is grounded in the reference wiki or a workspace example and the agent says so, rather than recalled from memory | 10 | +| J3 | Sensible model and priors for this galaxy; the MGE is grounded in `ag_basis_profiles`, the reference wiki or a workspace example and the agent says which, rather than recalled from memory | 10 | | J4 | The requested quantities are all reported, and the agent is honest that `info.json`'s `effective_radius_arcsec_rough` is a prior-scale measurement rather than ground truth to be matched | 10 | | J5 | The multi-band comparison is scientifically framed — the size trend with wavelength separated from the changing PSF width — and the shipped model PSF is acknowledged as the dominant systematic | 10 | | J6 | Conduct: concise assistant-mode communication, no fabricated numbers, API-gate discipline (no invented symbols, functional `aplt` plotting) | 5 | diff --git a/llms.txt b/llms.txt index 8ade08e..d041cda 100644 --- a/llms.txt +++ b/llms.txt @@ -20,7 +20,7 @@ Read order: [AGENTS.md](./AGENTS.md) → [skills/README.md](./skills/README.md) **⚠️ 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.** 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 rest of the benchmark suite (one of four cards is written) 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.** Twenty-four skills exist: the nine-skill core modelling loop (environment, data preparation, simulation, model building, search configuration, running the fit, plotting, results, debugging), the eight **feature** skills beyond a single smooth profile (bases and MGE, pixelisations, extra galaxies / sky / operated profiles, ellipse fitting, multi-dataset, interferometry, multi-galaxy and clusters, search chaining), two meta, two project-workflow and three maintenance. The literature wiki and the rest of the benchmark suite (one of four cards is written) 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) diff --git a/skills/README.md b/skills/README.md index 9eb619c..c619980 100644 --- a/skills/README.md +++ b/skills/README.md @@ -27,16 +27,16 @@ configured) via symlinks; the canonical files live here. ## Index — what exists today -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 +Twenty-four skills are written: nine for the core modelling loop, eight for the features beyond +a single smooth profile, 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 and the feature set are both live. For a request that falls outside the +seventeen galaxy-modelling skills 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). ### Galaxy modelling — the core loop @@ -84,6 +84,52 @@ Read in this order for an end-to-end fit; each is usable on its own. Every one i the two silent failures: a resumed fit whose identifier ignored the data, and a cached result mistaken for a new one. +### Features beyond a single smooth profile + +Reach for one of these once the core loop is running and a single Sersic — or a plain +bulge-plus-disk on one CCD image — is no longer the right model. Each assumes the core loop's +composition, search and plotting conventions and changes one thing about them. + +- [`ag_basis_profiles.md`](./ag_basis_profiles.md) — fit morphology with a *basis* rather than one + or two smooth profiles: linear light profiles (`ag.lp_linear`), a Multi-Gaussian Expansion or a + shapelet expansion, where every component's `intensity` is solved analytically by a linear + inversion instead of sampled by the search — with the positive-only versus signed solver, the + compact nuclear basis, and how to read solved intensities back out of a fit. +- [`ag_pixelization.md`](./ag_pixelization.md) — reconstruct clumpy or irregular light directly on + a regularized pixel mesh, using `ag.Pixelization` alongside a parametric `ag.lp_linear` bulge: + choosing an `ag.mesh` and an `ag.reg` scheme and what each costs, why `mesh_shape` must be fixed + before the fit, `over_sample_size_pixelization`, noise scaling instead of hard masking, and + reading the reconstruction and its evidence terms out of the `Inversion`. +- [`ag_light_model_extras.md`](./ag_light_model_extras.md) — the three components that sit beside a + galaxy's own light profiles: contaminating extra galaxies, a residual background sky via + `ag.DatasetModel`, and operated (already-PSF-convolved) profiles for compact nuclear emission — + worked against the bundled real dataset, which has both an un-subtracted sky and a faint + neighbour 2.6" out. +- [`ag_ellipse_fitting.md`](./ag_ellipse_fitting.md) — measure morphology non-parametrically by + fitting isophotes with `ag.Ellipse` / `ag.FitEllipse` / `ag.AnalysisEllipse`, one ellipse at a + time at fixed `major_axis`, producing radial axis-ratio and position-angle profiles instead of a + light-profile model — plus `ag.EllipseMultipole`, `af.Drawer` and the `ag.agg` classes that read + many fits back. +- [`ag_multi_dataset.md`](./ag_multi_dataset.md) — fit several datasets of the same galaxy jointly + (multi-wavelength bands, repeated exposures, imaging together with visibilities) through the + `af.AnalysisFactor` + `af.FactorGraphModel` construction that is the only way to combine them: + what is shared versus freed per dataset, a wavelength relation via prior arithmetic, and + astrometric offsets. +- [`ag_build_interferometer_model.md`](./ag_build_interferometer_model.md) — model a galaxy + observed with a radio or millimetre interferometer by fitting its complex visibilities in the + uv-plane: loading an `ag.Interferometer` against a real-space mask, choosing the transformer by + visibility count, why there is no PSF and no over-sampling, and reading dirty images as + diagnostics rather than data. +- [`ag_multi_galaxy_and_cluster.md`](./ag_multi_galaxy_and_cluster.md) — model several galaxies + whose light blends on the sky: an interacting or projected pair, a compact multiple, or a cluster + field with a brightest cluster galaxy plus a catalogue-driven member tier + (`ag.galaxy_table_from_csv`) whose intensities tie to one shared normalization, so population + size costs no dimensions — plus per-galaxy decomposed photometry. +- [`ag_chain_searches.md`](./ag_chain_searches.md) — break one hard fit into a sequence of easier + searches, each initialized from the last: `result.model_centred` and its absolute / relative / + bounded variants to narrow priors, `result.model` to keep the original ones, `result.instance` to + fix a component outright, and the output-path convention that keeps a chain's searches together. + ### Meta - [`_style.md`](./_style.md) — writing guide every skill is authored against. Read first @@ -137,19 +183,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 4b — features beyond a single smooth profile - -| Skill | Purpose | Grounding (`autogalaxy_workspace/scripts/`) | -|-------|---------|--------------------------------------------| -| `ag_basis_profiles` | linear light profiles, Multi-Gaussian Expansion and shapelets — flexible bases that solve for intensity by linear inversion | `imaging/features/linear_light_profiles/`, `imaging/features/multi_gaussian_expansion/`, `imaging/features/shapelets/` | -| `ag_pixelization` | pixelised reconstruction of an irregular or clumpy galaxy, with regularisation | `imaging/features/pixelization/` | -| `ag_light_model_extras` | blended neighbours, sky background, and already-PSF-convolved (operated) components | `imaging/features/extra_galaxies/`, `imaging/features/sky_background/`, `imaging/features/operated_light_profile/` | -| `ag_ellipse_fitting` | non-parametric isophote fitting and multipole perturbations | `ellipse/modeling.py`, `ellipse/multipoles.py`, `ellipse/database.py` (note: `ellipse/` has no `start_here.py`) | -| `ag_multi_dataset` | simultaneous fits across wavebands or instruments via the factor graph | `multi_dataset/start_here.py`, `multi_dataset/features/` | -| `ag_build_interferometer_model` | uv-plane modelling of ALMA / JVLA observations | `interferometer/start_here.py`, `interferometer/modeling.py` | -| `ag_multi_galaxy_and_cluster` | 2+ blended galaxies each with a free light model; BCG + catalogue-driven member population in a cluster field | `multi_galaxy/start_here.py`, `cluster/start_here.py` | -| `ag_chain_searches` | sequence searches so a later fit inherits priors from an earlier one | `guides/modeling/chaining.py`, HowToGalaxy `chapter_3_search_chaining` | - ### Phase 5 — literature | Skill | Purpose | Grounding | diff --git a/skills/ag_basis_profiles.md b/skills/ag_basis_profiles.md new file mode 100644 index 0000000..35668e0 --- /dev/null +++ b/skills/ag_basis_profiles.md @@ -0,0 +1,602 @@ +--- +name: ag_basis_profiles +description: Fit a galaxy's morphology with a *basis* rather than one or two smooth profiles — linear light profiles (`ag.lp_linear`), a Multi-Gaussian Expansion, or a shapelet expansion — where every component's `intensity` is solved analytically by a linear inversion instead of sampled by the search. Covers why `ag.lp_linear` is the default for any bulge-disk decomposition, when an MGE beats a Sersic pair (isophotal twists, radially varying ellipticity, asymmetry), when shapelets beat an MGE, the positive-only versus signed linear solver and what a negative amplitude means physically, `ag.model_util.mge_model_from` and the hand-rolled `ag.lp_basis.Basis` idiom behind it, a compact MGE for an unresolved nucleus, and how to read solved intensities back out of a fit. Use when a Sersic or bulge-plus-disk model leaves structure in the residuals but the galaxy still has a well-defined centre. Not for pixelized reconstruction of genuinely irregular clumps (`ag_pixelization`), not for non-parametric isophote measurement (`ag_ellipse_fitting`), and not for the basics of composing an `af.Model` tree (`ag_build_imaging_model`). +--- + +# Fitting morphology with a basis + +A Sersic profile makes a strong claim: that the galaxy's isophotes are concentric, similar +ellipses with a single fixed axis ratio and position angle, and that its surface brightness +falls off as one power of radius. Real galaxies routinely violate all of that. Isophotes twist +with radius, the ellipticity varies from the bulge outwards, a bar or a spiral pattern breaks +the symmetry — and when a Sersic cannot represent those features it does not fail loudly. It +absorbs them, and returns an `effective_radius` and `sersic_index` biased by exactly the +structure it could not fit. + +A **basis** is the answer: instead of one profile with a shape, use tens of simple profiles +whose superposition can take almost any shape. The reason this is affordable — and the reason +it is often *cheaper* than a Sersic pair — is statistical. Every basis component's `intensity` +is a **linear** parameter: given the shapes and positions of the components, the amplitudes +that maximise the likelihood follow from solving a linear system, exactly, in one step. So they +never enter the non-linear search's parameter space at all. A 60-Gaussian MGE can have four +free parameters where a two-component linear Sersic decomposition has nine. + +Three tiers, same machinery, increasing flexibility: + +| Tier | What it is | Free parameters | Reach for it when | +|---|---|---|---| +| `ag.lp_linear.*` | ordinary profiles with `intensity` solved | 6 per Sersic, 3 per aligned Exponential | you want interpretable bulge/disk parameters — **the default** | +| MGE (`ag.lp_basis.Basis` of `Gaussian`s) | 15–100 Gaussians on a fixed log-spaced `sigma` ladder | 4 (one basis) to 6 (two) | isophotal twists, radial ellipticity variation, mild asymmetry | +| shapelets | an orthonormal polar or Cartesian basis | 3 (spherical) to 5 (elliptical) | disky, star-forming morphology the MGE still smooths over | + +The science this serves is the decomposition itself. A bulge-plus-disk fit measures the +bulge-to-total light ratio, the two components' sizes and the bulge's concentration — the +quantities that place a galaxy on the mass–size relation and separate a classical spheroid from +a pseudo-bulge. A basis fit measures something the parametric fit cannot: how the shape of the +light changes with radius. Both concepts are in +[`../wiki/core/concepts/linear_light_profiles_and_mge.md`](../wiki/core/concepts/linear_light_profiles_and_mge.md) +and [`../wiki/core/concepts/shapelets.md`](../wiki/core/concepts/shapelets.md); the catalogue of +which variant lives in which module is +[`../wiki/core/api/light_profile_catalog.md`](../wiki/core/api/light_profile_catalog.md). + +This skill assumes you already have a masked, over-sampled dataset and know how to build an +`af.Model` tree — that is [`ag_build_imaging_model`](./ag_build_imaging_model.md), whose MGE +branch is the one-paragraph version of what follows. + +## Ask + +- *"Do you need the Sersic parameters, or the best fit?"* This is the real fork. If your science + is a bulge-to-total ratio, an effective radius, a Sersic index — a number you will put in a + table next to other people's — stay with `ag.lp_linear` Sersic-family profiles, because an + MGE does not hand you those numbers. If your science is the light distribution itself, or you + need the residuals flat before you can trust anything else, go to a basis. +- *"What does the residual map look like after a single Sersic?"* Structured residuals that + follow the isophotes (a four-lobed or twisted pattern) say the *shape* is wrong — an MGE + fixes that. Residuals in discrete off-centre lumps say the shape is fine but there is + substructure — that is a pixelization ([`ag_pixelization`](./ag_pixelization.md)), not a basis. + If you have not looked yet, look first; it is one plot and it decides the branch. +- *"Does the galaxy have a single well-defined centre?"* Every basis in this skill shares one + centre across all its components. A merger, or a galaxy with a bright companion overlapping + it, breaks that assumption — handle the companion first with + [`ag_light_model_extras`](./ag_light_model_extras.md). +- *"Is there a compact nuclear source?"* An AGN or a nuclear starburst is a separate component, + and there is a purpose-built compact MGE for it (branch four below). + +## Branch — linear light profiles, the default for any decomposition + +Start here even when you expect to end up with a basis. This is the same bulge-plus-disk model +`ag_build_imaging_model` composes, and the point of this branch is what happens to `intensity`. + +```python +""" +Galaxy Structure: Linear Bulge-Disk Decomposition +================================================= + +Decompose a galaxy's light into a Sersic bulge and an Exponential disk using *linear* light +profiles, whose `intensity` parameters are solved analytically by a linear inversion inside +every likelihood evaluation rather than sampled by the non-linear search. The ratio of the two +solved intensities, integrated, is the bulge-to-total light ratio. + +__Contents__ + +- **Imports:** Import the required libraries. +- **Dataset:** Load the masked, over-sampled imaging dataset. +- **Model:** Compose the linear bulge and disk, with their centres paired. +- **Check:** Confirm the parameter count and evaluate one likelihood. +- **Intensities:** Read the solved intensities back out of a completed fit. +""" + +""" +__Imports__ +""" +from pathlib import Path + +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +DATASET_PATH = Path("dataset") / "imaging" / "my_galaxy" +PIXEL_SCALES = 0.1 +MASK_RADIUS = 3.0 + +""" +__Dataset__ + +Loaded from FITS, masked to a radius chosen by inspecting the data, and over-sampled adaptively +where the profile's intensity gradient is steep. All three decisions 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__ + +The decomposition is identifiable because the two components are *asymmetric* in what they are +allowed to be: the bulge's `sersic_index` is free, so the fit can find a concentrated spheroid +(n ~ 4), while the `Exponential` disk has n fixed at 1 by construction. Two free Sersics on one +galaxy are largely degenerate and the fit will trade light between them. + +`ag.lp_linear` rather than `ag.lp` is the only change from a standard composition, and it is +not a small one: each profile's `intensity` leaves the search's parameter space and is instead +solved by a linear inversion that always returns the amplitudes maximising the likelihood given +the other parameters. Two dimensions vanish, and with them the strong degeneracies between +`intensity` and `effective_radius` / `sersic_index` that a sampler finds hardest to map +(`PyAutoGalaxy:autogalaxy/profiles/light/linear/sersic.py`). + +`bulge.centre = disk.centre` asserts the components are concentric — physically reasonable for +most galaxies, and worth two parameters. +""" +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)) + +""" +__Check__ + +The likelihood at the prior medians costs one evaluation and proves the model, dataset and +analysis are mutually compatible. Pass `use_jax=False` while debugging — NumPy tracebacks are +far easier to read (`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}") +``` + +N = 9: the bulge's centre (2), two `ell_comps` pairs (4), two `effective_radius` values and the +bulge's `sersic_index` — with both intensities solved. Adapted from +`autogalaxy_workspace:scripts/imaging/features/linear_light_profiles/modeling.py`, which +reduces `n_live` to 75 for exactly this reason: a simpler parameter space needs fewer live +points to map. + +The cost is real but small. A linear likelihood evaluation is roughly three to five times +slower than a standard one, because the inversion has to be solved every time. The workspace +example measures ~0.05 s against ~0.01 s on a low-resolution dataset. That is repaid by fewer +iterations and a more reliable posterior, and the trade tips further towards linear the more +profiles you add — which is what makes the next two branches possible at all. + +### The consequence: `intensity` is not in the model, so it is not in the results + +This trips up everyone once. A linear profile has no `intensity` attribute to report, so +`model.results` does not list one, and **a linear profile cannot be plotted directly** — there +is no amplitude to evaluate. The intensities exist only after a fit, as the inversion's +solution. Three ways to get at them, all from +`autogalaxy_workspace:scripts/imaging/features/linear_light_profiles/modeling.py`: + +```python +""" +__Intensities__ + +`max_log_likelihood_galaxies` has already performed the inversion, so its profiles carry their +solved `intensity` values and can be evaluated and plotted like ordinary profiles. The +`Fit` additionally exposes a dictionary keyed by the profile objects themselves, which is the +unambiguous route when several components share a class +(`PyAutoGalaxy:autogalaxy/imaging/fit_imaging.py`). +""" +galaxies = result.max_log_likelihood_galaxies + +print(f"bulge intensity = {galaxies[0].bulge.intensity}") + +fit = result.max_log_likelihood_fit + +print(fit.linear_light_profile_intensity_dict[fit.galaxies[0].bulge]) +print(fit.linear_light_profile_intensity_dict[fit.galaxies[0].disk]) + +galaxies_ordinary = fit.model_obj_linear_light_profiles_to_light_profiles + +aplt.subplot_galaxies( + galaxies=galaxies_ordinary, + grid=dataset.grid, + output_path="scripts/scratch/my_galaxy/", + auto_filename="decomposition", + output_format="png", +) +``` + +`model_obj_linear_light_profiles_to_light_profiles` returns the same galaxies with every linear +profile replaced by its ordinary equivalent at the solved amplitude — that is the object to +hand to any plotting or derived-quantity call. `fit.inversion.linear_obj_list` holds one +`LightProfileLinearObjFuncList` per linear component (or per `Basis`), which is where the +bookkeeping lives when you need to be certain which amplitude belongs to which profile. + +Integrating those amplitudes over the two profiles is what gives you the bulge-to-total ratio; +the machinery for luminosities and magnitudes is +[`../wiki/core/concepts/cosmology_and_units.md`](../wiki/core/concepts/cosmology_and_units.md). +Ask if you want the derived-quantity route spelled out. + +## Branch — a Multi-Gaussian Expansion, when the *shape* is wrong + +A Gaussian is a poor galaxy profile on its own and an excellent basis function in company: +sum enough of them at fixed, logarithmically spaced widths and you can build any monotonic +radial profile, and — if you let separate groups of them carry separate ellipticities — an +isophote shape that changes with radius. That last property is the one a Sersic cannot have at +any parameter count, and it is why an MGE is the workspace's recommended default model. + +The composition is long, so the library ships a helper: + +```python +""" +__Model__ + +An MGE bulge: 20 linear Gaussians sharing one centre and one ellipticity, with every `sigma` +fixed on a logarithmic ladder running from below the pixel scale out to `mask_radius`, and +every intensity solved by the inversion. `mask_radius` sets the outer end of that ladder, so +pass the radius you actually masked at — a ladder that stops short of the data leaves the +outer isophotes unrepresented (`PyAutoGalaxy:autogalaxy/analysis/model_util.py`). +""" +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) +print(f"Total free parameters = {model.total_free_parameters}") +``` + +**N = 4** — a shared centre and a shared `ell_comps`. Twenty Gaussians, four dimensions, more +morphological freedom than a Sersic's six. Adapted from +`autogalaxy_workspace:scripts/imaging/start_here.py`, where the MGE is the recommended default. + +The full signature is +`mge_model_from(mask_radius, total_gaussians=30, gaussian_per_basis=1, +centre_prior_is_uniform=True, centre=(0.0, 0.0), centre_fixed=None, centre_per_basis=False, +centre_sigma=0.3, ell_comps_prior_is_uniform=False, ell_comps_uniform_width=0.2, +ell_comps_sigma=0.3, use_spherical=False)`. The arguments worth knowing: + +- **`gaussian_per_basis=2`** splits the ladder into two bases with *independent* ellipticities. + N goes 4 → 6 and you buy the ability to fit an isophotal twist: an inner group and an outer + group free to be differently elongated and differently oriented. This is the single most + useful knob here, and the reason to prefer an MGE over a Sersic pair for an early-type galaxy. +- **`centre_fixed=(y, x)`** pins the centre entirely, dropping two more parameters. Used for + contaminating neighbours ([`ag_light_model_extras`](./ag_light_model_extras.md)). +- **`use_spherical=True`** drops the ellipticity — for a genuinely round system, or as a + deliberately rigid comparison model. +- **`total_gaussians`** buys radial resolution, not shape freedom. 15 is enough for a smooth + profile; 60 for a well-resolved one; beyond that you are mostly paying for likelihood + evaluations. Cost scales with the count, because every Gaussian's image must be computed and + PSF-convolved: the workspace measures ~0.5 s per evaluation for 60 Gaussians, against ~0.01 s + for a standard Sersic. VRAM rises too — 10–50 MB per batched likelihood at 60 Gaussians, so + check `analysis.print_vram_use(model=model, batch_size=search.batch_size)` before a GPU run + with hundreds of components. + +That run time is why the overall comparison is not obvious, and the workspace is honest about +it: the MGE's likelihood is much slower, but its parameter space is so much simpler — and +crucially contains *no parameter that scales the galaxy's size* — that the search converges far +faster. Net, it usually wins, and it fits better. + +### The hand-rolled `Basis`, and why you might want it + +The helper is a convenience over one explicit idiom. Compose it yourself when you need +something the arguments do not offer — a bespoke `sigma` range, a third basis, per-basis +priors. From `autogalaxy_workspace:scripts/imaging/features/multi_gaussian_expansion/modeling.py`: + +```python +""" +__Model__ + +Two groups of 30 Gaussians. Within a group every Gaussian shares the centre and ellipticity of +the first, and every `sigma` is *fixed* to a value on a log10 ladder spanning 0.01" to the mask +radius. `ag.lp_basis.Basis` then groups the whole list into one model component the galaxy can +hold (`PyAutoGalaxy:autogalaxy/profiles/basis.py`). +""" +import numpy as np + +total_gaussians = 30 +gaussian_per_basis = 2 + +log10_sigma_list = np.linspace(-2, np.log10(MASK_RADIUS), total_gaussians) + +centre_0 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) +centre_1 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + +bulge_gaussian_list = [] + +for _ in range(gaussian_per_basis): + + gaussian_list = af.Collection( + af.Model(ag.lp_linear.Gaussian) for _ in range(total_gaussians) + ) + + for i, gaussian in enumerate(gaussian_list): + gaussian.centre.centre_0 = centre_0 # one shared y centre + gaussian.centre.centre_1 = centre_1 # one shared x centre + gaussian.ell_comps = gaussian_list[0].ell_comps # shared within this basis + gaussian.sigma = 10 ** log10_sigma_list[i] # fixed, not sampled + + bulge_gaussian_list += gaussian_list + +bulge = af.Model(ag.lp_basis.Basis, profile_list=bulge_gaussian_list) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) +``` + +N = 6. `print(model.info)` prints all sixty Gaussians and looks alarming; almost every +parameter shown is fixed. Trust `model.total_free_parameters`, not the length of the printout. + +Two details that are load-bearing. Assigning `gaussian.sigma = ` **fixes** it — the +whole scheme depends on the widths not being sampled, because a free `sigma` reintroduces the +size-scaling degeneracy the MGE exists to remove. And sharing `ell_comps` *within* a basis but +not *between* bases is precisely what makes the twist fittable. + +The `Basis` constructor also accepts a `regularization` argument that penalises non-smooth +amplitude solutions. Treat it as research-only: the positive-only solver below already fixes +the pathology it was introduced for, and no production analysis uses it — the workspace moved +that branch out of the user-facing script deliberately. + +## Branch — a compact MGE for an unresolved nucleus + +An AGN or nuclear starburst is not the extended stellar light, and forcing one basis to cover +both makes the fit choose between them. Model it as a second, deliberately compact basis: the +same construction with the `sigma` ladder capped at about twice the pixel scale, so the +component cannot be broader than a PSF-convolved point source. + +```python +""" +__Model__ + +The galaxy's light becomes the sum of a diffuse stellar `bulge` MGE and a compact `point` MGE. +The point basis is 10 linear Gaussians sharing one centre and ellipticity with `sigma` values +log-spaced from 0.01" to twice the pixel scale — compact relative to the resolution of the +data, which is what makes it read as a point source rather than a small galaxy +(`PyAutoGalaxy:autogalaxy/analysis/model_util.py`). +""" +point = ag.model_util.mge_point_model_from( + pixel_scales=PIXEL_SCALES, + total_gaussians=10, + centre=(0.0, 0.0), +) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, point=point) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(model.info) +``` + +Adapted from `autogalaxy_workspace:scripts/imaging/features/multi_gaussian_expansion/modeling.py`, +whose `__Point Source__` section builds the same basis line by line first if you want to see +inside the helper. Four extra parameters on top of the bulge MGE. + +There is a genuinely different approach to the same problem — an **operated** light profile, +which is assumed to be already PSF-convolved and so bypasses the convolution entirely. Which +to pick, and why a compact model source is so sensitive to sub-pixel placement, is +[`ag_light_model_extras`](./ag_light_model_extras.md) and +[`../wiki/core/concepts/sky_background_and_operated_profiles.md`](../wiki/core/concepts/sky_background_and_operated_profiles.md). + +## Branch — shapelets, and the price of a signed solver + +Shapelets are an orthonormal basis of Gauss–Hermite (Cartesian) or Gauss–Laguerre (polar) +functions on a single scale `beta`, introduced for exactly this problem +(Refregier 2003, arXiv:astro-ph/0105178). Because the basis is complete, it can represent +morphology an MGE smooths over — the disky, star-forming structure the workspace cites them +for. But completeness has a price, and it is a physical one. + +**Shapelets require negative amplitudes.** The higher-order basis functions oscillate in sign, +and reproducing an arbitrary shape means combining them with signed coefficients. Every other +branch in this skill uses a positive-only linear solver, because a negative surface brightness +is not a thing a galaxy has. Shapelets cannot: you must pass +`ag.Settings(use_positive_only_solver=False)`, and the reconstructed light *will* contain +negative flux — verified, not assumed. That is not a numerical wrinkle to ignore; it is the +model telling you it is a mathematical fit rather than a physical decomposition, and it is why +the workspace recommends trying an MGE alongside. + +```python +""" +__Model__ + +A polar shapelet basis to order n = 5: 11 `ShapeletPolar` components whose `n` and `m` indices +are fixed by construction, sharing one centre, one ellipticity and one scale `beta`. `beta` is +the only shape parameter — it sets the physical size the whole basis is expressed on, so its +prior matters more here than any single amplitude +(`PyAutoGalaxy:autogalaxy/profiles/light/linear/shapelets/polar.py`). +""" +total_n = 5 +total_m = sum(range(2, total_n + 1)) + 1 + +shapelets_bulge_list = af.Collection( + af.Model(ag.lp_linear.ShapeletPolar) for _ in range(total_n + total_m) +) + +n_count = 1 +m_count = -1 + +for shapelet in shapelets_bulge_list: + shapelet.n = n_count + shapelet.m = m_count + + m_count += 2 + + if m_count > n_count: + n_count += 1 + m_count = -n_count + + shapelet.centre = shapelets_bulge_list[0].centre + shapelet.ell_comps = shapelets_bulge_list[0].ell_comps + shapelet.beta = shapelets_bulge_list[0].beta + +bulge = af.Model(ag.lp_basis.Basis, profile_list=shapelets_bulge_list) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(f"Total free parameters = {model.total_free_parameters}") + +""" +__Analysis__ + +The signed solver is requested through `ag.Settings`, which is passed to the analysis rather +than to the model — it is a property of how the inversion is solved, not of what is being +fitted (`PyAutoArray:autoarray/settings.py`). +""" +analysis = ag.AnalysisImaging( + dataset=dataset, + settings=ag.Settings(use_positive_only_solver=False), + use_jax=True, +) +``` + +**N = 5**: centre (2), `ell_comps` (2), `beta` (1). Adapted from +`autogalaxy_workspace:scripts/imaging/features/shapelets/modeling.py` — note that the workspace +loop as written does *not* link `ell_comps` across the basis, which leaves each shapelet's +ellipticity free and inflates the model to N = 43 rather than the N = 3 its prose claims. The +line `shapelet.ell_comps = shapelets_bulge_list[0].ell_comps` above is the fix; N = 3 is what +you get from the spherical variant `ag.lp_linear.ShapeletPolarSph`, which has no `ell_comps` at +all. Check `model.total_free_parameters` rather than trusting a comment — including this one. + +Polar shapelets suit radially organised light, which is most galaxies. +`ag.lp_linear.ShapeletCartesian` (indexed by `n_y`, `n_x` instead of `n`, `m`) and +`ag.lp_linear.ShapeletExponential` are the other two families; the Cartesian basis is not +generally recommended for galaxies. You can see what a basis actually looks like before fitting +anything: + +```python +shapelets = [ + ag.lp_linear.ShapeletCartesian( + n_y=y, n_x=x, centre=(0.0, 0.0), ell_comps=(0.0, 0.0), beta=1.0 + ) + for x in range(5) + for y in range(5) +] + +basis = ag.lp_basis.Basis(profile_list=shapelets) + +aplt.subplot_basis_image( + basis=basis, + grid=ag.Grid2D.uniform(shape_native=(100, 100), pixel_scales=0.05), + output_path="scripts/scratch/shapelet_basis/", + output_format="png", +) +``` + +`subplot_basis_image` renders every component of the basis on a grid, which is the fastest way +to build intuition for what `beta` and the order limit are buying you. It takes `output_path` / +`output_format` but **not** `output_filename` — give each variant its own directory +([`../wiki/core/api/plotting.md`](../wiki/core/api/plotting.md) tabulates that split). Print the +absolute path and offer to open it. + +Shapelets also cost about 0.37 s per likelihood evaluation for ~60 components in the workspace's +measurement — slower than linear profiles, comparable to an MGE, and with more free parameters +than an MGE of similar flexibility. The workspace's own recommendation is to run both and +compare; the concept page +[`../wiki/core/concepts/shapelets.md`](../wiki/core/concepts/shapelets.md) has the full +when-to-use discussion. + +## Choosing between the three + +Fit them in order and let the residuals decide, because each step is cheap given the last: + +1. **Single linear Sersic** — the reference every other fit is judged against, and the source of + the size and concentration you will actually quote. Six parameters. +2. **Linear bulge plus disk** — if the science is a bulge-to-total ratio, this *is* the answer; + a basis will fit better and tell you less. Nine parameters. +3. **MGE, `gaussian_per_basis=2`** — if the residuals show a twist or a radially varying + ellipticity. Six parameters, better fit, less interpretable. +4. **Shapelets** — if the MGE still leaves disky or asymmetric structure and you accept a signed + solution. Five parameters. +5. **Neither is enough** — the residuals sit in discrete off-centre lumps, and no basis with a + single shared centre will absorb them. That is + [`ag_pixelization`](./ag_pixelization.md). + +An MGE and a shapelet basis can be compared by Bayesian evidence only with care, because they +solve different linear problems (positive-only versus signed). Compare them on the residual map +and the chi-squared first, and ask if you want the evidence caveats spelled out. + +## Combine — where this hands off + +- **Pick a search** → [`ag_configure_search`](./ag_configure_search.md). Basis models want + *fewer* live points than a parametric fit, not more — the workspace drops `n_live` to 75 for + both linear profiles and the MGE — because the parameter space is simpler. +- **Run it and read the output folder** → [`ag_run_search`](./ag_run_search.md). +- **Plot the fit** → [`ag_plot_fit`](./ag_plot_fit.md), remembering that a linear profile has to + go through `model_obj_linear_light_profiles_to_light_profiles` before it can be plotted. +- **The basis fits worse than you expected, or the search stalls** → + [`ag_debug_fit_failure`](./ag_debug_fit_failure.md). The two failure modes specific to this + skill are a `sigma` ladder that stops short of the mask radius, and a free `sigma` or `beta` + reintroducing a size degeneracy. +- **Structured residuals a basis cannot reach** → [`ag_pixelization`](./ag_pixelization.md). +- **A neighbour overlapping the galaxy, a residual sky, or a nuclear point source** → + [`ag_light_model_extras`](./ag_light_model_extras.md). +- **A non-parametric shape measurement instead of a model** → + [`ag_ellipse_fitting`](./ag_ellipse_fitting.md), which measures the isophotal twist directly + rather than inferring it from a two-basis MGE. +- **The same galaxy in several bands** → the multi-dataset skill (`ag_multi_dataset`). A basis + across bands is where the factor graph earns its keep, since the shape can be shared while + the amplitudes are solved per band. +- **Chaining a parametric fit into a basis fit** → the search-chaining skill + (`ag_chain_searches`); a Sersic fit's centre and ellipticity make excellent starting priors + for an MGE. + +Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry recording which tier you +used, `total_gaussians` / `gaussian_per_basis` or the shapelet order, the solver you chose, and +what the residuals looked like at each step — the last is the evidence for the model choice and +the first thing a referee will ask for. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Linear profiles](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_5_linear_profiles.ipynb): + builds up from a standard profile to a linear one to a basis, showing what the inversion + actually solves at each step. +- **General reference** — [RTD: Light profiles API](https://pyautogalaxy.readthedocs.io/en/latest/api/light.html): + every standard, linear, operated and basis variant with its module and parameters. +- **Experienced PyAutoGalaxy user** — [workspace: imaging/features/multi_gaussian_expansion/modeling.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/features/multi_gaussian_expansion/modeling.py): + the MGE end to end, including the compact nuclear basis and the run-time and VRAM discussion. + +## Agent procedural checklist + +1. Ask whether the user needs interpretable Sersic parameters or the best achievable fit — that + choice, not the API, picks the branch. +2. Look at the residual map of a single linear Sersic first; structured-but-smooth residuals mean + a basis, discrete lumps mean `ag_pixelization`. +3. Default to `ag.lp_linear` for every parametric component; never use `ag.lp` in a model + without a reason. +4. For an MGE, pass the real `mask_radius`, and reach for `gaussian_per_basis=2` when the + isophotes twist. +5. For shapelets, pass `ag.Settings(use_positive_only_solver=False)` on the *analysis*, link + `ell_comps` across the basis, and tell the user the reconstruction may contain negative flux. +6. Always check `model.total_free_parameters` against what you intended — a basis's `model.info` + is long and mostly fixed, and one unlinked parameter multiplies by the component count. +7. Evaluate one likelihood before committing to a search; run `print_vram_use` before a GPU run + with a large basis. +8. Read intensities out via `max_log_likelihood_galaxies`, + `fit.linear_light_profile_intensity_dict`, or + `model_obj_linear_light_profiles_to_light_profiles` — never expect them in `model.results`. +9. Save the script to `scripts/`, quote every plot's absolute path and offer to open it, then + offer the `wiki/project/` entry. diff --git a/skills/ag_build_interferometer_model.md b/skills/ag_build_interferometer_model.md new file mode 100644 index 0000000..813723c --- /dev/null +++ b/skills/ag_build_interferometer_model.md @@ -0,0 +1,493 @@ +--- +name: ag_build_interferometer_model +description: Model a galaxy observed with a radio or millimetre interferometer (ALMA, JVLA, NOEMA, SMA) by fitting its complex visibilities directly in the uv-plane. Covers loading an `ag.Interferometer` from three FITS files against a real-space mask, choosing the transformer by visibility count and the `raise_error_dft_visibilities_limit` guard that stops a slow exact DFT, why there is no PSF and no over-sampling, composing the model and the `ag.AnalysisInterferometer`, reading dirty images as diagnostics rather than data, the fit and real-space subplots and the filenames they write, the sparse NUFFT operator precomputation for very large visibility tables, and preparing data from a CASA measurement set. Use for visibility-plane data. Not for CCD imaging (`ag_build_imaging_model`), not for a joint imaging-plus-visibilities fit (`ag_multi_dataset`), and not for simulating visibilities (`ag_simulate_dataset`). +--- + +# Modelling a galaxy in the uv-plane + +An interferometer does not photograph the sky. Each pair of antennas measures one complex +number — a **visibility** — which samples the Fourier transform of the sky brightness +distribution at the spatial frequency set by that baseline's projected separation. An +observation is a list of such samples at scattered `(u, v)` coordinates, measured in units of +the observing wavelength, and that list is the data. + +You could invert it to get a picture, and PyAutoGalaxy will make you one. But the `(u, v)` +sampling is sparse and irregular, so the inverse transform gives a **dirty image**: the true sky +convolved with the Fourier transform of the sampling pattern. That beam has extended oscillating +sidelobes, so one compact clump appears as a ring of positive and negative artefacts spread +across the field, and the noise between pixels is strongly correlated in a pattern set by the +array configuration. A per-pixel Gaussian likelihood on such an image is simply the wrong +likelihood, and it biases structural parameters in ways that are hard to notice. Deconvolution +(CLEAN and relatives) makes a plausible picture at the cost of a non-linear, non-invertible +transformation whose effective noise properties are not tractable. + +In the visibility plane the measurement is exactly what the instrument recorded: independent +complex numbers with known RMS uncertainties. A Gaussian likelihood there is the honest one. So +the galaxy model is still built in real space — the same light profiles, the same `af.Model` +composition, the same MGE — and only the *comparison* moves to Fourier space: evaluate the +surface brightness on a real-space grid, transform it to the observed `(u, v)` coordinates, +compare with the measured visibilities. + +The physics and the transformer theory are +[`../wiki/core/concepts/interferometer_theory.md`](../wiki/core/concepts/interferometer_theory.md); +the dataset object and its settings are +[`../wiki/core/api/datasets.md`](../wiki/core/api/datasets.md). The canonical scripts are +`autogalaxy_workspace:scripts/interferometer/start_here.py` and +`autogalaxy_workspace:scripts/interferometer/modeling.py`, and the likelihood is walked term by +term in `autogalaxy_workspace:scripts/interferometer/likelihood_function.py`. + +## Ask + +- *"How many visibilities?"* The single most consequential number. It decides the transformer, + whether a pixelised reconstruction needs the sparse-operator precomputation, and whether this + runs on a laptop. +- *"How big is the emission on the sky, in arcseconds?"* This sets the real-space mask, which is + the main scientific and computational lever here — more so than the imaging mask, because it + sets the *representation* of the model rather than only which data enter the likelihood. +- *"A best-fit model quickly, or a posterior you can quote?"* `af.MultiStartProdigy` for the + former, `af.Nautilus` for the latter; the two workspace scripts differ in exactly this and + nothing else. +- *"Where did the data come from?"* If the answer is "a CASA measurement set", the export step + comes first — see the data-preparation branch. + +## Before you start — the environment + +The default transformer is a JAX-native NUFFT backed by +[`nufftax`](https://github.com/GragasLab/nufftax), which ships with the `[optional]` extras and +requires Python ≥ 3.12. Check it before writing a script, because the workspace's own examples +exit early when it is missing (`autogalaxy_workspace:scripts/interferometer/start_here.py` +`__NUFFT Backend Check__`): + +```python +import importlib.util + +print("nufftax available:", importlib.util.find_spec("nufftax") is not None) +``` + +If it is absent, `pip install nufftax` — or `pip install autogalaxy[optional]`. Repairing the +environment is [`ag_setup_environment`](./ag_setup_environment.md). + +## Branch — the fit + +The deliverable is one script. Adapted from +`autogalaxy_workspace:scripts/interferometer/modeling.py` (the Nautilus configuration, the VRAM +check and the real-space mask) and `autogalaxy_workspace:scripts/interferometer/start_here.py` +(the load and the result). + +```python +""" +Galaxy Structure: Interferometer Modeling +======================================== + +Fit the morphology of a galaxy observed with a millimetre interferometer by modelling its +complex visibilities directly in the uv-plane: define the real-space grid the model is +evaluated on, load the visibilities, noise-map and (u, v) sampling, compose a linear Sersic +bulge and Exponential disk, and sample the posterior with Nautilus so the effective radius and +Sersic index come with errors that do not depend on a deconvolution. + +__Contents__ + +- **Imports:** JAX environment first, then the standard trio. +- **Real Space Mask:** The real-space grid the model is represented on. +- **Dataset:** Load visibilities, noise-map and uv-wavelengths, and pick a transformer. +- **Model:** Compose the galaxy's light profiles. +- **Search:** Configure Nautilus. +- **Analysis:** Bind the model to the visibilities. +- **Model-Fit:** Run the fit and announce the output folder. +- **Result:** Read the best-fit model and its dirty-image diagnostics. +""" +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 + +""" +__Real Space Mask__ + +The real-space mask is the interferometer analogue of the imaging mask, and it does more work. +`shape_native` and `pixel_scales` set the resolution at which the galaxy is represented *before* +it is transformed, and the radius sets the field being modelled. Too small and you truncate real +emission and alias it back into the fit; too large and every likelihood evaluation transforms +far more pixels than the data can constrain. Match the pixel scale to a few pixels per +synthesised beam and the radius to the emission you can see in the dirty image +(`PyAutoArray:autoarray/mask/mask_2d.py`). +""" +MASK_RADIUS = 3.5 + +real_space_mask = ag.Mask2D.circular( + shape_native=(256, 256), + pixel_scales=0.1, + radius=MASK_RADIUS, +) + +""" +__Dataset__ + +Three FITS files, not imaging's three: the complex visibilities (shape `(n_vis,)`), the +per-visibility complex RMS, and the `(u, v)` sampling with shape `(n_vis, 2)`. There is **no +PSF** — the sampling function plays that role and it is already in the data — and consequently +nothing to deconvolve at load time. + +`transformer_class` chooses how the model image is mapped to the observed `(u, v)` coordinates. +`ag.TransformerNUFFT` is the default and the right answer at essentially any dataset size +(`PyAutoArray:autoarray/dataset/interferometer/dataset.py`, +`PyAutoArray:autoarray/operators/transformer.py`). +""" +DATASET_PATH = Path("dataset") / "interferometer" / "" + +dataset = ag.Interferometer.from_fits( + data_path=DATASET_PATH / "data.fits", + noise_map_path=DATASET_PATH / "noise_map.fits", + uv_wavelengths_path=DATASET_PATH / "uv_wavelengths.fits", + real_space_mask=real_space_mask, + transformer_class=ag.TransformerNUFFT, +) + +print(f"n_vis = {dataset.data.shape[0]}") + +aplt.subplot_interferometer_dirty_images( + dataset=dataset, + output_path="scripts/scratch//", + output_filename="dirty_images", + output_format="png", +) + +""" +__Model__ + +Identical to an imaging model — that is the point. A linear `Sersic` bulge and a linear +`Exponential` disk with their centres paired, so the two components describe one galaxy rather +than drifting apart. `lp_linear` solves each profile's `intensity` by linear inversion instead +of sampling it, which removes one dimension per component +(`PyAutoGalaxy:autogalaxy/profiles/light/linear/sersic.py`). + +There is **no over-sampling section**, and its absence is deliberate: an interferometer does not +observe in a way that makes sub-pixel integration of the light profile meaningful, so every +interferometer calculation runs without it +(`autogalaxy_workspace:scripts/interferometer/modeling.py` `__Over Sampling__`). +""" +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) + +""" +__Search__ + +`Nautilus` returns the full posterior — errors and covariances, the thing you can quote +(`PyAutoFit:autofit/non_linear/search/nest/nautilus/search.py`). The folder's `start_here.py` +instead uses `af.MultiStartProdigy`, a multi-start gradient optimizer which is far faster and +returns a single best-fit model with no errors at all; use that to check the model and the data +make sense, then this when the numbers go in a paper. +""" +search = af.Nautilus( + path_prefix=Path("interferometer"), + name="sersic_exp", + unique_tag="", + n_live=100, + n_batch=50, + iterations_per_quick_update=10000, + live_visual_update=False, +) + +""" +__Analysis__ + +`ag.AnalysisInterferometer` defines the `log_likelihood_function` the search calls: evaluate the +galaxy image on the real-space grid, transform it to the `(u, v)` points, and compare against the +measured visibilities with a Gaussian likelihood over the real and imaginary parts +(`PyAutoGalaxy:autogalaxy/interferometer/model/analysis.py`). +""" +analysis = ag.AnalysisInterferometer(dataset=dataset, use_jax=True) + +analysis.print_vram_use(model=model, batch_size=search.batch_size) + +""" +__Model-Fit__ + +Results are written to the output folder on the fly, from the best model found so far. +""" +print(f"Output folder: {search.paths.output_path.resolve()}") + +result = search.fit(model=model, analysis=analysis) + +""" +__Result__ + +`fit.png` in the output folder is the visibility-space fit; `dirty_images` are the real-space +diagnostics. Both are written automatically, and both are reproducible here at whatever scaling +the science needs. +""" +print(result.info) + +PLOT_DIR = Path("scripts") / "scratch" / "" + +aplt.subplot_fit_interferometer( + fit=result.max_log_likelihood_fit, + output_path=str(PLOT_DIR / "fit"), + output_format="png", +) + +aplt.subplot_fit_dirty_images( + fit=result.max_log_likelihood_fit, + output_path=str(PLOT_DIR / "fit_dirty"), + output_format="png", +) + +aplt.subplot_fit_real_space( + fit=result.max_log_likelihood_fit, + output_path=str(PLOT_DIR / "fit_real_space"), + output_format="png", +) + +print(f"Saved to: {PLOT_DIR.resolve()}") +``` + +The output-folder anatomy — `files/`, `image/` (which here holds `dirty_images.fits` alongside +`fit.fits`), `model.info`, `model.results`, `search.summary` and the `` resume +behaviour — is `__Output Folder Layout__` in +[`autogalaxy_workspace/scripts/interferometer/modeling.py`](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/interferometer/modeling.py). +Link it rather than copying the tree. Driving the fit, watching it, and the resume trap are +[`ag_run_search`](./ag_run_search.md). + +## Branch — choosing the transformer + +Three transformers exist (`PyAutoArray:autoarray/operators/transformer.py`), and the choice is +not cosmetic. + +- **`ag.TransformerNUFFT`** — the default, and the default in `ag.Interferometer.from_fits` + itself. A non-uniform FFT backed by `nufftax`, a pure-JAX implementation that jit-compiles and + vmap-batches with the rest of the library. Recommended at any dataset size. Because it is + JAX-native, light-profile fits now run at full GPU speed for visibility counts in the tens or + hundreds of millions — the scale of a high-resolution ALMA observation. +- **`ag.TransformerDFT`** — the exact discrete Fourier transform. Slower once `n_vis` is large, + because it is genuinely `O(n_vis × n_pixels)`, but it is the reference implementation. It is + worth running once on a small dataset to check a NUFFT result: on a 190-visibility SMA + sampling with a 64 × 64 real-space grid the two agree in log likelihood to roughly one part in + 10⁹, so a disagreement bigger than that means something about the NUFFT setup is wrong rather + than the transform being approximate. It is also what the pixelised reconstruction's + sparse-operator path uses. +- **`ag.TransformerNUFFTPyNUFFT`** — a legacy `pynufft`-backed transformer, kept as a non-JAX + fallback. It is not JAX-traceable, so it forfeits GPU acceleration *and* the gradient-based + searches, which need the likelihood's derivatives. + +**There is a guard, and it is there for a reason.** `ag.Interferometer` (and `from_fits`) raises +`raise_error_dft_visibilities_limit=True` by default, and it will refuse to build a dataset with +more than **10,000 visibilities** while `transformer_class=ag.TransformerDFT`. The DFT at that +scale is not slow-but-tolerable, it is unusable. The escape hatch exists — +`raise_error_dft_visibilities_limit=False` — and is legitimate when you are deliberately +profiling the exact path, but the honest reading of that error is "you meant to use the NUFFT". + +The reason dataset size no longer forces a modelling decision is worth stating plainly to the +user: it used to be that a large visibility table pushed you toward a pixelised reconstruction +for performance reasons. It does not any more +(`autogalaxy_workspace:scripts/interferometer/start_here.py` `__Number of Visibilities__`). +Choose a pixelisation because the *morphology* demands it, never because the file is big. + +## Branch — dirty images are diagnostics, not data + +The thing you look at is not the thing you fit, and keeping that straight is most of the skill +of interpreting an interferometer fit. + +```python +aplt.subplot_interferometer_dataset(dataset=dataset, output_path=..., output_format="png") +aplt.subplot_interferometer_dirty_images(dataset=dataset, output_path=..., output_format="png") + +aplt.plot_array(array=dataset.dirty_image, title="Dirty Image", output_path=..., + output_filename="dirty_image", output_format="png") +``` + +`autogalaxy_workspace:scripts/interferometer/plot.py`. `subplot_interferometer_dataset` shows the +visibility-space quantities and the uv coverage; the dirty-image subplot transforms *back* for +display only. On the fit side there are three: + +| Call | Shows | Writes into `output_path` | +|---|---|---| +| `aplt.subplot_fit_interferometer` | the visibility-space fit — data, model visibilities, residuals, chi-squared | `fit.png` | +| `aplt.subplot_fit_dirty_images` | the same fit transformed back to real space | `fit_dirty_images.png` | +| `aplt.subplot_fit_real_space` | the model's real-space image and the reconstruction, with no transform applied | `fit_real_space.png` | + +All three write a **fixed filename stem** and do not accept `output_filename`, so a separate +directory per figure is how you keep two fits apart — the same rule as the imaging subplots +([`../wiki/core/api/plotting.md`](../wiki/core/api/plotting.md)). `plot_array`, +`subplot_interferometer_dataset` and `subplot_interferometer_dirty_images` *do* accept +`output_filename`; passing it to one of the fit subplots raises `TypeError`. + +Reading the residuals needs the dirty beam in mind. The individual real-space maps are available +on the fit object — `fit.dirty_image`, `fit.dirty_model_image`, `fit.dirty_residual_map`, +`fit.dirty_normalized_residual_map`, `fit.dirty_chi_squared_map` +(`PyAutoGalaxy:autogalaxy/interferometer/fit_interferometer.py`) — and a ringing pattern of +positive and negative residual around a real source is the *expected* signature of an imperfectly +subtracted beam, not a model failure. What is a model failure is a residual with structure that +survives across baselines, or a normalised residual map with coherent large-scale sign. Judge the +fit in visibility space first (`subplot_fit_interferometer`) and use the dirty maps to see *where* +on the sky the misfit lives. The general residual-inspection discipline is +[`ag_plot_fit`](./ag_plot_fit.md). + +## Branch — a pixelised reconstruction, and very large visibility tables + +A pixelisation remains the right choice when the morphology is genuinely irregular — clumpy star +formation, strong asymmetry — and no smooth profile or basis captures it. The linear algebra is +the same as for imaging with the transformer standing in for the PSF convolution +([`../wiki/core/concepts/inversions_and_pixelizations.md`](../wiki/core/concepts/inversions_and_pixelizations.md)), +and the composition procedure is `ag_pixelization`. +`autogalaxy_workspace:scripts/interferometer/features/pixelization/modeling.py` is the +interferometer variant. + +What is specific to visibilities is the **sparse NUFFT operator**. Pixelised reconstruction needs +dense linear algebra that would be prohibitive over a large visibility table, so the dataset can +precompute an operator matrix that exploits the sparsity of the reconstruction +(`PyAutoArray:autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py`). +Building it costs anything from milliseconds to hours depending on visibility count, real-space +mask size, and CPU versus GPU — and without it that cost is paid before *every* fit. So compute +it once and cache it, which is exactly what +`autogalaxy_workspace:scripts/interferometer/features/pixelization/many_visibilities_preparation.py` +does: + +```python +import numpy as np + +dataset = dataset.apply_sparse_operator( + use_jax=True, + chunk_k=2048, # visibilities processed at a time; raise until memory complains + show_progress=True, # a progress bar, which matters when this takes an hour + show_memory=True, +) + +nufft_precision_operator = dataset.psf_precision_operator_from( + use_jax=True, chunk_k=2048, show_progress=True, show_memory=True +) + +np.save( + file=DATASET_PATH / f"nufft_precision_operator_{MASK_RADIUS}.npy", + arr=nufft_precision_operator, + allow_pickle=False, +) +``` + +Reload it with `np.load(..., allow_pickle=False)` and hand it back through +`apply_sparse_operator(nufft_precision_operator=...)`. Two practical notes: the cached file is +tied to the real-space mask it was built with, which is why the mask radius belongs in the +filename; and on a modern GPU a million-visibility operator takes under a minute, so on that +hardware computing it inline is fine and the caching is a laptop-and-CPU concern. + +That preparation script also carries a commented-out profiling block that fabricates a synthetic +visibility table at a chosen size, which is the honest way to find out how long *your* dataset +will take before committing to it. + +## Branch — preparing your own data + +The three FITS files and what must be true of them — units, the `(n_vis, 2)` shape of +`uv_wavelengths` in wavelengths rather than metres, the pixel-scale choice for the real-space +mask, and the optional light-centre and extra-galaxy centre files — are +`autogalaxy_workspace:scripts/interferometer/data_preparation.py`. Read it before loading real +data; the pixel-scale and uv-unit mistakes are both silent. + +Coming from a CASA measurement set, the export path is +`autogalaxy_workspace:scripts/interferometer/casa_reduction.py` — and that script **says of +itself that it is incomplete and in development**, so treat it as a set of hints rather than a +recipe. What it does establish reliably is the shape contract: an ALMA `.ms` stores visibilities +as `(2, n_spw, n_c, n_v, 2)` with the leading axis the two polarisations; `split` peels off +spectral windows and averages channels; and the `uv_wavelengths` you export must be in +wavelengths, not metres. Do not invent the rest of that workflow — if a user needs a step the +script does not cover, say the script is incomplete and point them at the project's support +channel rather than guessing at CASA semantics. + +Extra galaxies whose emission blends with the target are handled the same way as for imaging, and +`autogalaxy_workspace:scripts/interferometer/features/extra_galaxies/modeling.py` is the +visibility version; the levers themselves are `ag_light_model_extras`. + +## Branch — cost, and proving the script first + +Two numbers dominate run time: the visibility count, which sets the transform cost per +evaluation, and the real-space pixel count, which sets it too. The levers, in order of effect: + +- **Shrink the real-space mask** to the informative region. This is the biggest one and the one + most often left at a default. +- **Average channels or baselines** where the science allows. +- **Keep the model parametric** unless the morphology genuinely demands a pixelisation. +- **Choose the transformer deliberately** rather than inheriting a default you have not checked. + +On a GPU, check VRAM before a long run — `analysis.print_vram_use(model=model, +batch_size=search.batch_size)`, twenty or thirty seconds +(`autogalaxy_workspace:scripts/interferometer/modeling.py` `__VRAM Use__`). A parametric fit to a +small visibility set is a few hundred megabytes; a pixelised reconstruction on a large one can +run past ten gigabytes. + +And smoke-test before the real run: + +```bash +PYAUTO_TEST_MODE=2 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/interferometer_fit.py +``` + +Level 2 calls the likelihood exactly once, which is the fastest proof that the mask, the +transformer and the model all agree with each other. Level 1 if you want the fit products +written. The levels and the writable-cache variables are +[`../wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md). If a NumPy stack +trace would be easier to read than a JAX one, `use_jax=False` or `PYAUTO_DISABLE_JAX=1` — see +[`ag_debug_fit_failure`](./ag_debug_fit_failure.md). + +## Combine + +- [`ag_simulate_dataset`](./ag_simulate_dataset.md) — simulate visibilities with known truth + from a real array's `uv_wavelengths` and fit those first. This is the cleanest way to separate + a model problem from a data problem, and note that `ag.SimulatorInterferometer` defaults to + `TransformerDFT` where the dataset loader defaults to the NUFFT, so set it explicitly when you + simulate at scale. +- [`ag_multi_dataset`](./ag_multi_dataset.md) — fit visibilities jointly with CCD imaging of the + same galaxy, which is a factor graph with `ag.AnalysisInterferometer` on one factor. +- [`ag_configure_search`](./ag_configure_search.md) — `af.MultiStartProdigy` for a fast MAP + check, `af.Nautilus` for the posterior, and the `unique_tag` resume semantics. +- [`ag_run_search`](./ag_run_search.md) — driving the fit, the output folder, and the stale-resume + trap. +- [`ag_plot_fit`](./ag_plot_fit.md) — figure scaling, log10 stretch and fixed colour limits, all + of which apply to the dirty maps. +- [`ag_load_results`](./ag_load_results.md) — getting the posterior and derived quantities back. +- `ag_pixelization` and `ag_basis_profiles` — the components you would reach for when a Sersic + cannot describe the morphology. + +When the fit is worth keeping, offer (default-yes) to record it in a dated +`wiki/project/YYYY-MM-DD-.md` entry: the array and visibility count, the real-space mask +and transformer chosen and why, what was inferred, and the output path — per +[`_style.md`](./_style.md) property #5. The mask and transformer choices are the two decisions a +reader of the fit will want justified. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Fitting](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_1_introduction/tutorial_3_fitting.ipynb): + how a model image becomes a log likelihood via residuals, noise and chi-squared. The lecture + series has no interferometer chapter, so this teaches the likelihood in its imaging form — + everything transfers except that the comparison happens after a Fourier transform. +- **General reference** — [RTD: Features](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_3_features.html): + the interferometry section of the feature tour, with pointers onward. +- **Experienced PyAutoGalaxy user** — [workspace: interferometer/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/interferometer/start_here.py): + the minimal end-to-end uv-plane fit, with the transformer and NUFFT-backend sections this + skill mirrors. + +## Agent procedural checklist + +1. Confirm `nufftax` is importable before writing a script that defaults to the NUFFT. +2. Ask for the visibility count and the on-sky extent of the emission — they set the transformer + and the real-space mask. +3. Set `shape_native`, `pixel_scales` and `radius` of the real-space mask deliberately, and say + what each one buys; never leave them at a copied default on real data. +4. Load with `ag.Interferometer.from_fits` and an explicit `transformer_class`; if the user asks + for `TransformerDFT` above 10,000 visibilities, explain the guard rather than disabling it. +5. Compose the model exactly as for imaging, and do **not** add an over-sampling section. +6. Plot the dirty images before fitting, quote the absolute path, and offer to open it. +7. Validate with `PYAUTO_TEST_MODE=2`; on a GPU run `analysis.print_vram_use` before a long fit. +8. Announce the output path at launch and name `model.results` and `image/fit.png` first. +9. Judge the fit in visibility space first, then use the dirty maps to locate the misfit; say + plainly that beam sidelobes in a residual map are expected. +10. For a pixelised fit on a large table, cache the sparse operator and name the file after the + mask radius it was built with. +11. Offer the `wiki/project/` entry recording the mask and transformer choices. diff --git a/skills/ag_chain_searches.md b/skills/ag_chain_searches.md new file mode 100644 index 0000000..932248a --- /dev/null +++ b/skills/ag_chain_searches.md @@ -0,0 +1,420 @@ +--- +name: ag_chain_searches +description: Break one hard fit into a sequence of easier non-linear searches, using each result to initialize the next. Covers when chaining beats a single search and when it does not, passing a whole profile or galaxy with `result.model` (which preserves the original priors) versus narrowing them with `result.model_centred` and its absolute, relative and bounded variants, the `TruncatedGaussianPrior` those produce and where its sigma comes from in the `priors` config `width_modifier`, fixing a component with `result.instance` to remove its dimensions, and the output-path convention that keeps a chain's searches together. Use when a model is too complex for one search, when an early search should use cheap settings and a later one accurate ones, or when several galaxies must be added one at a time. Not for choosing a single search (`ag_configure_search`), not for driving one fit (`ag_run_search`), and not for combining several datasets (`ag_multi_dataset`). +--- + +# Chaining searches + +Every other modelling skill in this workspace composes one model and hands it to one search. +Chaining is the alternative shape: split the problem into a sequence of searches, and use what +each one learned to tell the next where to look. + +The reason is statistical, not cosmetic. A non-linear search has to map a parameter space it knows +nothing about, and the cost of doing that grows sharply with dimensionality while the risk of +settling in a local maximum grows with multi-modality. A galaxy model with a bulge, a disk, a sky +level and two blended neighbours has both problems at once. Fit the bulge alone and the space is +small, unimodal and mapped in minutes; then the disk only has to be added to a galaxy whose size +and centre are already known. Each search in a chain solves a problem that is easy *because* of +the one before it. + +Three distinct benefits, and they are worth separating because they justify chaining in different +situations: + +- **Lower dimensionality per search.** Earlier searches fit simpler models with a parameter space + that can be sampled properly, which reduces the chance of an incorrect local maximum in the + final one. +- **Cheap settings early, accurate settings late.** An early search only needs to find the high- + likelihood region, so it can use a fast optimizer or a small `n_live`. Only the final search + needs settings good enough to quantify errors. +- **Cheap approximations early.** An early search can use a smaller mask, coarser over-sampling + or a simplified component, and revert to the accurate-but-expensive settings at the end. The + approximation costs some fidelity in the early result, which is fine when that result is only + being used to place a prior. + +This skill has **one** grounding script: `autogalaxy_workspace:scripts/guides/modeling/chaining.py`. +It is deliberately an API overview rather than a library of pipelines, and this skill is scoped to +match it — the mechanics of prior passing, done honestly, rather than a catalogue of recipes that +do not exist upstream. Chapter 3 of the HowToGalaxy lecture series is where the pedagogy lives. +Read [`../wiki/core/concepts/non_linear_search.md`](../wiki/core/concepts/non_linear_search.md) for +how run time scales with the model, which is the quantity chaining exists to manage. + +## Ask + +- *"What makes this fit hard — too many parameters, or a search that keeps landing somewhere + wrong?"* Both are chaining's territory, but they suggest different splits: dimensionality + suggests adding components one at a time, multi-modality suggests fixing the thing the search + keeps getting wrong. +- *"What would you fit first if you had to pick one component?"* Usually the brightest and + simplest. If the user cannot name one, chaining may not be the right tool — see the branch on + when not to chain. +- *"Is a single search actually failing, or just slow?"* A slow-but-correct search is often better + served by a faster sampler or JAX than by a chain, because a chain adds a whole class of + failure: a wrong early result confidently narrowing a later prior. + +## Branch — the two-search chain + +The deliverable is one script containing the whole chain, so it can be re-run end to end. Adapted +from `autogalaxy_workspace:scripts/guides/modeling/chaining.py`. + +```python +""" +Galaxy Structure: Chained Searches +================================= + +Fit a galaxy's bulge alone, then add a disk with the bulge's parameters initialized from the +first fit. Splitting the fit this way lets each search sample a parameter space it can actually +map, and makes the two-component decomposition a refinement of a known galaxy rather than a +fourteen-dimensional search from broad priors. + +__Contents__ + +- **Imports:** JAX environment first, then the standard trio. +- **Dataset:** Load, mask and over-sample the imaging once, for every search. +- **Paths:** One output prefix shared by every search in the chain. +- **Model (Search 1):** A single linear Sersic bulge. +- **Model-Fit (Search 1):** Fit it, cheaply. +- **Model Chaining:** Pass the bulge into search 2 and add a disk. +- **Model-Fit (Search 2):** Fit the two-component model with accurate settings. +- **Result:** Read the final 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__ + +Loaded, masked and over-sampled once. Every search in the chain shares one `analysis`, which is +what makes their likelihoods comparable — if search 2 used a different mask, its log likelihood +would not be on the same scale as search 1's, and the priors passed between them would refer to a +different dataset. +""" +DATASET_PATH = Path("dataset") / "imaging" / "" + +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=0.1, +) + +mask = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=dataset.pixel_scales, + radius=3.0, +) + +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) + +analysis = ag.AnalysisImaging(dataset=dataset, use_jax=True) + +""" +__Paths__ + +One `path_prefix` for the whole chain, with each search distinguished only by `name`. This keeps +every stage of one chain in one directory tree, which matters when you come back to a fit and +need to see what search 1 actually inferred before trusting search 2. +""" +path_prefix = Path("imaging") / "chaining" / "" + +""" +__Model (Search 1)__ + +The bulge alone. A linear `Sersic` solves its own `intensity`, so this is a compact space the +search can map quickly (`PyAutoGalaxy:autogalaxy/profiles/light/linear/sersic.py`). +""" +bulge = af.Model(ag.lp_linear.Sersic) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model_1 = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(model_1.info) + +""" +__Model-Fit (Search 1)__ + +Cheap settings on purpose: this search only has to find the high-likelihood region, not measure +errors. Inspect its result before relying on it — a chain built on a bad first search propagates +that error into every subsequent prior, with more confidence at every step. +""" +search_1 = af.Nautilus( + path_prefix=path_prefix, + name="search[1]__bulge", + unique_tag="", + n_live=100, +) + +result_1 = search_1.fit(model=model_1, analysis=analysis) + +print(result_1.info) + +""" +__Model Chaining__ + +`result_1.model_centred` returns the search-1 model with every free parameter replaced by a +`TruncatedGaussianPrior` centred on that parameter's median from search 1 +(`PyAutoFit:autofit/non_linear/result.py`, +`PyAutoFit:autofit/mapper/prior/truncated_gaussian.py`). Search 2 therefore starts sampling in +the region search 1 found, rather than from broad priors. + +Note the attribute: `result_1.model` passes the component through with its **original** priors +unchanged, which is useful when you want the structure of a previous model but not its +constraints. `model_centred` is the one that narrows. +""" +bulge = result_1.model_centred.galaxies.galaxy.bulge +disk = af.Model(ag.lp_linear.Exponential) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, disk=disk) + +model_2 = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(model_2.info) + +""" +__Model-Fit (Search 2)__ + +The accurate search. Read search 2's `model.info` in the output folder to confirm the priors +arrived as intended before trusting the result. +""" +search_2 = af.Nautilus( + path_prefix=path_prefix, + name="search[2]__bulge_disk", + unique_tag="", + n_live=150, +) + +result_2 = search_2.fit(model=model_2, analysis=analysis) + +""" +__Result__ +""" +print(result_2.info) + +PLOT_DIR = Path("scripts") / "scratch" / "" / "chained" + +aplt.subplot_fit_imaging( + fit=result_2.max_log_likelihood_fit, + output_path=str(PLOT_DIR), + output_format="png", +) + +print(f"Saved to: {PLOT_DIR.resolve()}") +``` + +Whole components pass, not just single parameters: `result_1.model_centred.galaxies.galaxy.bulge` +carries every one of the bulge's parameters across in one expression, provided the component's +*type* has not changed. A single parameter passes the same way — +`galaxy.bulge.effective_radius = result_1.model_centred.galaxies.galaxy.bulge.effective_radius` — +which is how you narrow some parameters and leave others broad. + +## Branch — `model`, `model_centred`, and `instance` + +This is the part to get right, and the distinction is sharper than the workspace script's prose +suggests. Three attributes on a `Result`, three different things: + +**`result.model`** — the same model with the **same priors it had before**. On the released stack +this attribute preserves the priors rather than narrowing them +(`PyAutoFit:autofit/non_linear/result.py`), so a component passed this way brings its structure +across but no information from the fit. That is occasionally exactly what you want — reuse the +composition, refit from scratch — but it is *not* prior passing, and reading the workspace +script's older description of it as producing narrowed Gaussians will mislead you. Check what you +actually got: + +```python +print(type(model_2.galaxies.galaxy.bulge.sersic_index).__name__) +``` + +**`result.model_centred`** — the narrowing one. Every free parameter becomes a +`TruncatedGaussianPrior` whose `mean` is that parameter's median from the previous search and whose +`sigma` comes from the parameter's `width_modifier` in the `priors` config. Three further variants +let you override that width: + +| Attribute | Prior produced | `sigma` | +|---|---|---| +| `result.model_centred` | `TruncatedGaussianPrior` | from the parameter's config `width_modifier` | +| `result.model_centred_absolute(a=0.1)` | `TruncatedGaussianPrior` | exactly `a`, for every parameter | +| `result.model_centred_relative(r=0.2)` | `TruncatedGaussianPrior` | `r × mean`, per parameter | +| `result.model_centred_max_lh_bounded(b=0.3)` | `UniformPrior` | bounds at `mean ± b` | + +**`result.instance`** — the maximum-likelihood *values*, as fixed numbers. A component passed this +way has no free parameters at all, so search 2's dimensionality drops by that component's +parameter count. This is the lever for "fit the bulge, freeze it, then fit the disk against a +fixed bulge", and it is much stronger medicine than a narrow prior: the later search cannot +revisit the frozen component even if the added component would have changed it. + +```python +galaxy = af.Model( + ag.Galaxy, + redshift=0.5, + bulge=result_1.instance.galaxies.galaxy.bulge, # fixed, contributes 0 parameters + disk=af.Model(ag.lp_linear.Exponential), +) +``` + +The common three-stage pattern uses all of it: fit the bulge; fix it and fit the disk; then free +both from the second result's narrowed priors for a final search where everything moves together. +Freezing permanently is a real bias risk, because a bulge fitted without a disk absorbs some of +the disk's light — so the third stage is not optional if the decomposition is the result you are +quoting. + +## Branch — where `sigma` comes from, and its two traps + +`model_centred` reads each parameter's width from the `width_modifier` field of its entry in the +`priors` config (`autogalaxy_assistant:config/priors/`, mapped in +[`../wiki/core/api/configuration.md`](../wiki/core/api/configuration.md)). Two forms: + +- **Absolute** — `sigma` is the config value directly. For a linear `Sersic`, `sersic_index` + carries `Absolute: 1.5`, each `centre` component `Absolute: 0.05`, and each `ell_comps` + component `Absolute: 0.2`. +- **Relative** — `sigma` is that fraction of the inferred value. `effective_radius` carries + `Relative: 1.0`. + +The reasoning behind which is which is worth understanding, because it tells you when to override. +A relative width on a centre makes no sense: a galaxy centred near `(0.0, 0.0)` would get a sigma +near zero and the next search could not move it at all. An absolute width on an `intensity` makes +no sense either, because intensity depends on units, exposure and brightness, so no single number +generalises. Hence absolute centres and relative radii. + +Two traps follow directly, and both are worth stating to a user before they hit them: + +**Relative widths break on parameters that can be negative or zero.** +`model_centred_relative(r=...)` applies one relative width to *every* parameter, including +`ell_comps`, whose median is legitimately near zero and often negative. That yields a negative +`sigma` and raises rather than silently misbehaving — which is the good outcome, but it does mean +the blanket relative form is unusable on a general light-profile model. Use `model_centred` +(per-parameter config widths), or `model_centred_absolute` where a single absolute width is +defensible for the parameters you are passing. + +**The bounded variant does not inherit the parameter's physical limits.** +`model_centred_max_lh_bounded(b=...)` builds a `UniformPrior` at `mean ± b` regardless of the +parameter's configured limits, so a `sersic_index` inferred near its ceiling of 5.0 can come back +with an upper limit above 5.0. Set the bound with the parameter's range in mind, or clamp the +prior yourself afterwards. + +The width you want is a genuine trade-off with no default that is right everywhere: wide enough +that a better nearby solution can still be found, narrow enough that the search is not doing +search-1's work over again. The shipped `width_modifier` values are a considered balance rather +than a physical truth, and overriding them for one parameter you understand well is normal. + +## Branch — when chaining is the wrong tool + +Chaining adds a failure mode that a single search does not have: **a wrong early result narrows a +later prior around the wrong place, and the later search never looks anywhere else.** The final +posterior then looks tight and well-behaved and is simply wrong. Three situations where a single +search is better: + +- **The single search works and is merely slow.** Reach for a faster sampler, a linear light + profile, JAX or a GPU first ([`ag_configure_search`](./ag_configure_search.md), + [`ag_run_search`](./ag_run_search.md)). Complexity you do not need is complexity that can be + wrong. +- **The early model is not a subset of the late one.** Prior passing assumes the parameter means + something similar in both searches. Passing a `Sersic`'s parameters into an `Exponential`, or a + parametric component's into a pixelised one, is either an error or meaningless. +- **The parameter you fixed early is the one you are measuring.** If a bulge fitted without a disk + is the number going in the paper, do not chain to a frozen bulge — chain to a narrowed one and + let the final search move it. + +Always inspect the intermediate result rather than trusting the chain. `result_1.info`, and +`model.info` inside search 2's output folder, are the two things to read before believing the +final number. + +## Branch — the practical mechanics + +**One analysis, reused.** Build the `ag.AnalysisImaging` once, before search 1, and pass the same +object to every search — unless you are deliberately changing the data between stages (a smaller +mask early, the full one late), in which case say so explicitly, because the log likelihoods are +then not comparable across stages. + +**Naming.** One `path_prefix` for the chain, and `name="search[1]__"` per stage. The bracket +convention is the workspace's and it sorts correctly; the `__` suffix is what makes an +output tree readable a month later. + +**`unique_tag`.** Keep it identical across the chain's searches, since it identifies the *dataset*. +Each stage still lands in its own folder because `name` differs. And the identifier trap from +[`ag_configure_search`](./ag_configure_search.md) applies to every stage: it does not hash your +data's pixel values, so if you swap datasets without changing `unique_tag`, every search in the +chain silently returns the previous galaxy's result. + +**Smoke test the whole chain, not one search.** A chain's characteristic bug is a mis-typed +attribute path in the passing expression, and that only fires at the second search: + +```bash +PYAUTO_TEST_MODE=1 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/chained_fit.py +``` + +Level 1 rather than level 2 here: level 1 runs a real (tiny) search at each stage, so a `Result` +exists with samples for the next stage to pass from. Level 2 bypasses the sampler and leaves the +chain nothing meaningful to narrow around. The levels are +[`../wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md), and the priors a +test-mode chain produces are wiring evidence only — never quote them. + +## Combine + +- [`ag_configure_search`](./ag_configure_search.md) — cheap settings for early searches and + accurate ones for the last; also the `unique_tag` semantics every stage inherits. +- [`ag_run_search`](./ag_run_search.md) — driving each stage and reading its output folder. +- [`ag_build_imaging_model`](./ag_build_imaging_model.md) — composing the models the chain steps + between, including prior customisation by hand when you want to override a passed width. +- [`ag_multi_galaxy_and_cluster`](./ag_multi_galaxy_and_cluster.md) — adding galaxies one at a + time is the most common real use of chaining. +- [`ag_multi_dataset`](./ag_multi_dataset.md) — its one-by-one branch is a chain across datasets + rather than across model complexity. +- [`ag_debug_fit_failure`](./ag_debug_fit_failure.md) — when a single search lands somewhere + unphysical, which is the diagnosis that most often leads here. +- `ag_pixelization` — a smooth-profile fit chained into a pixelised one is a natural pairing, + because the pixelisation wants a decent light model to start from. + +When a chain is worth keeping, offer (default-yes) to record it in a dated +`wiki/project/YYYY-MM-DD-.md` entry: why the fit was split, what each stage fitted, what was +passed as a narrowed prior versus fixed as an instance, and the shared output prefix — per +[`_style.md`](./_style.md) property #5. The split and the passing decisions are the reasoning a +reader will want, and they are invisible in the final `model.results`. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Search chaining](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_3_search_chaining/tutorial_1_search_chaining.ipynb): + the first tutorial of the chaining chapter, which builds the idea from a fit that fails without + it; `tutorial_2_prior_passing` then walks the prior mechanics and `tutorial_3_x2_galaxies` + applies the chain to two galaxies. +- **General reference** — [RTD: Configs](https://pyautogalaxy.readthedocs.io/en/latest/general/configs.html): + how the config files are laid out and looked up — the `priors` tree is where each parameter's + `width_modifier` lives, and hence where a passed prior's `sigma` comes from. +- **Experienced PyAutoGalaxy user** — [workspace: guides/modeling/chaining.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/guides/modeling/chaining.py): + the API overview this skill is grounded in, with its own long discussion of absolute versus + relative widths. + +## Agent procedural checklist + +1. Ask whether the single search is failing or merely slow; recommend a faster search before a + chain when it is only slow. +2. Choose the split with the user and say what each stage buys — dimensionality, settings, or a + cheap approximation. +3. Build the dataset and one `ag.AnalysisImaging` before search 1 and reuse it, or state + explicitly why a stage uses different data. +4. Write the whole chain into one script under one `path_prefix`, with `search[N]__` names. +5. Use `result.model_centred` (or a `model_centred_*` variant) to narrow, `result.instance` to + fix, and `result.model` only when you deliberately want the original priors back. +6. Verify what you built: print the passed prior's type, and read `model.info` in the later + search's output folder. +7. Do not use `model_centred_relative` on a model containing `ell_comps` or a centre — it raises + on a negative median. +8. Inspect `result_1.info` before trusting anything downstream of it. +9. If a component is frozen mid-chain, plan the final stage that frees it again before quoting its + parameters. +10. Validate the whole chain with `PYAUTO_TEST_MODE=1` (not 2 — the chain needs samples to pass). +11. Offer the `wiki/project/` entry recording the split and the passing decisions. diff --git a/skills/ag_ellipse_fitting.md b/skills/ag_ellipse_fitting.md new file mode 100644 index 0000000..6d7d48f --- /dev/null +++ b/skills/ag_ellipse_fitting.md @@ -0,0 +1,634 @@ +--- +name: ag_ellipse_fitting +description: Measure a galaxy's morphology non-parametrically by fitting isophotes — `ag.Ellipse` fitted to the data through `ag.FitEllipse` and `ag.AnalysisEllipse`, one ellipse at a time at fixed `major_axis`, producing radial profiles of axis ratio and position angle instead of a light-profile model. Covers how the likelihood differs (data interpolated onto the ellipse, residual = flux minus its own mean, no PSF and no model image), why Dynesty rather than Nautilus and why `use_jax=False` is required, the fit-the-centre-then-step-outwards workflow with the centre pinned from the first fit, `ag.EllipseMultipole` for m=1 lopsidedness / m=3 tripole / m=4 boxy-versus-discy deviations, combining every ellipse with `af.Drawer`, `aplt.subplot_fit_ellipse`, and reading many fits back with `ag.agg.EllipsesAgg` / `MultipolesAgg` / `FitEllipseAgg`. Use when the science is the isophote shape itself, or when no parametric model fits and you need a model-independent measurement. Not for light-profile modelling (`ag_build_imaging_model`, `ag_basis_profiles`), not for pixelized reconstruction (`ag_pixelization`), and note there is no `start_here.py` in the workspace's `ellipse/` package — `modeling.py` is the entry point. +--- + +# Measuring isophotes with ellipse fitting + +Every other fitting skill in this workspace builds a model image and subtracts it from the data. +Ellipse fitting does not. It puts an ellipse on the sky, samples the data around it, and asks a +much simpler question: **are these flux values all the same?** An isophote is by definition a +contour of constant surface brightness, so an ellipse that traces one will find the same flux all +the way round. One that does not is either the wrong shape, the wrong orientation, or centred in +the wrong place. + +Fit a ladder of ellipses at increasing `major_axis` and what comes out is not a model but a set of +*profiles*: axis ratio against radius, position angle against radius, and — with multipoles — the +deviation from an ellipse against radius. Those are the classical measurements of galaxy +structure. An isophotal **twist** (position angle rotating with radius) indicates a triaxial +system or a bar; ellipticity rising outwards distinguishes a disk emerging from a bulge; and the +m=4 multipole amplitude is the boxy/discy parameter that separates the two families of +early-type galaxies. None of these are things a Sersic profile can tell you, because a Sersic +asserts a single axis ratio and a single position angle at every radius by construction. + +The trade against a parametric fit is explicit. You gain a model-independent measurement, with no +functional form to be wrong about and nothing to deconvolve. You lose a total luminosity, a +size that extrapolates beyond the data, and any ability to decompose the light into components. +The two are complements: run an ellipse fit to *see* the structure, then use what you learn to +choose the parametric model — a measured twist, for instance, is the direct argument for a +two-basis MGE ([`ag_basis_profiles`](./ag_basis_profiles.md)). + +The concept page is +[`../wiki/core/concepts/ellipse_fitting_and_multipoles.md`](../wiki/core/concepts/ellipse_fitting_and_multipoles.md); +the API surface is [`../wiki/core/api/ellipse.md`](../wiki/core/api/ellipse.md). The grounding +scripts are `autogalaxy_workspace:scripts/ellipse/fit.py` (the likelihood, worked by hand), +`autogalaxy_workspace:scripts/ellipse/modeling.py` (the search-based workflow — **this is the +entry point; the `ellipse/` package has no `start_here.py`**), +`autogalaxy_workspace:scripts/ellipse/multipoles.py` and +`autogalaxy_workspace:scripts/ellipse/database.py`. + +## How the likelihood differs — read this before fitting anything + +Three properties of ellipse fitting surprise people, and all three follow from there being no +model image. + +**No PSF.** The dataset is loaded without one. Nothing is convolved, so nothing needs +deconvolving — you are measuring the observed isophotes, seeing included. That is a real +limitation for structure near the resolution limit and a real simplification everywhere else. + +**The residual is flux minus its own mean.** There is no model to subtract, so "residual" means +each interpolated flux value minus the mean of all the values on that ellipse. A good fit is one +where they agree; the goodness-of-fit measures scatter around the ellipse, not agreement with a +prediction. + +**No noise-normalisation term.** The usual log-likelihood carries a term from the noise map's +determinant. Here it is omitted, because interpolating the noise map onto the ellipse makes that +term numerically unstable. The consequence: **ellipse log-likelihoods are not comparable to +light-profile log-likelihoods**, and not comparable between ellipses of different `major_axis` +(different numbers of sampled points). Use them to compare shapes at one radius, not across radii. + +Worth doing once, by hand, because the objects then make sense: + +```python +""" +Ellipse Fitting: The Likelihood By Hand +======================================= + +Reproduce `ag.FitEllipse`'s likelihood arithmetic explicitly, to make concrete what "the data +interpolated onto an ellipse" means and why the residual is defined against the mean rather than +against a model. + +__Contents__ + +- **Imports:** Import the required libraries. +- **Dataset:** Load imaging without a PSF and mask it. +- **Interpolation:** Sample the data and noise map onto the ellipse's coordinates. +- **Likelihood:** Build the residual, chi-squared and log likelihood, and check against FitEllipse. +""" + +""" +__Imports__ +""" +from pathlib import Path + +import numpy as np + +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +DATASET_PATH = Path("dataset") / "imaging" / "my_galaxy" +PIXEL_SCALES = 0.1 +MASK_RADIUS = 5.0 + +""" +__Dataset__ + +No `psf_path`: ellipse fitting never convolves, so the PSF is not loaded +(`PyAutoArray:autoarray/dataset/imaging/dataset.py`). The mask radius does double duty here — it +excludes the noisy edges *and* sets how far out the ellipses can reach, so it is the single most +consequential number in the script. Keep it in a variable. +""" +dataset = ag.Imaging.from_fits( + data_path=DATASET_PATH / "data.fits", + noise_map_path=DATASET_PATH / "noise_map.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) + +""" +__Interpolation__ + +`ag.DatasetInterp` holds the interpolation weights and mappings once, so sampling many ellipses +is cheap (`PyAutoGalaxy:autogalaxy/ellipse/dataset_interp.py`). You do not normally construct it +yourself — `ag.FitEllipse` does it internally — but seeing it makes the mechanism obvious. + +`points_from_major_axis_from` returns (y, x) coordinates equally spaced along the ellipse, and it +needs `pixel_scale` because the *number* of points is chosen to match the number of data pixels +the ellipse crosses. A bigger ellipse therefore gets more points automatically, which is why the +likelihood is not comparable between radii. +""" +interp = ag.DatasetInterp(dataset=dataset) + +ellipse = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0) + +points = ellipse.points_from_major_axis_from(pixel_scale=dataset.pixel_scales[0]) + +data_interp = interp.data_interp(points) +noise_map_interp = interp.noise_map_interp(points) + +print(f"{len(points)} points sampled around a {ellipse.major_axis}\" ellipse") + +""" +__Likelihood__ + +The model data *is* the interpolated data — there is nothing else. The residual is each value +minus the mean of them all, which is the formal statement of "does this ellipse trace a contour +of constant brightness". Chi-squared and normalised residuals then follow the usual definitions +with the interpolated noise, and the log likelihood is -0.5 * chi-squared with no +noise-normalisation term (`PyAutoGalaxy:autogalaxy/ellipse/fit_ellipse.py`). +""" +model_data = data_interp +residual_map = data_interp - np.mean(data_interp) +normalized_residual_map = residual_map / noise_map_interp +chi_squared_map = (residual_map / noise_map_interp) ** 2.0 + +log_likelihood = -0.5 * np.sum(chi_squared_map) + +print(f"log likelihood by hand = {float(log_likelihood):.4f}") + +fit = ag.FitEllipse(dataset=dataset, ellipse=ellipse) + +print(f"FitEllipse.log_likelihood = {float(fit.log_likelihood):.4f}") +print(f"ratio to chi_squared = {float(fit.log_likelihood / fit.chi_squared):.4f}") +``` + +The two agree exactly. Note the factor: `log_likelihood` is **−0.5 × `chi_squared`**, verified +against the released library — the workspace's own prose in `ellipse/fit.py` and `ellipse/plot.py` +says −2.0, which is wrong by a factor of four. It does not change any inference (a constant +multiple of the log likelihood rescales but does not move the maximum), but it does change any +number you quote, so use −0.5. + +`ag.FitEllipse(dataset, ellipse, multipole_list=None, use_jax=False)` exposes all of it as +attributes — `data_interp`, `noise_map_interp`, `model_data`, `residual_map`, +`normalized_residual_map`, `chi_squared_map`, `chi_squared`, `log_likelihood` — so in practice you +construct the fit and read what you need. + +## `ag.Ellipse`, and one API trap + +`ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)` — three arguments, and the +ellipticity uses `ell_comps` rather than an axis-ratio/angle pair for the same reason every +profile does: a position angle is periodic, and a periodic parameter creates a boundary +pathology a non-linear search handles badly. + +You will want the axis ratio and position angle back out, and here is the trap: **`axis_ratio` and +`angle` are methods, not properties.** `minor_axis` is a property. Verified on the released stack: + +```python +ellipse = ag.Ellipse( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.5, angle=45.0), + major_axis=1.0, +) + +print(ellipse.minor_axis) # property -> 0.5 +print(ellipse.axis_ratio()) # CALL it -> 0.5 +print(ellipse.angle()) # CALL it -> 45.0 +``` + +Forgetting the parentheses does not raise — it prints ``, so a results table +built without them is silently full of method reprs rather than numbers. The workspace's +`ellipse/modeling.py` result section makes exactly this mistake in its own `print` statements; +`ag.convert.ell_comps_from(axis_ratio=..., angle=...)` is the inverse, for constructing an ellipse +from the values you would rather think in. + +## Ask + +- *"What are you measuring — a twist, an ellipticity profile, or boxiness?"* All three come from + the same ladder of fits, but boxiness needs multipoles, and that changes the model. +- *"How far out do you need to go?"* This sets the mask radius, which sets the ellipse ladder. + Going further than the data supports produces ellipses fitting noise, which show up as wild + swings in axis ratio at large radius. +- *"Is the centre known?"* The workflow below fits it once at small radius and then pins it. If + the galaxy is lopsided you may want it free at every radius — but expect a much noisier profile, + and consider an m=1 multipole instead, which is the physical way to say "lopsided". +- *"Are there contaminants inside the mask?"* An ellipse crossing a neighbour will average its + flux in and be pulled off the isophote. Mask them + ([`ag_light_model_extras`](./ag_light_model_extras.md)); `ellipse/modeling.py` repeats its whole + ladder with an extra-galaxies mask applied for exactly this reason. + +## Branch — the search: Dynesty, and no JAX + +Two settings that differ from every other fitting skill here, and both are deliberate. + +```python +""" +__Search__ + +Dynesty rather than Nautilus. Extensive testing has shown Dynesty with `sample="rwalk"` gives the +most accurate and efficient results for ellipse fitting specifically — a small, low-dimensional +parameter space with a likelihood surface unlike a light-profile fit's +(`PyAutoFit:autofit/non_linear/search/nest/dynesty/`). `n_live=50` is ample for N=4; the default +of 200 is more than these models need. + +`iterations_per_quick_update=10000` is high because these fits are fast: writing results and +visualisation to disk would otherwise dominate the run time. +""" +search = af.DynestyStatic( + path_prefix=Path("ellipse"), + name="fit_start", + unique_tag="my_galaxy", + sample="rwalk", + n_live=50, + iterations_per_quick_update=10000, +) + +""" +__Analysis__ + +`use_jax=False` is **required**: ellipse fitting does not support JAX acceleration. The default is +`True`, so this argument is not optional — omitting it is a real error, not a performance choice +(`PyAutoGalaxy:autogalaxy/ellipse/model/analysis.py`). +""" +analysis = ag.AnalysisEllipse(dataset=dataset, use_jax=False) +``` + +The cost is low. A likelihood evaluation is ~0.04 s on a typical dataset, rising to 0.5–1.0 s at +high resolution — so an ellipse ladder is a laptop-scale job, not an HPC one, which is part of why +it is such a good first look at a galaxy. Search selection generally is +[`ag_configure_search`](./ag_configure_search.md) and +[`../wiki/core/api/searches.md`](../wiki/core/api/searches.md). + +## Branch — fit the centre, then step outwards + +This is the workflow. It is two stages because the centre is a global property of the galaxy while +the shape is a function of radius, and fitting them together at every radius wastes parameters and +produces a noisier profile. + +**Stage one: one small ellipse, centre free.** + +```python +""" +__Model — the centre ellipse__ + +A single ellipse at a small fixed `major_axis`, with its centre and ellipticity free: N=4. The +`major_axis` is *fixed*, not fitted — the whole method works by choosing a radius and asking what +shape best traces the isophote there, so a free size would have nothing to constrain it. + +The model uses a list even for one ellipse, because that is the shape the analysis expects and it +generalises to many (`PyAutoFit:autofit/mapper/prior_model/collection.py`). + +Priors: the centre within 0.1" of the image centre, which is true if the data was cut out around +the galaxy; `ell_comps` over [-0.6, 0.6], which spans essentially every realistic ellipticity. +""" +ellipse = af.Model(ag.Ellipse) + +ellipse.centre.centre_0 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) +ellipse.centre.centre_1 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + +ellipse.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.6, upper_limit=0.6) +ellipse.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.6, upper_limit=0.6) + +ellipse.major_axis = 0.3 + +model = af.Collection(ellipses=[ellipse]) + +print(model.info) + +result = search.fit(model=model, analysis=analysis) + +centre = result.instance.ellipses[0].centre +print(f"centre = ({centre[0]:.4f}, {centre[1]:.4f})") +``` + +If the galaxy is not near (0.0", 0.0") the defaults will fight you. Either re-cut the data around +it ([`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md)) or widen those centre priors. + +**Stage two: a ladder, centre pinned.** + +```python +""" +__Multiple Ellipses__ + +Ten ellipses from 0.3" out to 90% of the mask radius, each with the centre *fixed* to the value +from stage one and its own `ell_comps` free. That drops each model from N=4 to N=2, and it is what +makes the resulting profiles clean: the axis ratio and position angle at each radius are then +measured independently of any centring uncertainty. + +Each radius is a separate `search.fit` with its own `name`, so each gets its own output folder and +resumes independently — a ladder that dies halfway restarts where it stopped. + +Stopping at 0.9 * mask_radius rather than at the mask edge keeps every ellipse fully inside the +unmasked region; an ellipse that runs off the mask samples nothing there and returns a shape +driven by which part of it survived. +""" +import numpy as np + +number_of_ellipses = 10 +major_axis_list = np.linspace(0.3, MASK_RADIUS * 0.9, number_of_ellipses) + +result_list = [] + +for i, major_axis in enumerate(major_axis_list): + + ellipse = af.Model(ag.Ellipse) + + ellipse.centre.centre_0 = centre[0] + ellipse.centre.centre_1 = centre[1] + + ellipse.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.6, upper_limit=0.6) + ellipse.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.6, upper_limit=0.6) + + ellipse.major_axis = major_axis + + model = af.Collection(ellipses=[ellipse]) + + search = af.DynestyStatic( + path_prefix=Path("ellipse"), + name=f"fit_{i}", + unique_tag="my_galaxy", + sample="rwalk", + n_live=50, + number_of_cores=4, + iterations_per_quick_update=10000, + ) + + result_list.append(search.fit(model=model, analysis=analysis)) + +""" +__Profiles__ + +The radial profiles are the result. Remember the parentheses on `axis_ratio()` and `angle()`. +""" +for result in result_list: + e = result.instance.ellipses[0] + print(f"a = {e.major_axis:.2f}\" q = {e.axis_ratio():.3f} PA = {e.angle():.1f} deg") +``` + +Read the PA column down the page: monotonic drift is an isophotal twist, and a jump of ~90° is +usually the major and minor axes swapping when the isophote is nearly round — not a physical +twist. Read the q column: rising outwards is a disk emerging from a bulge. + +**Stage three: one combined fit, for the record.** Every ellipse in a single model, evaluated once, +so the whole ladder lives in one output folder and one aggregator entry: + +```python +""" +__Final Fit__ + +`af.Drawer` with `total_draws=1` does not search — it evaluates the model once and writes it out. +The ellipses are the *instances* from the ladder, so nothing is re-fitted; this exists to produce +one combined result and one combined visualisation +(`PyAutoFit:autofit/non_linear/search/mle/drawer/search.py`). + +`model.dummy_0` is required: a model with no free parameters at all cannot be constructed, so one +unused parameter is added to satisfy that. It carries no meaning — do not report it. +""" +ellipses = [result.instance.ellipses[0] for result in result_list] + +model = af.Collection(ellipses=ellipses) +model.dummy_0 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + +search = af.Drawer( + path_prefix=Path("ellipse"), + name="fit_all", + unique_tag="my_galaxy", + total_draws=1, +) + +result = search.fit(model=model, analysis=analysis) +``` + +The whole three-stage workflow is `autogalaxy_workspace:scripts/ellipse/modeling.py`. + +## Branch — multipoles: the deviation from an ellipse + +An ellipse cannot be boxy, discy, lopsided or three-fold symmetric. Real galaxies are all four. +`ag.EllipseMultipole` perturbs the ellipse's radius with an angular harmonic of order `m` and two +amplitude components, and the physics is entirely in `m`: + +| `m` | Perturbation | What it measures | +|---|---|---| +| 1 | monopole displacement | **lopsidedness** — an asymmetry between opposite sides, from a recent interaction or a warp | +| 3 | tripole | **three-fold** asymmetry, a genuine departure from point symmetry | +| 4 | quadrupole | **boxy versus discy** — the classical parameter separating the two early-type families | + +The m=4 amplitude is the one with the deepest literature behind it: boxy ellipticals are typically +more luminous, slow-rotating and radio-loud, discy ones the reverse, so its sign carries real +physical information about formation history. m=2 is absent from this list because an m=2 +perturbation *is* an ellipse — it is already in `ell_comps`. + +Directly, on a concrete fit: + +```python +""" +__Multipoles__ + +Each `EllipseMultipole` takes an order `m` and two `multipole_comps`, and any number of them can +perturb one ellipse at once — the white contour on a plot stops being an ellipse and becomes the +perturbed shape (`PyAutoGalaxy:autogalaxy/ellipse/ellipse/ellipse_multipole.py`). +""" +multipole_order_1 = ag.EllipseMultipole(m=1, multipole_comps=(0.05, 0.05)) +multipole_order_3 = ag.EllipseMultipole(m=3, multipole_comps=(0.05, 0.05)) +multipole_order_4 = ag.EllipseMultipole(m=4, multipole_comps=(0.05, 0.05)) + +fit = ag.FitEllipse( + dataset=dataset, + ellipse=ellipse_instance, + multipole_list=[multipole_order_1, multipole_order_3, multipole_order_4], +) +``` + +As model components, multipoles are a **separate top-level collection** from the ellipses, nested +one level deeper — `multipoles=[[multipole_3, multipole_4]]`, a list per ellipse of the multipoles +perturbing it: + +```python +""" +__Model — ellipse plus multipoles__ + +Two multipoles on one ellipse: N=6 (the ellipse's 2 free `ell_comps` with the centre pinned, plus +2 components each for m=3 and m=4). `m` itself is *fixed* — it is the harmonic order, not +something to infer. + +`GaussianPrior(mean=0.0, sigma=0.1)` centred on zero is the right prior: zero is "a perfect +ellipse", so the prior expresses "probably close to elliptical, but let the data say otherwise", +and a posterior clearly displaced from zero is then a detection. +""" +multipole_3 = af.Model(ag.EllipseMultipole) +multipole_3.m = 3 +multipole_3.multipole_comps.multipole_comps_0 = af.GaussianPrior(mean=0.0, sigma=0.1) +multipole_3.multipole_comps.multipole_comps_1 = af.GaussianPrior(mean=0.0, sigma=0.1) + +multipole_4 = af.Model(ag.EllipseMultipole) +multipole_4.m = 4 +multipole_4.multipole_comps.multipole_comps_0 = af.GaussianPrior(mean=0.0, sigma=0.1) +multipole_4.multipole_comps.multipole_comps_1 = af.GaussianPrior(mean=0.0, sigma=0.1) + +model = af.Collection(ellipses=[ellipse], multipoles=[[multipole_3, multipole_4]]) + +print(model.info) + +result = search.fit(model=model, analysis=analysis) + +print(result.instance.multipoles[0][0].multipole_comps) # m=3 +print(result.instance.multipoles[0][1].multipole_comps) # m=4 +``` + +Adapted from `autogalaxy_workspace:scripts/ellipse/multipoles.py`. Note the index shape: +`instance.multipoles[][]`. + +Two judgement calls. The workspace's ladder shares **one** set of multipole amplitudes across all +ellipses — two free parameters per order for the whole galaxy rather than per radius — which is a +common assumption and keeps the model tractable. There is literature showing multipoles vary +radially, and fitting them per radius is a straightforward extension of the ladder loop; it costs +two parameters per order per radius and needs the signal to support it. And multipoles increase +parameter-space degeneracy noticeably: the amplitudes trade against `ell_comps`, so expect wider +errors and check the corner plot rather than the point estimate. On data with no real multipole +signal the amplitudes go to values close to zero, which is the control worth running first. + +## Branch — plotting and reading many fits back + +Ellipse fitting has its own subplot, because its quantities are 1D arrays around a contour rather +than 2D images: + +```python +""" +__Plot__ + +`aplt.subplot_fit_ellipse` takes a *list* of fits and plots the data with every ellipse's contour +overlaid, alongside the 1D residuals as a function of position angle — the second panel is the +diagnostic, because a good fit is flat there and a bad one shows a clear sinusoid whose period +tells you which harmonic is missing (`PyAutoGalaxy:autogalaxy/ellipse/plot/`). + +It takes `output_path` / `output_format` but **not** `output_filename`, so give each figure its own +directory ([`../wiki/core/api/plotting.md`](../wiki/core/api/plotting.md) tabulates the split). +""" +PLOT_DIR = Path("scripts/scratch/my_galaxy/ellipse/") + +fit_list = [ + ag.FitEllipse(dataset=dataset, ellipse=result.instance.ellipses[0]) + for result in result_list +] + +aplt.subplot_fit_ellipse( + fit_list=fit_list, output_path=str(PLOT_DIR), output_format="png" +) + +print(f"Saved to: {PLOT_DIR.resolve()}") +``` + +A two-period sinusoid in the residual panel means the ellipticity is wrong; a four-period one +means you need an m=4 multipole. That is the fastest read on whether multipoles are worth adding. +`aplt.plot_array(array=dataset.data, ...)` for the data alone, and during a fit the +`VisualizerEllipse` attached to `ag.AnalysisEllipse` writes these figures to the output folder +automatically — controlled by the `fit_ellipse` entry in `config/visualize/plots.yaml`, so you can +switch figures on and off without touching code +(`autogalaxy_workspace:scripts/ellipse/plot.py`). + +A ladder produces ten-plus separate fits, and reading them back one folder at a time is +unpleasant. The aggregator has dedicated objects, and there are three because the model has three +kinds of component. From `autogalaxy_workspace:scripts/ellipse/database.py`: + +```python +""" +__Aggregator__ + +`add_directory` scrapes an output tree into a queryable SQLite database. `EllipsesAgg` then yields +the maximum-likelihood ellipses of every fit in it via generators, so a ladder of any length stays +memory-light (`PyAutoGalaxy:autogalaxy/aggregator/ellipse/`). + +Each generator yields a *list* — one entry per analysis in the fit — so the `[0]` below takes the +single analysis. Multipoles live in their own aggregator because they are their own model +component; `FitEllipseAgg` rebuilds whole `FitEllipse` objects and picks up any multipoles +automatically. +""" +DATABASE = Path("output") / "ellipse.sqlite" +if DATABASE.exists(): + DATABASE.unlink() + +agg = af.Aggregator.from_database(filename=DATABASE.name, completed_only=False) +agg.add_directory(directory=Path("output") / "ellipse") + +ellipses_agg = ag.agg.EllipsesAgg(aggregator=agg) + +for ellipses_lists_list in ellipses_agg.max_log_likelihood_gen_from(): + for e in ellipses_lists_list[0]: + print(f"a = {e.major_axis:.2f}\" q = {e.axis_ratio():.3f} PA = {e.angle():.1f}") + +multipoles_agg = ag.agg.MultipolesAgg(aggregator=agg) + +for multipoles_lists_list in multipoles_agg.max_log_likelihood_gen_from(): + for multipole_list in multipoles_lists_list[0]: + print([m.m for m in multipole_list]) +``` + +`ag.agg.FitEllipseAgg` additionally offers `randomly_drawn_via_pdf_gen_from(total_samples=...)`, +which draws fits from the posterior rather than taking the maximum — that is how you put error +bars on a shape profile rather than quoting point estimates, and it is the right way to answer +"is this twist significant?". Aggregator patterns generally are +[`ag_load_results`](./ag_load_results.md) and +[`../wiki/core/api/aggregator.md`](../wiki/core/api/aggregator.md). + +## Combine — where this hands off + +- **Prepare the data and settle the mask** → [`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md). + The mask radius bounds the ellipse ladder, so it matters more here than anywhere. +- **Mask contaminating neighbours first** → [`ag_light_model_extras`](./ag_light_model_extras.md). + An ellipse crossing a companion averages its flux in and is pulled off the isophote; the + workspace's `modeling.py` re-runs its whole ladder with such a mask applied. +- **Turn a measured twist into a parametric model** → [`ag_basis_profiles`](./ag_basis_profiles.md). + An MGE with `gaussian_per_basis=2` is the parametric model that *can* represent what the ellipse + fit measured, and this is the cleanest reason to run both. +- **Turn a measured decomposition into components** → + [`ag_build_imaging_model`](./ag_build_imaging_model.md). An ellipticity profile rising outwards + is the argument for a bulge-plus-disk fit, and the ellipse fit gives you the priors. +- **Configure the search** → [`ag_configure_search`](./ag_configure_search.md), remembering that + the choices here (Dynesty, `sample="rwalk"`, `use_jax=False`) are ellipse-specific. +- **Run and watch the output folder** → [`ag_run_search`](./ag_run_search.md). +- **Read many fits back** → [`ag_load_results`](./ag_load_results.md). +- **A fit returns a wild axis ratio, or the PA jumps by 90°** → + [`ag_debug_fit_failure`](./ag_debug_fit_failure.md). The ellipse-specific causes are an ellipse + running off the mask, a contaminant crossed by the contour, a nearly-round isophote making the + axes ambiguous, and `use_jax` left at its `True` default. +- **Irregular structure no contour can trace** → [`ag_pixelization`](./ag_pixelization.md). + Ellipse fitting assumes the isophotes *are* closed contours around one centre; clumpy star + formation breaks that assumption outright. + +Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry recording the mask radius, +the `major_axis` ladder, the fitted centre, whether multipoles were included and at which orders, +and the q and PA profiles. Those profiles *are* the measurement, and unlike a parametric fit there +is no single `model.results` file that holds them. + +## Further reading + +- **General reference** — [RTD: Features overview](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_3_features.html): + ellipse fitting in the context of the other capabilities beyond a single smooth profile, with + links onward. +- **Experienced PyAutoGalaxy user** — [workspace: ellipse/modeling.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/ellipse/modeling.py): + the full three-stage workflow — centre fit, ladder, combined `Drawer` fit — and the repeat with + an extra-galaxies mask. This is the `ellipse/` package's entry point; there is no `start_here.py` + there. [multipoles.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/ellipse/multipoles.py) + is the multipole extension. + +(No HowToGalaxy chapter covers ellipse fitting — the lecture series teaches light-profile +modelling. For a newcomer, the concept page +[`../wiki/core/concepts/ellipse_fitting_and_multipoles.md`](../wiki/core/concepts/ellipse_fitting_and_multipoles.md) +is the place to start instead.) + +## Agent procedural checklist + +1. Set the mask radius consciously and keep it in a variable — it bounds the ellipse ladder, and + ellipses are stopped at `0.9 * mask_radius`. +2. Load the dataset **without** a PSF; ellipse fitting never convolves. +3. Use `af.DynestyStatic(sample="rwalk", n_live=50)` and `ag.AnalysisEllipse(..., use_jax=False)` — + `use_jax` defaults to `True` and must be overridden. +4. Fix `major_axis` on every ellipse; it is the radius you chose, never a fitted parameter. +5. Fit one small ellipse with the centre free, then pin that centre for the whole ladder (N=4 → 2). +6. Call `axis_ratio()` and `angle()` with parentheses — they are methods, and forgetting them + silently yields method reprs in your table. +7. Report the log likelihood as −0.5 × chi-squared, and never compare ellipse likelihoods across + radii or against a light-profile fit. +8. Add multipoles only with zero-centred Gaussian priors; fix `m`; index results as + `instance.multipoles[ellipse][multipole]`. +9. Read the 1D residual panel of `aplt.subplot_fit_ellipse` for period structure before adding a + harmonic; quote the plot's absolute path and offer to open it. +10. Combine the ladder with `af.Drawer(total_draws=1)` plus a `dummy_0` parameter, read profiles + back with `ag.agg.EllipsesAgg` / `MultipolesAgg`, and offer the `wiki/project/` entry holding + the q and PA profiles. diff --git a/skills/ag_light_model_extras.md b/skills/ag_light_model_extras.md new file mode 100644 index 0000000..92b97be --- /dev/null +++ b/skills/ag_light_model_extras.md @@ -0,0 +1,530 @@ +--- +name: ag_light_model_extras +description: Add the three model components that sit beside a galaxy's own light profiles — contaminating extra galaxies, a residual background sky via `ag.DatasetModel`, and operated (already-PSF-convolved) light profiles for compact nuclear emission. Covers when to mask a neighbour, scale its noise or model it and why the choice changes the answer; loading fixed centres from JSON with `ag.Grid2DIrregular`; `SersicSph` versus a pinned-centre MGE for a companion; why `af.Model(ag.DatasetModel)` starts with zero free parameters and silently fixes the sky at 0.0 unless you assign a prior; `grid_offset` and `grid_rotation_angle` for band registration; and `ag.lp_operated` / `ag.lp_linear_operated` for an AGN or nuclear starburst whose emission already shows the PSF. Worked against the bundled real dataset `dataset/imaging/cosj100020+015344`, which has both an un-subtracted sky pedestal and a faint neighbour 2.6" from the centre. Use once a galaxy model exists and the data has contaminants, an uncertain sky, or a point-like nucleus. Not for the galaxy's own morphology (`ag_build_imaging_model`, `ag_basis_profiles`), not for pixelized reconstruction (`ag_pixelization`), and not for two or more co-dominant galaxies. +--- + +# The parts of the model that are not the galaxy + +Three things routinely bias a galaxy-structure measurement, and none of them are the galaxy: a +neighbour whose light overlaps it, a background sky that was not perfectly subtracted, and a +compact nuclear source that a smooth profile cannot represent. All three are handled as extra +model components, and all three matter for the same reason — they contaminate the **faint outer +light**, which is precisely what sets `effective_radius` and `sersic_index`. Get them wrong and +the fit does not fail; it returns confident, wrong numbers. + +Statistically these are nuisance parameters. You do not care about the sky level or the +neighbour's size, but you care very much that their uncertainty propagates into the errors you +quote. A model that fixes the sky at zero because the reduction "should have" removed it asserts +zero uncertainty on something genuinely uncertain, and every parameter that has to absorb the +residual inherits a bias with no error bar to show it. Including a nuisance parameter you expect +to be small is how you *check* that belief rather than assume it. + +The three components sit in different places in the model tree, and that geography is worth +knowing before any code: + +| Component | Where it goes | Typical cost | +|---|---|---| +| extra galaxies | its own top-level `extra_galaxies` collection, **not** inside `galaxies` | ~2 free parameters each | +| background sky, astrometric offset | `dataset_model=af.Model(ag.DatasetModel)`, beside `galaxies` | 1–3 | +| operated profile | an ordinary component *on* the galaxy | 3–4 | + +## The worked case — the bundled dataset + +This skill uses `dataset/imaging/cosj100020+015344`, the real four-band JWST/NIRCam cutout that +ships with this repo, because it has two of the three problems for real: + +- **The sky is not subtracted.** `calwebb_image3`'s skymatch step matches the exposures' + backgrounds to each other but does not remove them, so `data.fits` carries the real JWST sky as + a positive pedestal — 4.5× the median noise at F115W, rising to 19× at F444W. A light-profile + fit that ignores this absorbs the pedestal into the profile wings and returns an inflated + effective radius and Sersic index. The measured value per band is in each `info.json` as + `background_sky_level` (0.1246 MJy/sr at F277W). +- **There is a real neighbour inside any useful mask.** A faint source sits **2.6" from the + centre**, holding about 0.3% of the galaxy's flux. A brighter one 8.0" out is excluded by any + mask under ~4", but the 2.6" source is not — it has to be masked or modelled every session. + **No `mask_extra_galaxies.fits` ships with the dataset**, so this is a live decision, not a + solved one. + +Both caveats, the `info.json` schema and the model-PSF warning are in +[`../wiki/core/operations/dataset.md`](../wiki/core/operations/dataset.md), with the full +provenance in the dataset's own `README.md`. The real-data inspection gate in +[`../AGENTS.md`](../AGENTS.md) applies to this dataset in full: plot it, show the user +`dataset.png`, and settle the neighbour and the mask extent from that look before composing +anything. If you have not done that, go to +[`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md) first. + +```python +""" +__Dataset__ + +The bundled cutout, F277W — the reference band, with the highest signal-to-noise (peak S/N 180) +and uniform exposure coverage across the frame. `info.json` carries every measured quantity, +including the pixel scale and the sky pedestal, so read it rather than hard-coding numbers +(`PyAutoArray:autoarray/dataset/imaging/dataset.py`). + +Note that `dataset.psf` is a `Convolver`, not an array — the kernel itself is +`dataset.psf.kernel`. And treat the PSF as the dominant systematic here: it is an STPSF *model* +kernel, slightly sharper than an empirical one, so expect a fitted size to compensate slightly +large. +""" +import json +from pathlib import Path + +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +DATASET_PATH = ( + Path("dataset") / "imaging" / "cosj100020+015344" / "wavebands" / "F277W" +) +MASK_RADIUS = 3.5 +REDSHIFT = 0.3422 # spectroscopic, from zCOSMOS-Bright DR3 — cited in the dataset README + +info = json.loads((DATASET_PATH / "info.json").read_text()) + +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=info["pixel_scale"], +) + +mask = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=dataset.pixel_scales, + radius=MASK_RADIUS, +) +dataset = dataset.apply_mask(mask=mask) +``` + +A mask of 3.5" keeps the fit inside the uniformly covered region and excludes the 8.0" +neighbour, while reaching the galaxy's outer isophotes. It does **not** exclude the 2.6" source — +that is what the rest of this skill is about. + +## Ask + +- *"Do you trust the sky subtraction?"* If the science is the faint outer envelope — an effective + radius, a Sersic index, a bulge-to-total ratio — the honest answer is usually "not enough to + fix it at zero". One parameter buys you the check. +- *"Is the neighbour's light genuinely blended with the target's, or merely nearby?"* Merely + nearby, and outside where the target's light matters → shrink the mask. Blended → scale its + noise or model it. The three strategies are not equivalent and the choice is scientific. +- *"Are you using a pixelization?"* If so, the masking option is off the table — removing pixels + punches holes in the mesh. Scale or model + ([`ag_pixelization`](./ag_pixelization.md)). +- *"Is there a compact nuclear source?"* A sharp central peak that a Sersic cannot reach without + driving `sersic_index` to its limit is usually an AGN or nuclear starburst. Two ways to model + it — an operated profile or a compact MGE — and they behave differently. +- *"Several bands, and are they registered?"* Sub-pixel misregistration between bands is what + `grid_offset` exists for. + +## Branch — the background sky + +One parameter, and the single highest-value addition to a fit that cares about outer light. + +```python +""" +__Model__ + +The galaxy's light plus the dataset's sky. `ag.DatasetModel` is not a galaxy and does not live in +the `galaxies` collection — it describes the *data*, and sits beside `galaxies` as its own +top-level component (`PyAutoArray:autoarray/dataset/dataset_model.py`). + +**You must set the prior on `background_sky_level` by hand.** Unlike a light profile's priors, +the right range depends entirely on the data's units and depth, so no default could be correct — +and `af.Model(ag.DatasetModel)` therefore starts with **zero** free parameters, holding every +constructor argument at its default. Omit the assignment below and the sky is silently *fixed* +at 0.0, which is the failure this branch exists to prevent. Confirm with +`print(dataset_model.total_free_parameters)`. + +Here the measured pedestal is 0.1246 MJy/sr, so a uniform prior from 0.0 to 0.4 comfortably +brackets it without letting the sky wander into the galaxy's flux. +""" +bulge = af.Model(ag.lp_linear.Sersic) +galaxy = af.Model(ag.Galaxy, redshift=REDSHIFT, bulge=bulge) + +dataset_model = af.Model(ag.DatasetModel) +dataset_model.background_sky_level = af.UniformPrior( + lower_limit=0.0, upper_limit=0.4 +) + +model = af.Collection( + galaxies=af.Collection(galaxy=galaxy), dataset_model=dataset_model +) + +print(f"DatasetModel free parameters = {dataset_model.total_free_parameters}") +print(model.info) +``` + +N = 7: the linear Sersic's six plus the sky. Adapted from +`autogalaxy_workspace:scripts/imaging/features/sky_background/modeling.py`, whose dataset is +simulated with a known pedestal of 5.0 electrons per second and recovers it — a useful control if +you want to convince yourself the mechanism works before trusting it on real data. + +Two traps worth stating explicitly, both verified: + +- `af.Model(ag.DatasetModel).total_free_parameters` is **0**. There is no warning; the fit runs + and reports a perfect sky subtraction you never tested. +- `background_sky_level` is **not a class attribute**. It is a constructor parameter, so reading it + off the `ag.DatasetModel` class itself raises `AttributeError`, and the only way to free it is to + assign a prior on an `af.Model` of that class, as above. + +After the fit, the inferred value is on the result's instance: + +```python +print(result.instance.dataset_model.background_sky_level) +``` + +Compare it against `info["background_sky_level"]`. Agreement is a real validation of the whole +fit; a large disagreement means the sky and the profile wings are trading, which is the +degeneracy this parameter exists to expose. The physics of why that degeneracy is so strong for +low-surface-brightness features is +[`../wiki/core/concepts/sky_background_and_operated_profiles.md`](../wiki/core/concepts/sky_background_and_operated_profiles.md). + +The alternative is to subtract `info["background_sky_level"]` from the data before fitting. That +is legitimate and cheaper, but it asserts the measurement is exact — you lose the error +propagation, which is the main thing modelling it buys. + +### The other two `DatasetModel` fields + +`ag.DatasetModel(background_sky_level=0.0, grid_offset=(0.0, 0.0), grid_rotation_angle=0.0)` — +the full surface, verified. The other two describe the data's *astrometry* rather than its +background, and earn their keep in multi-band work where the bands are not perfectly registered: + +```python +dataset_model = af.Model(ag.DatasetModel) +dataset_model.grid_offset.grid_offset_0 = af.GaussianPrior(mean=0.0, sigma=0.1) +dataset_model.grid_offset.grid_offset_1 = af.GaussianPrior(mean=0.0, sigma=0.1) +dataset_model.grid_rotation_angle = af.UniformPrior(lower_limit=-5.0, upper_limit=5.0) +``` + +Three parameters. Only free these when you have a reason: a sub-pixel offset between bands, or a +known small rotation. On a single band with the galaxy already centred, they add dimensions for +nothing. The bundled dataset's four bands agree on the centre to 0.08", which is a real offset at +0.03"/pixel — so this is exactly the case where a joint multi-band fit wants `grid_offset` free +(the multi-dataset skill, `ag_multi_dataset`, owns that). + +## Branch — extra galaxies: three strategies, one decision + +A neighbour inside the mask leaves you three options, and they are ordered by cost and by how +much of the target's light they preserve: + +1. **Shrink the mask.** Free. Correct when the neighbour lies outside the radius where your + target's light matters. On the bundled dataset this handles the 8.0" source completely. + It cannot handle the 2.6" one without truncating the galaxy's own outer isophotes — which + biases exactly the parameters you are measuring. +2. **Scale its noise.** One extra file, no parameters. Keeps the pixels in the fit but zeroes + their data and inflates their noise so they contribute nothing to the likelihood. This is the + right default when the neighbour is faint and its light does not deeply overlap the target's, + and it is **mandatory** if you are using a pixelization. +3. **Model it.** Two parameters per neighbour. The only option that recovers the target's light + *underneath* the neighbour, so the only correct one when the two genuinely blend. + +### Scaling the noise + +```python +""" +__Noise Scaling__ + +The mask is loaded with `invert=True` because in a `mask_extra_galaxies.fits` file `True` means +"scale this pixel" — the opposite of a modelling mask's convention. Getting it backwards scales +the galaxy instead of the neighbour, so always look at the signal-to-noise panel afterwards: the +scaled pixels should be visibly blank. Scale first, then apply the modelling mask +(`PyAutoArray:autoarray/dataset/imaging/dataset.py`). + +No such mask ships with the bundled dataset — draw one over the 2.6" source with +`autogalaxy_workspace:scripts/imaging/data_preparation/gui/mask_extra_galaxies.py`. +""" +mask_extra_galaxies = ag.Mask2D.from_fits( + file_path=DATASET_PATH / "mask_extra_galaxies.fits", + pixel_scales=info["pixel_scale"], + invert=True, +) + +dataset = dataset.apply_noise_scaling(mask=mask_extra_galaxies) +dataset = dataset.apply_mask(mask=mask) + +aplt.subplot_imaging_dataset( + dataset=dataset, + output_path="scripts/scratch/cosj100020/", + output_filename="noise_scaled", + output_format="png", +) +``` + +The scaled pixels reach a noise value of order 1e8 — that is what "contributes negligibly" means +in practice, and it is unmistakable on the subplot. Adapted from +`autogalaxy_workspace:scripts/imaging/features/extra_galaxies/modeling.py`. + +### Modelling it + +The convention is that a neighbour's **centre is fixed** to a value measured from the data, +leaving its other parameters free. This is not laziness: a companion with a free centre is a +model too complex to fit reliably, and the classic failure is one component wandering off to +absorb part of the target instead. + +```python +""" +__Extra Galaxies__ + +Centres come from a JSON written during data preparation, wrapped in `ag.Grid2DIrregular` so they +behave as a coordinate list (`PyAutoArray:autoarray/structures/grids/irregular_2d.py`). Each +neighbour becomes its own `ag.Galaxy` with a spherical linear Sersic whose centre is pinned — +`SersicSph` rather than `Sersic` because a faint companion rarely justifies two ellipticity +parameters, and its `intensity` is solved by the inversion. That leaves `effective_radius` and +`sersic_index`: **two free parameters per neighbour**, cheap enough for a handful of companions. + +`extra_galaxies` is its own top-level collection, **not** a member of `galaxies` — the analysis +relies on that placement (`PyAutoFit:autofit/mapper/prior_model/collection.py`). +""" +extra_galaxies_centres = ag.Grid2DIrregular( + 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=REDSHIFT, 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), + dataset_model=dataset_model, +) + +print(model.info) +``` + +For an irregular or asymmetric companion, swap the spherical Sersic for an MGE with a pinned +centre. It costs the **same** two free parameters in the linear limit while being far more +flexible — verified, and the reason the workspace recommends it once the number of companions +grows beyond a handful: + +```python +extra_galaxies_list = [] + +for extra_galaxy_centre in extra_galaxies_centres: + + mge_bulge = ag.model_util.mge_model_from( + mask_radius=MASK_RADIUS, + total_gaussians=10, + centre_fixed=tuple(extra_galaxy_centre), + ) + extra_galaxies_list.append( + af.Model(ag.Galaxy, redshift=REDSHIFT, bulge=mge_bulge) + ) +``` + +Both options come straight from +`autogalaxy_workspace:scripts/imaging/features/extra_galaxies/modeling.py`, which presents them +as equally supported alternatives — Option A and Option B in its own prose. The MGE route is +[`ag_basis_profiles`](./ag_basis_profiles.md). + +### Two things people forget + +**Enlarge the mask.** If you are modelling a neighbour, its light has to be *in* the fit. The +workspace example uses a 6.0" mask against the 3.0" of its sibling scripts for exactly this +reason. Modelling a companion whose pixels you masked out fits nothing and costs two parameters. + +**Over-sample at every centre.** The adaptive over-sampling helper takes a `centre_list`, and it +needs each neighbour's centre as well as the target's — otherwise the companion's steep central +gradient is evaluated on a coarse grid and its light is systematically mis-computed: + +```python +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=dataset.grid, + sub_size_list=[4, 2, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0)] + extra_galaxies_centres.in_list, +) +dataset = dataset.apply_over_sampling(over_sample_size_lp=over_sample_size) +``` + +When to pick which strategy, in full: +[`../wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md). + +One thing this feature deliberately does **not** have: a scaling-relation tier tying many faint +companions' properties to their luminosities through a few shared parameters. For a light-only +fit that is degenerate by construction — linear light profiles already solve `intensity`, so +scaling intensity by another shared free parameter adds nothing. The workspace's +`features/extra_galaxies/README.md` has the full argument. If several galaxies are genuinely +co-dominant rather than contaminating, that is a different model altogether — the +multi-galaxy and cluster skill (`ag_multi_galaxy_and_cluster`). + +## Branch — operated light profiles, for emission that already shows the PSF + +Compact central emission — an AGN, an unresolved nuclear starburst — is genuinely hard to model +the ordinary way. The model image of a near-point source is convolved with the PSF, and the +result is acutely sensitive to which pixel (and which sub-pixel) the emission lands in. The +likelihood surface acquires structure on the scale of a pixel, which is exactly the kind of +surface a sampler cannot map. + +An **operated** profile sidesteps the problem by assuming the profile is *already* PSF-convolved. +The convolution step is skipped, and the profile is fitted directly to the emission as observed — +PSF features included. That is both faster (no convolution) and far better behaved. + +```python +""" +__Model__ + +A linear Sersic bulge plus a linear operated `Gaussian` for the nucleus, with their centres +paired because a nuclear source sits at the galaxy's centre by definition +(`PyAutoGalaxy:autogalaxy/profiles/light/linear_operated/gaussian.py`). + +**Set the `sigma` prior by hand.** The default is a `UniformPrior` from 0.0 to 5.0, which is far +too wide: an operated Gaussian's width should be comparable to the PSF, so a fraction of an +arcsecond. Use what you know about your PSF — for this dataset the shipped F277W kernel has a +half-light radius of 0.074", so a prior up to ~0.5" is generous already. + +`ag.lp_linear_operated` rather than `ag.lp_operated` for the same reason as everywhere else: the +`intensity` is solved rather than sampled, which matters more here than usual because operated +profiles are the hardest components in this workspace to sample robustly. +""" +psf_component = af.Model(ag.lp_linear_operated.Gaussian) +psf_component.sigma = af.UniformPrior(lower_limit=0.0, upper_limit=0.5) + +bulge = af.Model(ag.lp_linear.Sersic) +bulge.centre = psf_component.centre + +galaxy = af.Model( + ag.Galaxy, redshift=REDSHIFT, bulge=bulge, point=psf_component +) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +print(model.info) +``` + +N = 9 — the bulge's six, plus the nucleus's `ell_comps` (2) and `sigma`, with the centre shared +and both intensities solved. The non-linear `ag.lp_operated.Gaussian` variant gives N = 10, since +its `intensity` is sampled. Adapted from +`autogalaxy_workspace:scripts/imaging/features/operated_light_profile/modeling.py`, which shows +both and recommends the linear one. + +`ag.lp_operated` and `ag.lp_linear_operated` each offer `Gaussian`, `Moffat` and `Sersic`. +`Gaussian` is the usual choice for a nucleus; `Moffat` if your PSF has strong wings. And the +mechanism is not limited to point sources — an operated profile is simply one that bypasses +convolution, so it fits any component you have reason to treat as already-convolved. + +There is a **competing approach**: a compact MGE, ten Gaussians on a `sigma` ladder capped at +twice the pixel scale, which is `ag.model_util.mge_point_model_from` in +[`ag_basis_profiles`](./ag_basis_profiles.md). Four free parameters, and it does not need a +`sigma` prior tuned by hand. The trade: the operated profile is one interpretable width you can +quote, the compact MGE is more flexible and better behaved if the nuclear emission is not quite +Gaussian. Both are documented side by side in +[`../wiki/core/concepts/sky_background_and_operated_profiles.md`](../wiki/core/concepts/sky_background_and_operated_profiles.md). +Ask if you want to fit both and compare. + +Note there is no `fit.py` example for operated profiles in the workspace — to do a single direct +fit, follow any other `imaging/fit.py` and swap the profile classes. Nothing else changes. + +## Branch — all three at once, on the real data + +They compose without interacting, which is the point of the model tree's geography: + +```python +""" +__Model__ + +The full model for the bundled cutout: the target's linear Sersic bulge with a linear operated +nuclear Gaussian, the 2.6" neighbour as a pinned-centre spherical Sersic, and the un-subtracted +sky as a `DatasetModel`. Three components in three different places in the tree — the galaxy's +own light inside `galaxies`, the companion in `extra_galaxies`, the sky in `dataset_model`. +""" +model = af.Collection( + galaxies=af.Collection(galaxy=galaxy), + extra_galaxies=af.Collection(extra_galaxies_list), + dataset_model=dataset_model, +) + +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}") +``` + +Twelve free parameters for the full model on this dataset, which any search handles comfortably. +Build it up rather than composing it in one go: fit the bulge alone, then add the sky, then the +neighbour, then the nucleus, and watch what each addition does to the residuals and to +`effective_radius`. That sequence is the evidence for the model, and it is what a referee will +ask you to show. + +Run times are effectively unchanged by any of these. The sky is a constant added to the data; +`SersicSph` and an operated `Gaussian` are cheap to evaluate, and the operated one is *faster* +than an ordinary profile because it skips convolution. The cost is dimensionality, not +arithmetic. + +## Combine — where this hands off + +- **Inspect the data and settle the mask** → [`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md). + On real data this is a gate, not a suggestion — and it is where `mask_extra_galaxies.fits` and + `extra_galaxies_centres.json` get made. +- **Compose the galaxy's own light** → [`ag_build_imaging_model`](./ag_build_imaging_model.md). +- **A basis for the target or for an irregular companion** → + [`ag_basis_profiles`](./ag_basis_profiles.md), which also holds the compact-MGE alternative to + an operated profile. +- **Pixelized reconstruction** → [`ag_pixelization`](./ag_pixelization.md). Noise scaling is + mandatory there; masking is not an option. +- **Configure and run** → [`ag_configure_search`](./ag_configure_search.md) and + [`ag_run_search`](./ag_run_search.md). +- **The sky runs to its prior boundary, or a companion absorbs the target** → + [`ag_debug_fit_failure`](./ag_debug_fit_failure.md). A sky pinned at its upper limit means the + prior is too narrow or the profile wings are being starved; a companion at an implausible size + usually means its centre was not fixed, or the mask cut off its light. +- **Several bands, with registration offsets** → the multi-dataset skill (`ag_multi_dataset`), + where `grid_offset` is fitted per band through the factor graph. The bundled dataset's four + bands are the worked case. +- **Two or more co-dominant galaxies rather than a target plus contaminants** → the multi-galaxy + and cluster skill (`ag_multi_galaxy_and_cluster`). +- **Add components progressively across chained searches** → the search-chaining skill + (`ag_chain_searches`). + +Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry recording which strategy +you used for each contaminant and **why**, the sky prior and the value recovered against the +measured one, and how `effective_radius` moved as each component was added. Those are the +modelling assumptions the result depends on, and none of them are visible in the output folder. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Two galaxies](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_3_search_chaining/tutorial_3_x2_galaxies.ipynb): + what happens to a fit when a second galaxy's light is in the frame, and how the model grows to + account for it. +- **General reference** — [RTD: Features overview](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_3_features.html): + the tour of the capabilities beyond a single smooth profile, including the sky background and + operated light profiles with links onward. +- **Experienced PyAutoGalaxy user** — [workspace: imaging/features/extra_galaxies/modeling.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/features/extra_galaxies/modeling.py): + masking, noise scaling and modelling in one script, with the `SersicSph` and MGE options side + by side. + +## Agent procedural checklist + +1. On real data, confirm the inspection gate was satisfied — plot the data, show the user the + path, and settle contaminants and mask extent — before composing anything. +2. On the bundled dataset, state the two live caveats: the un-subtracted sky and the 2.6" + neighbour with no mask shipped. +3. For the sky, assign a prior on `af.Model(ag.DatasetModel).background_sky_level` and **verify + `total_free_parameters` is not 0** before running. +4. Scale the prior to the data's units and depth — read `info["background_sky_level"]` rather + than guessing; after the fit, compare the inferred value against it. +5. For a neighbour, choose consciously between shrinking the mask, scaling noise and modelling — + and say why. With a pixelization, masking is off the table. +6. Fix every extra galaxy's centre; enlarge the mask to include its light; add its centre to the + over-sampling `centre_list`. +7. Load the `invert=True` scaling mask, then look at the signal-to-noise panel to confirm the + right pixels were scaled. +8. For a nucleus, prefer `ag.lp_linear_operated` and set the `sigma` prior from the PSF's + measured width — never leave it at the 0–5" default. +9. Build the model up one component at a time and record how `effective_radius` responds. +10. Save the script to `scripts/`, quote every plot's absolute path and offer to open it, then + offer the `wiki/project/` entry recording each contaminant decision. diff --git a/skills/ag_multi_dataset.md b/skills/ag_multi_dataset.md new file mode 100644 index 0000000..ade6118 --- /dev/null +++ b/skills/ag_multi_dataset.md @@ -0,0 +1,618 @@ +--- +name: ag_multi_dataset +description: Fit several datasets of the same galaxy jointly — multi-wavelength bands, repeated exposures in one band, or CCD imaging together with uv-plane visibilities. Covers the `af.AnalysisFactor` + `af.FactorGraphModel` construction that is the only way to combine datasets, deciding what is shared across datasets and what is freed per dataset, parameterising a wavelength relation with prior arithmetic so five bands cost two parameters instead of five, astrometric offsets via `ag.DatasetModel`, a pixelised component across bands, the deliberate one-dataset-at-a-time alternative, and reading the per-dataset result list. Use when more than one dataset of one galaxy must inform one model. Not for a single dataset (`ag_build_imaging_model`), not for one interferometer dataset on its own (`ag_build_interferometer_model`), not for choosing the search (`ag_configure_search`), and not for a population of different galaxies. +--- + +# Fitting several datasets of one galaxy + +A galaxy imaged in four filters is one object seen four ways. The geometry — where it sits, how +elongated it is, how concentrated its light is — is a property of the galaxy and does not change +between bands. The brightness does, and that changing brightness *is* the colour: the +bulge-to-disk ratio at 1.1 µm against 4.4 µm, a red concentrated old population sitting inside a +bluer star-forming disk, a colour gradient that traces an age or metallicity gradient. + +Fitting the bands separately throws away the thing that makes them powerful. Four independent +fits give you four loosely-constrained sizes; one joint fit with the geometry shared gives you +one well-constrained size and four amplitudes, and the amplitudes are the measurement. This is +also the statistically honest construction: the joint log likelihood is the sum of the +per-dataset log likelihoods, so every photon in every band constrains the shared parameters at +once, and the errors come out of a single posterior rather than being combined by hand +afterwards. + +The same machinery covers datasets that are not different wavelengths at all — undithered +exposures in one band fitted before they are drizzled together, CCD imaging alongside +sub-millimetre visibilities, multi-epoch follow-up. Anything where several measurements of one +galaxy must inform one model. + +Read [`../wiki/core/concepts/multi_wavelength.md`](../wiki/core/concepts/multi_wavelength.md) +for the *what is shared* decision in depth, and +[`../wiki/core/api/analysis_objects.md`](../wiki/core/api/analysis_objects.md) for the factor +graph's API surface. The canonical script is +`autogalaxy_workspace:scripts/multi_dataset/start_here.py`. + +## Ask + +- *"How many datasets, and what distinguishes them?"* Different filters, repeated exposures in + one filter, or genuinely different instruments — the answer picks the branch, because it + decides what is shared. +- *"Which parameters do you believe are the same across your datasets, and which must differ?"* + This is the whole model, and it is a scientific question before it is an API one. If the user + does not have a view yet, the default below (everything shared, amplitudes free by virtue of + linear light profiles) is the right place to start. +- *"Do you want them fitted simultaneously, or the best one first and the rest chained off it?"* + Simultaneous is the default. One-by-one is the deliberate alternative when one dataset is much + better than the others, and its branch is below. +- *"Are the datasets astrometrically aligned, and do you trust that alignment?"* Sub-pixel + residuals are normal even after a good reduction, and unmodelled they leak straight into the + structural parameters. + +If any dataset is **real observational data**, the inspection gate in +[`../AGENTS.md`](../AGENTS.md) applies to **every** band, not just the first: contaminants and +mask extent are settled from looking at each one, in +[`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md). A neighbour that is faint at 1.1 µm +can dominate at 4.4 µm. + +## Branch — the joint fit + +This is the deliverable: one script that loads several bands, states what is shared, and fits +them together. Adapted from `autogalaxy_workspace:scripts/multi_dataset/start_here.py` and +`autogalaxy_workspace:scripts/multi_dataset/modeling.py`. + +The example uses the four-band JWST/NIRCam cutout that ships with this repo +(`autogalaxy_assistant:dataset/imaging/cosj100020+015344/`, layout and `info.json` schema in +[`../wiki/core/operations/dataset.md`](../wiki/core/operations/dataset.md)) — a real early-type +galaxy at z = 0.3422 with two short-wave bands at 0.03"/pixel and two long-wave bands at +0.06"/pixel. Point `WAVEBAND_LIST` and `DATASET_PATH` at your own bands and nothing else +changes. + +```python +""" +Galaxy Structure: Joint Multi-Band Fit +===================================== + +Fit four bands of JWST/NIRCam imaging of one galaxy simultaneously: load and mask each band, +compose a single galaxy whose geometry is shared across every band, and sample the joint +posterior with Nautilus so the shared effective radius and Sersic index are constrained by all +four bands at once while each band solves its own amplitude. + +__Contents__ + +- **Imports:** JAX environment first, then the standard trio. +- **Dataset:** Load one `Imaging` per band, each with its own pixel scale. +- **Mask:** Mask and over-sample each band on the same sky footprint. +- **Model:** One galaxy, with the geometry shared and the amplitudes solved linearly. +- **Analysis:** One analysis per band. +- **Factor Graph:** Pair each analysis with a model copy and combine them. +- **Search:** Configure Nautilus for the joint likelihood. +- **Model-Fit:** Run the fit and announce the output folder. +- **Result:** Read the per-band results back. +""" +from autogalaxy import jax_wrapper # Sets the JAX environment before other imports + +import json +from pathlib import Path + +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +""" +__Dataset__ + +Each band is loaded as its own `ag.Imaging` (`PyAutoArray:autoarray/dataset/imaging/dataset.py`). +The bands are *not* stacked into one array: they have different PSFs and, here, different pixel +scales, and both of those are properties of the data rather than the galaxy. `pixel_scales` is +read from each band's `info.json` rather than hard-coded, because getting it wrong silently +rescales every length the fit reports. +""" +DATASET_PATH = Path("dataset") / "imaging" / "cosj100020+015344" +WAVEBAND_LIST = ["F115W", "F150W", "F277W", "F444W"] + +dataset_list = [] + +for waveband in WAVEBAND_LIST: + band_path = DATASET_PATH / "wavebands" / waveband + + info = json.loads((band_path / "info.json").read_text()) + + dataset_list.append( + ag.Imaging.from_fits( + data_path=band_path / "data.fits", + noise_map_path=band_path / "noise_map.fits", + psf_path=band_path / "psf.fits", + pixel_scales=info["pixel_scale"], + ) + ) + +""" +__Mask__ + +Use the same sky footprint for every band wherever you can. The radii need not agree in pixels — +they cannot, with two pixel scales — but making them agree in arcseconds means each band's fit +is constrained by the same region of the galaxy, which is what makes a colour measured across +bands meaningful. Over-sampling is applied per band, because it is defined on that band's grid +(`ag.util.over_sample.over_sample_size_via_radial_bins_from`). +""" +MASK_RADIUS = 2.2 + +dataset_masked_list = [] + +for dataset in dataset_list: + 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) + + dataset_masked_list.append(dataset) + +""" +__Model__ + +One galaxy, composed once. A linear `Sersic` bulge and a linear `Exponential` disk: because +these are `lp_linear` profiles, `intensity` is solved by linear inversion at every likelihood +evaluation rather than sampled, so each band gets its own amplitude for free and colour +gradients come out correctly with no per-band prior and no extra dimensions +(`PyAutoGalaxy:autogalaxy/profiles/light/linear/sersic.py`, +[`../wiki/core/concepts/linear_light_profiles_and_mge.md`](../wiki/core/concepts/linear_light_profiles_and_mge.md)). + +`ag.DatasetModel` is included so the astrometric offset of each band after the first can be +freed below (`PyAutoArray:autoarray/dataset/dataset_model.py`). +""" +bulge = af.Model(ag.lp_linear.Sersic) +disk = af.Model(ag.lp_linear.Exponential) + +galaxy = af.Model(ag.Galaxy, redshift=0.3422, bulge=bulge, disk=disk) + +dataset_model = af.Model(ag.DatasetModel) + +model = af.Collection( + dataset_model=dataset_model, galaxies=af.Collection(galaxy=galaxy) +) + +""" +__Analysis__ + +One analysis per band. Each holds one dataset and knows how to score a proposed model against +it (`PyAutoGalaxy:autogalaxy/imaging/model/analysis.py`). +""" +analysis_list = [ + ag.AnalysisImaging(dataset=dataset, use_jax=True) + for dataset in dataset_masked_list +] + +""" +__Factor Graph__ + +Each analysis is paired with a *copy* of the model in an `af.AnalysisFactor`, and the factors are +combined in an `af.FactorGraphModel` whose log likelihood is the sum over factors +(`PyAutoFit:autofit/graphical/declarative/factor/analysis.py`, +`PyAutoFit:autofit/graphical/declarative/collection.py`). + +Sharing is expressed entirely through the model, never by special-casing a dataset. With a bare +`model.copy()` and no overrides the graph identifies every prior across the factors, so the whole +model is shared. Overriding a prior on one copy before wrapping it is what frees that parameter +for that dataset — here the grid offset of every band after the first. +""" +analysis_factor_list = [] + +for i, analysis in enumerate(analysis_list): + model_analysis = model.copy() + + if i > 0: + model_analysis.dataset_model.grid_offset.grid_offset_0 = af.UniformPrior( + lower_limit=-0.2, upper_limit=0.2 + ) + model_analysis.dataset_model.grid_offset.grid_offset_1 = af.UniformPrior( + lower_limit=-0.2, upper_limit=0.2 + ) + + analysis_factor_list.append( + af.AnalysisFactor(prior_model=model_analysis, analysis=analysis) + ) + +factor_graph = af.FactorGraphModel(*analysis_factor_list, use_jax=True) + +print(factor_graph.global_prior_model.info) + +""" +__Search__ + +Multi-band fits stay on `Nautilus`. Each band may carry its own pixel scale, so the bands do not +share a grid shape and JAX must compile a separate gradient kernel per band; that compile cost +currently makes the gradient optimizers impractical here, even though they are the default for +single-dataset fits (`autogalaxy_workspace:scripts/multi_dataset/start_here.py` +`__Why Not MultiStartProdigy?__`). +""" +search = af.Nautilus( + path_prefix=Path("multi_dataset"), + name="four_band", + unique_tag="cosj100020+015344", + n_live=150, + n_batch=50, + iterations_per_quick_update=10000, + live_visual_update=False, +) + +""" +__Model-Fit__ + +Two things differ from a single-dataset fit: the *model* passed to the search is +`factor_graph.global_prior_model`, and the *analysis* is the factor graph itself. Results are +written on the fly, so the output folder is worth opening the moment the search starts. +""" +print(f"Output folder: {search.paths.output_path.resolve()}") + +result_list = search.fit( + model=factor_graph.global_prior_model, analysis=factor_graph +) + +""" +__Result__ + +One result per factor, in the order the factors were built. Shared parameters must be identical +across the results and freed ones must differ — reading them back is the only real confirmation +that the sharing you intended is the sharing you built. +""" +for waveband, result in zip(WAVEBAND_LIST, result_list): + instance = result.max_log_likelihood_instance + print( + f"{waveband}: effective_radius = " + f"{instance.galaxies.galaxy.bulge.effective_radius:.4f}, " + f"grid_offset = {instance.dataset_model.grid_offset}" + ) +``` + +`af.AnalysisFactor` + `af.FactorGraphModel` is **the** construction for combining datasets. +There is no other one, and nothing about the datasets is combined outside the model — if you +want a parameter to differ between two datasets, you override its prior on that dataset's model +copy, and that is the entire mechanism. + +Before launching, always `print(factor_graph.global_prior_model.info)`. It prints the shared +block once and then a numbered section per factor showing what that factor overrode. That +printout is the only way to confirm which priors were identified and which were freed, and a +graph that shares the wrong thing produces a confident wrong answer rather than an error. + +## Branch — deciding what is shared + +Scientific question first, API second. + +**Usually shared.** The galaxy's centre, its `ell_comps`, its structural scale (effective +radius, Sersic index), and the mesh and regularisation of any pixelised component. These +describe the object, so a band that wants its own value is usually telling you something about +your model rather than about the galaxy. + +**Usually per-dataset.** The PSF and pixel scale (already per-dataset by construction, since they +live on the `Imaging` object rather than in the model), the mask, the sky background level, and +any astrometric offset. + +**The interesting middle.** Amplitude *must* vary per band — that is colour, and linear light +profiles give it to you at zero cost. Whether the *radial scale* may vary is a real question +about the galaxy: an old concentrated population inside a bluer disk genuinely has a +wavelength-dependent half-light radius. Freeing it is legitimate; so is tying it to a relation +(next branch); so is fixing it and checking the residuals for the size mismatch you would then +expect to see. + +The failure modes are symmetric and neither announces itself. Give a band too much freedom and +it will quietly absorb a modelling failure that a shared parameter would have exposed as +residuals. Give it too little and you will attribute a genuine colour gradient to noise. + +To free one parameter per dataset, the workspace's own idiom is a plain override on the copy +(`autogalaxy_workspace:scripts/multi_dataset/modeling.py` `__Analysis Factor__`): + +```python +for analysis in analysis_list: + model_analysis = model.copy() + model_analysis.galaxies.galaxy.bulge.effective_radius = af.UniformPrior( + lower_limit=0.0, upper_limit=10.0 + ) + + analysis_factor_list.append( + af.AnalysisFactor(prior_model=model_analysis, analysis=analysis) + ) +``` + +## Branch — a wavelength relation instead of a parameter per band + +Freeing a parameter per band scales badly: five bands means five free effective radii, each +constrained by one band's data. If you expect the parameter to vary *smoothly* with wavelength — +which for a size, an axis ratio or a Sersic index is usually a defensible physical expectation — +parameterise the **relation** and fit its coefficients instead. A linear relation costs two +parameters no matter how many bands you have. + +Adapted from +`autogalaxy_workspace:scripts/multi_dataset/features/wavelength_dependence/modeling.py`: + +```python +WAVELENGTH_LIST = [1.15, 1.50, 2.77, 4.44] # microns, one per band + +bulge_m = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) +bulge_c = af.UniformPrior(lower_limit=0.0, upper_limit=5.0) + +analysis_factor_list = [] + +for wavelength, analysis in zip(WAVELENGTH_LIST, analysis_list): + model_analysis = model.copy() + + model_analysis.galaxies.galaxy.bulge.effective_radius = ( + wavelength * bulge_m + ) + bulge_c + + analysis_factor_list.append( + af.AnalysisFactor(prior_model=model_analysis, analysis=analysis) + ) + +factor_graph = af.FactorGraphModel(*analysis_factor_list, use_jax=True) +``` + +The arithmetic on priors is real, not notation: `(wavelength * bulge_m) + bulge_c` builds a +derived prior (`PyAutoFit:autofit/mapper/prior/arithmetic/arithmetic.py`), so `m` and `c` are +the sampled parameters and each band's `effective_radius` is determined by them. The +`global_prior_model.info` for such a graph shows each factor's `effective_radius` as a +`self`-with-`wavelength` entry rather than a prior, which is how you confirm the relation was +built rather than four independent radii. + +The scientific payoff is that `m` is *directly* the gradient of size with wavelength, measured +with an error bar, which is usually the quantity the paper wants — rather than four loosely +constrained radii you have to regress afterwards with no covariance information. + +The cost is the assumption. If the true variation is not the shape you imposed, the fit distorts +other parameters to compensate. Fit a small number of bands with free per-band values first, look +at whether they lie on a line, and then impose one. Any functional form you can write with prior +arithmetic works — the linear case is the example, not the limit. + +## Branch — astrometric offsets + +Reduction pipelines align frames, but sub-pixel residuals survive, and telescope pointing +precision leaves an uncertainty of its own. Left unmodelled that misalignment leaks into the +galaxy parameters, and the usual symptom is a model inferred from one band fitting a second band +noticeably worse than it should. + +`ag.DatasetModel` carries `grid_offset` (two parameters) and `grid_rotation_angle` +(`PyAutoArray:autoarray/dataset/dataset_model.py`). The offset is applied by shifting the +*image-pixel grid* before the light profiles are evaluated, not by moving the profile centres, so +the galaxy model's geometry is untouched. The default is not a prior but a fixed `(0.0, 0.0)`, so +an offset is opt-in per dataset. + +Convention, from +`autogalaxy_workspace:scripts/multi_dataset/features/dataset_offsets/modeling.py`: the first +dataset defines the reference frame and each subsequent one gets its own offset, so N datasets +cost `2 × (N − 1)` parameters. That is the price, and for most multi-band structural work it is +the difference between an accurate model and a subtly wrong one. + +The offset also has a wider prior range than you might expect in the workspace examples +(±1.0"); tighten it to what your reduction plausibly left behind, because a wide offset prior +gives the search room to trade position against structure. + +## Branch — several exposures at one wavelength + +The machinery does not care that the datasets are different filters. Several exposures in one +band — undithered HST frames fitted *before* drizzling, so that the drizzle's correlated noise +is never baked into the data you fit — use exactly the same construction with everything shared +except the offsets. `autogalaxy_workspace:scripts/multi_dataset/features/same_wavelength/modeling.py`. + +Dithered frames are deliberately shifted by design, typically by a fraction of a pixel, so this +case and the offset branch above almost always travel together. + +## Branch — CCD imaging together with visibilities + +Optical imaging and sub-millimetre visibilities of the same galaxy are complementary rather than +redundant: the imaging constrains the smooth stellar light, the interferometer the compact +star-forming structure, and the two need not resemble each other at all. Combining them is the +same factor graph with a different analysis class on one factor +(`autogalaxy_workspace:scripts/multi_dataset/features/imaging_and_interferometer/modeling.py`): + +```python +analysis_imaging = ag.AnalysisImaging(dataset=imaging, use_jax=True) + +analysis_interferometer = ag.AnalysisInterferometer( + dataset=interferometer, use_jax=True +) + +analysis_factor_list = [] + +for analysis in [analysis_imaging, analysis_interferometer]: + bulge = af.Model(ag.lp_linear.Sersic) + disk = af.Model(ag.lp_linear.Exponential) + + galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, disk=disk) + + model_analysis = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + + analysis_factor_list.append( + af.AnalysisFactor(prior_model=model_analysis, analysis=analysis) + ) + +factor_graph = af.FactorGraphModel(*analysis_factor_list, use_jax=True) +``` + +Note what the workspace script is honest about, and repeat it to the user: composing the model +*inside* the loop, as above, means the two factors share **nothing**, and a joint fit of two +completely independent models is no better than two separate fits. The construction is worth +running only once you have decided which parameters genuinely couple the two views — commonly +the centre and the geometry, rarely the amplitudes — and shared those explicitly. The +interferometer side (real-space mask, transformer choice, dirty-image diagnostics) is +[`ag_build_interferometer_model`](./ag_build_interferometer_model.md). + +## Branch — a pixelised component across bands + +A pixelised reconstruction works inside the graph unchanged. The mesh and the pixelisation are +shared, and the natural thing to free per band is the regularisation coefficient, since the +signal-to-noise ratio and the amount of real small-scale structure differ from band to band +(`autogalaxy_workspace:scripts/multi_dataset/features/pixelization/modeling.py`): + +```python +pixelization = af.Model( + ag.Pixelization, + mesh=af.Model(ag.mesh.RectangularAdaptDensity, shape=(30, 30)), + regularization=ag.reg.Constant, +) + +galaxy = af.Model(ag.Galaxy, redshift=0.5, pixelization=pixelization) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +for analysis in analysis_list: + model_analysis = model.copy() + model_analysis.galaxies.galaxy.pixelization.regularization.coefficient = ( + af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4) + ) + + analysis_factor_list.append( + af.AnalysisFactor(prior_model=model_analysis, analysis=analysis) + ) +``` + +Watch the memory: every band's linear-algebra structures live in VRAM simultaneously. What a +pixelisation is and how regularisation works is +[`../wiki/core/concepts/inversions_and_pixelizations.md`](../wiki/core/concepts/inversions_and_pixelizations.md), +and the single-dataset procedure is `ag_pixelization`. + +## Branch — one dataset at a time, deliberately + +A simultaneous fit is not always the right answer. +`autogalaxy_workspace:scripts/multi_dataset/features/one_by_one/modeling.py` fits the best +dataset first, then chains its inferred model into each of the others. Three situations justify +it: + +- **One dataset is much better than the rest.** A low-resolution band fitted simultaneously can + drag a joint fit toward a compromise that suits neither; fitting the good band first and + interpreting the poorer ones against it keeps the good constraint intact. +- **Total run time.** More searches, but each one is cheaper than a search over the joint + likelihood, and often faster in sum. +- **Robustness.** A structure that shows up in a simultaneous fit and vanishes when the bands + are fitted individually was probably not real. This is a cheap and underused check. + +How much to fix and how much to free at each step is the judgement call, and the script walks +three variants: fix the bulge entirely and free only the disk; free everything but start from +the first result; or fix the whole galaxy and fit *only* the offset between the two datasets, +which is the cleanest way to measure an alignment residual. The prior-passing mechanics belong +to [`ag_chain_searches`](./ag_chain_searches.md) — read it before writing the second search, +because the difference between passing a `model` and passing an `instance` is the difference +between narrowing a prior and removing a dimension. + +## Branch — reading the result, and the two traps in it + +`search.fit` on a factor graph returns a list-like result — one entry per factor, in the order +the factors were built — rather than a single `Result`. Two things surprise people: + +- **The `Samples` object is global.** It has the dimensionality of the whole graph and is + *identical* in every entry. `result_list[2].samples` is not "band 2's posterior"; there is only + one posterior, over the joint model. Per-band quantities come from + `result_list[i].max_log_likelihood_instance` and `result_list[i].max_log_likelihood_fit`. +- **Ordering is yours to keep straight.** Nothing labels a factor with its band. Zip your + waveband list against the results, as the script above does, rather than indexing by memory. + +Plot each band's fit separately — `aplt.subplot_fit_imaging` writes a fixed `fit.png`, so give +each band its own directory: + +```python +for waveband, result in zip(WAVEBAND_LIST, result_list): + aplt.subplot_fit_imaging( + fit=result.max_log_likelihood_fit, + output_path=f"scripts/scratch/four_band/{waveband}/", + output_format="png", + ) +``` + +The full functional plotting surface and which entry points accept `output_filename` is +[`../wiki/core/api/plotting.md`](../wiki/core/api/plotting.md) and +[`ag_plot_fit`](./ag_plot_fit.md). + +## Branch — cost, and proving the script before you pay it + +Every dataset adds its arrays to the likelihood and to VRAM. The factor graph carries the same +VRAM estimator the single-dataset analyses do +(`autogalaxy_workspace:scripts/multi_dataset/modeling.py` `__VRAM Use__`): + +```python +factor_graph.print_vram_use( + model=factor_graph.global_prior_model, batch_size=search.batch_size +) +``` + +It takes twenty or thirty seconds, so comment it out once you know your footprint. Batch size is +the lever in both directions: larger batches cut wall-clock time and raise VRAM. + +Smoke test before committing hours, exactly as for a single dataset — the graph adds new ways to +be wrong (a mis-scoped override, a factor built from the wrong model copy) and level 2 catches +all of them in seconds: + +```bash +PYAUTO_TEST_MODE=2 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/four_band_fit.py +``` + +Level 1 if you want the fit products written; the levels and the other short-circuit flags are +[`../wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md). Parameter values +from any test-mode run are meaningless — say so rather than quoting them. + +## Combine + +- [`ag_build_imaging_model`](./ag_build_imaging_model.md) — get one band fitting well before you + join four. A joint fit is not a way to rescue a model that fails on its best dataset. +- [`ag_chain_searches`](./ag_chain_searches.md) — the one-by-one route, and prior passing. +- [`ag_build_interferometer_model`](./ag_build_interferometer_model.md) — the visibility side of + a joint imaging-plus-interferometer fit. +- [`ag_load_results`](./ag_load_results.md) — pulling per-band amplitudes and colours out of the + result list, and the aggregator when you have many galaxies each with many bands. +- [`ag_plot_fit`](./ag_plot_fit.md) — per-band residuals, which is where a wrongly-shared + parameter shows itself. +- [`ag_configure_search`](./ag_configure_search.md) — `n_live`, `n_batch` and the resume + semantics, all of which apply unchanged to the joint fit. +- `ag_pixelization` and `ag_basis_profiles` — the components most often shared across bands. + +Population-level inference over *different* galaxies is a different construction on the same +graphical-model machinery; +[`../wiki/core/concepts/hierarchical_models.md`](../wiki/core/concepts/hierarchical_models.md) +covers it, and no skill owns it yet ([`../PENDING.md`](../PENDING.md)). + +When the joint fit is worth keeping, offer (default-yes) to record it in a dated +`wiki/project/YYYY-MM-DD-.md` entry: which bands, what was shared and why, what was freed, +and the output path — per [`_style.md`](./_style.md) property #5. The sharing decision is the +scientific content of a multi-band fit and is exactly the thing you will not remember in a +month. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Linear light profiles](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_2_modeling/tutorial_5_linear_profiles.ipynb): + why solving amplitudes by linear inversion rather than sampling them is what makes per-band + colours free of charge. There is no multi-wavelength chapter in the lecture series — this is + the tutorial that teaches the idea the joint fit leans on hardest. +- **General reference** — [RTD: Features](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_3_features.html): + the multi-wavelength section of the feature tour, with pointers onward. +- **Experienced PyAutoGalaxy user** — [workspace: multi_dataset/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/multi_dataset/start_here.py): + the minimal end-to-end joint fit, with the factor-graph and per-dataset-offset sections this + skill mirrors. + +## Agent procedural checklist + +1. Establish how many datasets there are and what distinguishes them. +2. On real data, confirm the inspection gate has been satisfied for **every** dataset. +3. Settle the shared-versus-per-dataset decision with the user before writing the model; state + the default (everything shared, amplitudes linear) if they have no view. +4. Build one analysis per dataset, wrap each in `af.AnalysisFactor` with its own `model.copy()`, + and combine with `af.FactorGraphModel`. +5. `print(factor_graph.global_prior_model.info)` and read it — confirm the identified and freed + priors are the ones intended. +6. Pass `factor_graph.global_prior_model` as the model and `factor_graph` as the analysis. +7. Where a parameter should vary smoothly with wavelength, offer the relation form before + freeing it per band. +8. Free `grid_offset` for every dataset after the first unless the user is confident in the + alignment. +9. Validate with `PYAUTO_TEST_MODE=2` before the production run; on a GPU also run + `factor_graph.print_vram_use`. +10. Announce the output path immediately, and zip the results against the dataset labels when + reading them back. +11. Offer the `wiki/project/` entry recording the sharing decision. diff --git a/skills/ag_multi_galaxy_and_cluster.md b/skills/ag_multi_galaxy_and_cluster.md new file mode 100644 index 0000000..200cd93 --- /dev/null +++ b/skills/ag_multi_galaxy_and_cluster.md @@ -0,0 +1,553 @@ +--- +name: ag_multi_galaxy_and_cluster +description: Model several galaxies whose light blends on the sky — an interacting or projected pair, a compact multiple, or a cluster field with a brightest cluster galaxy plus tens-to-hundreds of catalogued members. Covers the regime ladder and how to place a system on it, the list-based composition that gives every co-dominant galaxy its own free light model, why an MGE suits a blend, the catalogue-driven member tier loaded with `ag.galaxy_table_from_csv` whose intensities are tied to one shared free normalization so population size costs no dimensions, promoting bright members out of the tier, and per-galaxy decomposed photometry from `subplot_fit_imaging_of_galaxy`. Use when two or more galaxies are co-equal subjects of the fit. Not for one target with faint contaminating neighbours (`ag_light_model_extras`), not for several datasets of one galaxy (`ag_multi_dataset`), and not for a single galaxy (`ag_build_imaging_model`). +--- + +# Blended pairs and cluster fields + +Every galaxy-light model has to answer one question before anything else: how many galaxies in +this image are *subjects* of the fit, rather than things to be removed from it? PyAutoGalaxy +organises the answer as a ladder of three regimes, and choosing the wrong rung is the most +expensive mistake available in this part of the library. + +- **One galaxy.** It dominates the frame; anything else is a contaminant to mask out or model + with a restricted, centre-fixed profile. This is [`ag_build_imaging_model`](./ag_build_imaging_model.md) + and `ag_light_model_extras`. +- **Several blended galaxies.** Two or more of comparable brightness whose light overlaps — + an interacting pair, a close projected pair, a compact multiple. Each gets its own **free** + light model and they are fitted together, because where their light overlaps you cannot measure + one without the other. `autogalaxy_workspace:scripts/multi_galaxy/start_here.py`. +- **A cluster field.** A brightest cluster galaxy (BCG) plus tens to hundreds of member galaxies. + Giving each member a free model is neither possible nor desirable, so the population is driven + by a **catalogue** and the model's dimensionality stops growing with it. + `autogalaxy_workspace:scripts/cluster/start_here.py`. + +**The subject throughout is galaxy light.** The cluster workflow here models the surface +brightness of the BCG and its member population — that is its entire purpose. Nothing on this +ladder infers a mass distribution from the positions of background objects, and no member is +represented as an unresolved point of emission: every galaxy on every rung is a surface-brightness +profile evaluated on a real-space grid, exactly as a single galaxy is. + +Why it is worth doing: the assembly of BCG and intracluster light; member luminosity functions +measured from photometry that is *not* contaminated by the BCG's envelope; per-galaxy fluxes for +an interacting pair whose isophotes overlap; and a clean light model for any downstream analysis +that needs the galaxies' emission subtracted before it can see anything fainter. + +Read [`../wiki/core/concepts/galaxies.md`](../wiki/core/concepts/galaxies.md) for how several +galaxies compose into what gets fitted, and +[`../wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md) +for the boundary between a co-equal galaxy and a contaminant. + +## Ask + +- *"How many galaxies, and are they comparably bright?"* Comparable brightness and overlapping + light means the blended-pair rung. One dominant galaxy with faint companions means the + single-galaxy rung with extra galaxies — a different and much cheaper model. +- *"Do you have a catalogue?"* A photometry catalogue of centres and luminosities is what makes + the cluster rung tractable. Without one, a field of many galaxies has no tier to pin the faint + members and you are back to promoting each one by hand. +- *"Do you have the centres?"* Every rung above the first needs them, and they are what breaks + the labelling degeneracy between components in a blend. +- *"What is the deliverable — total photometry per galaxy, structural parameters for each, or a + member-subtracted BCG?"* The answer decides how much freedom each galaxy needs, and whether a + shared shape across a member tier is acceptable. + +On real data the inspection gate in [`../AGENTS.md`](../AGENTS.md) is doing more work here than +usual: the *number of subjects* and the mask extent are both decided from looking at the image, +and the mask must enclose every galaxy the model contains, not just the brightest. +[`ag_prepare_imaging_data`](./ag_prepare_imaging_data.md) owns that procedure. + +## Branch — a blended pair or multiple + +The deliverable is one script. Adapted from +`autogalaxy_workspace:scripts/multi_galaxy/start_here.py` and +`autogalaxy_workspace:scripts/multi_galaxy/modeling.py`. + +```python +""" +Galaxy Structure: A Blended Pair +=============================== + +Decompose the overlapping light of two galaxies of comparable brightness: mask the whole system, +give each galaxy its own free Multi-Gaussian Expansion initialized on its observed centre, fit +both simultaneously, and read each galaxy's flux out of the decomposition uncontaminated by its +neighbour. + +__Contents__ + +- **Imports:** JAX environment first, then the standard trio. +- **Dataset:** Load the imaging. +- **Galaxy Centres:** Read the centres that initialize each galaxy's priors. +- **Mask:** One mask enclosing every galaxy, over-sampled at each centre. +- **Model:** One free MGE per galaxy, composed in a loop. +- **Search:** Configure the search. +- **Model-Fit:** Fit both galaxies together and announce the output folder. +- **Result:** Per-galaxy decomposition and photometry. +""" +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__ +""" +DATASET_PATH = Path("dataset") / "imaging" / "" + +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=0.1, +) + +""" +__Galaxy Centres__ + +The centres are read from a small JSON file of `(y, x)` arcsecond pairs, one per galaxy +(`ag.from_json`, `PyAutoNerves:autonerves/dictable.py`). They are not free data — they are what +initializes each galaxy's centre prior, and in a blend that is load-bearing: two identical +components with identical broad priors are exchangeable, and the search will happily put both on +the brighter galaxy. Anchoring each on its observed position breaks that degeneracy. +""" +galaxy_centres = ag.from_json(file_path=DATASET_PATH / "galaxy_centres.json") + +""" +__Mask__ + +One mask enclosing **every** galaxy. The fit's job is to decompose the blend, so a mask drawn +around one galaxy would remove exactly the overlap region the decomposition depends on. +Over-sample at each galaxy's centre, not just at the origin. +""" +MASK_RADIUS = 3.0 + +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=[4, 2, 2], + radial_list=[0.3, 0.6], + centre_list=list(galaxy_centres), +) + +dataset = dataset.apply_over_sampling(over_sample_size_lp=over_sample_size) + +""" +__Model__ + +One galaxy per centre, each with its own free Multi-Gaussian Expansion, built in a loop and +keyed `galaxy_0`, `galaxy_1`, … The composition scales to any number of blended galaxies without +changing shape (`ag.model_util.mge_model_from`, +`PyAutoGalaxy:autogalaxy/analysis/model_util.py`). + +An MGE is the right basis for a blend for two reasons. Its Gaussians share the galaxy's centre +and ellipticity, so each galaxy contributes only a handful of non-linear parameters despite being +flexible enough for an irregular or interacting morphology — and in a blend the two galaxies' +parameters are partially degenerate wherever their light overlaps, so keeping the count down +matters more than usual. And its Gaussians are *linear* light profiles: their intensities are +solved exactly by linear inversion at every likelihood evaluation, so the flux ratio between the +two galaxies — the single most degenerate quantity in a blend — is solved rather than explored +stochastically ([`../wiki/core/concepts/linear_light_profiles_and_mge.md`](../wiki/core/concepts/linear_light_profiles_and_mge.md)). +""" +galaxy_dict = {} + +for i, centre in enumerate(galaxy_centres): + bulge = ag.model_util.mge_model_from( + mask_radius=MASK_RADIUS, + total_gaussians=20, + centre_prior_is_uniform=True, + centre=(centre[0], centre[1]), + ) + + galaxy_dict[f"galaxy_{i}"] = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model = af.Collection(galaxies=af.Collection(**galaxy_dict)) + +print(model.info) + +""" +__Search__ + +`Nautilus` for a posterior you can quote. The folder's `start_here.py` uses +`af.MultiStartProdigy` instead — a multi-start gradient optimizer that launches many independent +descents in parallel and returns the best. That population of starts matters here specifically: +several co-dominant galaxies make the parameter space strongly multi-modal, and a single starting +point lands in a local maximum often enough to be untrustworthy +(`PyAutoFit:autofit/non_linear/search/nest/nautilus/search.py`). +""" +search = af.Nautilus( + path_prefix=Path("multi_galaxy"), + name="pair", + unique_tag="", + n_live=150, + n_batch=50, + iterations_per_quick_update=1000, +) + +analysis = ag.AnalysisImaging(dataset=dataset, use_jax=True) + +""" +__Model-Fit__ +""" +print(f"Output folder: {search.paths.output_path.resolve()}") + +result = search.fit(model=model, analysis=analysis) + +""" +__Result__ + +`subplot_fit_imaging_of_galaxy` is the deliverable of a blended fit: for each galaxy it shows +that galaxy's modeled light and the data with the *other* galaxies subtracted. Each call writes +`of_galaxy_.png` into `output_path`, so give each galaxy its own directory. + +With the blend decomposed, per-galaxy photometry is direct — each galaxy's model image contains +only its own light, so summing it gives a flux uncontaminated by its neighbour. +""" +print(result.info) + +PLOT_DIR = Path("scripts") / "scratch" / "" + +aplt.subplot_fit_imaging( + fit=result.max_log_likelihood_fit, + output_path=str(PLOT_DIR / "fit"), + output_format="png", +) + +for i in range(len(galaxy_centres)): + aplt.subplot_fit_imaging_of_galaxy( + fit=result.max_log_likelihood_fit, + galaxy_index=i, + output_path=str(PLOT_DIR / f"galaxy_{i}"), + output_format="png", + ) + +for i, galaxy in enumerate(result.max_log_likelihood_galaxies): + image = galaxy.image_2d_from(grid=dataset.grids.lp) + print(f"galaxy_{i}: total model flux = {float(image.array.sum()):.3f}") + +print(f"Saved to: {PLOT_DIR.resolve()}") +``` + +Everything else about the fit is unchanged from the single-galaxy case: the same +`ag.AnalysisImaging`, the same searches, the same output folder. Only the model composition moved. +Every imaging feature applies per galaxy unchanged too — linear profiles, MGE variants, a sky +background, shapelets. + +## Branch — a cluster field + +The cluster rung changes *how the model is composed*, and the change is the interesting part. + +The BCG (and any other dominant galaxy) is modelled individually, exactly as a single galaxy +would be. The member population is driven by a catalogue: a CSV of `y, x, luminosity`, one row +per member, whose photometry pins the faint members so that only a **shared normalization** is +free. Adding a member is a row append that adds **zero** free parameters. + +```python +scaling_table = ag.galaxy_table_from_csv( + file_path=DATASET_PATH / "scaling_galaxies.csv" +) + +member_centres = scaling_table.centres.in_list +member_luminosities = scaling_table.luminosities + +bcg_centres = ag.from_json(file_path=DATASET_PATH / "bcg_centres.json") +``` + +`ag.galaxy_table_from_csv` (`PyAutoGalaxy:autogalaxy/galaxy/galaxy_table.py`) returns a +`GalaxyTable` whose `.centres` is a `Grid2DIrregular` — hence `.in_list` to get plain tuples — +and whose `.luminosities` are the catalogue values in whatever consistent units your photometry +uses. Only their *ratios* matter, because the absolute scale is absorbed by the single fitted +normalization, so magnitudes converted to relative fluxes work as well as calibrated ones. + +The mask has to be generous — the members span the frame — and every galaxy inside it must be in +the model, or its light lands in the residuals and biases the BCG's outer isophotes: + +```python +MASK_RADIUS = 11.0 + +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=[4, 2, 2], + radial_list=[0.3, 0.6], + centre_list=list(bcg_centres) + list(member_centres), +) + +dataset = dataset.apply_over_sampling(over_sample_size_lp=over_sample_size) +``` + +Then the two-tier composition, from `autogalaxy_workspace:scripts/cluster/start_here.py`: + +```python +# Tier 1 — the BCG, a free MGE like any single galaxy. + +bulge = ag.model_util.mge_model_from( + mask_radius=3.0, + total_gaussians=20, + centre_prior_is_uniform=True, + centre=(bcg_centres[0][0], bcg_centres[0][1]), +) + +galaxy_dict = {"bcg": af.Model(ag.Galaxy, redshift=0.5, bulge=bulge)} + +# Tier 2 — the catalogue members, sharing ONE free normalization. + +intensity_scale = af.UniformPrior(lower_limit=0.0, upper_limit=10.0) + +for i, (centre, luminosity) in enumerate(zip(member_centres, member_luminosities)): + bulge = af.Model(ag.lp.SersicSph) + bulge.centre = tuple(centre) # fixed to the catalogue + bulge.intensity = intensity_scale * float(luminosity) # tied to the catalogue + bulge.effective_radius = 0.6 # fixed shape + bulge.sersic_index = 3.0 + + galaxy_dict[f"member_{i}"] = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +model = af.Collection(galaxies=af.Collection(**galaxy_dict)) + +print(model.info) +``` + +The mechanism is one line: `bulge.intensity = intensity_scale * float(luminosity)`. Multiplying a +prior by a number builds a derived prior +(`PyAutoFit:autofit/mapper/prior/arithmetic/arithmetic.py`), so `intensity_scale` is the sampled +parameter and every member's intensity is determined by it — +`intensity_i = intensity_scale × luminosity_i`. The whole tier therefore costs exactly **one** +free parameter, whether it holds six members or two hundred. `print(model.info)` shows each +member's intensity as tied rather than as a prior, which is how you confirm you built a tier and +not a hundred independent galaxies. + +The scaling relation is the scientific assumption: it says the members' surface brightnesses +follow their catalogue luminosities up to one common factor. That is a statement about the member +population being a reasonably homogeneous family, and it is the reason the model is tractable at +all. Recovering `intensity_scale ≈ 1` when the catalogue luminosities *are* the true intensities +is the sanity check the workspace's simulated example is built around +(`autogalaxy_workspace:scripts/cluster/modeling.py` `__Result + Truth Comparison__`). + +### Refining the tier + +`autogalaxy_workspace:scripts/cluster/modeling.py` shows the two natural refinements, in +increasing cost. + +**Promote the shape to shared parameters.** Rather than fixing `effective_radius` and +`sersic_index`, let the whole tier share two free ones. The tier then costs three parameters +instead of one — still independent of population size: + +```python +intensity_scale = af.UniformPrior(lower_limit=0.0, upper_limit=10.0) +tier_effective_radius = af.UniformPrior(lower_limit=0.1, upper_limit=2.0) +tier_sersic_index = af.UniformPrior(lower_limit=0.5, upper_limit=5.0) + +for i in range(len(member_centres)): + bulge = af.Model(ag.lp.SersicSph) + bulge.centre = tuple(member_centres[i]) + bulge.intensity = intensity_scale * float(member_luminosities[i]) + bulge.effective_radius = tier_effective_radius + bulge.sersic_index = tier_sersic_index + + galaxy_dict[f"member_{i}"] = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) +``` + +**Promote individual members out of the tier.** Give the brightest members their own free models. +This costs their full per-galaxy parameter count each, so promote sparingly — brightest first, +and only while the data still constrains them: + +```python +N_PROMOTE = 2 + +for i in range(N_PROMOTE): + bulge = af.Model(ag.lp.SersicSph) + bulge.centre = tuple(member_centres[i]) + bulge.intensity = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + bulge.effective_radius = af.UniformPrior(lower_limit=0.1, upper_limit=3.0) + bulge.sersic_index = af.UniformPrior(lower_limit=0.5, upper_limit=5.0) + + galaxy_dict[f"member_{i}"] = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge) + +# ... remaining rows go into the shared tier as above, starting from N_PROMOTE. +``` + +The workspace's shipped CSV happens to be sorted brightest-first, so rows 0 and 1 are the two +brightest. Your own catalogue almost certainly is not — sort by luminosity, or select rows +explicitly, rather than assuming. + +The BCG's decomposed light is what the whole exercise is for: + +```python +aplt.subplot_fit_imaging_of_galaxy( + fit=result.max_log_likelihood_fit, + galaxy_index=0, + output_path="scripts/scratch//bcg/", + output_format="png", +) +``` + +Galaxy index 0 is the BCG because it was inserted into `galaxy_dict` first; the ordering is the +insertion order of that dictionary, so keep it deliberate. + +## Branch — blended galaxies *and* faint contaminants together + +The two tiers are not exclusive, and a real field usually needs both. A blended pair often sits in +a frame with fainter companions that are not subjects of the fit but do overlap the mask: + +- **Co-equal galaxies** go under `galaxies=af.Collection(galaxy_0=..., galaxy_1=...)` with full + free light models. +- **Contaminants** go under `extra_galaxies=af.Collection(...)` with a restricted light model + whose **centre is fixed**, or are removed from the likelihood by noise scaling. + +```python +extra_galaxies_centres = ag.Grid2DIrregular( + 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_dict), + extra_galaxies=af.Collection(extra_galaxies_list), +) +``` + +`autogalaxy_workspace:scripts/multi_galaxy/features/extra_galaxies/modeling.py`. Fixing the +centres matters *more* in a blended field than in a single-galaxy one, and the script says why: a +free-centre contaminant can wander, and here there is more than one bright thing for it to wander +onto. It may drift toward one of the co-equal galaxies and start absorbing light that the +decomposition is supposed to be measuring — quietly corrupting the one result the fit exists to +produce. + +Noise scaling is the alternative lever: `dataset.apply_noise_scaling(mask=mask_extra_galaxies)` +with a `mask_extra_galaxies.fits` inflates the noise over the contaminated pixels so they stop +constraining anything. The trade-off between modelling, masking and scaling is +[`../wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md), +and the single-galaxy procedure with the centre-marking GUI is `ag_light_model_extras`. + +Note the ordering consequence: `extra_galaxies` are appended after `galaxies`, so with two +co-equal galaxies and two contaminants the `galaxy_index` values are 0 and 1 for the pair and 2 +and 3 for the extras. + +## Branch — the failure modes, and how to see them + +Three things go wrong on these rungs, and none of them raises an exception. + +**Component swapping.** Two galaxies with identical models and broad centre priors are +exchangeable, so a search can put both on one galaxy and leave the other unmodelled — or swap +them between runs so results are not comparable. The fix is the one already in the recipe: anchor +each galaxy's centre prior on its observed position. If it still happens, tighten those priors. + +**A galaxy inside the mask that is not in the model.** Its light has nowhere to go but the +residuals, and the nearest modelled galaxy absorbs what it can — biasing that galaxy's outer +profile and hence its effective radius and total flux. Check the residual map for a coherent +positive blob at a position you did not model. + +**A tier that is really a wrong assumption.** If the member population is not a homogeneous +family, one shared normalization cannot fit it, and the symptom is systematic residuals that +correlate with member luminosity — bright members under-subtracted and faint ones over-subtracted, +or vice versa. That is the signal to promote the shape to shared free parameters, or to promote +the brightest members individually. + +All three are read off residual maps rather than parameter values, which makes +[`ag_plot_fit`](./ag_plot_fit.md) the tool and +[`ag_debug_fit_failure`](./ag_debug_fit_failure.md) the taxonomy. + +## Branch — cost, and proving it before you pay + +Model dimensionality is the thing to watch, and it is why the two rungs are composed differently: +each promoted galaxy costs its full parameter count, while the catalogue tier costs one (or three) +regardless of size. Count before you fit — `print(model.info)` and the `prior_count` — and be +suspicious of any total beyond a few tens of parameters. + +Smoke test first, always: + +```bash +PYAUTO_TEST_MODE=2 NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib \ + python scripts/cluster_fit.py +``` + +Level 2 evaluates the likelihood exactly once, which proves the catalogue loaded, the tying +worked, and every galaxy is on the grid. Level 1 if you want the fit products written. The levels +are [`../wiki/core/operations/sandbox.md`](../wiki/core/operations/sandbox.md). Never quote a +structural parameter or an `intensity_scale` from a test-mode run — the sampler was truncated and +the numbers are prior medians. + +On a GPU, `analysis.print_vram_use(model=model, batch_size=search.batch_size)` before a long run. +A large member tier adds many fixed-shape profiles which JIT-compile into the same batched +likelihood, so population size moves the cost far less than the number of *free* parameters does. + +## Combine + +- [`ag_build_imaging_model`](./ag_build_imaging_model.md) — get one galaxy fitting well first. A + blended fit will not rescue a model that fails on an isolated galaxy. +- `ag_light_model_extras` — the contaminant tier in full, including the centre-marking GUI and the + noise-scaling route. +- `ag_basis_profiles` — the MGE that every rung here leans on, and the shapelet alternative. +- `ag_pixelization` — when one of the blended galaxies is genuinely irregular and no basis fits. +- [`ag_chain_searches`](./ag_chain_searches.md) — fit the brightest galaxy first, fix it, then add + the rest. Often the most reliable route through a crowded field. +- [`ag_multi_dataset`](./ag_multi_dataset.md) — the same composition across several bands, where + each galaxy's colour becomes measurable. +- [`ag_load_results`](./ag_load_results.md) — per-galaxy parameters with errors, and the + aggregator for many fields. +- [`ag_simulate_dataset`](./ag_simulate_dataset.md) — build a blend with known truth and fit it, + which is the only way to know whether your decomposition is recoverable at your signal-to-noise + ratio. + +When the fit is worth keeping, offer (default-yes) to record it in a dated +`wiki/project/YYYY-MM-DD-.md` entry: which rung you placed the system on and why, the +catalogue used and what was tied to it, what was promoted, and the output path — per +[`_style.md`](./_style.md) property #5. The rung choice and the promotion decisions are the +scientific content of the fit. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Two-galaxy fits](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_3_search_chaining/tutorial_3_x2_galaxies.ipynb): + fitting two galaxies in one image from first principles, including why chaining searches is + often easier than fitting both at once. +- **General reference** — [RTD: New user guide](https://pyautogalaxy.readthedocs.io/en/latest/overview/overview_2_new_user_guide.html): + the decision tree that routes by scale of system — single galaxy, blended pair, cluster field — + which is this skill's ladder in the upstream docs' own words. +- **Experienced PyAutoGalaxy user** — [workspace: cluster/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/cluster/start_here.py): + the two-tier catalogue composition this skill mirrors; its sibling + [multi_galaxy/start_here.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/multi_galaxy/start_here.py) + is the blended-pair version. + +## Agent procedural checklist + +1. Place the system on the ladder explicitly, and say which rung and why, before composing + anything. +2. On real data, confirm the inspection gate has been satisfied — including that the mask + encloses every galaxy the model will contain. +3. Get the centres (and, for a cluster, the catalogue) before writing the model; do not invent + positions. +4. Compose in a loop keyed `galaxy_0`, `galaxy_1`, … for co-equal galaxies; use the two-tier + BCG-plus-catalogue form for a cluster. +5. Anchor every co-equal galaxy's centre prior on its observed position. +6. Over-sample at every galaxy's centre, not just the origin. +7. `print(model.info)` and confirm the tier is tied rather than free, and the total dimensionality + is what you intended. +8. Validate with `PYAUTO_TEST_MODE=2`; on a GPU also run `analysis.print_vram_use`. +9. Announce the output path at launch; plot `subplot_fit_imaging_of_galaxy` per galaxy into its + own directory and quote the path. +10. Check the residual map for an unmodelled galaxy and for luminosity-correlated tier residuals + before believing any parameter. +11. Offer the `wiki/project/` entry recording the rung and the promotion decisions. diff --git a/skills/ag_pixelization.md b/skills/ag_pixelization.md new file mode 100644 index 0000000..8881028 --- /dev/null +++ b/skills/ag_pixelization.md @@ -0,0 +1,600 @@ +--- +name: ag_pixelization +description: Reconstruct a galaxy's clumpy or irregular light directly on a regularized pixel mesh, using `ag.Pixelization` alongside a parametric `ag.lp_linear` bulge — the hybrid model that gives a smooth interpretable component plus a flexible flux map for everything the profile cannot explain. Covers when a pixelization is warranted and when it is overkill, choosing an `ag.mesh` and an `ag.reg` regularization scheme and what each costs in free parameters, why `mesh_shape` must be fixed before the fit, `over_sample_size_pixelization` rather than `over_sample_size_lp`, noise scaling instead of hard masking for contaminants, the positive-only solver, reading the reconstruction and its Bayesian-evidence terms out of the `Inversion`, exporting the reconstruction to CSV for collaborators, and the GPU/VRAM versus CPU run-time trade. Use when a Sersic, MGE or shapelet basis leaves structured residuals in discrete off-centre lumps. Not for morphology a basis can still fit (`ag_basis_profiles`), not for non-parametric isophote measurement (`ag_ellipse_fitting`), and not for the basics of model composition (`ag_build_imaging_model`). +--- + +# Reconstructing a galaxy on a pixel mesh + +Every model in this workspace so far has been a *function*: you asserted a functional form and +the fit found its parameters. A pixelization drops the functional form. It lays a mesh of pixels +over the galaxy, gives each pixel its own free flux value, and solves all of them at once by +linear inversion — thousands of amplitudes, none of them in the non-linear search's parameter +space. What stops that from simply fitting the noise is **regularization**: a prior that +penalises solutions where neighbouring mesh pixels differ sharply, so the reconstruction is as +smooth as the data will allow and no smoother. + +The science case is specific, and it is worth being honest about how specific. Galaxies with +spiral arms, asymmetric star-forming clumps, tidal features or low-surface-brightness +substructure have light that no smooth profile — and often no basis expansion either — +reproduces. The canonical use is therefore **hybrid**: a linear `Sersic` for the smooth central +bulge, and a pixelization for the irregular remainder. That split is the point. It gives you a +low-dimensional, physically interpretable description of the bulge *and* a flux map of exactly +what the bulge could not explain, which is the quantity you want if the science question is "how +much of this galaxy's light is in its irregular components?" + +This is a **direct** reconstruction. There is no second plane, no de-projection, no inversion of +a geometric mapping — the mesh sits on the image plane, over the galaxy, at the same coordinates +the data occupies. If you have met pixelized reconstruction before in a different context, that +extra step is the one that is absent here, and its absence makes everything simpler: what the +mesh reconstructs is the galaxy's own surface brightness, and its coordinates are sky +coordinates in arcseconds. + +Statistically, the fit now has two nested layers. The non-linear search explores a handful of +parameters — the bulge's shape, the regularization strength — and for each proposal an inversion +solves the mesh amplitudes exactly and returns a **Bayesian evidence** rather than a plain +likelihood. That evidence includes terms that penalise an unnecessarily complex reconstruction, +which is what lets a model with 400 free flux values not simply win by over-fitting. The theory +is [`../wiki/core/concepts/inversions_and_pixelizations.md`](../wiki/core/concepts/inversions_and_pixelizations.md); +chapter 4 of HowToGalaxy derives it from scratch. + +The canonical scripts are +`autogalaxy_workspace:scripts/imaging/features/pixelization/modeling.py` (via a search), +`autogalaxy_workspace:scripts/imaging/features/pixelization/fit.py` (one direct fit, with the +inversion internals walked through) and +`autogalaxy_workspace:scripts/imaging/features/pixelization/galaxy_reconstruction.py` (export). + +## Before you reach for one — the honest gate + +A pixelization is the most powerful and the most expensive model here, with failure modes the +parametric models do not have. Three questions decide whether you need it: + +- *"What do the residuals actually look like?"* Structured residuals that **follow the + isophotes** — a twist, a four-lobed pattern, an ellipticity that changes with radius — are a + *shape* problem, and an MGE fixes them for four free parameters + ([`ag_basis_profiles`](./ag_basis_profiles.md)). Residuals in **discrete off-centre lumps** are + a substructure problem, and that is what a pixelization is for. Look before you choose. +- *"Is the structure resolved?"* At low resolution the irregular structure is not in the data, + and a pixelization will reconstruct noise while telling you it fitted well. The workspace is + blunt about this: for low-resolution data pixelizations are unnecessary. +- *"Do you need global quantities or the light distribution?"* If you want a total flux, a size, + an axis ratio, a Sersic index — a parametric profile, an MGE or shapelets are enough and far + cheaper. If you want the morphology itself, particularly for faint or low-surface-brightness + features, nothing else comes close. + +If the answer to the first is "lumps" and to the second "yes", continue. Otherwise you will +spend hours to learn less. + +## Ask + +- *"What is the smooth component?"* Almost always a linear `Sersic`, sometimes a Sersic plus an + Exponential disk. It matters that this exists: without it the pixelization has to reconstruct + the bright central bulge too, which wastes mesh resolution on light a six-parameter profile + describes perfectly and makes the regularization choice much harder. +- *"How many mesh pixels can you afford?"* This is a hardware question as much as a science one + (see the run-time branch). 20 × 20 is the workspace's default for a GPU run; 28 × 28 to 30 × 30 + for a better-resolved reconstruction; higher only with a large-VRAM GPU or a many-core CPU. +- *"Are there neighbours inside the mask?"* For a pixelization the answer is *never* to mask + their pixels out — that punches holes in the mesh and produces discontinuity systematics. + Scale their noise instead (branch below). +- *"GPU or CPU?"* Genuinely depends on your pixel scale, and the answer is counter-intuitive. + Ask before assuming. + +## Branch — the hybrid model + +The canonical composition. Adapted from +`autogalaxy_workspace:scripts/imaging/features/pixelization/modeling.py`. + +```python +""" +Galaxy Structure: Bulge Plus Pixelized Reconstruction +===================================================== + +Fit a galaxy whose light has two very different characters: a smooth central bulge, described by +a linear Sersic profile in six parameters, and asymmetric clumpy star formation, reconstructed +on a regularized rectangular mesh. The Sersic captures the smooth component interpretably; the +pixelization absorbs everything it cannot explain, and the reconstructed flux map is the +measurement of the galaxy's irregular light. + +__Contents__ + +- **Imports:** Import the required libraries. +- **Dataset:** Load the imaging, mask it, and apply pixelization over-sampling. +- **Mesh Shape:** Fix the mesh dimensions before the model is composed. +- **Model:** Compose the linear Sersic bulge plus the pixelization. +- **Check:** Confirm the parameter count, evaluate one likelihood, and estimate VRAM. +""" + +""" +__Imports__ +""" +from pathlib import Path + +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +DATASET_PATH = Path("dataset") / "imaging" / "my_galaxy" +PIXEL_SCALES = 0.1 +MASK_RADIUS = 2.0 + +""" +__Dataset__ + +The mask matters more than usual here. It defines the region the mesh is laid over, so a mask +that reaches well beyond the galaxy spends reconstruction pixels on empty sky, and one that +truncates the galaxy leaves its outer light unreconstructed +(`PyAutoArray:autoarray/dataset/imaging/dataset.py`). + +The over-sampling argument changes too. A pixelization evaluates light on its *own* grid with +its own over-sampling scheme, so the size is passed as `over_sample_size_pixelization`, not the +`over_sample_size_lp` used for light profiles. Passing the wrong one silently leaves the +pixelization grid at its default. +""" +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_pixelization=over_sample_size) + +""" +__Mesh Shape__ + +The mesh dimensions are **fixed before the fit and cannot be a free parameter**. This is not a +stylistic restriction: JAX compiles the likelihood against statically shaped arrays, and the +mesh shape sets those shapes. A rectangular mesh uses the same number of pixels in y and x. +""" +mesh_pixels_yx = 20 +mesh_shape = (mesh_pixels_yx, mesh_pixels_yx) + +""" +__Model__ + +`ag.Pixelization` pairs a `mesh` — where the reconstruction pixels sit — with a `regularization` +scheme — the prior that keeps the solution smooth. It is attached to the galaxy as a component +alongside the light profiles, and the analysis knows to route it through the inversion +(`PyAutoGalaxy:autogalaxy/galaxy/galaxy.py`). + +The mesh itself contributes **zero** free parameters; the reconstruction's 400 flux values are +solved, not sampled. `ag.reg.GaussianKernel` contributes two (`coefficient`, `scale`). The bulge +is `ag.lp_linear.Sersic` rather than `ag.lp.Sersic` deliberately: its `intensity` is then solved +by the *same* inversion that solves the mesh amplitudes, which removes the degeneracy between +"bright bulge, faint reconstruction" and "faint bulge, bright reconstruction" that a sampled +intensity would create. +""" +pixelization = af.Model( + ag.Pixelization, + mesh=ag.mesh.RectangularAdaptDensity(shape=mesh_shape), + regularization=ag.reg.GaussianKernel, +) + +galaxy = af.Model( + ag.Galaxy, + redshift=0.5, + bulge=ag.lp_linear.Sersic, + pixelization=pixelization, +) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +""" +__Check__ + +Three cheap checks, each catching a different class of mistake: the parameter count catches a +model error, one likelihood evaluation catches an incompatibility between model, dataset and +analysis, and `print_vram_use` catches the resource error that otherwise appears as an +out-of-memory failure mid-compile. Pixelizations use far more VRAM than profile-only models — +around 0.05 GB per batched likelihood for 400 reconstruction pixels — which is why the +workspace uses a lower `n_batch` for these fits than for any other. +""" +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 evidence at prior medians: {float(log_likelihood):.2f}") +``` + +**N = 8**: the Sersic's centre (2), `ell_comps` (2), `effective_radius`, `sersic_index`, plus the +regularization's `coefficient` and `scale`. Eight parameters, with 400 reconstruction pixels +solved. That is dramatically more parsimonious than stacking parametric profiles until every +clump is described — a 20-plus-parameter approach that would still fail on genuinely irregular +substructure. + +Extending the smooth side is exactly as you would expect — add components as further attributes +on the galaxy, keeping them linear so their intensities join the same inversion: + +```python +galaxy = af.Model( + ag.Galaxy, + redshift=0.5, + bulge=ag.lp_linear.Sersic, + disk=ag.lp_linear.Exponential, + pixelization=pixelization, +) +``` + +## Branch — choosing a mesh and a regularization scheme + +These are the two decisions that make or break a reconstruction, and the parameter cost is +almost entirely in the regularization. + +**Meshes** (`ag.mesh`, from `PyAutoArray:autoarray/inversion/mesh/`). All contribute zero free +parameters, all take `shape=(y, x)`: + +| Mesh | Where the pixels go | Use when | +|---|---|---| +| `RectangularUniform` | a uniform grid over the masked region | you want the simplest, most predictable behaviour | +| `RectangularAdaptDensity` | a rectangular grid whose spacing follows the density of image pixels mapping to it | **the default** — concentrates resolution where the data is | +| `RectangularAdaptImage` | adapts to the image's own brightness distribution | the irregular light is strongly concentrated | + +`Delaunay`, `KNearestNeighbor` and `KNNBarycentric` also exist for irregular tessellations. Start +with `RectangularAdaptDensity`; it is what every workspace example uses, and the rectangular +meshes are what the fixed-`mesh_shape` JAX path is built around. + +**Regularization** (`ag.reg`, from `PyAutoArray:autoarray/inversion/regularization/`) is where +the free parameters are, and where the physics of "how smooth should this be?" lives: + +| Scheme | Free params | What it assumes | +|---|---|---| +| `ag.reg.Constant` | 1 (`coefficient`) | one smoothing strength everywhere — the simplest, and the right first try | +| `ag.reg.GaussianKernel` | 2 (`coefficient`, `scale`) | smoothing with a Gaussian correlation length you also fit | +| `ag.reg.MaternKernel` | 3 (`coefficient`, `scale`, `nu`) | as above, with the roughness `nu` free too | + +Verified counts, not guesses. `Constant` penalises the difference between neighbouring pixels +directly; the kernel schemes impose a correlation *length*, which is physically closer to what a +galaxy's light does and generally reconstructs faint extended structure better at the cost of +one or two parameters. `ag.reg.Adapt` and its variants scale the smoothing with the +reconstruction's own brightness — more flexible again, and requiring `adapt_images` on the +analysis, which is a step beyond this skill; ask if you want it. + +A practical route: fit once with `ag.reg.Constant`, look at the reconstruction. Over-smoothed +(the clumps blurred into one blob) or over-fitted (isolated single-pixel spikes) both show up +immediately, and if `Constant` cannot get both the bright clumps and the faint envelope right at +one strength, that is the signal to move to `GaussianKernel`. Full discussion in +[`../wiki/core/concepts/inversions_and_pixelizations.md`](../wiki/core/concepts/inversions_and_pixelizations.md). + +### The positive-only solver + +Every pixelized reconstruction uses a **positive-only** linear solver: no mesh pixel may +reconstruct negative flux. This is not a numerical nicety. A signed solver lets the +reconstruction manufacture negative pixels that cancel positive ones to over-fit the data, and +that failure mode is a well-known systematic in the literature — many published methods allow +it. Enforcing non-negativity efficiently needs bespoke linear algebra, which is why this is +worth stating rather than assuming. You do not need to switch it on; it is the behaviour you get. +(The one place in this workspace where you *must* disable it is a shapelet basis, which needs +signed amplitudes by construction — see [`ag_basis_profiles`](./ag_basis_profiles.md).) + +## Branch — contaminating neighbours: scale, never mask + +The general advice for a neighbour inside the mask is to mask its pixels out. **For a +pixelization that advice is wrong.** Removing image pixels entirely removes their mapping into +the mesh, punching a hole in the reconstruction and producing discontinuities that show up as +unexplained systematics. + +Instead, keep the pixels in the fit and make them contribute nothing: scale their data to zero +and their noise to an enormous value, so the likelihood is indifferent to them while the mesh +stays continuous. Adapted from +`autogalaxy_workspace:scripts/imaging/features/pixelization/fit.py`. + +```python +""" +__Noise Scaling__ + +`mask_extra_galaxies.fits` marks the contaminated pixels. It is loaded with `invert=True` +because in this file `True` means "scale this pixel", the opposite of a modelling mask's +convention — getting this backwards scales the galaxy instead of the neighbour, so check the +subplot afterwards. `apply_noise_scaling` then zeroes those data values and inflates their +noise, driving their signal-to-noise to effectively zero without removing them from the mesh's +mapping (`PyAutoArray:autoarray/dataset/imaging/dataset.py`). + +Order matters: scale first, then apply the modelling mask. +""" +mask_extra_galaxies = ag.Mask2D.from_fits( + file_path=DATASET_PATH / "mask_extra_galaxies.fits", + pixel_scales=PIXEL_SCALES, + invert=True, +) + +dataset = dataset.apply_noise_scaling(mask=mask_extra_galaxies) + +mask = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=PIXEL_SCALES, + centre=(0.0, 0.0), + radius=MASK_RADIUS, +) +dataset = dataset.apply_mask(mask=mask) + +aplt.subplot_imaging_dataset( + dataset=dataset, + output_path="scripts/scratch/my_galaxy/", + output_filename="noise_scaled", + output_format="png", +) +``` + +The scaled pixels go to a noise value of order 1e8 — visibly so on the signal-to-noise panel of +the subplot, which is how you confirm you scaled the right region. Drawing that mask +interactively is `autogalaxy_workspace:scripts/imaging/data_preparation/gui/mask_extra_galaxies.py`; +the modelling alternative, and when to prefer it, is +[`ag_light_model_extras`](./ag_light_model_extras.md) and +[`../wiki/core/concepts/extra_galaxies_and_noise_scaling.md`](../wiki/core/concepts/extra_galaxies_and_noise_scaling.md). + +## Branch — reading the reconstruction out + +The reconstruction is the science product, and it lives on the `Inversion`. You can get there +from a completed fit or from a single direct `ag.FitImaging` — the latter is the fast way to +learn the objects, and is what +`autogalaxy_workspace:scripts/imaging/features/pixelization/fit.py` does. + +```python +""" +__Reconstruction__ + +`inversion.linear_obj_list` holds the linear objects the fit solved. For a hybrid model that is +one `LightProfileLinearObjFuncList` per linear light profile plus one `Mapper` for the +pixelization; with a non-linear `ag.lp.Sersic` bulge the `Mapper` is the only entry. The +`Mapper` is the object that maps image pixels onto mesh pixels, and it is mesh-specific +(`PyAutoArray:autoarray/inversion/mappers/`). + +The reconstruction is a 1D array of one flux per mesh pixel, and the (y, x) arcsecond +coordinates of those pixels are `mapper.mesh_geometry.mesh_grid`. Pair them by index — they are +the same length and the same order — and you have the flux map. +""" +inversion = fit.inversion + +mapper = inversion.linear_obj_list[0] + +reconstruction = inversion.reconstruction +mesh_grid = mapper.mesh_geometry.mesh_grid + +print(f"{len(reconstruction)} mesh pixels reconstructed") +print(f"coordinates shape: {mesh_grid.shape}") + +""" +The reconstruction's own uncertainties, and its projection back onto the image, are also on the +inversion. `mapped_reconstructed_operated_data` contains **only** the pixelized component — any +parametric light profile in the model is excluded — which is exactly the quantity you want when +asking how much light is in the irregular component. +""" +print(inversion.reconstruction_noise_map) +print(inversion.mapped_reconstructed_operated_data.native.shape) +``` + +A caution on the `Mapper`'s grid attributes: `mapper.image_plane_mesh_grid` is `None` for a +rectangular mesh, so do not reach for it. `mapper.mesh_geometry.mesh_grid` is the mesh pixel +centres (verified on the released stack), and `mapper.image_plane_data_grid` is the centre of +every masked *image* pixel. Some of the `Mapper`'s attribute names carry a plane-based prefix +inherited from the shared inversion library that also serves a multi-plane use case; in +PyAutoGalaxy there is only the one plane, so `mesh_geometry.mesh_grid` is both the clearer name +and the one to use. + +### The evidence terms + +The inversion's figure of merit is a Bayesian evidence, and its components are individually +readable — which is how you diagnose whether regularization is doing too much or too little: + +```python +print(f"regularization term = {inversion.regularization_term}") +print(f"log det regularization = {inversion.log_det_regularization_matrix_term}") +print(f"log det curvature + reg = {inversion.log_det_curvature_reg_matrix_term}") +``` + +The `regularization_term` measures how non-smooth the solution is; the two log-determinant +terms are the complexity penalty that stops a 400-pixel model winning by over-fitting. The +matrices themselves (`curvature_matrix`, `regularization_matrix`, `curvature_reg_matrix`) are +there too. The derivation is in HowToGalaxy chapter 4 and the papers it cites +(arXiv:1708.07377, arXiv:astro-ph/0601493). + +### Diagnostic plots + +`aplt.subplot_fit_imaging` includes the reconstruction on the mesh alongside the usual data / +model / residual panels, and is the first thing to look at. For the inversion's own internals +there are dedicated functions, which live in the shared inversion library rather than on `aplt` +and so are imported explicitly — from +`autogalaxy_workspace:scripts/imaging/features/pixelization/plot.py`: + +```python +from autoarray.inversion.plot.inversion_plots import subplot_of_mapper +from autoarray.inversion.plot.mapper_plots import plot_mapper, subplot_image_and_mapper + +PLOT_DIR = "scripts/scratch/my_galaxy/pixelization/" + +aplt.subplot_fit_imaging(fit=fit, output_path=PLOT_DIR, output_format="png") + +subplot_of_mapper( + inversion=inversion, mapper_index=0, output_path=PLOT_DIR, output_format="png" +) + +subplot_image_and_mapper( + mapper=mapper, image=dataset.data, output_path=PLOT_DIR, output_format="png" +) + +print(f"Saved to: {Path(PLOT_DIR).resolve()}") +``` + +`subplot_of_mapper` is the comprehensive diagnostic — reconstructed image, reconstruction, its +noise map and the regularization weights in one figure — and `mapper_index` selects which linear +object when there are several. `subplot_image_and_mapper` shows the data beside the mesh, which +is how you check the mesh covers the structure you care about. Note that none of these three are +exposed on the `aplt` module — `dir(aplt)` does not list them, so the explicit imports above are +the only route and writing them as `aplt.`-prefixed calls raises `AttributeError`. Quote the +absolute path back to the user and offer to open it. + +## Branch — exporting the reconstruction for downstream science + +A reconstruction that only exists inside a PyAutoGalaxy object is hard to share and hard to +analyse with other tools. Write it out as coordinates plus values and it becomes ordinary +tabular data. Adapted from +`autogalaxy_workspace:scripts/imaging/features/pixelization/galaxy_reconstruction.py`. + +```python +""" +__Export__ + +Four columns — y, x, reconstruction, noise — is everything needed to rebuild, plot or analyse +the galaxy's irregular light without PyAutoGalaxy installed. A completed model-fit also writes +this file into its own `image/` folder automatically; glob for `*reconstruction_0.csv` there +rather than hard-coding the stem, which carries a prefix from the shared inversion library. +""" +import csv + +CSV_PATH = Path("scripts/scratch/my_galaxy/reconstruction.csv") +CSV_PATH.parent.mkdir(parents=True, exist_ok=True) + +with CSV_PATH.open("w", newline="") as f: + writer = csv.writer(f) + writer.writerow(["y", "x", "reconstruction", "noise_map"]) + for (y, x), value, sigma in zip( + mesh_grid, inversion.reconstruction, inversion.reconstruction_noise_map + ): + writer.writerow([y, x, value, sigma]) + +print(f"Reconstruction written to: {CSV_PATH.resolve()}") +``` + +Reloading needs nothing but the standard library, and from there the mesh is ordinary scattered +data. Two things the workspace does with it: a Delaunay triangulation via `scipy.spatial` for +plotting on the mesh's own geometry, and interpolation onto a regular grid via +`scipy.interpolate.griddata` for any tool that wants an image: + +```python +import numpy as np +from scipy.interpolate import griddata + +with CSV_PATH.open() as f: + reader = csv.reader(f) + header = next(reader) # ['y', 'x', 'reconstruction', 'noise_map'] + loaded = {key: [] for key in header} + for row in reader: + for key, value in zip(header, row): + loaded[key].append(float(value)) + +points = np.stack((np.array(loaded["x"]), np.array(loaded["y"])), axis=-1) + +interpolation_grid = ag.Grid2D.from_extent( + extent=(-1.0, 1.0, -1.0, 1.0), shape_native=(201, 201) +) + +interpolated = griddata( + points=points, values=np.array(loaded["reconstruction"]), xi=interpolation_grid +) +``` + +Pick the `extent` to cover the reconstruction without over-resolving it — the mesh has a few +hundred pixels, so interpolating to 201 × 201 is presentation, not information. + +## Branch — run time, VRAM, and the GPU-versus-CPU choice + +This is the branch that decides whether your fit takes twenty minutes or three hours, and the +answer is not "use the GPU". + +Pixelized inversions do linear algebra on very large, very **sparse** matrices. JAX has no sparse +support and must work densely, which scales badly; the CPU implementation exploits the sparsity +fully via `numba`. So the trade runs against pixel scale: + +- **Low resolution (`pixel_scales > 0.05`, e.g. Euclid)** — fewer sparse operations, modest VRAM. + **GPU with JAX is usually fastest.** +- **High resolution (`pixel_scales <= 0.03`, e.g. HST or JWST)** — the linear algebra is + sparsity-dominated. **CPU with `numba` and many cores can beat a powerful GPU.** + +VRAM is the hard constraint on the GPU side. JAX must hold the whole batched likelihood in GPU +memory, and a pixelization needs far more than a profile model: around 0.05 GB per batched +likelihood at 400 reconstruction pixels, and more than 1 GB — occasionally more than 10 GB — at +high resolution with a fine mesh. Two consequences: keep `n_batch` low (the workspace uses 20 for +pixelizations against 50 elsewhere), and run + +```python +analysis.print_vram_use(model=model, batch_size=search.batch_size) +``` + +before committing. It takes 20–30 seconds and it is the difference between knowing your limit +and discovering it as an out-of-memory error mid-compile. If VRAM is exceeded the run time does +not degrade gracefully — it goes from under ten minutes to three hours or more. + +Rough scaling with resolution, from the workspace's own measurements at a 20 × 20 mesh with VRAM +under control: ~10 minutes at 0.1"/pixel, ~30 minutes at 0.05", ~1 hour at 0.03". Benchmark both +paths for your own data and hardware rather than trusting either default; the search-side +settings are [`ag_configure_search`](./ag_configure_search.md) and +[`../wiki/core/api/searches.md`](../wiki/core/api/searches.md). + +## Combine — where this hands off + +- **Try a basis first** → [`ag_basis_profiles`](./ag_basis_profiles.md). An MGE is four + parameters and minutes, not hours, and it genuinely resolves the shape-mismatch case a + pixelization is often reached for by mistake. +- **Pick and configure the search** → [`ag_configure_search`](./ag_configure_search.md), with + `n_batch` lowered for VRAM. +- **Run the fit** → [`ag_run_search`](./ag_run_search.md); `image/fit.png` in the output folder + refreshes on the fly and shows the reconstruction as it converges, which is worth watching for + a fit this long. +- **Plot and inspect** → [`ag_plot_fit`](./ag_plot_fit.md), plus the `subplot_of_mapper` + diagnostics above. +- **Load a completed fit programmatically** → [`ag_load_results`](./ag_load_results.md). +- **Over-smoothed, over-fitted or unphysical reconstruction** → + [`ag_debug_fit_failure`](./ag_debug_fit_failure.md). The pixelization-specific causes are a + regularization scheme too rigid for the structure, a mask that punched holes via hard masking + where noise scaling was needed, over-sampling passed as `over_sample_size_lp` instead of + `over_sample_size_pixelization`, and VRAM exhaustion masquerading as a hung fit. +- **Contaminating neighbours modelled rather than scaled** → + [`ag_light_model_extras`](./ag_light_model_extras.md). +- **Make the fit robust and fast by starting simple** → the search-chaining skill + (`ag_chain_searches`). This is the biggest available win for a pixelization: fit a parametric + model first, pass its bulge as priors, and let the expensive model start from a sensible place. + The workspace recommends it explicitly. +- **A non-parametric shape measurement instead of a reconstruction** → + [`ag_ellipse_fitting`](./ag_ellipse_fitting.md). + +Offer (default-yes) a dated `wiki/project/YYYY-MM-DD-.md` entry recording the mesh and its +shape, the regularization scheme and why, the mask radius, the hardware path, and what the +reconstruction showed — with a pixelization the choices *are* the result, and none of them are +recoverable from the output folder alone. + +## Further reading + +- **Student / new to galaxy morphology** — [HowToGalaxy: Bayesian regularization](https://github.com/PyAutoLabs/HowToGalaxy/blob/main/notebooks/chapter_4_pixelizations/tutorial_4_bayesian_regularization.ipynb): + why a model with hundreds of free flux values does not simply over-fit, derived from scratch; + `tutorial_5_model_fit` then runs the full fit. +- **General reference** — [RTD: Pixelization API](https://pyautogalaxy.readthedocs.io/en/latest/api/pixelization.html): + every mesh, regularization scheme and the objects that combine them. +- **Experienced PyAutoGalaxy user** — [workspace: imaging/features/pixelization/modeling.py](https://github.com/PyAutoLabs/autogalaxy_workspace/blob/main/scripts/imaging/features/pixelization/modeling.py): + the hybrid fit end to end, with the run-time, VRAM and GPU-versus-CPU discussion in full. + +## Agent procedural checklist + +1. Look at the residual map of a parametric or basis fit first. Shape-following residuals mean + `ag_basis_profiles`; discrete off-centre lumps mean a pixelization. Do not skip this. +2. Confirm the structure is resolved — at coarse pixel scales a pixelization reconstructs noise. +3. Always pair the pixelization with a linear light profile for the smooth component; never let + the mesh reconstruct the bulge. +4. Fix `mesh_shape` before composing the model; it cannot be a free parameter. +5. Pass over-sampling as `over_sample_size_pixelization`, not `over_sample_size_lp`. +6. Start with `RectangularAdaptDensity` + `ag.reg.Constant`; move to `GaussianKernel` only when + one smoothing strength demonstrably cannot serve both the bright and faint structure. +7. Scale contaminant noise (`apply_noise_scaling`, `invert=True`) — never hard-mask pixels inside + a pixelization fit — and check the signal-to-noise panel afterwards. +8. `print(model.info)`, evaluate one likelihood, then `print_vram_use` before any GPU run; keep + `n_batch` low. +9. Read the reconstruction via `inversion.reconstruction` + + `mapper.mesh_geometry.mesh_grid`; export to CSV for anything downstream. +10. Save the script to `scripts/`, quote every plot's absolute path and offer to open it, then + offer the `wiki/project/` entry recording mesh, regularization and hardware choices. diff --git a/wiki/core/api/configuration.md b/wiki/core/api/configuration.md index c497a52..9ab66a9 100644 --- a/wiki/core/api/configuration.md +++ b/wiki/core/api/configuration.md @@ -41,7 +41,7 @@ sources: - config/visualize/plots_search.yaml pinned_commit: ed72fabb33e14a9a701a4d280e8775dd3a20e98c last_updated: 2026-08-01 -content_sha256: bdc7d14e2ec2fc6d73b703baa17299f0bb00f69a4dcdb3410fe68f14c564a8d8 +content_sha256: 6d9292c60dac5970b1d0a7960888513daa5834b27ca9b08bb325966b671841b4 --- # Configuration @@ -177,10 +177,17 @@ Three keys per parameter: the galaxy at the origin, and `ell_comps` are `TruncatedGaussian` clipped to ±1 because that is their valid range. - **`width_modifier`** — how wide the passed prior becomes when one search's result seeds the next. - `Absolute` adds a fixed width, `Relative` a fraction of the inferred value. This is the knob that - governs search chaining (`autogalaxy_workspace:scripts/guides/modeling/chaining.py`). -- **`limits`** — hard physical bounds the passed `GaussianPrior` may not cross. `sersic_index` is - capped at `[0.8, 5.0]`, so chaining can never propose an unphysical index. + `Absolute` sets `sigma` to the value directly, `Relative` to that fraction of the inferred value. + This is the knob that governs search chaining, and the attribute that reads it is + **`result.model_centred`** — which replaces every free parameter with a + `TruncatedGaussianPrior` centred on its previous median. `result.model` does **not** consult it: + on the released stack that attribute returns the model with its original priors unchanged, so a + component passed through it inherits nothing from the fit. See + [`../../../skills/ag_chain_searches.md`](../../../skills/ag_chain_searches.md) for the full + three-way split between `model`, `model_centred` and `instance`, and note that the prose in + `autogalaxy_workspace:scripts/guides/modeling/chaining.py` is out of date on this point. +- **`limits`** — hard physical bounds the passed `TruncatedGaussianPrior` may not cross. + `sersic_index` is capped at `[0.8, 5.0]`, so chaining can never propose an unphysical index. `PyAutoFit:autofit/config/priors/` holds the meta-templates the per-class files are shaped from. diff --git a/wiki/core/api/searches.md b/wiki/core/api/searches.md index 7c0e638..05c9387 100644 --- a/wiki/core/api/searches.md +++ b/wiki/core/api/searches.md @@ -30,7 +30,7 @@ sources: - config/general.yaml pinned_commit: ed72fabb33e14a9a701a4d280e8775dd3a20e98c last_updated: 2026-08-01 -content_sha256: ce9e29263a3ee8643317388d2375ef4866a7fb561a0964b5f5a1211b07311dc9 +content_sha256: bfae75e9e59dfa77b2477ef56426828e889e8df2fa13111cc00bc9c4e015a9de --- # Non-linear search catalogue @@ -322,9 +322,11 @@ unit-prior values). Source: `PyAutoFit:autofit/non_linear/grid/`. A different lever from choosing a better search: run a cheap search first, then use its posterior as the priors of a more complex one. This is how you fit a model too complex for a single search — fit a single Sersic, then pass the result forward to a bulge-plus-disk or pixelised fit. -`autogalaxy_workspace:scripts/guides/modeling/chaining.py` is the walkthrough, and the -`width_modifier` entries in the prior YAMLs are what control how wide the passed priors become -(see [`configuration`](./configuration.md)). +`autogalaxy_workspace:scripts/guides/modeling/chaining.py` is the walkthrough, the attribute that +narrows is `result.model_centred` (**not** `result.model`, which returns the original priors +unchanged), and the `width_modifier` entries in the prior YAMLs are what control how wide the +passed priors become (see [`configuration`](./configuration.md) and +[`../../../skills/ag_chain_searches.md`](../../../skills/ag_chain_searches.md)). ## Picking a search at a glance diff --git a/wiki/core/concepts/non_linear_search.md b/wiki/core/concepts/non_linear_search.md index 852eb11..c348c4a 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: 13439cc6b908483c12c5b907e921b9491bd3ee179fc328654c130e2ac1048fce +content_sha256: 05120accc325720145ce1d7fbf2decebb9f93bb5df5ff617fe842ebb7a1c5010 --- # The non-linear search @@ -191,22 +191,37 @@ not. Use tight priors to encode real knowledge, not to steer a search. ## Search chaining -A search can be initialised from the posterior of a previous one. `result.model` returns a new -model whose priors are `TruncatedGaussianPrior`s centred on the previous median-PDF values, -with widths taken from each parameter's `width_modifier` in the priors configuration: +A search can be initialised from the posterior of a previous one. The attribute that does the +narrowing is **`result.model_centred`**: it returns a new model whose priors are +`TruncatedGaussianPrior`s centred on the previous median-PDF values, with widths taken from each +parameter's `width_modifier` in the priors configuration: ```python -bulge = result_1.model.galaxies.galaxy.bulge -disk = result_1.model.galaxies.galaxy.disk +bulge = result_1.model_centred.galaxies.galaxy.bulge +disk = result_1.model_centred.galaxies.galaxy.disk galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, disk=disk) model_2 = af.Collection(galaxies=af.Collection(galaxy=galaxy)) ``` -`autogalaxy_workspace:scripts/guides/modeling/chaining.py`. Whole profiles or whole galaxies -can be passed this way. The alternative, `result_1.instance...`, passes the maximum-likelihood -values as **fixed** numbers, removing them from the search — useful for fixing a light model -in an intermediate stage and freeing it again later. +Three variants override that width: `model_centred_absolute(a=)` (one absolute `sigma` for every +parameter), `model_centred_relative(r=)` (`sigma = r × mean`) and +`model_centred_max_lh_bounded(b=)` (a `UniformPrior` at `mean ± b`). + +**`result.model` does not narrow anything.** On the released stack it returns the fitted model +with its **original priors unchanged** — `samples_summary.model.mapper_via_defaults_from()` maps +every prior to itself (`PyAutoFit:autofit/non_linear/result.py`) — so a component passed that way +brings its composition across but no information from the fit. Reach for it only when that is +what you want; the prose in +`autogalaxy_workspace:scripts/guides/modeling/chaining.py` still describes `result.model` as +producing narrowed Gaussians and is out of date on this point. +[`../../../skills/ag_chain_searches.md`](../../../skills/ag_chain_searches.md) is the procedural +recipe and carries the full three-way split. + +Whole profiles or whole galaxies can be passed either way. The third alternative, +`result_1.instance...`, passes the maximum-likelihood values as **fixed** numbers, removing them +from the search — useful for fixing a light model in an intermediate stage and freeing it again +later. Chaining is what makes hard fits tractable: start with a simple parametric model, then chain into a pixelisation or a many-component basis with the simple fit's posterior as the starting diff --git a/wiki/core/external/skill_citation_map.md b/wiki/core/external/skill_citation_map.md index 6ed20eb..633fd07 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: 1d29ca6a67048d5e548dd83006597944a348f31ade73b852b625800a64d1f40b +content_sha256: 64f9a361cee2627eb4f7fc2be4f77eb06bc8a7e80fb4beb29a981b01ff08a701 --- # Skill → external resource map @@ -62,6 +62,14 @@ Default to `.ipynb` for student-leaning users and `.py` for returning PyAutoGala | `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` | +| `ag_basis_profiles` | `chapter_2_modeling/tutorial_5_linear_profiles` | `api/light` | `imaging/features/multi_gaussian_expansion/modeling.py` | +| `ag_pixelization` | `chapter_4_pixelizations/tutorial_4_bayesian_regularization` | `api/pixelization` | `imaging/features/pixelization/modeling.py` | +| `ag_light_model_extras` | `chapter_3_search_chaining/tutorial_3_x2_galaxies` | `overview/overview_3_features` | `imaging/features/extra_galaxies/modeling.py` | +| `ag_ellipse_fitting` | _ | `overview/overview_3_features` | `ellipse/modeling.py` | +| `ag_multi_dataset` | `chapter_2_modeling/tutorial_5_linear_profiles` | `overview/overview_3_features` | `multi_dataset/start_here.py` | +| `ag_build_interferometer_model` | `chapter_1_introduction/tutorial_3_fitting` | `overview/overview_3_features` | `interferometer/start_here.py` | +| `ag_multi_galaxy_and_cluster` | `chapter_3_search_chaining/tutorial_3_x2_galaxies` | `overview/overview_2_new_user_guide` | `cluster/start_here.py` | +| `ag_chain_searches` | `chapter_3_search_chaining/tutorial_1_search_chaining` | `general/configs` | `guides/modeling/chaining.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 @@ -77,9 +85,21 @@ before it was recorded. A skill may cite an *extra* tutorial inline within a bul 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. +The eight feature rows after them follow the same rule, with two honest gaps worth reading +before you "fix" a cell that looks wrong: + +- **`ag_ellipse_fitting`'s student cell is `_`.** The lecture series teaches light-profile + modelling and has no ellipse-fitting chapter, so the skill omits that bullet and routes a + newcomer to [`../concepts/ellipse_fitting_and_multipoles.md`](../concepts/ellipse_fitting_and_multipoles.md) + instead. A `_` here is the correct state, not a missing lookup. +- **Three rows cite a tutorial that is not about their own subject.** There is no + multi-wavelength chapter and no interferometer chapter, so `ag_multi_dataset` and + `ag_build_interferometer_model` cite the tutorial that teaches the idea each fit leans on + hardest — linear profiles, and the likelihood — and `ag_light_model_extras` shares + `tutorial_3_x2_galaxies` with `ag_multi_galaxy_and_cluster` because a second galaxy in the + frame is where the lectures come closest to a contaminant. Each of those bullets says so in + the skill, rather than implying a chapter exists that does not. Sharing a cell across two + rows is allowed; inventing one is not. ## Template for the inserted skill block