Source code for eispac.core.fit_spectra_astropy

__all__ = ['fit_spectra_astropy']

import copy
from datetime import datetime
import numpy as np
from astropy.modeling import fitting
import eispac.core.fitting_functions as fit_fns
from eispac.core.eiscube import EISCube
from eispac.core.eisfittemplate import EISFitTemplate
from eispac.core.eisfitresult import EISFitResult
from eispac.core.scale_guess import scale_guess
from eispac.extern.mpfit import mpfit
from eispac.core.generate_astropy_model import generate_astropy_model

[docs] def fit_spectra_astropy(inten, template, parinfo=None, wave=None, errs=None, min_points=10): """Fit one or more EIS line spectra using mpfit. Parameters ---------- inten : EISCube object or array_like One or more intensity profiles to be fit. The code will loop over the data according to its dimensionality. 3D data is assumed to be a full EIS raster (or a sub region), 2D data is assumed to be a single EIS slit, and 1D data is assumed to be a single profile. template : EISFitTemplate object or dict Either an EISFitTemplate or the just a 'template' dictionary. parinfo : dict, optional Dictionary of fit parameters formatted for use with mpfit. Required if the 'template' parameter is given just a dict and ignored otherwise. wave : array_like, optional Associated wavelength values for the spectra. Required if 'inten' is given an array and ignored otherwise. errs : array_like, optional Intensity error values for the spectra. Required if 'inten' is given an array and ignored otherwise. min_points : int, optional Minimum number of good quality data points (i.e. non-zero values & errs) to be used in each fit. Spectra with fewer data points will be skipped. Default is 10. Returns ------- fit_res :EISFitResult class instance An EISFitResult object containing the output fit paramaters. """ # If given an EISCube, extract data arrays and zero out masked values. # Otherwise, make local copies and just zero out negative values if isinstance(inten, EISCube): wave_cube = inten.wavelength.copy() errs_cube = inten.uncertainty.array.copy() inten_cube = inten.data.copy() metadata = inten.meta loc_masked = np.where(inten.mask == True) inten_cube[loc_masked] = 0 errs_cube[loc_masked] = 0 else: wave_cube = wave.copy() errs_cube = errs.copy() inten_cube = inten.copy() metadata = {'filename_data':'unknown', 'index':{}, 'pointing':{}} loc_neg = np.where(inten_cube < 0) inten_cube[loc_neg] = 0 errs_cube[loc_neg] = 0 # Check contents of template and parinfo if parinfo is None and isinstance(template, EISFitTemplate): parinfo_copy = copy.deepcopy(template.parinfo) template_copy = copy.deepcopy(template.template) elif parinfo is None: raise TypeError('Please input either a "parinfo" dictionary or a full' ' EISFitTemplate object.') else: parinfo_copy = copy.deepcopy(parinfo) template_copy = copy.deepcopy(template) # Check the dimensions of input data. # If the the arrays are not 3D, add shallow dimensions of size 1 num_dims = inten_cube.ndim dims_size = inten_cube.shape if num_dims == 1: n_pxls = 1 n_steps = 1 wave_cube = wave_cube[np.newaxis, np.newaxis, :] inten_cube = inten_cube[np.newaxis, np.newaxis, :] errs_cube = errs_cube[np.newaxis, np.newaxis, :] elif num_dims == 2: n_pxls = dims_size[0] n_steps = 1 wave_cube = wave_cube[:, np.newaxis, :] inten_cube = inten_cube[:, np.newaxis, :] errs_cube = errs_cube[:, np.newaxis, :] elif num_dims == 3: n_pxls = dims_size[0] n_steps = dims_size[1] wave_cube = wave_cube inten_cube = inten_cube errs_cube = errs_cube # Initalize output object which will contain the fit results fit_res = EISFitResult(wave_cube, template_copy, parinfo_copy, func_name='multigaussian') fit_res.meta = metadata fit_res.fit_module = 'astropy' fit_res.fit_method = 'LevMarLSQ' # get the indices of the function parameters n_gauss = fit_res.n_gauss n_poly = fit_res.n_poly # peaks = np.arange(n_gauss)*3 # cents = np.arange(n_gauss)*3+1 # wdths = np.arange(n_gauss)*3+2 # backs = n_gauss*3 # timer t1 = datetime.now() # oldguess = template_copy['fit'] fitter = fitting.LevMarLSQFitter() tmplt_model = generate_astropy_model(template.filename_temp) tmplt_init_param_vals = copy.deepcopy(tmplt_model.parameters) # Get indices of parameter that are NOT tied tied_mask = np.zeros(tmplt_model.parameters.size, dtype=bool) for p in range(tmplt_model.parameters.size): if tmplt_model.__dict__[tmplt_model.param_names[p]].tied != False: tied_mask[p] = True loc_untied = np.where(tied_mask == False) # loop over pixel positions and slit steps in the entire raster for ii in range(n_pxls): for jj in range(n_steps): # Extract a single profile from the raster wave_ij = wave_cube[ii,jj,::] inten_ij = inten_cube[ii,jj,::] errs_ij = errs_cube[ii,jj,::] # Only use good data loc_good = np.where(errs_ij > 0) if len(loc_good[0]) < min_points: fit_res.fit['status'][ii,jj] = -1 continue wave_ij = wave_ij[loc_good] inten_ij = inten_ij[loc_good] errs_ij = errs_ij[loc_good] # Do the fit using Astropy fit_model = fitter(tmplt_model, wave_ij, inten_ij, weights=1.0/errs_ij, maxiter=2000) # check convergence status if fitter.fit_info['ierr'] in [1, 2, 3, 4]: # update initial parameter values setattr(tmplt_model, 'parameters', fit_model.parameters) # assemble fit structure # TO-DO: figure out a more compact way to do this that is still # easy to read and debug fit_dof = wave_ij.size - loc_untied[0].size chi_sq = sum(((inten_ij - fit_model(wave_ij))/errs_ij)**2) if fitter.fit_info['param_cov'] is not None: untied_perror = np.sqrt(np.diag(fitter.fit_info['param_cov'])) else: untied_perror = np.zeros(loc_untied[0].size) fit_res.fit['status'][ii,jj] = fitter.fit_info['ierr'] fit_res.fit['chi2'][ii,jj] = chi_sq/fit_dof fit_res.fit['wavelength'][ii,jj,::] = wave_cube[ii,jj,::] # fit_res.fit['peak'][ii,jj,::] = out.params[peaks] # fit_res.fit['err_peak'][ii,jj,::] = out.perror[peaks] # fit_res.fit['centroid'][ii,jj,::] = out.params[cents] # fit_res.fit['err_centroid'][ii,jj,::] = out.perror[cents] # fit_res.fit['width'][ii,jj,::] = out.params[wdths] # fit_res.fit['err_width'][ii,jj,::] = out.perror[wdths] # fit_res.fit['background'][ii,jj,::] = out.params[backs] # fit_res.fit['err_background'][ii,jj,::] = out.perror[backs] fit_res.fit['params'][ii,jj,::] = fit_model.parameters fit_res.fit['perror'][ii,jj,loc_untied] = untied_perror # fit_res.fit['int'][ii,jj,::] = int # fit_res.fit['err_int'][ii,jj,::] = err_int else: print(' ! fit did not converge!') fit_res.fit['status'][ii,jj] = fitter.fit_info['ierr'] # print status print(' + fit completed!') # timer t2 = datetime.now() print(' + fit runtime : {}'.format(t2-t1)) return fit_res