text: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function
  import os
  import copy
  import shutil
  import collections
  import logging
  import tempfile
  import filecmp
  import time
  import json
  import numpy as np
  from astropy.io import fits
  from astropy.table import Table, Column, vstack
  from gammapy.maps import Map, HpxGeom, WcsGeom, MapAxis, WcsNDMap, HpxNDMap
  import fermipy
  import fermipy.defaults as defaults
  import fermipy.utils as utils
  import fermipy.wcs_utils as wcs_utils
  import fermipy.fits_utils as fits_utils
  import fermipy.gtutils as gtutils
  import fermipy.srcmap_utils as srcmap_utils
  import fermipy.skymap as skymap
  import fermipy.plotting as plotting
  import fermipy.irfs as irfs
  import fermipy.sed as sed
  import fermipy.lightcurve as lightcurve
  from fermipy.residmap import ResidMapGenerator
  from fermipy.tsmap import TSMapGenerator, TSCubeGenerator
  from fermipy.sourcefind import SourceFind
  from fermipy.extension import ExtensionFit
  from fermipy.utils import merge_dict
  from fermipy.utils import create_hpx_disk_region_string
  from fermipy.utils import resolve_file_path
  from fermipy.roi_model import ROIModel
  from fermipy.ltcube import LTCube
  from fermipy.plotting import AnalysisPlotter
  from fermipy.logger import Logger, log_level
  from fermipy.config import ConfigSchema
  from fermipy.timing import Timer
  from fermipy.docstring_utils import DocstringMeta
  from fermipy.fitcache import FitCache
  from fermipy.data_struct import MutableNamedTuple
  # pylikelihood
  import GtApp
  import FluxDensity
  from LikelihoodState import LikelihoodState
  from fermipy.gtutils import BinnedAnalysis, SummedLikelihood
  import BinnedAnalysis as ba
  import pyLikelihood as pyLike
  
  norm_parameters = {
      'ConstantValue': ['Value'],
      'PowerLaw': ['Prefactor'],
      'PowerLaw2': ['Integral'],
      'BrokenPowerLaw': ['Prefactor'],
      'SmoothBrokenPowerLaw': ['Prefactor'],
      'LogParabola': ['norm'],
      'PLSuperExpCutoff': ['Prefactor'],
      'PLSuperExpCutoff2': ['Prefactor'],
      'ExpCutoff': ['Prefactor'],
      'FileFunction': ['Normalization'],
      'DMFitFunction': ['sigmav'],
      'Gaussian': ['Prefactor'],
      'PiecewisePowerLaw': ['dNdE0'],
  }
  
  shape_parameters = {
      'ConstantValue': [],
      'PowerLaw': ['Index'],
      'PowerLaw2': ['Index'],
      'BrokenPowerLaw': ['Index1', 'Index2'],
      'SmoothBrokenPowerLaw': ['Index1', 'Index2'],
      'LogParabola': ['alpha', 'beta'],
      'PLSuperExpCutoff': ['Index1', 'Cutoff'],
      'PLSuperExpCutoff2': ['Index1', 'Expfactor'],
      'ExpCutoff': ['Index1', 'Cutoff'],
      'FileFunction': [],
      'DMFitFunction': ['mass'],
      'Gaussian': ['Mean', 'Sigma'],
  }
  
  index_parameters = {
      'ConstantValue': [],
      'PowerLaw': ['Index'],
      'PowerLaw2': ['Index'],
      'BrokenPowerLaw': ['Index1', 'Index2'],
      'SmoothBrokenPowerLaw': ['Index1', 'Index2'],
      'LogParabola': ['alpha', 'beta'],
      'PLSuperExpCutoff': ['Index1', 'Index2'],
      'PLSuperExpCutoff2': ['Index1', 'Index2'],
      'ExpCutoff': ['Index1'],
      'FileFunction': [],
      'DMFitFunction': [],
      'Gaussian': [],
  }
  
  
  def set_edisp_kwargs(app, config, kw):
      gtapp = GtApp.GtApp(app)
      if 'edisp_bins' in gtapp.pars.keys():
          if config['gtlike']['edisp']:
              edisp_bins = config['gtlike']['edisp_bins']
          else:
              edisp_bins = 0
          kw['edisp_bins'] = edisp_bins
      elif 'edisp' in gtapp.pars.keys():
          kw['edisp'] = config['gtlike']['edisp']
  
  
  def make_scaled_srcmap(roi, srcmap0,
                         bexp_file0, bexproi_file0,
                         bexp_file1, bexproi_file1,
                         ccubefile1,
                         outfile):
  
      bexp0 = Map.read(bexp_file0)
      bexp1 = Map.read(bexp_file1)
      bexproi0 = Map.read(bexproi_file0)
      bexproi1 = Map.read(bexproi_file1)
      bexp_ratio = WcsNDMap(bexp0.geom,
                            bexp1.data / bexp0.data)
      bexproi_ratio = WcsNDMap(bexproi0.geom,
                               bexproi1.data / bexproi0.data)
      hdulist = fits.open(srcmap0)
      for src in roi.sources:
  
          if src.diffuse:
              ratio = bexproi_ratio.data
              hdulist[src.name].data *= ratio
          else:
              coord = (src.skydir,
                       bexp_ratio.geom.axes[0].center)
              ratio = bexp_ratio.interp_by_coord(coord)
              hdulist[src.name].data *= ratio[:, None, None]
  
      hdulist_ccube = fits.open(ccubefile1)
      hdulist['PRIMARY'] = hdulist_ccube['PRIMARY']
      hdulist['EBOUNDS'] = hdulist_ccube['EBOUNDS']
      hdulist['GTI'] = hdulist_ccube['GTI']
      hdulist.writeto(outfile, overwrite=True)
      hdulist_ccube.close()
      hdulist.close()
  
  
  def create_sc_table(scfile, colnames=None):
      """Load an FT2 file from a file or list of files."""
  
      if utils.is_fits_file(scfile) and colnames is None:
          return create_table_from_fits(scfile, 'SC_DATA')
  
      if utils.is_fits_file(scfile):
          files = [scfile]
      else:
          files = [line.strip() for line in open(scfile, 'r')]
  
      tables = [create_table_from_fits(f, 'SC_DATA', colnames)
                for f in files]
  
      return vstack(tables)
  
  
  def create_table_from_fits(fitsfile, hduname, colnames=None):
      """Memory efficient function for loading a table from a FITS
      file."""
  
      if colnames is None:
          return Table.read(fitsfile, hduname)
  
      cols = []
      with fits.open(fitsfile, memmap=True) as h:
          for k in colnames:
              data = h[hduname].data.field(k)
              cols += [Column(name=k, data=data)]
      return Table(cols)
  
  
  def get_spectral_index(src, egy):
      """Compute the local spectral index of a source."""
      delta = 1E-5
      f0 = src.spectrum()(pyLike.dArg(egy * (1 - delta)))
      f1 = src.spectrum()(pyLike.dArg(egy * (1 + delta)))
  
      if f0 > 0 and f1 > 0:
          gamma = np.log10(f0 / f1) / np.log10((1 - delta) / (1 + delta))
      else:
          gamma = np.nan
  
      return gamma
  
  
  def run_gtapp(appname, logger, kw, **kwargs):
  
      loglevel = kwargs.get('loglevel', logging.INFO)
  
      logger.log(loglevel, 'Running %s.', appname)
      t0 = time.time()
      filter_dict(kw, None)
      kw = utils.unicode_to_str(kw)
      gtapp = GtApp.GtApp(str(appname))
  
      for k, v in kw.items():
  
          if (appname == 'gtbin' and k == 'scfile' and
              not v.startswith('@') and
                  not utils.is_fits_file(v)):
              v = '@' + v
          gtapp[k] = v
  
      logger.log(loglevel, gtapp.command())
      stdin, stdout = gtapp.runWithOutput(print_command=False)
  
      for line in stdout:
          logger.log(loglevel, line.strip())
  
          # Capture return code?
  
      t1 = time.time()
      logger.log(loglevel, 'Finished %s. Execution time: %.2f s',
                 appname, t1 - t0)
  
  
  def filter_dict(d, val):
      for k, v in d.items():
          if v == val:
              del d[k]
  
  
  class GTAnalysis(fermipy.config.Configurable, sed.SEDGenerator,
                   ResidMapGenerator, TSMapGenerator, TSCubeGenerator,
                   SourceFind, ExtensionFit, lightcurve.LightCurve):
      """High-level analysis interface that manages a set of analysis
      component objects.  Most of the functionality of the Fermipy
      package is provided through the methods of this class.  The class
      constructor accepts a dictionary that defines the configuration
      for the analysis.  Keyword arguments to the constructor can be
      used to override parameters in the configuration dictionary.
      """
  
      __metaclass__ = DocstringMeta
  
      _docstring_registry = {
          'extension': defaults.extension,
          'sed': defaults.sed,
          'localize': defaults.localize,
          'tsmap': defaults.tsmap,
          'residmap': defaults.residmap,
          'lightcurve': defaults.lightcurve,
          'find_sources': defaults.sourcefind,
      }
  
      defaults = {'logging': defaults.logging,
                  'fileio': defaults.fileio,
                  'optimizer': defaults.optimizer,
                  'binning': defaults.binning,
                  'selection': defaults.selection,
                  'model': defaults.model,
                  'data': defaults.data,
                  'ltcube': defaults.ltcube,
                  'gtlike': defaults.gtlike,
                  'mc': defaults.mc,
                  'residmap': defaults.residmap,
                  'tsmap': defaults.tsmap,
                  'tscube': defaults.tscube,
                  'sourcefind': defaults.sourcefind,
                  'sed': defaults.sed,
                  'lightcurve': defaults.lightcurve,
                  'extension': defaults.extension,
                  'localize': defaults.localize,
                  'roiopt': defaults.roiopt,
                  'plotting': defaults.plotting,
                  'components': (None, '', list)}
  
      def __init__(self, config, roi=None, **kwargs):
  
          # Setup directories
          self._rootdir = os.getcwd()
          self._outdir = None
          validate = kwargs.pop('validate', True)
          self._loglevel = kwargs.pop('loglevel', logging.INFO)
  
          super(GTAnalysis, self).__init__(config, validate=validate,
                                           **kwargs)
  
          self._projtype = self.config['binning']['projtype']
          self._tmin = self.config['selection']['tmin']
          self._tmax = self.config['selection']['tmax']
          self._lck_params = {}
  
          # Set random seed
          np.random.seed(self.config['mc']['seed'])
  
          # Destination directory for output data products
          if self.config['fileio']['outdir'] is not None:
              self._outdir = os.path.join(self._rootdir,
                                          self.config['fileio']['outdir'])
              utils.mkdir(self._outdir)
          else:
              raise Exception('Save directory not defined.')
  
          if self.config['fileio']['logfile'] is None:
              self.config['fileio']['logfile'] = os.path.join(self.outdir,
                                                              'fermipy.log')
  
          self._logger_name = self.config['logging']['prefix'] + \
              self.__class__.__name__
          Logger.configure(self._logger_name,
                           self.config['fileio']['logfile'],
                           log_level(self.config['logging']['verbosity']))
  
          self.logger.log(self.loglevel, '\n' + '-' * 80 + '\n' +
                          "fermipy version {} ".
                          format(fermipy.__version__) + '\n' +
                          'ScienceTools version %s', fermipy.get_st_version())
          self.print_config(self.logger, loglevel=logging.DEBUG)
  
          # Working directory (can be the same as savedir)
          if self.config['fileio']['usescratch']:
              self.config['fileio']['workdir'] = tempfile.mkdtemp(
                  prefix=os.environ['USER'] + '.',
                  dir=self.config['fileio']['scratchdir'])
              self.logger.info('Created working directory: %s',
                               self.config['fileio']['workdir'])
          else:
              self._config['fileio']['workdir'] = self._outdir
  
          if 'FERMIPY_WORKDIR' not in os.environ:
              os.environ['FERMIPY_WORKDIR'] = self.config['fileio']['workdir']
  
          # put pfiles into savedir
          os.environ['PFILES'] = \
              self.workdir + ';' + os.environ['PFILES'].split(';')[-1]
  
          # Create Plotter
          self._plotter = AnalysisPlotter(self.config['plotting'],
                                          fileio=self.config['fileio'],
                                          logging=self.config['logging'])
  
          self._like = None
          self._components = []
          configs = self._create_component_configs()
  
          # Setup the ROI definition
          if roi is None:
              self._roi = ROIModel.create(configs[0]['selection'],
                                          configs[0]['model'],
                                          fileio=configs[0]['fileio'],
                                          coordsys=configs[0]['binning']['coordsys'])
          else:
              self._roi = roi
  
          for cfg in configs:
              comp = self._create_component(cfg, loglevel=self.loglevel)
              self._components.append(comp)
  
          self._files = {}
          self._files['ccube'] = os.path.join(self.workdir, 'ccube.fits')
  
          ebin_edges = np.zeros(0)
          roiwidths = np.zeros(0)
          binsz = np.zeros(0)
          for c in self.components:
              ebin_edges = np.concatenate((ebin_edges, c.log_energies))
              roiwidths = np.insert(roiwidths, 0, c.roiwidth)
              binsz = np.insert(binsz, 0, c.binsz)
  
          self._ebin_edges = np.sort(np.unique(ebin_edges.round(5)))
          self._enumbins = len(self._ebin_edges) - 1
          self._loge_bounds = np.array([self._ebin_edges[0],
                                        self._ebin_edges[-1]])
  
          self._roi_data = {
              'loglike': np.nan,
              'npred': 0.0,
              'npred_wt': 0.0,
              'counts': np.zeros(self.enumbins),
              'counts_wt': np.zeros(self.enumbins),
              'model_counts': np.zeros(self.enumbins),
              'model_counts_wt': np.zeros(self.enumbins),
              'energies': np.copy(self.energies),
              'log_energies': np.copy(self.log_energies),
              'loge_bounds': np.copy(self.loge_bounds),
              'components': []
          }
  
          for c in self._components:
              comp_model = [{'loglike': np.nan,
                             'npred': 0.0,
                             'npred_wt': 0.0,
                             'counts': np.zeros(c.enumbins),
                             'counts_wt': np.zeros(c.enumbins),
                             'model_counts': np.zeros(c.enumbins),
                             'model_counts_wt': np.zeros(c.enumbins),
                             'energies': np.copy(c.energies),
                             'log_energies': np.copy(c.log_energies),
                             'src_expscale': copy.deepcopy(c.src_expscale),
                             }]
  
              self._roi_data['components'] += comp_model
  
          self._roiwidth = max(roiwidths)
          self._binsz = min(binsz)
          self._npix = int(np.round(self._roiwidth / self._binsz))
  
          axes = [MapAxis.from_edges(self.energies, interp='log',
                                     name='energy', unit='MeV')]
          if self.projtype == 'HPX':
              is_nested = self.config['binning']['hpx_ordering_scheme'] == "NESTED"
              self._geom = HpxGeom.create(2**self.config['binning']['hpx_order'],
                                          nest=is_nested,
                                          coordsys=self.config['binning']['coordsys'],
                                          width=self.config['binning']['roiwidth'],
                                          skydir=self.roi.skydir,
                                          axes=axes)
          else:
  
              self._geom = WcsGeom.create(npix=self.npix, binsz=self._binsz,
                                          coordsys=self.config['binning']['coordsys'],
                                          proj=self.config['binning']['proj'],
                                          skydir=self.roi.skydir,
                                          axes=axes)
              
              # Update projection of ROI object
              self._roi.set_geom(self.geom)
  
          if self.config['fileio']['usescratch']:
              self.stage_input()
  
      def __del__(self):
          self.stage_output()
          self.cleanup()
  
      @property
      def logger(self):
          """Return the default loglevel."""
          return logging.getLogger(self._logger_name)
  
      @property
      def loglevel(self):
          """Return the default loglevel."""
          return self._loglevel
  
      @property
      def workdir(self):
          """Return the analysis working directory."""
          return self.config['fileio']['workdir']
  
      @property
      def outdir(self):
          """Return the analysis output directory."""
          return self._outdir
  
      @property
      def roi(self):
          """Return the ROI object."""
          return self._roi
  
      @property
      def plotter(self):
          """Return the plotter instance."""
          return self._plotter
  
      @property
      def like(self):
          """Return the global likelihood object."""
          return self._like
  
      @property
      def components(self):
          """Return the list of analysis components."""
          return self._components
  
      @property
      def energies(self):
          """Return the energy bin edges in MeV."""
          return 10**self._ebin_edges
  
      @property
      def log_energies(self):
          """Return the energy bin edges in log10(E/MeV)."""
          return self._ebin_edges
  
      @property
      def enumbins(self):
          """Return the number of energy bins."""
          return self._enumbins
  
      @property
      def npix(self):
          """Return the number of energy bins."""
          return self._npix
  
      @property
      def loge_bounds(self):
          """Current analysis energy bounds in log10(E/MeV)."""
          return self._loge_bounds
  
      @property
      def projtype(self):
          """Return the type of projection to use"""
          return self._projtype
  
      @property
      def geom(self):
          """ROI geometry."""
          return self._geom
  
      @property
      def tmin(self):
          """Return the MET time for the start of the observation."""
          return self._tmin
  
      @property
      def tmax(self):
          """Return the MET time for the end of the observation."""
          return self._tmax
  
      @property
      def files(self):
          return self._files
  
      @classmethod
      def create(cls, infile, config=None, params=None, mask=None):
          """Create a new instance of GTAnalysis from an analysis output file
          generated with `~fermipy.GTAnalysis.write_roi`.  By default
          the new instance will inherit the configuration of the saved
          analysis instance.  The configuration may be overriden by
          passing a configuration file path with the ``config``
          argument.
  
          Parameters
          ----------
  
          infile : str
              Path to the ROI results file.
  
          config : str
              Path to a configuration file.  This will override the
              configuration in the ROI results file.
  
          params : str
              Path to a yaml file with updated parameter values
          
          mask : str
              Path to a fits file with an updated mask
  
          """
  
          infile = os.path.abspath(infile)
          roi_file, roi_data = utils.load_data(infile)
  
          if config is None:
              config = roi_data['config']
              validate = False
          else:
              validate = True
  
          gta = cls(config, validate=validate)
          gta.setup(init_sources=False)
          gta.load_roi(infile, params=params, mask=mask)
          return gta
  
      def clone(self, config, **kwargs):
          """Make a clone of this analysis instance."""
          gta = GTAnalysis(config, **kwargs)
          gta._roi = copy.deepcopy(self.roi)
          return gta
  
      def set_random_seed(self, seed):
          """Set the seed for the random number generator"""
          self.config['mc']['seed'] = seed
          np.random.seed(seed)
  
      def set_log_level(self, level):
          self.logger.handlers[1].setLevel(level)
          for c in self.components:
              c.logger.handlers[1].setLevel(level)
  
      def _update_roi(self):
  
          rm = self._roi_data
  
          rm['loglike'] = -self.like()
  
          rm['model_counts'].fill(0)
          try:
              rm['model_counts_wt'].fill(0)
          except:
              rm['model_counts_wt'] = rm['model_counts'].copy()
          rm['npred'] = 0
          rm['npred_wt'] = 0
          for i, c in enumerate(self.components):
              rm['components'][i]['loglike'] = -c.like()
              rm['components'][i]['model_counts'].fill(0)
              try:
                  rm['components'][i]['model_counts_wt'].fill(0)
              except:
                  rm['components'][i]['model_counts_wt'] = rm['components'][i]['model_counts'].copy()
              rm['components'][i]['npred'] = 0
              rm['components'][i]['npred_wt'] = 0
  
          for name in self.like.sourceNames():
  
              # EAC, this is one of the only things we need from setup()
              # self.update_source(name)
  
              src = self.roi.get_source_by_name(name)
              rm['model_counts'] += src['model_counts']
              rm['npred'] += np.sum(src['model_counts'])
              try:
                  rm['model_counts_wt'] += src['model_counts_wt']
                  rm['npred_wt'] += np.sum(src['model_counts_wt'])
              except:
                  rm['model_counts_wt'] += src['model_counts']
                  rm['npred'] += np.sum(src['model_counts'])
  
              mc = self.model_counts_spectrum(name)
              mc_wt = self.model_counts_spectrum(name, weighted=True)
  
              for i, c in enumerate(self.components):
                  rm['components'][i]['model_counts'] += mc[i]
                  rm['components'][i]['model_counts_wt'] += mc_wt[i]
                  rm['components'][i]['npred'] += np.sum(mc[i])
                  rm['components'][i]['npred_wt'] += np.sum(mc_wt[i])
  
      def _update_srcmap(self, name, src, **kwargs):
  
          for c in self.components:
              c._update_srcmap(name, src, **kwargs)
  
          if self._fitcache is not None:
              self._fitcache.update_source(name)
  
      def _create_srcmap_cache(self, name, src):
          for c in self.components:
              c._create_srcmap_cache(name, src)
  
      def _clear_srcmap_cache(self):
          for c in self.components:
              c._srcmap_cache.clear()
  
      def reload_source(self, name, init_source=True):
          """Delete and reload a source in the model.  This will update
          the spatial model of this source to the one defined in the XML
          model."""
  
          for c in self.components:
              c.reload_source(name)
  
          if init_source:
              self._init_source(name)
  
          self.like.model = self.like.components[0].model
  
      def reload_sources(self, names, init_source=True):
  
          for c in self.components:
              c.reload_sources(names)
  
          if init_source:
              for name in names:
                  self._init_source(name)
  
          self.like.model = self.like.components[0].model
  
      def set_source_morphology(self, name, **kwargs):
          """Set the spatial model of a source.
  
          Parameters
          ----------
          name : str
             Source name.
  
          spatial_model : str
             Spatial model name (PointSource, RadialGaussian, etc.).
  
          spatial_pars : dict
             Dictionary of spatial parameters (optional).
  
          use_cache : bool        
             Generate the spatial model by interpolating the cached source
             map.
  
          use_pylike : bool
  
          """
  
          name = self.roi.get_source_by_name(name).name
          src = self.roi[name]
  
          spatial_model = kwargs.get('spatial_model', src['SpatialModel'])
          spatial_pars = kwargs.get('spatial_pars', {})
          use_pylike = kwargs.get('use_pylike', True)
          psf_scale_fn = kwargs.get('psf_scale_fn', None)
          update_source = kwargs.get('update_source', False)
  
          if hasattr(pyLike.BinnedLikelihood, 'setSourceMapImage') and not use_pylike:
              src.set_spatial_model(spatial_model, spatial_pars)
              self._update_srcmap(src.name, src, psf_scale_fn=psf_scale_fn)
          else:
              src = self.delete_source(name, loglevel=logging.DEBUG,
                                       save_template=False)
              src.set_spatial_model(spatial_model, spatial_pars)
              self.add_source(src.name, src, init_source=False,
                              use_pylike=use_pylike, loglevel=logging.DEBUG)
  
          if update_source:
              self.update_source(name)
  
      def set_source_spectrum(self, name, spectrum_type='PowerLaw',
                              spectrum_pars=None, update_source=True):
          """Set the spectral model of a source.  This function can be
          used to change the spectral type of a source or modify its
          spectral parameters.  If called with
          spectrum_type='FileFunction' and spectrum_pars=None, the
          source spectrum will be replaced with a FileFunction with the
          same differential flux distribution as the original spectrum.
  
          Parameters
          ----------
  
          name : str
             Source name.
  
          spectrum_type : str
             Spectrum type (PowerLaw, etc.).
  
          spectrum_pars : dict
             Dictionary of spectral parameters (optional).
  
          update_source : bool
             Recompute all source characteristics (flux, TS, NPred)
             using the new spectral model of the source.
  
          """
          name = self.roi.get_source_by_name(name).name
          src = self.roi[name]
  
          spectrum_pars = {} if spectrum_pars is None else spectrum_pars
          if (self.roi[name]['SpectrumType'] == 'PowerLaw' and
                  spectrum_type == 'LogParabola'):
              spectrum_pars.setdefault('beta', {'value': 0.0, 'scale': 1.0,
                                                'min': 0.0, 'max': 1.0})
              spectrum_pars.setdefault('Eb', src.spectral_pars['Scale'])
              spectrum_pars.setdefault('norm', src.spectral_pars['Prefactor'])
              if 'alpha' not in spectrum_pars:
                  spectrum_pars['alpha'] = src.spectral_pars['Index']
                  spectrum_pars['alpha']['value'] *= -1.0
                  if spectrum_pars['alpha']['scale'] == -1.0:
                      spectrum_pars['alpha']['value'] *= -1.0
                      spectrum_pars['alpha']['scale'] *= -1.0
  
          if spectrum_type == 'FileFunction':
              self._create_filefunction(name, spectrum_pars)
          else:
              fn = gtutils.create_spectrum_from_dict(spectrum_type,
                                                     spectrum_pars)
              self.like.setSpectrum(str(name), fn)
  
          # Get parameters
          src = self.components[0].like.logLike.getSource(str(name))
          pars_dict = gtutils.get_function_pars_dict(src.spectrum())
  
          self.roi[name]['SpectrumType'] = spectrum_type
          self.roi[name].set_spectral_pars(pars_dict)
          for c in self.components:
              c.roi[name]['SpectrumType'] = spectrum_type
              c.roi[name].set_spectral_pars(pars_dict)
  
          if update_source:
              self.update_source(name)
  
      def set_source_dnde(self, name, dnde, update_source=True):
          """Set the differential flux distribution of a source with the
          FileFunction spectral type.
  
          Parameters
          ----------
          name : str
             Source name.
  
          dnde : `~numpy.ndarray`
             Array of differential flux values (cm^{-2} s^{-1} MeV^{-1}).
          """
          name = self.roi.get_source_by_name(name).name
  
          if self.roi[name]['SpectrumType'] != 'FileFunction':
              msg = 'Wrong spectral type: %s' % self.roi[name]['SpectrumType']
              self.logger.error(msg)
              raise Exception(msg)
  
          xy = self.get_source_dnde(name)
  
          if len(dnde) != len(xy[0]):
              msg = 'Wrong length for dnde array: %i' % len(dnde)
              self.logger.error(msg)
              raise Exception(msg)
  
          for c in self.components:
              src = c.like.logLike.getSource(str(name))
              spectrum = src.spectrum()
              file_function = pyLike.FileFunction_cast(spectrum)
              file_function.setSpectrum(10**xy[0], dnde)
  
          if update_source:
              self.update_source(name)
  
      def get_source_dnde(self, name):
          """Return differential flux distribution of a source.  For
          sources with FileFunction spectral type this returns the
          internal differential flux array.
  
          Returns
          -------
          loge : `~numpy.ndarray`
             Array of energies at which the differential flux is
             evaluated (log10(E/MeV)).
  
          dnde : `~numpy.ndarray`
             Array of differential flux values (cm^{-2} s^{-1} MeV^{-1})
             evaluated at energies in ``loge``.
  
          """
          name = self.roi.get_source_by_name(name).name
  
          if self.roi[name]['SpectrumType'] != 'FileFunction':
  
              src = self.components[0].like.logLike.getSource(str(name))
              spectrum = src.spectrum()
              file_function = pyLike.FileFunction_cast(spectrum)
              loge = file_function.log_energy()
              logdnde = file_function.log_dnde()
  
              loge = np.log10(np.exp(loge))
              dnde = np.exp(logdnde)
  
              return loge, dnde
  
          else:
              ebinsz = (self.log_energies[-1] -
                        self.log_energies[0]) / self.enumbins
              loge = utils.extend_array(self.log_energies, ebinsz, 0.5, 6.5)
  
              dnde = np.array([self.like[name].spectrum()(pyLike.dArg(10 ** egy))
                               for egy in loge])
  
              return loge, dnde
  
      def _create_filefunction(self, name, spectrum_pars):
          """Replace the spectrum of an existing source with a
          FileFunction."""
  
          spectrum_pars = {} if spectrum_pars is None else spectrum_pars
  
          if 'loge' in spectrum_pars:
              loge = spectrum_pars.get('loge')
          else:
              ebinsz = (self.log_energies[-1] -
                        self.log_energies[0]) / self.enumbins
              loge = utils.extend_array(self.log_energies, ebinsz, 0.5, 6.5)
  
          # Get the values
          dnde = np.zeros(len(loge))
          if 'dnde' in spectrum_pars:
              dnde = spectrum_pars.get('dnde')
          else:
              dnde = np.array([self.like[name].spectrum()(pyLike.dArg(10 ** egy))
                               for egy in loge])
  
          filename = \
              os.path.join(self.workdir,
                           '%s_filespectrum.txt' % (name.lower().replace(' ', '_')))
  
          # Create file spectrum txt file
          np.savetxt(filename, np.vstack((10**loge, dnde)).T)
          self.like.setSpectrum(name, str('FileFunction'))
  
          self.roi[name]['Spectrum_Filename'] = filename
          # Update
          for c in self.components:
              src = c.like.logLike.getSource(str(name))
              spectrum = src.spectrum()
  
              spectrum.getParam(str('Normalization')).setBounds(1E-3, 1E3)
  
              file_function = pyLike.FileFunction_cast(spectrum)
              file_function.readFunction(str(filename))
              c.roi[name]['Spectrum_Filename'] = filename
  
      def _create_component_configs(self):
          configs = []
  
          components = self.config['components']
  
          common_config = GTBinnedAnalysis.get_config()
          common_config = merge_dict(common_config, self.config,
                                     add_new_keys=True)
  
          if components is None or len(components) == 0:
              cfg = copy.copy(common_config)
              cfg['file_suffix'] = '_00'
              cfg['name'] = '00'
              configs.append(cfg)
          elif isinstance(components, dict):
              for i, k in enumerate(sorted(components.keys())):
                  cfg = copy.copy(common_config)
                  cfg = merge_dict(cfg, components[k], add_new_keys=True)
                  cfg['file_suffix'] = '_' + k
                  cfg['name'] = k
                  configs.append(cfg)
          elif isinstance(components, list):
              for i, c in enumerate(components):
                  cfg = copy.copy(common_config)
                  cfg = merge_dict(cfg, c, add_new_keys=True)
                  cfg['file_suffix'] = '_%02i' % i
                  cfg['name'] = '%02i' % i
                  configs.append(cfg)
          else:
              raise Exception('Invalid type for component block.')
  
          return configs
  
      def _create_component(self, cfg, **kwargs):
  
          self.logger.debug("Creating Analysis Component: " + cfg['name'])
  
          cfg['fileio']['workdir'] = self.config['fileio']['workdir']
  
          for k in cfg.keys():
              if not k in GTBinnedAnalysis.defaults:
                  cfg.pop(k)
  
          comp = GTBinnedAnalysis(
              cfg, self.roi, logging=self.config['logging'], **kwargs)
  
          return comp
  
      def stage_output(self):
          """Copy data products to final output directory."""
  
          if self.workdir == self.outdir:
              return
          elif not os.path.isdir(self.workdir):
              self.logger.error('Working directory does not exist.')
              return
  
          regex = self.config['fileio']['outdir_regex']
          savefits = self.config['fileio']['savefits']
          files = os.listdir(self.workdir)
          self.logger.info('Staging files to %s', self.outdir)
  
          fitsfiles = []
          for c in self.components:
              for f in c.files.values():
  
                  if f is None:
                      continue
  
                  fitsfiles += [os.path.basename(f)]
  
          for f in files:
  
              wpath = os.path.join(self.workdir, f)
              opath = os.path.join(self.outdir, f)
  
              if not utils.match_regex_list(regex, os.path.basename(f)):
                  continue
              if os.path.isfile(opath) and filecmp.cmp(wpath, opath, False):
                  continue
              if not savefits and f in fitsfiles:
                  continue
  
  
              self.logger.debug('Copying ' + f)
              self.logger.info('Copying ' + f)
              shutil.copy(wpath, self.outdir)
  
          self.logger.info('Finished.')
  
      def stage_input(self):
          """Copy input files to working directory."""
  
          if self.workdir == self.outdir:
              return
          elif not os.path.isdir(self.workdir):
              self.logger.error('Working directory does not exist.')
              return
  
          self.logger.info('Staging files to %s', self.workdir)
  
          files = [os.path.join(self.outdir, f)
                   for f in os.listdir(self.outdir)]
  
          regex = copy.deepcopy(self.config['fileio']['workdir_regex'])
  
          for f in files:
  
              if not os.path.isfile(f):
                  continue
              if not utils.match_regex_list(regex, os.path.basename(f)):
                  continue
  
              self.logger.debug('Copying ' + os.path.basename(f))
              shutil.copy(f, self.workdir)
  
          for c in self.components:
              for f in c.files.values():
  
                  if f is None:
                      continue
  
                  wpath = os.path.join(self.workdir, os.path.basename(f))
                  opath = os.path.join(self.outdir, os.path.basename(f))
  
                  if os.path.isfile(wpath):
                      continue
                  elif os.path.isfile(opath):
                      self.logger.debug('Copying ' + os.path.basename(f))
                      shutil.copy(opath, self.workdir)
  
          self.logger.info('Finished.')
  
      def setup(self, init_sources=True, overwrite=False, **kwargs):
          """Run pre-processing for each analysis component and
          construct a joint likelihood object.  This function performs
          the following tasks: data selection (gtselect, gtmktime),
          data binning (gtbin), and model generation (gtexpcube2,gtsrcmaps).
  
          Parameters
          ----------
  
          init_sources : bool
  
             Choose whether to compute properties (flux, TS, etc.) for
             individual sources.
  
          overwrite : bool
  
             Run all pre-processing steps even if the output file of
             that step is present in the working directory.  By default
             this function will skip any steps for which the output file
             already exists.
  
          """
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          self.logger.log(loglevel, 'Running setup.')
  
          # Make spatial maps for extended sources
          for s in self.roi.sources:
              if s.diffuse:
                  continue
              if not s.extended:
                  continue
              self.make_template(s)
  
          # Run setup for each component
          for i, c in enumerate(self.components):
              c.setup(overwrite=overwrite)
  
          # Create likelihood
          self._create_likelihood()
  
          # Determine tmin, tmax
          for i, c in enumerate(self._components):
              self._tmin = (c.tmin if self._tmin is None
                            else min(self._tmin, c.tmin))
              self._tmax = (c.tmax if self._tmax is None
                            else min(self._tmax, c.tmax))
  
          if init_sources:
  
              self.logger.log(loglevel, 'Initializing source properties')
              for name in self.like.sourceNames():
                  self.logger.debug('Initializing source %s', name)
                  self._init_source(name)
              self._update_roi()
  
          self.logger.log(loglevel, 'Finished setup.')
  
      def _create_likelihood(self, srcmdl=None):
          """Instantiate the likelihood object for each component and
          create a SummedLikelihood."""
  
          self._like = SummedLikelihood()
          for c in self.components:
              c._create_binned_analysis(srcmdl)
              self._like.addComponent(c.like)
  
          self.like.model = self.like.components[0].model
          self._fitcache = None
          self._init_roi_model()
  
      def _init_roi_model(self):
  
          rm = self._roi_data
  
          rm['counts'] = np.zeros(self.enumbins)
          rm['counts_wt'] = np.zeros(self.enumbins)
          rm['loglike'] = -self.like()
  
          cmaps = []
          wmaps = []
          proj_type = 0
          for c, crd in zip(self.components, rm['components']):
              cm = c.counts_map()
              wm = c.weight_map()
              cmaps += [cm]
              wmaps += [wm]
              if isinstance(cm, WcsNDMap):
                  crd['counts'] = np.sum(cm.data, axis=(1, 2), keepdims=False)
                  crd['counts_wt'] = np.sum(cm.data * wm.data,
                                            axis=(1, 2), keepdims=False)
              elif isinstance(cm, HpxNDMap):
                  proj_type = 1
                  crd['counts'] = np.sum(cm.data, axis=(1,))
                  crd['counts_wt'] = np.sum(cm.data * wm.data, axis=(1,))
              else:
                  raise ValueError
              crd['loglike'] = -c.like()
  
          if proj_type == 0:
              shape = (self.enumbins, self.npix, self.npix)
          elif proj_type == 1:
              shape = (self.enumbins, np.max(self.geom.npix))
  
          self._coadd_maps(cmaps, shape, rm, wmaps)
  
      def _init_source(self, name):
  
          src = self.roi.get_source_by_name(name)
          src.update_data({'sed': None})
          sd = self.get_src_model(name, paramsonly=True)
          src.update_data(sd)
          for c in self.components:
              src = c.roi.get_source_by_name(name)
              src.update_data(sd)
  
          return src
  
      def cleanup(self):
  
          if self.workdir == self.outdir:
              return
          elif os.path.isdir(self.workdir):
              self.delete_workdir()
  
      def delete_workdir(self):
          self.logger.info('Deleting working directory: ' +
                           self.workdir)
          shutil.rmtree(self.workdir)
  
      def generate_model(self, model_name=None):
          """Generate model maps for all components.  model_name should
          be a unique identifier for the model.  If model_name is None
          then the model maps will be generated using the current
          parameters of the ROI."""
  
          for i, c in enumerate(self._components):
              c.generate_model(model_name=model_name)
  
              # If all model maps have the same spatial/energy binning we
              # could generate a co-added model map here
  
      def set_energy_range(self, logemin, logemax):
          """Set the energy bounds of the analysis.  This restricts the
          evaluation of the likelihood to the data that falls in this
          range.  Input values will be rounded to the closest bin edge
          value.  If either argument is None then the lower or upper
          bound of the analysis instance will be used.
  
          Parameters
          ----------
  
          logemin : float
             Lower energy bound in log10(E/MeV).
  
          logemax : float
             Upper energy bound in log10(E/MeV).
  
          Returns
          -------
  
          eminmax : array
             Minimum and maximum energy in log10(E/MeV).
  
          """
  
          if logemin is None:
              logemin = self.log_energies[0]
          else:
              imin = int(utils.val_to_edge(self.log_energies, logemin)[0])
              logemin = self.log_energies[imin]
  
          if logemax is None:
              logemax = self.log_energies[-1]
          else:
              imax = int(utils.val_to_edge(self.log_energies, logemax)[0])
              logemax = self.log_energies[imax]
  
          self._loge_bounds = np.array([logemin, logemax])
          self._roi_data['loge_bounds'] = np.copy(self.loge_bounds)
          for c in self.components:
              c.set_energy_range(logemin, logemax)
  
          return self._loge_bounds
  
      def counts_map(self):
          """Return a `~fermipy.skymap.Map` representation of the counts map.
  
          Returns
          -------
  
          map : `~fermipy.skymap.Map`
  
          """
          return self._ccube
  
      def weight_map(self):
          """Return a `~fermipy.skymap.Map` representation of the weights map.
  
          Returns
          -------
  
          map : `~fermipy.skymap.Map`
  
          """
          return self._wcube
  
      def model_counts_map(self, name=None, exclude=None, use_mask=False):
          """Return the model counts map for a single source, a list of
          sources, or for the sum of all sources in the ROI.  The
          exclude parameter can be used to exclude one or more
          components when generating the model map.
  
          Parameters
          ----------
          name : str or list of str
  
             Parameter controlling the set of sources for which the
             model counts map will be calculated.  If name=None the
             model map will be generated for all sources in the ROI.
  
          exclude : str or list of str
  
             List of sources that will be excluded when calculating the
             model map.
  
          use_mask : bool
  
             Parameter that specifies in the model counts map should include
             mask pixels (i.e., ones whose weights are <= 0)
  
          Returns
          -------
          map : `~gammapy.maps.Map`
          """        
          maps = [c.model_counts_map(name, exclude, use_mask=use_mask)
                  for c in self.components]
          return skymap.coadd_maps(self.geom, maps)
  
      def model_counts_spectrum(self, name, logemin=None, logemax=None,
                                summed=False, weighted=False):
          """Return the predicted number of model counts versus energy
          for a given source and energy range.  If summed=True return
          the counts spectrum summed over all components otherwise
          return a list of model spectra.   If weighted=True return 
          the weighted version of the counts spectrum
          """
  
          if logemin is None:
              logemin = self.log_energies[0]
          if logemax is None:
              logemax = self.log_energies[-1]
  
          if summed:
              cs = np.zeros(self.enumbins)
              imin = utils.val_to_bin_bounded(self.log_energies,
                                              logemin + 1E-7)[0]
              imax = utils.val_to_bin_bounded(self.log_energies,
                                              logemax - 1E-7)[0] + 1
  
              for c in self.components:
                  ecenter = 0.5 * (c.log_energies[:-1] + c.log_energies[1:])
                  counts = c.model_counts_spectrum(name, self.log_energies[0],
                                                   self.log_energies[-1], weighted)
                  cs += np.histogram(ecenter,
                                     weights=counts,
                                     bins=self.log_energies)[0]
  
              return cs[imin:imax]
          else:
              cs = []
              for c in self.components:
                  cs += [c.model_counts_spectrum(name, logemin,
                                                 logemax, weighted=weighted)]
              return cs
  
      def get_sources(self, cuts=None, distance=None, skydir=None,
                      minmax_ts=None, minmax_npred=None, exclude=None,
                      square=False):
          """Retrieve list of sources in the ROI satisfying the given
          selections.
  
          Returns
          -------
          srcs : list
              A list of `~fermipy.roi_model.Model` objects.
  
          """
  
          coordsys = self.config['binning']['coordsys']
          return self.roi.get_sources(skydir, distance, cuts,
                                      minmax_ts, minmax_npred,
                                      exclude, square,
                                      coordsys=coordsys)
  
      def add_source(self, name, src_dict, free=None, init_source=True,
                     save_source_maps=True, use_pylike=True,
                     use_single_psf=False, **kwargs):
          """Add a source to the ROI model.  This function may be called
          either before or after `~fermipy.gtanalysis.GTAnalysis.setup`.
  
          Parameters
          ----------
          name : str
              Source name.
  
          src_dict : dict or `~fermipy.roi_model.Source` object
              Dictionary or source object defining the source properties
              (coordinates, spectral parameters, etc.).
  
          free : bool
              Initialize the source with a free normalization parameter.
  
          use_pylike : bool
              Create source maps with pyLikelihood.
  
          use_single_psf : bool 
              Use the PSF model calculated for the ROI center.  If false
              then a new model will be generated using the position of
              the source.
  
          """
  
          if self.roi.has_source(name):
              msg = 'Source %s already exists.' % name
              self.logger.error(msg)
              raise Exception(msg)
  
          loglevel = kwargs.pop('loglevel', self.loglevel)
  
          self.logger.log(loglevel, 'Adding source ' + name)
  
          src = self.roi.create_source(name, src_dict, rescale=True)
          self.make_template(src)
  
          for c in self.components:
              c.add_source(name, src_dict, free=free,
                           save_source_maps=save_source_maps,
                           use_pylike=use_pylike,
                           use_single_psf=use_single_psf)
  
          if self._like is None:
              return
  
          use_edisp = self.config['gtlike']['edisp'] and self.config['gtlike']['edisp_bins'] != 0
          if use_edisp and src.name not in \
                  self.config['gtlike']['edisp_disable']:
              self.set_edisp_flag(src.name, True)
  
          self.like.syncSrcParams(str(name))
          self.like.model = self.like.components[0].model
          # if free is not None:
          #    self.free_norm(name, free, loglevel=logging.DEBUG)
  
          if init_source:
              self._init_source(name)
              self._update_roi()
  
          if self._fitcache is not None:
              self._fitcache.update_source(name)
  
      def add_sources_from_roi(self, names, roi, free=False, **kwargs):
          """Add multiple sources to the current ROI model copied from another ROI model.
  
          Parameters
          ----------
  
          names : list
              List of str source names to add.
  
          roi : `~fermipy.roi_model.ROIModel` object
              The roi model from which to add sources.
  
          free : bool
              Initialize the source with a free normalization paramter.
  
          """
  
          for name in names:
              self.add_source(name, roi[name].data, free=free, **kwargs)
  
      def delete_source(self, name, save_template=True, delete_source_map=False,
                        build_fixed_wts=True, **kwargs):
          """Delete a source from the ROI model.
  
          Parameters
          ----------
          name : str
              Source name.
  
          save_template : bool
              Keep the SpatialMap FITS template associated with this
              source.
  
          delete_source_map : bool
              Delete the source map associated with this source from the
              source maps file.
  
          Returns
          -------
          src : `~fermipy.roi_model.Model`
              The deleted source object.
  
          """
  
          if not self.roi.has_source(name):
              self.logger.error('No source with name: %s', name)
              return
  
          loglevel = kwargs.pop('loglevel', self.loglevel)
  
          self.logger.log(loglevel, 'Deleting source %s', name)
  
          # STs require a source to be freed before deletion
          if self.like is not None:
              self.free_norm(name, loglevel=logging.DEBUG)
  
          for c in self.components:
              c.delete_source(name, save_template=save_template,
                              delete_source_map=delete_source_map,
                              build_fixed_wts=build_fixed_wts)
  
          src = self.roi.get_source_by_name(name)
          self.roi.delete_sources([src])
          if self.like is not None:
              self.like.model = self.like.components[0].model
              self._update_roi()
          return src
  
      def delete_sources(self, cuts=None, distance=None,
                         skydir=None, minmax_ts=None, minmax_npred=None,
                         exclude=None, square=False, names=None):
          """Delete sources in the ROI model satisfying the given
          selection criteria.
  
          Parameters
          ----------
          cuts : dict
              Dictionary of [min,max] selections on source properties.
  
          distance : float
              Cut on angular distance from ``skydir``.  If None then no
              selection will be applied.
  
          skydir : `~astropy.coordinates.SkyCoord`
              Reference sky coordinate for ``distance`` selection.  If
              None then the distance selection will be applied with
              respect to the ROI center.
  
          minmax_ts : list
              Select sources that have TS in the range [min,max].  If
              either min or max are None then only a lower (upper) bound
              will be applied.  If this parameter is none no selection
              will be applied.
  
          minmax_npred : list
              Select sources that have npred in the range [min,max].  If
              either min or max are None then only a lower (upper) bound
              will be applied.  If this parameter is none no selection
              will be applied.
  
          square : bool
              Switch between applying a circular or square (ROI-like)
              selection on the maximum projected distance from the ROI
              center.
  
          names : list            
              Select sources matching a name in this list.
  
          Returns
          -------
          srcs : list
              A list of `~fermipy.roi_model.Model` objects.
  
          """
  
          srcs = self.roi.get_sources(skydir=skydir, distance=distance, cuts=cuts,
                                      minmax_ts=minmax_ts, minmax_npred=minmax_npred,
                                      exclude=exclude, square=square,
                                      coordsys=self.config[
                                          'binning']['coordsys'],
                                      names=names)
  
          for s in srcs:
              self.delete_source(s.name, build_fixed_wts=False)
  
          if self.like is not None:
              # Build fixed model weights in one pass
              for c in self.components:
                  c.like.logLike.buildFixedModelWts()
  
              self._update_roi()
  
          return srcs
  
      def set_weights_map(self, wmap, update_roi=True):
  
          if self.like is None:
              raise RuntimeError(
                  "GTAnalysis.set_weights_map must be run after GTAnalysis.setup")
  
          likes = []
          if isinstance(self.like, gtutils.SummedLikelihood):
              for c in self.like.components:
                  likes += [c.logLike]
          else:
              likes += [self.like.logLike]
          projmaps = []
          if isinstance(wmap, list):
              for w in wmap:
                  if w is None:
                      projmap = None
                  else:
                      projmap = pyLike.WcsMapLibrary.instance().wcsmap(w, "")
                  projmaps += [projmap]
              if len(projmaps) != len(likes):
                  raise ValueError("Length of weights maps list is not equal to number of components %i %i" % (
                      len(projmaps), len(likes)))
              for l, projmap in zip(likes, projmaps):
                  l.setWeightsMap(projmap)
          else:
              if wmap is None:
                  projmap = None
              else:
                  projmap = pyLike.WcsMapLibrary.instance().wcsmap(wmap, "")
              for l in likes:
                  l.setWeightsMap(projmap)
  
          if update_roi is False:
              return
          self._init_roi_model()
          for name in self.like.sourceNames():
              self.logger.debug('Initializing source %s', name)
              self._init_source(name)
          self._update_roi()
  
      def make_template(self, src):
  
          if src['SpatialType'] != 'SpatialMap':
              return
  
          if src['SpatialModel'] in ['RadialGaussian']:
  
              template_file = '%s_template_gauss_%05.3f.fits' % (src.name,
                                                                 src['SpatialWidth'])
              template_file = os.path.join(self.workdir, template_file)
  
              sigma = src['SpatialWidth'] / 1.5095921854516636
              srcmap_utils.make_gaussian_spatial_map(src.skydir, sigma,
                                                     template_file)
              src['Spatial_Filename'] = template_file
          elif src['SpatialModel'] in ['RadialDisk']:
  
              template_file = '%s_template_disk_%05.3f.fits' % (src.name,
                                                                src['SpatialWidth'])
              template_file = os.path.join(self.workdir, template_file)
  
              radius = src['SpatialWidth'] / 0.8246211251235321
              srcmap_utils.make_disk_spatial_map(src.skydir, radius,
                                                 template_file)
              src['Spatial_Filename'] = template_file
  
      def free_sources_by_name(self, names, free=True, pars=None,
                               **kwargs):
          """Free all sources with names matching ``names``.
  
          Parameters
          ----------
          names : list
              List of source names.
  
          free : bool
              Choose whether to free (free=True) or fix (free=False)
              source parameters.
  
          pars : list
              Set a list of parameters to be freed/fixed for each
              source.  If none then all source parameters will be
              freed/fixed.  If pars='norm' then only normalization
              parameters will be freed.
  
          Returns
          -------
          srcs : list
              A list of `~fermipy.roi_model.Model` objects.            
          """
  
          if names is None:
              return
  
          names = [names] if not isinstance(names, list) else names
          names = [self.roi.get_source_by_name(t).name for t in names]
          srcs = [s for s in self.roi.sources if s.name in names]
          for s in srcs:
              self.free_source(s.name, free=free, pars=pars, **kwargs)
  
          return srcs
  
      def free_sources(self, free=True, pars=None, cuts=None,
                       distance=None, skydir=None, minmax_ts=None, minmax_npred=None,
                       exclude=None, square=False, **kwargs):
          """Free or fix sources in the ROI model satisfying the given
          selection.  When multiple selections are defined, the selected
          sources will be those satisfying the logical AND of all
          selections (e.g. distance < X && minmax_ts[0] < ts <
          minmax_ts[1] && ...).
  
          Parameters
          ----------
  
          free : bool
              Choose whether to free (free=True) or fix (free=False)
              source parameters.
  
          pars : list
              Set a list of parameters to be freed/fixed for each
              source.  If none then all source parameters will be
              freed/fixed.  If pars='norm' then only normalization
              parameters will be freed.
  
          cuts : dict
              Dictionary of [min,max] selections on source properties.
  
          distance : float
              Cut on angular distance from ``skydir``.  If None then no
              selection will be applied.
  
          skydir : `~astropy.coordinates.SkyCoord`
              Reference sky coordinate for ``distance`` selection.  If
              None then the distance selection will be applied with
              respect to the ROI center.
  
          minmax_ts : list
              Free sources that have TS in the range [min,max].  If
              either min or max are None then only a lower (upper) bound
              will be applied.  If this parameter is none no selection
              will be applied.
  
          minmax_npred : list
              Free sources that have npred in the range [min,max].  If
              either min or max are None then only a lower (upper) bound
              will be applied.  If this parameter is none no selection
              will be applied.
  
          exclude : list
              Names of sources that will be excluded from the selection.
  
          square : bool
              Switch between applying a circular or square (ROI-like)
              selection on the maximum projected distance from the ROI
              center.
  
          Returns
          -------
          srcs : list
              A list of `~fermipy.roi_model.Model` objects.
  
          """
  
          srcs = self.roi.get_sources(skydir=skydir, distance=distance,
                                      cuts=cuts, minmax_ts=minmax_ts,
                                      minmax_npred=minmax_npred, exclude=exclude,
                                      square=square,
                                      coordsys=self.config['binning']['coordsys'])
  
          for s in srcs:
              self.free_source(s.name, free=free, pars=pars, **kwargs)
  
          return srcs
  
      def set_edisp_flag(self, name, flag=True):
          """Enable or disable the energy dispersion correction for the
          given source."""
          src = self.roi.get_source_by_name(name)
          name = src.name
  
          for c in self.components:
              c.like[name].src.set_edisp_flag(flag)
  
      def scale_parameter(self, name, par, scale):
  
          idx = self.like.par_index(name, par)
          self.like[idx].setScale(self.like[idx].getScale() * scale)
          self._sync_params(name)
  
      def _set_value_bounded(self, idx, value):
          bounds = list(self.like.model[idx].getBounds())
          value = max(value, bounds[0])
          value = min(value, bounds[1])
          self.like[idx].setValue(value)
  
      def set_parameter(self, name, par, value, true_value=True, scale=None,
                        bounds=None, error=None, update_source=True):
          """
          Update the value of a parameter.  Parameter bounds will
          automatically be adjusted to encompass the new parameter
          value.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          par : str
              Parameter name.
  
          value : float
              Parameter value.  By default this argument should be the
              unscaled (True) parameter value.
  
          scale : float
              Parameter scale (optional).  Value argument is interpreted
              with respect to the scale parameter if it is provided.
  
          error : float
              Parameter error (optional).  By default this argument should be the
              unscaled (True) parameter value.
  
          update_source : bool
              Update the source dictionary for the object.
  
          """
          name = self.roi.get_source_by_name(name).name
          idx = self.like.par_index(name, par)
          current_bounds = list(self.like.model[idx].getBounds())
  
          if scale is not None:
              self.like[idx].setScale(scale)
          else:
              scale = self.like.model[idx].getScale()
  
          if true_value:
              current_bounds[0] = min(current_bounds[0], value / scale)
              current_bounds[1] = max(current_bounds[1], value / scale)
              if error is not None:
                  error = error / scale
          else:
              current_bounds[0] = min(current_bounds[0], value)
              current_bounds[1] = max(current_bounds[1], value)
  
          # update current bounds to encompass new value
          self.like[idx].setBounds(*current_bounds)
  
          if true_value:
              for p in self.like[idx].pars:
                  p.setTrueValue(value)
          else:
              self.like[idx].setValue(value)
  
          if bounds is not None:
              if true_value:
                  bounds[0] = min(bounds[0], value / scale)
                  bounds[1] = max(bounds[1], value / scale)
              else:
                  bounds[0] = min(bounds[0], value)
                  bounds[1] = max(bounds[1], value)
  
          # For some reason the numerical accuracy is causing this to throw exceptions.
          try:
              if bounds is not None:
                  self.like[idx].setBounds(*bounds)
          except RuntimeError:
              self.logger.warning(
                  "Caught failure on setBounds for %s::%s." % (name, par))
              pass
  
          if error is not None:
              self.like[idx].setError(error)
  
          self._sync_params(name)
  
          if update_source:
              self.update_source(name)
  
      def set_parameter_scale(self, name, par, scale):
          """Update the scale of a parameter while keeping its value constant."""
          name = self.roi.get_source_by_name(name).name
          idx = self.like.par_index(name, par)
          current_bounds = list(self.like.model[idx].getBounds())
          current_scale = self.like.model[idx].getScale()
          current_value = self.like[idx].getValue()
  
          self.like[idx].setScale(scale)
          self.like[idx].setValue(current_value * current_scale / scale)
          self.like[idx].setBounds(current_bounds[0] * current_scale / scale,
                                   current_bounds[1] * current_scale / scale)
          self._sync_params(name)
  
      def set_parameter_bounds(self, name, par, bounds):
          """Set the bounds on the scaled value of a parameter.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          par : str
              Parameter name.
  
          bounds : list
              Upper and lower bound.
  
          """
          idx = self.like.par_index(name, par)
          self.like[idx].setBounds(*bounds)
          self._sync_params(name)
  
      def set_parameter_error(self, name, par, error):
          """Set the error on the value of a parameter.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          par : str
              Parameter name.
  
          error : float
              The value for the parameter error
          """
          idx = self.like.par_index(name, par)
          self.like[idx].setError(error)
          self._sync_params(name)
  
      def lock_parameter(self, name, par, lock=True):
          """Set parameter to locked/unlocked state.  A locked parameter
          will be ignored when running methods that free/fix sources or
          parameters.
  
          Parameters
          ----------
          name : str
              Source name.
  
          par : str
              Parameter name.
  
          lock : bool       
              Set parameter to locked (True) or unlocked (False) state.
          """
          name = self.roi.get_source_by_name(name).name
          lck_params = self._lck_params.setdefault(name, [])
  
          if lock:
              self.free_parameter(name, par, False)
              if not par in lck_params:
                  lck_params += [par]
          else:
              if par in lck_params:
                  lck_params.remove(par)
  
      def free_parameter(self, name, par, free=True):
          """Free/Fix a parameter of a source by name.
  
          Parameters
          ----------
          name : str
              Source name.
  
          par : str
              Parameter name.
          """
          name = self.get_source_name(name)
          if par in self._lck_params.get(name, []):
              return
          idx = self.like.par_index(name, par)
          self.like[idx].setFree(free)
          self._sync_params(name)
  
      def lock_source(self, name, lock=True):
          """Set all parameters of a source to a locked/unlocked state.
          Locked parameters will be ignored when running methods that
          free/fix sources or parameters.
  
          Parameters
          ----------
          name : str
              Source name.
  
          lock : bool        
              Set source parameters to locked (True) or unlocked (False)
              state.
          """
          name = self.get_source_name(name)
          if lock:
              par_names = self.get_source_params(name)
              self.free_source(name, False, pars=par_names)
              self._lck_params[name] = par_names
          else:
              self._lck_params[name] = []
  
      def free_source(self, name, free=True, pars=None, **kwargs):
          """Free/Fix parameters of a source.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          free : bool
              Choose whether to free (free=True) or fix (free=False)
              source parameters.
  
          pars : list
              Set a list of parameters to be freed/fixed for this source.  If
              none then all source parameters will be freed/fixed with the
              exception of those defined in the skip_pars list.
  
          """
  
          free_pars = self.get_free_param_vector()
  
          loglevel = kwargs.pop('loglevel', self.loglevel)
  
          # Find the source
          src = self.roi.get_source_by_name(name)
          name = src.name
  
          if pars is None or (isinstance(pars, list) and not pars):
              pars = []
              pars += norm_parameters.get(src['SpectrumType'], [])
              pars += shape_parameters.get(src['SpectrumType'], [])
          elif pars == 'norm':
              pars = []
              pars += norm_parameters.get(src['SpectrumType'], [])
          elif pars == 'shape':
              pars = []
              pars += shape_parameters.get(src['SpectrumType'], [])
          elif isinstance(pars, list):
              pass
          else:
              raise Exception('Invalid parameter list.')
  
          # Remove locked parameters
          lck_params = self._lck_params.get(name, [])
          pars = [p for p in pars if p not in lck_params]
  
          # Deduce here the names of all parameters from the spectral type
          src_par_names = pyLike.StringVector()
          self.like[name].src.spectrum().getParamNames(src_par_names)
  
          par_indices = []
          par_names = []
          for p in src_par_names:
              if pars is not None and p not in pars:
                  continue
  
              idx = self.like.par_index(name, p)
              if free == free_pars[idx]:
                  continue
  
              par_indices.append(idx)
              par_names.append(p)
  
          if len(par_names) == 0:
              return
  
          if free:
              self.logger.log(loglevel, 'Freeing parameters for %-22s: %s',
                              name, par_names)
          else:
              self.logger.log(loglevel, 'Fixing parameters for %-22s: %s',
                              name, par_names)
  
          for (idx, par_name) in zip(par_indices, par_names):
              self.like[idx].setFree(free)
          self._sync_params_state(name)
  
      def set_norm_scale(self, name, value):
          name = self.get_source_name(name)
          normPar = self.like.normPar(name)
          normPar.setScale(value)
          self._sync_params(name)
  
      def set_norm(self, name, value, update_source=True):
          name = self.get_source_name(name)
          par = self.like.normPar(name).getName()
          self.set_parameter(name, par, value, true_value=False,
                             update_source=update_source)
  
      def set_norm_bounds(self, name, bounds):
          name = self.get_source_name(name)
          par = self.like.normPar(name).getName()
          self.set_parameter_bounds(name, par, bounds)
  
      def free_norm(self, name, free=True, **kwargs):
          """Free/Fix normalization of a source.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          free : bool
              Choose whether to free (free=True) or fix (free=False).
  
          """
          name = self.get_source_name(name)
          normPar = self.like.normPar(name).getName()
          self.free_source(name, pars=[normPar], free=free, **kwargs)
  
      def free_index(self, name, free=True, **kwargs):
          """Free/Fix index of a source.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          free : bool
              Choose whether to free (free=True) or fix (free=False).
  
          """
          src = self.roi.get_source_by_name(name)
          self.free_source(name, free=free,
                           pars=index_parameters.get(src['SpectrumType'], []),
                           **kwargs)
  
      def free_shape(self, name, free=True, **kwargs):
          """Free/Fix shape parameters of a source.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          free : bool
              Choose whether to free (free=True) or fix (free=False).
          """
          src = self.roi.get_source_by_name(name)
          self.free_source(name, free=free,
                           pars=shape_parameters[src['SpectrumType']],
                           **kwargs)
  
      def _sync_params(self, name):
          self.like.syncSrcParams(str(name))
          src = self.components[0].like.logLike.getSource(str(name))
          spectral_pars = gtutils.get_function_pars_dict(src.spectrum())
          self.roi[name].set_spectral_pars(spectral_pars)
  
      def _sync_params_state(self, name=None):
  
          if name is None:
              self.like.syncSrcParams()
              names = self.like.sourceNames()
          else:
              self.like.syncSrcParams(str(name))
              names = [name]
  
          for name in names:
  
              src = self.like[name].src
              pars = gtutils.get_function_pars(src.spectrum())
              for p in pars:
                  self.roi[name].spectral_pars[p['name']]['free'] = p['free']
  
      def get_norm(self, name):
          name = self.get_source_name(name)
          return self.like.normPar(name).getValue()
  
      def _get_param(self, name, par_name):
  
          idx = self.like.par_index(name, par_name)
          par = self.like.model[idx]
          bounds = par.getBounds()
          return {'value': par.getValue(),
                  'error': par.error(),
                  'scale': par.getScale(),
                  'free': par.isFree(),
                  'min': bounds[0], 'max': bounds[1]}
  
      def get_params(self, freeonly=False):
  
          params = {}
          for srcName in self.like.sourceNames():
  
              par_names = pyLike.StringVector()
              src = self.components[0].like.logLike.getSource(str(srcName))
              src.spectrum().getParamNames(par_names)
  
  #            for parName in self.get_free_source_params(srcName):
              for parName in par_names:
                  idx = self.like.par_index(srcName, parName)
                  par = self.like.model[idx]
  
                  if freeonly and not par.isFree():
                      continue
  
                  is_norm = parName == self.like.normPar(srcName).getName()
                  params[idx] = {'src_name': srcName,
                                 'par_name': parName,
                                 'is_norm': is_norm,
                                 'idx': idx}
                  params[idx].update(self._get_param(srcName, parName))
  
          return [params[k] for k in sorted(params.keys())]
  
      def get_free_param_vector(self):
          free = []
          for p in self.like.params():
              free.append(p.isFree())
          return free
  
      def set_free_param_vector(self, free):
          for i, t in enumerate(free):
              if t:
                  self.like.thaw(i)
              else:
                  self.like.freeze(i)
  
          self._sync_params_state()
  
      def _latch_free_params(self):
          self._free_params = self.get_free_param_vector()
  
      def _restore_free_params(self):
          self.set_free_param_vector(self._free_params)
  
      def _latch_state(self):
          self._saved_state = LikelihoodState(self.like)
          return self._saved_state
  
      def _restore_state(self):
          if self._saved_state is None:
              return
          self._saved_state.restore()
  
      def get_source_params(self, name):
          name = self.get_source_name(name)
          spectrum = self.like[name].src.spectrum()
          parNames = pyLike.StringVector()
          spectrum.getParamNames(parNames)
          return [str(p) for p in parNames]
  
      def get_free_source_params(self, name):
          name = self.get_source_name(name)
          spectrum = self.like[name].src.spectrum()
          parNames = pyLike.StringVector()
          spectrum.getFreeParamNames(parNames)
          return [str(p) for p in parNames]
  
      def get_source_name(self, name):
          """Return the name of a source as it is defined in the
          pyLikelihood model object."""
          if name not in self.like.sourceNames():
              name = self.roi.get_source_by_name(name).name
          return name
  
      def zero_source(self, name, **kwargs):
          normPar = self.like.normPar(name).getName()
          self.scale_parameter(name, normPar, 1E-10)
          self.free_source(name, free=False, **kwargs)
  
      def unzero_source(self, name, **kwargs):
          normPar = self.like.normPar(name).getName()
          self.scale_parameter(name, normPar, 1E10)
  
      def optimize(self, **kwargs):
          """Iteratively optimize the ROI model.  The optimization is
          performed in three sequential steps:
  
          * Free the normalization of the N largest components (as
            determined from NPred) that contain a fraction ``npred_frac``
            of the total predicted counts in the model and perform a
            simultaneous fit of the normalization parameters of these
            components.
  
          * Individually fit the normalizations of all sources that were
            not included in the first step in order of their npred
            values.  Skip any sources that have NPred <
            ``npred_threshold``.
  
          * Individually fit the shape and normalization parameters of
            all sources with TS > ``shape_ts_threshold`` where TS is
            determined from the first two steps of the ROI optimization.
  
          To ensure that the model is fully optimized this method can be
          run multiple times.
  
          Parameters
          ----------
  
          npred_frac : float
              Threshold on the fractional number of counts in the N
              largest components in the ROI.  This parameter determines
              the set of sources that are fit in the first optimization
              step.
  
          npred_threshold : float
              Threshold on the minimum number of counts of individual
              sources.  This parameter determines the sources that are
              fit in the second optimization step.
  
          shape_ts_threshold : float
              Threshold on source TS used for determining the sources
              that will be fit in the third optimization step.
  
          max_free_sources : int
              Maximum number of sources that will be fit simultaneously
              in the first optimization step.
  
          skip : list
              List of str source names to skip while optimizing.
  
          optimizer : dict
              Dictionary that overrides the default optimizer settings.
  
          """
  
          loglevel = kwargs.pop('loglevel', self.loglevel)
          timer = Timer.create(start=True)
          self.logger.log(loglevel, 'Starting')
  
          loglike0 = -self.like()
          self.logger.debug('LogLike: %f' % loglike0)
  
          # Extract options from kwargs
          config = copy.deepcopy(self.config['roiopt'])
          config['optimizer'] = copy.deepcopy(self.config['optimizer'])
          fermipy.config.validate_config(kwargs, config)
          config = merge_dict(config, kwargs)
  
          # Extract options from kwargs
          npred_frac_threshold = config['npred_frac']
          npred_threshold = config['npred_threshold']
          shape_ts_threshold = config['shape_ts_threshold']
          max_free_sources = config['max_free_sources']
          skip = copy.deepcopy(config['skip'])
  
          o = defaults.make_default_dict(defaults.roiopt_output)
          o['config'] = config
          o['loglike0'] = loglike0
  
          # preserve free parameters
          free = self.get_free_param_vector()
  
          # Fix all parameters
          self.free_sources(free=False, loglevel=logging.DEBUG)
  
          # Free norms of sources for which the sum of npred is a
          # fraction > npred_frac of the total model counts in the ROI
          npred_sum = 0
          skip_sources = skip if skip != None else []
          joint_norm_fit = []
  
          # EAC, we need this try block by older version of the ST don't have has_weights function
          try:
              if self.like.logLike.has_weights():
                  npred_str = 'npred_wt'
              else:
                  npred_str = 'npred'
          except AttributeError:
              npred_str = 'npred'
  
          # FIXME, EAC, use npred_wt here
          for s in sorted(self.roi.sources, key=lambda t: t[npred_str],
                          reverse=True):
  
              npred_sum += s[npred_str]
              npred_frac = npred_sum / self._roi_data[npred_str]
  
              if s.name in skip_sources:
                  continue
  
              self.free_norm(s.name, loglevel=logging.DEBUG)
              joint_norm_fit.append(s.name)
  
              if npred_frac > npred_frac_threshold:
                  break
              if s[npred_str] < npred_threshold:
                  break
              if len(joint_norm_fit) >= max_free_sources:
                  break
  
          print ("Joint fit ", joint_norm_fit)
          self.fit(loglevel=logging.DEBUG, **config['optimizer'])
          self.free_sources(free=False, loglevel=logging.DEBUG)
  
          # Step through remaining sources and re-fit normalizations
          # FIXME, EAC, use npred_wt here
          for s in sorted(self.roi.sources, key=lambda t: t[npred_str],
                          reverse=True):
  
              if s.name in skip_sources or s.name in joint_norm_fit:
                  continue
  
              if s[npred_str] < npred_threshold:
                  self.logger.debug(
                      'Skipping %s with npred %10.3f', s.name, s[npred_str])
                  continue
  
              self.logger.debug('Fitting %s npred: %10.3f TS: %10.3f',
                                s.name, s[npred_str], s['ts'])
              self.free_norm(s.name, loglevel=logging.DEBUG)
              self.fit(loglevel=logging.DEBUG, **config['optimizer'])
              self.logger.debug('Post-fit Results npred: %10.3f TS: %10.3f',
                                s[npred_str], s['ts'])
              self.free_norm(s.name, free=False, loglevel=logging.DEBUG)
  
          # Refit spectral shape parameters for sources with TS >
          # shape_ts_threshold
          for s in sorted(self.roi.sources,
                          key=lambda t: t['ts'] if np.isfinite(t['ts']) else 0,
                          reverse=True):
  
              if s.name in skip_sources:
                  continue
  
              if s['ts'] < shape_ts_threshold or not np.isfinite(s['ts']):
                  continue
  
              print ('Fitting shape %s TS: %10.3f' % (s.name, s['ts']))
              self.logger.debug('Fitting shape %s TS: %10.3f', s.name, s['ts'])
              self.free_source(s.name, loglevel=logging.DEBUG)
              self.fit(loglevel=logging.DEBUG, **config['optimizer'])
              self.free_source(s.name, free=False, loglevel=logging.DEBUG)
  
          self.set_free_param_vector(free)
  
          loglike1 = -self.like()
  
          o['loglike1'] = loglike1
          o['dloglike'] = loglike1 - loglike0
  
          self.logger.log(loglevel, 'Finished')
          self.logger.log(loglevel, 'LogLike: %f Delta-LogLike: %f',
                          loglike1, loglike1 - loglike0)
          self.logger.log(loglevel, 'Execution time: %.2f s', timer.elapsed_time)
          return o
  
      def profile_norm(self, name, logemin=None, logemax=None, reoptimize=False,
                       xvals=None, npts=None, fix_shape=True, savestate=True,
                       **kwargs):
          """Profile the normalization of a source.
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          reoptimize : bool
              Re-optimize free parameters in the model at each point
              in the profile likelihood scan.
  
          """
  
          self.logger.debug('Profiling %s', name)
  
          if savestate:
              saved_state = LikelihoodState(self.like)
  
          if fix_shape:
              self.free_sources(False, pars='shape', loglevel=logging.DEBUG)
  
          if npts is None:
              npts = self.config['gtlike']['llscan_npts']
  
          # Find the source
          name = self.roi.get_source_by_name(name).name
          parName = self.like.normPar(name).getName()
  
          loge_bounds = self.loge_bounds
          if logemin is not None or logemax is not None:
              self.set_energy_range(logemin, logemax)
  
          # Find a sequence of values for the normalization scan
          if xvals is None:
              if reoptimize:
                  xvals = self._find_scan_pts_reopt(name, npts=npts,
                                                    **kwargs)
              else:
                  xvals = self._find_scan_pts(name, npts=9)
                  lnlp = self.profile(name, parName,
                                      reoptimize=False, xvals=xvals)
                  lims = utils.get_parameter_limits(lnlp['xvals'],
                                                    lnlp['dloglike'],
                                                    cl_limit=0.99)
  
                  if not np.isfinite(lims['ul']):
                      self.logger.warning('Upper limit not found.  '
                                          'Refitting normalization.')
                      self.like.optimize(0)
                      xvals = self._find_scan_pts(name, npts=npts)
                      lnlp = self.profile(name, parName,
                                          reoptimize=False,
                                          xvals=xvals)
                      lims = utils.get_parameter_limits(lnlp['xvals'],
                                                        lnlp['dloglike'],
                                                        cl_limit=0.99)
  
                  if np.isfinite(lims['ll']):
                      xhi = np.linspace(lims['x0'], lims['ul'], npts - npts // 2)
                      xlo = np.linspace(lims['ll'], lims['x0'], npts // 2)
                      xvals = np.concatenate((xlo[:-1], xhi))
                      xvals = np.insert(xvals, 0, 0.0)
                  elif np.abs(lnlp['dloglike'][0] - lims['lnlmax']) > 0.1:
                      lims['ll'] = 0.0
                      xhi = np.linspace(lims['x0'], lims['ul'],
                                        (npts + 1) - (npts + 1) // 2)
                      xlo = np.linspace(lims['ll'], lims['x0'], (npts + 1) // 2)
                      xvals = np.concatenate((xlo[:-1], xhi))
                  else:
                      xvals = np.linspace(0, lims['ul'], npts)
  
          o = self.profile(name, parName,
                           reoptimize=reoptimize, xvals=xvals,
                           savestate=savestate, **kwargs)
  
          if savestate:
              saved_state.restore()
  
          if logemin is not None or logemax is not None:
              self.set_energy_range(*loge_bounds)
  
          self.logger.debug('Finished')
  
          return o
  
      def _find_scan_pts(self, name, logemin=None, logemax=None, npts=20):
  
          par = self.like.normPar(name)
  
          loge_bounds = [self.loge_bounds[0] if logemin is None else logemin,
                         self.loge_bounds[1] if logemax is None else logemax]
  
          val = par.getValue()
  
          if val == 0:
              par.setValue(1.0)
              self.like.syncSrcParams(str(name))
              cs = self.model_counts_spectrum(name,
                                              loge_bounds[0],
                                              loge_bounds[1],
                                              summed=True, weighted=True)
              npred = np.sum(cs)
              # EAC, protect against npred == 0
              val = 1. / max(npred, 1.)
              npred = 1.0
              par.setValue(0.0)
              self.like.syncSrcParams(str(name))
          else:
              cs = self.model_counts_spectrum(name,
                                              loge_bounds[0],
                                              loge_bounds[1],
                                              summed=True, weighted=True)
              npred = np.sum(cs)
  
          if npred < 10:
              val *= 1. / max(1e-16, npred)
              xvals = val * 10 ** np.linspace(-1.0, 3.0, npts - 1)
              xvals = np.insert(xvals, 0, 0.0)
          else:
  
              npts_lo = npts // 2
              npts_hi = npts - npts_lo
              xhi = np.linspace(0, 1, npts_hi)
              xlo = np.linspace(-1, 0, npts_lo)
              xvals = val * 10 ** np.concatenate((xlo[:-1], xhi))
              xvals = np.insert(xvals, 0, 0.0)
  
          return xvals
  
      def _find_scan_pts_reopt(self, name, logemin=None, logemax=None, npts=20,
                               dloglike_thresh=2.7059, **kwargs):
  
          parName = self.like.normPar(name).getName()
  
          cl_limit = utils.onesided_dlnl_to_cl(dloglike_thresh)
  
          npts = max(npts, 5)
          xvals = self._find_scan_pts(name, logemin=logemin, logemax=logemax,
                                      npts=20)
  
          # Generate likelihood profile w/ fixed nuisance pars
          lnlp0 = self.profile(name, parName, logemin=logemin, logemax=logemax,
                               reoptimize=False, xvals=xvals, **kwargs)
          xval0 = self.like.normPar(name).getValue()
          lims0 = utils.get_parameter_limits(lnlp0['xvals'], lnlp0['dloglike'],
                                             cl_limit=cl_limit)
  
          # Difference in lnL at maximum with respect to lnL(x=0)
          dlnlmax = np.abs(lnlp0['dloglike'][0] - lims0['lnlmax'])
  
          if dlnlmax < 0.05:
              xvals = np.array([0.0, lims0['x0'] + lims0['err_hi'], lims0['ul']])
          else:
              xvals = np.array([lims0['ll'],
                                lims0['x0'] - lims0['err_lo'], lims0['x0'],
                                lims0['x0'] + lims0['err_hi'], lims0['ul']])
              xvals = xvals[np.isfinite(xvals)]
  
          # Generate likelihood profile w/ free nuisance pars
          if np.isnan(xvals).any():
              raise ValueError(
                  "Parameter scan points for %s::%s include infinite value." % (name, parName))
  
          lnlp1 = self.profile(name, parName, logemin=logemin, logemax=logemax,
                               reoptimize=True, xvals=xvals, **kwargs)
  
          loglike = copy.deepcopy(lnlp1['loglike'])
          dloglike = copy.deepcopy(lnlp1['dloglike'])
          dloglike0 = dloglike[-1]
          xup = xvals[-1]
  
          # Search for value at which likelihood changes by a given amount
          for i in range(20):
  
              lims1 = utils.get_parameter_limits(xvals, dloglike,
                                                 cl_limit=cl_limit)
              delta = np.abs(dloglike0 - lims1['lnlmax'])
  
              if np.abs(delta - dloglike_thresh) < 0.1:
                  break
  
              if not np.isfinite(lims1['ul']) or delta < 1.0:
                  xup = 2.0 * xvals[-1]
              else:
                  xup = lims1['ul']
  
              lnlp = self.profile(name, parName, logemin=logemin,
                                  logemax=logemax,
                                  reoptimize=True, xvals=[xup], **kwargs)
              dloglike0 = lnlp['dloglike']
              loglike0 = lnlp['loglike']
  
              loglike = np.concatenate((loglike, loglike0))
              dloglike = np.concatenate((dloglike, dloglike0))
              xvals = np.concatenate((xvals, [xup]))
              isort = np.argsort(xvals)
              dloglike = dloglike[isort]
              loglike = loglike[isort]
              xvals = xvals[isort]
  
          if np.isfinite(lims1['ll']):
              xlo = np.concatenate(
                  ([0.0], np.linspace(lims1['ll'], xval0, (npts + 1) // 2 - 1)))
          elif np.abs(dloglike[0] - lims1['lnlmax']) > 0.1:
              xlo = np.linspace(0.0, xval0, (npts + 1) // 2)
          else:
              xlo = np.array([0.0, xval0])
  
          if np.isfinite(lims1['ul']):
              xhi = np.linspace(xval0, lims1['ul'], npts + 1 - len(xlo))[1:]
          else:
              xhi = np.linspace(xval0, lims0['ul'], npts + 1 - len(xlo))[1:]
  
          xvals = np.sort(np.concatenate((xlo, xhi)))
          return xvals
  
      def profile(self, name, parName, logemin=None, logemax=None,
                  reoptimize=False,
                  xvals=None, npts=None, savestate=True, **kwargs):
          """Profile the likelihood for the given source and parameter.
  
          Parameters
          ----------
  
          name : str
             Source name.
  
          parName : str
             Parameter name.
  
          reoptimize : bool
             Re-fit nuisance parameters at each step in the scan.  Note
             that enabling this option will only re-fit parameters that
             were free when the method was executed.
  
          Returns
          -------
  
          lnlprofile : dict
             Dictionary containing results of likelihood scan.
  
          """
          # Find the source
          name = self.roi.get_source_by_name(name).name
  
          par = self.like.normPar(name)
          parName = self.like.normPar(name).getName()
          idx = self.like.par_index(name, parName)
          bounds = self.like.model[idx].getBounds()
          value = self.like.model[idx].getValue()
          loge_bounds = self.loge_bounds
  
          optimizer = kwargs.get('optimizer', self.config['optimizer'])
  
          if savestate:
              saved_state = self._latch_state()
  
          # If parameter is fixed temporarily free it
          par.setFree(True)
          if optimizer['optimizer'] == 'NEWTON':
              self._create_fitcache()
  
          if logemin is not None or logemax is not None:
              loge_bounds = self.set_energy_range(logemin, logemax)
          else:
              loge_bounds = self.loge_bounds
  
          loglike0 = -self.like()
  
          if xvals is None:
  
              err = par.error()
              val = par.getValue()
              if err <= 0 or val <= 3 * err:
                  xvals = 10 ** np.linspace(-2.0, 2.0, 51)
                  if val < xvals[0]:
                      xvals = np.insert(xvals, val, 0)
              else:
                  xvals = np.linspace(0, 1, 25)
                  xvals = np.concatenate((-1.0 * xvals[1:][::-1], xvals))
                  xvals = val * 10 ** xvals
  
          if np.isnan(xvals).any():
              raise RuntimeError(
                  "Parameter scan points for %s::%s include infinite value." % (name, parName))
  
          # Update parameter bounds to encompass scan range
          try:
              self.like[idx].setBounds(min(min(xvals), value, bounds[0]),
                                       max(max(xvals), value, bounds[1]))
          except RuntimeError:
              self.logger.warning(
                  "Caught failure on setBounds for %s::%s." % (name, parName))
  
          o = {'xvals': xvals,
               'npred': np.zeros(len(xvals)),
               'npred_wt': np.zeros(len(xvals)),
               'dnde': np.zeros(len(xvals)),
               'flux': np.zeros(len(xvals)),
               'eflux': np.zeros(len(xvals)),
               'dloglike': np.zeros(len(xvals)),
               'loglike': np.zeros(len(xvals))
               }
  
          if reoptimize and hasattr(self.like.components[0].logLike,
                                    'setUpdateFixedWeights'):
  
              for c in self.components:
                  c.like.logLike.setUpdateFixedWeights(False)
  
          for i, x in enumerate(xvals):
  
              try:
                  self.like[idx] = x
              except RuntimeError:
                  self.logger.warning(
                      "Caught failure on set for %s::%s: %.2f" % (name, parName, x))
  
              if self.like.nFreeParams() > 1 and reoptimize:
                  # Only reoptimize if not all frozen
                  self.like.freeze(idx)
                  fit_output = self._fit(errors=False, **optimizer)
                  loglike1 = fit_output['loglike']
                  self.like.thaw(idx)
              else:
                  loglike1 = -self.like()
  
              flux = self.like[name].flux(10 ** loge_bounds[0],
                                          10 ** loge_bounds[1])
              eflux = self.like[name].energyFlux(10 ** loge_bounds[0],
                                                 10 ** loge_bounds[1])
              prefactor = self.like[idx]
  
              o['dloglike'][i] = loglike1 - loglike0
              o['loglike'][i] = loglike1
              o['dnde'][i] = prefactor.getTrueValue()
              o['flux'][i] = flux
              o['eflux'][i] = eflux
  
              cs = self.model_counts_spectrum(name,
                                              loge_bounds[0],
                                              loge_bounds[1], summed=True)
              o['npred'][i] += np.sum(cs)
              cs_wt = self.model_counts_spectrum(name,
                                                 loge_bounds[0],
                                                 loge_bounds[1], summed=True, weighted=True)
              o['npred_wt'][i] += np.sum(cs_wt)
  
          self.like[idx] = value
  
          if reoptimize and hasattr(self.like.components[0].logLike,
                                    'setUpdateFixedWeights'):
  
              for c in self.components:
                  c.like.logLike.setUpdateFixedWeights(True)
  
          # Restore model parameters to original values
          if savestate:
              saved_state.restore()
  
          self.like[idx].setBounds(*bounds)
          if logemin is not None or logemax is not None:
              self.set_energy_range(*loge_bounds)
  
          return o
  
      def constrain_norms(self, srcNames, cov_scale=1.0):
          """Constrain the normalizations of one or more sources by
          adding gaussian priors with sigma equal to the parameter
          error times a scaling factor."""
  
          # Get the covariance matrix
  
          for name in srcNames:
              par = self.like.normPar(name)
  
              err = par.error()
              val = par.getValue()
  
              if par.error() == 0.0 or not par.isFree():
                  continue
  
              self.add_gauss_prior(name, par.getName(),
                                   val, err * cov_scale)
  
      def add_gauss_prior(self, name, parName, mean, sigma):
  
          par = self.like[name].funcs["Spectrum"].params[parName]
          par.addGaussianPrior(mean, sigma)
  
      def remove_prior(self, name, parName):
  
          par = self.like[name].funcs["Spectrum"].params[parName]
          par.removePrior()
  
      def remove_priors(self):
          """Clear all priors."""
  
          for src in self.roi.sources:
  
              for par in self.like[src.name].funcs["Spectrum"].params.values():
                  par.removePrior()
  
      def _create_optObject(self, **kwargs):
          """ Make MINUIT or NewMinuit type optimizer object """
  
          optimizer = kwargs.get('optimizer',
                                 self.config['optimizer']['optimizer'])
  
          if optimizer.upper() == 'MINUIT':
              optObject = pyLike.Minuit(self.like.logLike)
          elif optimizer.upper() == 'NEWMINUIT':
              optObject = pyLike.NewMinuit(self.like.logLike)
          else:
              optFactory = pyLike.OptimizerFactory_instance()
              optObject = optFactory.create(str(optimizer), self.like.logLike)
          return optObject
  
      def _fit_optimizer_iter(self, **kwargs):
  
          min_fit_quality = kwargs.get('min_fit_quality', 3)
          retries = kwargs.get('retries', 3)
          covar = kwargs.get('covar', True)
  
          num_free = self.like.nFreeParams()
          o = {'values': np.ones(num_free) * np.nan,
               'errors': np.ones(num_free) * np.nan,
               'covariance': np.ones((num_free, num_free)) * np.nan,
               'correlation': np.ones((num_free, num_free)) * np.nan,
               'indices': np.zeros(num_free, dtype=int),
               'is_norm': np.empty(num_free, dtype=bool),
               'src_names': num_free * [None],
               'par_names': num_free * [None],
               'fit_quality': 3,
               'fit_status': 0,
               'fit_success': True,
               'edm': 0,
               'niter': 0,
               'loglike': np.nan}
  
          if num_free == 0:
              o['loglike'] = -self.like()
              return o
  
          quality = 0
          niter = 0
          while niter < retries:
              niter += 1
              quality, status, edm, loglike = self._fit_optimizer(**kwargs)
  
              if quality >= min_fit_quality and status == 0:
                  break
  
              self.logger.debug("Retry fit iter: %i quality: %i edm: %8.4f loglike: %12.3f",
                                niter, quality, edm, loglike)
  
          o['fit_quality'] = quality
          o['fit_status'] = status
          o['edm'] = edm
          o['niter'] = niter
          o['loglike'] = loglike
  
          if quality < min_fit_quality or o['fit_status']:
              o['fit_success'] = False
  
          if (covar and self.like.covariance is not None and
                  len(self.like.covariance) == num_free):
              o['covariance'] = np.array(self.like.covariance)
              o['errors'] = np.diag(o['covariance'])**0.5
              errinv = 1. / o['errors']
              o['correlation'] = \
                  o['covariance'] * errinv[:, np.newaxis] * errinv[np.newaxis, :]
  
          free_params = self.get_params(True)
          for i, p in enumerate(free_params):
              o['values'][i] = p['value']
              o['errors'][i] = p['error']
              o['indices'][i] = p['idx']
              o['src_names'][i] = p['src_name']
              o['par_names'][i] = p['par_name']
              o['is_norm'][i] = p['is_norm']
  
          return o
  
      def _fit_optimizer(self, **kwargs):
  
          errors = kwargs.get('errors', True)
          optObject = self._create_optObject(optimizer=kwargs.get('optimizer',
                                                                  'MINUIT'))
  
          kw = {}
  
          quality = 3
          status = 0
          edm = 0
          loglike = 0
  
          try:
              if errors:
                  kw['verbosity'] = kwargs.get('verbosity', 0)
                  kw['tol'] = kwargs.get('tol', None)
                  kw['covar'] = kwargs.get('covar', True)
                  kw['optObject'] = optObject
                  self.like.fit(**kw)
              else:
                  kw['verbosity'] = kwargs.get('verbosity', 0)
                  kw['tol'] = kwargs.get('tol', None)
                  kw['optObject'] = optObject
                  self.like.optimize(**kw)
  
              status = optObject.getRetCode()
  
              if isinstance(optObject, pyLike.Minuit):
                  quality = optObject.getQuality()
  
              if isinstance(optObject, pyLike.Minuit) or \
                      isinstance(optObject, pyLike.NewMinuit):
                  edm = optObject.getDistance()
  
              loglike = optObject.stat().value()
  
          except Exception:
              self.logger.error('Likelihood optimization failed.', exc_info=True)
              quality = 0
              status = 1
  
          return quality, status, edm, loglike
  
      def _fit(self, **kwargs):
  
          optimizer = kwargs.get('optimizer',
                                 self.config['optimizer']['optimizer']).upper()
  
          if optimizer == 'NEWTON':
  
              params = self.get_params(True)
              for p in params:
  
                  if p['free'] and not p['is_norm']:
                      optimizer = 'MINUIT'
                      kwargs['optimizer'] = 'MINUIT'
                      self.logger.debug(
                          'Found non-norm free parameter.  Reverting to MINUIT.')
                      break
  
          if optimizer == 'NEWTON':
              return self._fit_newton(**kwargs)
          else:
              return self._fit_optimizer_iter(**kwargs)
  
      def fit(self, update=True, **kwargs):
          """Run the likelihood optimization.  This will execute a fit of all
          parameters that are currently free in the model and update the
          charateristics of the corresponding model components (TS,
          npred, etc.).  The fit will be repeated N times (set with the
          `retries` parameter) until a fit quality greater than or equal
          to `min_fit_quality` and a fit status code of 0 is obtained.
          If the fit does not succeed after N retries then all parameter
          values will be reverted to their state prior to the execution
          of the fit.
  
          Parameters
          ----------
  
          update : bool
             Update the model dictionary for all sources with free
             parameters.
  
          tol : float
             Set the optimizer tolerance.
  
          verbosity : int
             Set the optimizer output level.
  
          optimizer : str
             Set the likelihood optimizer (e.g. MINUIT or NEWMINUIT).
  
          retries : int
             Set the number of times to rerun the fit when the fit quality
             is < 3.
  
          min_fit_quality : int
             Set the minimum fit quality.  If the fit quality is smaller
             than this value then all model parameters will be
             restored to their values prior to the fit.
  
          reoptimize : bool
             Refit background sources when updating source properties
             (TS and likelihood profiles).
  
          Returns
          -------
  
          fit : dict
             Dictionary containing diagnostic information from the fit
             (fit quality, parameter covariances, etc.).
  
          """
  
          loglevel = kwargs.pop('loglevel', self.loglevel)
          self.logger.log(loglevel, "Starting fit.")
  
          # Extract options from kwargs
          config = copy.deepcopy(self.config['optimizer'])
          config.setdefault('covar', True)
          config.setdefault('reoptimize', False)
          config = utils.merge_dict(config, kwargs)
  
          num_free = self.like.nFreeParams()
  
          loglike0 = -self.like()
  
          # Initialize output dict
          o = {'fit_quality': 3,
               'fit_status': 0,
               'fit_success': True,
               'dloglike': 0.0,
               'edm': 0.0,
               'loglike': loglike0,
               'covariance': None,
               'correlation': None,
               'values': np.ones(num_free) * np.nan,
               'errors': np.ones(num_free) * np.nan,
               'indices': np.zeros(num_free, dtype=int),
               'is_norm': np.empty(num_free, dtype=bool),
               'src_names': num_free * [None],
               'par_names': num_free * [None],
               'config': config
               }
  
          if not num_free:
              self.logger.log(loglevel, "Skipping fit.  No free parameters.")
              return o
  
          saved_state = LikelihoodState(self.like)
  
          fit_output = self._fit(**config)
          o.update(fit_output)
  
          self.logger.debug("Fit complete.")
          o['dloglike'] = o['loglike'] - loglike0
  
          if not o['fit_success']:
              self.logger.error('%s failed with status code %i fit quality %i',
                                config['optimizer'], o['fit_status'],
                                o['fit_quality'])
              saved_state.restore()
              return o
  
          if update:
  
              free_params = self.get_params(True)
              self._extract_correlation(o, free_params)
              for name in self.like.sourceNames():
                  freePars = self.get_free_source_params(name)
                  if len(freePars) == 0:
                      continue
                  self.update_source(name, reoptimize=config['reoptimize'])
  
              # Update roi model counts
              self._update_roi()
  
          self.logger.log(loglevel, "Fit returned successfully. " +
                          "Quality: %3i Status: %3i",
                          o['fit_quality'], o['fit_status'])
          self.logger.log(loglevel, "LogLike: %12.3f DeltaLogLike: %12.3f ",
                          o['loglike'], o['dloglike'])
          return o
  
      def _create_fitcache(self, **kwargs):
  
          create_fitcache = kwargs.get('create_fitcache', False)
          if self._fitcache is not None and not create_fitcache:
              return self._fitcache
  
          tol = kwargs.get('tol', self.config['optimizer']['tol'])
          max_iter = kwargs.get('max_iter',
                                self.config['optimizer']['max_iter'])
          init_lambda = kwargs.get('init_lambda',
                                   self.config['optimizer']['init_lambda'])
          use_reduced = kwargs.get('use_reduced', True)
  
          self.logger.debug('Creating FitCache')
          self.logger.debug('\ntol: %.5g\nmax_iter: %i\ninit_lambda: %.5g',
                            tol, max_iter, init_lambda)
  
          params = self.get_params()
          self._fitcache = FitCache(self.like, params,
                                    tol, max_iter, init_lambda, use_reduced)
  
          return self._fitcache
  
      def _fit_newton(self, fitcache=None, ebin=None, **kwargs):
          """Fast fitting method using newton fitter."""
  
          tol = kwargs.get('tol', self.config['optimizer']['tol'])
          max_iter = kwargs.get('max_iter',
                                self.config['optimizer']['max_iter'])
          init_lambda = kwargs.get('init_lambda',
                                   self.config['optimizer']['init_lambda'])
          use_reduced = kwargs.get('use_reduced', True)
  
          free_params = self.get_params(True)
          free_norm_params = [p for p in free_params if p['is_norm'] is True]
  
          if len(free_params) != len(free_norm_params):
              msg = 'Executing Newton fitter with one ' + \
                  'or more free shape parameters.'
              self.logger.error(msg)
              raise Exception(msg)
  
          verbosity = kwargs.get('verbosity', 0)
          if fitcache is None:
              fitcache = self._create_fitcache(**kwargs)
  
          fitcache.update(self.get_params(),
                          tol, max_iter, init_lambda, use_reduced)
  
          logemin = self.loge_bounds[0]
          logemax = self.loge_bounds[1]
          imin = int(utils.val_to_edge(self.log_energies, logemin)[0])
          imax = int(utils.val_to_edge(self.log_energies, logemax)[0])
  
          if ebin is not None:
              fitcache.fitcache.setEnergyBin(ebin)
          elif imin == 0 and imax == self.enumbins:
              fitcache.fitcache.setEnergyBin(-1)
          else:
              fitcache.fitcache.setEnergyBins(imin, imax)
  
          num_free = len(free_norm_params)
  
          o = {'fit_status': 0,
               'fit_quality': 3,
               'fit_success': True,
               'edm': 0,
               'loglike': None,
               'values': np.ones(num_free) * np.nan,
               'errors': np.ones(num_free) * np.nan,
               'indices': np.zeros(num_free, dtype=int),
               'is_norm': np.empty(num_free, dtype=bool),
               'src_names': num_free * [None],
               'par_names': num_free * [None],
               }
  
          if num_free == 0:
              return o
  
          ref_vals = np.array(fitcache.fitcache.refValues())
          free = np.array(fitcache.fitcache.currentFree())
          norm_vals = ref_vals[free]
  
          norm_idxs = []
          for i, p in enumerate(free_norm_params):
  
              norm_idxs += [p['idx']]
              o['indices'][i] = p['idx']
              o['src_names'][i] = p['src_name']
              o['par_names'][i] = p['par_name']
              o['is_norm'][i] = p['is_norm']
  
          o['fit_status'] = fitcache.fit(verbose=verbosity)
          o['edm'] = fitcache.fitcache.currentEDM()
  
          pars, errs, cov = fitcache.get_pars()
          pars *= norm_vals
          errs *= norm_vals
          cov = cov * np.outer(norm_vals, norm_vals)
  
          o['values'] = pars
          o['errors'] = errs
          o['covariance'] = cov
          errinv = np.zeros_like(o['errors'])
          m = o['errors'] > 0
          errinv[m] = 1. / o['errors'][m]
          o['correlation'] = o['covariance'] * np.outer(errinv, errinv)
  
          if o['fit_status'] in [-2, 0]:
              for idx, val, err in zip(norm_idxs, pars, errs):
                  self._set_value_bounded(idx, val)
                  self.like[idx].setError(err)
              self.like.syncSrcParams()
              o['fit_success'] = True
          else:
              o['fit_success'] = False
  
          if o['fit_status']:
              self.logger.error('Error in NEWTON fit. Fit Status: %i',
                                o['fit_status'])
  
          # FIXME: Figure out why currentLogLike gets out of sync
          #loglike = fitcache.fitcache.currentLogLike()
          #prior_vals, prior_errs, has_prior = gtutils.get_priors(self.like)
          #loglike -= np.sum(has_prior) * np.log(np.sqrt(2 * np.pi))
          loglike = -self.like()
          o['loglike'] = loglike
  
          return o
  
      def fit_correlation(self):
  
          saved_state = LikelihoodState(self.like)
          self.free_sources(False, loglevel=logging.DEBUG)
          self.free_sources(pars='norm', loglevel=logging.DEBUG)
          fit_results = self.fit(loglevel=logging.DEBUG, min_fit_quality=2)
          free_params = self.get_params(True)
          self._extract_correlation(fit_results, free_params)
          saved_state.restore()
  
      def _extract_correlation(self, fit_results, free_params):
  
          for i, p0 in enumerate(free_params):
              if not p0['is_norm']:
                  continue
  
              src = self.roi[p0['src_name']]
              for j, p1 in enumerate(free_params):
  
                  if not p1['is_norm']:
                      continue
  
                  src['correlation'][p1['src_name']] = fit_results[
                      'correlation'][i, j]
  
      def load_xml(self, xmlfile):
          """Load model definition from XML.
  
          Parameters
          ----------
          xmlfile : str
              Name of the input XML file.
  
          """
  
          self.logger.info('Loading XML')
  
          for c in self.components:
              c.load_xml(xmlfile)
  
          for name in self.like.sourceNames():
              self.update_source(name)
  
          self._fitcache = None
  
          self.logger.info('Finished Loading XML')
  
      def write_xml(self, xmlfile):
          """Save current model definition as XML file.
  
          Parameters
          ----------
          xmlfile : str
              Name of the output XML file.
  
          """
          for c in self._components:
              c.write_xml(xmlfile)
  
      def load_parameters_from_yaml(self, yamlfile, update_sources=False):
          """Load model parameters from yaml
  
          Parameters
          ----------
          yamlfile : str
              Name of the input yaml file.
          """
          d = utils.load_yaml(yamlfile)
          for src, src_pars in d.items():
              for par_name, par_dict in src_pars.items():
                  if par_name in ['SpectrumType']:
                      continue
                  par_value = par_dict.get('value', None)
                  par_error = par_dict.get('error', None)
                  par_scale = par_dict.get('scale', None)
                  par_min = par_dict.get('min', None)
                  par_max = par_dict.get('max', None)
                  par_free = par_dict.get('free', None)
                  if par_min is not None and par_max is not None:
                      par_bounds = [par_min, par_max]
                  else:
                      par_bounds = None
                  try:
                      self.set_parameter(src, par_name, par_value, true_value=False,
                                         scale=par_scale, bounds=par_bounds, error=par_error,
                                         update_source=update_sources)
                  except RuntimeError as msg:
                      self.logger.warn(msg)
                      self.logger.warn("Did not set parameter %s:%s"%(src,par_name))                    
                      continue
                  except Exception as msg:
                      self.logger.warn(msg)
                      continue
                  if par_free is not None:
                      self.free_parameter(src, par_name, par_free)
          self._sync_params_state()
  
      def _restore_counts_maps(self):
          """
          Revert counts maps to their state prior to injecting any simulated
          components.
          """
  
          for c in self.components:
              c.restore_counts_maps()
  
          if hasattr(self.like.components[0].logLike, 'setCountsMap'):
              self._init_roi_model()
          else:
              self.write_xml('tmp')
              self._like = SummedLikelihood()
              for i, c in enumerate(self._components):
                  c._create_binned_analysis()
                  self._like.addComponent(c.like)
              self._init_roi_model()
              self.load_xml('tmp')
  
      def simulate_source(self, src_dict=None):
          """
          Inject simulated source counts into the data.
  
          Parameters
          ----------
          src_dict : dict
             Dictionary defining the spatial and spectral properties of
             the source that will be injected.
  
          """
  
          self._fitcache = None
  
          if src_dict is None:
              src_dict = {}
          else:
              src_dict = copy.deepcopy(src_dict)
  
          skydir = wcs_utils.get_target_skydir(src_dict, self.roi.skydir)
  
          src_dict.setdefault('ra', skydir.ra.deg)
          src_dict.setdefault('dec', skydir.dec.deg)
          src_dict.setdefault('SpatialModel', 'PointSource')
          src_dict.setdefault('SpatialWidth', 0.3)
          src_dict.setdefault('Index', 2.0)
          src_dict.setdefault('Prefactor', 1E-13)
  
          self.add_source('mcsource', src_dict, free=True,
                          init_source=False)
          for c in self.components:
              c.simulate_roi('mcsource', clear=False)
  
          self.delete_source('mcsource')
  
          if hasattr(self.like.components[0].logLike, 'setCountsMap'):
              self._init_roi_model()
          else:
              self.write_xml('tmp')
              self._like = SummedLikelihood()
              for i, c in enumerate(self._components):
                  c._create_binned_analysis('tmp.xml')
                  self._like.addComponent(c.like)
              self._init_roi_model()
              self.load_xml('tmp')
  
      def simulate_roi(self, name=None, randomize=True, restore=False):
          """Generate a simulation of the ROI using the current best-fit model
          and replace the data counts cube with this simulation.  The
          simulation is created by generating an array of Poisson random
          numbers with expectation values drawn from the model cube of
          the binned analysis instance.  This function will update the
          counts cube both in memory and in the source map file.  The
          counts cube can be restored to its original state by calling
          this method with ``restore`` = True.
  
          Parameters
          ----------
          name : str
             Name of the model component to be simulated.  If None then
             the whole ROI will be simulated.
  
          restore : bool
             Restore the data counts cube to its original state.
          """
  
          self.logger.info('Simulating ROI')
  
          self._fitcache = None
  
          if restore:
              self.logger.info('Restoring')
              self._restore_counts_maps()
              self.logger.info('Finished')
              return
  
          for c in self.components:
              c.simulate_roi(name=name, clear=True, randomize=randomize)
  
          if hasattr(self.like.components[0].logLike, 'setCountsMap'):
              self._init_roi_model()
          else:
              self.write_xml('tmp')
              self._like = SummedLikelihood()
              for i, c in enumerate(self._components):
                  c._create_binned_analysis('tmp.xml')
                  self._like.addComponent(c.like)
              self._init_roi_model()
              self.load_xml('tmp')
  
          self.logger.info('Finished')
  
      def write_model_map(self, model_name, name=None):
          """Save the counts model map to a FITS file.
  
          Parameters
          ----------
          model_name : str
              String that will be append to the name of the output file.
          name : str
              Name of the component.
  
          Returns
          -------
  
          """
          maps = [c.write_model_map(model_name, name) for c in self.components]
  
          outfile = os.path.join(self.workdir,
                                 'mcube_%s.fits' % (model_name))
  
          mmap = Map.from_geom(self.geom)
          for m in maps:
              mmap.coadd(m)
          mmap.write(outfile, overwrite=True, conv='fgst-ccube')
          return [mmap] + maps
  
      def write_weight_map(self, model_name):
          """Save the counts model map to a FITS file.
  
          Parameters
          ----------
          model_name : str
              String that will be append to the name of the output file.
  
          Returns
          -------
  
          """
          maps = [c.write_weight_map(model_name) for c in self.components]
  
          outfile = os.path.join(self.workdir,
                                 'wcube_%s.fits' % (model_name))
  
          wmap = Map.from_geom(self.geom)
          # FIXME: Should we average weights maps rather than coadding?
          for m in maps:
              wmap.coadd(m)
          wmap.write(outfile, overwrite=True, conv='fgst-ccube')
          return [wmap] + maps
  
      def print_roi(self, loglevel=logging.INFO):
          """Print information about the spectral and spatial properties
          of the ROI (sources, diffuse components)."""
          self.logger.log(loglevel, '\n' + str(self.roi))
  
      def print_params(self, allpars=False, loglevel=logging.INFO):
          """Print information about the model parameters (values,
          errors, bounds, scale)."""
  
          pars = self.get_params()
  
          o = '\n'
          o += '%4s %-20s%10s%10s%10s%10s%10s%5s\n' % (
              'idx', 'parname', 'value', 'error',
              'min', 'max', 'scale', 'free')
  
          o += '-' * 80 + '\n'
  
          src_pars = collections.OrderedDict()
          for p in pars:
  
              src_pars.setdefault(p['src_name'], [])
              src_pars[p['src_name']] += [p]
  
          free_sources = []
          for k, v in src_pars.items():
  
              for p in v:
                  if not p['free']:
                      continue
  
                  free_sources += [k]
  
          for k, v in src_pars.items():
  
              if not allpars and k not in free_sources:
                  continue
  
              o += '%s\n' % k
              for p in v:
  
                  o += '%4i %-20.19s' % (p['idx'], p['par_name'])
                  o += '%10.3g%10.3g' % (p['value'], p['error'])
                  o += '%10.3g%10.3g%10.3g' % (p['min'], p['max'],
                                               p['scale'])
  
                  if p['free']:
                      o += '    *'
                  else:
                      o += '     '
  
                  o += '\n'
  
          self.logger.log(loglevel, o)
  
      def print_model(self, loglevel=logging.INFO):
  
          o = '\n'
          o += '%-20s%8s%8s%7s%10s%10s%12s%5s\n' % (
              'sourcename', 'offset', 'norm', 'eflux', 'index',
              'ts', 'npred', 'free')
          o += '-' * 80 + '\n'
  
          for s in sorted(self.roi.sources, key=lambda t: t['offset']):
              if s.diffuse:
                  continue
  
              normVal = self.like.normPar(s.name).getValue()
              fixed = self.like[s.name].fixedSpectrum()
  
              if fixed:
                  free_str = ' '
              else:
                  free_str = '*'
  
              if s['SpectrumType'] == 'PowerLaw':
                  index = s['dnde1000_index']
              else:
                  index = 0.5 * (s['dnde1000_index'] +
                                 s['dnde10000_index'])
  
              o += '%-20.19s%8.3f%8.3f%10.3g%7.2f%10.2f%12.1f%5s\n' % (
                  s['name'], s['offset'], normVal, s['eflux'], index,
                  s['ts'], s['npred'], free_str)
  
          for s in sorted(self.roi.sources, key=lambda t: t['offset']):
              if not s.diffuse:
                  continue
  
              normVal = self.like.normPar(s.name).getValue()
              fixed = self.like[s.name].fixedSpectrum()
  
              if fixed:
                  free_str = ' '
              else:
                  free_str = '*'
  
              if s['SpectrumType'] == 'PowerLaw':
                  index = s['dnde1000_index']
              else:
                  index = 0.5 * (s['dnde1000_index'] +
                                 s['dnde10000_index'])
  
              o += '%-20.19s%8s%8.3f%10.3g%7.2f%10.2f%12.1f%5s\n' % (
                  s['name'],
                  '---', normVal, s['eflux'], index, s['ts'], s['npred'], free_str)
  
          self.logger.log(loglevel, o)
  
      def load_roi(self, infile, reload_sources=False, params=None, mask=None):
          """This function reloads the analysis state from a previously
          saved instance generated with
          `~fermipy.gtanalysis.GTAnalysis.write_roi`.
  
          Parameters
          ----------
  
          infile : str
  
          reload_sources : bool
             Regenerate source maps for non-diffuse sources.
  
          params : str
              Path to a yaml file with updated parameter values
          
          mask : str
              Path to a fits file with an updated mask
  
          """
  
          infile = utils.resolve_path(infile, workdir=self.workdir)
          roi_file, roi_data = utils.load_data(infile, workdir=self.workdir)
  
          self.logger.info('Loading ROI file: %s', roi_file)
  
          key_map = {'dfde': 'dnde',
                     'dfde100': 'dnde100',
                     'dfde1000': 'dnde1000',
                     'dfde10000': 'dnde10000',
                     'dfde_index': 'dnde_index',
                     'dfde100_index': 'dnde100_index',
                     'dfde1000_index': 'dnde1000_index',
                     'dfde10000_index': 'dnde10000_index',
                     'e2dfde': 'e2dnde',
                     'e2dfde100': 'e2dnde100',
                     'e2dfde1000': 'e2dnde1000',
                     'e2dfde10000': 'e2dnde10000',
                     'Npred': 'npred',
                     'Npred_wt': 'npred_wt',
                     'logLike': 'loglike',
                     'dlogLike': 'dloglike',
                     'emin': 'e_min',
                     'ectr': 'e_ctr',
                     'emax': 'e_max',
                     'logemin': 'loge_min',
                     'logectr': 'loge_ctr',
                     'logemax': 'loge_max',
                     'ref_dfde': 'ref_dnde',
                     'ref_e2dfde': 'ref_e2dnde',
                     'ref_dfde_emin': 'ref_dnde_e_min',
                     'ref_dfde_emax': 'ref_dnde_e_max',
                     }
  
          self._roi_data = utils.update_keys(roi_data['roi'], key_map)
  
          if 'erange' in self._roi_data:
              self._roi_data['loge_bounds'] = self._roi_data.pop('erange')
  
          self._loge_bounds = self._roi_data.setdefault('loge_bounds',
                                                        self.loge_bounds)
  
          sources = roi_data.pop('sources')
          sources = utils.update_keys(sources, key_map)
          for k0, v0 in sources.items():
              for k, v in defaults.source_flux_output.items():
                  if k not in v0:
                      continue
                  if v[2] == float and isinstance(v0[k], np.ndarray):
                      sources[k0][k], sources[k0][k + '_err'] \
                          = v0[k][0], v0[k][1]
  
          self.roi.load_sources(sources.values())
          for i, c in enumerate(self.components):
              if 'src_expscale' in self._roi_data['components'][i]:
                  c._src_expscale = copy.deepcopy(self._roi_data['components']
                                                  [i]['src_expscale'])
  
          self._create_likelihood(infile)
          self.set_energy_range(self.loge_bounds[0], self.loge_bounds[1])
  
          if params is not None:
              self.load_parameters_from_yaml(params)
  
          if mask is not None:
              self.set_weights_map(mask, update_roi=False)
  
          if reload_sources:
              names = [s.name for s in self.roi.sources if not s.diffuse]
              self.reload_sources(names, False)
  
          self.logger.info('Finished Loading ROI')
  
      def write_roi(self, outfile=None,
                    save_model_map=False, **kwargs):
          """Write current state of the analysis to a file.  This method
          writes an XML model definition, a ROI dictionary, and a FITS
          source catalog file.  A previously saved analysis state can be
          reloaded from the ROI dictionary file with the
          `~fermipy.gtanalysis.GTAnalysis.load_roi` method.
  
          Parameters
          ----------
  
          outfile : str
              String prefix of the output files.  The extension of this
              string will be stripped when generating the XML, YAML and
              npy filenames.
  
          make_plots : bool
              Generate diagnostic plots.
  
          save_model_map : bool
              Save the current counts model to a FITS file.
  
          """
          # extract the results in a convenient format
  
          make_plots = kwargs.get('make_plots', False)
          save_weight_map = kwargs.get('save_weight_map', False)
  
          if outfile is None:
              pathprefix = os.path.join(self.config['fileio']['workdir'],
                                        'results')
          elif not os.path.isabs(outfile):
              pathprefix = os.path.join(self.config['fileio']['workdir'],
                                        outfile)
          else:
              pathprefix = outfile
  
          pathprefix = utils.strip_suffix(pathprefix,
                                          ['fits', 'yaml', 'npy'])
  #        pathprefix, ext = os.path.splitext(pathprefix)
          prefix = os.path.basename(pathprefix)
  
          xmlfile = pathprefix + '.xml'
          fitsfile = pathprefix + '.fits'
          npyfile = pathprefix + '.npy'
  
          self.write_xml(xmlfile)
          self.write_fits(fitsfile)
  
          if not self.config['gtlike']['use_external_srcmap']:
              for c in self.components:
                  c.like.logLike.saveSourceMaps(str(c.files['srcmap']))
  
          if save_model_map:
              self.write_model_map(prefix)
  
          if save_weight_map:
              self.write_weight_map(prefix)
  
          o = {}
          o['roi'] = copy.deepcopy(self._roi_data)
          o['config'] = copy.deepcopy(self.config)
          o['version'] = fermipy.__version__
          o['stversion'] = fermipy.get_st_version()
          o['sources'] = {}
  
          for s in self.roi.sources:
              o['sources'][s.name] = copy.deepcopy(s.data)
  
          for i, c in enumerate(self.components):
              o['roi']['components'][i][
                  'src_expscale'] = copy.deepcopy(c.src_expscale)
  
          self.logger.info('Writing %s...', npyfile)
          np.save(npyfile, o)
  
          if make_plots:
              self.make_plots(prefix, None,
                              **kwargs.get('plotting', {}))
  
      def write_fits(self, fitsfile):
  
          self.logger.info('Writing %s...', fitsfile)
  
          tab = self.roi.create_table()
          hdu_data = fits.table_to_hdu(tab)
          hdu_data.name = 'CATALOG'
  
          tab_srcs = []
          tab_params = []
          for i, c in enumerate(self.components):
              tab = self.roi.create_source_table()
              tab['component'] = i
              tab['expscale'] = np.nan
  
              for k, v in c.src_expscale.items():
                  m = tab['source_name'] == k
                  tab['expscale'][m] = v
  
              tab_srcs += [tab]
  
              tab = self.roi.create_param_table()
              tab['component'] = i
              tab_params += [tab]
  
          from astropy.table import vstack
  
          hdu_srcs = fits.table_to_hdu(vstack(tab_srcs))
          hdu_srcs.name = 'SOURCES'
          hdu_params = fits.table_to_hdu(vstack(tab_params))
          hdu_params.name = 'PARAMS'
          hdu_roi = fits.table_to_hdu(self.create_roi_table())
          hdu_roi.name = 'ROI'
  
          hdus = [fits.PrimaryHDU(), hdu_data, hdu_roi, hdu_srcs, hdu_params]
          hdus[0].header['CONFIG'] = json.dumps(self.config)
          hdus[1].header['CONFIG'] = json.dumps(self.config)
          fits_utils.write_hdus(hdus, fitsfile)
  
      def create_roi_table(self):
  
          rd = copy.deepcopy(self._roi_data)
          tab = fits_utils.dict_to_table(rd)
          tab['component'] = -1
          row_dict = {}
          for i, c in enumerate(rd['components']):
              c['component'] = i
              c['loge_bounds'] = rd['loge_bounds']
  
              row = []
              for k in tab.columns:
  
                  shape = tab.columns[k].shape
                  ndim = tab.columns[k].ndim
                  if ndim == 1:
                      val = c[k]
                  else:
                      val = np.ones(shape[1:]) * np.nan
                      val[:len(c[k])] = c[k]
                  row += [val]
              tab.add_row(row)
          return tab
  
      def make_plots(self, prefix, mcube_map=None, **kwargs):
          """Make diagnostic plots using the current ROI model."""
  
          #mcube_maps = kwargs.pop('mcube_maps', None)
          if mcube_map is None:
              mcube_map = self.model_counts_map()
  
          plotter = plotting.AnalysisPlotter(self.config['plotting'],
                                             fileio=self.config['fileio'],
                                             logging=self.config['logging'])
          plotter.run(self, mcube_map, prefix=prefix, **kwargs)
  
      def curvature(self, name, **kwargs):
          """Test whether a source shows spectral curvature by comparing
          the likelihood ratio of PowerLaw and LogParabola spectral
          models.
  
          Parameters
          ----------
          name : str
              Source name.
  
          """
  
          name = self.roi.get_source_by_name(name).name
  
          saved_state = LikelihoodState(self.like)
          source = self.components[0].like.logLike.getSource(str(name))
  
          old_spectrum = source.spectrum()
          old_pars = copy.deepcopy(self.roi[name].spectral_pars)
          old_type = self.roi[name]['SpectrumType']
  
          if old_type != 'PowerLaw':
  
              dnde = self.like[name].spectrum()(pyLike.dArg(1000.))
              value, scale = utils.scale_parameter(dnde)
              pars0 = {
                  'Prefactor':
                      {'value': value, 'scale': scale,
                          'min': 1E-5, 'max': 1000., 'free': True},
                  'Index':
                      {'value': 2.0, 'scale': -1.0, 'min': 0.0,
                          'max': 5.0, 'free': False},
                  'Scale':
                      {'value': 1E3, 'scale': 1.0, 'min': 1.,
                          'max': 1E6, 'free': False},
              }
              self.set_source_spectrum(str(name), 'PowerLaw',
                                       spectrum_pars=pars0,
                                       update_source=False)
  
          self.free_source(name, loglevel=logging.DEBUG)
          fit_pl = self._fit(loglevel=logging.DEBUG)
  
          prefactor = self._get_param(name, 'Prefactor')
          index = self._get_param(name, 'Index')
          scale = self._get_param(name, 'Scale')
  
          pars1 = {
              'norm': copy.deepcopy(prefactor),
              'alpha': copy.deepcopy(index),
              'Eb': copy.deepcopy(scale),
          }
          pars1['alpha']['scale'] *= -1
          pars1['alpha']['min'] = -5.0
          pars1['alpha']['max'] = 5.0
  
          self.set_source_spectrum(str(name), 'LogParabola',
                                   spectrum_pars=pars1,
                                   update_source=False)
  
          self.free_source(name, loglevel=logging.DEBUG)
          fit_lp = self._fit(loglevel=logging.DEBUG)
  
          self.set_source_spectrum(str(name), old_type,
                                   spectrum_pars=old_pars,
                                   update_source=False)
  
          pars2 = {
              'Prefactor': copy.deepcopy(prefactor),
              'Index1': copy.deepcopy(index),
              'Cutoff': {'value': 1000.0, 'scale': 1E3,
                         'min': 10.0, 'max': 1E4, 'free': True},
              'Index2': {'value': 1.0, 'scale': 1.0,
                         'min': 1.0, 'max': 1.0, 'free': False},
              'Scale': copy.deepcopy(scale)
          }
  
          self.set_source_spectrum(str(name), 'PLSuperExpCutoff',
                                   spectrum_pars=pars2,
                                   update_source=False)
  
          self.free_source(name, loglevel=logging.DEBUG)
          fit_ple = self._fit(loglevel=logging.DEBUG)
  
          # Revert to initial spectral model
          self.set_source_spectrum(str(name), old_type,
                                   spectrum_pars=old_pars,
                                   update_source=False)
          saved_state.restore()
  
          lp_ts_curv = 2.0 * (fit_lp['loglike'] - fit_pl['loglike'])
          ple_ts_curv = 2.0 * (fit_ple['loglike'] - fit_pl['loglike'])
          o = MutableNamedTuple(ts_curv=lp_ts_curv,
                                lp_ts_curv=lp_ts_curv,
                                ple_ts_curv=ple_ts_curv,
                                loglike_pl=fit_pl['loglike'],
                                loglike_lp=fit_lp['loglike'],
                                loglike_ple=fit_ple['loglike'])
  
          self.logger.info('LogLike_PL: %12.3f LogLike_LP: %12.3f LogLike_PLE: %12.3f',
                           o.loglike_pl, o.loglike_lp, o.loglike_ple)
          self.logger.info('TS_curv:        %.3f (LP)', o.lp_ts_curv)
          self.logger.info('TS_curv:        %.3f (PLE)', o.ple_ts_curv)
          return o
  
      def bowtie(self, name, fd=None, loge=None):
          """Generate a spectral uncertainty band (bowtie) for the given
          source.  This will create an uncertainty band on the
          differential flux as a function of energy by propagating the
          errors on the global fit parameters.  Note that this band only
          reflects the uncertainty for parameters that are currently
          free in the model.
  
          Parameters
          ----------
  
          name : str
             Source name.
  
          fd : FluxDensity
             Flux density object.  If this parameter is None then one
             will be created.
  
          loge : array-like
             Sequence of energies in log10(E/MeV) at which the flux band
             will be evaluated.
  
          """
  
          if loge is None:
              logemin = self.log_energies[0]
              logemax = self.log_energies[-1]
              loge = np.linspace(logemin, logemax, 50)
  
          o = {'energies': 10**loge,
               'log_energies': loge,
               'dnde': np.zeros(len(loge)) * np.nan,
               'dnde_lo': np.zeros(len(loge)) * np.nan,
               'dnde_hi': np.zeros(len(loge)) * np.nan,
               'dnde_err': np.zeros(len(loge)) * np.nan,
               'dnde_ferr': np.zeros(len(loge)) * np.nan,
               'pivot_energy': np.nan}
  
          try:
              if fd is None:
                  fd = FluxDensity.FluxDensity(self.like, name)
          except RuntimeError:
              self.logger.error('Failed to create FluxDensity',
                                exc_info=True)
              return o
  
          dnde = [fd.value(10 ** x) for x in loge]
          dnde_err = [fd.error(10 ** x) for x in loge]
  
          dnde = np.array(dnde)
          dnde_err = np.array(dnde_err)
          m = dnde > 0
  
          fhi = np.zeros_like(dnde)
          flo = np.zeros_like(dnde)
          ferr = np.zeros_like(dnde)
  
          fhi[m] = dnde[m] * (1.0 + dnde_err[m] / dnde[m])
          flo[m] = dnde[m] / (1.0 + dnde_err[m] / dnde[m])
          ferr[m] = 0.5 * (fhi[m] - flo[m]) / dnde[m]
          fhi[~m] = dnde_err[~m]
  
          o['dnde'] = dnde
          o['dnde_lo'] = flo
          o['dnde_hi'] = fhi
          o['dnde_err'] = dnde_err
          o['dnde_ferr'] = ferr
  
          try:
              o['pivot_energy'] = 10 ** utils.interpolate_function_min(loge, o[
                                                                       'dnde_ferr'])
          except Exception:
              self.logger.error('Failed to compute pivot energy',
                                exc_info=True)
  
          return o
  
      def _coadd_maps(self, cmaps, shape, rm, wmaps):
          """
          """
          self._ccube = skymap.coadd_maps(self.geom, cmaps)
          self._wcube = skymap.coadd_maps(self.geom, wmaps)
          self._ccube.write(self.files['ccube'], overwrite=True, conv='fgst-ccube')
  
          if self.projtype == "WCS":
              rm['counts'] += np.sum(self._ccube.data,
                                     axis=(1, 2), keepdims=False)
              rm['counts_wt'] += np.sum(self._ccube.data * self._wcube.data,
                                        axis=(1, 2), keepdims=False)
          elif self.projtype == "HPX":
              rm['counts'] += np.sum(self._ccube.data, axis=(1,), keepdims=False)
              rm['counts_wt'] += np.sum(self._ccube.data * self._wcube.data,
                                        axis=(1,), keepdims=False)
          else:
              raise Exception(
                  "Did not recognize projection type %s", self.projtype)
  
      def update_source(self, name, paramsonly=False, reoptimize=False, **kwargs):
          """Update the dictionary for this source.
  
          Parameters
          ----------
  
          name : str
  
          paramsonly : bool
  
          reoptimize : bool
             Re-fit background parameters in likelihood scan.
  
          """
  
          npts = self.config['gtlike']['llscan_npts']
          optimizer = kwargs.get('optimizer', self.config['optimizer'])
  
          sd = self.get_src_model(name, paramsonly, reoptimize, npts,
                                  optimizer=optimizer)
          src = self.roi.get_source_by_name(name)
          src.update_data(sd)
  
      def get_src_model(self, name, paramsonly=False, reoptimize=False,
                        npts=None, **kwargs):
          """Compose a dictionary for a source with the current best-fit
          parameters.
  
          Parameters
          ----------
  
          name : str
  
          paramsonly : bool
             Skip computing TS and likelihood profile.
  
          reoptimize : bool
             Re-fit background parameters in likelihood scan.
  
          npts : int
             Number of points for likelihood scan.
  
          Returns
          -------
          src_dict : dict
  
          """
  
          self.logger.debug('Generating source dict for ' + name)
  
          optimizer = kwargs.get('optimizer', self.config['optimizer'])
          if npts is None:
              npts = self.config['gtlike']['llscan_npts']
  
          name = self.get_source_name(name)
          source = self.like[name].src
          spectrum = source.spectrum()
          normPar = self.like.normPar(name)
  
          src_dict = defaults.make_default_dict(defaults.source_flux_output)
          src_dict.update({'name': name,
                           'pivot_energy': 1000.,
                           'ts': np.nan,
                           'loglike': np.nan,
                           'npred': 0.0,
                           'npred_wt': 0.0,
                           'loglike_scan': np.nan * np.ones(npts),
                           'dloglike_scan': np.nan * np.ones(npts),
                           'eflux_scan': np.nan * np.ones(npts),
                           'flux_scan': np.nan * np.ones(npts),
                           'norm_scan': np.nan * np.ones(npts),
                           })
  
          src_dict.update(gtutils.gtlike_spectrum_to_vectors(spectrum))
          src_dict['spectral_pars'] = gtutils.get_function_pars_dict(spectrum)
  
          # Get Counts Spectrum
          src_dict['model_counts'] = self.model_counts_spectrum(
              name, summed=True)
          src_dict['model_counts_wt'] = self.model_counts_spectrum(
              name, summed=True, weighted=True)
  
          # Get NPred
          src_dict['npred'] = self.like.NpredValue(str(name))
          # EAC, we need this b/c older version of the ST don't have the right signature
          try:
              src_dict['npred_wt'] = self.like.NpredValue(str(name), True)
          except (TypeError, NotImplementedError):
              src_dict['npred_wt'] = src_dict['npred']
  
          # Get the Model Fluxes
          try:
              thesrc = self.like[name]
              src_dict['flux'] = self.like.flux(name, self.energies[0],
                                                self.energies[-1])
              src_dict['flux100'] = self.like.flux(name, 100., 10 ** 5.5)
              src_dict['flux1000'] = self.like.flux(name, 1000., 10 ** 5.5)
              src_dict['flux10000'] = self.like.flux(name, 10000., 10 ** 5.5)
              src_dict['eflux'] = self.like.energyFlux(name,
                                                       self.energies[0],
                                                       self.energies[-1])
              src_dict['eflux100'] = self.like.energyFlux(name, 100.,
                                                          10 ** 5.5)
              src_dict['eflux1000'] = self.like.energyFlux(name, 1000.,
                                                           10 ** 5.5)
              src_dict['eflux10000'] = self.like.energyFlux(name, 10000.,
                                                            10 ** 5.5)
              src_dict['dnde'] = self.like[name].spectrum()(
                  pyLike.dArg(src_dict['pivot_energy']))
              src_dict['dnde100'] = self.like[name].spectrum()(
                  pyLike.dArg(100.))
              src_dict['dnde1000'] = self.like[name].spectrum()(
                  pyLike.dArg(1000.))
              src_dict['dnde10000'] = self.like[name].spectrum()(
                  pyLike.dArg(10000.))
              if normPar.getValue() == 0:
                  normPar.setValue(1.0)
  
                  dnde_index = -get_spectral_index(self.like[name],
                                                   src_dict['pivot_energy'])
  
                  dnde100_index = -get_spectral_index(self.like[name],
                                                      100.)
                  dnde1000_index = -get_spectral_index(self.like[name],
                                                       1000.)
                  dnde10000_index = -get_spectral_index(self.like[name],
                                                        10000.)
  
                  normPar.setValue(0.0)
              else:
                  dnde_index = -get_spectral_index(self.like[name],
                                                   src_dict['pivot_energy'])
  
                  dnde100_index = -get_spectral_index(self.like[name],
                                                      100.)
                  dnde1000_index = -get_spectral_index(self.like[name],
                                                       1000.)
                  dnde10000_index = -get_spectral_index(self.like[name],
                                                        10000.)
              src_dict['dnde_index'] = dnde_index
              src_dict['dnde100_index'] = dnde100_index
              src_dict['dnde1000_index'] = dnde1000_index
              src_dict['dnde10000_index'] = dnde10000_index
  
          except Exception:
              self.logger.error('Failed to update source parameters.',
                                exc_info=True)
  
          # Only compute TS, errors, and ULs if the source was free in
          # the fit
          if not self.get_free_source_params(name) or paramsonly:
              return src_dict
  
          emax = 10 ** 5.5
  
          try:
              src_dict['flux_err'] = self.like.fluxError(name,
                                                         self.energies[0],
                                                         self.energies[-1])
              src_dict['flux100_err'] = self.like.fluxError(name, 100., emax)
              src_dict['flux1000_err'] = self.like.fluxError(name, 1000., emax)
              src_dict['flux10000_err'] = self.like.fluxError(name, 10000., emax)
              src_dict['eflux_err'] = \
                  self.like.energyFluxError(name, self.energies[0],
                                            self.energies[-1])
              src_dict['eflux100_err'] = self.like.energyFluxError(name, 100.,
                                                                   emax)
              src_dict['eflux1000_err'] = self.like.energyFluxError(name, 1000.,
                                                                    emax)
              src_dict['eflux10000_err'] = self.like.energyFluxError(name, 10000.,
                                                                     emax)
  
          except Exception:
              pass
          # self.logger.error('Failed to update source parameters.',
          #  exc_info=True)
          lnlp = self.profile_norm(name, savestate=True,
                                   reoptimize=reoptimize, npts=npts,
                                   optimizer=optimizer)
  
          src_dict['loglike_scan'] = lnlp['loglike']
          src_dict['dloglike_scan'] = lnlp['dloglike']
          src_dict['eflux_scan'] = lnlp['eflux']
          src_dict['flux_scan'] = lnlp['flux']
          src_dict['norm_scan'] = lnlp['xvals']
          src_dict['loglike'] = np.max(lnlp['loglike'])
  
          flux_ul_data = utils.get_parameter_limits(
              lnlp['flux'], lnlp['dloglike'])
          eflux_ul_data = utils.get_parameter_limits(
              lnlp['eflux'], lnlp['dloglike'])
  
          if normPar.getValue() == 0:
              normPar.setValue(1.0)
              flux = self.like.flux(name, self.energies[0], self.energies[-1])
              flux100 = self.like.flux(name, 100., emax)
              flux1000 = self.like.flux(name, 1000., emax)
              flux10000 = self.like.flux(name, 10000., emax)
              eflux = self.like.energyFlux(name, self.energies[0],
                                           self.energies[-1])
              eflux100 = self.like.energyFlux(name, 100., emax)
              eflux1000 = self.like.energyFlux(name, 1000., emax)
              eflux10000 = self.like.energyFlux(name, 10000., emax)
  
              flux100_ratio = flux100 / flux
              flux1000_ratio = flux1000 / flux
              flux10000_ratio = flux10000 / flux
              eflux100_ratio = eflux100 / eflux
              eflux1000_ratio = eflux1000 / eflux
              eflux10000_ratio = eflux10000 / eflux
              normPar.setValue(0.0)
          else:
              flux100_ratio = src_dict['flux100'] / src_dict['flux']
              flux1000_ratio = src_dict['flux1000'] / src_dict['flux']
              flux10000_ratio = src_dict['flux10000'] / src_dict['flux']
  
              eflux100_ratio = src_dict['eflux100'] / src_dict['eflux']
              eflux1000_ratio = src_dict['eflux1000'] / src_dict['eflux']
              eflux10000_ratio = src_dict['eflux10000'] / src_dict['eflux']
  
          src_dict['flux_ul95'] = flux_ul_data['ul']
          src_dict['flux100_ul95'] = flux_ul_data['ul'] * flux100_ratio
          src_dict['flux1000_ul95'] = flux_ul_data['ul'] * flux1000_ratio
          src_dict['flux10000_ul95'] = flux_ul_data['ul'] * flux10000_ratio
  
          src_dict['eflux_ul95'] = eflux_ul_data['ul']
          src_dict['eflux100_ul95'] = eflux_ul_data['ul'] * eflux100_ratio
          src_dict['eflux1000_ul95'] = eflux_ul_data['ul'] * eflux1000_ratio
          src_dict['eflux10000_ul95'] = eflux_ul_data['ul'] * eflux10000_ratio
  
          # Extract covariance matrix
          fd = None
          try:
              fd = FluxDensity.FluxDensity(self.like, name)
              src_dict['covar'] = fd.covar
          except RuntimeError:
              pass
          # if ex.message == 'Covariance matrix has not been
          # computed.':
  
          # Extract bowtie
          if fd and len(src_dict['covar']) and src_dict['covar'].ndim >= 1:
              loge = np.linspace(self.log_energies[0],
                                 self.log_energies[-1], 50)
              src_dict['model_flux'] = self.bowtie(name, fd=fd, loge=loge)
              src_dict['dnde100_err'] = fd.error(100.)
              src_dict['dnde1000_err'] = fd.error(1000.)
              src_dict['dnde10000_err'] = fd.error(10000.)
  
              src_dict['pivot_energy'] = src_dict['model_flux']['pivot_energy']
  
              e0 = src_dict['pivot_energy']
              src_dict['dnde'] = self.like[name].spectrum()(pyLike.dArg(e0))
              src_dict['dnde_err'] = fd.error(e0)
  
          if not reoptimize:
              src_dict['ts'] = self.like.Ts2(name, reoptimize=reoptimize)
          else:
              src_dict['ts'] = -2.0 * lnlp['dloglike'][0]
  
          return src_dict
  
      def compute_srcprob(self,xmlfile=None, overwrite=False):
          """Run the gtsrcprob app with the current model or a user provided xmlfile"""
  
          for i,c in enumerate(self.components):
              # compute diffuse response, necessary for srcprob
              c._diffrsp_app(xmlfile=xmlfile)
              # compute srcprob
              c._srcprob_app(xmlfile = xmlfile, overwrite = overwrite)
  
  
  class GTBinnedAnalysis(fermipy.config.Configurable):
      defaults = dict(selection=defaults.selection,
                      binning=defaults.binning,
                      ltcube=defaults.ltcube,
                      gtlike=defaults.gtlike,
                      data=defaults.data,
                      model=defaults.model,
                      logging=defaults.logging,
                      fileio=defaults.fileio,
                      name=('00', '', str),
                      file_suffix=('', '', str))
  
      def __init__(self, config, roi, **kwargs):
  
          self._loglevel = kwargs.pop('loglevel', logging.INFO)
  
          super(GTBinnedAnalysis, self).__init__(config, **kwargs)
  
          self._projtype = self.config['binning']['projtype']
          self._logger_name = self.config['logging']['prefix'] + \
              self.__class__.__name__
  
          Logger.configure(self._logger_name, self.config['fileio']['logfile'],
                           log_level(self.config['logging']['verbosity']))
  
          self._roi = roi
          workdir = self.config['fileio']['workdir']
          self._name = self.config['name']
  
          search_dirs = [workdir]
  
          self._files = {}
  
          # Set default file names
          self._files['ft1'] = 'ft1%s.fits'
          self._files['ft1_filtered'] = 'ft1_filtered%s.fits'
          self._files['ccube'] = 'ccube%s.fits'
          self._files['ccubemc'] = 'ccubemc%s.fits'
          self._files['srcmap'] = 'srcmap%s.fits'
          self._files['bexpmap'] = 'bexpmap%s.fits'
          self._files['bexpmap_roi'] = 'bexpmap_roi%s.fits'
          self._files['srcmdl'] = 'srcmdl%s.xml'
          self._files['ltcube'] = 'ltcube%s.fits'
          self._files = {
              k: os.path.join(workdir, v % self.config['file_suffix'])
              for k, v in self._files.items()}
  
          # Override files defined in gtlike
          for k in ['srcmap', 'bexpmap', 'bexpmap_roi']:
  
              if self.config['gtlike'].get(k, None) is None:
                  continue
  
              self._files[k] = resolve_file_path(self.config['gtlike'][k],
                                                 search_dirs=search_dirs,
                                                 expand=True)
  
          self._data_files = {}
          self._data_files['evfile'] = self.config['data']['evfile']
          self._data_files['scfile'] = self.config['data']['scfile']
  
          for k, v in self._data_files.items():
              if v is None:
                  continue
              if not os.path.isfile(v):
                  v = os.path.join(workdir, v)
              if not os.path.isfile(v):
                  continue
              if not utils.is_fits_file(v):
                  self._data_files[k] = \
                      utils.resolve_file_path_list(v, workdir,
                                                   prefix=k + '_' + self.name)
              else:
                  self._data_files[k] = v
  
          self._srcmap_cache = {}
          self._srcmap = {}
  
          # Fill dictionary of exposure corrections
          self._src_expscale = {}
          if self.config['gtlike']['expscale'] is not None:
              self._src_expscale = {
                  src.name: self.config['gtlike']['expscale']
                  for src in self.roi}
  
          if self.config['gtlike']['src_expscale']:
              for k, v in self.config['gtlike']['src_expscale'].items():
                  self._src_expscale[k] = v
  
          # Setup external LT cube
          self._ext_ltcube = resolve_file_path(self.config['data']['ltcube'],
                                               search_dirs=search_dirs,
                                               expand=True)
  
          if self._ext_ltcube is not None:
              self.files['ltcube'] = self._ext_ltcube
  
          # Setup weights map
          self._files['wmap'] = resolve_file_path(self.config['gtlike']['wmap'],
                                                  search_dirs=search_dirs,
                                                  expand=True)
  
          try:
              emin = self.config['selection']['emin']
              emax = self.config['selection']['emax']
              logemin = np.log10(emin)
              logemax = np.log10(emax)
          except AttributeError:
              logemin = self.config['selection']['logemin']
              logemax = self.config['selection']['logemax']
              emin = np.power(10., logemin)
              emax = np.power(10., logemax)
  
          self.config['selection']['logemin'] = logemin
          self.config['selection']['logemax'] = logemax
          self.config['selection']['emin'] = emin
          self.config['selection']['emax'] = emax
  
          if self.config['binning']['enumbins'] is not None:
              self._enumbins = int(self.config['binning']['enumbins'])
          else:
              self._enumbins = np.round(self.config['binning']['binsperdec'] *
                                        np.log10(emax / emin))
              self._enumbins = int(self._enumbins)
  
          self._ebin_edges = np.linspace(logemin, logemax,
                                         self._enumbins + 1)
          self._ebin_center = 0.5 * \
              (self._ebin_edges[1:] + self._ebin_edges[:-1])
  
          if self.config['binning']['npix'] is None:
              roiwidth = (self.config['binning']['roiwidth']
                          if self.config['binning']['roiwidth'] is not None else 360.)
              self._npix = int(np.round(roiwidth /
                                        self.config['binning']['binsz']))
          else:
              self._npix = self.config['binning']['npix']
  
          if self.config['selection']['radius'] is None:
              self._config['selection']['radius'] = float(
                  np.sqrt(2.) * 0.5 * self.npix *
                  self.config['binning']['binsz'] + 0.5)
              self.logger.debug(
                  'Automatically setting selection radius to %s deg',
                  self.config['selection']['radius'])
  
          if self.config['binning']['coordsys'] == 'CEL':
              self._xref = float(self.roi.skydir.ra.deg)
              self._yref = float(self.roi.skydir.dec.deg)
          elif self.config['binning']['coordsys'] == 'GAL':
              self._xref = float(self.roi.skydir.galactic.l.deg)
              self._yref = float(self.roi.skydir.galactic.b.deg)
          else:
              raise Exception('Unrecognized coord system: ' +
                              self.config['binning']['coordsys'])
  
          self._like = None
          self._coordsys = self.config['binning']['coordsys']
          self._tmin = self.config['selection']['tmin']
          self._tmax = self.config['selection']['tmax']
          axes = [MapAxis.from_edges(self.energies, interp='log',
                                     name='energy', unit='MeV')]
          if self.projtype == 'HPX':
              is_nested = self.config['binning']['hpx_ordering_scheme'] == "NESTED"
              self._geom = HpxGeom.create(2**self.config['binning']['hpx_order'],
                                          nest=is_nested,
                                          coordsys=self.config['binning']['coordsys'],
                                          width=self.config['binning']['roiwidth'],
                                          skydir=self.roi.skydir,
                                          axes=axes)
          elif self.projtype == "WCS":
              self._geom = WcsGeom.create(npix=self.npix, binsz=self.binsz,
                                          coordsys=self.config['binning']['coordsys'],
                                          proj=self.config['binning']['proj'],
                                          skydir=self.roi.skydir,
                                          axes=axes)
          else:
              raise Exception(
                  "Did not recognize projection type %s", self.projtype)
  
          self.print_config(self.logger, loglevel=logging.DEBUG)
  
      @property
      def logger(self):
          """Return the default loglevel."""
          return logging.getLogger(self._logger_name)
  
      @property
      def loglevel(self):
          """Return the default loglevel."""
          return self._loglevel
  
      @property
      def workdir(self):
          """Return the analysis working directory."""
          return self.config['fileio']['workdir']
  
      @property
      def roi(self):
          return self._roi
  
      @property
      def like(self):
          return self._like
  
      @property
      def name(self):
          return self._name
  
      @property
      def energies(self):
          """Return the energy bin edges in MeV."""
          return 10**self._ebin_edges
  
      @property
      def log_energies(self):
          """Return the energy bin edges in log10(E/MeV)."""
          return self._ebin_edges
  
      @property
      def enumbins(self):
          return len(self._ebin_edges) - 1
  
      @property
      def npix(self):
          return self._npix
  
      @property
      def binsz(self):
          return self.config['binning']['binsz']
  
      @property
      def roiwidth(self):
          return self._npix * self.config['binning']['binsz']
  
      @property
      def projtype(self):
          """Return the type of projection to use"""
          return self._projtype
  
      @property
      def geom(self):
          """ROI geometry."""
          return self._geom
  
      @property
      def tmin(self):
          """Return the MET time for the start of the observation."""
          return self._tmin
  
      @property
      def tmax(self):
          """Return the MET time for the end of the observation."""
          return self._tmax
  
      @property
      def geom(self):
          return self._geom
  
      @property
      def wcs(self):
          if self.projtype == "WCS":
              return self.geom.wcs
          return None
  
      @property
      def hpx(self):
          if self.projtype == "HPX":
              return self.geom
          return None
  
      @property
      def coordsys(self):
          return self._coordsys
  
      @property
      def files(self):
          return self._files
  
      @property
      def data_files(self):
          return self._data_files
  
      @property
      def src_expscale(self):
          return self._src_expscale
  
      def reload_source(self, name):
          """Recompute the source map for a single source in the model.
          """
  
          src = self.roi.get_source_by_name(name)
  
          if hasattr(self.like.logLike, 'loadSourceMap'):
              self.like.logLike.loadSourceMap(str(name), True, False)
              srcmap_utils.delete_source_map(self.files['srcmap'], name)
              self.like.logLike.saveSourceMaps(str(self.files['srcmap']))
              self._scale_srcmap(self._src_expscale, check_header=False,
                                 names=[name])
              self.like.logLike.buildFixedModelWts()
          else:
              self.write_xml('tmp')
              src = self.delete_source(name)
              self.add_source(name, src, free=True)
              self.load_xml('tmp')
  
      def reload_sources(self, names):
          """Recompute the source map for a list of sources in the model.
          """
  
          try:
              self.like.logLike.loadSourceMaps(names, True, True)
              # loadSourceMaps doesn't overwrite the header so we need
              # to ignore EXPSCALE by setting check_header=False
              self._scale_srcmap(self._src_expscale, check_header=False,
                                 names=names)
          except:
              for name in names:
                  self.reload_source(name)
  
      def add_source(self, name, src_dict, free=None, save_source_maps=True,
                     use_pylike=True, use_single_psf=False):
          """Add a new source to the model.  Source properties
          (spectrum, spatial model) are set with the src_dict argument.
  
  
          Parameters
          ----------
  
          name : str
              Source name.
  
          src_dict : dict or `~fermipy.roi_model.Source` object
              Dictionary or Source object defining the properties of the
              source.
  
          free : bool
              Initialize the source with the normalization parameter free.
  
          save_source_maps : bool
              Write the source map for this source to the source maps file.
  
          use_pylike : bool
  
          use_single_psf : bool
  
          """
  
          # if self.roi.has_source(name):
          #    msg = 'Source %s already exists.' % name
          #    self.logger.error(msg)
          #    raise Exception(msg)
  
          srcmap_utils.delete_source_map(self.files['srcmap'], name)
  
          src = self.roi[name]
          if self.config['gtlike']['expscale'] is not None and \
                  name not in self._src_expscale:
              self._src_expscale[name] = self.config['gtlike']['expscale']
  
          if self._like is None:
              return
  
          if not use_pylike:
              self._update_srcmap_file([src], True)
  
          pylike_src = self._create_source(src)
  
          # Initialize source as free/fixed
          if free is not None:
              pylike_src.spectrum().normPar().setFree(free)
  
          if hasattr(pyLike, 'PsfIntegConfig') and \
                  hasattr(pyLike.PsfIntegConfig, 'set_use_single_psf'):
              config = pyLike.BinnedLikeConfig(self.like.logLike.config())
              config.psf_integ_config().set_use_single_psf(use_single_psf)
              self.like.addSource(pylike_src, config)
          else:
              self.like.addSource(pylike_src)
  
          self.like.syncSrcParams(str(name))
          self.like.logLike.buildFixedModelWts()
          if save_source_maps and \
                  not self.config['gtlike']['use_external_srcmap']:
              self.like.logLike.saveSourceMaps(str(self.files['srcmap']))
  
          self.set_exposure_scale(name)
  
      def _create_source(self, src):
          """Create a pyLikelihood Source object from a
          `~fermipy.roi_model.Model` object."""
  
          if src['SpatialType'] == 'SkyDirFunction':
              pylike_src = pyLike.PointSource(self.like.logLike.observation())
              pylike_src.setDir(src.skydir.ra.deg, src.skydir.dec.deg, False,
                                False)
          elif src['SpatialType'] == 'SpatialMap':
              filepath = str(utils.path_to_xmlpath(src['Spatial_Filename']))
              sm = pyLike.SpatialMap(filepath)
              pylike_src = pyLike.DiffuseSource(sm,
                                                self.like.logLike.observation(),
                                                False)
          elif src['SpatialType'] == 'RadialProfile':
              filepath = str(utils.path_to_xmlpath(src['Spatial_Filename']))
              sm = pyLike.RadialProfile(filepath)
              sm.setCenter(src['ra'], src['dec'])
              pylike_src = pyLike.DiffuseSource(sm,
                                                self.like.logLike.observation(),
                                                False)
          elif src['SpatialType'] == 'RadialGaussian':
              sm = pyLike.RadialGaussian(src.skydir.ra.deg, src.skydir.dec.deg,
                                         src.spatial_pars['Sigma']['value'])
              pylike_src = pyLike.DiffuseSource(sm,
                                                self.like.logLike.observation(),
                                                False)
  
          elif src['SpatialType'] == 'RadialDisk':
              sm = pyLike.RadialDisk(src.skydir.ra.deg, src.skydir.dec.deg,
                                     src.spatial_pars['Radius']['value'])
              pylike_src = pyLike.DiffuseSource(sm,
                                                self.like.logLike.observation(),
                                                False)
  
          elif src['SpatialType'] == 'MapCubeFunction':
              filepath = str(utils.path_to_xmlpath(src['Spatial_Filename']))
              mcf = pyLike.MapCubeFunction2(filepath)
              pylike_src = pyLike.DiffuseSource(mcf,
                                                self.like.logLike.observation(),
                                                False)
          else:
              raise Exception('Unrecognized spatial type: %s',
                              src['SpatialType'])
  
          if src['SpectrumType'] == 'FileFunction':
              fn = gtutils.create_spectrum_from_dict(src['SpectrumType'],
                                                     src.spectral_pars)
  
              file_function = pyLike.FileFunction_cast(fn)
              filename = str(os.path.expandvars(src['Spectrum_Filename']))
              file_function.readFunction(filename)
          elif src['SpectrumType'] == 'DMFitFunction':
  
              fn = pyLike.DMFitFunction()
              fn = gtutils.create_spectrum_from_dict(src['SpectrumType'],
                                                     src.spectral_pars, fn)
              filename = str(os.path.expandvars(src['Spectrum_Filename']))
              fn.readFunction(filename)
          else:
              fn = gtutils.create_spectrum_from_dict(src['SpectrumType'],
                                                     src.spectral_pars)
  
          pylike_src.setSpectrum(fn)
          pylike_src.setName(str(src.name))
  
          return pylike_src
  
      def delete_source(self, name, save_template=True, delete_source_map=False,
                        build_fixed_wts=True):
  
          src = self.roi.get_source_by_name(name)
  
          self.logger.debug('Deleting source %s', name)
  
          if self.like is not None:
  
              normPar = self.like.normPar(name)
              isFree = normPar.isFree()
  
              if str(src.name) in self.like.sourceNames():
                  self.like.deleteSource(str(src.name))
                  self.like.logLike.eraseSourceMap(str(src.name))
  
              if not isFree and build_fixed_wts:
                  self.like.logLike.buildFixedModelWts()
  
          if (not save_template and
              src['Spatial_Filename'] is not None and
              os.path.isfile(src['Spatial_Filename']) and
                  os.path.dirname(src['Spatial_Filename']) == self.config['fileio']['workdir']):
              os.remove(src['Spatial_Filename'])
  
          if delete_source_map:
              srcmap_utils.delete_source_map(self.files['srcmap'], name)
  
          return src
  
      def set_exposure_scale(self, name, scale=None):
          """Set the exposure correction of a source.
  
          Parameters
          ----------
          name : str
              Source name.
  
          scale : factor
              Exposure scale factor (1.0 = nominal exposure).
  
          """
          name = self.roi.get_source_by_name(name).name
          if scale is None and name not in self._src_expscale:
              return
          elif scale is None:
              scale = self._src_expscale.get(name, 1.0)
          else:
              self._src_expscale[name] = scale
          self._scale_srcmap({name: scale})
  
      def set_edisp_flag(self, name, flag=True):
          """Enable/Disable the energy dispersion correction for a
          source."""
          src = self.roi.get_source_by_name(name)
          name = src.name
          self.like[name].src.set_edisp_flag(flag)
  
      def set_energy_range(self, logemin, logemax):
          """Set the energy range of the analysis.
  
          Parameters
          ----------
          logemin: float
             Lower end of energy range in log10(E/MeV).
  
          logemax : float
             Upper end of energy range in log10(E/MeV).
  
          """
  
          if logemin is None:
              logemin = self.log_energies[0]
  
          if logemax is None:
              logemax = self.log_energies[-1]
  
          imin = int(utils.val_to_edge(self.log_energies, logemin)[0])
          imax = int(utils.val_to_edge(self.log_energies, logemax)[0])
  
          if imin - imax == 0:
              imin = int(len(self.log_energies) - 1)
              imax = int(len(self.log_energies) - 1)
  
          klims = self.like.logLike.klims()
          if imin != klims[0] or imax != klims[1]:
              self.like.selectEbounds(imin, imax)
  
          return np.array([self.log_energies[imin], self.log_energies[imax]])
  
      def counts_map(self):
          """Return 3-D counts map for this component as a Map object.
  
          Returns
          -------
          map : `~fermipy.skymap.MapBase`
  
          """
          try:
              if isinstance(self.like, gtutils.SummedLikelihood):
                  cmap = self.like.components[0].logLike.countsMap()
                  p_method = cmap.projection().method()
              else:
                  cmap = self.like.logLike.countsMap()
                  p_method = cmap.projection().method()
          except Exception:
              p_method = 0
  
          if p_method == 0:  # WCS
              z = cmap.data()
              z = np.array(z).reshape(self.enumbins, self.npix, self.npix)
              return WcsNDMap(copy.deepcopy(self.geom), z)
          elif p_method == 1:  # HPX
              z = cmap.data()
              z = np.array(z).reshape(self.enumbins, np.max(self.geom.npix))
              return HpxNDMap(copy.deepcopy(self.geom), z)
          else:
              self.logger.error('Did not recognize CountsMap type %i' % p_method,
                                exc_info=True)
          return None
  
      def weight_map(self):
          """Return 3-D weights map for this component as a Map object.
  
          Returns
          -------
          map : `~fermipy.skymap.MapBase`
  
          """
          # EAC we need the try blocks b/c older versions of the ST don't have some of these functions
          if isinstance(self.like, gtutils.SummedLikelihood):
              cmap = self.like.components[0].logLike.countsMap()
              try:
                  p_method = cmap.projection().method()
              except AttributeError:
                  p_method = 0
              try:
                  if self.like.components[0].logLike.has_weights():
                      wmap = self.like.components[0].logLike.weightMap()
                  else:
                      wmap = None
              except Exception:
                  wmap = None
          else:
              cmap = self.like.logLike.countsMap()
              try:
                  p_method = cmap.projection().method()
              except AttributeError:
                  p_method = 0
              try:
                  if self.like.logLike.has_weights():
                      wmap = self.like.logLike.weightMap()
                  else:
                      wmap = None
              except Exception:
                  wmap = None
  
          if p_method == 0:  # WCS
              if wmap is None:
                  z = np.ones((self.enumbins, self.npix, self.npix))
              else:
                  z = wmap.model()
                  z = np.array(z).reshape(self.enumbins, self.npix, self.npix)
              return WcsNDMap(copy.deepcopy(self._geom), z)
          elif p_method == 1:  # HPX
              nhpix = np.max(self.geom.npix)
              if wmap is None:
                  z = np.ones((self.enumbins, nhpix))
              else:
                  z = wmap.model()
                  z = np.array(z).reshape(self.enumbins, nhpix)
              return HpxNDMap(self.geom, z)
          else:
              self.logger.error('Did not recognize CountsMap type %i' % p_method,
                                exc_info=True)
          return None
  
      def model_counts_map(self, name=None, exclude=None, use_mask=False):
          """Return the model expectation map for a single source, a set
          of sources, or all sources in the ROI.  The map will be
          computed using the current model parameters.
  
          Parameters
          ----------
          name : str
  
             Parameter that defines the sources for which the model map
             will be calculated.  If name=None a model map will be
             generated for all sources in the model.  If name='diffuse'
             a map for all diffuse sources will be generated.
  
          exclude : list
  
             Source name or list of source names that will be excluded
             from the model map.
  
          use_mask : bool
  
             Parameter that specifies in the model counts map should include
             mask pixels (i.e., ones whose weights are <= 0)
  
          Returns
          -------
          map : `~fermipy.skymap.Map`
  
             A map object containing the counts and WCS projection.
  
          """
          if self.projtype == "WCS":
              v = pyLike.FloatVector(self.npix ** 2 * self.enumbins)
          elif self.projtype == "HPX":
              v = pyLike.FloatVector(np.max(self.geom.npix) * self.enumbins)
          else:
              raise Exception("Unknown projection type %s", self.projtype)
  
          exclude = utils.arg_to_list(exclude)
          names = utils.arg_to_list(name)
  
          excluded_names = []
          for i, t in enumerate(exclude):
              srcs = self.roi.get_sources_by_name(t)
              excluded_names += [s.name for s in srcs]
  
          if not hasattr(self.like.logLike, 'loadSourceMaps'):
              # Update fixed model
              self.like.logLike.buildFixedModelWts()
              # Populate source map hash
              self.like.logLike.buildFixedModelWts(True)
          elif (name is None or name == 'all') and not exclude:
              self.like.logLike.loadSourceMaps()
  
          src_names = []
          if (name is None) or (name == 'all'):
              src_names = [src.name for src in self.roi.sources]
          elif name == 'diffuse':
              src_names = [src.name for src in self.roi.sources if
                           src.diffuse]
          else:
              srcs = [self.roi.get_source_by_name(t) for t in names]
              src_names = [src.name for src in srcs]
  
          # Remove sources in exclude list
          src_names = [str(t) for t in src_names if t not in excluded_names]
  
          # EAC we need the try blocks b/c older versions of the ST don't have some of these functions
          if len(src_names) == len(self.roi.sources):
              try:
                  self.like.logLike.computeModelMap(v, use_mask)
              except (TypeError, NotImplementedError):
                  self.like.logLike.computeModelMap(v)
          elif not hasattr(self.like.logLike, 'setSourceMapImage'):
              for s in src_names:
                  model = self.like.logLike.sourceMap(str(s))
                  try:
                      self.like.logLike.updateModelMap(v, model, use_mask)
                  except (TypeError, NotImplementedError):
                      self.like.logLike.updateModelMap(v, model)
          else:
              try:
                  if hasattr(self.like.logLike, 'has_weights'):
                      self.like.logLike.computeModelMap(src_names, v, use_mask)
                  else:
                      self.like.logLike.computeModelMap(src_names, v)
              except:
                  vsum = np.zeros(v.size())
                  for s in src_names:
                      vtmp = pyLike.FloatVector(v.size())
                      if hasattr(self.like.logLike, 'has_weights'):
                          self.like.logLike.computeModelMap(
                              str(s), vtmp, use_mask)
                      else:
                          self.like.logLike.computeModelMap(str(s), vtmp)
                      vsum += vtmp
                  v = pyLike.FloatVector(vsum)
  
          if self.projtype == "WCS":
              z = np.array(v).reshape(self.enumbins, self.npix, self.npix)
              return WcsNDMap(copy.deepcopy(self._geom), z)
          elif self.projtype == "HPX":
              z = np.array(v).reshape(self.enumbins, np.max(self._geom.npix))
              return HpxNDMap(copy.deepcopy(self._geom), z)
          else:
              raise Exception(
                  "Did not recognize projection type %s", self.projtype)
  
      def model_counts_spectrum(self, name, logemin, logemax, weighted=False):
          """Return the model counts spectrum of a source.
  
          Parameters
          ----------
          name : str
             Source name.
  
          """
          # EAC, we need this b/c older version of the ST don't have the right signature
          try:
              cs = np.array(self.like.logLike.modelCountsSpectrum(
                  str(name), weighted))
          except (TypeError, NotImplementedError):
              cs = np.array(self.like.logLike.modelCountsSpectrum(str(name)))
          imin = utils.val_to_edge(self.log_energies, logemin)[0]
          imax = utils.val_to_edge(self.log_energies, logemax)[0]
          if imax <= imin:
              raise Exception('Invalid energy range.')
          return cs[imin:imax]
  
      def setup(self, overwrite=False, **kwargs):
          """Run pre-processing step for this component.  This will
          generate all of the auxiliary files needed to instantiate a
          likelihood object.  By default this function will skip any
          steps for which the output file already exists.
  
          Parameters
          ----------
          overwrite : bool
  
              Run all pre-processing steps even if the output file of
              that step is present in the working directory.
          """
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          self.logger.log(loglevel, 'Running setup for component %s',
                          self.name)
  
          use_external_srcmap = self.config['gtlike']['use_external_srcmap']
  
          # Run data selection
          if not use_external_srcmap:
              self._select_data(overwrite=overwrite, **kwargs)
  
          # Create LT Cube
          if self._ext_ltcube is not None:
              self.logger.log(loglevel, 'Using external LT cube.')
          else:
              self._create_ltcube(overwrite=overwrite, **kwargs)
  
          self.logger.debug('Loading LT Cube %s', self.files['ltcube'])
          self._ltc = LTCube.create(self.files['ltcube'])
  
          # Extract tmin, tmax from LT cube
          self._tmin = self._ltc.tstart
          self._tmax = self._ltc.tstop
  
          self.logger.debug('Creating PSF model')
          self._psf = irfs.PSFModel.create(self.roi.skydir, self._ltc,
                                           self.config['gtlike']['irfs'],
                                           self.config['selection']['evtype'],
                                           self.energies)
  
          # Bin data and create exposure cube
          if not use_external_srcmap:
              self._bin_data(overwrite=overwrite, **kwargs)
              self._create_expcube(overwrite=overwrite, **kwargs)
  
          # This is needed in case the exposure map is in HEALPix
          hpxhduname = "HPXEXPOSURES"
          try:
              self._bexp = Map.read(self.files['bexpmap'], hdu=hpxhduname)
          except KeyError:
              self._bexp = Map.read(self.files['bexpmap'])
  
          # Write ROI XML
          self.roi.write_xml(self.files['srcmdl'], self.config['model'])
  
          # Create source maps file
          if not use_external_srcmap:
              self._create_srcmaps(overwrite=overwrite)
  
          if not self.config['data']['cacheft1'] and os.path.isfile(self.files['ft1']):
              self.logger.debug('Deleting FT1 file.')
              os.remove(self.files['ft1'])
  
          self.logger.log(loglevel, 'Finished setup for component %s',
                          self.name)
  
      def _select_data(self, overwrite=False, **kwargs):
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          if (os.path.isfile(self.files['ft1']) and
              os.path.isfile(self.files['ccube']) and
                  not overwrite):
              self.logger.log(loglevel, 'Skipping data selection.')
              return
  
          # Run gtselect and gtmktime
          kw_gtselect = dict(infile=self.data_files['evfile'],
                             outfile=self.files['ft1'],
                             ra=self.roi.skydir.ra.deg,
                             dec=self.roi.skydir.dec.deg,
                             rad=self.config['selection']['radius'],
                             convtype=self.config['selection']['convtype'],
                             phasemin=self.config['selection']['phasemin'],
                             phasemax=self.config['selection']['phasemax'],
                             evtype=self.config['selection']['evtype'],
                             evclass=self.config['selection']['evclass'],
                             tmin=self.config['selection']['tmin'],
                             tmax=self.config['selection']['tmax'],
                             emin=self.config['selection']['emin'],
                             emax=self.config['selection']['emax'],
                             zmax=self.config['selection']['zmax'],
                             chatter=self.config['logging']['chatter'])
  
          kw_gtmktime = dict(evfile=self.files['ft1'],
                             outfile=self.files['ft1_filtered'],
                             scfile=self.data_files['scfile'],
                             roicut=self.config['selection']['roicut'],
                             filter=self.config['selection']['filter'])
  
          run_gtapp('gtselect', self.logger, kw_gtselect,
                    loglevel=loglevel)
          if self.config['selection']['roicut'] == 'yes' or \
                  self.config['selection']['filter'] is not None:
              run_gtapp('gtmktime', self.logger, kw_gtmktime,
                        loglevel=loglevel)
              os.system('mv %s %s' % (self.files['ft1_filtered'],
                                      self.files['ft1']))
  
      def _bin_data(self, overwrite=False, **kwargs):
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          # Run gtbin
          if self.projtype == "WCS":
              kw = dict(algorithm='ccube',
                        nxpix=self.npix, nypix=self.npix,
                        binsz=self.config['binning']['binsz'],
                        evfile=self.files['ft1'],
                        outfile=self.files['ccube'],
                        scfile=self.data_files['scfile'],
                        xref=self._xref,
                        yref=self._yref,
                        axisrot=0,
                        proj=self.config['binning']['proj'],
                        ebinalg='LOG',
                        emin=self.config['selection']['emin'],
                        emax=self.config['selection']['emax'],
                        enumbins=self._enumbins,
                        coordsys=self.config['binning']['coordsys'],
                        chatter=self.config['logging']['chatter'])
          elif self.projtype == "HPX":
              if self.config['binning']['roiwidth'] is None or self.config['binning']['roiwidth'] >= 180.:
                  hpx_region = None
              else:
                  hpx_region = "DISK(%.3f,%.3f,%.3f)" % (
                      self._xref, self._yref, 0.5 * self.config['binning']['roiwidth'])
              kw = dict(algorithm='healpix',
                        evfile=self.files['ft1'],
                        outfile=self.files['ccube'],
                        scfile=self.data_files['scfile'],
                        xref=self._xref,
                        yref=self._yref,
                        proj=self.config['binning']['proj'],
                        hpx_ordering_scheme=self.config['binning'][
                            'hpx_ordering_scheme'],
                        hpx_order=self.config['binning']['hpx_order'],
                        hpx_ebin=self.config['binning']['hpx_ebin'],
                        hpx_region=hpx_region,
                        ebinalg='LOG',
                        emin=self.config['selection']['emin'],
                        emax=self.config['selection']['emax'],
                        enumbins=self._enumbins,
                        coordsys=self.config['binning']['coordsys'],
                        chatter=self.config['logging']['chatter'])
          else:
              self.logger.error(
                  'Unknown projection type, %s. Choices are WCS or HPX',
                  self.projtype)
  
          if not os.path.isfile(self.files['ccube']) or overwrite:
              run_gtapp('gtbin', self.logger, kw, loglevel=loglevel)
          else:
              self.logger.debug('Skipping gtbin.')
  
      def _create_ltcube(self, overwrite=False, **kwargs):
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          if os.path.isfile(self.files['ltcube']) and not overwrite:
              self.logger.log(loglevel, 'Skipping LT Cube.')
              return
  
          # Run gtltcube
          kw = dict(evfile=self.files['ft1'],
                    scfile=self.data_files['scfile'],
                    outfile=self.files['ltcube'],
                    binsz=self.config['ltcube']['binsz'],
                    dcostheta=self.config['ltcube']['dcostheta'],
                    phibins=self.config['ltcube']['phibins'],
                    zmax=self.config['selection']['zmax'])
  
          if self.config['ltcube']['use_local_ltcube']:
              self.logger.info('Generating local LT cube.')
              colnames = ['START', 'STOP', 'LIVETIME',
                          'RA_SCZ', 'DEC_SCZ',
                          'RA_ZENITH', 'DEC_ZENITH']
              tab_sc = create_sc_table(self.data_files['scfile'],
                                       colnames=colnames)
              tab_gti = Table.read(self.files['ft1'], 'GTI')
              radius = self.config['selection']['radius'] + 10.0
              ltc_new = LTCube.create_from_gti(self.roi.skydir, tab_sc, tab_gti,
                                               self.config['selection']['zmax'],
                                               radius=radius)
              ltc_new.write(self.files['ltcube'])
          else:
              run_gtapp('gtltcube', self.logger, kw, loglevel=loglevel)
  
      def _create_expcube(self, overwrite=False, **kwargs):
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          if os.path.isfile(self.files['bexpmap']) and not overwrite:
              self.logger.log(loglevel, 'Skipping gtexpcube.')
              return
  
          if self.config['gtlike']['irfs'] == 'CALDB':
              if self.projtype == "HPX":
                  cmap = None
              else:
                  cmap = self.files['ccube']
          else:
              cmap = 'none'
  
          # Run gtexpcube2
          kw = dict(infile=self.files['ltcube'], cmap=cmap,
                    ebinalg='LOG',
                    emin=self.config['selection']['emin'],
                    emax=self.config['selection']['emax'],
                    enumbins=self._enumbins,
                    outfile=self.files['bexpmap'], proj='CAR',
                    nxpix=360, nypix=180, binsz=1,
                    xref=0.0, yref=0.0,
                    evtype=self.config['selection']['evtype'],
                    irfs=self.config['gtlike']['irfs'],
                    coordsys=self.config['binning']['coordsys'],
                    chatter=self.config['logging']['chatter'])
  
          set_edisp_kwargs('gtexpcube2', self.config, kw)
          
          run_gtapp('gtexpcube2', self.logger, kw, loglevel=loglevel)
  
          if self.projtype == "WCS":
              kw = dict(infile=self.files['ltcube'], cmap='none',
                        ebinalg='LOG',
                        emin=self.config['selection']['emin'],
                        emax=self.config['selection']['emax'],
                        enumbins=self._enumbins,
                        outfile=self.files['bexpmap_roi'], proj='CAR',
                        nxpix=self.npix, nypix=self.npix,
                        binsz=self.config['binning']['binsz'],
                        xref=self._xref, yref=self._yref,
                        evtype=self.config['selection']['evtype'],
                        irfs=self.config['gtlike']['irfs'],
                        coordsys=self.config['binning']['coordsys'],
                        chatter=self.config['logging']['chatter'])
  
              run_gtapp('gtexpcube2', self.logger, kw, loglevel=loglevel)
          elif self.projtype == "HPX":
              self.logger.debug('Skipping local gtexpcube for HEALPix')
          else:
              raise Exception(
                  "Did not recognize projection type %s", self.projtype)
  
      def _create_srcmaps(self, overwrite=False, **kwargs):
  
          loglevel = kwargs.get('loglevel', self.loglevel)
          use_scaled_srcmap = self.config['gtlike']['use_scaled_srcmap']
  
          # Run gtsrcmaps
          kw = dict(scfile=self.data_files['scfile'],
                    expcube=self.files['ltcube'],
                    cmap=self.files['ccube'],
                    srcmdl=self.files['srcmdl'],
                    bexpmap=self.files['bexpmap'],
                    outfile=self.files['srcmap'],
                    irfs=self.config['gtlike']['irfs'],
                    wmap=self.files['wmap'],
                    evtype=self.config['selection']['evtype'],
                    rfactor=self.config['gtlike']['rfactor'],
                    #                   resample=self.config['resample'],
                    minbinsz=self.config['gtlike']['minbinsz'],
                    chatter=self.config['logging']['chatter'],
                    emapbnds='no')
  
          set_edisp_kwargs('gtsrcmaps', self.config, kw)
  
          if os.path.isfile(self.files['srcmap']) and not overwrite:
              self.logger.log(loglevel, 'Skipping gtsrcmaps.')
          elif use_scaled_srcmap:
              make_scaled_srcmap(self.roi,
                                 self.config['gtlike']['srcmap_base'],
                                 self.config['gtlike']['bexpmap_base'],
                                 self.config['gtlike']['bexpmap_roi_base'],
                                 self.files['bexpmap'],
                                 self.files['bexpmap_roi'],
                                 self.files['ccube'],
                                 self.files['srcmap'])
          else:
              run_gtapp('gtsrcmaps', self.logger, kw, loglevel=loglevel)
  
      def _create_binned_analysis(self, xmlfile=None, **kwargs):
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          self.logger.log(loglevel, 'Creating BinnedAnalysis for component %s.',
                          self.name)
  
          srcmdl_file = self.files['srcmdl']
          if xmlfile is not None:
              srcmdl_file = self.get_model_path(xmlfile)
  
          # Create BinnedObs
          self.logger.debug('Creating BinnedObs')
          kw = dict(srcMaps=self.files['srcmap'], expCube=self.files['ltcube'],
                    binnedExpMap=self.files['bexpmap'],
                    irfs=self.config['gtlike']['irfs'])
          self.logger.debug(kw)
  
          self._obs = ba.BinnedObs(**utils.unicode_to_str(kw))
  
  
          # Create BinnedAnalysis
          self.logger.debug('Creating BinnedAnalysis')
          kw = dict(srcModel=srcmdl_file,
                    optimizer='MINUIT',
                    wmap=self._files['wmap'],
                    convolve=self.config['gtlike']['convolve'],
                    resample=self.config['gtlike']['resample'],
                    minbinsz=self.config['gtlike']['minbinsz'],
                    resamp_fact=self.config['gtlike']['rfactor'])
  
          set_edisp_kwargs('gtlike', self.config, kw)
          self.logger.debug(kw)
  
          self._like = BinnedAnalysis(binnedData=self._obs,
                                      **utils.unicode_to_str(kw))
  
  #        print(self.like.logLike.use_single_fixed_map())
  #        self.like.logLike.set_use_single_fixed_map(False)
  #        print(self.like.logLike.use_single_fixed_map())
  
          for s in self.config['gtlike']['edisp_disable']:
  
              if not self.roi.has_source(s):
                  continue
  
              self.logger.debug('Disabling energy dispersion for %s', s)
              self.set_edisp_flag(s, False)
  
          if hasattr(self._like.logLike, 'set_save_all_srcmaps'):
              self._like.logLike.set_save_all_srcmaps(True)
  
          # Recompute fixed model weights
          self.logger.debug('Computing fixed weights')
          self.like.logLike.buildFixedModelWts()
          self.logger.debug('Updating source maps')
          if not self.config['gtlike']['use_external_srcmap']:
              self.like.logLike.saveSourceMaps(str(self.files['srcmap']))
  
          # Apply exposure corrections
          self._scale_srcmap(self._src_expscale)
  
      def _scale_srcmap(self, scale_map, check_header=True, names=None):
          """Apply exposure corrections to the source map file.
  
          Parameters
          ----------
          scale_map : dict
              Dictionary of exposure corrections.
  
          check_header : bool
              Check EXPSCALE header keyword to see if an exposure
              correction has already been applied to this source.
  
          names : list, optional
              Names of sources to which the exposure correction will be
              applied.  If None then all sources will be corrected.
          """
  
          srcmap = fits.open(self.files['srcmap'])
  
          for hdu in srcmap[1:]:
              if hdu.name not in scale_map:
                  continue
              if names is not None and hdu.name not in names:
                  continue
  
              scale = scale_map[hdu.name]
              if scale < 1e-20:
                  self.logger.warning(
                      "The expscale parameter was zero, setting it to 1e-8")
                  scale = 1e-8
              if 'EXPSCALE' in hdu.header and check_header:
                  old_scale = hdu.header['EXPSCALE']
              else:
                  old_scale = 1.0
              hdu.data *= scale / old_scale
              hdu.header['EXPSCALE'] = (scale,
                                        'Exposure correction applied to this map')
  
          srcmap.writeto(self.files['srcmap'], overwrite=True)
          srcmap.close()
  
          # Force reloading the map from disk
          for name in scale_map.keys():
              self.like.logLike.eraseSourceMap(str(name))
          self.like.logLike.buildFixedModelWts()
  
      def _make_scaled_srcmap(self):
          """Make an exposure cube with the same binning as the counts map."""
  
          self.logger.info('Computing scaled source map.')
  
          bexp0 = fits.open(self.files['bexpmap_roi'])
          bexp1 = fits.open(self.config['gtlike']['bexpmap'])
          srcmap = fits.open(self.config['gtlike']['srcmap'])
  
          if bexp0[0].data.shape != bexp1[0].data.shape:
              raise Exception('Wrong shape for input exposure map file.')
  
          bexp_ratio = bexp0[0].data / bexp1[0].data
  
          self.logger.info(
              'Min/Med/Max exposure correction: %f %f %f' % (np.min(bexp_ratio),
                                                             np.median(
                                                                 bexp_ratio),
                                                             np.max(bexp_ratio)))
  
          for hdu in srcmap[1:]:
  
              if hdu.name == 'GTI':
                  continue
              if hdu.name == 'EBOUNDS':
                  continue
              hdu.data *= bexp_ratio
  
          srcmap.writeto(self.files['srcmap'], overwrite=True)
  
      def restore_counts_maps(self):
  
          cmap = Map.read(self.files['ccube'])
  
          if hasattr(self.like.logLike, 'setCountsMap'):
              self.like.logLike.setCountsMap(np.ravel(cmap.data.astype(float)))
  
          srcmap_utils.update_source_maps(self.files['srcmap'],
                                          {'PRIMARY': cmap.data},
                                          logger=self.logger)
  
      def simulate_roi(self, name=None, clear=True, randomize=True):
          """Simulate the whole ROI or inject a simulation of one or
          more model components into the data.
  
          Parameters
          ----------
          name : str
             Name of the model component to be simulated.  If None then
             the whole ROI will be simulated.
  
          clear : bool
             Zero the current counts map before injecting the simulation.
  
          randomize : bool
             Fill with each pixel with random values drawn from a
             poisson distribution.  If false then fill each pixel with
             the counts expectation value.
  
          """
          cm = self.counts_map()
          data = cm.data
          m = self.model_counts_map(name) 
  
          if clear:
              data.fill(0.0)
  
          if randomize:
              if m.data.min()<0.:
                  self.logger.warning('At least on negative value found in model map.'
                                      ' Changing it/them to 0')
                  indexcond = np.where( m.data <0. )
                  m.data[indexcond]=np.zeros(len(m.data[indexcond]))
                  
              data += np.random.poisson(m.data).astype(float)
          else:
              data += m.data
  
          if hasattr(self.like.logLike, 'setCountsMap'):
              self.like.logLike.setCountsMap(np.ravel(data))
  
          srcmap_utils.update_source_maps(self.files['srcmap'],
                                          {'PRIMARY': data},
                                          logger=self.logger)
          cm.write(self.files['ccubemc'], overwrite=True, conv='fgst-ccube')
  
      def write_model_map(self, model_name=None, name=None):
          """Save counts model map to a FITS file.
  
          """
  
          if model_name is None:
              suffix = self.config['file_suffix']
          else:
              suffix = '_%s%s' % (model_name, self.config['file_suffix'])
  
          self.logger.info('Generating model map for component %s.', self.name)
  
          outfile = os.path.join(self.config['fileio']['workdir'],
                                 'mcube%s.fits' % (suffix))
  
          cmap = self.model_counts_map(name, use_mask=False)
          cmap.write(outfile, overwrite=True, conv='fgst-ccube')
          return cmap
  
      def write_weight_map(self, model_name=None):
          """Save counts model map to a FITS file.
  
          """
  
          if model_name is None:
              suffix = self.config['file_suffix']
          else:
              suffix = '_%s%s' % (model_name, self.config['file_suffix'])
  
          self.logger.info('Generating model map for component %s.', self.name)
  
          outfile = os.path.join(self.config['fileio']['workdir'],
                                 'wcube%s.fits' % (suffix))
  
          wmap = self.weight_map()
          wmap.write(outfile, overwrite=True, conv='fgst-ccube')
          return wmap
  
      def _update_srcmap_file(self, sources, overwrite=True):
          """Check the contents of the source map file and generate
          source maps for any components that are not present."""
  
          if not os.path.isfile(self.files['srcmap']):
              return
  
          hdulist = fits.open(self.files['srcmap'])
          hdunames = [hdu.name.upper() for hdu in hdulist]
  
          srcmaps = {}
  
          for src in sources:
  
              if src.name.upper() in hdunames and not overwrite:
                  continue
              self.logger.debug('Creating source map for %s', src.name)
              srcmaps[src.name] = self._create_srcmap(src.name, src)
  
          if srcmaps:
              self.logger.debug(
                  'Updating source map file for component %s.', self.name)
              srcmap_utils.update_source_maps(self.files['srcmap'], srcmaps,
                                              logger=self.logger)
          hdulist.close()
  
      def _create_srcmap_cache(self, name, src, **kwargs):
  
          from fermipy.srcmap_utils import SourceMapCache
  
          skydir = src.skydir
          spatial_model = src['SpatialModel']
          spatial_width = src['SpatialWidth']
          xpix, ypix = self.geom.to_image().coord_to_pix(skydir)
          exp = self._bexp.interp_by_coord(
              (skydir, self._bexp.geom.axes[0].center))
          rebin = min(int(np.ceil(self.binsz / 0.01)), 8)
          shape_out = (self.enumbins + 1, self.npix, self.npix)
          cache = SourceMapCache.create(self._psf, exp, spatial_model,
                                        spatial_width, shape_out,
                                        self.config['binning']['binsz'],
                                        rebin=rebin)
          self._srcmap_cache[name] = cache
  
      def _create_srcmap(self, name, src, **kwargs):
          """Generate the source map for a source."""
  
          psf_scale_fn = kwargs.get('psf_scale_fn', None)
          skydir = src.skydir
          spatial_model = src['SpatialModel']
          spatial_width = src['SpatialWidth']
          xpix, ypix = self.geom.to_image().coord_to_pix(skydir)
          exp = self._bexp.interp_by_coord(
              (skydir, self._bexp.geom.axes[0].center))
          cache = self._srcmap_cache.get(name, None)
          if cache is not None:
              k = cache.create_map([ypix, xpix])
          else:
              k = srcmap_utils.make_srcmap(self._psf, exp, spatial_model,
                                           spatial_width,
                                           npix=self.npix, xpix=xpix, ypix=ypix,
                                           cdelt=self.config['binning']['binsz'],
                                           psf_scale_fn=psf_scale_fn,
                                           sparse=True)
  
          return k
  
      def _update_srcmap(self, name, src, **kwargs):
          """Update the source map for an existing source in memory."""
  
          k = self._create_srcmap(name, src, **kwargs)
          scale = self._src_expscale.get(name, 1.0)
          k *= scale
  
          # Force the source map to be cached
          # FIXME: No longer necessary to force cacheing in ST after 11-05-02
          self.like.logLike.sourceMap(str(name)).model()
          self.like.logLike.setSourceMapImage(str(name), np.ravel(k))
          self.like.logLike.sourceMap(str(name)).model()
  
          normPar = self.like.normPar(name)        
          if not normPar.isFree():
              self.like.logLike.buildFixedModelWts()
  
      def generate_model(self, model_name=None, outfile=None):
          """Generate a counts model map from an XML model file using
          gtmodel.
  
          Parameters
          ----------
          model_name : str
              Name of the model.  If no name is given it will use the
              baseline model.
  
          outfile : str
              Override the name of the output model file.
          """
  
          if model_name is not None:
              model_name = os.path.splitext(model_name)[0]
  
          if model_name is None or model_name == '':
              srcmdl = self.files['srcmdl']
          else:
              srcmdl = self.get_model_path(model_name)
  
          if not os.path.isfile(srcmdl):
              raise Exception("Model file does not exist: %s", srcmdl)
  
          if model_name is None:
              suffix = self.config['file_suffix']
          else:
              suffix = '_%s%s' % (model_name, self.config['file_suffix'])
  
          outfile = os.path.join(self.config['fileio']['workdir'],
                                 'mcube%s.fits' % (suffix))
  
          # May consider generating a custom source model file
          if not os.path.isfile(outfile):
  
              kw = dict(srcmaps=self.files['srcmap'],
                        srcmdl=srcmdl,
                        bexpmap=self.files['bexpmap'],
                        outfile=outfile,
                        expcube=self.files['ltcube'],
                        irfs=self.config['gtlike']['irfs'],
                        evtype=self.config['selection']['evtype'],
                        outtype='ccube',
  
                        chatter=self.config['logging']['chatter'])
  
              set_edisp_kwargs('gtmodel', self.config, kw)
              run_gtapp('gtmodel', self.logger, kw)
          else:
              self.logger.info('Skipping gtmodel')
  
      def load_xml(self, xmlfile):
  
          xmlfile = self.get_model_path(xmlfile)
          self.logger.info('Loading %s' % xmlfile)
          self.like.logLike.reReadXml(str(xmlfile))
          if not self.like.logLike.fixedModelUpdated():
              self.like.logLike.buildFixedModelWts()
  
      def write_xml(self, xmlfile):
          """Write the XML model for this analysis component."""
  
          xmlfile = self.get_model_path(xmlfile)
          self.logger.info('Writing %s...', xmlfile)
          self.like.writeXml(str(xmlfile))
  
      def get_model_path(self, name):
          """Infer the path to the XML model name."""
  
          name, ext = os.path.splitext(name)
          ext = '.xml'
          xmlfile = name + self.config['file_suffix'] + ext
          xmlfile = utils.resolve_path(xmlfile,
                                       workdir=self.config['fileio']['workdir'])
  
          return xmlfile
  
      def _tscube_app(self, xmlfile):
          """Run gttscube as an application."""
  
          xmlfile = self.get_model_path(xmlfile)
  
          outfile = os.path.join(self.config['fileio']['workdir'],
                                 'tscube%s.fits' % (self.config['file_suffix']))
  
          kw = dict(cmap=self.files['ccube'],
                    expcube=self.files['ltcube'],
                    bexpmap=self.files['bexpmap'],
                    irfs=self.config['gtlike']['irfs'],
                    evtype=self.config['selection']['evtype'],
                    srcmdl=xmlfile,
                    nxpix=self.npix, nypix=self.npix,
                    binsz=self.config['binning']['binsz'],
                    xref=float(self.roi.skydir.ra.deg),
                    yref=float(self.roi.skydir.dec.deg),
                    proj=self.config['binning']['proj'],
                    stlevel=0,
                    coordsys=self.config['binning']['coordsys'],
                    outfile=outfile)
  
          run_gtapp('gttscube', self.logger, kw)
  
      def _diffrsp_app(self,xmlfile=None, **kwargs):
          """
          Compute the diffuse response 
          """
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          self.logger.log(loglevel, 'Computing diffuse repsonce for component %s.',
                          self.name)
  
          # set the srcmdl
          srcmdl_file = self.files['srcmdl']
          if xmlfile is not None:
              srcmdl_file = self.get_model_path(xmlfile)
  
          kw = dict(evfile=self.files['ft1'],
              scfile=self.data_files['scfile'],
              irfs = self.config['gtlike']['irfs'],
              evtype = self.config['selection']['evtype'],
              srcmdl = srcmdl_file)
  
          run_gtapp('gtdiffrsp', self.logger, kw, loglevel=loglevel)
          return 
  
      def _srcprob_app(self,xmlfile=None, overwrite=False, **kwargs):
          """
          Run srcprob for an analysis component as an application
          """
  
          loglevel = kwargs.get('loglevel', self.loglevel)
  
          self.logger.log(loglevel, 'Computing src probability for component %s.',
                          self.name)
  
          # set the srcmdl
          srcmdl_file = self.files['srcmdl']
          if xmlfile is not None:
              srcmdl_file = self.get_model_path(xmlfile)
  
          # set the outfile
          # it's defined here and not in self.files dict
          # so that it is copied with the stage_output module
          # even if savefits is False
          outfile = os.path.join(self.workdir, 
                      'ft1_srcprob{0[file_suffix]:s}.fits'.format(self.config))
  
          kw = dict(evfile=self.files['ft1'],
              scfile=self.data_files['scfile'],
              outfile= outfile,
              irfs = self.config['gtlike']['irfs'],
              srcmdl = srcmdl_file)
          self.logger.debug(kw)
  
          # run gtapp for the srcprob
          if os.path.isfile(outfile) and not overwrite:
              self.logger.info('Skipping gtsrcprob')
          else:
              run_gtapp('gtsrcprob', self.logger, kw, loglevel=loglevel)
  