Friday, May 21, 2010

A very quick note on 2D user models

Our documentation currently doesn't help people out who want to create a 2D user model in Sherpa. Here's a quick run example of using load_user_model() to create a flat model in Python (so a replacement for Sherpa's const2d model).

% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.2 Monday, November 30, 2009

sherpa-1> def flatmod(params, x, y):
print("Lower left = {0},{1}".format(x[0],y[0]))
print("Upper right = {0},{1}".format(x[-1],y[-1]))
return x*0 + params[0]
sherpa-2> load_user_model(flatmod, "reallyflat")
sherpa-3> add_user_pars("reallyflat", ["ampl"])
sherpa-4> print(reallyflat)
usermodel.reallyflat
Param Type Value Min Max Units
----- ---- ----- --- --- -----
reallyflat.ampl thawed 0 -3.40282e+38 3.40282e+38
sherpa-5> reallyflat.ampl = 2

So here I have created the model reallyflat which sets every pixel to the value of 2 (i.e. it ignores the coordinates of the input). As a side effect it prints out the first and last coordinate sent to it (which will be the bottom-left and top-right pixels if there are no image filters). I created the model on the fly from within Sherpa for demonstration purposes; in normal use I would write the model definition in an external file and load it in using execfile() or import.

Now we load in an image and use the model:

sherpa-6> load_data("img.fits")
sherpa-7> set_source(reallyflat)
sherpa-8> image_data()
sherpa-9> image_model()
Lower left = 1.0,1.0
Upper right = 421.0,367.0

Displaying the model causes it to be evaluated, as shown above. If we use image_source() we see the same result (we have no PSF set up so there is no actual difference between the source and model views in this case).

sherpa-10> image_source()
Lower left = 1.0,1.0
Upper right = 421.0,367.0

The default coordinate system is "logical", which matches the FITS convention, where (1,1) is the center of the bottom-left pixel and each pixel is square with a width of 1. Below we change to the "physical" system, which uses the ACIS SKY coordinates (as this is an ACIS image). Note that the coordinate values sent to the model have changed.

sherpa-11> set_coord("physical")
sherpa-12> image_source()
Lower left = 4064.0,3511.0
Upper right = 4484.0,3877.0
sherpa-13> image_model()
Lower left = 4064.0,3511.0
Upper right = 4484.0,3877.0
sherpa-14> image_fit()
Lower left = 4064.0,3511.0
Upper right = 4484.0,3877.0
Lower left = 4064.0,3511.0
Upper right = 4484.0,3877.0

The image_fit() output above shows that the model is evaluated twice; once for the model display and once to create the residual image.

An example of a slightly-more complex 2D model would be:

def mygauss2d(pars,x,y):
"My version of a gaussian"
(sigma, ampl, xpos, ypos) = pars
r2 = (x - xpos)**2 + (y - ypos)**2
return ampl * np.exp(-r2/(sigma*sigma))

load_user_model(mygauss2d, "g2d")
add_user_pars("g2d",
["sigma", "ampl", "xpos", "ypos"],
[5,1,1,1], [0.1,0,0,0])

where I use the optional third and fourth arguments of add_user_pars() to set up both the starting and minimum values of the parameters.

2 comments:

  1. Hi Doug,

    thanks for the nice examples.

    It seems it is possible to implement 3d user models just as easily:
    def mygauss3d(pars, x, y, z):
    "My version of a gaussian"
    (sigma, ampl, xpos, ypos, zpos) = pars
    r2 = (x - xpos)**2 + (y - ypos)**2 + + (z - zpos)**2
    return ampl * np.exp(-r2/(sigma*sigma))

    load_user_model(mygauss3d, "g3d")
    add_user_pars("g3d",
    ["sigma", "ampl", "xpos", "ypos", "zpos"],
    [10, 10, 50, 50, 50])
    print g3d
    print g3d([50],[50],[50])

    Is there a way to make a 3d data set?
    This doesn't work:
    dataspace2d([10, 10, 2])
    show_data()

    Thanks!
    Christoph

    ReplyDelete
  2. Christoph,

    Internally Sherpa is designed so as to support > 2D datasets, it is just that we haven't provided top-level support for such beasts. So, your user model looks correct, but unfortunately we do not have a "dataspace3d" command, and I am not sure how easy it would be to provide support for 3D datasets. Hopefully Brian can clarify this.

    I should also add a post to comment on "integrated" vs "non-integrated" models, since this can change the number of arguments the user model gets. Maybe next week.

    Doug

    ReplyDelete