Fix to ePSF algorithm changes in v0.7

[Onoddil]

This PR fixes the bug in EPSFBuilder in its v0.7 update. The bug, as best I pinned it down, is due to a subtle interaction between a "halfway" algorithm update to the original Anderson & King (2000) paper, from its form in v0.6 which was, as I understand it, implemented via discussions with Jay, unbeknownst to me at the time. I have subsequently found what I believe to be a citable version of the updates (and referenced them accordingly), and thus re-aligned the ePSF algorithm with its updated version. Ironically if I had made more code changes instead of fewer, and further aligned the code with the workflow of AK00 it would likely have been ok, but there we are...!

The main differences should therefore lie in the removal of the residual smooth and recenter loop of AK00, to follow ISR WFC3 2016-12; for continuity I simply reversed the recentering algorithm completely, keeping the original iteratively-find-the-center-via-centroiding function. There should hopefully be very few differences in the v0.6 and v0.7 results -- some preliminary testing on very poor datasets (16 randomly placed Gaussian stars with sigma of 0.5 pixels, noisy images and poor initial centroids for an oversampling factor of 4) gives results that are not complete gibberish as with the previous examples in issues raised, so hopefully this should successfully revert the bug.

This PR also implements the minor change to include a custom SigmaClip, allowing for custom sigma clip functions to be passed to EPSFBuilder.

This PR does not include changes to one major difference between the v0.6 algorithm and both AK00 and ISR 2016-12: whether the residual sampling to grid point mapping is a nearest neighbour, "box of half a grid spacing length" function or, as is the case in the two cited works, residual samplings are placed in all grid points within a grid spacing (i.e., a box of a full grid spacing length, 0.25 pixels of oversampling=4) for the purposes of sigma-clipping residuals for ePSF grid interpolation point updating. Given the current issue being one of non-obvious changes to cited works being made through offline, internal discussions I felt perhaps it was wise to hold off on this change before confirming whether this change was also a case of out-of-band 'unofficial' updates/suggestions. Is this the case, or should I make a separate PR to implement the swap from a one-to-one sampling:grid point mapping to a each-sampling-is-put-in-four-grid-points mapping? I can see arguments for reducing the distance a sampling point can be from any grid point, to reduce the effects of data points from 'extreme' distances affecting the sigma-clipped residual point.

Finally, this PR has come about through a fairly roundabout route, with several rebases a long the way, so it would be good to get a few pairs of eyes on the code to ensure no further bugs are introduced trying to fix the first bug. I believe the changes to at least be robust to some level of my poking the code with a stick, but cc @eteq @larrybradley for some fresh perspective!

Fixes #969

[pep8speaks]

Hello @Onoddil! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file photutils/psf/epsf.py:

Line 315:80: E501 line too long (116 > 79 characters)

Comment last updated at 2019-11-14 22:06:12 UTC

[Onoddil]

I also haven't included the changelog yet, holding off so I can get some opinions on what type of changes various things are, whether they go in 0.7.1 or 0.8, etc. I assume the bug fix of the EPSFBuilder algorithm is 0.7.1, and therefore the API change of SigmaClip goes along with it, but the bug fix also contains API changes (with the re-introduction of center_accuracy and recentering_boxsize, e.g.)

Finally, this reversal of changes throws into question the now-default centroid_epsf as a) the default centroiding algorithm, and b) useful at all, since it's just for EPSFBuilder currently. The re-inclusion of box_size in _recenter_epsf makes it a bit awkward to handle the AK00 centroiding function, as it just involves the few pixels either side of a given shift value. Opinions on whether that has been handled correctly in the code, and if we should change the default again, are appreciated.

[Onoddil]

Attached is a screengrab of running the notebook in #969 with these changes, showing qualitatively the same PSF as shown in the issue for v0.6's algorithm.

[eteq]

A few minor things to suggest but mostly good now!

Also re: API, since this is mostly just adding back API options that were in 0.6, I think it's fine to do this in 0.7.2 . sigclip is new but it's an easy/simple thing so best just to not delay this further and keep it here (and therefore in 0.7.2).

So you can add it to the changelog in the 0.7.2 section.

[eteq]

Suggested change

print(epsf.data, truth_epsf, np.amax(epsf.data), np.amax(truth_epsf))

(print not needed in the test)

[eteq]

Suggested change

	epsf_builder = EPSFBuilder(oversampling=oversampling, maxiters=20,
	epsf_builder = EPSFBuilder(oversampling=oversampling, maxiters=5,

I checked and if it's 20 it takes >~60 sec in the current master, but if we change it to 5 is still works with this fix but finishes (and correctly fails) in a few seconds on master. So we should change this so it doesn't accidentally hang the tests

[eteq]

I think this name needs to be more specific - maybe flux_residual_sigclip. And the docstring should say a bit more about how it enters into the algorithm as opposed to mechanically where it goes in the code. E.g. "A ... to determine which pixels are ignored based on the flux residual."

[eteq]

There are several issues this should fix - @Onoddil pointed out #969 is fixed by this, and I just checked @larrybradley's notebook from #951 (comment) , with proof that this fix addresses that here:

Note I believe this should fix the other people who mentioned this in #951, as well as #860, but cannot test those. I think we should plan on merging this since it fixes at last part of that, and then ask the folks there to test.

[Onoddil]

I've made the changes as requested by @eteq -- changing sigclip's name, and minor edits to the new test. I had to add the 0.7.2 release to the changelog, so please let me know if I messed up the syntax for the changelog at all and I'll fix that.

The one outstanding issue I have, which might require a quick update after feedback, is whether we should also roll back the default centroiding algorithm. Currently the change in v0.7 holds, where the default centroiding algorithm is centroid_epsf, the algorithm detailed in AK00; however, if this is a larger "roll back miscommunication re: changes to the algorithm since 2000" then perhaps the default centroid algorithm shouldn't be the one laid out in the 2000 paper?

The centre-of-mass algorithm is certainly more flexible as a default, and I do use it in the new testing (although offline tests show that the symmetry centroiding of centroid_epsf also works), so should I put the default back to centroid_com?

[eteq]

After some out-of-band discussion with @larrybradley, the conclusion was that it is indeed better to add the flux_residual_sigclip in 0.8. But so as to make this easy to merge quickly, I've just pushed up a few commits "hiding" the functionality to be merged in this PR for 0.7.2, and then #984 brings it back (which will then be milestoned 0.8).

So hopefully if the tests are happy this is now good to go (but should probably get a glance from @larrybradley before merging)

[larrybradley]

Many thanks, @Onoddil.

[emirkmo]

I don't know if this is the right place to write this but 0.72 still gives me a diverging PSF fit instead of converging to a good fit when using EPSFBuilder. I.e., in 0.71 and 0.72 I get the same behavior as mentioned in: #969

All versions >= 0.7 give the same diverging result. However, 0.6X versions correctly converge. So for now, I had to downgrade to that. This also means Astropy needs to be downgraded to a version before 4.0.

Let me know if an example would help. In my tests, the issue was not solved.