Activity

1. Jammy2211 commented on Aug 28, 2026

   

   Jammy2211

   on Aug 28, 2026

   CollaboratorAuthor

   More actions

   Step 0 checkpoint — instrumented split + measured geometry

   Harness commit 3f69b34 (autolens_profiling, feature/numba-hst-curvature-matrix-phase2).
   All three cells run with OMP_NUM_THREADS=1 AUTOARRAY_NUMBA_OPERATED_MEMO=0, n_repeats 10,
   against the task worktree's PyAutoArray (PYTHONPATH-pinned, autoarray.__file__ asserted).
   First post-rebuild run discarded (numba cache: true).

   Per-step split (s/eval)

   Step	hst rect (bilinear, 784)	euclid rect (bilinear, 784)	hst delaunay (Hilbert, 1250)
   Blurred image (FFT convolve)	0.010794	0.005689	0.011387
   Inversion build (trace+mesh+mapper)	0.026765	0.012052	0.028456
   PSF-weighted data [numba]	0.009156	0.001931	0.009175
   MGE: image_2d_from(grid) x60	0.116026	0.030895	0.116086
   MGE: image_2d_from(blurring_grid) x60	0.045998	0.024407	0.046184
   MGE: batched PSF convolution	0.054078	0.014001	0.064275
   MGE residual (dict assembly)	-0.007544	0.001098	-0.003476
   (MGE row total, un-split)	0.208558	0.070401	0.223069
   Mapper sparse triplets	0.004870	0.001153	0.019793
   Data vector D [numba]	0.004511	0.000750	0.005662
   F: mapper x mapper [sparse-op]	0.255055	0.049054	0.130403
   F: mapper x linear-func [FFT conv + scatter]	0.049128	0.014918	0.054698
   F: linear-func x linear-func [BLAS]	0.018941	0.003284	0.019782
   F residual (mirror + diag-add)	-0.006710	0.000521	-0.013915
   Regularization matrix H	0.001925	0.002463	0.072216
   Reconstruction solve (BLAS)	0.009757	0.006513	0.044785
   log det (F+H) + log det H	0.016614	0.013514	0.058052
   Log evidence (figure of merit)	0.015803	0.012274	0.059074
   TOTAL (step-by-step)	0.633439	0.197839	0.734423
   Direct log_likelihood_function	0.585031	0.175982	0.663736

   The four MGE rows sum to the raw un-split MGE step, as the four F rows do (both
   residuals are differences of averaged timings and can come out slightly negative;
   recorded as measured, never clipped). mge_split_reproduces_step: true on all three
   cells — the three timed pieces reconstruct the step's dict bit-identically, so
   the attribution is not an inference.

   Geometry constants (model-independent, recorded once into the JSON)

   	hst rect	euclid rect	hst delaunay
   PSF shape (native)	21x21	21x21	21x21
   image pixels (masked)	15361	3841	15361
   sparse_operator stored pairs = lengths.sum()	11,080,573	2,318,869	11,080,573
   lengths.mean() / .max()	721.34 / 841	603.71 / 841	721.34 / 841
   data_to_pix_unique.shape	(15361, 4)	(3841, 4)	(15361, 3)
   pix_lengths.mean() / .max() / .sum()	4.0 / 4 / 61444	4.0 / 4 / 15364	1.5523 / 3 / 23845
   mapper.params (n_source)	784	784	1250
   mapper x mapper inner ops	1.7729e8	3.7102e7	3.2692e7
   two-stage stage-1 scatters	4.4322e7	9.2755e6	1.8187e7
   two-stage stage-2 AXPY	4.8172e7	1.2045e7	2.9806e7
   two-stage accumulator clear	1.2043e7	3.0113e6	1.9201e7
   two-stage total ops	1.0454e8	2.4332e7	6.7195e7

   mapper x mapper inner ops is the exact count of
   curvature_matrix_via_sparse_operator_from's innermost accumulation,
   sum over stored (data_0, data_1) pairs of u0(data_0) * u1(data_1) — not an estimate.

   Verdict

   (a) MGE — profile evaluation, not convolution. At HST the MGE row is 0.2086 s, of
   which 0.1620 s (78 %) is the 120 numpy profile evaluations (0.1160 s on the
   over-sampled data grid + 0.0460 s on the blurring grid) and 0.0541 s (26 %) is the single
   batched real-space convolution; the residual is -0.0075 s. The plan's inference from #505
   ("~0.17 of the 0.22 s is the 120 per-profile image evaluations") is confirmed within
   5 %. Step 3's shared-geometry fast path in PyAutoGalaxy is therefore aimed at the right
   0.16 s, and the ConvolverState.fft_kernel idea correctly stayed out of scope. The 2.5x
   asymmetry between the grid and blurring-grid evaluations is the 4x4 radial over-sampling
   of the data grid, not a per-call overhead.

   (b) mapper x mapper op count — confirmed, but the two-stage margin is dataset-dependent.
   HST rectangular measures 1.773e8 innermost accumulations against the plan's ~2e8: right
   to within 12 %. At 0.2551 s that is 1.44 ns/accumulation, ~4.6 cycles — the signature of an
   irregular read-modify-write into a 4.9 MB matrix, exactly the cost the two-stage
   reformulation targets. Its op count on the same geometry is 1.045e8, a 1.70x reduction
   before any locality or vectorisation gain, so the plan's 1.7-2.5x estimate for HST
   rectangular stands.

   The Delaunay cell is the caveat, and it is a new finding from this measurement: at
   1250 vertices with pix_lengths.mean() = 1.55 (barycentric, vs 4.0 for bilinear) the
   reference kernel does only 3.27e7 ops while the two-stage form does 6.72e7 — 2.05x
   more. The crossover is not a function of n_source alone: it depends on n_source,
   pix_lengths.mean() (u0) and lengths.mean() (P) together, because stage 2 costs
   u0 x n_source per data pixel regardless of how sparse the mapping is. Step 2 will
   measure whether vectorised dense AXPYs still win at that op ratio, and will report the
   measured crossover; the branch will be an explicit, documented threshold, never a silent
   heuristic.

   Proceeding to step 1 (oracle + bit-identical hoist).
2. Jammy2211

   added a commit that references this issue on Aug 28, 2026

   style: ruff format the two numba breakdown harnesses (#507)

   ...

   7cc4bbd
3. Jammy2211 commented on Aug 28, 2026

   

   Jammy2211

   on Aug 28, 2026

   CollaboratorAuthor

   More actions

   Phase 2 complete — steps 3 and 4, and the result

   Continuing from the step-0 checkpoint and the step-1/2
   commits. The goal is met: HST rectangular 0.6214 -> 0.3334 s per evaluation
   (step total), 0.6184 -> 0.3013 s directly timed — 1.86x / 2.05x, against the
   target of ~0.35 s.

   Step 3 — the MGE operated matrix

   3a, PyAutoArray d8bc3bac — the OverSampler divisor.
   binned_array_2d_from recomputed np.bincount(segment_ids) and the
   sub_is_uniform np.isclose sweep on every call — 120 calls per evaluation for
   a 60-Gaussian MGE, and once more for every ordinary numpy light profile. Both
   depend only on the constructor's sub_size; both are now cached_property. On
   the HST over sampler (17980 sub-pixels into 15361 pixels) the divisor cost
   38.8 us per call and the uniformity check 61.3 us, against 81.5 us for the whole
   cached call. Bit-identical. The cached divisor is writeable = False and built
   with np.maximum(counts, 1), because the old zero-guard mutated in place and
   the array is now shared by every caller. OverSampler is pickled to Nautilus
   workers, so a pickle round-trip is asserted.

   3b, PyAutoGalaxy 6031aa0a — one eccentric-radius grid per MGE basis. This
   is where the MGE time actually was. Each of the 60 Gaussians re-ran the same
   reference-frame transform (arctan2 + sin + cos per coordinate) and formed
   the same eccentric radii from it; cProfile on the HST cell put the transform
   alone at 951 us of the 1.3 ms each profile took.
   LightProfileLinearObjFuncList._image_slim_list_from now groups profiles by
   (class, centre, ell_comps) and computes each group's transform and radii once.
   Bit-identical. Two decisions worth recording:

   • Grouping, not an all-or-nothing check. The workspace's canonical MGE
     recipe stacks two sets of 30 Gaussians which share a centre but carry their own
     ell_comps, and GalaxiesToInversion puts both in one func list. An
     all-or-nothing test would have fallen back to 60 independent evaluations on
     exactly the model users are told to write.

   • Sersic is deliberately excluded, even in a same-geometry Sersic basis: it
     already avoids the polar transform via _eccentric_radii_grid_from_cartesian,
     a branch taken only when the grid handed to it is not pre-transformed — so
     hoisting a transformed grid into it would be both slower and a different
     floating-point expression. The membership test is the identity of the
     image_2d_from function object against Gaussian.image_2d_from, which also
     correctly excludes GaussianMultipole.

     Follow-up worth filing: giving Gaussian that same Cartesian shortcut would
     remove the polar transform for every Gaussian evaluation, not just an MGE
     basis. It moves values, so it is a separate task.

   Binning stays a per-profile call (binned_array_2d_from takes a 1-D array), and
   so does the exp — measured, that arithmetic is transcendental-bound rather than
   overhead-bound, so the plan's "vectorise over sigma into one block" would buy
   nothing and would have to be sliced column-wise to bin.

   Step 4 — measured

   Three arms, paired B/A/B/A in one session (base = PyAutoArray main @
   1b89404b + PyAutoGalaxy main @ 0fbe863d; after steps 1-2; final), three
   rounds per cell in rotating arm order, first discarded for the numba recompile,
   mean of the last two. OMP_NUM_THREADS=1, AUTOARRAY_NUMBA_OPERATED_MEMO=0,
   n_repeats=10.

   row	base	step 2	final	step 3	whole phase
   hst rect, step total	0.6214	0.4324	0.3334	1.30x	1.86x
   hst rect, direct eval	0.6184	0.4207	0.3013	1.40x	2.05x
   hst rect, F mapper x mapper	0.2581	0.0819	0.0825	0.99x	3.13x
   hst rect, MGE total	0.1970	0.1937	0.0844	2.29x	2.33x
   hst rect, MGE 60x image_2d_from(grid)	0.1042	0.0995	0.0233	4.28x	4.48x
   hst rect, MGE 60x image_2d_from(blurring)	0.0436	0.0420	0.0115	3.67x	3.81x
   hst rect, MGE batched PSF convolution	0.0512	0.0504	0.0523	0.96x	0.98x
   hst rect, Blurred image (FFT)	0.0101	0.0097	0.0099	0.98x	1.02x
   euclid rect, step total	0.2226	0.1870	0.1417	1.32x	1.57x
   euclid rect, F mapper x mapper	0.0531	0.0188	0.0191	0.98x	2.78x
   euclid rect, MGE total	0.0802	0.0797	0.0316	2.52x	2.54x
   hst delaunay-1250, step total	0.6331	0.6076	0.4884	1.24x	1.30x
   hst delaunay, F mapper x mapper	0.1088	0.0697	0.0676	1.03x	1.61x
   hst delaunay, MGE total	0.1881	0.1932	0.0839	2.30x	2.24x

   The last two hst-rect rows are the control: nothing outside the two terms step 0
   identified moves. The Delaunay cell gains least because its evaluation is now
   spread across terms this issue did not touch — after phase 2 its largest steps
   are the ConstantSplit regularization matrix (0.060 s), the batched MGE
   convolution (0.057 s) and the mapper x linear-func block (0.049 s).

   hst rectangular RTU was re-run once for currency, unpaired: it is a 7.1 s cell
   dominated by its 6.73 s mapper-sparse-triplets step and is GPU-only by the
   2026-08-28 decision.

   Correctness

   Pinned log likelihoods at the explicit rtol=1e-6, PASSED on every arm and
   unchanged to every recorded digit: hst bilinear 27661.91013366411 (pinned
   27661.910133665442), hst RTU 27180.70471569685, hst Delaunay
   29090.527210448134. euclid, unpinned, 6213.306873885871 on all three arms.

   Array-level, same instance on each arm:

   comparison	MGE operated matrix	inversion.curvature_matrix
   step 3 alone (step2 -> final), hst	bit-identical	bit-identical
   step 3 alone, euclid	bit-identical	bit-identical
   whole phase (base -> final), hst	bit-identical	max rel 1.24e-14, allclose(rtol=1e-12) True
   whole phase, euclid	bit-identical	max rel 2.63e-15, allclose(rtol=1e-12) True

   Step 3 is bit-identical by construction in both libraries; the only
   floating-point movement in the whole phase is step 2's reassociation, twelve
   orders inside the pins.

   Tests: pytest test_autoarray 1326 passed (1296 at the phase-1 baseline),
   pytest test_autogalaxy 1149 passed (1144 before). Control runs recorded for
   every rtol test and for every branch of the new step-3b membership test
   (dropping the transform, elliptical instead of eccentric radii, transforming the
   un-over-sampled grid, removing the image_2d_from-identity guard — each fails
   the new tests).

   Pool run

   cpu_fast_modeling.py's first fit, 8 cores, PYAUTO_TEST_MODE=1,
   PYAUTO_SMALL_DATASETS unset (2828 masked pixels), 60-Gaussian linear MGE bulge
   added to the lens. Per evaluation:

   	pool of 8 (median of 3)	serial
   base	0.2000	0.4099
   final	0.1697	0.3634

   Parallel speed-up ratio 2.05x -> 2.14x — flat to slightly up, which is the
   pass condition; it is the number that would fall if any lever had introduced
   hidden threads. End-to-end the pool improves 15 % and the single process 11 % on
   this small workspace dataset, where F and the MGE are a much smaller share of the
   evaluation than in the HST harness cells. Read it as a regression check, not a
   speed-up measurement.

   One harness bug found and fixed

   The MGE: image_2d_from(blurring_grid) x60 row hand-rolled the per-profile loop
   that operated_mapping_matrix_override used to contain. After step 3b that loop
   no longer exists in the library, so the harness would have kept timing the old
   code on the new library and mis-attributed the entire win to the residual row —
   the first final run showed exactly that (a -0.030 s residual). It now calls the
   same path the override calls, dispatching on the attribute so both libraries stay
   measurable, and records which branch it took in the result JSON as
   mge_blurring_stack_path.

   PRs (do not merge out of order)

   • perf(numba): the mapper x mapper block and the MGE operated matrix (#507) #508 — steps 1, 2, 3a
   • perf(linear): share one eccentric-radius grid across an MGE basis (PyAutoArray#507) PyAutoGalaxy#590 — step 3b
   • profiling: phase-2 instrumentation and the measured result (#507) autolens_profiling#190 — steps 0 and 4 (harness, note, artifacts)

   Merge order: #508 -> #590 -> #190. PyAutoGalaxy's CI clones PyAutoArray from
   source and checks out the same-named branch when it exists, so #590 is presently
   testing against the PyAutoArray branch and falls back to main once #508 merges
   and that branch is deleted. No PyAutoGalaxy test depends on the new OverSampler
   behaviour, so #590 is green either way.

   Heart is RED at ship time for a pre-existing, unrelated reason
   ("release validation FAILED (stage integrate)"); the human authorisation to
   ship over it is quoted verbatim in each PR body.
4. Jammy2211

   added a commit that references this issue on Aug 28, 2026

   Merge pull request #508 from PyAutoLabs/feature/numba-hst-curvature-m…

   ...

   f4b8e75
5. Jammy2211 commented on Aug 28, 2026

   

   Jammy2211

   on Aug 28, 2026

   CollaboratorAuthor

   More actions

   Shipped

   Phase 2 is merged. The issue's goal — HST rectangular numba evaluation ~0.60 s → ≤ ~0.35 s with the log likelihood unchanged to pinned tolerance — is met.

   Merged PRs (library-first gate honoured, in this order):

   PR	merge commit
   #508 — steps 1, 2, 3a	f4b8e758 (head d8bc3bac)
   PyAutoLabs/PyAutoGalaxy#590 — step 3b	90317122 (head 6031aa0a)
   PyAutoLabs/autolens_profiling#190 — steps 0 and 4	ed1e265a (head 7cc4bbd7)

   Result (paired B/A/B/A in one session, OMP_NUM_THREADS=1, AUTOARRAY_NUMBA_OPERATED_MEMO=0, n_repeats=10):

   cell	base (1b89404b)	final	whole phase
   hst rectangular, step total	0.6214 s	0.3334 s	1.86x
   hst rectangular, direct eval	0.6184 s	0.3013 s	2.05x
   euclid rectangular	0.2226 s	0.1417 s	1.57x
   hst Delaunay-1250	0.6331 s	0.4884 s	1.30x
   hst: F mapper × mapper	0.2581 s	0.0825 s	3.13x
   hst: MGE total	0.1970 s	0.0844 s	2.33x

   The two levers, both chosen by step 0's instrumentation rather than the prompt's candidate list: a two-stage source-space accumulator for the mapper × mapper block (w0 does not depend on data_1, so the sum factorises into a per-data-pixel dense accumulator plus contiguous AXPYs), and one shared eccentric-radius grid per MGE group in PyAutoGalaxy (78% of the MGE row was profile evaluation, not the PSF convolution — the transform and radii were recomputed for all 60 Gaussians that share centre/ell_comps). Plus a cached OverSampler binning divisor and uniformity check, called 120× per evaluation.

   Verification: test_autoarray 1326 passed, test_autogalaxy 1149 passed. All three pinned log likelihoods unchanged to every digit at rtol=1e-6 (hst bilinear 27661.91013366411, hst RTU 27180.70471569685, hst Delaunay 29090.527210448134). curvature_matrix bit-identical for steps 1 and 3, allclose(rtol=1e-12) across the phase. Nautilus pool ratio 2.05x → 2.14x on 8 cores — no oversubscription.

   CI green on every run and every matrix leg before each merge: PyAutoArray Tests [pull_request] (unittest (3.12), unittest (3.13), unittest-nojax); PyAutoGalaxy Tests [pull_request] (same three legs) + Docs [pull_request] (docs-build); autolens_profiling lint [pull_request].

   Heart RED at merge — human override. Merged over the pre-existing Heart RED (pyauto-heart readiness --json 2026-08-28T21:34:42Z, score 45; red_reasons: "release validation FAILED (stage integrate)"). Human authorisation 2026-08-28, verbatim: "prm and then kick off phase 2, I authorize the heart RED thing" — scoped to this task only; release stays human. The RED is pre-existing on main and unrelated: the failing workspace scripts sit on the JAX likelihood path, while this diff touches only the numba imaging-inversion path plus additive-only changes.

   Record: PyAutoMind/complete/2026/08/numba-hst-curvature-matrix-phase2.md.
6. Jammy2211

   closed this as completedon Aug 28, 2026