Convolution consistency - #5782

This test is wrong: https://github.com/astropy/astropy/blob/master/astropy/convolution/tests/test_convolve.py#L120
tests should be failing now but are not. I don't get it yet.

(it's wrong because this is a 1D test and the output is 3D - that can not happen.)

EDIT: it was wrong because it was broadcasting the comparison operation!

@astrofrog - I think the behavior of convolve simply doesn't make sense in some cases. We need to chat offline about this. For example, this convolution:

x = np.array([[0., 0., 3.],
              [1., np.nan, 0.],
              [0., 2., 0.]], dtype='float64')

y = np.array([[1., 1., 1.],
              [1., 1., 1.],
              [1., 1., 1.]], dtype='float64')

convolve gives:

In [3]: convolve(x,y)
Out[3]:
array([[ 1.75,  4.75,  3.75],
       [ 3.75,  6.75,  5.75],
       [ 3.75,  3.75,  2.75]])

which... I can't reconcile with any interpretation of convolution.

Just as problematic, that is not the value being tested for in the tests, but the test is not failing.

(please tell me if you get a different local result)

I believe, after re-reading #926, that the convolve_fft behavior is correct and convolve behavior is incorrect in the presence of NaNs. But, the keyword interpolate_nan may be misleading. The intent of convolve_fft's interpolate_nan functionality is to implement the "normalized convolution" method @danielwe described in #926. By contrast, convolve is performing a kernel-weighted average at the position of the NaN, which is fine for that value, but results in unexpected and weird behavior for other surrounding pixels.

In the example above, the middle NaN is assigned a temporary value of 0.75 = 6/8 = (1+3+2) / 8, where 8 is the number of non-NaN values being averaged together. Unfortunately, the kernel is not an averaging (normalized) kernel, but an additive non-normalized kernel. That results in this 0.75 being added to every neighbor. I don't think this is ever a useful operation.

So, my plan is to replace this weighted average with a simple zeroing of the NaN value... which I think will have the desired effect, but it will take some guess and check or a deeper investigation of the code to get it right.

I recently wrote a bunch of (numba-based) convolution functions because I disliked the interpolation before convolution (I remarked on that on astropy-dev a long time ago).

However there are only three options when you have invalid values (nan or masked):

• interpolate them before convolution
• ignore them during convolution and re-normalize the header
• replace them by a fill value

I'd rather think that ignoring and re-normalizing the header is the best approach but that depends on the size of the kernel (accuracy) and the applied convolution method (sum/average/median). I don't know what would be the "least astonishment" option here but there will always be differences in opinion. But changing how these are handled is definetly breaking backwards-compatibility so should not be taken lightly.

Or am I totally wrong and there is a consensus on the "right way"?

Note: I haven't looked at the actual changes, that's just picking up your last comment.

@MSeifert04 I don't think there is a 'consensus right way' - both "ignore and renormalize" and "replace with fill" are valid options. I'm not convinced there's ever a case where 'interpolate-then-convolve' is a reasonable behavior for a convolution function.

I also agree that this will break backward compatibility and is not a small change. It is probably necessary to preserve all three behaviors in some way, though I would like to deprecate the interpolate-first option. I don't want to provide the interpolate-first option as part of convolve_fft either - if someone really wants that behavior, it's straightforward to produce anyway by doing one extra convolution.

Is it possible for us to include numba functions in astropy? Or should we stick with cython?

And looking at that astropy-dev conversation.... well, you and I were right, and the current convolve behavior does not make much sense. I don't even remember having that conversation.

I think part of the issue is that prior to that discussion, we had been using the term 'interpolate' to mean two different things - in convolve, it meant interpolate-then-convolve, while in convolve_fft, it meant interpolate-by-normalized-convolution.

I think numba isn't a good dependency. It's a mess to install LLVM and numba without conda.

In case anyone is following along, I'm completely changing the behavior of convolve to match that of convolve_fft. This will absolutely be a backwards-incompatible change because of the added keywords and a different default behavior, but I think the change is critical. The most important component of this, which I haven't even started on, will be improving the documentation to illustrate what we mean by the various options, because anyone who has not coded this up themselves will have a hard time understanding the behaviors just by skimming the docstring.

One question has come up that warrants some feedback: is interpolate_nan too confusing a keyword name? Should we rename interpolate_nan to normalized_convolution or something similar? I had never heard of the term prior to #926, which is why the current keyword is there, but perhaps using an 'official' term might be better?

@keflavich - I need to catch up with all this, but was wondering whether it wouldn't make more sense to change the default convolve_fft behavior - specifically I find the 'classic' result here to be the most sensible default: #926 (comment) - what do you think?

@astrofrog - Yes, probably the default behavior of convolve_fft should also change if that's even possible (I don't remember right now, but it's possible that convolve_fft cannot interpolate without a normalized kernel).

However, I now firmly believe that the default behavior of convolve in the presence of NaNs is incorrect, in the sense that is not behavior users should ever expect. If you'd like to talk about this more, let me know, but I think the reasons show up in this and other PRs.

@keflavich - I'd like to discuss this more - I'm traveling for the next couple of weeks but maybe we can plan to talk during the week of March 13th? (if so I can send you a doodle)

Side note: some users clearly do expect and like the NaN behavior in convolve: astropy/astropy.github.com#117 (comment) :)

Just to put this out there, one thing we could do to not break backward compatibility in any bad way is to actually rename the functions and keep the old ones as wrappers of the new ones with the correct defaults. Since the main benefit of these convolution functions is to deal with NaNs in some way (compared to the convolve functions in scipy) we could always make the functions e.g. nan_friendly_convolve or something like that and point users to the scipy convolve in cases where NaNs are not an issue.

I'm traveling a lot over the next few weeks as well, but let's find a time to talk.

We could split out the nan-fixing interpolation, but I think we need to always provide a way to access it via keywords from the default convolve functions.

And, regarding users' liking the code: astropy's convolve (both variants) is the only version I know of that can handle missing data at all, so it's no surprise users like it, but we should be delivering correct results!

I've completed the code fixes. There are some new open questions, and some more tests need to be written.

Tests - checking convolve - convolve_fft consistency.

Question: Do we change the default parameters to interpolate_nan=True, normalize_kernel=True? If so, how do we do the transition?

circleci tests behave different from my own machine; they seem to reset the context between script calls within convolutions/index.rst.

EDIT: travis does this too. I'm very sad now, this was a good solution that worked nicely, now I have to ditch everything.

@astrofrog how did you make the docs in wcsaxes not try to access remote data on appveyor? My additions are causing this error: https://ci.appveyor.com/project/Astropy/astropy/build/1.0.7154/job/c8bwnq9v73aufi2p

problem line in question says remote_data: none instead of remote_data: astropy:
https://ci.appveyor.com/project/Astropy/astropy/build/1.0.7154/job/c8bwnq9v73aufi2p#L2496

There are no failed tests on travis-ci, only timeout errors. These seem to be unrelated to this PR.

Appveyor is failing because the documentation tests (i.e., the tests of the documentation, not the embedded docs test) are failing as noted above

This PR is ready for final review if someone can help figure out the doc failures on travis and the illegal remote access on appveyor.

For the Appveyor error, it looks like your doc example tries to grab remote data when it is not marked to do so.

@pllim - Either both @keflavich and I were both very tired last week or something is fundamentally not right. The docs seems to use the same logic and thus remote data as the wcsdocs wcsaxes, yet the latter works while this one isn't.
Also I feel that the timeouts are related to this PR, we don't see them anywhere else.

Are you sure the WCS docs actually run through doctest? A lot of FITS and WCS stuff get a "docskip all" treatment, right?

I meant wcsaxes...

Looks like this is now ready to go as soon as the tests pass. It's probably good already but so many rebases... ®️⚾️

All green, so merging. Thanks a lot @keflavich for finishing off this marathon PR.

Thanks @keflavich!