Using a private DS9 Session in Sherpa

Currently, Sherpa uses the default DS9 session 'ds9' and this can potentially conflict with other open DS9 windows that are unrelated to Sherpa analysis. To enable private DS9 sessions in Sherpa, consider the following work-around.

import sherpa.image
sherpa.image.backend.DS9._DefTemplate="sherpa"

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.

Simple user-defined models using a function

Typically, users have a particular model in mind for use in Sherpa. Here is a simple user-model example of creating a gaussian plus a constant.

I define a function 'mymodel' that accepts a Python list of parameter values and a Numpy ndarray of x values. It calculates and returns the sum of a gaussian and constant as a NumPy ndarray.

I define a 1D dataspace from [0.1, 10] with a bin size of 0.01 and load the user-model function identified as 'mdl'. Using the identifier, I populate the parameter names and initial values.

Setup the model, fake the data, and plot!

#!/usr/bin/env python
from sherpa.astro.ui import *
import numpy as np


def mymodel(pars, x, *args, **kwargs):
    """
    Evaluate a 1D Gaussian plus constant
    pars = { fwhm, pos, ampl, const } 
    """
    gfactor = 4.0*np.log(2.0)
    ratio = (x - pars[1])/pars[0]
    return pars[2] * np.exp(-gfactor * ratio * ratio) + pars[3]


# create simple grid [0.1, 10] with a bin size of 0.01
dataspace1d(0.1,10,0.01,dstype=Data1D)

# load the user-model mymodel() with parameter names into
# variable 'mdl'
load_user_model(mymodel,'mdl')
add_user_pars('mdl',['fwhm', 'pos', 'ampl', 'const'])

# initialize parameter values
mdl.pos = 5.0
mdl.fwhm = 2.5
mdl.ampl = 75
mdl.const = 0.1

# set the identifier as the associated model
set_model(mdl)

# evaluate, add poisson noise, and populate the
# dataspace
fake()

# plot faked data
plot_data()

So, can I shut Sherpa up?

Say you want to run a bunch of Sherpa commands (maybe in a script, or in an interactive session) but don't want to see the output from the commands. What should you do?

Well you can take advantage of the not-yet-documented-in-CIAO-but-really-should-be Python logging infrastructure used by Sherpa. Start with

import logging

logger = logging.getLogger("sherpa")

and then you turn off the output using one of

logger.setLevel(logging.WARN)

logger.setLevel(logging.ERROR)

and turn it back on again with

logger.setLevel(logging.INFO)

The argument determines what severity of messages are displayed (so logging.WARN means warning and above, which excludes the INFO level used for most output from Sherpa).

As an example, compare the screen output for the three load_data() calls below:

% sherpa

-----------------------------------------------------

Welcome to Sherpa: CXC's Modeling and Fitting Package

-----------------------------------------------------

CIAO 4.2 Monday, November 30, 2009

sherpa-1> import logging

sherpa-2> logger = logging.getLogger("sherpa")

sherpa-3> load_pha("3c273.pi")

WARNING: systematic errors were not found in file '3c273.pi'

statistical errors were found in file '3c273.pi' but not used; to use them, re-read with use_errors=True

read ARF file 3c273.arf

read RMF file 3c273.rmf

WARNING: systematic errors were not found in file '3c273_bg.pi'

statistical errors were found in file '3c273_bg.pi' but not used; to use them, re-read with use_errors=True

read background file 3c273_bg.pi

sherpa-4> logger.setLevel(logging.WARN)

sherpa-5> load_pha("3c273.pi")

WARNING: systematic errors were not found in file '3c273.pi'

WARNING: systematic errors were not found in file '3c273_bg.pi'

sherpa-6> logger.setLevel(logging.ERROR)

sherpa-7> load_pha("3c273.pi")

sherpa-8> logger.setLevel(logging.INFO)

This can be particularly useful for hiding the on-screen results of a call to fit() or one of the error routines such as conf() and proj(); if you do this though you will probably want to use the corresponding get_xxx_results routine to get at the calculated data!

You can do a lot more with this "logger" object, such as re-directing the output to a file (or even to both the screen and a file). See 'help logging' or this logging tutorial (one of many that are out there for you intrepid readers) for more information.

Fake data FAST!

Sometimes, users just what to simulate a simple model and plot. This can be done quickly with five easy steps. Define the grid, the model, the initial model parameter values, add noise, and plot.

Example:

#!/usr/bin/env python
from sherpa.astro.ui import *

# create simple grid [0.1, 10] with a
# bin size of 0.01
dataspace1d(0.1,10,0.01, dstype=Data1D)

# define a gaussian model plus a constant
set_model(gauss1d.g1+const1d.c1)

# initialize parameter values
g1.pos = 5
g1.fwhm = 2.5
g1.ampl = 75
c1.c0 = 0.1

# evaluate the model, add poisson noise,
# and populate the dataspace
fake()

# plot faked data
plot_data()

New Sherpa this Summer

It looks like we will be releasing a new Sherpa version this summer. Standby for updates.

Installing Standalone Sherpa for Python.

Tom has put Instructions on Astropython  to install standalone Sherpa for X-ray analysis. Note that this is mainly useful for Python users who want to import Sherpa to their Python environment.

Combine Sherpa Models with Deprojected Models in Sherpa

Yesterday, I used a nice trick in Sherpa to combine standard models and deprojected models. Using  Deproject Python package in Sherpa one can fit the emission models to X-ray clusters accounting for the 3D to 2D projection effects in the data (standard "onion-peeling" technique described by Fabian et al 1981, and Kriss et al 1983). This allows for obtaining deprojected temperature, density and entropy distributions in the cluster. Sometimes, a combination of deprojected and standard model components needs to be used and  I could use get_source() and set_source() to make such combinations. The trick was to setup the deprojected model components first and then use get_source() to parse the model expression into a new set_source() expression where I added an extra model component. 

          dep.set_source("xsmekal")

         set_source(0,get_source(0)+powlaw1d.p1)

         set_source(1,get_source(1)+powlaw1d.p1)     

         set_source(2,get_source(2)+powlaw1d.p2)

         set_source(3,get_source(3)+powlaw1d.p2)

This is a part of my script where I insert the expression for each data set by hand. Of course there must be a nicer way to do this with loops etc. 
The deprojected source expression is set by dep.set_source() for all data sets as explained in the Deproject documentation. Here I used just a standard XSPEC mekal model.  Note that there are 2 annuli and each of them has 2 data sets. I add the same power law model component (p1 and p2) to each annulus. 


get_model() allows you to inspect the full model expression.