image_dynamic_range cannot detect over-subtraction: wrong image, wrong location, wrong direction

[Athanaseus]

Summary

image_dynamic_range already returns exactly the statistics needed to spot an
over-subtracted source:

{"deepest_negative": peak_flux / abs(min_flux),
 "local_rms":        peak_flux / local_std,
 "global_rms":       peak_flux / global_std}

But of a DI/DD image pair where the DD run over-subtracted a peeled source into a
−24.7 mJy hole, it reports the DD image as better:

	DI	DD	verdict
deepest_negative (current behaviour)	83.5	93.4	DD looks better

The same formula applied to the restored image, field-wide, gives:

	DI	DD	verdict
peak / abs(min) on restored, field-wide	137.9	1.7	DD 81× worse

So the metric is right; what it is measured on is not.

Why it misses it — four independent reasons

Current implementation (aimfast.py ~874-892):

min_flux   = residual_data[imslice].min()    # RESIDUAL, box at the restored peak
local_std  = target_area.std()               # RESIDUAL, same box
global_std = residual_data[0,0,...].std()    # RESIDUAL, whole image
peak_flux  = restored_data.max()             # RESTORED

1. Wrong image. The hole lives in the restored image (−24.74 mJy). The DD residual
minimum inside that box is only −0.45 mJy.

2. Wrong location. The box is centred on each image's own brightest pixel — and peeling
moves the peak:

DI peak: pixel (4564,2986)  RA 355.2791 DEC +0.3096   162.24 mJy
DD peak: pixel (4613,2445)  RA 355.2587 DEC +0.0841    42.21 mJy
                                             0.226 deg apart

Each image nominates a different reference source, with peak values differing by 3.8×. That
is not a comparison: it measures noise near source A in one image and near source B in the
other. The damage at RA 355.2791 is never inside DD's box.

3. Wrong direction. Over-subtraction flattens residuals. Residual-based metrics
therefore reward the failure mode.

4. The "box" is not a box — it is four whole rows

The intended region is area_factor beams around the peak (32 px ≈ 48″ with BMAJ 7.65″,
area_factor=6, 1.5″ pixels). But imslice is built as four edge coordinates and then
used as arr[imslice], which numpy interprets as fancy indexing selecting four rows:

imslice = np.array([pix_coord[2] - ra_num_pix/2, pix_coord[2] + ra_num_pix/2,
                    pix_coord[3] - dec_num_pix/2, pix_coord[3] + dec_num_pix/2])
imslice = np.array(list(map(int, imslice)))
...
target_area = residual_data[0, frq_ax, :, :][imslice]

On a 6075² image:

imslice = [2970 3002 4548 4580]
target_area.shape = (4, 6075)      # expected ~(32, 32)

So local_rms and deepest_negative are computed over four full-width rows spanning the
entire image — 24,300 pixels — rather than a box at the source. Worse, np.argwhere returns
[stokes, chan, y, x], so pix_coord[3] is the x coordinate being used as a row index:
two of the four rows (4548, 4580) sit at arbitrary declinations unrelated to the source.

Consequences: local_rms is not local, and deepest_negative finds the minimum anywhere
along those rows rather than near the reference source.

Fix: slice a real 2D box, clipped to the image bounds, with dec along y and ra along x
(and abs() on the pixel counts, since CDELT1 is normally negative).

This changes reported DR values for all existing users. local_rms and
deepest_negative will differ from previous runs — previously they described four arbitrary
image rows, now they describe the intended neighbourhood of the reference position.

Which factor dominates

Same metric, three variants, same data:

what is measured	DI	DD	verdict
residual, "box" at own peak (current)	83.5	93.4	DD better
restored, real box at own peak	322.5	44.2	DD 7.3× worse
restored, real box at pinned position	322.5	1.00	DD 323× worse
restored, field-wide min	137.9	1.7	DD 81× worse

With the fixes applied, pinning DD to DI's reference position reports
ref_peak_flux = 24.74 mJy at RA 355.2791 DEC +0.3096 — i.e. it finds the over-subtracted
hole, which is exactly the intended behaviour.

Switching image helps; pinning the reference position helps far more; field-wide most.

Proposed fix (four independent parts, each useful alone)

0. Allow the reference position to be pinned, and always report it. Cheapest change, and
the only one that makes DR an actual comparison rather than two unrelated measurements.

image_dynamic_range sees exactly one image, and results are stored per image
(output_dict[restored_label]), so a comparison happens by reading the results afterwards.
That per-image design is deliberate and should stay. The minimal fix is therefore to make the
per-image result self-describing:

def image_dynamic_range(fitsname, residual, area_factor=6, ref_position=None)

• ref_position=None → current behaviour exactly (locate the peak in this image), so existing
  callers and self-cal DR tracking are unaffected
• ref_position=(ra, dec) in degrees → measure there instead
• always return ref_position and ref_peak_flux in the result dict, so the position lands
  in fidelity_results.json
• raise if a supplied position falls outside the image or lands on blanked pixels —
  deliberately not falling back to this image's own peak, since silently measuring somewhere
  other than the position asked for is precisely the failure this argument exists to prevent

The silent failure then becomes visible: two DR entries carry the positions they were measured
at, and 0.226° apart on different sources is obvious on inspection. A warning when comparing
entries measured at different positions closes it completely.

"Use image A's peak for image B" needs no extra parameter — it is composition:

dr_a = image_dynamic_range(img_a, res_a)
dr_b = image_dynamic_range(img_b, res_b, ref_position=dr_a["ref_position"])

Deliberately not accepting an image as the reference argument: a parameter taking either a
tuple or a path is awkward to validate and document, and a second ref_image parameter would
need mutual-exclusion logic for a case already served by the above.

This also repairs the existing use case. Tracking DR across self-cal rounds where peeling
occurs: when the brightest source is peeled, the metric silently switches to a different
source, so the DR curve jumps because the reference object changed, not the image quality —
and nothing reports it. At minimum, report which source/position the DR refers to so a switch
becomes visible.

1. Let the user select restored vs residual. Small, backwards-compatible.

2. Stop measuring only one box at the peak. Field-wide minimum, or per-source boxes —
which is where this meets the sign-blindness of _source_residual_results (separate issue).

Not a bug in the original design

DR was built to answer "how clean is it around the brightest source" — the classic
dynamic-range question, and it does that correctly. Finding damage elsewhere in the field was
never in scope. This is an extension, with the caveat that peeling silently invalidates
the implicit assumption that the brightest source is the same object in both images.

Workaround available today

image_dynamic_range(restored, restored) — no code change; gets the restored-image variant,
though still only the box at the peak.

Tests

• Two synthetic images, one with a known negative hole away from the peak: assert the
  restored/field-wide variant flags it and the current residual/peak-box variant does not.
• Assert that when comparing two images, both are measured at the same sky position.
• Assert the reported output names the reference source/position used.