Skip to content

models: bounds / tied issues #1064

New issue
New issue
Closed
Closed
models: bounds / tied issues#1064
Labels
modeling
Milestone
v0.3.0

Description

@adrn
adrn
opened on May 7, 2013

Writing this down so I don't forget...

I started using astropy.models in my own code, but have run in to some issues. First off, we need to sprint hard on the documentation -- no one will use this if we don't have good documentation explaining all of the nice features in here, which are currently buried in code!

I'm having trouble getting bounded / tied parameters to work with the Gaussian2DModel. Haven't tried this with others, and haven't yet tried to dig in to the code, but here's an example of what I'm trying to do -- fit a 2D Gaussian to some image data:

gauss = am.builtin_models.Gaussian2DModel(amplitude=10., x_mean=5., y_mean=5.,
                                      x_stddev=4., y_stddev=4., theta=0.5,
                                      bounds={"x_mean" : [0.,11.],
                                              "y_mean" : [0.,11.],
                                              "x_stddev" : [1.,4.],
                                              "y_stddev" : [1.,4.]})
gauss_fit = am.fitting.NonLinearLSQFitter(gauss)

but when I run the code, the parameters I get are definitely not bounded:

Model: Gaussian2DModel
    Dim:   2
    Degree: N/A
    Parameter sets: 1
    Parameters:
               amplitude: [1620.8534792444862]
               x_stddev: [71.600369674585707]
               y_stddev: [4.023002873594919]
               x_mean: [-7.4567104062070975]
               y_mean: [-17.633779415014097]
               theta: [4.204525599803083]

Also tried this with using tied:

gauss = am.builtin_models.Gaussian2DModel(amplitude=10., x_mean=5., y_mean=5.,
                                      x_stddev=4., y_stddev=4., theta=0.5)
gauss.x_stddev.tied = lambda g: g.y_stddev

but also does not seem to work.

cc @nden

Activity

  1. nden commented on May 9, 2013

    @nden
    nden
    on May 9, 2013
    Contributor

    @adrn Could you fetch upstream and try your example again. When given in the initializer the constraints were not propagated. This should be fixed by #1075.
    I agree about the documentation. May be this can become an example of constraint fitting?
    I am not closing this until I hear it works.

  2. nden commented on May 9, 2013

    @nden
    nden
    on May 9, 2013
    Contributor

    @adrn actually I missed one commit, so part of this won't work

  3. nden
    added 2 commits that reference this issue on May 16, 2013

    Address problems with constraints reported in astropy#1060 and astrop…

    bfae565

    Merge pull request #1107 from nden/constraints_issues

    dc6885c

  4. eteq commented on May 18, 2013

    @eteq
    eteq
    on May 18, 2013
    Member

    @adrn - was this fixed by #1107 ?

  5. adrn commented on May 19, 2013

    @adrn
    adrn
    on May 19, 2013
    MemberAuthor

    @eteq I don't think so. I fetched from master, but the bounds still don't work -- run this:

    import numpy as np
import astropy.models as am


data = np.array([505.0, 556.0, 630.0, 595.0, 561.0, 553.0, 543.0, 496.0, 460.0, 469.0, 426.0, 518.0, 684.0, 798.0, 830.0, 794.0, 649.0, 706.0, 671.0, 545.0, 479.0, 454.0, 505.0, 700.0, 1058.0, 1231.0, 1325.0, 997.0, 1036.0, 884.0, 610.0, 487.0, 453.0, 527.0, 780.0, 1094.0, 1983.0, 1993.0, 1809.0, 1525.0, 1056.0, 895.0, 604.0, 466.0, 510.0, 678.0, 1130.0, 1986.0, 2670.0, 2535.0, 1878.0, 1450.0, 1200.0, 663.0, 511.0, 474.0, 569.0, 848.0, 1670.0, 2611.0, 3129.0, 2507.0, 1782.0, 1211.0, 723.0, 541.0, 511.0, 518.0, 597.0, 1137.0, 1993.0, 2925.0, 2438.0, 1910.0, 1230.0, 738.0, 506.0, 461.0, 486.0, 597.0, 733.0, 1262.0, 1896.0, 2342.0, 1792.0, 1180.0, 667.0, 482.0, 454.0, 482.0, 504.0, 566.0, 789.0, 1194.0, 1545.0, 1361.0, 933.0, 562.0, 418.0, 463.0, 435.0, 466.0, 528.0, 487.0, 664.0, 799.0, 746.0, 550.0, 478.0, 535.0, 443.0, 416.0, 439.0, 472.0, 472.0, 492.0, 523.0, 569.0, 487.0, 441.0, 428.0])
data = data.reshape(11,11)
X,Y = np.meshgrid(np.arange(11),np.arange(11))

gauss = am.builtin_models.Gaussian2DModel(amplitude=10., x_mean=5., y_mean=5.,
                                          x_stddev=4., y_stddev=4., theta=0.5,
                                          bounds={"x_mean" : [0.,11.],
                                                  "y_mean" : [0.,11.],
                                                  "x_stddev" : [1.,4],
                                                  "y_stddev" : [1.,4]})
gauss_fit = am.fitting.NonLinearLSQFitter(gauss)
gauss_fit(X, Y, data)

print(gauss)

    I get negative mean values and huge standard deviations.

  6. adrn commented on May 19, 2013

    @adrn
    adrn
    on May 19, 2013
    MemberAuthor

    Oh! Oops, didn't notice the new import stuff -- hmm, I wonder why the code above worked...must have a phantom old version of astropy installed alongside the dev version...

    With the new astropy.modeling import, the above kind of works...

  7. adrn commented on May 19, 2013

    @adrn
    adrn
    on May 19, 2013
    MemberAuthor

    Spoke too soon -- looks like the bounds work for the stddev's but not the mean? strange.

  8. adrn commented on May 19, 2013

    @adrn
    adrn
    on May 19, 2013
    MemberAuthor

    @nden @eteq Here's some code that generates simulated data, then uses astropy and scipy.optimize to fit the data. For some reason, astropy.modeling is doing much worse:

    https://gist.github.com/adrn/5606738

  9. cdeil commented on May 19, 2013

    @cdeil
    cdeil
    on May 19, 2013
    Member

    @adrn Note that currently the astropy Gaussian2DModel eval gives incorrect results (see point 1 in #1077).
    Care to fix Gaussian2DModel and add a unit test for eval (and gderiv if possible)?

  10. cdeil commented on May 19, 2013

    @cdeil
    cdeil
    on May 19, 2013
    Member

    I just saw that @astrofrog already fixed the astropy Gaussian models in #1106 .

  11. adrn commented on May 19, 2013

    @adrn
    adrn
    on May 19, 2013
    MemberAuthor

    Hm, there is some deeper issue here. I tried out #1106 on the gist code I posted above, still doesn't work. I'm going to take a look at this...

  12. nden
    added a commit that references this issue on May 20, 2013

    Address problems with constraints reported in astropy#1060 and astrop…

    c42998e

  13. nden commented on May 22, 2013

    @nden
    nden
    on May 22, 2013
    Contributor

    @adrn In the gist the initial values p0 are compared with the output of optimize. If you change the true output to the output of the NonLinearLSQFitter the example works.

    print("true", ["{0:.3f}".format(x) for x in gauss.parameters])

    I will add this to the tests unless you want to add it.
    Were the original problems with bounds fixed?

  14. adrn commented on May 22, 2013

    @adrn
    adrn
    on May 22, 2013
    MemberAuthor

    @nden "true" are the parameters I used to generate the data, not meant to be the NonLinearLSQFitter parameters. I updated the gist code to also print those values, though...still getting weirdness. Here's the code again:
    https://gist.github.com/adrn/5606738

    It seems like the bounds work now, but when I run it with bounds the two stddev's get slammed to the upper bound (4.). If I delete the bounds, the y_stddev value goes to ~271 (true value is 2.).

  15. adrn commented on May 22, 2013

    @adrn
    adrn
    on May 22, 2013
    MemberAuthor

    Here are the results from removing the bounds on the stddev's:

    true ['137.000', '5.100', '5.400', '1.500', '2.000', '0.785']
    leastsq ['137.135', '5.094', '5.410', '1.498', '2.000', '0.786']
    astropy.modeling ['132.801', '5.094', '5.409', '-0.877', '271.430', '0.644']

    Not sure why it's doing so much worse than leastsq...

  16. nden commented on May 22, 2013

    @nden
    nden
    on May 22, 2013
    Contributor

    @adrn I don't think this test compares apples to apples. The gaussian2d function does not evaluate to the same thing as models.Gaussian2DModel. Here's a test which I think does what you intended but let me know if I misunderstood something (I copied the eval method from Gaussian2DModel):

    from astropy.modeling import models, fitting
import numpy as np
from scipy.optimize import leastsq

def test_gaussian2d_fitting():
    def gaussian2d(p, x, y):
        a = 0.5 * ((np.cos(p[5]) / p[3]) ** 2 +
                   (np.sin(p[5]) / p[4]) ** 2)
        b = 0.5 * (np.cos(p[5]) * np.sin(p[5]) *
                   (1. / p[3] ** 2 - 1. / p[4] ** 2))
        c = 0.5 * ((np.sin(p[5]) / p[3]) ** 2 +
                   (np.cos(p[5]) / p[4]) ** 2)

        return p[0] * np.exp(-(a * (x - p[1]) ** 2 +
                               b * (x - p[1]) * (y - p[2]) +
                               c * (y - p[2]) ** 2))

    def error_function(p, x, y, z):
        return z - gaussian2d(p, x, y)

    X,Y = np.meshgrid(np.arange(11),np.arange(11))
    # generate fake data
    p0 = [137., 5.1, 5.4, 1.5, 2., np.pi/4]
    data = gaussian2dm(p0, X.ravel(), Y.ravel()).reshape(X.shape)
    rsn=np.random.RandomState(1234567890)
    n=rsn.randn(data.size)
    n.shape = data.shape
    data += n 

    gauss = models.Gaussian2DModel(amplitude=10., x_mean=5., y_mean=5.,
                               x_stddev=4., y_stddev=4., theta=0.5)
    gauss_fit = fitting.NonLinearLSQFitter(gauss)
    gauss_fit(X, Y, data)
    p_opt, success = leastsq(error_function, 
                         [10., 5., 5., 4., 4., 0.5], 
                         args=(X.ravel(),Y.ravel(),data.ravel()), 
                         maxfev=10000)

    print("true", ["{0:.3f}".format(x) for x in p0])
    print("leastsq", ["{0:.3f}".format(x) for x in p_opt])
    print("modeling", ["{0:.3f}".format(x) for x in gauss.parameters])

    The result of this is:

    true ['137.000', '5.100', '5.400', '1.500', '2.000', '0.785']
leastsq ['137.396', '5.094', '5.391', '1.507', '1.983', '0.769']
modeling ['137.396', '5.094', '5.391', '1.507', '1.983', '0.769']

  17. cdeil commented on Aug 8, 2013

    @cdeil
    cdeil
    on Aug 8, 2013
    Member

    @adrn The Gaussian2DModel and contraints code has been updated in the meantime and is not covered by unit tests. Could you try again and let us know if you still find any issues?

    @nden Running the example from your latest comment I don't get identical results for leastsq and modeling!?

    In [4]: test_gaussian2d_fitting()
('true', ['137.000', '5.100', '5.400', '1.500', '2.000', '0.785'])
('leastsq', ['137.396', '5.094', '5.391', '1.507', '1.983', '0.769'])
('modeling', ['25.809', '4.540', '3.425', '3.261', '2.820', '0.500'])

    (I had to change the line

    data = gaussian2dm(p0, X.ravel(), Y.ravel()).reshape(X.shape)

    to

    data = gaussian2d(p0, X.ravel(), Y.ravel()).reshape(X.shape)

    to make it run.)

  18. adrn commented on Oct 8, 2013

    @adrn
    adrn
    on Oct 8, 2013
    MemberAuthor

    @nden what's the status of this? I'll check with some tests of my own but want to hear from you too..

  19. nden commented on Oct 8, 2013

    @nden
    nden
    on Oct 8, 2013
    Contributor

    @adrn After merging #1534 it will be ready to test.

  20. nden commented on Oct 11, 2013

    @nden
    nden
    on Oct 11, 2013
    Contributor

    @adrn #1534 was merged and I think this is ready to be closed. Can you confirm it works for you?

  21. nden commented on Oct 18, 2013

    @nden
    nden
    on Oct 18, 2013
    Contributor

    @adrn I am closing this as it was included as a test in test_models. Please open a different issue if you think something is wrong.

  22. nden
    closed this as completedon Oct 18, 2013
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    modeling

    Type

    No type

    Projects

    No projects

      Milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions