Skip to content

Commit

Permalink
Fix RMF handling for non-variable-length RMF files (AstroPy)
Browse files Browse the repository at this point in the history
PR sherpa#358 added AstroPy support for RMF files where the data was not
encoded using FITS variable-length arrays for the matrix. This means
that the data can not just be concatenated together. Instead each
row in the matrix (which, in this case, is a regular nchannels by
nbins array) has to be subset (maybe multiple times).

The commit I added to do this - cfad850
worked in the Swift case, where the F_CHAN/N_CHAN values were set up
to effectively remove the "space saving" mechanism they provide (that
is, all elements of each row in the MATRIX column are used). It does
not work for cases where F_CHAN/N_CHAN are not set to 0/nbins for
each row. The code was embarassingly-confused over what it was meant
to be doing, and only worked for the Swift case "by luck".

This commit should fix this.
  • Loading branch information
DougBurke committed Dec 7, 2017
1 parent de3338b commit d83681b
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 59 deletions.
41 changes: 17 additions & 24 deletions sherpa/astro/io/pyfits_backend.py
Original file line number Diff line number Diff line change
Expand Up @@ -741,23 +741,11 @@ def get_rmf_data(arg, make_copy=False):
finally:
rmf.close()

# ## For every row i of the response matrix, such that
# ## n_grp[i] == 0, we need to remove that row from the
# ## n_chan, f_chan, and matrix arrays we are constructing
# ## to be passed up to the DataRMF data structure.

# ## This is trivial for n_chan and f_chan. For the matrix
# ## array this can be more work -- can't just remove all
# ## zeroes, because some rows where n_grp[row] > 0 might
# ## still have zeroes in the matrix. I add new code first
# ## to deal with the matrix, then simpler code to remove zeroes
# ## from n_chan and f_chan.

# Read in MATRIX column with structure as-is -- i.e., as an array of
# arrays. Then flatten it, but include only those arrays that come from
# rows where n_grp[row] > 0. Zero elements can only be included from
# rows where n_grp[row] > 0. SMD 05/23/13

# Remove any rows from the MATRIX and F_CHAN/N_CHAN where N_GRP is 0,
# since they do not add any data. This is to match the Crates backend.
# Note that crates uses the sherpa.astro.utils.resp_init routine,
# but it's not clear why, so it is not used here for now.
#
good = (data['n_grp'] > 0)
data['matrix'] = data['matrix'][good]

Expand All @@ -781,14 +769,18 @@ def get_rmf_data(arg, make_copy=False):
n_chan = data['n_chan'][good]

rowdata = []
for mrow, ng, fcs, ncs in zip(matrix, n_grp, f_chan, n_chan):
for mrow, ng, ncs in zip(matrix, n_grp, n_chan):
# Need a RMF which ng>1 to test this with.
if ng == 1:
fcs = [fcs]
ncs = [ncs]
for fc, nc in zip(fcs, ncs):
# It appears that F_CHAN is 0-based
#
rowdata.append(mrow[fc: fc + nc])
start = 0
for nc in ncs:
# n_chan can be an unsigned integer. Adding a Python
# integer to a NumPy unsigned integer appears to return
# a float.
end = start + numpy.int(nc)
rowdata.append(mrow[start:end])
start = end

data['matrix'] = numpy.concatenate(rowdata)

Expand All @@ -797,7 +789,6 @@ def get_rmf_data(arg, make_copy=False):
# Flatten f_chan and n_chan vectors into 1D arrays as crates does
# according to group
#
# TODO: should the n_grp > 0 filter be applied before all this?
if data['f_chan'].ndim > 1 and data['n_chan'].ndim > 1:
f_chan = []
n_chan = []
Expand All @@ -807,6 +798,8 @@ def get_rmf_data(arg, make_copy=False):
f_chan.append(fch[i])
n_chan.append(nch[i])

# This automatically filters out rows where N_GRP is 0.
#
data['f_chan'] = numpy.asarray(f_chan, SherpaUInt)
data['n_chan'] = numpy.asarray(n_chan, SherpaUInt)
else:
Expand Down
38 changes: 3 additions & 35 deletions sherpa/astro/tests/test_astro_data_rosat_unit.py
Original file line number Diff line number Diff line change
Expand Up @@ -429,7 +429,6 @@ def test_read_rmf_crates_py3(make_data_path):
@requires_data
@requires_fits
@requires_backend('pyfits')
@pytest.mark.xfail
def test_read_rmf_pyfits(make_data_path):
"""Can we read in a PSPC RMF using the astropy back end.
"""
Expand Down Expand Up @@ -460,15 +459,16 @@ def test_read_rmf_fails_pha(make_data_path):

@requires_data
@requires_fits
@requires_backend('crates')
@pytest.mark.skip(not six.PY2, reason="Python 3 currently not supported")
def test_can_use_pspc_data(make_data_path):
"""A basic check that we can read in and use the ROSAT PSPC data.
Unlike the previous tests, that directly access the io module,
this uses the ui interface.
"""

if not six.PY2 and (backend == "crates"):
pytest.skip('Python3 and Crates: known to fail')

# The PSPC PHA file does not have the ANCRFILE/RESPFILE keywords
# set up, so the responses has to be manually added.
#
Expand Down Expand Up @@ -506,35 +506,3 @@ def test_can_use_pspc_data(make_data_path):
assert s.dof == 5

assert_allclose(s.statval, sexpected, rtol=0, atol=0.005)


@requires_data
@requires_fits
@requires_backend('pyfits')
def test_can_use_pspc_data_astropy(make_data_path):
"""A basic check that we can read in and use the ROSAT PSPC data.
Unlike the previous tests, that directly access the io module,
this uses the ui interface.
"""

# The PSPC PHA file does not have the ANCRFILE/RESPFILE keywords
# set up, so the responses has to be manually added.
#
ui.load_pha(make_data_path(PHAFILE))
assert ui.get_analysis() == 'channel'

ui.load_rmf(make_data_path(RMFFILE))
assert ui.get_analysis() == 'energy'

ui.set_source(ui.powlaw1d.pl)
ui.set_par('pl.ampl', 0.0003)

# Due to the filtering out of the 0 N_GRP bins in only some, not
# all, columns, we have an inconsisent RMF, so any attempt to
# use it will fail.
#
with pytest.raises(ValueError) as exc:
ui.calc_stat()

assert str(exc.value) == 'RMF data is invalid or inconsistent'

0 comments on commit d83681b

Please sign in to comment.