default template fitter is no running

[Former user]

Dear Pyspeckit team,

I have tried to use your default template fitter that you have as example in your web page but when I run it I get an error that seems to be related with the input parameter length. Below are the final lines of the error message but in the attached files you can find it complete in the file error_get. You can also find my code, the spectrum to fit and the template file.

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pyspeckit-0.1.16.dev-py2.7.egg/pyspeckit/spectrum/models/model.pyc in L(x) 292 lower_parind = jjself.npars+self.vheight 293 upper_parind = (jj+1)self.npars+self.vheight --> 294 v += self.modelfunc(x, parvals[lower_parind:upper_parind], *kwargs) 295 return v 296 return L

TypeError: 'NoneType' object is not callable

I will appreciate if you can help me with this issue

[Adam Ginsburg]

Thanks for the bug report. I'm looking into it.

[Adam Ginsburg]

OK, there were a few bugs affecting the template fitter. Your example should now run without too much difficulty using this commit: 53035c984cc5bee7536456b901ba5ffb89263a56

However, from a quick examination of your example, I don't expect it to work very well - the template looks like a bad match to the spectrum.

I should also warn you that you're the first serious user of the template fitter, so you may run into further issues. Please post additional issues as you encounter them.

[Julián Mejía-Restrepo]

Hi Adam,

Thank you very much for your help. Sorry for the delay in the answer, I was on vacations for two weeks and I did not see this reply until today.

The template fitter is now running well but as you said the fit is not quite well. I have been playing with the code and I realized that the scaling of the template is very well performed by the fitter. However, I figured out that the template fitter is not shifting appropriately.

I reach to this conclusion after manually shifting the template by 2.0, 5.0 and 10.0A and refitting and I have always get a shift value close to 0A and not much larger tha 1.0A. I then tried setting proper limits to the fitter:

dataspec.specfit(fittype='template',xmin=xmin_fe,xmax=xmax_fe,guesses=[1,0],limits=[(0.1,10.0),(-15.0,15.0)],limited=[(True,True),(True,True)])

But I get the same results: I very good scaling but very bad shifting (and even the same value that I get without setting limits). It seems to me that the code internally is not giving enough freedom to the shift parameter. Do you have any idea about the reason for this issue?

[Adam Ginsburg]

Your issue may be because of how the edges of the template are treated. If your template is smaller (covers a smaller wavelength range) than your data, the left and right side of the template are set to zero. If this is the case, you could try cropping your spectrum to only the part that is within the template's range. If not, we'll have to look in more detail.

[Julián Mejía-Restrepo]

Hi Adam, I tried your suggestion of cropping the spectrum to the template size and I still get the same results. For some reason the code do not give enough freedom to the shift parameter and the output result is always close to the initial guess given to the shift parameter .

[Adam Ginsburg]

Hi Julian, This seems to be a challenging problem, one that we might be addressing with the help of a student in a few months. In the meantime, I've done some diagnostic work:

# diagnostics
import pylab as pl
import numpy as np
shifts = np.arange(-100, 100, 1)
chi2 = np.array([((dataspec.data - template_fitter.modelfunc(dataspec.xarr,
                                                             1.18, shift,
                                                             xshift_units='angstroms'))**2).sum()
                 for shift in shifts])

dataspec_cropped = dataspec.slice(xmin_fe, xmax_fe, units='angstroms')
chi2_c = np.array([((dataspec_cropped.data -
                     template_fitter.modelfunc(dataspec_cropped.xarr, 1.18,
                                               shift,
                                               xshift_units='angstroms'))**2).sum()
                   for shift in shifts])

dataspec_flagged = dataspec_cropped.copy()
a,b = dataspec_flagged.xarr.x_to_pix(2766),dataspec_flagged.xarr.x_to_pix(2833)
dataspec_flagged.data[a:b] = np.nan
chi2_f = np.array([np.nansum((dataspec_flagged.data -
                              template_fitter.modelfunc(dataspec_flagged.xarr,
                                                        1.18, shift,
                                                        xshift_units='angstroms'))**2)
                   for shift in shifts])

pl.figure(3)
pl.clf()
pl.plot(shifts, chi2/len(dataspec), label='raw')
pl.plot(shifts, chi2_c/len(dataspec_cropped), label='cropped')
pl.plot(shifts, chi2_f/len(dataspec_flagged-abs(b-a)), label='flagged')
pl.legend(loc='best')


dataspec.specfit.selectregion(xmin=xmin_fe, xmax=xmax_fe, exclude=[2766, 2833])
dataspec.plotter(xmin=xmin_fe-25, xmax=xmax_fe+25)
dataspec.specfit.highlight_fitregion()
dataspec.specfit(fittype='template', xmin=xmin_fe, xmax=xmax_fe, guesses=[1,-15],
                 limits=[(0.1, 10.0), (-150.0, 150.0)], limited=[(True, True),
                                                                 (True, True)],
                 exclude=[2766, 2833])
dataspec.specfit.highlight_fitregion()

The best I can come up with now is that either the cost function is being incorrectly computed (perhaps the chi^2 value is ignoring the edges or including too much of the edges) or the step size within mpfit needs to be manually adjusted. I'll keep looking, but I'd also appreciate any help you can offer since template fitting is somewhat outside my expertise.

[Julián Mejía-Restrepo]

Hi Adam,

Thank you very much for this diagnostic. I just checked it and I also checked the chi2 values for different shifts obtained directly by specfit and the distribution is not similar at all to any of the raw, cropped and flagged approaches.

chi2_specfit=np.zeros_like(shifts)
dataspec.specfit.selectregion(xmin=xmin_fe, xmax=xmax_fe, exclude=[2766, 2833])
dataspec.plotter(xmin=xmin_fe-25, xmax=xmax_fe+25)
dataspec.specfit.highlight_fitregion()
i=0
for shift in shifts:

    dataspec.specfit(fittype='template', xmin=xmin_fe, xmax=xmax_fe, guesses=[1,shift],
                     limits=[(0.1, 10.0), (shift-0.5, shift+0.5)], limited=[(True, True),
                                                                     (True, True)],
                     exclude=[2766, 2833])
    #dataspec.specfit.highlight_fitregqion()
    chi2_specfit[i]=dataspec.specfit.chi2/6000
    i=i+1

pl.plot(shifts,chi2_specfit,label='specfit')
pl.legend(loc='best')

I then think that the problem is more related, as you suggested, to an incorrect computation of the chi2 because of problems at the edges maybe because of the drastical drop of the template to zero at the edges.

Which is (are) there any particular routine(s) to check the way the chi2 is computed?

[Adam Ginsburg]

chi^2 isn't directly computed, exactly. The optimizers, mpfit and lmfit, use the returned residuals scaled by the error, then take the sum of squares themselves.

The residuals are computed in this line: https://bitbucket.org/pyspeckit/pyspeckit/src/13facf82258b942e496f5577976276627257a7cc/pyspeckit/spectrum/models/model.py?at=master#cl-306

The tricky thing is, it is not possible to modify the length of the spectrum or the value of the spectrum or error during the fitting process. Doing so doesn't make sense in the context of the optimization algorithm anyway. I suppose we could try to have the template use the spectrum data outside of its own range such that the residual is exactly zero, but that isn't trivial.