audit: Delaunay pixel areas and magnification paths (euclid-dr1-prep phase 8)

[Jammy2211]

Overview

Phase 8 of the euclid-dr1-prep epic (was 6c; renumbered in the 2026-09-01 Cortex split). The Delaunay source model's magnification has never been validated and the epic's owner suspects the pixel areas. This is a source-code audit of every Delaunay area/magnification path in PyAutoArray, with known-answer tests landed regardless of outcome. It is the software half of the question; the empirical half is Cortex phase 7 (PyAutoCortex/phases/euclid/magnification_robustness.md), which this audit must brief. Any real defect is filed as a separate bug prompt — this task does not grow into the fix.

Plan

• Map every Delaunay area/magnification code path — Voronoi cell areas behind areas_for_magnification, barycentric dual areas behind the split-point weights, the tripcolor plotting path, the workspace source_science.py consumer, and the Euclid pipeline's magnification latent — and state what each computes.
• Prove or drop each suspicion with known-answer tests (analytic lattices, convex-hull totals, exact integrals of a linear field), not resemblance.
• Land the missing direct tests for areas_for_magnification and barycentric_dual_area_from, and correct the areas_for_magnification docstring to what is actually zeroed.
• Check whether the Euclid pipeline's magnification latent survives a Delaunay vis_pix / full_model stage at all.
• Write the verdict here and in the epic ledger, with an explicit statement for Cortex phase 7: are Delaunay magnification numbers trustworthy, and which area must they use?
• Any proven defect → separate draft/bug/autoarray/ prompt with the reproduced failing case.

Detailed implementation plan

Affected Repositories

• PyAutoArray (primary — tests + docstring)
• PyAutoLens (read-only: autolens/analysis/latent.py::magnification reachability check; no edits unless a defect is found, then a separate prompt)

Branch Survey

Repository	Current Branch	Dirty?
./PyAutoArray	main	clean
./PyAutoLens	main	clean

Suggested branch: feature/delaunay-area-magnification-audit

What the pre-plan survey established

1. Two different "areas" exist for a Delaunay mesh on two different paths:
   • MeshGeometryDelaunay.areas_for_magnification (autoarray/inversion/mesh/mesh_geometry/delaunay.py:195) = voronoi_areas_numpy (shoelace over scipy Voronoi cells). Only unbounded cells are marked -1 and zeroed; bounded-but-huge boundary cells are kept (the existing test__voronoi_areas_via_delaunay_from shows index 3 = 29.8 kept, index 4 = −1). The docstring says boundary pixels are zeroed — it overstates.
   • barycentric_dual_area_from (autoarray/inversion/mesh/interpolator/delaunay.py:342, plus an in-graph JAX scatter-add copy at ~:309) = Σ(triangle area / 3) per vertex; used only to weight split-regularisation points via areas_factor * sqrt(areas). Not used for flux or magnification.
2. Plotting (autoarray/plot/inversion.py::_plot_delaunay) uses tripcolor Gouraud — linear interpolation over triangles, whose integral is the barycentric dual area, not the Voronoi cell. Picture and number come from different geometry.
3. areas_for_magnification has zero library callers; it serves only the four pixelized source_science.py workspace scripts (Σ reconstruction × mesh_areas as source-plane flux). Those scripts are what Cortex phase 7's Delaunay leg would use.
4. The Euclid pipeline's magnification latent (autolens/analysis/latent.py::magnification) is a flux ratio via galaxies[-1].image_2d_from(...) and never touches areas_for_magnification. Whether it is finite for a Pixelization source galaxy is unverified.
5. sqrt is not on the magnification path; the known NaN-gradient hazard is confined to the interpolator's split-point weights.
6. The cluster-epic prompt draft/test/workspaces/mesh_magnification_correctness.md owns the library-level magnification API, non-Delaunay meshes and the simulate-and-recover campaign. This task stays Delaunay-only.

Implementation Steps

1. test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py — add:
   • test__areas_for_magnification__uniform_lattice: n×n unit lattice → every interior Voronoi cell area == 1.0; every hull point unbounded → 0.0; Σ == (n−2)². Pins the real semantics.
   • test__areas_for_magnification__bounded_boundary_cells_are_kept: the existing 6-point configuration → index 3 (≈29.8) is kept, only index 4 zeroed — documents the bias candidate explicitly.
   • test__areas_for_magnification__repeat_calls_agree: two calls agree; the -1 sentinel is not written back into any cached state.
2. test_autoarray/inversion/mesh/interpolator/test_delaunay.py (create if absent, mirroring the existing interpolator test layout) — add:
   • test__barycentric_dual_area__single_triangle: each vertex == A/3.
   • test__barycentric_dual_area__sums_to_convex_hull_area: random points → Σ dual == scipy.spatial.ConvexHull(points).volume; assert Σ areas_for_magnification differs for the same points (records the divergence).
   • test__linear_field_integral__dual_areas_exact: f = a + b·y + c·x at the vertices; exact integral over the hull equals Σ f_i·dual_i to rounding; Σ f_i·voronoi_i differs on a non-uniform mesh.
   • NumPy vs JAX parity: barycentric_dual_area_from(xp=np) equals the in-graph scatter-add on the same padded simplices.
3. autoarray/inversion/mesh/mesh_geometry/delaunay.py — rewrite the areas_for_magnification docstring to the proven semantics (unbounded cells → 0; bounded boundary cells kept; these are Voronoi cells, not the interpolant's dual areas). No behavioural change.
4. Audit evidence (scratch script, numbers reported here, not shipped): a Gaussian source strictly inside the hull of a Delaunay mesh; compare Σ reconstruction·voronoi areas, Σ reconstruction·dual areas, and the mapped-reconstruction image-plane flux against the analytic source flux. Report which denominator recovers μ and the sign/size of the boundary effect.
5. Euclid reachability (read-only): look for magnification in the latent outputs of the vis_pix / full_model stages under euclid_strong_lens_modeling_pipeline/output/ral_*; if absent (NaN latents are dropped from latent_summary.json), confirm by evaluating the latent on a stored Delaunay fit in test mode. Outcome is a finding.
6. If a defect is proven (expected candidate: the magnification denominator should integrate the interpolant with barycentric dual areas, or bounded boundary cells bias the Voronoi sum), file draft/bug/autoarray/<name>.md via /intake with the reproduced case and Epic: euclid-dr1-prep; cross-link from the ledger and here. Do not fix in this PR.
7. Write the per-path audit table (path, what it computes, verdict, evidence) and the Cortex-phase-7 statement into this issue and the epic ledger.

Key Files

• autoarray/inversion/mesh/mesh_geometry/delaunay.py — voronoi_areas_numpy, voronoi_areas, areas_for_magnification
• autoarray/inversion/mesh/interpolator/delaunay.py — barycentric_dual_area_from, in-graph dual areas, split-point weights
• autoarray/plot/inversion.py — _plot_delaunay
• test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py — existing test__voronoi_areas_via_delaunay_from
• autolens_workspace/scripts/imaging/features/pixelization/source_science.py — the only consumer of areas_for_magnification
• PyAutoLens/autolens/analysis/latent.py — magnification, total_source_flux

Verification

• pytest test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py test_autoarray/inversion/mesh -q green in the worktree; full pytest test_autoarray -q before the PR.
• Every Delaunay area path either matched to an analytic answer in a landed test or flagged with a reproduced failing case here.

Original Prompt

Click to expand starting prompt

Audit the Delaunay pixel-area and magnification source code for correctness bugs

Type: bug
Target: autoarray
Repos:

• PyAutoArray
• PyAutoLens
  Themes:
• pixelization
• euclid
  Difficulty: small-medium
  Autonomy: supervised
  Priority: normal
  Status: formalised
  Consequence: judge
  Review-minutes: 25
  Unattended: ready
  Epic: euclid-dr1-prep
  Phase: 8
  Parent: draft/feature/euclid/euclid_dr1_prep_epic.md
  Filed: 2026-08-28

Phase 8 of the Euclid DR1 preparation epic (was 6c; renumbered 2026-09-01 in the Cortex split,
where the old 6b became PyAutoCortex phases/euclid/magnification_robustness). Can run
alongside the magnification-robustness phase and does not gate on it — that phase is the
empirical half (does μ come out right on data?), this one is the
source-code half (does the code look right?). Either can find the defect first.

User request (verbatim, the follow-up clause of 6b):

"""
We have also never really validated or tested the magnifcitations using
the Delaunay source model, and this could even have bugs (E.g due to pixel areas not being quite right). So, for this
Euclid epic do all of the magnification comparisons you can on the 10 lenses in this spirit, but also have a follow up issue
which checks if the soruce code itself looks like it might have a bug or issue with the Delaunay area and thus magnification calculations.
"""

Where to look

• PyAutoArray mesh geometry: the Delaunay areas_for_magnification implementation
  (voronoi_areas, with boundary cells zeroed via a -1 sentinel) and its
  rectangular-adaptive sibling. The cluster epic's audit
  (draft/test/workspaces/mesh_magnification_correctness.md) already established that
  areas_for_magnification exists for only two mesh geometries, has zero library
  callers (it serves only workspace scripts), and has no direct test — only its
  delegates are tested. That is exactly the shape a latent bug hides in.
• The boundary-cell semantics. Zeroing boundary cells changes the summed source-plane
  area and therefore μ directly. Is the zeroing correct, or does it bias μ high?
• PyAutoArray/autoarray/plot/inversion.py::_plot_delaunay and the Delaunay mapper —
  check whether the areas used for plotting and the areas used for magnification
  come from the same computation. Divergence there is a classic source of "the picture
  looks right but the number is wrong".
• The interaction with sqrt on dual areas — there is a known NaN hazard in the Delaunay
  dual-area gradient path; check whether the same expression is on the magnification
  path.

Method

This is an audit, not a rewrite. Read the code, then prove each suspicion or drop it:

• Construct known-answer configurations (a uniform triangulation of known total area,
  a mesh with an analytically computable source-plane area) and check the code returns
  the analytic value. A resemblance verdict is not verification.
• Where a suspected bug is found, reproduce it in a minimal test before proposing a
  fix, and check whether it is reachable from the paths Euclid actually uses.
• Do not exempt a case because it "looks intentional" — if boundary zeroing is
  deliberate, find the commit or comment that says so.

Deliverables

1. A written audit: each area/magnification code path, what it computes, and whether it
   is correct.
2. Known-answer tests for Delaunay areas_for_magnification (the missing direct test),
   landed regardless of whether a bug is found — the absence of a test is itself the
   defect the cluster-epic audit flagged.
3. If a real defect is found: file a separate bug prompt for the fix rather than
   growing this one, and tell phase 6b immediately so its Delaunay leg is interpreted
   correctly.

Acceptance / gate

• Every Delaunay area/magnification path either verified against an analytic answer or
  flagged with a reproduced failing case.
• Direct tests exist where there were none.
• A clear statement to 6b: are its Delaunay magnification numbers trustworthy?

[Jammy2211]

Audit deliverable (phase 8)

PR: #523 (tests + docstring only, no behavioural change)
Reproduction scripts are in the session scratchpad (part1_flux_integral.py … part3b_fit_latent.py); they were not committed.

Phase 8 audit — Delaunay pixel areas and the magnification calculation

Read-only source audit + numerical experiments. No repository was modified.
All canonical checkouts on main (PyAutoArray, PyAutoLens, PyAutoGalaxy,
PyAutoFit, autolens_workspace, euclid_strong_lens_modeling_pipeline).
Everything below is reproducible from the scripts in this directory
(part1_flux_integral.py, part1b_irregular_hull.py, part1c_pedestal.py,
part1d_zeroed_ring.py, part1e_rectangular_control.py,
part2_plot_consistency.py, part3_latent_reach.py, part3b_fit_latent.py)
run as NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/matplotlib python <script>,
numpy path (xp=np) throughout.

Answer to the epic owner's question in one line: yes, the source code has a
real bug on the Delaunay magnification path. areas_for_magnification returns
Voronoi cell areas, but the Delaunay mapper is a barycentric-linear
interpolant whose exact quadrature weights are the barycentric dual areas.
The two differ by up to five orders of magnitude on individual boundary cells,
and the resulting magnification is biased by −13 % to −99 % on realistic
irregular source meshes.

---

1. Per-path table

Path	File:line	What it computes	Consumers	Verdict	Evidence
MeshGeometryDelaunay.areas_for_magnification	PyAutoArray/autoarray/inversion/mesh/mesh_geometry/delaunay.py:195-207	voronoi_areas with -1 entries set to 0.0	the four source_science.py scripts (autolens_workspace/scripts/{imaging:399,group:404,multi_galaxy:411,interferometer:387}/features/pixelization/source_science.py) — only if the user swaps the mesh to Delaunay; nothing in PyAutoLens/PyAutoGalaxy/Euclid	BIASED (wrong quadrature) + UNTESTED	μ error +0.03 %…+9.3 % (pinned-square hull), −0.06 %…−99.3 % (realistic irregular hull), −13 %…−53 % (adaptive-style mesh). grep areas_for_magnification over PyAutoArray/test_autoarray, PyAutoLens/test_autolens, PyAutoGalaxy/test_autogalaxy, autolens_workspace_test → zero hits
voronoi_areas_numpy	.../mesh_geometry/delaunay.py:13-91	scipy Voronoi + shoelace per cell; unbounded regions marked -1.0 at :89	voronoi_areas (:191) → areas_for_magnification	correct arithmetic, wrong quantity for this use	shoelace verified (abs, correct wraparound, padded terms contribute 0). Its only test, PyAutoArray/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py:41-57, pins the pathology: it asserts voronoi_areas[3] == 29.836324 on a mesh whose real cells are O(1)
Zeroing guard areas[areas == -1] = 0.0	.../mesh_geometry/delaunay.py:205	zeroes only strictly unbounded cells	as above	INSUFFICIENT	the pathological cells are bounded. On 6 seeds of a jittered 24×24 mesh: 15-21 unbounded cells zeroed, but 37-50 bounded cells with area > 5× their dual area, max ratio 1.2 × 10⁵, largest single cell 973 arcsec² on a hull of area 4.15 arcsec² (234× the whole mesh)
barycentric_dual_area_from	PyAutoArray/autoarray/inversion/mesh/interpolator/delaunay.py:342-398	Σ (triangle area / 3) per vertex	scipy_delaunay:39 and jax_delaunay:317 → only the areas_factor·sqrt(areas) split-point offsets	CORRECT — and it is the area the flux integral needs	F_dual / F_map = 1.000000 to ≤ 2 × 10⁻⁵ across 5 source configurations and 3 mesh resolutions (Part 1)
pixel_weights_delaunay_from	.../interpolator/delaunay.py:509-586	barycentric weights (:546-560); nearest-vertex [1,0,0] fallback outside the hull (:568-580)	Mapper.mapping_matrix	correct, but the out-of-hull fallback is a second, separate flux hazard	with the data grid extending past the hull, F_map exceeded the true flux by 4.7×–16.7× (first run of Part 1, before the hull was pinned)
Mapper.mapping_matrix	PyAutoArray/autoarray/inversion/mappers/abstract.py:255-279	data-pixel × source-pixel weights	mapped_reconstructed_(operated_)data, i.e. the magnification numerator	correct	row sums = 1 exactly; Σ colsum · pixel_area = hull area to 1 × 10⁻⁷
_plot_delaunay	PyAutoArray/autoarray/plot/inversion.py:269-322	ax.tripcolor(x, y, vals) — no triangles=, no shading=	every Delaunay source-plane figure	consistent with the mapper, inconsistent with the magnification denominator; docstring wrong	integral of the plotted surface / F_dual = 1.00000000; / F_vor = 18.47 (Part 2)
MeshGeometryRectangular.areas_for_magnification	.../mesh_geometry/rectangular.py:472-481 → areas_transformed	transformed rectangular cell areas	the actual shipped configuration of all four source_science.py scripts (RectangularBilinearAdaptDensity)	FLAGGED — not established, needs its own check	identity-lens μ = +6.2 % (RectangularUniform), +28.5 % (RectangularBilinearAdaptDensity); per-pixel effective/area is a tight scale factor (median 1.061 / 1.279), consistent with a guard-ring node-count mismatch — see caveat in §5
latent.magnification	PyAutoLens/autolens/analysis/latent.py:218-229	total_lensed_source_flux_mujy / total_source_flux_mujy	Euclid pipeline vis_lp, vis_pix, full_model	UNREACHABLE for any pixelized stage (never touches areas_for_magnification)	Delaunay fit on the shipped simulated dataset: total_source_flux = 0.0 → magnification = inf; 9/9 production vis_pix results record magnification = 0.0
Galaxy.image_2d_from	PyAutoGalaxy/autogalaxy/galaxy/galaxy.py:245-278	return xp.zeros(...) at :278 when the galaxy has no non-linear LightProfile	total_source_flux(_mujy)	the root cause of the above	measured: pixelization-only galaxy → sum = 0.000000e+00; Sersic galaxy → 7.941301e+01
areas_factor · sqrt(areas)	.../interpolator/delaunay.py:45 (numpy), :323 (JAX)	split-point offset distance	ConstantSplit / AdaptSplit regularisation matrices	NOT ON PATH — cannot change the flux integral	see §4

---

2. PART 1 — which area integrates the reconstruction?

Setup

Identity lens (source-plane data grid = image-plane data grid), so the true
magnification is exactly 1.0. A Delaunay mapper is built from a jittered
24×24 lattice; a Gaussian (or random positive vector) is sampled at the mesh
vertices as reconstruction. Four integrals:

• F_map = Σ over data pixels of (mapping_matrix @ s) × pixel_area — what the
  mapper paints on the sky, i.e. what mapped_reconstructed_data sums to
• F_vor = Σ s · areas_for_magnification — the current workspace denominator
• F_dual = Σ s · barycentric_dual_area_from(points, simplices)
• F_true = analytic Gaussian flux 2πσ²A

Two confounds were removed by construction so the numbers mean what they say:

1. Data pixels outside the convex hull take the nearest-vertex fallback
   (interpolator/delaunay.py:568-580) and dump flux onto boundary vertices.
   The outer ring of mesh vertices is therefore pinned un-jittered to the data
   boundary, and n_fallback == 0 is reported in every run.
2. F_map integrates the data-grid footprint while F_dual integrates the hull.
   Pinning the hull to the data square makes them the same region — reported as
   hull area 4.000000 vs data footprint 4.000000 (match to 0.0000 %).

Only interior vertices are jittered, so the triangulation stays irregular.

Result (part1_output.txt)

config	flux in outer ring	F_map	F_vor	F_dual	F_true	μ(vor)	μ(dual)
(i) σ=0.10″, hull 1.0″	0.0000	0.06271824	0.06266253	0.06271942	0.06283185	1.000889 (+0.089 %)	0.999981 (−0.002 %)
(i-b) σ=0.20″	0.0004	0.25135922	0.25128810	0.25135685	0.25132741	1.000283 (+0.028 %)	1.000009 (+0.001 %)
(ii) σ=0.50″	0.1823	1.42951462	1.40011211	1.42950317	1.57079633	1.021000 (+2.100 %)	1.000008 (+0.001 %)
(ii-b) σ=1.00″	0.3868	2.92569264	2.74590086	2.92567536	6.28318531	1.065476 (+6.548 %)	1.000006 (+0.001 %)
(iii) random positive	0.4723	2.33212568	2.13409902	2.33213093	—	1.092792 (+9.279 %)	0.999998 (−0.000 %)

Resolution check (σ=0.20″): F_dual/F_map = 0.99999083 / 0.99999959 / 1.00000751
at n_side = 16 / 24 / 40.

Which denominator gives 1.0: F_dual. It reproduces the identity-lens truth
to ≤ 2 × 10⁻⁵ in every case, including the random reconstruction — the identity
is geometric, not a property of the source. F_vor never does, and its error
grows monotonically with how much reconstruction weight sits near the boundary
(0.089 % → 9.3 % as the outer-ring weight goes 0.0000 → 0.4723).

F_map/F_true is a separate, honest discretisation error of the mesh itself
(0.998 → 0.466 as the source outgrows the mesh); it is not a denominator defect,
which is exactly why μ and not F_true is the right yardstick here.

The realistic case is far worse (part1b_irregular_hull.py, part1b_output.txt)

Because F_dual == F_map is now an established identity, F_dual / F_vor is
the magnification bias, with no matched-domain bookkeeping needed. On a fully
irregular hull (outer ring jittered too — what an adaptive source mesh looks
like after tracing), scipy's Voronoi hands some bounded boundary sites
enormous cells, because their far vertices are the circumcentres of very flat
boundary triangles. areas_for_magnification zeroes only the unbounded ones.

mesh / source	zeroed (unbounded)	bounded cells with area > 5× dual	max area/dual	largest cell	μ bias, 6 seeds
jittered 24×24, compact σ=0.15″ (flux deep inside)	15-21	37-50	5.9 × 10³ – 1.7 × 10⁴	36-137 arcsec² (hull 4.15)	−0.31 % … +0.41 %
same meshes, extended σ=0.80″ (flux at the boundary)	15-21	42-50	up to 1.2 × 10⁵	up to 973 arcsec² = 234× the entire hull	−99.30 % … −95.03 %
adaptive-style mesh (500 points drawn from the source's own brightness), σ=0.30″	8-11	2-7	7.8 – 286	0.21-3.31 arcsec² (hull 1.5-1.9)	−52.96 % … −13.35 %

Direction: negative — the inflated boundary cells swell the denominator, so
source_science.py under-estimates the magnification whenever the
reconstruction is non-zero at the mesh edge. (The positive branch, +0.03 %…+9.3 %,
occurs in the opposite regime, where the boundary cells are unbounded and get
zeroed, so their flux is discarded from the denominator instead.)

The "compact source is fine" case is a knife edge (part1c_pedestal.py, part1c_output.txt)

The ±0.4 % agreement for a compact source is an accident: the pathological cells
multiply a reconstruction value of ≈ 0. A real inversion does not put exactly 0
there. Adding a random pedestal to the vertices, as a fraction of the peak:

pedestal / peak	seed 0	seed 3	seed 5
0	+0.41 %	+0.11 %	−0.05 %
1e-4	−2.82 %	−5.05 %	−3.48 %
1e-3	−47.40 %	−32.24 %	−77.30 %
1e-2	−88.08 %	−81.51 %	−96.81 %
1e-1	−96.36 %	−96.15 %	−99.42 %

A pedestal of one part in ten thousand of the peak already costs several per
cent. There is no regime in which this denominator is safe by design.

zeroed_pixels does not rescue it (part1d_zeroed_ring.py, part1d_output.txt)

The Euclid pipeline builds its Delaunay source mesh with
zeroed_pixels=edge_pixels_total (30), a ring of image-plane points appended by
al.image_mesh.append_with_circle_edge_points
(euclid_strong_lens_modeling_pipeline/scripts/full_model.py:730-734, 786, 828-832, 863;
scripts/initial_lens_model.py:592) whose source values are fixed to exactly 0.
Peeling convex hulls shows why that is not enough:

seed 0, pathological bounded cells (area > 5x dual) by hull-ring depth:
   ring 0   :  21 vertices,   0 pathological   <- the convex hull: ALL already unbounded and zeroed
   ring 1   :  22 vertices,  14 pathological, max area/dual 17374.9, sum(vor) 161.6 vs sum(dual) 0.12
   ring 2   :  23 vertices,  12 pathological, max area/dual   777.1, sum(vor)   8.0 vs sum(dual) 0.14
   ring 3   :  23 vertices,   6 pathological, max area/dual    36.6
   interior : 487 vertices,   5 pathological

The pathology lives one ring inside the convex hull, which zeroed_pixels
never touches. Zeroing one hull ring moves the flat-source bias from −97.63 % to
−97.69 % (seed 0); −97.03 → −97.09 (seed 1); −97.40 → −97.46 (seed 2);
−97.93 → −97.97 (seed 3). For a compact source with a 0.1 % pedestal, zeroing one
ring moves −7.76 % → −13.65 % (seed 0), −19.38 % → −18.88 % (seed 1) — i.e. it
does nothing useful in either direction. Two rings are needed before the number
comes back (−1.68 %, −0.62 %, −0.90 %, −1.97 %), and three before it is clean.

Is the mapping barycentric-linear?

Yes, and this is what makes Σ sᵢ · dualᵢ the exact integral rather than an
approximation. pixel_weights_delaunay_from
(PyAutoArray/autoarray/inversion/mesh/interpolator/delaunay.py:509-586)
computes the three sub-triangle areas a0, a1, a2 at :546-552 and normalises
by their sum at :556-560 — textbook barycentric coordinates. For a piecewise
linear function on a triangulation,
∫_hull f = Σ_tri area_tri · mean(3 vertices) = Σ_i s_i · (Σ area/3) = Σ_i s_i · dual_i,
which is precisely barycentric_dual_area_from (:342-398). The one departure
is the out-of-hull nearest-vertex fallback at :568-580, which is piecewise
constant, not linear — a separate hazard, quantified above.

Incidental: scipy_delaunay pads simplices to (2N, 3) with -1 rows for
static JAX shapes (:24). -1 indexing wraps to the last vertex, giving a
degenerate zero-area triangle, so passing the padded array to
barycentric_dual_area_from is harmless — verified identical to the unpadded
result (part2_output.txt).

---

3. PART 2 — plotting vs magnification consistency

_plot_delaunay (PyAutoArray/autoarray/plot/inversion.py:269-322) calls
ax.tripcolor(x, y, vals, cmap=colormap, norm=norm) at :319 with neither a
triangles= nor a shading= argument. Two things follow. First, matplotlib
3.10's tripcolor defaults to shading='flat' — not Gouraud, contradicting
the docstring at :272; with point-valued C, flat shading paints each triangle
with the mean of its three vertices, whose area-weighted sum is
Σ_tri area·mean(3) = Σ_i s_i·dual_i — identical to the Gouraud integral by
linearity, so the docstring inaccuracy is cosmetic here but should not be relied
on. Second, without triangles=, matplotlib rebuilds its own Qhull
triangulation of the same points instead of reusing
mapper.interpolator.delaunay.simplices; on the test mesh the two simplex sets
were identical (780 vs 780 triangles, 0 differences), so this is a latent
fragility (co-circular points could diverge) rather than an active defect.
The measured consequence is unambiguous: the integral of the plotted surface
equals F_dual to 1.00000000 and exceeds F_vor by a factor of 18.47
(part2_output.txt). The picture the user looks at and the mapper's own flux
already agree with each other and disagree with the magnification denominator;
areas_for_magnification is the odd one out, not the plot.

---

4. PART 4 — the sqrt / NaN path

grep -n sqrt PyAutoArray/autoarray/inversion/mesh/mesh_geometry/delaunay.py
returns nothing — there is no square root anywhere on the
voronoi_areas_numpy → voronoi_areas → areas_for_magnification path, so no
NaN can be manufactured there by a negative argument (the shoelace at :83
takes np.abs before any use, and -1.0 at :89 is a sentinel that :205
maps to 0.0, never a sqrt input). The only sqrt in the Delaunay package is
areas_factor * np.sqrt(areas) at interpolator/delaunay.py:45 (numpy) and
:323 (JAX), and it is provably not on the flux path: its output is
area_weights, which split_points_from
(PyAutoArray/autoarray/inversion/regularization/regularization_util.py:33-84)
uses solely as an offset distance to place four cross points around each
vertex; those points feed splitted_mappings / splitted_weights, which are
consumed only by the split regularisation schemes
(autoarray/inversion/regularization/adapt_split.py:127-139, reg_split_np_from).
They build the regularisation matrix H, which shapes which reconstruction the
inversion prefers but never rescales a reconstruction that has already been
solved, and they appear nowhere in mapping_matrix,
mapped_reconstructed_data, or areas_for_magnification. So the split-point
weighting cannot change the integral of the reconstructed values — it can only
change the values themselves via the fit, which is a modelling choice, not an
arithmetic defect.

---

5. PART 3 — Euclid reachability

(a) Latent outputs

euclid_strong_lens_modeling_pipeline/output/ does not exist in the canonical
checkout — .gitignore:13-14 excludes output/ and output_*/, and the RAL
results are not present there. Production outputs were found on this machine
under /tmp/scoped_output, /tmp/scoped_output_sed, /tmp/scoped_323738
(read-only; 523 latent_summary.json files across the dr1_prelim_grade_ab
sample). Only two initial_lens_model stages exist in them —
there are no full_model outputs on this machine. Values are the
median_pdf_sample entry of latent_summary.json:

stage	results	magnification present	value
initial_lens_model/vis_lp (Sersic source)	9	9/9	6.2727, 10.0782, 18.4297, … — real numbers
initial_lens_model/vis_pix (Delaunay source)	9	9/9	0.0 in every tile, alongside total_source_flux = 0.0
sersic_lens_model/{vis,nir_y,nir_j,nir_h,decam_grizt}	455	455/455	real numbers (2.93 – 28.80)

So the vis_pix magnification is not merely absent — it is present as a
meaningless 0.0 sentinel that a downstream catalogue would silently ingest.
The vis_pix key set is otherwise complete (9 keys, same as vis_lp), so the
NaN-drop signature does not fire on these particular files.

(b) Analytic reading, confirmed empirically

total_source_flux_mujy (PyAutoLens/autolens/analysis/latent.py:189-215)
computes tracer.galaxies[-1].image_2d_from(grid=fit.dataset.grids.lp) on
fit.tracer_linear_light_profiles_to_light_profiles. Galaxy.image_2d_from
(PyAutoGalaxy/autogalaxy/galaxy/galaxy.py:245-278) is:

if len(self.cls_list_from(cls=LightProfile, cls_filtered=LightProfileLinear)) > 0:
    return sum(self.image_2d_list_from(grid=grid, xp=xp, operated_only=operated_only))
return xp.zeros((grid.shape[0],))        # galaxy.py:278

A Pixelization is not a LightProfile, so a source galaxy whose only light
component is a pixelization takes the return zeros branch. Measured
(part3_output.txt): pixelization-only galaxy → sum(image_2d_from) = 0.000000e+00;
Sersic galaxy → 7.941301e+01. Then ab_mag_via_flux_from(0.0, magzero=25) = inf
and flux_mujy_via_ab_mag_from(inf) = 0.0 exactly, so the denominator is a hard
zero and magnification = lensed / 0.0 = +inf.

I confirmed this empirically rather than resting on the reading.
part3b_fit_latent.py builds a real al.FitImaging on the pipeline's shipped
dataset/simulated/euclid_dr1_like (magzero 24.6), with the truth lens from
truth.json, and evaluates the library latents directly — no search:

source model	total_source_flux	total_source_flux_mujy	magnification
Sersic (the vis_lp control)	0.8679858287	0.4555254383	15.146299 — bit-identical to truth.json's 15.146298766165936
Delaunay pixelization + zeroed_pixels=30 (vis_pix / full_model)	0.0	0.0	inf

inf is non-finite, so PyAutoFit/autofit/non_linear/analysis/latent.py:258-263
(if np.isfinite(all_values[:, i]).any()) drops the whole magnification
column: absence of the key is the signature under the current code, and 0.0
under whatever version produced the archived vis_pix results. Either way
the Euclid magnification latent can never report a magnification for a
Delaunay stage, and it never touches areas_for_magnification at all — the two
Delaunay-magnification defects are independent and neither masks the other.

The Sersic control also validates the latent itself: it reproduces the
simulator's recorded truth exactly, so the definition is sound and only the
pixelized-source route is broken.

Note also scripts/tools/diagnose_latent_vis_pix.py's docstring claim that
"the library latents that do not depend on the source reconstruction
(magnification, …) are identical between the two stages" is incorrect for
magnification: it depends on the source model, which is exactly what changes
at vis_pix, and the production numbers (6.27/10.08/18.43 vs 0.0/0.0/0.0)
disprove it directly.

Caveat on the rectangular control (part1e_rectangular_control.py)

The shipped source_science.py scripts all use
al.mesh.RectangularBilinearAdaptDensity, whose areas_for_magnification is
areas_transformed (mesh_geometry/rectangular.py:472-481) — a different code
path from the Delaunay one this audit was scoped to. Running the same
identity-lens test there gives μ = +6.16 % (RectangularUniform) and
+28.47 % (RectangularBilinearAdaptDensity) at σ=0.15″, with the per-source-pixel
ratio effective_mapper_area / areas_for_magnification a tight scale factor
(median 1.0607 and 1.2785, spread ~0.3 % away from the boundary) rather than the
Delaunay path's wild per-cell scatter. That signature is consistent with a
guard-ring node-count mismatch — the rectangular docstring at :490-497 notes a
guard ring occupies the border, so an n vs n−3 denominator would give
(24/21)² ≈ 1.31. I did not establish a matched-domain identity for the
rectangular meshes and I am not claiming a bug there — I ran it only so the
Delaunay verdict is stated against a control rather than in isolation. It is a
lead for the cluster epic's mesh_magnification_correctness prompt, not a
finding of this audit.

---

6. Statement to Cortex phase 7 (magnification_robustness)

Delaunay magnification numbers produced by the source_science.py recipe are
not trustworthy as they stand, and the defect is in the library, not the
script. MeshGeometryDelaunay.areas_for_magnification
(PyAutoArray/autoarray/inversion/mesh/mesh_geometry/delaunay.py:195-207)
returns Voronoi cell areas, but the Delaunay mapper is a barycentric-linear
interpolant whose exact quadrature weights are the barycentric dual areas
already computed one module away in barycentric_dual_area_from
(inversion/mesh/interpolator/delaunay.py:342-398) and currently used only to
position regularisation split points. The area that should be used is the
barycentric dual area: with an identity lens it recovers μ = 1.0 to ≤ 2 × 10⁻⁵
for every source shape tested, including a random reconstruction, and it is also
exactly what the plotted tripcolor surface integrates to — so switching the
denominator makes the number, the picture and the mapper agree, where today only
the first disagrees. The expected bias direction is negative — magnification
systematically under-estimated — whenever the reconstruction is non-zero near
the mesh edge, because areas_for_magnification zeroes only the unbounded
boundary cells while the pathology sits in the first ring inside the convex
hull, where bounded Voronoi cells reach 10⁵× their dual area (one cell measured
at 234× the entire mesh hull). Magnitude: −13 % to −53 % on adaptive-style
meshes, −95 % to −99 % when the source fills the hull, and −3 % to −5 % even for
a compact central source once its edge pixels carry a pedestal of one part in
ten thousand of the peak; the apparent ±0.4 % accuracy for a perfectly compact
source is a knife edge that a real non-negative solver will not sit on. Three
operational consequences for the phase. (1) The Delaunay rung (rung 5) cannot be
scored against a latent.magnification number from the pipeline: for any
pixelized source total_source_flux is a hard 0.0
(PyAutoGalaxy/autogalaxy/galaxy/galaxy.py:278), so the latent evaluates to
inf and is dropped, and 9/9 archived vis_pix results record it as a
meaningless 0.0 — that leg needs the area-based A_img/A_src quantity
computed explicitly, and the phase's "state which magnification is compared at
each rung" instruction is therefore load-bearing, not bookkeeping. (2) The
Sersic rungs (1-4) are unaffected by this bug and their latent is sound — it
reproduced the simulator truth 15.146298766165936 bit-for-bit — so a
matched-vs-mismatched comparison across rungs 1-4 is interpretable today.
(3) areas_for_magnification has no test in any repository, and the only
test that touches voronoi_areas pins a 29.8 arcsec² boundary cell as expected
behaviour, so any fix needs a regression test built from the identity-lens
construction in part1_flux_integral.py rather than a re-pin of the existing
snapshot. A separate, un-established lead worth a line in the phase's runs: the
rectangular meshes that source_science.py actually ships with show a
consistent +6 %/+28 % scale offset with a guard-ring signature; that belongs to
the cluster epic's mesh_magnification_correctness prompt, and this audit
neither confirms nor rules it out.

(No fix is proposed or implemented here, per the audit's brief.)

[Jammy2211]

Shipped

PR #523 merged into main (merge commit 548ff1e, head 8c2e0d18) — CI green on all three legs (unittest 3.12, 3.13, nojax).

What landed: the audit's test deliverable and one docstring correction, no behaviour change.

• Three known-answer tests for MeshGeometryDelaunay.areas_for_magnification (unit lattice: interior cells 1.0, hull cells zeroed, sum (n-2)²; bounded boundary cells kept; repeat calls agree without writing the -1 sentinel back).
• Four tests for barycentric_dual_area_from (single triangle A/3; sum equals the convex-hull area; exact integration of a linear field; NumPy-vs-JAX parity).
• areas_for_magnification docstring now states the proven semantics: Voronoi cell areas, only unbounded cells zeroed, not the interpolant's barycentric dual areas.

Audit verdict (full table in the comment above): two real defects, deliberately not fixed here, filed as separate PyAutoMind bug prompts —

• draft/bug/autoarray/delaunay_magnification_uses_voronoi_not_dual_areas.md — the Voronoi denominator biases μ; the exact quadrature weight for the barycentric-linear mapper is the dual area.
• draft/bug/autolens/magnification_latent_zero_for_pixelized_source.md — the pipeline's magnification latent is a hard 0/0 for any pixelization-only source.

The two tests that pin the current Voronoi semantics carry comments so the follow-up fix flips them deliberately.

Closing as completed; the Mind record is PyAutoMind/complete/2026/09/delaunay-area-magnification-audit.md. pending-release stays on the PR until a PyAutoArray release publishes it.

---

Generated by Claude Code