Restrict "tweaked moments" for "free shape"

[HannoSpreeuw]

1. Consider this the equivalent of commit f3b10b1. That commit tightened the condition for applying "tweaked moments" for forced beams: apply "tweaked moments" only when the rounded barycenter position is part of the island. The current commit tightens the conditions for "tweaked moments" for a "free shape" in the same way in the sense that reverting to a maximum pixel value and corresponding position is no longer possible; i.e. "basevalue" and "basepos" can no longer be used as the base point for a "tweaked moments" extrapolation - unless it coincides with the rounded barycenter, of course. We noticed that in that case the extrapolations from "tweaked moments" could reach unrealistic values. They are bounded, that is why they did not show up as prominently as for forced beams.

2. Fixed a serious bug that occurs for "ill-formed" or large islands: "posx.size==posy.size!=no_pixels" is the source of error. Since we have indices "[0][0]" in the line "i = np.nonzero(mask)[0][0]" this will often still end up right, but we have encountered cases from a collection of 60 GRB images described here: Unrealistic values with --force-beam and --vectorize #180 (comment) where a run "pyse GRB201006A/*.fits --vectorized --back-size-x 50 --back-size-y 50" will reveal serious problems, i.e. unrealistic peak brightnesses. This can occur on one machine, while not causing any problems on another machine with exactly the same code. This comes from the way the "np.empty" fills values in these lines in the "extract" module:

   pyse/sourcefinder/extract.py

   Line 2229 in 591717f

   xpositions = np.empty((num_islands, max_pixels), dtype=np.int32)

   
   "xpositions = np.empty((num_islands, max_pixels), dtype=np.int32)"
   and

   pyse/sourcefinder/extract.py

   Line 2230 in 591717f

   ypositions = np.empty((num_islands, max_pixels), dtype=np.int32)

   
   "ypositions = np.empty((num_islands, max_pixels), dtype=np.int32)"
   What is called "xpositions" and "ypositions" in the "extract" module is propagated to "posx" and "posy" in "fitting.moments_enhanced", respectively. The problem can occur for any island smaller than the largest island, i.e. for "no_pixels < max_pixels" from this line

   pyse/sourcefinder/measure.py

   Line 498 in 591717f

   mask = (posx == rounded_barycenter[0]) & (posy == rounded_barycenter[1])

   "
   mask = (posx == rounded_barycenter[0]) & (posy == rounded_barycenter[1])
   "
   which generally causes "mask.size>no_pixels", i.e. "mask" generally becomes (much) larger than the island. This should never occur.

3. After the fix from 1) we end up with a maximum ratio peak/maxi of 1.62 from the 60 images described above, see this comment: Add regression test to cross-check against non-vectorized output. #182 (comment). That is why "upp_bound = 2.0 * basevalue" seems appropriate, in the sense that this bound will never be reached. This implies that we could replace that line by "upp_bound = np.inf", but that seems a bit risky at this stage; we have not explored enough pathological cases yet.

4. The comment about the meaning of the maxi input to measure.moments_enhanced should be clearer now.

[HannoSpreeuw]

Note: this does not yet fix #182 , since this does not include a regression test to detect the errors described under 1) and 2). The current images under test/data do not seem to be suitable for this and adding a new image could take too much GH bandwidth for CI.

[HannoSpreeuw]

• The equivalent code change for measure.moments to tackle 1) still needs to be implemented; 2) does not seem to be a problem for measure.moments.

• In terms of equivalence between measure.moments and measure.moments_enhanced one may seek an equivalent of forced beam for measure.moments.