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