Fixed Gaussian2D model

[larrybradley]

These are fixes to the Gaussian2D model. Specifically, the b parameter was incorrect in the code (but correct in the docs). Also, the a parameter was incorrect in the docs (but correct in the code). This PR also updates the derivatives based on the correct b value. The other changes represent some refactoring to hopefully make things a little clearer/cleaner.

[embray]

Kind of annoying, but I wonder if it might also be worth saving x - x_mean and y - y_mean in temporary variable since they're computed twice.

[nden]

@larrybradley Looks good to me. Just one minor comment. To be consistent with the rest of the code could you insert space between operators, things like np.cos(theta) ** 2 . Thanks!

[nden]

Should we consider this a bug fix and tag it 0.3.1 or is it too late for that?

[embray]

I don't think it's too late. I'm still a little stuck on 0.3.1 in getting all the tests to pass so the release has been delayed a bit by that. This should merge uncontroversially. Please just include a changelog entry under the 0.3.1 section.

[cdeil]

@larrybradley Can you add a unit test that would fail without your fix but now succeeds? (For me it's hard seeing which formulas are correct from just looking at the code.)

I think the travis-ci test failure in astropy.io.fits is unrelated to this pull request and might go away if someone re-starts the test or if you push a new commit to trigger a re-test.

[larrybradley]

@cdeil I've added the test.

[cdeil]

Thanks for adding the test!

With the hardcoded result numbers you have its still hard to see if the model now is correct, right?

An alternative (or addition) could be to test that a rotated model gives the same result at a rotated point ... something like this (untested):

from numpy.testing import assert_allclose
from astropy.coordinates import Angle
from astropy.modeling.models import MatrixRotation2D
from astropy.modeling.models import Gaussian2D

amplitude = 42
x_mean, y_mean = 43, 44
x_stddev, y_stddev = 2, 3
theta = Angle(10, 'deg')

pars = dict(amplitude=amplitude, x_mean=x_mean, y_mean=y_mean,
            x_stddev=x_stddev, y_stddev=y_stddev)
rotation_matrix = MatrixRotation2D(angle=theta.degree)

point1 = (x_mean + 2 * x_stddev, y_mean + 2 * y_stddev)
point2 = rotation_matrix(*point1)

g1 = Gaussian2D(theta=0, **pars)
g2 = Gaussian2D(theta=theta.radian, **pars)

value1 = g1(*point1)
value2 = g2(*point2)

assert_allclose(value1, value2)

Actually this could be easily extended to test all models that can be rotated and a few rotation angles (e.g. a negative one, and one larger than 360 deg).

@larrybradley If you don't want to add this here, I'm 👍 to merge this now and I can make a new pull request with this more generic test for rotated models.

[perrygreenfield]

I'd say not as an alternative, perhaps in addition to, though it is hard to see how it adds any real new test to the existing one. As an alternative, it is vulnerable to a change to the function that makes it not a correct 2D gaussian. It only tests that the rotation is applied correctly to the functional form. On the other hand, if the rotation is being computed incorrectly, the existing test is very unlikely to get the right answer, and it ensures that the functional form is right at the same time. Whether the answers in the test are obviously right is a different matter, but that's what external validation of the test is for.

[larrybradley]

I've added an additional, alternative test. Thanks, @cdeil. But note that since MatrixRotation2D rotates about the origin, the Gaussian2D x_mean and y_mean must be 0.0 for this test (i.e. both the 2D Gaussian and test point are then rotated about the origin).

Thanks, @perrygreenfield. Actually the way the b cross term was previously incorrectly calculated caused the 2D Gaussian to change its shape/ellipticity as a function of rotation angle (and the rotation angle was in the wrong direction). This alternative rotation test fails in that case. It ensures that at least the shape of Gaussian isn't changing as a function of rotation, so it's a good alternative test. But yes, bugs in the MatrixRotation2D rotation would also cause this new test to fail.

[cdeil]

@perrygreenfield Agreed that the test I suggested only tests the rotation part of the model evaluation. But I believe the other parts were already covered by this test.

Concerning the point of "external validation" you mention, here is a check of @larrybradley's numbers against Sherpa:

import numpy as np
from sherpa.models import Gauss2D
sigma_to_fwhm = np.sqrt(8 * np.log(2))

g = Gauss2D()
g.ampl = 100
g.xpos = 1.7
g.ypos = 3.1
g.fwhm = sigma_to_fwhm * 5.0
g.ellip = 1 - 3.3 / 5.0
g.theta = np.pi/6.
g.integrate = False

x, y = np.mgrid[0:5, 0:5]
g(x.flatten(), y.flatten()).reshape(x.shape)

I get these numbers which are different:

array([[ 71.3747051 ,  84.33994996,  92.1020545 ,  92.95061421,  86.69262967],
       [ 70.9461993 ,  85.73626166,  95.75179377,  98.82714765,  94.26545545],
       [ 66.88289735,  82.66028215,  94.41166469,  99.6555279 ,  97.21294495],
       [ 59.80013416,  75.58408495,  88.28878253,  95.3076285 ,  95.08166221],
       [ 50.70962213,  65.54883422,  78.30446686,  86.44802358,  88.20040255]])

Probably I made a mistake ... I'll sleep on it and have a fresh look tomorrow.
Still, @larrybradley, did you get your reference numbers from some external implementation such as scipy.stats.multivariate_normal.pdf or did you simply paste numbers from your updated formulas in astropy.modeling.Gaussian2D?

[cdeil]

In any case I'm wondering if "hard-coding" the rotation and stretching into the Gaussian2D is the best approach, or if maybe factoring that part out and re-using it for other 2D models would be better?

At the moment Gaussian2D is the only model that has this ellipticity in astropy, but as can be seen e.g. by the beta2d and gauss2d and a few more models in Sherpa, there are several 2D models where this will probably be added also in astropy, and implementing and testing it once would IMO be preferable.

[larrybradley]

@cdeil I quickly noticed that I exactly reproduce your Sherpa results if I use theta = np.pi/3. instead of theta = np.pi/6. (i.e. 90 deg - theta) in my test. I will take a look at this further tomorrow.

[cdeil]

@larrybradley OK, so one of us is confused about (x, y) versus (y, x) and / or the direction of rotation theta to get this effect. It's probably me, because I don't have time now to sit down and think through the formulas in the docstrings of Sherpa and astropy.

One thought I had: Shouldn't it be

y, x = np.mgrid[0:5, 0:5]

instead of

x, y = np.mgrid[0:5, 0:5]

because isn't this correct

In [8]: y, x = np.mgrid[0:3, 10:12]

In [9]: x
Out[9]: 
array([[10, 11],
       [10, 11],
       [10, 11]])

In [10]: y
Out[10]: 
array([[0, 0],
       [1, 1],
       [2, 2]])

In [11]: x.shape
Out[11]: (3, 2)

In [12]: x[0, 0]
Out[12]: 10

?

What would help IMO would be a description of the axis order and rotation angle sense with a few examples of how they relate to the numpy axis order in the astropy docs ... like what @kbarbary proposes here: astropy/photutils#32

[astrofrog]

@cdeil - I think the convention should probably be that the orientation should be 'correct' when output to a FITS file or displayed in a matplotlib image (which I think is the same).

[cdeil]

This should help understand what is going on (I used this branch to run this example)
http://nbviewer.ipython.org/gist/cdeil/9134816

My conclusion:

• The current Gaussian2D.eval as implemented in this PR is correct.
• The current test_Gaussian2D should be changed to call y, x = np.mgrid[0:5, 0:5] instead of x, y = np.mgrid[0:5, 0:5] and the hard-coded matrix g transposed.

This (x, y) versus (y, x) axis order "bug" (or call it axis order convention difference) needs to be fixed in ~ 10 other places in astropy and documented.
First of all: @larrybradley , @perrygreenfield , @astrofrog Do you agree?
If yes: @larrybradley Do you want to fix this in this PR or should we merge this now and you or I make a new PR?

[cdeil]

@larrybradley This needs to be rebased against master, the deriv model method has been renamed recently to fit_deriv.

[cdeil]

I was trying to find out if the deriv method is covered by the unit tests.

Turns out it is, but the only way I was able to find this coverage information for this pull request was via the travis-ci log, which seems overly complicated:
https://travis-ci.org/astropy/astropy/jobs/18914850#L1843
https://coveralls.io/files/137268781#L229
@astrofrog Is there a better way of browsing to coverage information for astropy pull requests?

[astrofrog]

@cdeil - you can run:

python setup.py test --coverage

or go directly to https://coveralls.io/r/astropy/astropy

[cdeil]

But the https://coveralls.io/r/astropy/astropy page just has the builds in order, so it's very hard to find the coverage report for a given pull request if the build isn't very recent, no?
(Am I missing something?)

[astrofrog]

Ah right yes - in that case there is no easy way. The issue numbers are listed on coveralls though, but I agree it'd be nice to have an easier way.

[larrybradley]

@cdeil Yes, I agree! I actually suspected last night that the mgrid order needed to be reversed in the test, but I wanted to check it. I will update the test to y, x = np.mgrid[0:5, 0:5] and transpose the test data. I think we should leave the other 10 axis order fixes in astropy to a separate, self-contained, PR. Thanks!

[larrybradley]

I've updated the test and rebased on the current master. Hopefully this is ready to be merged now.

[cdeil]

👍 to merge this now.