# Fitting the Data

Fitting of the spectra involves selecting a spectral line of interest (e.g. 195.119 Å) from one of the spectral windows of in the data, choosing a function (or combination of functions) to fit, and determining an initial guess for each parameter. The next ingredient for a fit is the selection of an optimization method. By default, EISPAC uses a Python implementation of the well-known IDL method MPFIT, which solves the non-linear least squares problem using the Levenberg- Marquardt algorithm. The Python module, mpfit.py, can be found on GitHub  and is included in EISPAC. Future versions of the code will include full support for other fitting packages such as the newer astropy.modeling framework.

## Template Files

In order to make fitting quick and easy, we’ve created a set of fit templates for different spectral lines. Each template consists of one or more component Gaussian functions and a constant background term. Template files are named according to the following pattern: {primary spectral line}.{number of Gaussians}c.template.h5. 119 templates are provided with EISPAC, covering all of the commonly observed spectral line combinations.

The easiest way to browse and copy fit template files is to use the eis_browse_templates GUI tool included with EISPAC. This tool can be used to display a list of all templates available for each spectral window in an EIS observation and copy the template files to your work directory (see the eis_browse_templates section for a few tips). A similar result can be obtained with the match_templates function, which takes either an EISCube or the path to a header file and outputs a list (or list of lists) of all fit templates available for a given spectral window (or data file).

>>> import eispac
>>> data_filename = 'eis_20190404_131513.data.h5'
>>> template_list = eispac.match_templates(data_cube)
>>> for file in template_list:
>>>     print(file.name)
fe_12_195_119.1c.template.h5
fe_12_195_119.2c.template.h5
fe_12_195_179.2c.template.h5


In the example above, notice how both one and two component fits are available for the same line. Furthermore, as with many multi-component templates, the third template file is actually a duplicate of the second template, just with the name of the secondary spectral line in the filename. This was designed to make it easier to see exactly which spectral lines can be fit. However, it also means you cannot simply loop over all template files available for a given spectral window, as you would be unnecessarily multiplying the amount of fitting performed (and resultant output files).

We recommend copying the template files you wish to use to your project directory so you have a more self-contained record of your fitting process. Assuming you are running Python from your project directory, the following code snippet illustrates one method you could use.

>>> import os
>>> import shutil
>>> current_dir = os.getcwd()
>>> shutil.copy(template_list, current_dir)


Attention

Take care to only copy template files. Moving a template from your EISPAC installation will lead to incomplete results in future calls to match_templates().

An h5dump  on one of the template files shows that it contains a /template group for the initial guess on the fit parameters and a /parinfo group containing constraints on the parameters for use by mpfit.

h5dump -n fe_12_195_119.2c.template.h5
HDF5 "fe_12_195_119.2c.template.h5" {
FILE_CONTENTS {
group      /
group      /parinfo
dataset    /parinfo/fixed
dataset    /parinfo/limited
dataset    /parinfo/limits
dataset    /parinfo/tied
dataset    /parinfo/value
group      /template
dataset    /template/component
dataset    /template/data_e
dataset    /template/data_x
dataset    /template/data_y
dataset    /template/fit
dataset    /template/fit_back
dataset    /template/fit_gauss
dataset    /template/line_ids
dataset    /template/n_gauss
dataset    /template/n_poly
dataset    /template/order
dataset    /template/wmax
dataset    /template/wmin
}


The read_template function can be used to read a template file and examine the contents.

>>> import eispac
>>> tmplt_filename = 'fe_12_195_119.2c.template.h5'


Use the command print(TEMPLATE) to view the initial parameter values and constraints in a nice format.

>>> print(tmplt)
--- FIT TEMPLATE PARAMETER CONSTRAINTS ---
*         Value      Fixed        Limited            Limits               Tied
p     57514.6647     0          1    0       0.0000       0.0000
p       195.1179     0          1    1     195.0778     195.1581
p         0.0289     0          1    1       0.0191       0.0510
p      8013.4013     0          1    0       0.0000       0.0000
p       195.1779     0          1    1     195.1378     195.2181          p+0.06
p         0.0289     0          1    1       0.0191       0.0510          p
p       664.3349     0          0    0       0.0000       0.0000


The structure of parinfo is specific to MPFIT and should be familiar to anyone who has used the original IDL version; please see the section on Constraining Parameter Values for more details. The templates provided with EISPAC consist of one or more Gaussian functions (with parameters in the order of peak, centroid, & width) followed by one or more background polynomial terms (usually just a single, constant value). The values .template['n_gauss'] and .template['n_poly'] indicate, respectively, the number of Gaussian functions and background polynomial terms in a given template.

Note

The multigaussian function is composed of generalized Gaussian functions of the form $$f(x) = A exp(-(x-b)^2/2c^2)$$, where A is the amplitude (peak value), b is the position of the center of the peak (centroid), and c is the standard deviation (width). This is consistent with the fit parameters used for EIS data in the IDL SolarSoftWare (SSW) analysis suite.

## Fitting Spectra

Once you’ve read in a template file, you can use the central wavelength to find the desired spectral window in the data using read_cube.

>>> data_filename = 'eis_20190404_131513.data.h5'


As mentioned in the previous chapter, read_cube automatically applies all of the pointing and wavelength corrections, bad data masking, and error estimations needed for scientific analysis. By default, the code also converts the data from photon counts to intensity units of erg cm$$^{-2}$$ s$$^{-1}$$ sr$$^{-1}$$ using the appropriate pre-flight calibration curve. This conversion can be disabled by setting the keyword apply_radcal=False, should you prefer to run your fits in count space.

On to the fitting! Now that you have a template and the data elements, you can perform a fit of the entire data cube by calling the top-level fitting routine, fit_spectra. The easiest way to use fit_spectra is to just give it both an EISCube and EISFitTemplate object (or filepaths to the data and template HDF5 files). You may slice your EISCube however you wish before fitting and the code will loop over the data appropriately (this includes fitting a single spectra or slit observation). Additionally, fit_spectra takes advantage of the multiprocessing package in the Python standard library to automatically parallelize the fitting process and minimize the run time. You may control the number of processing cores used for the fitting with ncpu keyword, or set it equal to “max” or None to use the maximum number of cores available. Please see the full doc string for fit_spectra for additional options and parameters.

Attention

Due to the specifics of how the multiprocessing library works, any statements that call fit_spectra() using ncpu > 1 MUST be wrapped in a "if __name__ == __main__:" statement in the top-level script or program. If such a “name guard” statement is not detected, fit_spectra() will fall back to using a single process. Unfortunately, this means you can not directly use parallel fitting from an interactive Python shell, you must first write a program that you save and run.

Here is a minimal example program that just loads and fits the data.

import eispac

if __name__ == '__main__':
# input data and template files
data_filepath = './eis_20190404_131513.data.h5'
template_filepath = './fe_12_195_119.2c.template.h5'

# Read spectral window into an EISCube

# Fit the data, then save it to disk and test loading it back in
fit_res = eispac.fit_spectra(data_cube, tmplt, ncpu='max')
save_filepaths = eispac.save_fit(fit_res, save_dir='cwd')
FITS_file = eispac.export_fits(fit_res, save_dir='cwd')


Note

The command line script eis_fit_files can be used to quickly fit a directory of files using one or more templates in another directory.

## EISFitResult Objects

fit_spectra outputs an EISFitResult object, which may be saved to an HDF5 file and read back in later using the save_fit and read_fit functions (as shown in the example above). The output fit parameters are stored in a dictionary of arrays.

>>> for key in fit_res.fit.keys():
...     print(f"{key:<15} {fit_res.fit[key].dtype} {fit_res.fit[key].shape}")

line_ids        <U14 (2,)
main_component  int16 ()
n_gauss         int16 ()
n_poly          int16 ()
wave_range      float64 (2,)
status          float64 (128, 32)
chi2            float64 (128, 32)
wavelength      float64 (128, 32, 24)
int             float64 (128, 32, 2)
err_int         float64 (128, 32, 2)
vel             float64 (128, 32, 2)
err_vel         float64 (128, 32, 2)
params          float64 (128, 32, 7)
perror          float64 (128, 32, 7)
component       int32 (7,)
param_names     <U32 (7,)
param_units     <U32 (7,)


We can extract an array of the fit parameters or intensity profile using the get_params and get_fit_profile methods. Both methods take optional keywords for selecting the component number and/or an individual pixel (using array coordinates). get_fit_profile also has a num_wavelengths keyword that allows us to interpolate the fit profile at a higher wavelength resolution than observed by EIS. The use of these methods are demonstrated in the longer example program below, which also shows one method for finding the indices and coordinates of maximum intensity.

import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.wcs.utils import wcs_to_celestial_frame
import eispac

if __name__ == '__main__':
# Read in the fit template and EIS observation
data_filepath = './eis_20190404_131513.data.h5'
template_filepath = './fe_12_195_119.2c.template.h5'

# Select a cutout of the raster
eis_frame = wcs_to_celestial_frame(data_cube.wcs)
lower_left = [None, SkyCoord(Tx=-25, Ty=225, unit=u.arcsec, frame=eis_frame)]
upper_right = [None, SkyCoord(Tx=175, Ty=425, unit=u.arcsec, frame=eis_frame)]
raster_cutout = data_cube.crop(lower_left, upper_right)

# Fit the data and save it to disk
fit_res = eispac.fit_spectra(raster_cutout, tmplt, ncpu='max')
save_filepaths = eispac.save_fit(fit_res, save_dir='cwd')

# Find indices and world coordinates of max intensity
sum_data_inten = raster_cutout.sum_spectra().data
iy, ix = np.unravel_index(sum_data_inten.argmax(), sum_data_inten.shape)
ex_world_coords = raster_cutout.wcs.array_index_to_world(iy, ix, 0)
y_arcsec, x_arcsec = ex_world_coords.Ty.value, ex_world_coords.Tx.value

# Extract data profile and interpolate fit at higher spectral resolution
data_x = raster_cutout.wavelength[iy, ix, :]
data_y = raster_cutout.data[iy, ix, :]
data_err = raster_cutout.uncertainty.array[iy, ix, :]
fit_x, fit_y = fit_res.get_fit_profile(coords=[iy,ix], num_wavelengths=100)
c0_x, c0_y = fit_res.get_fit_profile(0, coords=[iy,ix], num_wavelengths=100)
c1_x, c1_y = fit_res.get_fit_profile(1, coords=[iy,ix], num_wavelengths=100)
c2_x, c2_y = fit_res.get_fit_profile(2, coords=[iy,ix], num_wavelengths=100)

# Make a multi-panel figure with the cutout and example profile
fig = plt.figure(figsize=[10,5])

data_subplt.imshow(sum_data_inten, origin='lower', extent=cutout_extent)
data_subplt.scatter(x_arcsec, y_arcsec, color='r', marker='x')
data_subplt.set_title('Data Cutout\n'+raster_cutout.meta['mod_index']['date_obs'])
data_subplt.set_xlabel('Solar-X [arcsec]')
data_subplt.set_ylabel('Solar-Y [arcsec]')

profile_subplt.errorbar(data_x, data_y, yerr=data_err, ls='', marker='o', color='k')
profile_subplt.plot(fit_x, fit_y, color='b', label='Combined profile')
profile_subplt.plot(c0_x, c0_y, color='r', label=fit_res.fit['line_ids'])
profile_subplt.plot(c1_x, c1_y, color='r', ls='--', label=fit_res.fit['line_ids'])
profile_subplt.plot(c2_x, c2_y, color='g', label='Background')
profile_subplt.set_title(f'Cutout indices: iy = {iy}, ix = {ix}')
profile_subplt.set_xlabel('Wavelength [$\AA$]')
profile_subplt.set_ylabel('Intensity ['+raster_cutout.unit.to_string()+']')
profile_subplt.legend(loc='upper left', frameon=False)
plt.show() Fig. 5 Example data cutout (left) and fit profile (right) for the spectral window containing the Fe XII 195.119 Å line. The red X shows the location of the maximum summed intensity.

## EISMaps for Sunpy

The fit line intensities, velocities, and widths can be loaded into an EISMap, which is a subclass of sunpy.map.Map. This allow us to leverage the full power of Sunpy to do all sorts of cool science like comparing spacecraft locations, co-aligning images, reprojecting maps, and performing field extrapolations (see the Map sections of the SunPy User’s Guide and Example Gallery for some demonstrations). You can get an EISMap by either using the get_map method or saving the measurements to FITS files using export_fits and then loading them in with either eispac.EISMap(FILENAME) or even sunpy.map.Map(FILENAME) (assuming EISPAC is also imported in your program).

For now, we will just show you some examples of the quick-look plots.

>>> # Fit intensity (in a nice sunpy Map)
>>> inten_map = fit_res.get_map(component=0, measurement='intensity')
>>> inten_map.peek()

>>> # Fit velocity map