Fitting of models from GUI does not seem to be the same as the Astropy machinery

[eteq]

This may be a dupe #179 or #180, but I'm not sure, so I'm creating a new issue. @nmearl, feel free to close this immediately if you think it's the same underlying cause.

I noted strange/poor behavior with the model fitter in specviz, and thought I would try to track down what's going on by first trying to reproduce the behavior in a python terminal. However, it seems that something is critically different between how the models are evaluated in the way I use astropy.modeling and the way specviz seems to run. I'll illustrate the test case here and maybe we can get to the bottom of it by having @pllim or @nmearl tell me if there's some assumption I've made that's wrong.

I started by loading the example_spec_ascii.txt file in specviz, which came from a Trello card (I can try to find exactly which one if it seems important). I then selected a region using the span tool. I then set up a model with a Gaussian and a Const, with the arithmetic behavior Gaussian1+Const1. I hit "create layer", tweaked the defaults a bit to get a reasonably close answer (using "Update Layer" as necessary - although see #181 which came out of this), and then did "Perform fit". To my surprise, the fit ended up way off:

To try to diagnose this, I opened up a python session and tried to fit the data directly in Astropy:

from astropy import table
from astropy.modeling import models,fitting

t=table.Table.read('example_spec_ascii.txt',format='ascii')
x=t['wavelength']
f=t['flux']

# these are the same values as what I put in by hand into the specviz GUI
c=models.Const1D(amplitude=3e-15)
g=models.Gaussian1D(amplitude=9e-16,mean=1.74,stddev=0.1)
fr=fitting.LevMarLSQFitter()
fitmodel=fr(c+g, x, f)

from matplotlib import pyplot as plt
plt.plot(x,f)
msk=(x<2)&(x>1.5)  # this is roughly the fitting region in specviz as the screenshot above shows
plt.plot(x[msk],f[msk])
plt.plot(x, fitmodel(x))

To my even greater surprise I saw this:

The red line is the as-fit model. It appears to be exactly what I would have thought the fit should do, and was initially surprised to see specviz did not manage to do.

So it appears something in specviz is not using the astropy fitting machinery the way I would have thought. Any ideas what's going on?

[eteq]

While this is a shot in the dark, my first hunch is that somehow the gaussian and constant models are being fit separately instead of the composite model being fit. That would explain why the specviz fit is both too wide and too tall.

[nmearl]

Fitting uses the composite model.

Here is what I get when I run your script without any alterations:

Another thing I'd like to point out is that, in SpecViz, the data is masked before it's fed to the fitter. Making such changes to your script:

%matplotlib notebook

from astropy import table
from astropy.modeling import models,fitting

t=table.Table.read('/Users/nearl/Desktop/example_spec_ascii.txt',format='ascii')
x=t['wavelength']
f=t['flux']
msk=(x<2)&(x>1.5)  # this is roughly the fitting region in specviz as the screenshot above shows
x = x[msk]
f = f[msk]

# these are the same values as what I put in by hand into the specviz GUI
c=models.Const1D(amplitude=3e-15)
g=models.Gaussian1D(amplitude=9e-16,mean=1.74,stddev=0.1)
fr=fitting.LevMarLSQFitter()
fitmodel=fr(c+g, x, f)

from matplotlib import pyplot as plt
plt.plot(x,f)
plt.plot(x, fitmodel(x))

Which gives me:

If I feed your initial values in SpecViz and let the fitter try, I get:

Honestly, this feels more like the fitter is just really sensitive?

[pllim]

I am no help in the department, but perhaps @ibusko can comment further.

[ibusko]

I suspect the issue is in how specviz handles the masks when fitting. If I get the same data but extract just the region around the line (~ 1 to 2) I get a nice fit after just adding a Const + Gaussian model and running L-M about 4-5 times. The standard L-M call runs at most 100 iterations,a s far as I remember. So I am just saying that it needs ~400-500 iterations to converge.

[ibusko]

I didn't explain right: I extract the region of the spectrum into a new ascii file, then ingest that file. So no need to rely on ROIs.

[eteq]

@nmearl - oops, my bad in the initial issue report. I copied-and-pasted the wrong line there. Indeed you're right that the fitter is supposed to use the masked version of the wavelength/flux arrays. So in the original issue if you change the line fitmodel=fr(c+g, x, f) to fitmodel=fr(c+g, x[msk], f[msk]) (and move the msk= line above that of course), you get the plot I showed. So the corrected version is:

from astropy import table
from astropy.modeling import models,fitting

t=table.Table.read('example_spec_ascii.txt',format='ascii')
x=t['wavelength']
f=t['flux']
msk=(x<2)&(x>1.5)  # this is roughly the fitting region in specviz as the screenshot above shows

# these are the same values as what I put in by hand into the specviz GUI
c=models.Const1D(amplitude=3e-15)
g=models.Gaussian1D(amplitude=9e-16,mean=1.74,stddev=0.1)
fr=fitting.LevMarLSQFitter()
fitmodel=fr(c+g, x[msk], f[msk])

from matplotlib import pyplot as plt
plt.plot(x,f)
plt.plot(x[msk],f[msk])
plt.plot(x, fitmodel(x))

Which I think is functionally equivalent to what you did above.

[eteq]

Hmm, curious. I wonder if maybe I cut out a different size layer than I initially thought? When I do things in specviz in a more careful manner, I see that I can often get a decent fit if I just hit the "Perform fit" button a few times, occasionally tweaking the parameters a little towards sanity. That suggests that @ibusko's suggestion is the key: turning up the allowed number of iterations...? If so, this can probably be closed as a red herring and a more specific issue/PR created.

[ibusko]

The astropy fitters can be called with a keyword parameter that controls the maximum number of iterations. Maybe it would be useful to expose that parameter in the UI. They have other calling parameters as well, enough to give a good deal of control on the way they work. But, regardless, in the way the fitter-model-UI interaction works right now, just by re-doing the fit in sequence multiple times we can achieve the same result as passing a larges number of iterations in a single call.

[eteq]

Hmm, so maybe the short-term fix is to just note this in the documentation somewhere? That feels a little jankety, but it's also true that non-linear fitting by its nature is rather heuristic and "well lets just see what happens"-y...

[ibusko]

Indeed, documenting this would be nice.

These are the parameters accepted by astropy fitters. Maybe exposing 'maxiter' and 'acc' would be fine; I'm not sure about the others.

    x : array
       input coordinates
    y : array
       input coordinates
    z : array (optional)
       input coordinates
    weights : array (optional)
       weights
    maxiter : int
        maximum number of iterations
    acc : float
        Relative error desired in the approximate solution
    epsilon : float
        A suitable step length for the forward-difference
        approximation of the Jacobian (if model.fjac=None). If
        epsfcn is less than the machine precision, it is
        assumed that the relative errors in the functions are
        of the order of the machine precision.
    estimate_jacobian : bool
        If False (default) and if the model has a fit_deriv method,
        it will be used. Otherwise the Jacobian will be estimated.
        If True, the Jacobian will be estimated in any case.

[eteq]

epsilon could be useful too, but is rather model-specific in a way that acc and maxiter aren't. Maybe the thing to do is expose this stuff in an "advanced" panel option (I think @hcferguson advocated for this at one point). That's certainly not a high-priority item right now, though.

@nmearl or @pllim - if you like, I (or either of you) can close this and replace it with more targeted issues like "document tips for the fitter" and/or "add advanced options", as the initial title/topic doesn't really seem correct anymore.

[pllim]

To be continued in #121 and #184.