Repository navigation
Make convolve and convolve_fft more consistent #926
Description
Activity
I'm re-tagging this v0.3 as it's a behavior change, not a bugfix
Sounds good to me
Triaging this to 0.4, I won't have time to work on this in the near future.
@keflavich - I'm thinking we should change
interpolate_nanto True by default inconvolve_fftto make it consistent withconvolve. Do you have any objections to this?No, that seems reasonable and should be the expected response.
@keflavich - I'm having issues understanding how to make
convolve_fftbehave likeconvolve:import numpy as np image = np.random.random((128, 128)) image[30:60,40:90] = np.nan subplot(221) imshow(image) title('Original') from astropy.convolution import convolve_fft, Gaussian2DKernel, convolve kernel = Gaussian2DKernel(10) new1 = convolve(image, kernel) subplot(222) imshow(new1) title('Classic') new2 = convolve_fft(image, kernel) subplot(223) imshow(new2) title('FFT interpolate_nan=False') new3 = convolve_fft(image, kernel, interpolate_nan=True) subplot(224) imshow(new3) title('FFT interpolate_nan=True')gives:
@keflavich - any idea what is going on? The
interpolate_nanargument doesn't seem to have any effect.This looks like an error; will investigate
@keflavich - Ah, it looks like I need to set
normalize_kernel=Trueto get it to work.So I guess the question is, since
convolvedoesn't require the kernel to be normalized, is there a way that we can makeconvolve_fftalso not care?Valid question. I'm not 100% sure my reasoning for including the check for kernel normalization still holds; I think that is likely a holdover from old problems. keflavich@2e3450a removes it. Need to run the tests, but if you use that patch, it gets to "almost consistent" with convolve (subtle differences still visible)
Tests do not pass. I'll have to go back and check the math... maybe this afternoon.
5 remaining items
Removing the milestone since I don't have time to work on this and I don't think anyone else does either. There is no bug here, so not critical.
Hi everyone! I'd like to pitch in on this and the related #1647. I was delighted to find out that astropy could provide convolution of numpy arrays with missing values, but was then somewhat disappointed upon inspecting the implementation and finding that it performs a separate interpolation step before the convolution. Replacing NaNs with "fabricated" data using interpolation (coincidentally using the convolution kernel as interpolation weight) is not the same as treating NaNs as missing data.
Consistent handling of missing data in convolution/filtering can be done e.g. by the procedure of Normalized Convolution, which looks a lot like what the astropy convolution routines are trying to do. In fact, if the interpolation step was taken out, and the
if not npy_isnan(fixed[i, j])check removed from the "proper convolution" step, the convolve routines would be performing Normalized Convolution straight out of the book, as far as I can see.For a few intuitive and non-rigorous explanations of Normalized Convolution as applied to missing data, see http://www.mathworks.com/matlabcentral/fileexchange/41961-nanconv, http://www.scribd.com/doc/204528731/Normalized-convolution-for-image-interpolation, and http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PIRODDI1/NormConv/NormConv.html.
In general, one can write the procedure as
h[n] = sum(c[n - k] * f[n - k] * g[k]) / sum(c[n - k] * g[k]),
where f is the signal, g is the kernel, and c is a weight function which is 0 at missing data and 1 otherwise (in some applications, c can take any value in [0, 1] and be interpreted as the level of certainty about each value, but I don't see the relevance of this for astropy at the moment). As mentioned above, this is equivalent to how astropy performs the main convolution step in the convolve routines. I haven't looked at convolve_fft, but to implement identical behavior, one simply needs to realize that the expression above is the elementwise quotient of two ordinary discrete convolutions. In convolve_fft, the two convolutions could be computed separately in Fourier space, and the quotient taken afterwards (of course, there may exist a more efficient algorithm that I don't know about yet).As impied by the name, and also implicitly pointed out above and in #1647, Normalized Convolution doesn't make sense if one does not care about the normalization of the kernel. The procedure entails a separate renormalization of the kernel for each element in the output, which is of course the wrong thing to do if the convolution is just viewed as an unnormalized wieghted sum, and is particularly nonsensical when using a zero-sum kernel that cannot be normalized even in principle.
I would be happy to spend some time trying to adjust the behavior of the convolution routines such that they handle missing data in a consistent way using Normalized Convolution, while working properly with unnormalized/zero-sum kernels. However, I'd like to get some feedback from you guys before starting on this. Does this sound interesting to you? Are such modifications, which will neccesarily lead to some changes in behavior, likely to be accepted by the project? How to avoid breaking backward compatibility too badly? What should be the expected behavior in various cases?
@danielwe - Normalized convolution sounds like what I was going for in
convolve_fft. A great first step toward resolving this issue (consistency) and makingconvolvesupport normalized convolution would be to add a test case based off of, e.g., the matlab implementation. Would you be willing to do that?@keflavich:
I had a look at the convolution_consistency branch in your fork now. We're definitely thinking about the same way of handling NaNs. What's the status here -- do you think that your version could work, and only needs an appropriate test in test_convolve_fft to prove it? I'm not very familiar with these kinds of tests, but I could probably have a look at that. However, I still see a few things I'd suggest to change in the implementation: the handling of the edge behavior options could be simplified, the interpolate_nan and ignore_edge_zeros options could be made redundant (but remain for backward compatibility), etc.I'd also like to work on the convolve function, which should be very straightforward to adjust such that it behaves sensibly in all cases. This would involve only handling NaNs as missing when normalize_kernel=True, and otherwise setting NaN to 0 and computing naïvely (one could also add an option to treat NaN as 0 even if normalize_kernel=True). In some sense I think the natural first step here is make sure convolve does whatever we come to agree is the right thing to do in all cases, since it is the simpler function both in terms of available options and implementation technicalities. When this is achieved, convolve_fft can be modified for consistency with convolve.
Anyway, I'll probably make a fork and start tinkering this weekend if I find the time. If you can be a little more specific about what you would like me to do about the test case, I can look into that.
@danielwe:
Yes, I think my version can work. I don't think theinterpolate_nanandignore_edge_zerosare entirely redundant, but it's been a few years since I looks closely at that problem and maybe I missed something.It would be great if you could work on the
convolvefunction. The important thing is starting from an agreed-upon "right" behavior that should be matched betweenconvolve_fftandconvolvefor, e.g., the Lena example from the docs you linked. So you'd start by making a test that loads the undersampled version of the image and the "correctly" convolved image (say, from matlab), and check that bothconvolveandconvolve_fftreproduce it to some tolerance (withnp.testing.assert_allclose, for example).Normalized convolution would be a nice addition. @keflavich and @danielwe - are either of you still working on this?
@larrybradley - based on previous discussion, I'm fairly sure that
convolve_fftdoes implement "Normalized Convolution", but it's not as obvious whetherconvolvedoes. @astrofrog - let's chat about this briefly tomorrow.I thought that
convolveis interpolating over thenp.nans (?).I'd also like to see
convolvegeneralized such that it could take amaskfor missing/bad data (instead of treating onlynp.nanas missing).Hi,
While writing my earlier posts, I realized I can obtain normalized convolution by doing something like
convolve(data.filled(0.0), kernel) / convolve((~data.mask).astype(float), kernel)wheredatais a masked numpy array where the mask marks missing data,kernelis whatever you would like to convolvedatawith, andconvolveis any standard convolution routine using zero-padding.I still plan to keep my promise and look into the astropy implementations, but now that I don't urgently need the functionality myself I don't know when I'll get the time. I think the current status is that
convolve_ffttries to implement normalized convolution, but has some issues and needs more/better tests (tests should be trivial, given the relation above), whileconvolvejust interpolates.I also agree that the functions ought to take masked arrays.
At the moment, convolve and convolve_fft have slightly different behaviors - for instance:
There may be other differences too. We should try and make the behavior more consistent.