script for linear least squared

test_fit

from astropy.io import fits
from scipy import optimize
import numpy as np
from extract_spectral.trace import GaussianPolynomialTrace, ModelFrame2D, LogLikelihood, LinearLstSqExtraction

from DIRECT import solve
from astropy.modeling import models
#download from 
data = fits.getdata('IRCQ00467448.fits')
#data = fits.getdata('test_data/fors_many_stars.fits').T
#cutout = data
cutout = data[:,150:250]

#cutout[500:, 50:] = np.mean(cutout)

x, y = np.mgrid[:cutout.shape[0], :cutout.shape[1]]


bounds = ([(0.1, 5)] + [(-0.5, 0.5)] * 4)
bounds = np.array(bounds)
l = bounds[:,0]
u = bounds[:,1]

model = ModelFrame2D(cutout) | (GaussianPolynomialTrace(5, sigma=5,
domain=(0, cutout.shape[0]), window=(0, cutout.shape[1]))) | LinearLstSqExtraction(cutout)


llhood = model | LogLikelihood(cutout)


def fit_func(x, data=None):

    llhood_val = -2 * llhood.evaluate(x[0], *x[1:])
    #print x, llhood_val
    return llhood_val
    
# to run I use:

optimize.differential_evolution(fit_func, bounds, disp=True)