diff --git a/README.md b/README.md index 358712ef8..e38d0c223 100644 --- a/README.md +++ b/README.md @@ -62,7 +62,7 @@ Results are framed by **astronomy instrument** (HST, Euclid, JWST, …) rather t | `imaging/delaunay` | hst | hpc_a100_mp | dense (mapping) | 96.8 ms | v2026.7.6.649 | | `imaging/delaunay` | hst | hpc_a100_mp | sparse (w-tilde) | 95.5 ms | v2026.7.6.649 | | `imaging/delaunay_numba` | euclid | local_cpu_fp64 | sparse (numba) | 1.19 s | v2026.8.17.1 | -| `imaging/delaunay_numba` | hst | local_cpu_fp64 | sparse (numba) | 679.2 ms | v2026.8.17.1 | +| `imaging/delaunay_numba` | hst | local_cpu_fp64 | sparse (numba) | 482.4 ms | v2026.8.17.1 | | `imaging/mge` | hst | local_cpu_fp64 | dense (mapping) | 179.5 ms | v2026.7.6.649 | | `imaging/mge` | hst | hpc_a100_fp64 | dense (mapping) | 7.8 ms | v2026.7.6.649 | | `imaging/pixelization` | hst | local_cpu_fp64 | dense (mapping) | 8.65 s | v2026.7.6.649 | @@ -71,8 +71,8 @@ Results are framed by **astronomy instrument** (HST, Euclid, JWST, …) rather t | `imaging/pixelization` | hst | hpc_a100_fp64 | sparse (w-tilde) | 57.9 ms | v2026.7.6.649 | | `imaging/pixelization` | hst | hpc_a100_mp | dense (mapping) | 56.4 ms | v2026.7.6.649 | | `imaging/pixelization` | hst | hpc_a100_mp | sparse (w-tilde) | 55.5 ms | v2026.7.6.649 | -| `imaging/pixelization_numba` | euclid | local_cpu_fp64 | sparse (numba) | 288.6 ms | v2026.8.17.1 | -| `imaging/pixelization_numba` | hst | local_cpu_fp64 | sparse (numba) | 691.5 ms | v2026.8.17.1 | +| `imaging/pixelization_numba` | euclid | local_cpu_fp64 | sparse (numba) | 124.3 ms | v2026.8.17.1 | +| `imaging/pixelization_numba` | hst | local_cpu_fp64 | sparse (numba) | 304.5 ms | v2026.8.17.1 | The tables above are auto-generated by `scripts/misc/tooling/build_readme.py` from the artifacts under [`results/`](./results/README.md) — never edit them by hand; run `python scripts/misc/tooling/build_readme.py` after a profiling run and commit the result (CI checks idempotence via `--check`). Narrative context — per-cell "where to optimize next" recommendations and the mp-vs-fp64 verdicts — lives in [`scripts/misc/likelihood_runtime/OPTIMIZATION_NOTES.md`](./scripts/misc/likelihood_runtime/OPTIMIZATION_NOTES.md). diff --git a/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.json b/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.json index 4c02451d3..87a542f95 100644 --- a/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.json +++ b/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.json @@ -25,33 +25,68 @@ "omp_num_threads": "1" }, "steps": { - "FitImaging construct": 7.37200001367455e-05, - "Blurred image (FFT convolve)": 0.010467500000095243, - "Profile subtracted image": 5.2989999994679236e-05, - "Inversion build (trace+Delaunay+mapper)": 0.02722744000016064, - "PSF-weighted data [numba]": 0.008432200000015655, - "MGE operated mapping matrix (60 funcs)": 0.2072877299999618, - "Mapper sparse triplets": 0.018328900000051364, - "Data vector D [numba]": 0.004825830000118003, - "F: mapper\u00d7mapper block [sparse-op]": 0.12305467000014687, - "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.0515464300000076, - "F: linear-func\u00d7linear-func block [BLAS]": 0.01916231999989577, - "Curvature matrix F [residual: mirror + diag-add]": -0.010250559999940373, - "Regularization matrix H (ConstantSplit)": 0.0681292599999324, - "F + H": 0.002749279999989085, - "Reconstruction solve (BLAS)": 0.03926761000002443, - "Mapped reconstruction [numba+FFT]": 0.0034941299999900364, - "Model data + chi^2": 0.003210120000039751, - "log det (F+H) [Cholesky]": 0.02755216999998993, - "log det H [Cholesky]": 0.024576129999877593, - "Regularization term s'Hs": 0.0007812499999545252, - "Log evidence (figure of merit)": 0.0492773199999192 + "FitImaging construct": 7.01600003594649e-05, + "Blurred image (FFT convolve)": 0.009483080000063638, + "Profile subtracted image": 5.181000033189775e-05, + "Inversion build (trace+Delaunay+mapper)": 0.023930310000287135, + "PSF-weighted data [numba]": 0.008010229999672446, + "MGE: image_2d_from(grid) x60 [numpy]": 0.022226080000109505, + "MGE: image_2d_from(blurring_grid) x60 [numpy]": 0.01153678000009677, + "MGE: batched PSF convolution [real-space np]": 0.05661098999953538, + "MGE operated mapping matrix [residual: dict assembly]": 3.100000085395904e-05, + "Mapper sparse triplets": 0.017609500000071422, + "Data vector D [numba]": 0.004609609999715758, + "F: mapper\u00d7mapper block [sparse-op]": 0.06782932999958576, + "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.04873490999980277, + "F: linear-func\u00d7linear-func block [BLAS]": 0.017525710000336404, + "Curvature matrix F [residual: mirror + diag-add]": -0.005954410000049398, + "Regularization matrix H (ConstantSplit)": 0.06016452000058052, + "F + H": 0.001769830000193906, + "Reconstruction solve (BLAS)": 0.0378300099999251, + "Mapped reconstruction [numba+FFT]": 0.0035184800002753037, + "Model data + chi^2": 0.0034454800004823484, + "log det (F+H) [Cholesky]": 0.02494014000021707, + "log det H [Cholesky]": 0.023139239999909477, + "Regularization term s'Hs": 0.0006984200001170394, + "Log evidence (figure of merit)": 0.04456750999979704 + }, + "total_step_by_step": 0.4823787200022708, + "curvature_matrix_f_total_s": 0.12813553999967553, + "mge_operated_mapping_matrix_total_s": 0.09040485000059562, + "mge_operated_mapping_matrix_split_note": "The four 'MGE: ...' / 'MGE operated mapping matrix' rows sum to mge_operated_mapping_matrix_total_s, the raw un-split cost of the inversion.linear_func_operated_mapping_matrix_dict step. The three piece rows re-run, uncached, the three parts of LightProfileLinearObjFuncList.operated_mapping_matrix_override (60 profile images on the data grid, 60 on the blurring grid, one batched real-space PSF convolution); the residual row is the step minus those three. mge_split_reproduces_step records that the three pieces reconstruct the step's output bit-identically. PyAutoArray#507 step 0.", + "mge_split_reproduces_step": true, + "mge_blurring_stack_path": "shared_geometry", + "geometry": { + "psf_shape_native": [ + 21, + 21 + ], + "image_pixels_masked": 15361, + "psf_precision_pairs_stored": 11080573, + "psf_precision_lengths_sum": 11080573, + "psf_precision_lengths_mean": 721.3445088210403, + "psf_precision_lengths_max": 841, + "mappers": [ + { + "params": 1250, + "data_to_pix_unique_shape": [ + 15361, + 3 + ], + "pix_lengths_mean": 1.5523077924614284, + "pix_lengths_max": 3, + "pix_lengths_sum": 23845, + "mapper_mapper_inner_ops": 32691733, + "two_stage_stage_1_scatter_ops": 18187061, + "two_stage_stage_2_axpy_ops": 29806250, + "two_stage_clear_ops": 19201250, + "two_stage_total_ops": 67194561 + } + ] }, - "total_step_by_step": 0.6792464400003608, - "curvature_matrix_f_total_s": 0.18351286000010986, "curvature_matrix_f_split_note": "The four 'F: ...' / 'Curvature matrix F' rows sum to curvature_matrix_f_total_s, the raw un-split cost of the inversion.curvature_matrix step. The three block rows are direct measurements of the (uncached) per-block helpers in imaging_numba/sparse.py; the 'Curvature matrix F [residual: ...]' row is t(curvature_matrix) minus those three, i.e. the global mirror, the no-regularization diagonal add and assembly overhead, and it carries the combined noise of the four timings. PyAutoArray#505 step 0.", - "direct_log_likelihood_function_per_call": 0.7576715100000001, - "warmup_incl_numba_compile_s": 4.814360099999703, + "direct_log_likelihood_function_per_call": 0.4267223799999556, + "warmup_incl_numba_compile_s": 3.621436499999618, "log_likelihood": 29090.527210448134, "pinned_expected": 29090.527192092646, "pinned_drift": [] diff --git a/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.png b/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.png index 36e82c12b..4bf419a46 100644 Binary files a/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.png and b/results/breakdown/imaging/delaunay_numba_breakdown_hst_v2026.8.17.1.png differ diff --git a/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.json b/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.json index 5ad6c2511..22f7561ce 100644 --- a/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.json +++ b/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.json @@ -25,33 +25,68 @@ "omp_num_threads": "1" }, "steps": { - "FitImaging construct": 7.546999995611259e-05, - "Blurred image (FFT convolve)": 0.008001979999926334, - "Profile subtracted image": 3.681000011965807e-05, - "Inversion build (trace+mesh+mapper)": 0.01710581000006641, - "PSF-weighted data [numba]": 0.0026192299999820536, - "MGE operated mapping matrix (60 funcs)": 0.10346844000009696, - "Mapper sparse triplets": 0.0016650700000354845, - "Data vector D [numba]": 0.001205819999859159, - "F: mapper\u00d7mapper block [sparse-op]": 0.06530645000011645, - "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.022554969999964668, - "F: linear-func\u00d7linear-func block [BLAS]": 0.004474150000078225, - "Curvature matrix F [residual: mirror + diag-add]": 0.004506459999856813, - "Regularization matrix H": 0.004038749999972424, - "F + H": 0.0012899899999865739, - "Reconstruction solve (BLAS)": 0.0099663199999668, - "Mapped reconstruction [numba+FFT]": 0.0015609000000040396, - "Model data + chi^2": 0.001559590000033495, - "log det (F+H) [Cholesky]": 0.010603099999934784, - "log det H [Cholesky]": 0.00950673000002098, - "Regularization term s'Hs": 0.0004289500000595581, - "Log evidence (figure of merit)": 0.018604910000021846 + "FitImaging construct": 6.200999996508472e-05, + "Blurred image (FFT convolve)": 0.0056353500001932845, + "Profile subtracted image": 2.655000043887412e-05, + "Inversion build (trace+mesh+mapper)": 0.012572869999530667, + "PSF-weighted data [numba]": 0.002042729999811854, + "MGE: image_2d_from(grid) x60 [numpy]": 0.007493010000325739, + "MGE: image_2d_from(blurring_grid) x60 [numpy]": 0.006900779999705264, + "MGE: batched PSF convolution [real-space np]": 0.015458829999624868, + "MGE operated mapping matrix [residual: dict assembly]": -0.0011304999994536033, + "Mapper sparse triplets": 0.001228639999862935, + "Data vector D [numba]": 0.0007602999996379368, + "F: mapper\u00d7mapper block [sparse-op]": 0.017563919999702193, + "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.014607230000365234, + "F: linear-func\u00d7linear-func block [BLAS]": 0.0033264399999097805, + "Curvature matrix F [residual: mirror + diag-add]": 0.00015849000028538007, + "Regularization matrix H": 0.0016086299998278263, + "F + H": 0.0006751000000804197, + "Reconstruction solve (BLAS)": 0.0067057800002658045, + "Mapped reconstruction [numba+FFT]": 0.001095799999711744, + "Model data + chi^2": 0.0012144199999966077, + "log det (F+H) [Cholesky]": 0.007371300000340852, + "log det H [Cholesky]": 0.006293040000127803, + "Regularization term s'Hs": 0.00035453999989840666, + "Log evidence (figure of merit)": 0.012255779999759398 + }, + "total_step_by_step": 0.12428103999991434, + "curvature_matrix_f_total_s": 0.03565608000026259, + "mge_operated_mapping_matrix_total_s": 0.028722120000202266, + "mge_operated_mapping_matrix_split_note": "The four 'MGE: ...' / 'MGE operated mapping matrix' rows sum to mge_operated_mapping_matrix_total_s, the raw un-split cost of the inversion.linear_func_operated_mapping_matrix_dict step. The three piece rows re-run, uncached, the three parts of LightProfileLinearObjFuncList.operated_mapping_matrix_override (60 profile images on the data grid, 60 on the blurring grid, one batched real-space PSF convolution); the residual row is the step minus those three. mge_split_reproduces_step records that the three pieces reconstruct the step's output bit-identically. PyAutoArray#507 step 0.", + "mge_split_reproduces_step": true, + "mge_blurring_stack_path": "shared_geometry", + "geometry": { + "psf_shape_native": [ + 21, + 21 + ], + "image_pixels_masked": 3841, + "psf_precision_pairs_stored": 2318869, + "psf_precision_lengths_sum": 2318869, + "psf_precision_lengths_mean": 603.7149179901068, + "psf_precision_lengths_max": 841, + "mappers": [ + { + "params": 784, + "data_to_pix_unique_shape": [ + 3841, + 4 + ], + "pix_lengths_mean": 4.0, + "pix_lengths_max": 4, + "pix_lengths_sum": 15364, + "mapper_mapper_inner_ops": 37101904, + "two_stage_stage_1_scatter_ops": 9275476, + "two_stage_stage_2_axpy_ops": 12045376, + "two_stage_clear_ops": 3011344, + "two_stage_total_ops": 24332196 + } + ] }, - "total_step_by_step": 0.28857990000005884, - "curvature_matrix_f_total_s": 0.09684203000001616, "curvature_matrix_f_split_note": "The four 'F: ...' / 'Curvature matrix F' rows sum to curvature_matrix_f_total_s, the raw un-split cost of the inversion.curvature_matrix step. The three block rows are direct measurements of the (uncached) per-block helpers in imaging_numba/sparse.py; the 'Curvature matrix F [residual: ...]' row is t(curvature_matrix) minus those three, i.e. the global mirror, the no-regularization diagonal add and assembly overhead, and it carries the combined noise of the four timings. PyAutoArray#505 step 0.", - "direct_log_likelihood_function_per_call": 0.24940190999996048, - "warmup_incl_numba_compile_s": 2.1857472000001508, + "direct_log_likelihood_function_per_call": 0.11088282000000618, + "warmup_incl_numba_compile_s": 1.5452468999992561, "log_likelihood": 6213.306873885871, "pinned_expected": null, "pinned_drift": [] diff --git a/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.png b/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.png index 4e7f698d1..b27cdd6e8 100644 Binary files a/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.png and b/results/breakdown/imaging/pixelization_numba_breakdown_euclid_v2026.8.17.1.png differ diff --git a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.json b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.json index 3ebd3341a..6bbea8a22 100644 --- a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.json +++ b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.json @@ -25,33 +25,68 @@ "omp_num_threads": "1" }, "steps": { - "FitImaging construct": 8.33399999919493e-05, - "Blurred image (FFT convolve)": 0.01084737999999561, - "Profile subtracted image": 6.784000001971435e-05, - "Inversion build (trace+mesh+mapper)": 0.026095339999892532, - "PSF-weighted data [numba]": 0.008763999999928274, - "MGE operated mapping matrix (60 funcs)": 0.22401575999997475, - "Mapper sparse triplets": 0.005336110000007466, - "Data vector D [numba]": 0.0052595000000565054, - "F: mapper\u00d7mapper block [sparse-op]": 0.27679772000001324, - "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.05443964000005508, - "F: linear-func\u00d7linear-func block [BLAS]": 0.021487699999988764, - "Curvature matrix F [residual: mirror + diag-add]": 0.0066316599999936, - "Regularization matrix H": 0.0018421100000068692, - "F + H": 0.0008287500000278669, - "Reconstruction solve (BLAS)": 0.00876641000008931, - "Mapped reconstruction [numba+FFT]": 0.003643519999968703, - "Model data + chi^2": 0.0037507700000787734, - "log det (F+H) [Cholesky]": 0.008669489999920187, - "log det H [Cholesky]": 0.007879109999976208, - "Regularization term s'Hs": 0.00032316000001628706, - "Log evidence (figure of merit)": 0.015974480000113544 + "FitImaging construct": 6.414999952539802e-05, + "Blurred image (FFT convolve)": 0.00926681000037206, + "Profile subtracted image": 4.731000008177944e-05, + "Inversion build (trace+mesh+mapper)": 0.021709480000026815, + "PSF-weighted data [numba]": 0.007878119999440969, + "MGE: image_2d_from(grid) x60 [numpy]": 0.02258933999983128, + "MGE: image_2d_from(blurring_grid) x60 [numpy]": 0.010614800000075775, + "MGE: batched PSF convolution [real-space np]": 0.046522589999585764, + "MGE operated mapping matrix [residual: dict assembly]": -0.003177729998969883, + "Mapper sparse triplets": 0.004005030000007537, + "Data vector D [numba]": 0.003737139999975625, + "F: mapper\u00d7mapper block [sparse-op]": 0.0751765200000591, + "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.045131479999690784, + "F: linear-func\u00d7linear-func block [BLAS]": 0.017325829999936104, + "Curvature matrix F [residual: mirror + diag-add]": -0.0015314400001443518, + "Regularization matrix H": 0.0013350499997613952, + "F + H": 0.0006081399998947746, + "Reconstruction solve (BLAS)": 0.007613869999840972, + "Mapped reconstruction [numba+FFT]": 0.002937240000028396, + "Model data + chi^2": 0.003083939999669383, + "log det (F+H) [Cholesky]": 0.007987400000092747, + "log det H [Cholesky]": 0.007959590000427851, + "Regularization term s'Hs": 0.00029186999981902775, + "Log evidence (figure of merit)": 0.01328899000036472 + }, + "total_step_by_step": 0.30446551999939403, + "curvature_matrix_f_total_s": 0.13610238999954163, + "mge_operated_mapping_matrix_total_s": 0.07654900000052293, + "mge_operated_mapping_matrix_split_note": "The four 'MGE: ...' / 'MGE operated mapping matrix' rows sum to mge_operated_mapping_matrix_total_s, the raw un-split cost of the inversion.linear_func_operated_mapping_matrix_dict step. The three piece rows re-run, uncached, the three parts of LightProfileLinearObjFuncList.operated_mapping_matrix_override (60 profile images on the data grid, 60 on the blurring grid, one batched real-space PSF convolution); the residual row is the step minus those three. mge_split_reproduces_step records that the three pieces reconstruct the step's output bit-identically. PyAutoArray#507 step 0.", + "mge_split_reproduces_step": true, + "mge_blurring_stack_path": "shared_geometry", + "geometry": { + "psf_shape_native": [ + 21, + 21 + ], + "image_pixels_masked": 15361, + "psf_precision_pairs_stored": 11080573, + "psf_precision_lengths_sum": 11080573, + "psf_precision_lengths_mean": 721.3445088210403, + "psf_precision_lengths_max": 841, + "mappers": [ + { + "params": 784, + "data_to_pix_unique_shape": [ + 15361, + 4 + ], + "pix_lengths_mean": 4.0, + "pix_lengths_max": 4, + "pix_lengths_sum": 61444, + "mapper_mapper_inner_ops": 177289168, + "two_stage_stage_1_scatter_ops": 44322292, + "two_stage_stage_2_axpy_ops": 48172096, + "two_stage_clear_ops": 12043024, + "two_stage_total_ops": 104537412 + } + ] }, - "total_step_by_step": 0.6915037900001153, - "curvature_matrix_f_total_s": 0.3593567200000507, "curvature_matrix_f_split_note": "The four 'F: ...' / 'Curvature matrix F' rows sum to curvature_matrix_f_total_s, the raw un-split cost of the inversion.curvature_matrix step. The three block rows are direct measurements of the (uncached) per-block helpers in imaging_numba/sparse.py; the 'Curvature matrix F [residual: ...]' row is t(curvature_matrix) minus those three, i.e. the global mirror, the no-regularization diagonal add and assembly overhead, and it carries the combined noise of the four timings. PyAutoArray#505 step 0.", - "direct_log_likelihood_function_per_call": 0.5949823200000083, - "warmup_incl_numba_compile_s": 4.023859600000378, + "direct_log_likelihood_function_per_call": 0.2970425999999861, + "warmup_incl_numba_compile_s": 3.0814270000009856, "log_likelihood": 27661.91013366411, "pinned_expected": 27661.910133665442, "pinned_drift": [] diff --git a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.png b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.png index 2e974322a..40b4eeab6 100644 Binary files a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.png and b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1.png differ diff --git a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.json b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.json index f1727d9d3..a6e72ee2b 100644 --- a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.json +++ b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.json @@ -25,33 +25,68 @@ "omp_num_threads": "1" }, "steps": { - "FitImaging construct": 6.170000001475273e-05, - "Blurred image (FFT convolve)": 0.01020472999980484, - "Profile subtracted image": 5.66400000025169e-05, - "Inversion build (trace+mesh+mapper)": 0.024389229999997042, - "PSF-weighted data [numba]": 0.00845435000001089, - "MGE operated mapping matrix (60 funcs)": 0.20154515999997785, - "Mapper sparse triplets": 7.080621299999985, - "Data vector D [numba]": 0.00409297999995033, - "F: mapper\u00d7mapper block [sparse-op]": 0.24928207000002658, - "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.06102243000013914, - "F: linear-func\u00d7linear-func block [BLAS]": 0.019332219999932933, - "Curvature matrix F [residual: mirror + diag-add]": 0.004671089999874256, - "Regularization matrix H": 0.0016499499999554246, - "F + H": 0.0008134300001074735, - "Reconstruction solve (BLAS)": 0.010040690000050746, - "Mapped reconstruction [numba+FFT]": 0.0035483700000440875, - "Model data + chi^2": 0.0035280499999771563, - "log det (F+H) [Cholesky]": 0.008806140000069717, - "log det H [Cholesky]": 0.007731979999971372, - "Regularization term s'Hs": 0.00030322000002342976, - "Log evidence (figure of merit)": 0.014850420000084342 + "FitImaging construct": 6.83200003550155e-05, + "Blurred image (FFT convolve)": 0.01055562000037753, + "Profile subtracted image": 5.696000043826643e-05, + "Inversion build (trace+mesh+mapper)": 0.02428079000001162, + "PSF-weighted data [numba]": 0.00830354999998235, + "MGE: image_2d_from(grid) x60 [numpy]": 0.023540689999936148, + "MGE: image_2d_from(blurring_grid) x60 [numpy]": 0.012065089999668998, + "MGE: batched PSF convolution [real-space np]": 0.05100648000006913, + "MGE operated mapping matrix [residual: dict assembly]": 0.0073249800003395565, + "Mapper sparse triplets": 6.733646810000209, + "Data vector D [numba]": 0.004568759999710892, + "F: mapper\u00d7mapper block [sparse-op]": 0.08872232000012445, + "F: mapper\u00d7linear-func block [FFT conv + scatter]": 0.05873444999961066, + "F: linear-func\u00d7linear-func block [BLAS]": 0.02017093000049499, + "Curvature matrix F [residual: mirror + diag-add]": -0.016193320000274986, + "Regularization matrix H": 0.0016231399999014685, + "F + H": 0.0007177499999670544, + "Reconstruction solve (BLAS)": 0.008862000000408442, + "Mapped reconstruction [numba+FFT]": 0.0033353099999658296, + "Model data + chi^2": 0.003341770000224642, + "log det (F+H) [Cholesky]": 0.008235310000418394, + "log det H [Cholesky]": 0.007438530000035826, + "Regularization term s'Hs": 0.00031590000016876727, + "Log evidence (figure of merit)": 0.014587909999863768 + }, + "total_step_by_step": 7.075310050002008, + "curvature_matrix_f_total_s": 0.15143437999995513, + "mge_operated_mapping_matrix_total_s": 0.09393724000001384, + "mge_operated_mapping_matrix_split_note": "The four 'MGE: ...' / 'MGE operated mapping matrix' rows sum to mge_operated_mapping_matrix_total_s, the raw un-split cost of the inversion.linear_func_operated_mapping_matrix_dict step. The three piece rows re-run, uncached, the three parts of LightProfileLinearObjFuncList.operated_mapping_matrix_override (60 profile images on the data grid, 60 on the blurring grid, one batched real-space PSF convolution); the residual row is the step minus those three. mge_split_reproduces_step records that the three pieces reconstruct the step's output bit-identically. PyAutoArray#507 step 0.", + "mge_split_reproduces_step": true, + "mge_blurring_stack_path": "shared_geometry", + "geometry": { + "psf_shape_native": [ + 21, + 21 + ], + "image_pixels_masked": 15361, + "psf_precision_pairs_stored": 11080573, + "psf_precision_lengths_sum": 11080573, + "psf_precision_lengths_mean": 721.3445088210403, + "psf_precision_lengths_max": 841, + "mappers": [ + { + "params": 784, + "data_to_pix_unique_shape": [ + 15361, + 4 + ], + "pix_lengths_mean": 4.0, + "pix_lengths_max": 4, + "pix_lengths_sum": 61444, + "mapper_mapper_inner_ops": 177289168, + "two_stage_stage_1_scatter_ops": 44322292, + "two_stage_stage_2_axpy_ops": 48172096, + "two_stage_clear_ops": 12043024, + "two_stage_total_ops": 104537412 + } + ] }, - "total_step_by_step": 7.715006150000001, - "curvature_matrix_f_total_s": 0.33430780999997295, "curvature_matrix_f_split_note": "The four 'F: ...' / 'Curvature matrix F' rows sum to curvature_matrix_f_total_s, the raw un-split cost of the inversion.curvature_matrix step. The three block rows are direct measurements of the (uncached) per-block helpers in imaging_numba/sparse.py; the 'Curvature matrix F [residual: ...]' row is t(curvature_matrix) minus those three, i.e. the global mirror, the no-regularization diagonal add and assembly overhead, and it carries the combined noise of the four timings. PyAutoArray#505 step 0.", - "direct_log_likelihood_function_per_call": 7.479174180000018, - "warmup_incl_numba_compile_s": 11.205480999999963, + "direct_log_likelihood_function_per_call": 6.6226883100000125, + "warmup_incl_numba_compile_s": 9.757296899999346, "log_likelihood": 27180.70471569685, "pinned_expected": 27180.704715698186, "pinned_drift": [] diff --git a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.png b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.png index 7e42d15f9..a35dc1975 100644 Binary files a/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.png and b/results/breakdown/imaging/pixelization_numba_breakdown_hst_v2026.8.17.1_rtu.png differ diff --git a/results/notes/numba_curvature_matrix_f_split.md b/results/notes/numba_curvature_matrix_f_split.md index 4c79c5f6c..7fc6f90bb 100644 --- a/results/notes/numba_curvature_matrix_f_split.md +++ b/results/notes/numba_curvature_matrix_f_split.md @@ -247,3 +247,256 @@ This matches the code: `Convolver`'s FFT path uses `np.fft.rfft2`, which has no dataset is a 2828-pixel 0.1"/pixel simulation with a small rectangular mesh, so F is a much smaller share of the evaluation than in the HST harness cells and the end-to-end gain is correspondingly smaller than the 2.6x measured there. + + +--- + +# Phase 2 (#507): the mapper x mapper block + the MGE operated matrix + +Measured 2026-08-28 on the same local WSL host and settings as phase 1 +(`OMP_NUM_THREADS=1`, `AUTOARRAY_NUMBA_OPERATED_MEMO=0`, `n_repeats = 10`, +autolens `v2026.8.17.1`), branch `feature/numba-hst-curvature-matrix-phase2` +across PyAutoArray, PyAutoGalaxy and this repo. + +Phase 1 left the HST rectangular evaluation at 0.60 s with two terms carrying +it: the mapper x mapper sparse-operator block (46 %) and the MGE operated +mapping matrix (37 %). Phase 2 is those two. + +## Step 0 — the split, measured before anything was changed + +Two additions to both numba CPU breakdown harnesses (commit `3f69b34`): the +`MGE operated mapping matrix` row split into its three pieces plus a residual +(the same residual scheme the phase-1 F rows use, so the four rows sum to the +un-split step, recorded as `mge_operated_mapping_matrix_total_s`), and the +model-independent geometry constants written once into the result JSON. + +| cell | step total | mapper x mapper | MGE total | MGE = grid eval + blurring eval + convolution | +|---|---|---|---|---| +| hst, rectangular bilinear | 0.633 | **0.2551** | **0.2086** | 0.1160 + 0.0460 + 0.0541 (residual -0.0075) | +| euclid, rectangular bilinear | 0.198 | 0.0491 | 0.0704 | | +| hst, Delaunay-1250 | 0.734 | 0.1304 | 0.2229 | | + +**Profile evaluation is 78 % of the MGE row** at HST — not the convolution, +which phase 1's batching had already reduced to 0.054 s. And the geometry +constants turned the plan's "~2e8 inner operations" estimate into a +measurement: `sum over stored (data_0, data_1) pairs of u0 * u1` = **1.773e8** +at HST (11.08e6 stored pairs, mean row length 721, `u0 = u1 = 4.0` for the +bilinear mapper). + +That checkpoint fixed the two levers: restructure the mapper x mapper kernel, +and hoist the geometry the 60 Gaussians share. + +## What changed in the libraries + +**Step 1 (PyAutoArray `a583f1a6`) — oracle first, then hoist.** The kernel had +no direct unit test; it now has one pinned to `F = M.T @ W @ M` with a dense +`W`, plus a symmetry / halved-diagonal test, and the pre-change quadruple loop +is retained as `curvature_matrix_via_sparse_operator_reference_from`. The hoist +itself takes 1-D row views once per data pixel instead of re-gathering from +wide-stride 2-D arrays u0 times, leaving the accumulated expression +operand-for-operand identical. Bit-identical (`np.array_equal` on both the +kernel output and the assembled `curvature_matrix`), 1.04-1.12x. + +**Step 2 (PyAutoArray `88e14bc6`) — the two-stage reformulation.** `w0` does +not depend on `data_1`, so the sum factorises into a per-data-pixel dense +accumulator over source space followed by contiguous AXPYs over whole rows of +`F`. Replaces ~1.8e8 irregular read-modify-writes into a 4.9 MB matrix with +~4.4e7 L1 scatters plus vectorisable dense adds. Agrees with the direct form to +floating-point reassociation, not bit-identically. + +**Step 3a (PyAutoArray `d8bc3bac`) — the OverSampler divisor.** +`binned_array_2d_from` recomputed `np.bincount(segment_ids)` *and* the +`sub_is_uniform` check on every call — 120 calls per evaluation for a +60-Gaussian MGE. Both depend only on the constructor's `sub_size`, and are now +`cached_property`. On the HST over sampler (17980 sub-pixels into 15361 pixels) +the divisor cost 38.8 us and the uniformity check 61.3 us per call, against +81.5 us for the whole cached call. Bit-identical; the cached divisor is +read-only, because the old code guarded zeros by mutating in place. + +**Step 3b (PyAutoGalaxy `6031aa0a`) — one eccentric-radius grid per MGE +basis.** 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: the transform alone was 951 us of the 1.3 ms +each profile took, ~70 % of the MGE's profile-evaluation cost. +`LightProfileLinearObjFuncList._image_slim_list_from` now groups the profiles by +`(class, centre, ell_comps)` and computes each group's transform and radii once. +Bit-identical. + +Grouping rather than an all-or-nothing check is deliberate: the workspace's +canonical MGE recipe stacks two sets of 30 Gaussians which share a centre but +carry their own `ell_comps`, and `GalaxiesToInversion` lands both sets in a +single func list — an all-or-nothing test would fall back to 60 independent +evaluations on exactly the model users are told to write. + +The membership test is the identity of the `image_2d_from` function object +against `Gaussian.image_2d_from`, not `isinstance`. `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 into it would be +both slower and a different floating-point expression. Giving `Gaussian` that +same Cartesian shortcut is a larger, value-changing lever, and is left as a +follow-up. + +## The measured crossover (step 2) + +`pix_pixels` swept on the two production HST geometries — real 15361-pixel +sparse operator, real mapper mappings, only the source-space extent varied. +Speed-up of two-stage over direct: + +``` +pix_pixels 128 784 1250 2048 4096 6144 8192 +rectangular u0=4.00 3.10 2.93 2.47 2.10 1.36 1.12 1.04 +delaunay u0=1.55 1.79 1.72 1.53 1.33 1.04 1.01 0.98 +``` + +The crossover is geometry-dependent and sits at or beyond 8192 source pixels — +outside the range a PyAuto pixelization is run at. Two-stage also wins on the +Delaunay cell (1.53x at its native 1250) despite executing 2.05x *more* +arithmetic there, because a vectorised contiguous AXPY is several times cheaper +per operation than a scattered read-modify-write into a 12.5 MB matrix. The +guard `CURVATURE_TWO_STAGE_MAX_PIX_PIXELS = 4096` is therefore the conservative +end of the measured envelope, an explicit module-level constant carrying this +sweep, not a fitted crossover and not a silent heuristic. + +## Before / after (seconds per evaluation) + +Three arms measured **paired in one session, alternating B/A/B/A**, on this +host's 20-30 % run-to-run variance: `base` = PyAutoArray `main` @ `1b89404b` +with PyAutoGalaxy `main` @ `0fbe863d` (the phase-2 starting point, i.e. phase 1 +as merged); `step2` = after PyAutoArray steps 1-2; `final` = after step 3 in +both libraries. Three rounds were run per cell in rotating arm order; the first +round is discarded (numba `cache: true` recompiles on an arm switch) and the +table is the mean of rounds 2 and 3. + +**hst, rectangular bilinear (784 source pixels)** + +| row | base | step 2 | final | step 3 | whole phase | +|---|---|---|---|---|---| +| step total | 0.6214 | 0.4324 | **0.3334** | 1.30x | **1.86x** | +| direct `log_likelihood_function` | 0.6184 | 0.4207 | **0.3013** | 1.40x | **2.05x** | +| F total | 0.3246 | 0.1441 | 0.1494 | 0.96x | 2.17x | +| F: mapper x mapper | 0.2581 | 0.0819 | 0.0825 | 0.99x | **3.13x** | +| MGE total | 0.1970 | 0.1937 | **0.0844** | **2.29x** | 2.33x | +| MGE: 60x image_2d_from(grid) | 0.1042 | 0.0995 | 0.0233 | **4.28x** | 4.48x | +| MGE: 60x image_2d_from(blurring) | 0.0436 | 0.0420 | 0.0115 | 3.67x | 3.81x | +| MGE: batched PSF convolution | 0.0512 | 0.0504 | 0.0523 | 0.96x | 0.98x | +| Blurred image (FFT convolve) | 0.0101 | 0.0097 | 0.0099 | 0.98x | 1.02x | + +**euclid, rectangular bilinear (784 source pixels)** + +| row | base | step 2 | final | step 3 | whole phase | +|---|---|---|---|---|---| +| step total | 0.2226 | 0.1870 | **0.1417** | 1.32x | **1.57x** | +| direct `log_likelihood_function` | 0.2076 | 0.1644 | 0.1224 | 1.34x | 1.70x | +| F: mapper x mapper | 0.0531 | 0.0188 | 0.0191 | 0.98x | 2.78x | +| MGE total | 0.0802 | 0.0797 | 0.0316 | 2.52x | 2.54x | +| MGE: 60x image_2d_from(grid) | 0.0353 | 0.0353 | 0.0081 | 4.37x | 4.36x | +| MGE: 60x image_2d_from(blurring) | 0.0283 | 0.0280 | 0.0075 | 3.75x | 3.79x | + +**hst, Delaunay-1250** + +| row | base | step 2 | final | step 3 | whole phase | +|---|---|---|---|---|---| +| step total | 0.6331 | 0.6076 | **0.4884** | 1.24x | **1.30x** | +| direct `log_likelihood_function` | 0.5785 | 0.5640 | 0.4563 | 1.24x | 1.27x | +| F: mapper x mapper | 0.1088 | 0.0697 | 0.0676 | 1.03x | 1.61x | +| MGE total | 0.1881 | 0.1932 | 0.0839 | 2.30x | 2.24x | +| MGE: 60x image_2d_from(grid) | 0.0971 | 0.1015 | 0.0224 | 4.53x | 4.33x | + +**The issue's goal — 0.60 s -> ~0.35 s at HST rectangular — is met**: 0.6214 -> +0.3334 s step total, 0.6184 -> 0.3013 s directly timed. The mapper x mapper +block and the MGE profile evaluation, the two terms step 0 identified, are down +3.1x and 4.3x respectively. Nothing outside those two moves, which is the +control the unrelated rows above provide (`Blurred image` 0.98-1.02x, the +batched PSF convolution 0.96-1.00x). + +The Delaunay cell gains least (1.30x) because its evaluation is spread across +more terms: after phase 2 its largest steps are the regularization matrix H +(ConstantSplit, 0.060 s), the batched MGE convolution (0.057 s) and the +mapper x linear-func block (0.049 s), none of which this issue touched. + +### What was not paired + +`hst, rectangular RTU` is an 7.1 s cell dominated by its 6.73 s "mapper sparse +triplets" step and is GPU-only by the 2026-08-28 decision, so it was re-run once +on `final` for currency (`_rtu` artifact) rather than paired. Its pinned check +passes and its MGE row moves with the others. + +The harness's `MGE: image_2d_from(blurring_grid) x60` row was changed in this +phase (commit below): it hand-rolled the per-profile loop, which after step 3b +would have kept timing the *old* code on the new library and mis-attributed the +win to the residual row. It now calls the same path +`operated_mapping_matrix_override` calls, dispatching on the attribute so both +libraries are measurable, and records which branch it took in the result JSON as +`mge_blurring_stack_path` (`per_profile_loop` on `base`/`step2`, +`shared_geometry` on `final`, as the artifacts show). + +## Correctness + +Pinned log-likelihood checks (explicit `rtol=1e-6`) **PASSED** on all three +pinned cells, on every arm, to every recorded digit: + +| cell | pinned | measured (base / step2 / final) | +|---|---|---| +| hst bilinear | 27661.910133665442 | 27661.91013366411 on all three | +| hst RTU | 27180.704715698186 | 27180.70471569685 (final, unpaired) | +| hst Delaunay | 29090.527192092646 | 29090.527210448134 on all three | + +euclid has no pinned value by design and measured 6213.306873885871 on all +three arms — unchanged to every recorded digit. + +Array-level, on `inversion.curvature_matrix` and on the MGE operated mapping +matrix, computed from the same instance on each arm: + +| comparison | MGE operated matrix | curvature matrix F | +|---|---|---| +| step 2 -> final (step 3 alone), hst | bit-identical | **bit-identical** | +| step 2 -> final (step 3 alone), euclid | bit-identical | **bit-identical** | +| base -> final (whole phase), hst | bit-identical | max rel 1.24e-14, `allclose(rtol=1e-12)` True | +| base -> final (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 comes from step 2's reassociation, +and it is twelve orders inside the pins. + +## Pool run (multiprocessing oversubscription check) + +Same purpose and setup as phase 1's: the harness measures one process at +`OMP_NUM_THREADS=1`, but Nautilus runs one process per core, so a single-thread +win only counts if it survives the real pool. + +**Setup.** `autolens_workspace/scripts/imaging/features/pixelization/cpu_fast_modeling.py`, +first (non-SLaM) fit, reproduced from the workspace CWD on an 8-core host, +against `base` (PyAutoArray `1b89404b` + PyAutoGalaxy `0fbe863d`) and `final`, +with `autoarray.__file__` asserted on both. `PYAUTO_TEST_MODE=1` (Nautilus +`n_like_max = 1`, one exploration batch of 100 evaluations, so per-evaluation +time is wall/100), `PYAUTO_SMALL_DATASETS` unset (2828 masked image-pixels, the +simulator's full resolution), `number_of_cores` raised 2 -> 8. One deviation +from phase 1's config B: the lens light is +`al.model_util.mge_model_from(total_gaussians=60)` — a single 60-Gaussian basis +— rather than phase 1's 30 x 2 recipe, because the two-basis form gives its two +halves different `ell_comps` and this run needs the step-3b path exercised. +(Grouping means the 30 x 2 form takes it too, as two groups.) + +Per evaluation (s), alternating arms: + +| | run 1 | run 2 | run 3 | median | +|---|---|---|---|---| +| base, pool of 8 | 0.2115 | 0.1918 | 0.2000 | 0.2000 | +| final, pool of 8 | 0.1606 | 0.1816 | 0.1697 | 0.1697 | +| base, serial (`number_of_cores=1`) | 0.4099 | | | 0.4099 | +| final, serial (`number_of_cores=1`) | 0.3634 | | | 0.3634 | + +**Verdict — no oversubscription.** The pool improves by 15 % and the single +process by 11 %, so the pool gain tracks (slightly exceeds) the single-thread +gain rather than eroding it. The parallel speed-up ratio is flat to slightly up +across the change — 0.4099/0.2000 = **2.05x** on 8 cores before, 0.3634/0.1697 += **2.14x** after — which is the number that would fall if any of the three +levers had introduced hidden threads. None of them can: the two-stage kernel is +a single-threaded numba `njit` with no `prange`, and both step-3 changes remove +work rather than add machinery. + +**Read this as a regression check, not a speed-up measurement**, for the same +reason as phase 1: the workspace dataset is a 2828-pixel 0.1"/pixel simulation +with a small rectangular mesh, so both the F block and the MGE are a much +smaller share of the evaluation than in the HST harness cells. diff --git a/scripts/imaging/likelihood_breakdown/delaunay_numba.py b/scripts/imaging/likelihood_breakdown/delaunay_numba.py index 83cfa93c0..2760a3608 100644 --- a/scripts/imaging/likelihood_breakdown/delaunay_numba.py +++ b/scripts/imaging/likelihood_breakdown/delaunay_numba.py @@ -185,6 +185,7 @@ def _profiling_root() -> _Path: from autoarray.inversion.inversion.imaging_numba.sparse import ( # noqa: E402 InversionImagingSparseNumba, ) +from autoarray.inversion.mappers.abstract import Mapper # noqa: E402 n_image_pixels = dataset.data.shape[0] n_over_sampled_pixels = dataset.grids.lp.over_sampled.shape[0] @@ -259,6 +260,190 @@ def _f_block_scratch(inversion) -> np.ndarray: return scratch +# ------------------------------------------------------------------- +# MGE operated mapping matrix sub-block instrumentation (PyAutoArray#507 step 0) +# ------------------------------------------------------------------- +# +# ``inversion.linear_func_operated_mapping_matrix_dict`` builds, for the linear +# MGE lens light, the PSF-convolved image of every one of its 60 Gaussians. On +# the numpy path that is +# ``LightProfileLinearObjFuncList.operated_mapping_matrix_override`` +# (autogalaxy/profiles/light/linear/abstract.py), whose fast branch is three +# distinct pieces of work: +# +# 1. ``mapping_matrix`` -- 60x ``image_2d_from(grid)``, stacked +# 2. the blurring stack -- 60x ``image_2d_from(blurring_grid)``, stacked +# 3. one batched real-space PSF convolution of the two stacks +# +# The override is a ``cached_property`` on the linear-func object and the three +# pieces below are recomputed from scratch (``mapping_matrix`` is a plain +# ``property``), so timing them does *not* prime the step that follows. As with +# the F sub-blocks, the pre-existing row is therefore reported as a RESIDUAL, +# +# MGE residual = t(linear_func_operated_mapping_matrix_dict) - (t1 + t2 + t3) +# +# so the four MGE rows sum to exactly the un-instrumented step and the artifact +# stays comparable with pre-instrumentation runs. The raw (unsplit) MGE total is +# written to the result JSON as ``mge_operated_mapping_matrix_total_s``. +# +# ``mge_split_reproduces_step`` in the JSON records a one-off bit-identical +# check that the three timed pieces really do reconstruct the dict the step +# returns -- if a dataset ever takes a different branch of the override the +# split is flagged rather than silently mis-attributed. +# +# Caveat (as for F): the residual is a difference of averaged timings and can in +# principle come out slightly negative. It is recorded as measured, never clipped. + +MGE_PROFILE_IMAGE_LABEL = "MGE: image_2d_from(grid) x60 [numpy]" +MGE_BLURRING_IMAGE_LABEL = "MGE: image_2d_from(blurring_grid) x60 [numpy]" +MGE_CONVOLVE_LABEL = "MGE: batched PSF convolution [real-space np]" +MGE_RESIDUAL_LABEL = "MGE operated mapping matrix [residual: dict assembly]" + +_mge_scratch: dict[str, list] = {} + + +def _mge_linear_func_list(inversion) -> list: + """The inversion's linear-func objects (the MGE lens light), never its mappers.""" + return [obj for obj in inversion.linear_obj_list if hasattr(obj, "light_profile_list")] + + +def _mge_mapping_matrix(fit) -> list: + """Piece 1: the 60 unblurred profile images on the (over-sampled) data grid.""" + matrices = [linear_func.mapping_matrix for linear_func in _mge_linear_func_list(fit.inversion)] + _mge_scratch["mapping_matrix"] = matrices + return matrices + + +def _mge_blurring_stack(linear_func): + """The blurring stack of one linear-func object, built the way the override builds it. + + ``operated_mapping_matrix_override`` stacked the 60 blurring-grid images with an + inline per-profile loop until PyAutoArray#507 step 3, which moved both stacks behind + ``LightProfileLinearObjFuncList._image_slim_list_from`` so that an MGE basis shares one + reference-frame transform and one eccentric-radius grid. This harness is run against + both libraries (the paired before/after of that step), so it dispatches on the + attribute rather than hard-coding either shape -- a hand-rolled loop would keep timing + the *old* code on the new library and mis-attribute the win to the residual row. Which + branch was taken is recorded in the result JSON as ``mge_blurring_stack_path``. + """ + if hasattr(linear_func, "_image_slim_list_from"): + _mge_scratch["blurring_stack_path"] = "shared_geometry" + + return np.stack( + linear_func._image_slim_list_from(grid=linear_func.blurring_grid, xp=np), + axis=1, + ) + + _mge_scratch["blurring_stack_path"] = "per_profile_loop" + + return np.stack( + [ + light_profile.image_2d_from(grid=linear_func.blurring_grid, xp=np).slim.array + for light_profile in linear_func.light_profile_list + ], + axis=1, + ) + + +def _mge_blurring_mapping_matrix(fit) -> list: + """Piece 2: the same 60 profiles on the blurring grid (flux blurred in from + outside the mask), stacked exactly as the override stacks them.""" + matrices = [ + _mge_blurring_stack(linear_func) for linear_func in _mge_linear_func_list(fit.inversion) + ] + _mge_scratch["blurring_mapping_matrix"] = matrices + return matrices + + +def _mge_convolved(fit) -> list: + """Piece 3: the single batched real-space convolution of the two stacks.""" + return [ + linear_func.psf.convolved_mapping_matrix_via_real_space_np_from( + mapping_matrix=_mge_scratch["mapping_matrix"][index], + mask=linear_func.grid.mask, + blurring_mapping_matrix=_mge_scratch["blurring_mapping_matrix"][index], + blurring_mask=linear_func.blurring_grid.mask, + ) + for index, linear_func in enumerate(_mge_linear_func_list(fit.inversion)) + ] + + +def _mge_split_reproduces_step(fit) -> bool: + """Bit-identical check that pieces 1-3 reconstruct the timed step's output.""" + _mge_mapping_matrix(fit) + _mge_blurring_mapping_matrix(fit) + reconstructed = _mge_convolved(fit) + + actual = list(fit.inversion.linear_func_operated_mapping_matrix_dict.values()) + + if len(actual) != len(reconstructed) or len(actual) == 0: + return False + + return all(np.array_equal(np.asarray(a), np.asarray(b)) for a, b in zip(actual, reconstructed)) + + +def _geometry_constants(fit) -> dict: + """Model-independent geometry of the sparse operator and of each mapper. + + Recorded once, outside the timed loop, so that the complexity of the + mapper x mapper kernel (`curvature_matrix_via_sparse_operator_from`) is a + measurement rather than an estimate. `mapper_mapper_inner_ops` is the exact + count of that kernel's innermost accumulations, + + sum over stored (data_0, data_1) pairs of u0(data_0) * u1(data_1), + + and the `two_stage_*` counts are the same geometry costed for the phase-2 + reformulation (a dense n_source accumulator per data pixel, then one AXPY + per mapping). PyAutoArray#507 step 0. + """ + inversion = fit.inversion + + lengths = np.asarray(inversion.sparse_operator.lengths).astype("int64") + indexes = np.asarray(inversion.sparse_operator.indexes).astype("int64") + + geometry = { + "psf_shape_native": [int(s) for s in dataset.psf.kernel.shape_native], + "image_pixels_masked": int(n_image_pixels), + "psf_precision_pairs_stored": int(indexes.shape[0]), + "psf_precision_lengths_sum": int(lengths.sum()), + "psf_precision_lengths_mean": float(lengths.mean()), + "psf_precision_lengths_max": int(lengths.max()), + "mappers": [], + } + + for mapper in inversion.cls_list_from(cls=Mapper): + pix_lengths = np.asarray(mapper.unique_mappings.pix_lengths).astype("int64") + + pair_u0 = np.repeat(pix_lengths, lengths) + pair_u1 = pix_lengths[indexes] + + n_source = int(mapper.params) + data_pixels = int(lengths.shape[0]) + + stage_1_ops = int(pair_u1.sum()) + stage_2_ops = int(pix_lengths.sum()) * n_source + clear_ops = data_pixels * n_source + + geometry["mappers"].append( + { + "params": n_source, + "data_to_pix_unique_shape": [ + int(s) for s in np.asarray(mapper.unique_mappings.data_to_pix_unique).shape + ], + "pix_lengths_mean": float(pix_lengths.mean()), + "pix_lengths_max": int(pix_lengths.max()), + "pix_lengths_sum": int(pix_lengths.sum()), + "mapper_mapper_inner_ops": int((pair_u0 * pair_u1).sum()), + "two_stage_stage_1_scatter_ops": stage_1_ops, + "two_stage_stage_2_axpy_ops": stage_2_ops, + "two_stage_clear_ops": clear_ops, + "two_stage_total_ops": stage_1_ops + stage_2_ops + clear_ops, + } + ) + + return geometry + + STEP_ACCESSORS = [ ("FitImaging construct", None), # handled specially (constructor) ("Blurred image (FFT convolve)", lambda fit: fit.blurred_image), @@ -268,8 +453,11 @@ def _f_block_scratch(inversion) -> np.ndarray: lambda fit: fit.inversion, ), ("PSF-weighted data [numba]", lambda fit: fit.inversion.psf_weighted_data), + (MGE_PROFILE_IMAGE_LABEL, _mge_mapping_matrix), + (MGE_BLURRING_IMAGE_LABEL, _mge_blurring_mapping_matrix), + (MGE_CONVOLVE_LABEL, _mge_convolved), ( - "MGE operated mapping matrix (60 funcs)", + MGE_RESIDUAL_LABEL, lambda fit: fit.inversion.linear_func_operated_mapping_matrix_dict, ), ( @@ -351,6 +539,16 @@ def one_decomposed_evaluation() -> tuple[dict[str, float], float]: assert isinstance(fit_check.inversion, InversionImagingSparseNumba), ( f"Expected InversionImagingSparseNumba, got {type(fit_check.inversion).__name__}" ) + +# One-off, untimed (PyAutoArray#507 step 0). +mge_split_reproduces_step = _mge_split_reproduces_step(fit_check) +geometry_constants = _geometry_constants(fit_check) + +print("\n--- Geometry constants (model-independent) ---") +print(json.dumps(geometry_constants, indent=2)) +print(f" MGE split reproduces step bit-identically: {mge_split_reproduces_step}") +print(f" MGE blurring stack path: {_mge_scratch.get('blurring_stack_path')}") + del fit_check n_repeats = 10 @@ -375,9 +573,24 @@ def one_decomposed_evaluation() -> tuple[dict[str, float], float]: + _step_dict[F_MAPPER_FUNC_LABEL] + _step_dict[F_FUNC_FUNC_LABEL] ) + +# The same scheme for the MGE row (PyAutoArray#507 step 0): the three timed +# pieces of the operated-mapping-matrix override are subtracted from the raw +# step, so the four MGE rows sum to the raw, un-split MGE cost recorded here. +mge_operated_mapping_matrix_total = _step_dict[MGE_RESIDUAL_LABEL] +_mge_piece_total = ( + _step_dict[MGE_PROFILE_IMAGE_LABEL] + + _step_dict[MGE_BLURRING_IMAGE_LABEL] + + _step_dict[MGE_CONVOLVE_LABEL] +) + +_residual_of = { + F_RESIDUAL_LABEL: curvature_matrix_f_total - _f_block_total, + MGE_RESIDUAL_LABEL: mge_operated_mapping_matrix_total - _mge_piece_total, +} + likelihood_steps = [ - (label, curvature_matrix_f_total - _f_block_total if label == F_RESIDUAL_LABEL else per_call) - for label, per_call in likelihood_steps + (label, _residual_of.get(label, per_call)) for label, per_call in likelihood_steps ] start = time.perf_counter() @@ -433,6 +646,21 @@ def one_decomposed_evaluation() -> tuple[dict[str, float], float]: "steps": {label: per_call for label, per_call in likelihood_steps}, "total_step_by_step": step_total, "curvature_matrix_f_total_s": curvature_matrix_f_total, + "mge_operated_mapping_matrix_total_s": mge_operated_mapping_matrix_total, + "mge_operated_mapping_matrix_split_note": ( + "The four 'MGE: ...' / 'MGE operated mapping matrix' rows sum to " + "mge_operated_mapping_matrix_total_s, the raw un-split cost of the " + "inversion.linear_func_operated_mapping_matrix_dict step. The three " + "piece rows re-run, uncached, the three parts of " + "LightProfileLinearObjFuncList.operated_mapping_matrix_override (60 " + "profile images on the data grid, 60 on the blurring grid, one batched " + "real-space PSF convolution); the residual row is the step minus those " + "three. mge_split_reproduces_step records that the three pieces " + "reconstruct the step's output bit-identically. PyAutoArray#507 step 0." + ), + "mge_split_reproduces_step": mge_split_reproduces_step, + "mge_blurring_stack_path": _mge_scratch.get("blurring_stack_path"), + "geometry": geometry_constants, "curvature_matrix_f_split_note": ( "The four 'F: ...' / 'Curvature matrix F' rows sum to " "curvature_matrix_f_total_s, the raw un-split cost of the " diff --git a/scripts/imaging/likelihood_breakdown/pixelization_numba.py b/scripts/imaging/likelihood_breakdown/pixelization_numba.py index 38bf96684..fe7af79cd 100644 --- a/scripts/imaging/likelihood_breakdown/pixelization_numba.py +++ b/scripts/imaging/likelihood_breakdown/pixelization_numba.py @@ -19,7 +19,14 @@ 4. Inversion build — ray trace, adaptive mesh + interpolation, mapper 5. PSF-weighted data [numba] 6. MGE operated mapping matrix — the 60 linear-Gaussian lens-light images - + their PSF convolution (cached, reused by data vector + curvature) + + their PSF convolution (cached, reused by data vector + curvature), + split into its three pieces plus the assembly residual — see the "MGE + operated mapping matrix sub-block instrumentation" note beside + ``STEP_ACCESSORS``: + 6a. 60x image_2d_from(grid) [numpy profile evaluation] + 6b. 60x image_2d_from(blurring_grid) [numpy profile evaluation] + 6c. 1x batched real-space PSF convolution + 6d. MGE residual: dict assembly 7. Mapper sparse triplets — the sparse [data -> source-pixel] mapping weights (cached, reused by data vector + curvature) 8. Data vector D [numba] @@ -195,6 +202,7 @@ def _profiling_root() -> _Path: from autoarray.inversion.inversion.imaging_numba.sparse import ( # noqa: E402 InversionImagingSparseNumba, ) +from autoarray.inversion.mappers.abstract import Mapper # noqa: E402 n_image_pixels = dataset.data.shape[0] n_over_sampled_pixels = dataset.grids.lp.over_sampled.shape[0] @@ -270,14 +278,201 @@ def _f_block_scratch(inversion) -> np.ndarray: return scratch +# ------------------------------------------------------------------- +# MGE operated mapping matrix sub-block instrumentation (PyAutoArray#507 step 0) +# ------------------------------------------------------------------- +# +# ``inversion.linear_func_operated_mapping_matrix_dict`` builds, for the linear +# MGE lens light, the PSF-convolved image of every one of its 60 Gaussians. On +# the numpy path that is +# ``LightProfileLinearObjFuncList.operated_mapping_matrix_override`` +# (autogalaxy/profiles/light/linear/abstract.py), whose fast branch is three +# distinct pieces of work: +# +# 1. ``mapping_matrix`` -- 60x ``image_2d_from(grid)``, stacked +# 2. the blurring stack -- 60x ``image_2d_from(blurring_grid)``, stacked +# 3. one batched real-space PSF convolution of the two stacks +# +# The override is a ``cached_property`` on the linear-func object and the three +# pieces below are recomputed from scratch (``mapping_matrix`` is a plain +# ``property``), so timing them does *not* prime the step that follows. As with +# the F sub-blocks, the pre-existing row is therefore reported as a RESIDUAL, +# +# MGE residual = t(linear_func_operated_mapping_matrix_dict) - (t1 + t2 + t3) +# +# so the four MGE rows sum to exactly the un-instrumented step and the artifact +# stays comparable with pre-instrumentation runs. The raw (unsplit) MGE total is +# written to the result JSON as ``mge_operated_mapping_matrix_total_s``. +# +# ``mge_split_reproduces_step`` in the JSON records a one-off bit-identical +# check that the three timed pieces really do reconstruct the dict the step +# returns -- if a dataset ever takes a different branch of the override the +# split is flagged rather than silently mis-attributed. +# +# Caveat (as for F): the residual is a difference of averaged timings and can in +# principle come out slightly negative. It is recorded as measured, never clipped. + +MGE_PROFILE_IMAGE_LABEL = "MGE: image_2d_from(grid) x60 [numpy]" +MGE_BLURRING_IMAGE_LABEL = "MGE: image_2d_from(blurring_grid) x60 [numpy]" +MGE_CONVOLVE_LABEL = "MGE: batched PSF convolution [real-space np]" +MGE_RESIDUAL_LABEL = "MGE operated mapping matrix [residual: dict assembly]" + +_mge_scratch: dict[str, list] = {} + + +def _mge_linear_func_list(inversion) -> list: + """The inversion's linear-func objects (the MGE lens light), never its mappers.""" + return [obj for obj in inversion.linear_obj_list if hasattr(obj, "light_profile_list")] + + +def _mge_mapping_matrix(fit) -> list: + """Piece 1: the 60 unblurred profile images on the (over-sampled) data grid.""" + matrices = [linear_func.mapping_matrix for linear_func in _mge_linear_func_list(fit.inversion)] + _mge_scratch["mapping_matrix"] = matrices + return matrices + + +def _mge_blurring_stack(linear_func): + """The blurring stack of one linear-func object, built the way the override builds it. + + ``operated_mapping_matrix_override`` stacked the 60 blurring-grid images with an + inline per-profile loop until PyAutoArray#507 step 3, which moved both stacks behind + ``LightProfileLinearObjFuncList._image_slim_list_from`` so that an MGE basis shares one + reference-frame transform and one eccentric-radius grid. This harness is run against + both libraries (the paired before/after of that step), so it dispatches on the + attribute rather than hard-coding either shape -- a hand-rolled loop would keep timing + the *old* code on the new library and mis-attribute the win to the residual row. Which + branch was taken is recorded in the result JSON as ``mge_blurring_stack_path``. + """ + if hasattr(linear_func, "_image_slim_list_from"): + _mge_scratch["blurring_stack_path"] = "shared_geometry" + + return np.stack( + linear_func._image_slim_list_from(grid=linear_func.blurring_grid, xp=np), + axis=1, + ) + + _mge_scratch["blurring_stack_path"] = "per_profile_loop" + + return np.stack( + [ + light_profile.image_2d_from(grid=linear_func.blurring_grid, xp=np).slim.array + for light_profile in linear_func.light_profile_list + ], + axis=1, + ) + + +def _mge_blurring_mapping_matrix(fit) -> list: + """Piece 2: the same 60 profiles on the blurring grid (flux blurred in from + outside the mask), stacked exactly as the override stacks them.""" + matrices = [ + _mge_blurring_stack(linear_func) for linear_func in _mge_linear_func_list(fit.inversion) + ] + _mge_scratch["blurring_mapping_matrix"] = matrices + return matrices + + +def _mge_convolved(fit) -> list: + """Piece 3: the single batched real-space convolution of the two stacks.""" + return [ + linear_func.psf.convolved_mapping_matrix_via_real_space_np_from( + mapping_matrix=_mge_scratch["mapping_matrix"][index], + mask=linear_func.grid.mask, + blurring_mapping_matrix=_mge_scratch["blurring_mapping_matrix"][index], + blurring_mask=linear_func.blurring_grid.mask, + ) + for index, linear_func in enumerate(_mge_linear_func_list(fit.inversion)) + ] + + +def _mge_split_reproduces_step(fit) -> bool: + """Bit-identical check that pieces 1-3 reconstruct the timed step's output.""" + _mge_mapping_matrix(fit) + _mge_blurring_mapping_matrix(fit) + reconstructed = _mge_convolved(fit) + + actual = list(fit.inversion.linear_func_operated_mapping_matrix_dict.values()) + + if len(actual) != len(reconstructed) or len(actual) == 0: + return False + + return all(np.array_equal(np.asarray(a), np.asarray(b)) for a, b in zip(actual, reconstructed)) + + +def _geometry_constants(fit) -> dict: + """Model-independent geometry of the sparse operator and of each mapper. + + Recorded once, outside the timed loop, so that the complexity of the + mapper x mapper kernel (`curvature_matrix_via_sparse_operator_from`) is a + measurement rather than an estimate. `mapper_mapper_inner_ops` is the exact + count of that kernel's innermost accumulations, + + sum over stored (data_0, data_1) pairs of u0(data_0) * u1(data_1), + + and the `two_stage_*` counts are the same geometry costed for the phase-2 + reformulation (a dense n_source accumulator per data pixel, then one AXPY + per mapping). PyAutoArray#507 step 0. + """ + inversion = fit.inversion + + lengths = np.asarray(inversion.sparse_operator.lengths).astype("int64") + indexes = np.asarray(inversion.sparse_operator.indexes).astype("int64") + + geometry = { + "psf_shape_native": [int(s) for s in dataset.psf.kernel.shape_native], + "image_pixels_masked": int(n_image_pixels), + "psf_precision_pairs_stored": int(indexes.shape[0]), + "psf_precision_lengths_sum": int(lengths.sum()), + "psf_precision_lengths_mean": float(lengths.mean()), + "psf_precision_lengths_max": int(lengths.max()), + "mappers": [], + } + + for mapper in inversion.cls_list_from(cls=Mapper): + pix_lengths = np.asarray(mapper.unique_mappings.pix_lengths).astype("int64") + + pair_u0 = np.repeat(pix_lengths, lengths) + pair_u1 = pix_lengths[indexes] + + n_source = int(mapper.params) + data_pixels = int(lengths.shape[0]) + + stage_1_ops = int(pair_u1.sum()) + stage_2_ops = int(pix_lengths.sum()) * n_source + clear_ops = data_pixels * n_source + + geometry["mappers"].append( + { + "params": n_source, + "data_to_pix_unique_shape": [ + int(s) for s in np.asarray(mapper.unique_mappings.data_to_pix_unique).shape + ], + "pix_lengths_mean": float(pix_lengths.mean()), + "pix_lengths_max": int(pix_lengths.max()), + "pix_lengths_sum": int(pix_lengths.sum()), + "mapper_mapper_inner_ops": int((pair_u0 * pair_u1).sum()), + "two_stage_stage_1_scatter_ops": stage_1_ops, + "two_stage_stage_2_axpy_ops": stage_2_ops, + "two_stage_clear_ops": clear_ops, + "two_stage_total_ops": stage_1_ops + stage_2_ops + clear_ops, + } + ) + + return geometry + + STEP_ACCESSORS = [ ("FitImaging construct", None), # handled specially (constructor) ("Blurred image (FFT convolve)", lambda fit: fit.blurred_image), ("Profile subtracted image", lambda fit: fit.profile_subtracted_image), ("Inversion build (trace+mesh+mapper)", lambda fit: fit.inversion), ("PSF-weighted data [numba]", lambda fit: fit.inversion.psf_weighted_data), + (MGE_PROFILE_IMAGE_LABEL, _mge_mapping_matrix), + (MGE_BLURRING_IMAGE_LABEL, _mge_blurring_mapping_matrix), + (MGE_CONVOLVE_LABEL, _mge_convolved), ( - "MGE operated mapping matrix (60 funcs)", + MGE_RESIDUAL_LABEL, lambda fit: fit.inversion.linear_func_operated_mapping_matrix_dict, ), ( @@ -360,6 +555,16 @@ def one_decomposed_evaluation() -> tuple[dict[str, float], float]: assert isinstance(fit_check.inversion, InversionImagingSparseNumba), ( f"Expected InversionImagingSparseNumba, got {type(fit_check.inversion).__name__}" ) + +# One-off, untimed (PyAutoArray#507 step 0). +mge_split_reproduces_step = _mge_split_reproduces_step(fit_check) +geometry_constants = _geometry_constants(fit_check) + +print("\n--- Geometry constants (model-independent) ---") +print(json.dumps(geometry_constants, indent=2)) +print(f" MGE split reproduces step bit-identically: {mge_split_reproduces_step}") +print(f" MGE blurring stack path: {_mge_scratch.get('blurring_stack_path')}") + del fit_check n_repeats = 10 @@ -384,9 +589,24 @@ def one_decomposed_evaluation() -> tuple[dict[str, float], float]: + _step_dict[F_MAPPER_FUNC_LABEL] + _step_dict[F_FUNC_FUNC_LABEL] ) + +# The same scheme for the MGE row (PyAutoArray#507 step 0): the three timed +# pieces of the operated-mapping-matrix override are subtracted from the raw +# step, so the four MGE rows sum to the raw, un-split MGE cost recorded here. +mge_operated_mapping_matrix_total = _step_dict[MGE_RESIDUAL_LABEL] +_mge_piece_total = ( + _step_dict[MGE_PROFILE_IMAGE_LABEL] + + _step_dict[MGE_BLURRING_IMAGE_LABEL] + + _step_dict[MGE_CONVOLVE_LABEL] +) + +_residual_of = { + F_RESIDUAL_LABEL: curvature_matrix_f_total - _f_block_total, + MGE_RESIDUAL_LABEL: mge_operated_mapping_matrix_total - _mge_piece_total, +} + likelihood_steps = [ - (label, curvature_matrix_f_total - _f_block_total if label == F_RESIDUAL_LABEL else per_call) - for label, per_call in likelihood_steps + (label, _residual_of.get(label, per_call)) for label, per_call in likelihood_steps ] # Cross-check: a directly timed, undecomposed evaluation via the production @@ -443,6 +663,21 @@ def one_decomposed_evaluation() -> tuple[dict[str, float], float]: "steps": {label: per_call for label, per_call in likelihood_steps}, "total_step_by_step": step_total, "curvature_matrix_f_total_s": curvature_matrix_f_total, + "mge_operated_mapping_matrix_total_s": mge_operated_mapping_matrix_total, + "mge_operated_mapping_matrix_split_note": ( + "The four 'MGE: ...' / 'MGE operated mapping matrix' rows sum to " + "mge_operated_mapping_matrix_total_s, the raw un-split cost of the " + "inversion.linear_func_operated_mapping_matrix_dict step. The three " + "piece rows re-run, uncached, the three parts of " + "LightProfileLinearObjFuncList.operated_mapping_matrix_override (60 " + "profile images on the data grid, 60 on the blurring grid, one batched " + "real-space PSF convolution); the residual row is the step minus those " + "three. mge_split_reproduces_step records that the three pieces " + "reconstruct the step's output bit-identically. PyAutoArray#507 step 0." + ), + "mge_split_reproduces_step": mge_split_reproduces_step, + "mge_blurring_stack_path": _mge_scratch.get("blurring_stack_path"), + "geometry": geometry_constants, "curvature_matrix_f_split_note": ( "The four 'F: ...' / 'Curvature matrix F' rows sum to " "curvature_matrix_f_total_s, the raw un-split cost of the " diff --git a/scripts/misc/likelihood_breakdown/README.md b/scripts/misc/likelihood_breakdown/README.md index c37b0e753..f56648c68 100644 --- a/scripts/misc/likelihood_breakdown/README.md +++ b/scripts/misc/likelihood_breakdown/README.md @@ -26,7 +26,7 @@ For *how long the likelihood actually takes* on production hardware — i.e. a s | `imaging/delaunay` | hst | hpc_a100_mp | dense (mapping) | 96.8 ms | v2026.7.6.649 | | `imaging/delaunay` | hst | hpc_a100_mp | sparse (w-tilde) | 95.5 ms | v2026.7.6.649 | | `imaging/delaunay_numba` | euclid | local_cpu_fp64 | sparse (numba) | 1.19 s | v2026.8.17.1 | -| `imaging/delaunay_numba` | hst | local_cpu_fp64 | sparse (numba) | 679.2 ms | v2026.8.17.1 | +| `imaging/delaunay_numba` | hst | local_cpu_fp64 | sparse (numba) | 482.4 ms | v2026.8.17.1 | | `imaging/mge` | hst | local_cpu_fp64 | dense (mapping) | 179.5 ms | v2026.7.6.649 | | `imaging/mge` | hst | hpc_a100_fp64 | dense (mapping) | 7.8 ms | v2026.7.6.649 | | `imaging/pixelization` | hst | local_cpu_fp64 | dense (mapping) | 8.65 s | v2026.7.6.649 | @@ -35,8 +35,8 @@ For *how long the likelihood actually takes* on production hardware — i.e. a s | `imaging/pixelization` | hst | hpc_a100_fp64 | sparse (w-tilde) | 57.9 ms | v2026.7.6.649 | | `imaging/pixelization` | hst | hpc_a100_mp | dense (mapping) | 56.4 ms | v2026.7.6.649 | | `imaging/pixelization` | hst | hpc_a100_mp | sparse (w-tilde) | 55.5 ms | v2026.7.6.649 | -| `imaging/pixelization_numba` | euclid | local_cpu_fp64 | sparse (numba) | 288.6 ms | v2026.8.17.1 | -| `imaging/pixelization_numba` | hst | local_cpu_fp64 | sparse (numba) | 691.5 ms | v2026.8.17.1 | +| `imaging/pixelization_numba` | euclid | local_cpu_fp64 | sparse (numba) | 124.3 ms | v2026.8.17.1 | +| `imaging/pixelization_numba` | hst | local_cpu_fp64 | sparse (numba) | 304.5 ms | v2026.8.17.1 | Auto-generated by `scripts/build_readme.py` from the versioned artifacts under `results/breakdown/`.