def mwl_fit_low_level():
    """Use high-level Sherpa API.

    Low-level = no session, classes.

    Example: http://python4astronomers.github.io/fitting/low-level.html
    """
    fermi_data = FermiData().sherpa_data
    hess_data = IACTData().sherpa_data

    # spec_model = PowLaw1D('spec_model')
    spec_model = LogParabola('spec_model')
    spec_model.c1 = 0.5
    spec_model.c2 = 0.2
    spec_model.ampl = 5e-11

    data = DataSimulFit(name='global_data', datasets=[fermi_data, hess_data])
    # TODO: Figure out how to notice using the low-level API
    # data.notice(mins=1e-3, maxes=None, axislist=None)
    model = SimulFitModel(name='global_model', parts=[spec_model, spec_model])
    stat = FermiStat()
    method = LevMar()
    fit = Fit(data=data, model=model, stat=stat, method=method)
    result = fit.fit()

    # IPython.embed()
    return Bunch(results=result, model=spec_model)

    def fit(self, current):
        self._fit.model.thawedpars = current

        # nm = NelderMead()
        # nm.config['iquad'] = 0
        # nm.config['finalsimplex'] = 1

        #lm = LevMar()
        #lm.config['maxfev'] = 5

        cv = Covariance()

        # Use the fit method defined before called get_draws().  This way the user
        # does not have to pass in the fitting method and method options.
        fit = Fit(self._fit.data, self._fit.model, self._fit.stat, self._fit.method,
                  cv)

        fit_result = fit.fit()
        covar_result = fit.est_errors()
        sigma = np.array(covar_result.extra_output)

        if np.isnan(sigma).any():
            raise CovarError("NaNs found in covariance matrix")

        # cache the fitting scales
        self._sigma = sigma

 def test_same_cache(self):
     poly = Polynom1D()
     poly.pars[1].thaw()
     sdata = DataSimulFit('d1d2d3', (self.d1, self.d2, self.d3))
     smodel = SimulFitModel('same', (poly, poly, poly))
     sfit = Fit(sdata, smodel, method=NelderMead(), stat=Cash())
     result = sfit.fit()
     self.compare_results(self._fit_same_poly_bench, result)

    def fit(self):
        """Fit spectrum"""
        from sherpa.fit import Fit
        from sherpa.models import ArithmeticModel, SimulFitModel
        from sherpa.astro.instrument import Response1D
        from sherpa.data import DataSimulFit

        # Translate model to sherpa model if necessary
        if isinstance(self.model, models.SpectralModel):
            model = self.model.to_sherpa()
        else:
            model = self.model

        if not isinstance(model, ArithmeticModel):
            raise ValueError('Model not understood: {}'.format(model))

        # Make model amplitude O(1e0)
        val = model.ampl.val * self.FLUX_FACTOR ** (-1)
        model.ampl = val

        if self.fit_range is not None:
            log.info('Restricting fit range to {}'.format(self.fit_range))
            fitmin = self.fit_range[0].to('keV').value
            fitmax = self.fit_range[1].to('keV').value

        # Loop over observations
        pha = list()
        folded_model = list()
        nobs = len(self.obs_list)
        for ii in range(nobs):
            temp = self.obs_list[ii].to_sherpa()
            if self.fit_range is not None:
                temp.notice(fitmin, fitmax)
                if temp.get_background() is not None:
                    temp.get_background().notice(fitmin, fitmax)
            temp.ignore_bad()
            if temp.get_background() is not None:
                temp.get_background().ignore_bad()
            pha.append(temp)
            # Forward folding
            resp = Response1D(pha[ii])
            folded_model.append(resp(model) * self.FLUX_FACTOR)

        data = DataSimulFit('simul fit data', pha)
        fitmodel = SimulFitModel('simul fit model', folded_model)

        log.debug(fitmodel)
        fit = Fit(data, fitmodel, self.statistic)
        fitresult = fit.fit()
        log.debug(fitresult)
        # The model instance passed to the Fit now holds the best fit values
        covar = fit.est_errors()
        log.debug(covar)

        for ii in range(nobs):
            efilter = pha[ii].get_filter()
            shmodel = fitmodel.parts[ii]
            self.result[ii].fit = _sherpa_to_fitresult(shmodel, covar, efilter, fitresult)

 def test_gauss_gauss(self):
     g1, g2 = Gauss1D(), Gauss1D()
     g1.fwhm = 1.3
     g1.pos = 1.5
     g2.fwhm = 4.
     g2.pos = -2.0
     sdata = DataSimulFit('d4d5', (self.d4, self.d5))
     smodel = SimulFitModel('g1g2', (g1, g2))
     sfit = Fit(sdata, smodel, method=LevMar(), stat=LeastSq())
     result = sfit.fit()
     self.compare_results(self._fit_g2g2_bench, result)

 def test_diff_cache(self):
     poly1 = Polynom1D()
     poly2 = Polynom1D()
     poly3 = Polynom1D()
     poly1.pars[1].thaw()
     poly2.pars[1].thaw()
     poly3.pars[1].thaw()
     sdata = DataSimulFit('d123', (self.d1, self.d2, self.d3))
     smodel = SimulFitModel('diff', (poly1, poly2, poly3))
     sfit = Fit(sdata, smodel, method=NelderMead(), stat=Cash())
     result = sfit.fit()
     self.compare_results(self._fit_diff_poly_bench, result)

 def test_simul_stat_fit(self):
     data1 = self.data_pi2278
     data2 = self.data_pi2286
     model1 = self.model_mult
     model2 = self.model_mult
     data = DataSimulFit(name='data1data2', datasets=[data1, data2])
     model = SimulFitModel(name='model1model2', parts=[model1, model2])
     fit = Fit(data=data, model=model, stat=MyChiNoBkg(),
               method=NelderMead())
     result = fit.fit()
     self.compare_results(self._fit_simul_datavarstat_results_bench,
                          result)

def test_rsp1d_matrix_pha_zero_energy_bin():
    """What happens when the first bin starts at 0, with replacement.

    Unlike test_rsp1d_delta_pha_zero_energy_bin this directly
    calls Response1D to create the model.
    """

    ethresh = 1.0e-5

    rdata = create_non_delta_rmf()

    # hack the first bin to have 0 energy
    rdata.energ_lo[0] = 0.0

    # PHA and ARF have different exposure ties
    exposure_arf = 0.1
    exposure_pha = 2.4

    specresp = create_non_delta_specresp()

    with warnings.catch_warnings(record=True) as ws:
        warnings.simplefilter("always")
        adata = create_arf(rdata.energ_lo,
                           rdata.energ_hi,
                           specresp,
                           exposure=exposure_arf,
                           ethresh=ethresh)

    validate_zero_replacement(ws, 'ARF', 'user-arf', ethresh)

    nchans = rdata.e_min.size
    channels = np.arange(1, nchans + 1, dtype=np.int16)
    counts = np.ones(nchans, dtype=np.int16)
    pha = DataPHA('test-pha', channel=channels, counts=counts,
                  exposure=exposure_pha)
    pha.set_rmf(rdata)
    pha.set_arf(adata)

    pha.set_analysis('energy')

    mdl = MyPowLaw1D()

    rsp = Response1D(pha)
    wrapped = rsp(mdl)

    # Evaluate the statistic / model. The value was calculated using
    # commit a65fb94004664eab219cc09652172ffe1dad80a6 on a linux
    # system (Ubuntu 17.04).
    #
    f = Fit(pha, wrapped)
    ans = f.calc_stat()
    assert ans == pytest.approx(37971.8716151947)

    def _fit_sherpa(self):
        """Wrapper around sherpa minimizer."""
        from sherpa.fit import Fit
        from sherpa.data import Data1DInt
        from sherpa.optmethods import NelderMead
        from .sherpa_utils import SherpaModel, SherpaStat

        binning = self.obs_list[0].e_reco
        # The sherpa data object is not usued in the fit. It is set to the
        # first observation for debugging purposes, see below
        data = self.obs_list[0].on_vector.data.data.value
        data = Data1DInt('Dummy data', binning[:-1].value,
                         binning[1:].value, data)
        # DEBUG
        # from sherpa.models import PowLaw1D
        # from sherpa.stats import Cash
        # model = PowLaw1D('sherpa')
        # model.ref = 0.1
        # fit = Fit(data, model, Cash(), NelderMead())

        # NOTE: We cannot use the Levenbergr-Marquart optimizer in Sherpa
        # because it relies on the fvec return value of the fit statistic (we
        # return None). The computation of fvec is not straightforwad, not just
        # stats per bin. E.g. for a cash fit the sherpa stat computes it
        # according to cstat
        # see https://github.com/sherpa/sherpa/blob/master/sherpa/include/sherpa/stats.hh#L122

        self._sherpa_fit = Fit(data,
                               SherpaModel(self),
                               SherpaStat(self),
                               NelderMead())
        fitresult = self._sherpa_fit.fit()
        log.debug(fitresult)
        self._make_fit_result()

    def setUp(self):
        try:
            from sherpa.astro.io import read_pha
            from sherpa.astro.xspec import XSwabs, XSpowerlaw
        except:
            return
        #self.startdir = os.getcwd()
        self.old_level = logger.getEffectiveLevel()
        logger.setLevel(logging.CRITICAL)

        datadir = SherpaTestCase.datadir
        if datadir is None:
            return

        pha = os.path.join(datadir, "refake_0934_1_21_1e4.fak")
        rmf = os.path.join(datadir, "ccdid7_default.rmf")
        arf = os.path.join(datadir, "quiet_0934.arf")
        
        self.simarf = os.path.join(datadir, "aref_sample.fits")
        self.pcaarf = os.path.join(datadir, "aref_Cedge.fits")

        data = read_pha(pha)
        data.ignore(None,0.3)
        data.ignore(7.0,None)

        rsp = Response1D(data)
        self.abs1 = XSwabs('abs1')
        self.p1 = XSpowerlaw('p1')
        model = rsp(self.abs1*self.p1)
        
        self.fit = Fit(data, model, CStat(), NelderMead(), Covariance())

def test_sherpa_crab_fit():
    from sherpa.models import NormGauss2D, PowLaw1D, TableModel, Const2D
    from sherpa.stats import Chi2ConstVar
    from sherpa.optmethods import LevMar
    from sherpa.fit import Fit
    from ..sherpa_ import CombinedModel3D

    filename = gammapy_extra.filename('experiments/sherpa_cube_analysis/counts.fits.gz')
    # Note: The cube is stored in incorrect format
    counts = SkyCube.read(filename, format='fermi-counts')
    cube = counts.to_sherpa_data3d()

    # Set up exposure table model
    filename = gammapy_extra.filename('experiments/sherpa_cube_analysis/exposure.fits.gz')
    exposure_data = fits.getdata(filename)
    exposure = TableModel('exposure')
    exposure.load(None, exposure_data.ravel())

    # Freeze exposure amplitude
    exposure.ampl.freeze()

    # Setup combined spatial and spectral model
    spatial_model = NormGauss2D('spatial-model')
    spectral_model = PowLaw1D('spectral-model')
    source_model = CombinedModel3D(spatial_model=spatial_model, spectral_model=spectral_model)

    # Set starting values
    source_model.gamma = 2.2
    source_model.xpos = 83.6
    source_model.ypos = 22.01
    source_model.fwhm = 0.12
    source_model.ampl = 0.05

    model = 1E-9 * exposure * source_model  # 1E-9 flux factor

    # Fit
    fit = Fit(data=cube, model=model, stat=Chi2ConstVar(), method=LevMar())
    result = fit.fit()

    reference = [0.121556,
                 83.625627,
                 22.015564,
                 0.096903,
                 2.240989]

    assert_allclose(result.parvals, reference, rtol=1E-3)

def test_sherpa_crab_fit():
    from sherpa.models import NormGauss2D, PowLaw1D, TableModel, Const2D
    from sherpa.stats import Chi2ConstVar
    from sherpa.optmethods import LevMar
    from sherpa.fit import Fit
    from ..sherpa_ import Data3D, CombinedModel3D

    filename = gammapy_extra.filename('experiments/sherpa_cube_analysis/counts.fits.gz')
    counts = SkyCube.read(filename)
    cube = counts.to_sherpa_data3d()

    # Set up exposure table model
    filename = gammapy_extra.filename('experiments/sherpa_cube_analysis/exposure.fits.gz')
    exposure_data = fits.getdata(filename)
    exposure = TableModel('exposure')
    exposure.load(None, exposure_data.ravel())

    # Freeze exposure amplitude
    exposure.ampl.freeze()

    # Setup combined spatial and spectral model
    spatial_model = NormGauss2D('spatial-model')
    spectral_model = PowLaw1D('spectral-model')
    source_model = CombinedModel3D(spatial_model=spatial_model, spectral_model=spectral_model)

    # Set starting values
    source_model.gamma = 2.2
    source_model.xpos = 83.6
    source_model.ypos = 22.01
    source_model.fwhm = 0.12
    source_model.ampl = 0.05

    model = 1E-9 * exposure * (source_model) # 1E-9 flux factor

    # Fit
    fit = Fit(data=cube, model=model, stat=Chi2ConstVar(), method=LevMar())
    result = fit.fit()

    reference = (0.11925401159500593,
                 83.640630749333056,
                 22.020525848447541,
                 0.036353759774770608,
                 1.1900312815970555)

    assert_allclose(result.parvals, reference, rtol=1E-8)

def mwl_fit_low_level_calling_fermi():
    """Example how to do a Sherpa model fit,
    but use the Fermi ScienceTools to evaluate
    the likelihood for the Fermi dataset.
    """

    spec_model = LogParabola('spec_model')
    spec_model.c1 = 0.5
    spec_model.c2 = 0.2
    spec_model.ampl = 5e-11

    model = spec_model

    data = FermiDataShim()
    stat = FermiStatShim()
    method = LevMar()
    fit = Fit(data=data, model=model, stat=stat, method=method)
    result = fit.fit()

    # IPython.embed()
    return Bunch(results=result, model=spec_model)

    def setUp(self):
        data = Data1D('fake', self._x, self._y, self._err)

        g1 = Gauss1D('g1')
        g1.fwhm.set(1.0, _tiny, _max, frozen=False)
        g1.pos.set(1.0, -_max, _max, frozen=False)
        g1.ampl.set(1.0, -_max, _max, frozen=False)
        p1 = PowLaw1D('p1')
        p1.gamma.set(1.0, -10, 10, frozen=False)
        p1.ampl.set(1.0, 0.0, _max, frozen=False)
        p1.ref.set(1.0, -_max, _max, frozen=True)      
        model = p1 + g1

        method = LevMar()
        method.config['maxfev'] = 10000
        method.config['ftol'] = float(_eps)
        method.config['epsfcn'] = float(_eps)
        method.config['gtol'] = float(_eps)
        method.config['xtol'] = float(_eps)
        method.config['factor'] = float(100)

        self.fit = Fit(data, model, Chi2DataVar(), method, Covariance())
        results = self.fit.fit()

        for key in ["succeeded", "numpoints", "nfev"]:
            assert self._fit_results_bench[key] == int(getattr(results, key))

        for key in ["rstat", "qval", "statval", "dof"]:
            assert numpy.allclose(float(self._fit_results_bench[key]),
                                  float(getattr(results, key)),
                                  1.e-7, 1.e-7)

        for key in ["parvals"]:
            try:
                assert numpy.allclose(self._fit_results_bench[key],
                                      getattr(results, key),
                                      1.e-4, 1.e-4)
            except AssertionError:
                print 'parvals bench: ', self._fit_results_bench[key]
                print 'parvals fit:   ', getattr(results, key)
                print 'results', results
                raise

        covresults = self.fit.est_errors()
        self.dof = results.dof
        self.mu = numpy.array(results.parvals)
        self.cov = numpy.array(covresults.extra_output)
        self.num = 10

    def setUp(self):
        from sherpa.astro.io import read_pha
        from sherpa.astro.xspec import XSwabs, XSpowerlaw
        # self.startdir = os.getcwd()
        self.old_level = logger.getEffectiveLevel()
        logger.setLevel(logging.CRITICAL)

        pha = self.make_path("refake_0934_1_21_1e4.fak")
        # rmf = self.make_path("ccdid7_default.rmf")
        # arf = self.make_path("quiet_0934.arf")

        self.simarf = self.make_path("aref_sample.fits")
        self.pcaarf = self.make_path("aref_Cedge.fits")

        data = read_pha(pha)
        data.ignore(None, 0.3)
        data.ignore(7.0, None)

        rsp = Response1D(data)
        self.abs1 = XSwabs('abs1')
        self.p1 = XSpowerlaw('p1')
        model = rsp(self.abs1 * self.p1)

        self.fit = Fit(data, model, CStat(), NelderMead(), Covariance())

 def test_mycash_nobkgdata_modelhasbkg(self):
     data = self.bkg
     fit = Fit(data, self.model, MyCashWithBkg(), NelderMead())
     results = fit.fit()
     self.compare_results(self._fit_mycashnobkg_results_bench, results)

 def test_mycash_datahasbkg_modelhasnobkg(self):
     fit = Fit(self.data, self.model, MyCashNoBkg(), NelderMead())
     results = fit.fit()
     self.compare_results(self._fit_mycash_results_bench, results)

 def test_mychi_data_and_model_donothave_bkg(self):
     data = self.bkg
     fit = Fit(data, self.model, MyChiNoBkg(), LevMar())
     results = fit.fit()
     self.compare_results(self._fit_mychinobkg_results_bench, results, 1e-5)

 def test_mycash_data_and_model_donothave_bkg(self):
     data = self.bkg
     fit = Fit(data, self.model, MyCashNoBkg(), NelderMead())
     results = fit.fit()
     self.compare_results(self._fit_mycashnobkg_results_bench, results)

 def test_mychi_data_and_model_have_bkg(self):
     fit = Fit(self.data, self.model, MyChiWithBkg(), LevMar())
     results = fit.fit()
     self.compare_results(self._fit_mychi_results_bench, results)