Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixed Gaussian2D model #2038

Merged
merged 9 commits into from
Feb 21, 2014
Merged

Fixed Gaussian2D model #2038

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
fe8a9be
Fixed Gaussian2D model
larrybradley Feb 3, 2014
f7f8779
Stored several calculations for reuse
larrybradley Feb 3, 2014
2850e9a
insert spaces around all operators for internal consistency
larrybradley Feb 14, 2014
e6f5bda
updated changelog
larrybradley Feb 14, 2014
c67e8b9
updated changelog
larrybradley Feb 14, 2014
ad4565f
Added new Gaussian2D model test
larrybradley Feb 20, 2014
27df53d
doc fixes
larrybradley Feb 20, 2014
7a88c9f
added new Gaussian2D model test
larrybradley Feb 20, 2014
c5011c9
updated Gaussian2D test
larrybradley Feb 21, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 2 additions & 0 deletions CHANGES.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -313,6 +313,8 @@ Bug Fixes

- Raise a `NotImplementedError` when fitting composite models. [#1915]

- Fixed bug in computation of ``Gaussian2D`` model. [#2038]

- ``astropy.nddata``

- ``astropy.sphinx``
Expand Down
163 changes: 89 additions & 74 deletions astropy/modeling/functional_models.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,11 +34,11 @@ class Gaussian1D(Parametric1DModel):
Parameters
----------
amplitude : float
Amplitude of the gaussian
Amplitude of the Gaussian.
mean : float
Mean of the gaussian
Mean of the Gaussian.
stddev : float
Standard deviation of the gaussian
Standard deviation of the Gaussian.

Notes
-----
Expand Down Expand Up @@ -102,17 +102,18 @@ def __init__(self, amplitude, mean, stddev, **constraints):
super(Gaussian1D, self).__init__(param_dim=param_dim,
amplitude=amplitude, mean=mean,
stddev=stddev, **constraints)

@staticmethod
def eval(x, amplitude, mean, stddev):
"""
Model function Gauss1D
Gaussian1D model function.
"""
return amplitude * np.exp(- 0.5 * (x - mean) ** 2 / stddev ** 2)

@staticmethod
def fit_deriv(x, amplitude, mean, stddev):
"""
Model function derivatives Gauss1D
Gaussian1D model function derivatives.
"""

d_amplitude = np.exp(-0.5 / stddev ** 2 * (x - mean) ** 2)
Expand All @@ -128,22 +129,24 @@ class Gaussian2D(Parametric2DModel):
Parameters
----------
amplitude : float
Amplitude of the gaussian
Amplitude of the Gaussian.
x_mean : float
Mean of the gaussian in x
Mean of the Gaussian in x.
y_mean : float
Mean of the gaussian in y
Mean of the Gaussian in y.
x_stddev : float
Standard deviation of the gaussian in x
Either x_fwhm or x_stddev must be specified
Standard deviation of the Gaussian in x.
x_stddev and y_stddev must be specified unless a covariance
matrix (cov_matrix) is input.
y_stddev : float
Standard deviation of the gaussian in y
Either y_fwhm or y_stddev must be specified
theta : float
Rotation angle in radians. Note: increases clockwise.
cov_matrix : ndarray
A 2x2 covariance matrix. If specified, overrides stddev, fwhm, and
theta specification.
Standard deviation of the Gaussian in y.
x_stddev and y_stddev must be specified unless a covariance
matrix (cov_matrix) is input.
theta : float, optional
Rotation angle in radians. The rotation angle increases clockwise.
cov_matrix : ndarray, optional
A 2x2 covariance matrix. If specified, overrides the x_stddev,
y_stddev, and theta specification.

Notes
-----
Expand All @@ -157,14 +160,14 @@ class Gaussian2D(Parametric2DModel):
Using the following definitions:

.. math::
a = \\left(- \\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}} -
\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}}\\right)
a = \\left(\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}} +
\\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}}\\right)

b = \\left(\\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{y}^{2}} -
\\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{x}^{2}}\\right)
b = \\left(\\frac{-\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{x}^{2}} +
\\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{y}^{2}}\\right)

c = \\left(\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}} +
\\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}}\\right)
c = \\left(\\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}} +
\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}}\\right)


See Also
Expand All @@ -183,11 +186,11 @@ def __init__(self, amplitude, x_mean, y_mean, x_stddev=None, y_stddev=None,
theta=0.0, cov_matrix=None, **constraints):
if y_stddev is None and cov_matrix is None:
raise InputParameterError(
"Either y_stddev must be specified, or a "
"Either x/y_stddev must be specified, or a "
"covariance matrix.")
elif x_stddev is None and cov_matrix is None:
raise InputParameterError(
"Either x_stddev must be specified, or a "
"Either x/y_stddev must be specified, or a "
"covariance matrix.")
elif cov_matrix is not None and (x_stddev is not None or
y_stddev is not None):
Expand All @@ -211,54 +214,65 @@ def __init__(self, amplitude, x_mean, y_mean, x_stddev=None, y_stddev=None,
def eval(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta):
"""Two dimensional Gaussian function"""

a = 0.5 * ((np.cos(theta) / x_stddev) ** 2 +
(np.sin(theta) / y_stddev) ** 2)
b = 0.5 * (np.cos(theta) * np.sin(theta) *
(1. / x_stddev ** 2 - 1. / y_stddev ** 2))
c = 0.5 * ((np.sin(theta) / x_stddev) ** 2 +
(np.cos(theta) / y_stddev) ** 2)

return amplitude * np.exp(-(a * (x - x_mean) ** 2 +
b * (x - x_mean) * (y - y_mean) +
c * (y - y_mean) ** 2))
cost2 = np.cos(theta) ** 2
sint2 = np.sin(theta) ** 2
sin2t = np.sin(2. * theta)
xstd2 = x_stddev ** 2
ystd2 = y_stddev ** 2
xdiff = x - x_mean
ydiff = y - y_mean
a = 0.5 * ((cost2 / xstd2) + (sint2 / ystd2))
b = 0.5 * (-(sin2t / xstd2) + (sin2t / ystd2))
c = 0.5 * ((sint2 / xstd2) + (cost2 / ystd2))
return amplitude * np.exp(-((a * xdiff ** 2) + (b * xdiff * ydiff) +
(c * ydiff ** 2)))

@staticmethod
def fit_deriv(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta):
"""Two dimensional Gaussian function derivative with respect to parameters"""

# Helper quantities
a = 0.5 * ((np.cos(theta) / x_stddev) ** 2 +
(np.sin(theta) / y_stddev) ** 2)
b = 0.5 * (np.cos(theta) * np.sin(theta) *
(1. / x_stddev ** 2 - 1. / y_stddev ** 2))
c = 0.5 * ((np.sin(theta) / x_stddev) ** 2 +
(np.cos(theta) / y_stddev) ** 2)

da_dtheta = np.sin(theta) * np.cos(theta) * (1/y_stddev**2 - 1/x_stddev**2)
da_dx_stddev = -np.cos(theta)**2/x_stddev**3
da_dy_stddev = -np.sin(theta)**2/y_stddev**3
db_dtheta = 0.5*np.cos(2*theta) * (1/x_stddev**2 - 1/y_stddev**2)
db_dx_stddev = -np.cos(theta)*np.sin(theta)/x_stddev**3
db_dy_stddev = np.cos(theta)*np.sin(theta)/y_stddev**3
dc_dtheta = np.cos(theta)*np.sin(theta)*(1/x_stddev**2 - 1/y_stddev**2)
dc_dx_stddev = -np.sin(theta)**2/x_stddev**3
dc_dy_stddev = -np.cos(theta)**2/y_stddev**3

d_A = np.exp(- a * (x - x_mean) ** 2
- b * (x - x_mean) * (y - y_mean)
- c * (y - y_mean) ** 2)
d_theta = amplitude * (-(x - x_mean) ** 2 * da_dtheta - (x - x_mean) *
(y - y_mean) * db_dtheta - (y - y_mean) ** 2 * dc_dtheta) * d_A

d_x_stddev = amplitude * (-(x - x_mean) ** 2 *da_dx_stddev -
(x - x_mean) * (y - y_mean) * db_dx_stddev
- (y - y_mean) ** 2 * dc_dx_stddev) * d_A
d_y_stddev = amplitude * (-(x - x_mean) ** 2 *da_dy_stddev -
(x - x_mean) * (y - y_mean) * db_dy_stddev
- (y - y_mean) ** 2 * dc_dy_stddev) * d_A
d_y_mean = amplitude * (+(x - x_mean) * b + 2 * (y - y_mean) * c) * d_A
d_x_mean = amplitude * ((y - y_mean) * b + 2 * (x - x_mean) * a) * d_A
return [d_A, d_x_mean, d_y_mean, d_x_stddev, d_y_stddev, d_theta]
cost = np.cos(theta)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@cdeil - you can run:

python setup.py test --coverage

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

sint = np.sin(theta)
cost2 = np.cos(theta) ** 2
sint2 = np.sin(theta) ** 2
cos2t = np.cos(2. * theta)
sin2t = np.sin(2. * theta)
xstd2 = x_stddev ** 2
ystd2 = y_stddev ** 2
xstd3 = x_stddev ** 3
ystd3 = y_stddev ** 3
xdiff = x - x_mean
ydiff = y - y_mean
xdiff2 = xdiff ** 2
ydiff2 = ydiff ** 2
a = 0.5 * ((cost2 / xstd2) + (sint2 / ystd2))
b = 0.5 * (-(sin2t / xstd2) + (sin2t / ystd2))
c = 0.5 * ((sint2 / xstd2) + (cost2 / ystd2))
g = amplitude * np.exp(-((a * xdiff2) + (b * xdiff * ydiff) +
(c * ydiff2)))
da_dtheta = (sint * cost * ((1. / ystd2) - (1. / xstd2)))
da_dx_stddev = -cost2 / xstd3
da_dy_stddev = -sint2 / ystd3
db_dtheta = (-cos2t / xstd2) + (cos2t / ystd2)
db_dx_stddev = sin2t / xstd3
db_dy_stddev = -sin2t / ystd3
dc_dtheta = -da_dtheta
dc_dx_stddev = -sint2 / xstd3
dc_dy_stddev = -cost2 / ystd3
dg_dA = g / amplitude
dg_dx_mean = g * ((2. * a * xdiff) + (b * ydiff))
dg_dy_mean = g * ((b * xdiff) + (2. * c * ydiff))
dg_dx_stddev = g * (-(da_dx_stddev * xdiff2 +
db_dx_stddev * xdiff * ydiff +
dc_dx_stddev * ydiff2))
dg_dy_stddev = g * (-(da_dy_stddev * xdiff2 +
db_dy_stddev * xdiff * ydiff +
dc_dy_stddev * ydiff2))
dg_dtheta = g * (-(da_dtheta * xdiff2 +
db_dtheta * xdiff * ydiff +
dc_dtheta * ydiff2))
return [dg_dA, dg_dx_mean, dg_dy_mean, dg_dx_stddev, dg_dy_stddev,
dg_dtheta]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ditto here.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, I can do that. I suppose I should also save np.cos(theta) and np.sin(theta) since they are calculated many times.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh yes, definitely those too. At some point I still want to push for more efficient implementations of these functions, but for now a lot is gained simply by reusing intermediate results like those.



class Shift(Model):
Expand Down Expand Up @@ -423,7 +437,7 @@ class Linear1D(Parametric1DModel):

def __init__(self, slope, intercept, **constraints):
super(Linear1D, self).__init__(slope=slope, intercept=intercept,
**constraints)
**constraints)

@staticmethod
def eval(x, slope, intercept):
Expand Down Expand Up @@ -472,7 +486,7 @@ class Lorentz1D(Parametric1DModel):

def __init__(self, amplitude, x_0, fwhm, **constraints):
super(Lorentz1D, self).__init__(amplitude=amplitude, x_0=x_0,
fwhm=fwhm, **constraints)
fwhm=fwhm, **constraints)

@staticmethod
def eval(x, amplitude, x_0, fwhm):
Expand Down Expand Up @@ -738,7 +752,7 @@ def eval(x, amplitude, x_0, width):

return np.select([np.logical_and(x >= x_0 - width / 2.,
x <= x_0 + width / 2.)],
[amplitude], 0)
[amplitude], 0)

@classmethod
def fit_deriv(cls, x, amplitude, x_0, width):
Expand Down Expand Up @@ -801,8 +815,10 @@ def __init__(self, amplitude, x_0, y_0, x_width, y_width, **constraints):
@staticmethod
def eval(x, y, amplitude, x_0, y_0, x_width, y_width):
"""Two dimensional Box model function"""
x_range = np.logical_and(x >= x_0 - x_width / 2., x <= x_0 + x_width / 2.)
y_range = np.logical_and(y >= y_0 - y_width / 2., y <= y_0 + y_width / 2.)
x_range = np.logical_and(x >= x_0 - x_width / 2.,
x <= x_0 + x_width / 2.)
y_range = np.logical_and(y >= y_0 - y_width / 2.,
y <= y_0 + y_width / 2.)
return np.select([np.logical_and(x_range, y_range)], [amplitude], 0)


Expand Down Expand Up @@ -889,7 +905,6 @@ def __init__(self, amplitude, x_0, y_0, R_0, slope, **constraints):
x_0=x_0, y_0=y_0, R_0=R_0,
slope=slope, **constraints)


@staticmethod
def eval(x, y, amplitude, x_0, y_0, R_0, slope):
"""Two dimensional Trapezoid Disk model function"""
Expand Down
36 changes: 36 additions & 0 deletions astropy/modeling/tests/test_functional_models.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,15 @@
import numpy as np
from numpy.testing import assert_allclose
from .. import models
from astropy.coordinates import Angle

try:
from scipy import optimize
HAS_SCIPY = True
except ImportError:
HAS_SCIPY = False


def test_Trapezoid1D():
"""Regression test for
https://github.com/astropy/astropy/issues/1721
Expand All @@ -22,3 +24,37 @@ def test_Trapezoid1D():
yy = model(xx)
yy_ref = [0., 1.41428571, 3.12857143, 4.2, 4.2, 3.12857143, 1.41428571, 0.]
assert_allclose(yy, yy_ref, rtol=0, atol=1e-6)


def test_Gaussian2D():
"""
Test rotated elliptical Gaussian2D model.
https://github.com/astropy/astropy/pull/2038
"""
model = models.Gaussian2D(100, 1.7, 3.1, x_stddev=3.3, y_stddev=5.0,
theta=np.pi/6.)
y, x = np.mgrid[0:5, 0:5]
g = model(x, y)
g_ref = [[77.86787722, 79.8450774, 75.66323816, 66.26262987, 53.6289606],
[86.01743889, 90.20335997, 87.419011, 78.29536082, 64.80568981],
[90.11888627, 96.6492347, 95.79172412, 87.74139348, 74.27249818],
[89.54601432, 98.21442091, 99.55228375, 93.25543692, 80.73169335],
[84.38744554, 94.65710923, 98.12408063, 94.00369635, 83.22642276]]
assert_allclose(g, g_ref, rtol=0, atol=1e-6)


def test_Gaussian2DRotation():
amplitude = 42
x_mean, y_mean = 0, 0
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 = models.MatrixRotation2D(angle=theta.degree)
point1 = (x_mean + 2 * x_stddev, y_mean + 2 * y_stddev)
point2 = rotation_matrix(*point1)
g1 = models.Gaussian2D(theta=0, **pars)
g2 = models.Gaussian2D(theta=theta.radian, **pars)
value1 = g1(*point1)
value2 = g2(*point2)
assert_allclose(value1, value2)