Python mpfit.mpfit函数代码示例

def do_mpfit(x,y,covarin,guess,functname=thepolynomial):
    # check if covariance or error bars were given
    covar=covarin
    if np.size(np.shape(covarin)) == 1:
        err=covarin
        print('err')
        #Prepare arguments for mpfit
        fa={'x':double(x),'y':double(y),'err':double(err),'functname':functname}
        costfct=fdeviates
    else:
        print('covar')
        #Fist do a SVD decomposition
        u,s,v=np.linalg.svd(covarin)
        #Prepare arguments for mpfit
        fa={'x':double(x),'y':double(y),'svdvals':double(s),'v':double(v),'functname':functname}
        costfct=chi2svd

    #Run MPFIT
    mpf = mpfit.mpfit(costfct, guess, functkw=fa, autoderivative=1)
    print('Status of the Fit',mpf.status)
    print('Chi2=',mpf.fnorm)
    print('ndf=',mpf.dof)
    print('Fitted params:',mpf.params)

    return(mpf,mpf.params,mpf.perror,mpf.covar)

def test_quadfit():

    x = N.array([-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
          -2.8394297E-01,1.3416969E+00,1.3757038E+00,
          -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
          8.2065094E-01])
    y = N.array([2.3095947E+01,2.6449392E+01,1.0204468E+01,
          5.40507,1.5787588E+01,1.6520903E+01,
          1.5971818E+01,4.7668524E+00,4.9337711E+00,
          8.7348375E+00])
    ey=0.2*N.ones(y.shape,dtype='float64')
    p0=N.array([1.0,1.0,1.0],dtype='float64')  #initial conditions
    pactual=N.array([4.7, 0.0, 6.2]) #actual values used to make data
    parbase={'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}
    parinfo=[]
    for i in range(len(pactual)):
        parinfo.append(copy.deepcopy(parbase))
        parinfo[i]['value'] = p0[i]
    fa = {'x': x, 'y': y, 'err': ey}
    m = mpfit(myfunctquad, p0, parinfo = parinfo, functkw = fa)
    if (m.status <= 0):
        print 'status = ', m.status
    print m.perror
    assert N.allclose(m.params, N.array([ 4.703829, 0.062586,  6.163087], dtype='float64'))
    assert N.allclose(m.perror, N.array([ 0.097512, 0.054802,  0.054433], dtype='float64'))
    chisq = (myfunctquad(m.params, x=x, y=y, err=ey)[1]**2).sum()
    assert N.allclose(N.array([chisq], dtype='float64'), N.array([5.679323], dtype='float64'))
    assert m.dof == 7
    return

def test_gaussfit_fixed():

    x = N.array([-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
        -2.8394297E-01,1.3416969E+00,1.3757038E+00,
        -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
        8.2065094E-01])
    y = N.array([-4.4494256E-02,8.7324673E-01,7.4443483E-01,
        4.7631559E+00,1.7187297E-01,1.1639182E-01,
        1.5646480E+00,5.2322268E+00,4.2543168E+00,
        6.2792623E-01])
    ey=0.5*N.ones(y.shape,dtype='float64')
    p0=N.array([0.0, 1.0, 0.0, 0.1],dtype='float64')  #initial conditions          
    pactual=N.array([0.0, 4.70, 0.0, 0.5]) #actual values used to make data
    parbase={'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}
    parinfo=[]
    for i in range(len(pactual)):
        parinfo.append(copy.deepcopy(parbase))
        parinfo[i]['value'] = p0[i]
    
    parinfo[0]['fixed'] = 1
    parinfo[2]['fixed'] = 1

    fa = {'x': x, 'y': y, 'err': ey}
    m = mpfit(myfunctgauss, p0, parinfo = parinfo, functkw = fa)
    if (m.status <= 0):
        print 'status = ', m.status
    assert N.allclose(m.params, N.array([ 0., 5.059244,  0.,  0.479746 ], dtype='float64'))
    assert N.allclose(m.perror, N.array([ 0., 0.329307,  0.,  0.053804 ], dtype='float64'))
    chisq = (myfunctgauss(m.params, x=x, y=y, err=ey)[1]**2).sum()
    assert N.allclose(N.array([chisq], dtype='float64'), N.array([15.516134], dtype='float64'))
    assert m.dof == 8
    return

def zandtfit(Inclination, Planetaryradius_SI, Keplerlc, Teff, Logg, Metalicity, sigma, Orbitalperiod_SI, KeplerID, Planetnumberer):
	Planetaryradius_Earth = (Planetaryradius_SI) /  6.3675E6 #convert the planetary radius to earth radii
	p0 = [math.cos(1.570795), Planetaryradius_Earth, Teff, Logg, np.log10(Metalicity)] #inital guess (from kepler mast and headers)
	fa = {'Keplerlc':Keplerlc, 'sigma':sigma, 'Orbitalperiod_SI':Orbitalperiod_SI, 'KeplerID':KeplerID, 'Planetnumberer':Planetnumberer} #additional variables to be past to the function
	parinfo = [{'value':p0[0], 'fixed':0, 'limited':[1,1], 'limits':[-1,1], 'step':0.01}, {'value':p0[1], 'fixed':0, 'limited':[1,0], 'limits':[0.01,100], 'step':0},\
		{'value':p0[2], 'fixed':0, 'limited':[1,1], 'limits':[0,40000], 'step':0}, {'value':p0[3], 'fixed':0, 'limited':[1,1], 'limits':[0,5.0], 'step':0},\
		{'value':p0[4], 'fixed':0, 'limited':[1,1], 'limits':[-5,1], 'step':0} ]
	m = mpfit(minimisefunction, p0, functkw = fa, quiet = 0, maxiter = 35, parinfo = parinfo, ftol=0.0001) #run the mpfit for the best fit
	#print m #testing
	bestfitarray = m.params #extract the results 
	errors = m.perror #these need to be properly stored
	covar = m.covar #these need to be properly stored
	if errors == None: #i.e. if mpfits is a pain and decides to not return the errors or covariance properly
		errors = np.zeros(5) #create an empty errors array
	if covar == None: #if no covariance array is produced, create an array of zeros so that there are no problems later
		covar = np.zeros((5,5))
	#print bestfitarray
	Inclination_fit = math.acos(bestfitarray[0]) #Extract the fit inclination
	Planetaryradius_SI_fit =   bestfitarray[1] * 6.3675E6 #Extract the fit planetary radius and convert it to SI units
	Teff_fit = bestfitarray[2] #Extract the fit Teff
	Logg_fit = bestfitarray[3] #Extract the fit logg of the surface gravity
	Metalicity_fit = (10**bestfitarray[4]) #Extract the metalicity and remove the logg
	Mstar_SI, Rstar_SI, Stellardensity_SI = RMrhocalc(Teff_fit, Logg_fit, Metalicity_fit) #Calculate the stellar mass, radius and density
	semimajoraxis_SI = semimajoraxiscalc(Orbitalperiod_SI, Mstar_SI) #calculate the semi-major axis
	T_dur_SI = Transitdurationcalc(Rstar_SI, Planetaryradius_SI_fit, semimajoraxis_SI, Inclination_fit, Orbitalperiod_SI) #calculate the transit duration
	return Inclination_fit, Planetaryradius_SI_fit, Teff_fit, Logg_fit, Metalicity_fit, Mstar_SI, Rstar_SI, Stellardensity_SI, semimajoraxis_SI, T_dur_SI, errors, covar

    def _run_mpfit(self):
        """Run an MPFIT regression"""
        # Make up the PARINFO list

        parinfo = []
        for i in range(len(self._guess)):
            pdict = {
                "value": self._guess[i].copy(),
                "fixed": 0,
                "limited": self._ilimited[i],
                "limits": self._ilimits[i],
            }
            parinfo.append(pdict)

        quiet = 1
        if self.debug:
            quiet = 0

        m = mpfit.mpfit(self._residualsMPFIT, self._guess, parinfo=parinfo, quiet=quiet, debug=self.debug)

        self._result = m.params
        self._niter = m.niter

        if m.perror is not None:
            self._stdev = m.perror
        else:
            self._stdev = zeros(len(self._guess))

        if m.covar is not None:
            self._covar = m.covar
        else:
            self._covar = zeros((len(self._guess), len(self._guess)))

        self._mpfit = m

def newfit2( infile, outfile, xcol, ycol, p0, \
    constrain=[ {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0} ] ) :
  [x,y] = getdata( infile, xcol, ycol )
  err = numpy.ones( len(x), dtype=float )
  parinfo = [ {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0} ]
  fa = {'x':x, 'y':y, 'err':err}
  m = mpfit.mpfit(myfunc, p0, functkw=fa, parinfo=constrain)

  fout = open( outfile, "w" )
  fout.write("# %s\n" %  outfile )
  fout.write("# fit to spectrum %s, xcol=%d, ycol=%d\n" % ( infile, xcol, ycol ) )
  fout.write("# final fit:\n")
  fout.write("#   yavg= %0.5f   %d \n" %  (m.params[0], constrain[0]['fixed'])  )
  fout.write("#   slope = %0.5f   %d\n" % (m.params[1], constrain[1]['fixed'])  )
  fout.write("#   amp = %0.5f   %d\n" %  (m.params[2], constrain[2]['fixed'])  )
  fout.write("#   v0 = %0.3f   %d\n" %  (m.params[3], constrain[3]['fixed'])  )
  fout.write("#   FWHM = %0.3f   %d\n" % (m.params[4], constrain[4]['fixed'])  )
  fout.write("#\n")

  xavg = 0.5 * (x[0] + x[len(x)-1])
  fit = eval('m.params[0] + m.params[1]*(x-xavg) + \
     m.params[2]*numpy.exp(-2.772 * pow( (x-m.params[3])/m.params[4], 2.))')
     # fit to raw data, including slope
  yminusSlope = eval('y - m.params[1]*(x-xavg)')
     # raw data with linear slope removed
  gfit = eval('m.params[0] + m.params[2]*numpy.exp(-2.772 * pow( (x-m.params[3])/m.params[4], 2.))')
     # gaussian fit to data with slope removed
  yminusSlopeNorm = eval('yminusSlope/m.params[0]')
     # normalized absorption 
  gfitNorm = eval('gfit/m.params[0]')
  for i in range(0,len(x)) :
    fout.write(" %10.3f  %9.4f  %9.4f  %9.4f  %9.4f %9.5f %9.5f\n" % \
	  (x[i],y[i],fit[i],yminusSlope[i],gfit[i],yminusSlopeNorm[i],gfitNorm[i] ) )
  fout.close()

def onedgaborfit(xax, data, err=None,
        params=[0,1,0,1,1,0],fixed=[False,False,False,False,False,False],
        limitedmin=[False,False,False,True,True,True],
        limitedmax=[False,False,False,False,False,True], minpars=[0,0,0,0,0,0],
        maxpars=[0,0,0,0,0,360], quiet=True, shh=True,
        veryverbose=False):
    """
    Inputs:
       xax - x axis
       data - y axis
       err - error corresponding to data

       params - Fit parameters: Height of background, Amplitude, Shift, Width, Wavelength, Phase
       fixed - Is parameter fixed?
       limitedmin/minpars - set lower limits on each parameter (default: width>0)
       limitedmax/maxpars - set upper limits on each parameter
       quiet - should MPFIT output each iteration?
       shh - output final parameters?

    Returns:
       Fit parameters
       Model
       Fit errors
       chi2
    """

    def mpfitfun(x,y,err):
        if err is None:
            def f(p,fjac=None): return [0,(y-onedgabor(x,*p))]
        else:
            def f(p,fjac=None): return [0,(y-onedgabor(x,*p))/err]
        return f

    if xax == None:
        xax = np.arange(len(data))

    parinfo = [ {'n':0,'value':params[0],'limits':[minpars[0],maxpars[0]],'limited':[limitedmin[0],limitedmax[0]],'fixed':fixed[0],'parname':"HEIGHT",'error':0} ,
                {'n':1,'value':params[1],'limits':[minpars[1],maxpars[1]],'limited':[limitedmin[1],limitedmax[1]],'fixed':fixed[1],'parname':"AMPLITUDE",'error':0},
                {'n':2,'value':params[2],'limits':[minpars[2],maxpars[2]],'limited':[limitedmin[2],limitedmax[2]],'fixed':fixed[2],'parname':"SHIFT",'error':0},
                {'n':3,'value':params[3],'limits':[minpars[3],maxpars[3]],'limited':[limitedmin[3],limitedmax[3]],'fixed':fixed[3],'parname':"WIDTH",'error':0},
                {'n':4,'value':params[4],'limits':[minpars[4],maxpars[4]],'limited':[limitedmin[4],limitedmax[4]],'fixed':fixed[4],'parname':"WAVELENGTH",'error':0},
                {'n':5,'value':params[5],'limits':[minpars[5],maxpars[5]],'limited':[limitedmin[5],limitedmax[5]],'fixed':fixed[5],'parname':"PHASE",'error':0}]

    mp = mpfit(mpfitfun(xax,data,err),parinfo=parinfo,quiet=quiet)
    mpp = mp.params
    mpperr = mp.perror
    chi2 = mp.fnorm

    if mp.status == 0:
        raise Exception(mp.errmsg)

    if (not shh) or veryverbose:
        print "Fit status: ",mp.status
        for i,p in enumerate(mpp):
            parinfo[i]['value'] = p
            print parinfo[i]['parname'],p," +/- ",mpperr[i]
        print "Chi2: ",mp.fnorm," Reduced Chi2: ",mp.fnorm/len(data)," DOF:",len(data)-len(mpp)

    return mpp,onedgabor(xax,*mpp),mpperr,chi2

def ab_fit(src,xd,td,lguess,xmin_id,xmax_id,evmin1,evmax1,evmin2,evmax2):
	npar    =  len(lguess)
	guessp  = np.array(lguess, dtype='float64')
	plimd   = [[False,False]]*npar
	plims   = [[0.,0.]]*npar
	parbase = {'value': 0., 'fixed': 0, 'parname': '', 'limited': [0, 0], 'limits': [0., 0.]}
	pname   = ['bg']*npar
	pfix    = [False]*npar

	count   = 1
	for i in range(npar):
		if(i==0):
			pname[i] = 'Tcont'

		if( ((i-1)%3 == 0) and (i>0) ):
			pname[i] = 'tau'+str(count)

		if( ((i-2)%3 == 0) and (i>0) ):
			pname[i] = 'v0'+str(count)

		if( ((i-3)%3 == 0) and (i>0) ):
			pname[i] = 'wid'+str(count)
			count    = count + 1

	parinfo = []
	for i in range(len(guessp)):
		parinfo.append(copy.deepcopy(parbase))

	for i in range(len(guessp)):
		parinfo[i]['value']   = guessp[i]
		parinfo[i]['fixed']   = pfix[i]
		parinfo[i]['parname'] = pname[i]
		parinfo[i]['limited'] = plimd[i]

	##  data and 1-sigma uncertainty ##
	x    = xd[xmin_id:xmax_id]
	y    = td[xmin_id:xmax_id]
	sg,\
	erry = get_tb_sigma(xd, td, evmin1, evmax1, evmin2, evmax2)
	erry = erry[xmin_id:xmax_id]

	x    = x.astype(np.float64)
	y    = y.astype(np.float64)
	erry = erry.astype(np.float64)

	fa   = {'x':x, 'y':y, 'err':erry}
	mp   = mpfit(myabfunc, guessp, parinfo=parinfo, functkw=fa, quiet=True)

	## Fit values of Tau ##
	abp   = mp.params
	abper = mp.perror
	fit   = 0.
	for i in range(1, len(abp), 3):
		fit = fit + abp[i]*np.exp(- ( (x-abp[i+1])/(0.6005612*abp[i+2]))**2)  

	fit = np.exp(-fit)
	# fit = abp[1]*fit

	return x,fit,y/abp[0],abp,abper,npar,parbase,pname,parinfo

def delPuncorr_bestfit(k, delP_over_P): 
    '''
    '''
    fa = {'x': k, 'y': delP_over_P}
    
    param_guess = [-4.0]
    bestfit = mpfit.mpfit(nongauss_poly_mpfit, param_guess, functkw=fa, quiet=True)
    return bestfit.params

def mpfitexpr(func, x, y, err, start_params, check=True, full_output=False,
              imports=None, **kw):
    """Fit the used defined expression to the data
    Input:
    - func: string with the function definition
    - x: x vector
    - y: y vector
    - err: vector with the errors of y
    - start_params: the starting parameters for the fit
    Output:
    - The tuple (params, yfit) with best-fit params and the values of func evaluated at x
    Keywords:
    - check: boolean parameter. If true(default) the function will be checked for sanity
    - full_output: boolean parameter. If True(default is False) then instead of best-fit parameters the mpfit object is returned
    - imports: list of strings, of optional modules to be imported, required to evaluate the function
    Example:
    params,yfit=mpfitexpr('p[0]+p[2]*(x-p[1])',x,y,err,[0,10,1])

    If you need to use numpy and scipy functions in your function, then
            you must to use the full names of these functions, e.g.:
            numpy.sin, numpy.cos etc.

    This function is motivated by mpfitexpr() from wonderful MPFIT IDL package
            written by Craig Markwardt

    """

    hash = {}
    hash['numpy'] = numpy
    hash['scipy'] = scipy

    if imports is not None:
        for i in imports:
            #exec '%s=__import__("%s")'%(a,b) in globals(),locals()
            hash[i] = __import__(i)

    def myfunc(p, fjac=None, x=None, y=None, err=None):
        return [0, eval('(y-(%s))/err' % func, hash, locals())]

    myre = "[^a-zA-Z]p\[(\d+)\]"
    r = re.compile(myre)
    maxp = -1
    for m in re.finditer(r, func):
        curp = int(m.group(1))
        maxp = curp if curp > maxp else maxp
    if check:
        if maxp == -1:
            raise Exception("wrong function format")
        if maxp + 1 != len(start_params):
            raise Exception(
                "the length of the start_params != the length of the parameter verctor of the function")
    fa = {'x': x, 'y': y, 'err': err}
    res = mpfit.mpfit(myfunc, start_params, functkw=fa, **kw)
    yfit = eval(func, globals(), {'x': x, 'p': res.params})
    if full_output:
        return (res, yfit)
    else:
        return (res.params, yfit)

def main():

    # Generate a noisy polynomial
    #     [off,  x1,  x2,  x3,  y1,  y2,  y3]
    pIn = [2.0, 1.5, 0.1, 0.3, 1.0, 2.0, 0.05]
    pIn = [1.0, 0.2, 0.0, 0.0, 0.1, 0.0, 0.0]
    shape = (200, 200)
    X, Y, Z, xyData = genpolydata(pIn, shape, 300, 10.2)

    print X.shape
    print Y.shape
    print Z.shape
    print xyData.shape
    # raw_input()

    # Define an function to evaluate the residual
    def errFn(p, fjac=None):
        status = 0
        # poly_surface' returns the 'rfunc' function and the X,Y data is
        # inserted via argument unpacking.
        return status, poly_surface(p)(*[Y, X]) - Z
    
    # Fit the data starting from an initial guess
    mp = mpfit(errFn, parinfo=inParms, quiet=False)
    print 
    for i in range(len(inParms)):
        print "%s = %f +/- %f" % (inParms[i]['parname'],
                                  mp.params[i],
                                  mp.perror[i])
    p1 = mp.params

    #-------------------------------------------------------------------------#
  
    # Plot the original, fit & residual
    fig = pl.figure(figsize=(18,4.3))
    
    ax1 = fig.add_subplot(1,3,1)
    cax1 = ax1.imshow(xyData, origin='lower',cmap=mpl.cm.jet)
    cbar1=fig.colorbar(cax1, pad=0.0)
    ax1.scatter(X, Y, c=Z, s=40, cmap=mpl.cm.jet)
    ax1.set_title("Sampled Data")
    ax1.set_xlim(0, shape[-1]-1)
    ax1.set_ylim(0, shape[-2]-1)
    ax1.set_aspect('equal')
    
    ax2 = fig.add_subplot(1,3,2)
    xyDataFit = poly_surface(p1, shape)
    cax2 = ax2.imshow(xyDataFit, origin='lower', cmap=mpl.cm.jet)
    cbar2=fig.colorbar(cax2, pad=0.0)
    ax2.set_title("Model Fit")
    
    ax3 = fig.add_subplot(1,3,3)
    xyDataRes = xyData - xyDataFit
    cax3 = ax3.imshow(xyDataRes, origin='lower', cmap=mpl.cm.jet)
    cbar2=fig.colorbar(cax3, pad=0.0)
    ax3.set_title("Residual")
    
    pl.show()

def disk_bulge_AIM(Iobs,err,psf,guessP,quiet=0,ftol=1.e-10,return_image=False):

	#Dictionary to pass parameters to function `handleFunc`
	fa = {'y':Iobs,'x':psf,'err':err}

	#Possible constraints for each parameter
	parinfo= [{'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]} for i in range(len(guessP))]
	#Pass initial values to parinfo
	for i in range(len(guessP)):
		parinfo[i]['value'] = guessP[i]
	
	L = guessP[0]

	#Set parameter limits and whether fixed
	parinfo[0]['fixed'] = 1				#xmax
	parinfo[1]['fixed'] = 1				#N
	
	parinfo[2]['limited'] = [1,1]		#disk flux
	parinfo[2]['limits'] = [0.1,1.e10] #
	parinfo[3]['limited'] = [1,1]		#bulge flux
	parinfo[3]['limits'] = [0.1,1.e10] #
	parinfo[4]['limited'] = [1,1]		#q (q=10 corresponds to e = 0.82)
	parinfo[4]['limits'] = [1.0,20.0]
	parinfo[5]['limited'] = [1,1]		#phi (there is 2pi degeneracy, but keep to avoid boundaries)
	parinfo[5]['limits'] = [0.0,6.283185]
	parinfo[6]['limited'] = [1,1]		#disk re
	parinfo[6]['limits'] = [0.05,3.0*L] #
	parinfo[7]['limited'] = [1,1]		#bulge re
	parinfo[7]['limits'] = [0.05,3.0*L] #
	parinfo[8]['limited'] = [1,1]
	parinfo[8]['limits'] = [0.5,4.0]  #n_b
	

	parinfo[9]['fixed'] = 1				#kappa
	parinfo[10]['fixed'] = 1				#Shear 1
	parinfo[11]['fixed'] = 1				#Shear 2
	parinfo[12]['limited'] = [1,1]		#f1
	parinfo[12]['limits'] = [-1.0,1.0]
	parinfo[13]['limited'] = [1,1]		#f2
	parinfo[13]['limits'] = [-1.0,1.0]	
	parinfo[14]['limited'] = [1,1]		#g1
	parinfo[14]['limits'] = [-1.0,1.0]	
	parinfo[15]['limited'] = [1,1]		#g2
	parinfo[15]['limits'] = [-1.0,1.0]
	parinfo[16]['limited'] = [1,1]		#cx
	parinfo[16]['limits'] = [-L,L]
	parinfo[17]['limited'] = [1,1]		#cy
	parinfo[17]['limits'] = [-L,L]


	m = mpfit.mpfit(disk_bulge_func,guessP,parinfo=parinfo,functkw = fa,quiet=quiet,ftol=ftol)

	if return_image == True:
		fit_image = disk_bulge_model(m.params,psf)
		return m,fit_image
	else:
		return m

    def __init__(self,f, xdata, ydata, p0=None, sigma=None, fixed=None,limits=None, maxiter=200,quiet=1):
        args, varargs, varkw, defaults = inspect.getargspec(f)
       
        self.quiet = quiet
        if len(args) < 2:
            msg = "Unable to determine number of fit parameters."
            raise ValueError(msg)
       
        def general_function(params, fjac=None, xdata=None, ydata=None):
            return [0, f(xdata, *params) - ydata]

        def weighted_general_function(params,fjac=None, xdata=None, ydata=None, weights=None):
            return [0, weights * (f(xdata, *params) - ydata)]


        parinfo = [{'value':1., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}  for i in range(len(args)-1)] 
        

        if p0 != None:
            for x,y in zip(parinfo,p0):
                x['value']=y
        if fixed != None:
            for x,y in zip(parinfo,fixed):
                x['fixed']=y


        if limits != None:
            for x,y in zip(parinfo,limits):
                x['limited']=y[0]
                x['limits']=y[1]
            



    #args = (xdata, ydata, f)
        if sigma is None:
            func = general_function
           # print 'here'
          #  my_functkw = {'xdata':xdata,'ydata':ydata,'function':f}
            my_functkw = {'xdata':xdata,'ydata':ydata}
        else:
            func = weighted_general_function
            weights= 1.0/asarray(sigma)
           # my_functkw = {'xdata':xdata,'ydata':ydata,'function':f,'weights':weights}
            my_functkw = {'xdata':xdata,'ydata':ydata,'weights':weights}


        
       
        result = mpfit(func, functkw=my_functkw, parinfo = parinfo,quiet=self.quiet,maxiter=maxiter)
      
        self.params = result.params
        self.errors = result.perror
        self.dof = result.dof
        self.chi2 = result.fnorm
        self.covar = result.covar

def fit_pixel(inputs):
	wl,q_array,u_array,sigma_q,sigma_u=inputs
	pos=[1,1,10]
	funcargs= {'x':wl, 'y':[q_array,u_array],'err':[sigma_q,sigma_u]}
	fit=mpfit.mpfit(alpha_function,pos,functkw=funcargs,autoderivative=1,quiet=1)
	parms=fit.params
	parms[:2]*=1e-7
	dparms=fit.perror
	dparms[:2]*=1e-7
	return np.array([parms,dparms,fit.fnorm])

def dlosenv_peakfit_sigma_env_fit_test(cat_corr, n_NN=3, fit='gauss', **kwargs): 
    """ Linear fit of the best-fit sigma as a function of environment. 
    Bestfit sigma values are computed for the dLOS distribution in bins of 
    galaxy environment. Linear fit is done using MPFit.

    """

    env_low, env_high, fpeaks, fpeak_errs, sigmas, sigma_errs, amps, nbins =  \
            dlos_envbin_peakfit(
                    cat_corr, 
                    n_NN = n_NN,
                    fit = fit, 
                    **kwargs
                    )
        
    env_mid = np.array( 0.5 * (env_low + env_high) )
    p0 = [ -0.03, 4.0 ] # guess
    fa = {'x': env_mid, 'y': sigmas, 'err': sigma_errs}
    
    fit_param = mpfit.mpfit(mpfit_linear, p0, functkw=fa, quiet=1)
        
    best_slope = fit_param.params[0]
    best_yint = fit_param.params[1]

    print 'Entire range'
    print best_slope, best_yint
    
    for range_test in [20.0, 30.0, 40.0, 50.0]: 

        test_range = np.where(
                env_mid < range_test 
                )
        fa = {'x': env_mid[test_range], 'y': sigmas[test_range], 'err': sigma_errs[test_range]}
        
        fit_param = mpfit.mpfit(mpfit_linear, p0, functkw=fa, quiet=1)
            
        best_slope = fit_param.params[0]
        best_yint = fit_param.params[1]
        
        print 'dNN < ', str(range_test)
        print best_slope, best_yint

    return best_slope, best_yint