mkrprm: make radial profile response matrix

[kglotfelty]

New script to create the equivalent of a response matrix for radial profile; that is the fraction of the PSF that is included in adjacent radial bins. Based on notebook/discussion with @hamogu a while ago (like several years).

The output is a 2D array, ie image, ie response matrix.

Hopefully can replace/extent-to using mkpsfmap+dmimgadpt w/ marx when offset issue is fixed.

There is also a new sherpa convolution model: MatrixModel which will multiply (ie conovolve) with a matrix.

import sherpa.astro.ui as ui
import numpy as np

# Create an identity matrix
from sherpa_contrib.matrix_model import MatrixModel

from pycrates import read_file
# ui.load_data(1, "radial.prof[cols rmid,sur_bri,sur_bri_err]", Data1D)
ui.load_data(1,"radial_profile.fits",3,["RMID","SUR_FLUX","SUR_FLUX_ERR"])

response_matrix = read_file("frac.out").get_image().values*1.0
x_vals = ui.get_data().x
mymatrix = MatrixModel(response_matrix, x_vals, name="rprm")

pp = ui.polynom1d("poly")
ui.thaw(pp.c1)
pp.c1 = -1.0
pp.c0 = 12.5

ui.set_source(mymatrix(pp) )
ui.fit()
ui.plot_fit_delchi()

The MatrixModel is definitely rough and there are some oddities that will need to be worked out.

[kglotfelty]

Problem with MatrixModel sherpa model when using notice/ignore

The MatrixModel needs to evaluate the model over the entire dataspace and the data grid must match the response matrix grid. When you notice or ignore a subset of bins, it breaks.

import numpy as np
import sherpa.astro.ui as ui
from matrix_model import MatrixModel

npts = 10
xx = np.arange(npts)+10
yy = np.random.rand(npts)
ee = np.ones_like(xx)*0.01

ui.load_arrays(1, xx, yy, ee, ui.Data1D)

m_model = MatrixModel(np.identity(npts))

ui.set_model(m_model(ui.const1d.c0))
ui.fit()

ui.notice(12,17)
ui.fit()

The first fit() works, the 2nd fit() fails.

Not sure how to get the full dataspace grid?

I don't think a show stopper for release but would be nice to address.

[DougBurke]

So, the Sherpa instrument models that represent an OGIP-like response (RSP, RMF, or ARF [*]) deal with this since they are evaluated on a particular grid (in this case channels) but then evaluate the "wrapped" models - so the const1d in m_model(ui.const1d.c0) - with an energy grid, and then they deal with the conversion/filtering/... [actually, they generally don't bother with the filtering, but they could do].

So you'd want some way to associate the "unfiltered" grid when creating the MatrixModel - or, a la the PSF code, with a method like fold - and then the calc method can deal with all this unpleasantness. Of course, the way the instrument models are set up we have several classes all intertwined, but maybe we can do it with just one class.

[kglotfelty]

I updated the MatrixModel to take in the full X-coordinate grid as part of the constructor with various checks that the full x-coordinate grid is a superset of the x-values that are being evaluated. The model is evaluated over the entire grid and then only the noticed/ignored bins are returned.