Development Guidelines

The development strategy for the DAP has been to allow for high-level interaction with the software to make subtle changes to its execution parameters. For help with the high-level interaction with the DAP parameters, see the description of its Execution.

We have also recently abstracted the input to the DAP to allow for custom DataCube objects to be provided by the user, enabling the DAP to analyze non-MaNGA data; see How to Fit a non-MaNGA Datacube.

Additionally, we have constructed the low-level, core algorithms in a way that is largely agnostic to the data being analyzed. This should allow users to use DAP methods/classes in new scripts tailored to their own needs and analysis of non-MaNGA data. For example, we demonstrate How to Fit One Spectrum using the DAP fitting modules.

We are very interested in having the community building on the existing DAP algorithms. As a means of introducing some of the nuances of its structure, the following is list of some of the guidelines we tried to follow in the development of the code:

  • Minimize code repetition by isolating repeated functionality into a flexible function or class.

  • Provide utility level classes that ease interaction with the input and output data.

  • Provide utility level classes that ease inclusion of new spectral template libraries, emission-line parameters, and spectral-index parameters for use in the main analysis classes.

  • Compartmentalize analysis into low-level functions that are independent of the specifics of the MaNGA data.

  • Generalize the interfaces of these low-level functions such that they can be incorporated into user-level scripts.

  • Where possible and appropriate, ensure the low-level functions operate on arrays of spectra that can effectively be considered independent of one another.

  • Minimize hard-coded behavior in favor of parameter or configuration files.

  • Provide an infrastructure that includes hooks for user-provided functions to override the base-level DAP analysis functions while still maintaining the output DAP datamodel.

  • Use an astropy.io.fits.HDUList object as the primary interface to the high-level classes to generalize interaction and allow them to be directly written to their reference fits files with a well-formed datamodel.

TemplateLibrary example

As an example of what some of this looks like in practice, consider the Class usage examples for a TemplateLibrary.

The TemplateLibrary object is setup to read and process the template library used to model the stellar continuum. The survey-level execution of the DAP currently uses the MILESHC library for the stellar kinematics and the MASTARHC library to model the stellar continuum in the emission-line module. The way that the DAP knows where to find the templates and how to process them is via the configuration files located in $MANGADAP_DIR/mangadap/config/spectral_templates. The one specific to the MILESHC library is $MANGADAP_DIR/mangadap/config/spectral_templates/miles_hc.ini and looks like this:

[default]
 key              = MILESHC
 file_search      = miles_cluster/*.fits
 fwhm             = 2.50
 sres_ext
 in_vacuum        = False
 wave_limit       = 3575, 7400
 lower_flux_limit = 0.0
 log10            = False

This configuration (ini) file allows one to construct the MILESHC library as follows:

tpl = TemplateLibrary('MILESHC',
                      velscale_ratio=4,     # Set the pixel size to 1/4 the MaNGA step
                      spectral_step=1e-4,   # The MaNGA step is dlogLambda = 1e-4
                      log=True,             # Sample the templates logarithmically
                      hardcopy=False)       # Don't save a hardcopy of the library to disk

This is possible because the DAP will use the keyword to identify the appropriate ini file in the $MANGADAP_DIR/mangadap/config/spectral_templates directory. If you have a template library that you want to use that’s not provided by the DAP, you can use the Analysis Plans to define a new one.

The template library key for use with the stellar-continuum fitting is one of the parameters in the StellarContinuumModelDef parameter set; see Analysis Plans.

However, you can also write scripts that incorporate the DAP functionality and circumvents the normal mangadap workflow and does not use a plan file. Assuming you have a script that uses a TemplateLibrary object, you can define a new template library in the code itself using the TemplateLibraryDef object. The TemplateLibraryDef object is actually the product of the parsed configuration file within the main DAP code. For example:

# Imports
from mangadap.proc.templatelibrary import TemplateLibraryDef, TemplateLibrary

# Define the search string for the library
search_str = '/path/to/library/*.fits'
search_sres_str = '/path/to/library/with/sres/*.fits'

# Define the template library parameters
new_tpl_lst = TemplateLibraryDef(key='MYLIB',            # Unique keyword for the library
                                 file_search=search_str, # Search string
                                 fwhm=2.50,              # FWHM of resolution element
                                 in_vacuum=False,        # Wavelength in vacuum?
                                 wave_limit=numpy.array([ 3575., 7400. ]),   # Valid Range
                                 lower_flux_limit=0.0,   # Lower limit for valid flux
                                 log10=False)            # Log binned?

# Or if you there is an extension SPECRES in *all* the files with
# the spectral resolution:
new_tpl_list = [new_tpl_list,
                TemplateLibraryDef(key='MYLIB_SRES',       # Unique library keyword
                                   file_search=search_sres_str, # Search string
                                   sres_ext='SPECRES',     # Spectral Resolution Extension
                                   in_vacuum=False,        # Wavelength in vacuum?
                                   wave_limit=numpy.array([ 3575., 7400. ]),   # Valid range
                                   lower_flux_limit=0.0,   # Lower limit for valid flux
                                   log10=False)            # Log binned?
               ]

# Read and process template library
tpl = TemplateLibrary('MYLIB',
                      tpllib_list=new_tpl_lst,  # Available list of template libraries
                      velscale_ratio=4,     # Set the pixel size to 1/4 the MaNGA step
                      spectral_step=1e-4,   # The MaNGA step is dlogLambda = 1e-4
                      log=True,             # Sample the templates logarithmically
                      hardcopy=False)       # Don't save a hardcopy of the library to disk

Adding new functionality

The pairing of the defining parameters of a specific analysis method and the instantiation of the method itself, like what is shown above, is ubiquitous in the DAP, following from one of the main design principles. Apart from allowing one to alter the details of how the DAP proceeds via changing the analysis plan file, this also facilitates incorporating new algorithms within the existing infrastructure.

The implementation of hooks for including new algorithms into the DAP infrastructure exists; however, it is minimal in some respects and hasn’t been well tested. That means that, for anyone that tries this, there are sure to be some growing pains in making sure it works properly. At the moment, a primary limitation to incorporating new algorithms is that one has to alter a number of bits of code in the DAP itself. There are a few ways to do this, but testing of these methods has been limited. The following are two sketched out examples of including new functionality or algorithms.

Warning

The code examples below are from earlier versions of the DAP (version 3.1 and earlier). They’re still reasonably relevant to this discussion, but the details will have changed. Please Submit an issue if you’re attempting to do something similar to what’s discussed below!

1. Adding a new binning scheme

The class that constructs the spatially binned spectra is SpatiallyBinnedSpectra. The method used to construct a class instance is defined using SpatiallyBinnedSpectraDef, which has components that define a set of binning parameters, a binning class instance, and/or a binning function. The existing binning schemes are:

All of these classes provide a common interface that SpatiallyBinnedSpectra calls to determine which spaxels are assigned to each bin. The format of this function must be:

def binning_function(x, y, par=None):
    # Bin the data
    ...
    return bin_id

That is, the function must take in the on-sky x and y positions of each spaxel, accept some set of parameters provided by the par dictionary and return a bin ID number associated with each x and y position.

So let’s say that you wanted to bin all spectra in a set of apertures. You could define and implement a function that performs this binning, and then execute this binning approach within SpatiallyBinnedSpectra as follows (the code is untested!):

#!/usr/bin/env python3

import time
import warnings
import numpy
import astropy.constants

from mangadap.datacube import MaNGADataCube
from mangadap.proc.reductionassessments import ReductionAssessment
from mangadap.proc.spectralstack import SpectralStackPar, SpectralStack
from mangadap.proc.spatiallybinnedspectra import SpatiallyBinnedSpectra, SpatiallyBinnedSpectraDef
from mangadap.proc.stellarcontinuummodel import StellarContinuumModel
from mangadap.proc.emissionlinemoments import EmissionLineMoments
from mangadap.proc.emissionlinemodel import EmissionLineModel
from mangadap.proc.spectralindices import SpectralIndices
from mangadap.dapfits import construct_maps_file, construct_cube_file

#-----------------------------------------------------------------------------

class ApertureBinning():
    """
    Perform aperture binning

    Args:
        x (array-like):
            List of on-sky x coordinates for apertures
        y (array-like):
            List of on-sky y coordinates for apertures
        r (array-like):
            Single or list of radii of the apertures

    Attributes:
        n (:obj:`int`):
            Number of apertures
        x (`numpy.ndarray`_):
            On-sky x coordinates for apertures
        y (`numpy.ndarray`_):
            On-sky y coordinates for apertures
        r (`numpy.ndarray`_):
            Aperture radii
    """
    def __init__(self, x, y, r):
        self.x = numpy.asarray(x)
        if len(self.x.shape) != 1:
            raise ValueError('On-sky coordinates must be one-dimensional.')
        self.n = x.size
        if len(y) != self.n:
            raise ValueError('Input coordinates are of different lengths.')
        self.y = numpy.asarray(y)
        if len(r) != 1 and len(r) != self.n:
            raise ValueError('Radii must be common to all apertures or unique to each aperture.')
        self.r = numpy.full(self.n, r, dtype=float) if len(r) == 1 else numpy.asarray(r)

    def bin_spaxels(self, x, y, par=None):
        _x = numpy.asarray(x)
        if len(_x.shape) != 1:
            raise ValueError('On-sky coordinates must be one-dimensional.')
        nspaxels = _x.size
        if len(y) != nspaxels:
            raise ValueError('Input coordinates are of different lengths.')
        _y = numpy.asarray(y)

        # Find which spaxels land in each aperture
        indx = numpy.square(_x[:,None]-self.x[None,:]) + numpy.square(_y[:,None]-self.y[None,:]) \
                    < numpy.square(self.r[None,:])
        if numpy.any(numpy.sum(indx, axis=1) > 1):
            warnings.warn('Spaxels found in multiple apertures!')

        # Return the aperture index that each spaxel is within,
        # isolating only one aperture per spaxel; spaxels not in any
        # aperture have a bin ID of -1
        binid = numpy.full((nspaxels, self.n), -1, dtype=int)
        binid[indx] = numpy.array([numpy.arange(self.n)]*nspaxels)[indx]
        return numpy.amax(binid, axis=1)

#-----------------------------------------------------------------------------
if __name__ == '__main__':
    t = time.perf_counter()

    # Set the plate, ifu, and initial velocity/redshift
    plate = 7495
    ifu = 12704
    vel = 8675.5
    nsa_redshift = vel/astropy.constants.c.to('km/s').value

    # Read the DRP LOGCUBE file
    cube = MaNGADataCube.from_plateifu(plate, ifu)

    # Calculate the S/N and coordinates
    rdxqa = ReductionAssessment('SNRG', cube)

    # Setup the aperture binning class
    ax = numpy.array([0.0, 3.0, 6.0])
    ay = numpy.array([0.0, 0.0, 0.0])
    apbin = ApertureBinning(ax, ay, 2.5)

    # Setup the stacking operations
    stackpar = SpectralStackPar('mean',         # Operation for stack
                                False,          # Apply a velocity registration
                                None,           # Velocity offsets for registration
                                'channels',     # Covariance mode and parameters
                                SpectralStack.parse_covariance_parameters('channels', 11))
    stacker = SpectralStack()

    # Create a new binning method
    binning_method = SpatiallyBinnedSpectraDef('Aperture',      # Key for binning method
                                               'ODonnell',      # Galactic reddening function
                                               3.1,             # Rv for Galactic reddening
                                               0.0,             # Minimum S/N to include
                                               None,            # Object with binning pars
                                               None,            # Binning class instance
                                               apbin.bin_spaxels,   # Binning function
                                               stackpar,        # Object with stacking pars
                                               stacker,         # Stacking class instance
                                               stacker.stack_datacube)   # Stacking function

    # Bin the spectra using the new binning method
    binned_spectra = SpatiallyBinnedSpectra('Aperture',     # Key for binning method
                                            cube,           # DRP data to bin
                                            rdxqa,          # Cube coordinates and S/N
                                            method_list=binning_method) # Binning methods

    # The rest of this is just a single execution of the remaining
    # analysis steps in
    # $MANGADAP_DIR/python/mangadap/survey/manga_dap.py , with some
    # simplifications
    stellar_continuum = StellarContinuumModel('GAU-MILESHC', binned_spectra, guess_vel=vel,
                                              guess_sig=100.)

    emission_line_moments = EmissionLineMoments('EMOMF', binned_spectra,
                                                stellar_continuum=stellar_continuum,
                                                redshift=nsa_redshift)

    emission_line_model = EmissionLineModel('EFITF', binned_spectra,
                                            stellar_continuum=stellar_continuum,
                                            redshift=nsa_redshift, dispersion=100.0)

    spectral_indices = SpectralIndices('INDXEN', binned_spectra, redshift=nsa_redshift,
                                       stellar_continuum=stellar_continuum,
                                       emission_line_model=emission_line_model)

    construct_maps_file(cube, rdxqa=rdxqa, binned_spectra=binned_spectra,
                        stellar_continuum=stellar_continuum,
                        emission_line_moments=emission_line_moments,
                        emission_line_model=emission_line_model,
                        spectral_indices=spectral_indices, nsa_redshift=nsa_redshift)

    construct_cube_file(cube, binned_spectra=binned_spectra,
                        stellar_continuum=stellar_continuum,
                        emission_line_model=emission_line_model)

    print('Elapsed time: {0} seconds'.format(time.perf_counter() - t))

You’ll notice that there are some limitations in what one can implement. In this example, the limitation is that each spaxel must be assigned to a single unique bin, not multiple bins if the apertures are overlapping.

2. Adding a new emission-line fitter

Adding a new binning scheme is relatively straight-forward because all that’s required is to provide a new binning function that adheres to the specified form. Things become much more complicated when you want to replace a core algorithm that provides much of the content of the output data model. Still, it can be done, it’s just that one has to follow more requirements that can be more stringent.

Let’s say you want to add a new emission-line fitter (as we did for MPL-6 in changing from Elric to Sasuke). The class that constructs the parameterized emission-line models is EmissionLineModel. The method used to construct a class instance is defined using EmissionLineModelDef, which has components that define a set of model-fitting parameters, a model-fitting class instance, and/or a model-fitting function. The common function call that any emission-line fitter must provide looks like:

def fit(binned_spectra, par=None, loggers=None, quiet=False):
    # Fit the spectra
    ...
    return model_eml_flux, model_eml_base, model_eml_mask, model_fit_par, \
            model_eml_par, model_binid

where binned_spectra is a SpatiallyBinnedSpectra object and the returned arrays are:

  • model_eml_flux: Model emission-line flux only; shape is \((N_{\rm mod}, N_{\rm wave})\). The first axis is ordered by model ID number.

  • model_eml_base: Any baseline resulting from the emission-line fit such that the model fit to each spectrum is: stellar_continuum + model_eml_flux + model_eml_base; shape is \((N_{\rm mod}, N_{\rm wave})\). The first axis is ordered by model ID number.

  • model_eml_mask: Boolean or bit-mask array for fitted models; shape is \((N_{\rm mod}, N_{\rm wave})\). The first axis is ordered by model ID number.

  • model_fit_par: A numpy record array that provides the results of each fit. This can be None and the output maps file will still be successfully written

  • model_eml_par: A numpy record array that provides the output model parameters. The data type must be as returned by mangadap.proc.spectralfitting.EmissionLineFit._per_emission_line_dtype() and the shape must be \((N_{\rm mod},)\) with the parameters ordered by the model ID number.

  • model_binid: A 2D map of the ID numbers assigned to each spaxel with a fitted model. Any spaxel without an emission-line model should have model_binid = -1, and the number of IDs that are greater than -1 must be \(N_{\rm mod}\). The shape is \((N_x,N_y)\); i.e., it must match the spatial dimensions of the fitted DRP data cube. This can be {{{None}}}, which indicates that the bin IDS are the same as the bin IDs set in the binned_spectra object.

In addition to the high-level data model interfaces, like EmissionLineModel, the DAP attempts to provide a set of base classes that provided functionality common to a set of abstracted spectral-fitting routines. For the emission-line fitting, this is the EmissionLineFit class. This object provides the data table description that should be common to all emission-line model output for the construction of the output data model by EmissionLineModel (see model_eml_par above).

Warning

The code discussed/provided below was an initial go at sketching out the code for the emission-line module used in MPL-6. The code will not work out of the box and is meant to illustrate the solution to the problem. For the actual solution, see the main DAP interface class Sasuke and the primary fitting function written by Xihan Ji and Michele Cappellari (with some significant edits by Kyle Westfall), mangadap.contrib.xjmc.

So, let’s say you have a function first_second_iteration that applies the new fitting approach that we want for the emission lines. The following is untested code that could be used to implement this function as the emission-line model fitter while still providing the same DAP output data model:

#!/usr/bin/env python3

import time
import warnings
import numpy
import astropy.constants

from mangadap.datacube import MaNGADataCube
from mangadap.proc.reductionassessments import ReductionAssessment
from mangadap.proc.spatiallybinnedspectra import SpatiallyBinnedSpectra
from mangadap.proc.stellarcontinuummodel import StellarContinuumModel
from mangadap.proc.ppxffit import PPXFFit
from mangadap.util.instrument import spectrum_velocity_scale
from mangadap.util.fitsutil import DAPFitsUtil
from mangadap.util.fileio import init_record_array
from mangadap.proc.emissionlinemoments import EmissionLineMoments
from mangadap.proc.spectralfitting import EmissioneLineFit
from mangadap.proc.emissionlinemodel import EmissionLineModelDef, EmissionLineModel
from mangadap.par.emissionlinedb import EmissionLineDB
from mangadap.proc.spectralindices import SpectralIndices
from mangadap.dapfits import construct_maps_file, construct_cube_file

from mangadap.contrib.xjmc import first_second_iteration

#-----------------------------------------------------------------------------

class XJMCEmissionLineFitter(EmissionLineFit):
    def __init__(self, par=None):
        if par is None:
            # Set the default parameter set.  The guess_redshift,
            # stellar_continuum, and emission_lines values can be
            # filled by EmissionLineModel._fill_method_par()
            par = { 'guess_redshift': None,     # The guess redshift for each binned spectrum
                    'stellar_continuum': None,  # The StellarContinuumModel object
                    'emission_lines': None,     # The EmissionLineDB object
                    'degree': 8,                # Additive polynomial order
                    'mdegree': 0 }              # Multiplicative polynomial order
        EmissionLineFit.__init__(self, 'XJMC', None, par=par)


    def fit(self, binned_spectra, par=None, loggers=None, quiet=False):
        if par is not None:
            self.par = par

        # Check the parameter keys
        required_keys = [ 'guess_redshift', 'stellar_continuum', 'emission_lines', 'degree',
                          'mdegree' ]
        if numpy.any([ reqk not in self.par.keys() for reqk in required_keys ]):
            raise ValueError('Parameter dictionary does not have all the required keys.')

        # Wavelengths are in vacuum
        wave = binned_spectra['WAVE'].data.copy()
        # Velocity step per pixel
        velscale = spectrum_velocity_scale(wave)
        # Flux and noise masked arrays; shape is (Nspaxels,Nwave)
        # where Nspaxels is Nx*Ny
        flux = binned_spectra.cube.copy_to_masked_array(flag=['DONOTUSE', 'FORESTAR'])
        noise = numpy.ma.power(binned_spectra.cube.copy_to_masked_array(
                                    attr='ivar', flag=['DONOTUSE', 'FORESTAR']), -0.5)
        # Spaxel coordinates; shape is (nspaxels,)
        x = binned_spectra.rdxqa['SPECTRUM'].data['SKY_COO'][:,0]
        y = binned_spectra.rdxqa['SPECTRUM'].data['SKY_COO'][:,1]
        # Binned flux and binned noise masked arrays; shape is (nbins,nwave)
        flux_binned = binned_spectra.copy_to_masked_array(
                                                    flag=binned_spectra.do_not_fit_flags())
        noise_binned = numpy.ma.power(binned_spectra.copy_to_masked_array(ext='IVAR',
                                                flag=binned_spectra.do_not_fit_flags()) , -0.5)
        # Bin coordinates; shape is (nbins,)
        x_binned = binned_spectra['BINS'].data['SKY_COO'][:,0]
        y_binned = binned_spectra['BINS'].data['SKY_COO'][:,1]
        # Set initial guesses for the velocity and velocity
        # dispersion
        if self.par['guess_redshift'] is not None:
            # Use guess_redshift if provided
            vel = self.par['guess_redshift'] * astropy.constants.c.to('km/s').value
            # And set default velocity dispersion to 100 km/s
            sig = numpy.full(vel.size, 100, dtype=float)
        elif self.par['stellar_continuum'] is not None:
            # Otherwise use the stellar-continuum result
            vel, sig = self.par['stellar_continuum'].matched_guess_kinematics(binned_spectra,
                                                                              cz=True)
        else:
            # TODO: Set default guess kinematics
            vel, sig = None, None

        # Get the stellar templates;
        # shape is (Ntemplates, Nwave_templates)
        if self.par['stellar_continuum'] is not None:
            stars_templates = self.par['stellar_continuum'].method['fitpar']['template_library']
            stars_templates_wave = stars_templates['WAVE'].data.copy()
            stars_templates = stars_templates['FLUX'].data.copy()
            velscale_ratio = self.par['stellar_continuum'].method['fitpar']['velscale_ratio']
            dv = -PPXFFit.ppxf_tpl_obj_voff(stars_templates_wave, wave, velscale,
                                            velscale_ratio=velscale_ratio)
        else:
            # TODO: Default construction of stellar templates
            stars_templates = None
            velscale_ratio = None
            dv = None

        # TODO: Construct gas templates
        gas_templates = None
        gas_names = None

        # TODO: Default polynomial orders
        degree = 8 if self.par['degree'] is None else self.par['degree']
        mdegree = 0 if self.par['mdegree'] is None else self.par['mdegree']

        # Output is:
        #   - model_flux: stellar-continuum + emission-line model;
        #     shape is (Nmod, Nwave); first axis is ordered by model
        #     ID number
        #   - model_eml_flux: model emission-line flux only; shape
        #     is (Nmod, Nwave); first axis is ordered by model ID
        #     number
        #   - model_mask: boolean or bit mask for fitted models;
        #     shape is (Nmod, Nwave); first axis is ordered by model
        #     ID number
        #   - model_binid: ID numbers assigned to each spaxel with a
        #     fitted model; any spaxel without a model should have
        #     model_binid = -1; the number of >-1 IDs must be Nmod;
        #     shape is (Nx,Ny) which is equivalent to:
        #       flux[:,0].reshape((numpy.sqrt(Nspaxels).astype(int),)*2).shape
        #   - eml_flux: Flux of each emission line; shape is
        #     (Nmod,Neml)
        #   - eml_fluxerr: Error in emission-line fluxes; shape is
        #     (Nmod, Neml)
        #   - eml_kin: Kinematics (velocity and velocity dispersion)
        #     of each emission line; shape is (Nmod,Neml,Nkin)
        #   - eml_kinerr: Error in the kinematics of each emission
        #     line
        #   - eml_sigmacorr: Quadrature corrections required to
        #     obtain the astrophysical velocity dispersion; shape is
        #     (Nmod,Neml); corrections are expected to be applied as
        #     follows:
        #       sigma = numpy.ma.sqrt( numpy.square(eml_kin[:,:,1])
        #                               - numpy.square(eml_sigmacorr))
        model_flux, model_eml_flux, model_mask, model_binid, eml_flux, eml_fluxerr, \
                eml_kin, eml_kinerr, eml_sigmacorr \
                        = first_second_iteration(wave, flux, noise, flux_binned, noise_binned,
                                                 velscale, velscale_ratio, dv, vel, sig,
                                                 stars_templates, gas_templates, gas_names,
                                                 degree, mdegree, x, y, x_binned, y_binned)

        # The ordered indices in the flatted bin ID map with/for
        # each model
        model_srt = numpy.argsort(model_binid.ravel())[model_binid.ravel() > -1]

        # Construct the output emission-line database.  The data
        # type defined by
        # EmissionLineFit._per_emission_line_dtype(); shape is
        # (Nmod,); parameters must be ordered by model ID number
        nmod = len(model_srt)
        neml = eml_flux.shape[1]
        nkin = eml_kin.shape[-1]
        model_eml_par = init_record_array(nmod,
                            EmissionLineFit._per_emission_line_dtype(neml, nkin, numpy.int16))
        model_eml_par['BINID'] = model_binid.ravel()[model_srt]
        model_eml_par['BINID_INDEX'] = numpy.arange(nmod)
        model_eml_par['MASK'][:,:] = 0
        model_eml_par['FLUX'] = eml_flux
        model_eml_par['FLUXERR'] = eml_fluxerr
        model_eml_par['KIN'] = eml_kin
        model_eml_par['KINERR'] = eml_kinerr
        model_eml_par['SIGMACORR'] = eml_sigmacorr

        # Include the equivalent width measurements
        if self.par['emission_lines'] is not None:
            EmissionLineFit.measure_equivalent_width(wave, flux[model_srt,:],
                                                     par['emission_lines'], model_eml_par)

        # Calculate the "emission-line baseline" as the difference
        # between the stellar continuum model determined for the
        # kinematics and the one determined by the optimized
        # stellar-continuum + emission-line fit:
        if self.par['stellar_continuum'] is not None:
            # Construct the full 3D cube for the stellar continuum
            # models
            sc_model_flux, sc_model_mask \
                    = DAPFitsUtil.reconstruct_cube(binned_spectra.shape,
                                                   self.par['stellar_continuum']['BINID'].data,
                                               [ self.par['stellar_continuum']['FLUX'].data,
                                                 self.par['stellar_continuum']['MASK'].data ])
            # Set any masked pixels to 0
            sc_model_flux[sc_model_mask>0] = 0.0

            # Construct the full 3D cube of the new stellar
            # continuum from the combined stellar-continuum +
            # emission-line fit
            el_continuum = DAPFitsUtil.reconstruct_cube(binned_spectra.shape, model_binid,
                                                        model_flux - model_eml_flux)
            # Get the difference, restructure it to match the shape
            # of the emission-line models, and zero any masked
            # pixels
            model_eml_base = (el_model_flux - sc_model_flux).reshape(-1,wave.size)[model_srt,:]
            if model_mask is not None:
                model_eml_base[model_mask>0] = 0.0
        else:
            model_eml_base = numpy.zeros(model_flux.shape, dtype=float)

        # Returned arrays are:
        #   - model_eml_flux: model emission-line flux only; shape
        #     is (Nmod, Nwave); first axis is ordered by model ID
        #     number
        #   - model_eml_base: difference between the combined fit
        #     and the stars-only fit; shape is (Nmod, Nwave); first
        #     axis is ordered by model ID number
        #   - model_mask: boolean or bit mask for fitted models;
        #     shape is (Nmod, Nwave); first axis is ordered by model
        #     ID number
        #   - model_fit_par: This provides the results of each fit;
        #     TODO: The is set to None.  Provide metrics of the ppxf
        #     fit to each spectrum?
        #   - model_eml_par: output model parameters; data type must
        #     be EmissionLineFit._per_emission_line_dtype(); shape
        #     is (Nmod,); parameters must be ordered by model ID
        #     number
        #   - model_binid: ID numbers assigned to each spaxel with a
        #     fitted model; any spaxel with a model should have
        #     model_binid = -1; the number of >-1 IDs must be Nmod;
        #     shape is (Nx,Ny)
        return model_eml_flux, model_eml_base, model_mask, None, model_eml_par, model_binid


#-----------------------------------------------------------------------------
if __name__ == '__main__':
    t = time.perf_counter()

    # Set the plate, ifu, and initial velocity/redshift
    plate = 7495
    ifu = 12704
    vel = 8675.5
    nsa_redshift = vel/astropy.constants.c.to('km/s').value

    # Read the DRP LOGCUBE file
    cube = MaNGADataCube.from_plateifu(plate, ifu)

    # Calculate the S/N and coordinates
    rdxqa = ReductionAssessment('SNRG', cube)

    # Peform the Voronoi binning to S/N>~10
    binned_spectra = SpatiallyBinnedSpectra('VOR10', cube, rdxqa)

    # Fit the stellar kinematics
    stellar_continuum = StellarContinuumModel('GAU-MILESHC', binned_spectra, guess_vel=vel,
                                              guess_sig=100.)

    # Get the emission-line moments
    emission_line_moments = EmissionLineMoments('EMOMF', binned_spectra,
                                                stellar_continuum=stellar_continuum,
                                                redshift=nsa_redshift)

    # Get an estimate of the redshift of each bin using the first
    # moment of the H-alpha emission line:
    el_init_redshift = numpy.full(binned_spectra.nbins, nsa_redshift, dtype=float)
    # HARDCODED FOR A SPECIFIC EMISSION-LINE MOMENT DATABASE
    halpha_channel = 7
    halpha_mom1_masked = emission_line_moments['ELMMNTS'].data['MASK'][:,halpha_channel] > 0
    # - Use the 1st moment of the H-alpha line
    el_init_redshift[ emission_line_moments['ELMMNTS'].data['BINID_INDEX'] ] \
                = emission_line_moments['ELMMNTS'].data['MOM1'][:,halpha_channel] \
                                / astropy.constants.c.to('km/s').value
    # - For missing bins in the moment measurements and bad H-alpha
    #   moment measurements, use the value for the nearest good bin
    bad_bins = numpy.append(emission_line_moments.missing_bins,
                emission_line_moments['ELMMNTS'].data['BINID'][halpha_mom1_masked]).astype(int)
    if len(bad_bins) > 0:
        nearest_good_bin_index = binned_spectra.find_nearest_bin(bad_bins, indices=True)
        bad_bin_index = binned_spectra.get_bin_indices(bad_bins)
        el_init_redshift[bad_bin_index] = el_init_redshift[nearest_good_bin_index]

    # Setup the new emission-line fitter
    fitter = XJMCEmissionLineFitter()

    # Setup the new fitting method
    fit_method = EmissionLineModelDef('XJMC',       # Key for the fitting method
                                      0.0,          # Minimum S/N of the binned spectra
                                      None,         # Keyword for an artifact mask
                                      None,         # Keyword for an emission-line database
                                      fitter.par,   # Object with fit parameters
                                      fitter,       # Fitting class instance
                                      fitter.fit)   # Fitting function

    # Fit the emission lines
    emission_line_model = EmissionLineModel('XJMC',
                                            binned_spectra,
                                            stellar_continuum=stellar_continuum,
                                            redshift=el_init_redshift, dispersion=100.0,
                                            method_list=fit_method)

    # The rest of this is just a single execution of the remaining
    # analysis steps in
    # $MANGADAP_DIR/python/mangadap/survey/manga_dap.py , with some
    # simplifications
    spectral_indices = SpectralIndices('INDXEN', binned_spectra, redshift=nsa_redshift,
                                       stellar_continuum=stellar_continuum,
                                       emission_line_model=emission_line_model)

    construct_maps_file(cube, rdxqa=rdxqa, binned_spectra=binned_spectra,
                        stellar_continuum=stellar_continuum,
                        emission_line_moments=emission_line_moments,
                        emission_line_model=emission_line_model,
                        spectral_indices=spectral_indices, nsa_redshift=nsa_redshift)

    construct_cube_file(cube, binned_spectra=binned_spectra,
                        stellar_continuum=stellar_continuum,
                        emission_line_model=emission_line_model)

    print('Elapsed time: {0} seconds'.format(time.perf_counter() - t))