class test_filter_energy_grid(SherpaTestCase):

    _notice = numpy.ones(46, dtype=bool)
    _notice[44:46]=False

    _ignore = numpy.zeros(46, dtype=bool)
    _ignore[14:33]=True

    _emin = numpy.array([
        1.46000006e-03,   2.48199999e-01,   3.06600004e-01,   4.67200011e-01,
        5.69400012e-01,   6.42400026e-01,   7.00800002e-01,   7.44599998e-01,
        7.88399994e-01,   8.17600012e-01,   8.61400008e-01,   8.90600026e-01,
        9.49000001e-01,   9.92799997e-01,   1.03659999e+00,   1.09500003e+00,
        1.13880002e+00,   1.19719994e+00,   1.28480005e+00,   1.40160000e+00,
        1.47459996e+00,   1.60599995e+00,   1.69360006e+00,   1.81040001e+00,
        1.89800000e+00,   1.94180000e+00,   2.02940011e+00,   2.08780003e+00,
        2.19000006e+00,   2.27760005e+00,   2.39439988e+00,   2.58419991e+00,
        2.71560001e+00,   2.86159992e+00,   3.08060002e+00,   3.38720012e+00,
        3.56240010e+00,   3.79600000e+00,   4.02960014e+00,   4.24860001e+00,
        4.71579981e+00,   5.02239990e+00,   5.37279987e+00,   5.89839983e+00,
        6.57000017e+00,   9.86960030e+00], numpy.float)

    _emax = numpy.array([
        0.2482    ,   0.3066    ,   0.46720001,   0.56940001,   0.64240003,
        0.7008    ,   0.7446    ,   0.78839999,   0.81760001,   0.86140001,
        0.89060003,   0.949     ,   0.9928    ,   1.03659999,   1.09500003,
        1.13880002,   1.19719994,   1.28480005,   1.4016    ,   1.47459996,
        1.60599995,   1.69360006,   1.81040001,   1.898     ,   1.9418    ,
        2.02940011,   2.08780003,   2.19000006,   2.27760005,   2.39439988,
        2.58419991,   2.71560001,   2.86159992,   3.08060002,   3.38720012,
        3.5624001 ,   3.796     ,   4.02960014,   4.24860001,   4.71579981,
        5.0223999 ,   5.37279987,   5.89839983,   6.57000017,   9.8696003 ,
        14.95040035], numpy.float)

    def setUp(self):
        self.old_level = logger.getEffectiveLevel()
        logger.setLevel(logging.ERROR)
        self.pha = DataPHA('', numpy.arange(46, dtype=float)+1.,
                           numpy.zeros(46),
                           bin_lo = self._emin, bin_hi = self._emax )
        self.pha.units="energy"

    def tearDown(self):
        logger.setLevel(self.old_level)

    def test_notice(self):
        #clear mask
        self.pha.notice()        
        self.pha.notice(0.0, 6.0)
        #self.assertEqual(self._notice, self.pha.mask)
        assert (self._notice==numpy.asarray(self.pha.mask)).all()


    def test_ignore(self):
        #clear mask
        self.pha.notice()
        self.pha.ignore(0.0, 1.0)
        self.pha.ignore(3.0, 15.0)
        #self.assertEqual(self._ignore, self.pha.mask)
        assert (self._ignore==numpy.asarray(self.pha.mask)).all()

class test_filter_wave_grid(SherpaTestCase):

    _notice = np.ones(16384, dtype=bool)
    _notice[8465:16384] = False

    _ignore = np.zeros(16384, dtype=bool)
    _ignore[14064:15984] = True

    _emin = np.arange(205.7875, 0.9875, -0.0125)

    _emax = _emin + 0.0125

    def setUp(self):
        self.old_level = logger.getEffectiveLevel()
        logger.setLevel(logging.ERROR)
        self.pha = DataPHA('', np.arange(16384, dtype=float) + 1,
                           np.zeros(16384),
                           bin_lo=self._emin,
                           bin_hi=self._emax)

    def tearDown(self):
        logger.setLevel(self.old_level)

    def test_notice(self):
        self.pha.units = 'wavelength'
        self.pha.notice()
        self.pha.notice(100.0, 225.0)
        assert (self._notice == np.asarray(self.pha.mask)).all()

    def test_ignore(self):
        self.pha.units = 'wavelength'
        self.pha.notice()
        self.pha.ignore(30.01, 225.0)
        self.pha.ignore(0.1, 6.0)
        assert (self._ignore == np.asarray(self.pha.mask)).all()

def test_rsp_normf_error(analysis):
    """Check that an error is raised on set_analysis

    """

    exposure = 200.1

    # rdata is only used to define the grids
    rdata = create_non_delta_rmf()
    specresp = create_non_delta_specresp()
    adata = create_arf(rdata.energ_lo,
                       rdata.energ_hi,
                       specresp,
                       exposure=exposure)

    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.set_arf(adata)

    with pytest.raises(DataErr) as exc:
        pha.set_analysis(analysis)

    emsg = "response incomplete for dataset test-pha, " + \
           "check the instrument model"
    assert str(exc.value) == emsg

def calc_wstat_sherpa(mu_sig, n_on, n_off, alpha):
    import sherpa.stats as ss
    from sherpa.astro.data import DataPHA
    from sherpa.models import Const1D
    wstat = ss.WStat()

    model = Const1D()
    model.c0 = mu_sig
    data = DataPHA(counts=np.atleast_1d(n_on),
                   name='dummy',
                   channel=np.atleast_1d(1),
                   backscal=1,
                   exposure=1)
    background = DataPHA(counts=np.atleast_1d(n_off),
                         name='dummy background',
                         channel=np.atleast_1d(1),
                         backscal=np.atleast_1d(1. / alpha),
                         exposure=1)

    data.set_background(background, 1)

    # Docstring for ``calc_stat``
    # https://github.com/sherpa/sherpa/blob/fe8508818662346cb6d9050ba676e23318e747dd/sherpa/stats/__init__.py#L219

    stat = wstat.calc_stat(model=model, data=data)
    print("Sherpa stat: {}".format(stat[0]))
    print("Sherpa fvec: {}".format(stat[1]))

def test_rsp_no_arf_matrix_call(analysis, phaexp):
    """Check out Response1D with matrix but no ARF

    analysis is the analysis setting
    arfexp determines whether the arf has an exposure time
    phaexp determines whether the PHA has an exposure time
    """

    if phaexp:
        pha_exposure = 220.9
    else:
        pha_exposure = None

    if phaexp:
        exposure = pha_exposure
        mdl_label = '({} * flat)'.format(exposure)
    else:
        exposure = 1.0
        mdl_label = 'flat'

    rdata = create_non_delta_rmf()

    constant = 2.3
    mdl = Const1D('flat')
    mdl.c0 = constant

    # Turn off integration on this model, so that it is not integrated
    # across the bin width.
    #
    mdl.integrate = False

    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=pha_exposure)

    pha.set_rmf(rdata)

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

    assert isinstance(wrapped, ArithmeticModel)

    expname = 'apply_rmf({})'.format(mdl_label)
    assert wrapped.name == expname

    modvals = exposure * constant * np.ones(rdata.energ_lo.size)
    matrix = get_non_delta_matrix()
    expected = np.matmul(modvals, matrix)

    pha.set_analysis(analysis)
    out = wrapped([4, 5])
    assert_allclose(out, expected)

 def setUp(self):
     self.old_level = logger.getEffectiveLevel()
     logger.setLevel(logging.ERROR)
     self.pha = DataPHA('', np.arange(16384, dtype=float) + 1,
                        np.zeros(16384),
                        bin_lo=self._emin,
                        bin_hi=self._emax)

 def setUp(self):
     self.old_level = logger.getEffectiveLevel()
     logger.setLevel(logging.ERROR)
     self.pha = DataPHA('', numpy.arange(46, dtype=float)+1.,
                        numpy.zeros(46),
                        bin_lo = self._emin, bin_hi = self._emax )
     self.pha.units="energy"

def test_arfmodelpha_call(ignore):
    """What happens calling an arf with a pha?

    The ignore value indicates what channel to ignore (0 means
    nothing is ignored). The aim is to check edge effects,
    and as there are only a few channels, it was decided to
    test all channels.
    """

    # Note: the exposure is set in the PHA and ARF, but should not be
    #       used when evaluating the model; it's value has been
    #       set to a value that the test will fail it it is.
    #
    exposure = 200.1
    estep = 0.01
    egrid = np.arange(0.01, 0.06, estep)
    svals = [1.1, 1.2, 1.3, 1.4]
    specresp = np.asarray(svals)
    adata = create_arf(egrid[:-1], egrid[1:], specresp,
                       exposure=exposure)

    constant = 2.3
    mdl = Const1D('flat')
    mdl.c0 = constant

    channels = np.arange(1, 5, dtype=np.int16)
    counts = np.asarray([10, 5, 12, 7], dtype=np.int16)
    pha = DataPHA('test-pha', channel=channels, counts=counts,
                  exposure=exposure)
    pha.set_arf(adata)

    # force energy units (only needed if ignore is set)
    pha.set_analysis('energy')

    if ignore is not None:
        de = estep * 0.9
        e0 = egrid[ignore]
        pha.notice(lo=e0, hi=e0 + de, ignore=True)

        # The assert are intended to help people reading this
        # code rather than being a useful check that the code
        # is working.
        mask = [True, True, True, True]
        mask[ignore] = False
        assert (pha.mask == mask).all()

    wrapped = ARFModelPHA(adata, pha, mdl)

    # The model is evaluated on the ARF grid, not whatever
    # is sent in. It is also integrated across the bins,
    # which is why there is a multiplication by the
    # grid width (for this constant model).
    #
    # Note that the filter doesn't change the grid.
    #
    de = egrid[1:] - egrid[:-1]
    expected = constant * np.asarray(svals) * de
    out = wrapped([4, 5])
    assert_allclose(out, expected)

def test_rspmodelpha_delta_call_wave():
    """What happens calling a rsp with a pha (RMF is a delta fn)? Wavelength.

    Unlike the energy case no bins are ignored, as this code path
    has already been tested.
    """

    exposure = 200.1
    estep = 0.025
    egrid = np.arange(0.1, 0.8, estep)
    elo = egrid[:-1]
    ehi = egrid[1:]
    specresp = 2.4 * np.ones(elo.size, dtype=np.float32)
    specresp[2:5] = 0.0
    specresp[16:19] = 3.2
    adata = create_arf(elo, ehi, specresp, exposure=exposure)
    rdata = create_delta_rmf(elo, ehi, e_min=elo, e_max=ehi)
    nchans = elo.size

    constant = 2.3
    mdl = Const1D('flat')
    mdl.c0 = constant

    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.set_rmf(rdata)

    pha.set_analysis('wave')

    wrapped = RSPModelPHA(adata, rdata, pha, mdl)

    # Note that this is a Sherpa model, so it's normalization is
    # per unit x axis, so when integrated here the bins are in
    # Angstroms, so the bin width to multiply by is
    # Angstroms, not keV.
    #
    dl = (DataPHA._hc / elo) - (DataPHA._hc / ehi)
    expected = constant * specresp * dl

    out = wrapped([4, 5])
    assert_allclose(out, expected)

def test_rspmodelpha_matrix_call_xspec():
    """Check XSPEC constant is invariant to wavelength/energy setting.

    As XSPEC models internally convert from Angstrom to keV,
    do a simple check here.
    """

    exposure = 200.1
    rdata = create_non_delta_rmf()
    specresp = create_non_delta_specresp()
    adata = create_arf(rdata.energ_lo,
                       rdata.energ_hi,
                       specresp,
                       exposure=exposure)

    constant = 2.3
    mdl = XSconstant('flat')
    mdl.factor = constant

    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)

    # The set_arf call isn't necessary, but leave in
    pha.set_arf(adata)
    pha.set_rmf(rdata)

    # The XSPEC models are evaluated on an energy grid, even when
    # the analysis setting is wavelength. Also, unlike the Sherpa
    # Constant model, the XSPEC XSconstant model is defined
    # over the integrated bin, so no correction is needed for the
    # bin width.
    #
    modvals = constant * specresp
    matrix = get_non_delta_matrix()
    expected = np.matmul(modvals, matrix)

    wrapped = RSPModelPHA(adata, rdata, pha, mdl)

    pha.set_analysis('wave')
    out_wl = wrapped([4, 5])
    assert_allclose(out_wl, expected)

    pha.set_analysis('energy')
    out_en = wrapped([4, 5])
    assert_allclose(out_en, expected)

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

    ethresh = 1.0e-10

    # Note: the two exposures are different to check which is
    #       used (the answer is neither, which seems surprising)
    #
    exposure1 = 0.1
    egrid = np.asarray([0.0, 0.1, 0.2, 0.4, 0.5, 0.7, 0.8])
    elo = egrid[:-1]
    ehi = egrid[1:]
    specresp = np.asarray([10.2, 9.8, 10.0, 12.0, 8.0, 10.0])

    with warnings.catch_warnings(record=True) as ws:
        warnings.simplefilter("always")
        adata = create_arf(elo, ehi, specresp, exposure=exposure1,
                           ethresh=ethresh)

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

    arf = ARF1D(adata)

    exposure2 = 2.4
    channels = np.arange(1, 7, dtype=np.int16)
    counts = np.ones(6, dtype=np.int16)
    pha = DataPHA('test-pha', channel=channels, counts=counts,
                  exposure=exposure2)
    pha.set_arf(adata)

    pha.set_analysis('energy')

    mdl = MyPowLaw1D()
    tmdl = PowLaw1D()

    wrapped = ARFModelPHA(arf, pha, mdl)

    out = wrapped([0.1, 0.2])
    elo[0] = ethresh
    expected = specresp * tmdl(elo, ehi)

    assert_allclose(out, expected)
    assert not np.isnan(out[0])

def test_rspmodelpha_delta_call(ignore):
    """What happens calling a rsp with a pha (RMF is a delta fn)?

    The ignore value gives the channel to ignore (counting from 0).
    """

    exposure = 200.1
    estep = 0.025
    egrid = np.arange(0.1, 0.8, estep)
    elo = egrid[:-1]
    ehi = egrid[1:]
    specresp = 2.4 * np.ones(elo.size, dtype=np.float32)
    specresp[2:5] = 0.0
    specresp[16:19] = 3.2
    adata = create_arf(elo, ehi, specresp, exposure=exposure)
    rdata = create_delta_rmf(elo, ehi, e_min=elo, e_max=ehi)
    nchans = elo.size

    constant = 2.3
    mdl = Const1D('flat')
    mdl.c0 = constant

    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.set_rmf(rdata)

    # force energy units (only needed if ignore is set)
    pha.set_analysis('energy')

    if ignore is not None:
        de = estep * 0.9
        e0 = egrid[ignore]
        pha.notice(lo=e0, hi=e0 + de, ignore=True)

        # The assert are intended to help people reading this
        # code rather than being a useful check that the code
        # is working.
        mask = [True] * nchans
        mask[ignore] = False
        assert (pha.mask == mask).all()

    wrapped = RSPModelPHA(adata, rdata, pha, mdl)

    # The model is evaluated on the RMF grid, not whatever
    # is sent in. It is also integrated across the bins,
    # which is why there is a multiplication by the
    # grid width (for this constant model).
    #
    # Note that the filter doesn't change the grid.
    #
    de = egrid[1:] - egrid[:-1]
    expected = constant * specresp * de
    out = wrapped([4, 5])
    assert_allclose(out, expected)

def test_rspmodelpha_delta_call_channel():
    """What happens calling a rsp with a pha (RMF is a delta fn)? Channels.

    I am not convinced I understand the bin width calculation here,
    as it doesn't seem to match the wavelength case.
    """

    exposure = 200.1
    estep = 0.025
    egrid = np.arange(0.1, 0.8, estep)
    elo = egrid[:-1]
    ehi = egrid[1:]
    specresp = 2.4 * np.ones(elo.size, dtype=np.float32)
    specresp[2:5] = 0.0
    specresp[16:19] = 3.2
    adata = create_arf(elo, ehi, specresp, exposure=exposure)
    rdata = create_delta_rmf(elo, ehi, e_min=elo, e_max=ehi)
    nchans = elo.size

    constant = 2.3
    mdl = Const1D('flat')
    mdl.c0 = constant

    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.set_rmf(rdata)

    pha.set_analysis('channel')

    wrapped = RSPModelPHA(adata, rdata, pha, mdl)

    # Since this is channels you might expect the bin width to be 1,
    # but it is actually still dE.
    #
    de = ehi - elo
    expected = constant * specresp * de

    out = wrapped([4, 5])
    assert_allclose(out, expected)

def test_rmfmodelpha_matrix_mismatch(analysis):
    """Check that an error is raised if there's a mismatch.

    """

    exposure = 200.1
    rdata = create_non_delta_rmf()

    # nchans should be rdata.e_min.size for the sizes to match
    nchans = rdata.energ_lo.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.set_rmf(rdata)

    with pytest.raises(DataErr) as exc:
        pha.set_analysis(analysis)

    emsg = "RMF 'non-delta-rmf' is incompatible with PHA dataset 'test-pha'"
    assert str(exc.value) == emsg

def test_rspmodelpha_matrix_call(ignore):
    """What happens calling a rsp with a pha (RMF is a matrix)?

    The ignore value gives the channel to ignore (counting from 0).
    """

    exposure = 200.1
    rdata = create_non_delta_rmf()
    specresp = create_non_delta_specresp()
    elo = rdata.energ_lo
    ehi = rdata.energ_hi

    adata = create_arf(elo, ehi, specresp, exposure=exposure)
    nchans = rdata.e_min.size

    constant = 22.3
    slope = -1.2
    mdl = Polynom1D('sloped')
    mdl.c0 = constant
    mdl.c1 = slope

    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.set_rmf(rdata)

    # force energy units (only needed if ignore is set)
    pha.set_analysis('energy')

    if ignore is not None:
        e0 = rdata.e_min[ignore]
        e1 = rdata.e_max[ignore]
        de = 0.9 * (e1 - e0)
        pha.notice(lo=e0, hi=e0 + de, ignore=True)

        # The assert are intended to help people reading this
        # code rather than being a useful check that the code
        # is working.
        mask = [True] * nchans
        mask[ignore] = False
        assert (pha.mask == mask).all()

    wrapped = RSPModelPHA(adata, rdata, pha, mdl)

    # The filter does not change the grid
    modvals = specresp * mdl(rdata.energ_lo, rdata.energ_hi)
    matrix = get_non_delta_matrix()
    expected = np.matmul(modvals, matrix)

    out = wrapped([4, 5])
    assert_allclose(out, expected)

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 test_rmfmodelpha_matrix_call(ignore):
    """What happens calling an rmf (matrix) with a pha?

    The ignore value gives the channel to ignore (counting from 0).
    """

    exposure = 200.1
    rdata = create_non_delta_rmf()
    elo = rdata.e_min
    ehi = rdata.e_max
    nchans = elo.size

    constant = 12.2
    slope = 0.01
    mdl = Polynom1D('not-flat')
    mdl.c0 = constant
    mdl.c1 = slope

    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.set_rmf(rdata)

    # force energy units (only needed if ignore is set)
    pha.set_analysis('energy')

    if ignore is not None:
        e0 = elo[ignore]
        e1 = ehi[ignore]
        de = 0.9 * (e1 - e0)
        pha.notice(lo=e0, hi=e0 + de, ignore=True)

        # The assert are intended to help people reading this
        # code rather than being a useful check that the code
        # is working.
        mask = [True] * nchans
        mask[ignore] = False
        assert (pha.mask == mask).all()

    wrapped = RMFModelPHA(rdata, pha, mdl)

    # Note that the evaluation ignores any filter we've applied.
    # and the exposure time is not used.
    #
    modvals = mdl(rdata.energ_lo, rdata.energ_hi)
    matrix = get_non_delta_matrix()
    expected = np.matmul(modvals, matrix)

    out = wrapped([4, 5])
    assert_allclose(out, expected)

def get_data(mu_sig, n_on, n_off, alpha):
    from sherpa.astro.data import DataPHA
    from sherpa.models import Const1D

    model = Const1D()
    model.c0 = mu_sig
    data = DataPHA(
        counts=np.atleast_1d(n_on),
        name="dummy",
        channel=np.atleast_1d(1),
        backscal=1,
        exposure=1,
    )
    background = DataPHA(
        counts=np.atleast_1d(n_off),
        name="dummy background",
        channel=np.atleast_1d(1),
        backscal=np.atleast_1d(1. / alpha),
        exposure=1,
    )

    data.set_background(background, 1)

    return model, data

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

    ethresh = 2e-7

    egrid = np.asarray([0.0, 0.1, 0.2, 0.4, 0.5, 0.7, 0.8])
    elo = egrid[:-1]
    ehi = egrid[1:]

    with warnings.catch_warnings(record=True) as ws:
        warnings.simplefilter("always")
        rdata = create_delta_rmf(elo, ehi, ethresh=ethresh)

    validate_zero_replacement(ws, 'RMF', 'delta-rmf', ethresh)

    exposure = 2.4
    channels = np.arange(1, 7, dtype=np.int16)
    counts = np.ones(6, dtype=np.int16)
    pha = DataPHA('test-pha', channel=channels, counts=counts,
                  exposure=exposure)
    pha.set_rmf(rdata)

    pha.set_analysis('energy')

    mdl = MyPowLaw1D()
    tmdl = PowLaw1D()

    wrapped = RMFModelPHA(rdata, pha, mdl)

    out = wrapped([0.1, 0.2])

    elo[0] = ethresh
    expected = tmdl(elo, ehi)

    assert_allclose(out, expected)
    assert not np.isnan(out[0])

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

    ethresh = 2.0e-7

    # PHA and ARF have different exposure ties
    exposure1 = 0.1
    exposure2 = 2.4
    egrid = np.asarray([0.0, 0.1, 0.2, 0.4, 0.5, 0.7, 0.8])
    elo = egrid[:-1]
    ehi = egrid[1:]
    specresp = np.asarray([10.2, 9.8, 10.0, 12.0, 8.0, 10.0])

    with warnings.catch_warnings(record=True) as ws:
        warnings.simplefilter("always")
        adata = create_arf(elo, ehi, specresp, exposure=exposure1,
                           ethresh=ethresh)

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

    with warnings.catch_warnings(record=True) as ws:
        warnings.simplefilter("always")
        rdata = create_delta_rmf(elo, ehi, ethresh=ethresh)

    validate_zero_replacement(ws, 'RMF', 'delta-rmf', ethresh)

    channels = np.arange(1, 7, dtype=np.int16)
    counts = np.ones(6, dtype=np.int16)
    pha = DataPHA('test-pha', channel=channels, counts=counts,
                  exposure=exposure2)
    pha.set_rmf(rdata)
    pha.set_arf(adata)

    pha.set_analysis('energy')

    mdl = MyPowLaw1D()
    tmdl = PowLaw1D()

    wrapped = RSPModelPHA(adata, rdata, pha, mdl)

    out = wrapped([0.1, 0.2])

    elo[0] = ethresh
    expected = specresp * tmdl(elo, ehi)

    assert_allclose(out, expected)
    assert not np.isnan(out[0])