1

I'm trying to reproduce something in python that works correctly in IDL. (This is part of a larger program, so just using IDL is not an option.)

What I'm trying to do is fit a spectrum which has not been wavelength-calibrated, to another which already has, and thereby get a wavelength dispersion solution for the uncalibrated spectrum.

IDL Code:

function fit_wav, x, p, wav=w, flux=f ;our function for mpfitfun
    ;eval a polynomial function for wavelength and interpolate
    ;to the wavelength of the fiducial spectrum
    wav = poly(w, p)
    return, interpol(f, wav, x)
end

function normalize, arr ;scale an array into [0,1]
    return, (arr - min(arr)) / (max(arr) - min(arr))
end

readcol, '~/Desktop/fiducial.txt', wavelength, flux, format='F,F'
readcol, '~/Desktop/to_fit.txt', uncalibrated, format='F'

npix = n_elements(uncalibrated)
pix = [0d:npix-1]
p_init = [1.95, 0.5/npix, 0.]
coef = mpfitfun('fit_wav', wavelength, normalize(flux), 1., $
                p_init, functargs={wav:pix, flux:normalize(uncalibrated)})
print, coef

;plot the results
calib_wav = poly(pix, coef)
cgplot, calib_wav, normalize(uncalibrated), xra=minmax([calib_wav, wavelength])
cgplot, /overplot, wavelength, normalize(flux), color='red'
end

The output for this looks like (fiducial is red, fit is black):

    1.9652186   0.00070232586  -1.5183943e-07

Fit with MPFITFUN in IDL

The python code is essentially the same:

from astropy.modeling import Fittable1DModel, Parameter, fitting
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

class WavModel(Fittable1DModel):
    '''Our model for fitting the wavelength dispersion solution
    based on an existing calibrated spectrum'''
    a = Parameter(name='a',default=0.)
    b = Parameter(name='b',default=1.)
    c = Parameter(name='c',default=0.)

    def __init__(self, a, b, c, spec=None, **kwargs):
        self.spec = spec
        self.pix = np.arange(spec.size, dtype=float)
        super(WavModel, self).__init__(a, b, c, **kwargs)

    def evaluate(self, x, a, b, c):
        wav = a + self.pix * b + self.pix * self.pix * c
        return interp1d(wav, self.spec, bounds_error=False, fill_value=0.)(x)

    fit_deriv = None

normalize = lambda arr: (arr - arr.min()) / (arr.max() - arr.min()) #scale arr into [0,1]

wavelength, flux = np.loadtxt('fiducial.txt').T
uncalibrated = np.loadtxt('to_fit.txt')

npix = uncalibrated.size
pix = np.arange(npix, dtype=float)
fitter = fitting.LevMarLSQFitter()
p_init = WavModel(1.95, 0.5/npix, 0., spec=normalize(uncalibrated))
coef = fitter(p_init, wavelength, normalize(flux)).parameters

print coef

#now we plot!
calib_wav = coef[0] + coef[1] * pix + coef[2] * pix * pix
plt.plot(calib_wav, normalize(uncalibrated), 'b')
plt.plot(wavelength, normalize(flux), 'r')
plt.savefig('python_test.png')

The results of THIS, however, are somewhat different:

[  1.98065438e+00   3.70817592e-04   1.06795564e-06]

Python fit to dispersion solution.

The IDL version is definitely the correct one, since the wavelength limit of the detector is somewhere around 2.5 microns...

I've pasted the data here and here. Any suggestions on how to get the python fit to look like the IDL fit? Thank you so much!

sirpercival
  • 236
  • 2
  • 9
  • You're using very different scales for both when it comes to plotting. Other than that, these seem to be pretty identical. – Slater Victoroff Nov 01 '15 at 15:46
  • The scales seem a lot different because the wavelength fit in the python case was wrong, so the x-axis range has to be larger to accommodate the incorrect calibration. I agree that they otherwise seem identical, so I'm wondering why (when the code is otherwise identical) the results are so different. – sirpercival Nov 01 '15 at 15:55
  • 1
    can you plot them with the same x-axis? It's very difficult to tell exactly what is different here other than the way you've chosen to plot the results. – Slater Victoroff Nov 01 '15 at 16:50
  • 1
    It would help to see some other diagnostics (e.g., chi^2) and, as Slater Tyranus said, the plots with the same X-axis. It's really not obvious that the two fits are different except in the dispersion is halved in the astropy fit. This question may be more appropriate for the astropy-users list (http://www.astropy.org/help.html) – keflavich Nov 02 '15 at 07:09
  • Do we know what algorithm mpfitfun is using? – Iguananaut Nov 03 '15 at 11:08
  • By the way, in this particular case the way you're using astropy.modeling it's little more than an OO wrapper around [scipy.optimize.leastsq](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html#scipy.optimize.leastsq). If you use that directly your code will look even more similar to the IDL version and it will be easier to identify where there might be a difference, and to rule out anything specifically having to do with Astropy. – Iguananaut Nov 03 '15 at 11:29
  • This isn't really an answer, but just to let you know, in IDL 8.5 we released the IDL-Python Bridge, which lets you call IDL code from Python (and vice versa). I'm wondering whether you could simply call your IDL code from your Python application, and not even worry about trying to get the Python code to work properly. Hope this helps! – Chris Torrence Nov 06 '15 at 17:53
  • @Iguananaut - `mpfitfun.pro` is a wrapper for `mpfit.pro`, both written by Craig Markwardt. They use the Levenberg-Marquardt algorithm. – honeste_vivere Jan 21 '16 at 18:52
  • @sirpercival - Is the `fitter` routine the python equivalent of `mpfit`? – honeste_vivere Jan 21 '16 at 18:58

0 Answers0