Trying out Sherpa and ChIPS

Let's try running some simple Sherpa/ChIPS routines from the IPython notebook. First we load in the modules and do some set up to make things a bit nicer (there are a bunch of other notebooks - which I should learn how to link to - that explain what this does):

Set up

In [1]:
from sherpa.astro.ui import *
from pycrates import *
from pychips import *
from pychips.hlui import *
In [2]:
# Sherpa settings
import logging
logging.getLogger('sherpa').propagate = False

import sys
sys.tracebacklimit = None
In [3]:
# ChIPS settings
set_preference('window.display', 'false')

from IPython.core.display import Image
import IPython.core.display as display

import tempfile
import os

def my_print(fname=None):
    "Display the current ChIPS window in IPython as a PNG"
    use_temp = False
    if fname is None:
        fname = tempfile.mkstemp(prefix='chips_temp', suffix='.png')[1]

    # Get ChIPS to create a PNG of the current window
    print_window(fname, 'clobber=true')
    img = Image(filename=fname)
    display.display_png(img)
    if use_temp:
        os.remove(fname)
        
def show_gui():
    print("show_gui has been disabled from the notebook")
In [4]:
# HTML display of Sherpa models

from IPython.core.display import HTML

import collections

def parameter_header_raw():
    """Create a TR row for the header items."""
    def th(n): return "<th>{}</th>".format(n)
    return "<tr class='parvals'>" + "".join([th(n) for n in ["Component","Name","Frozen","Value","Min","Max","Units"]]) + "</tr>"

# Note that this code needs to be extended to support properly displaying linked parameters.
def par_to_tr_raw(par, bgcol=False, cptcount=1, newcpt=False):
    """Convert a parameter to a TR.
    
    If bgcol is set True then the model class will be added to the tr
    component (bad parameter name).
    
    The cptcount argument gives the number of parameters there are
    for this component and newcpt is True if it is the first
    parameter of the component.
    """
    def td(n, cname=None): 
        ostr = "<td"
        if cname != None:
            ostr += " class='{}'".format(cname)
        if newcpt and cptcount > 1 and n == "modelname":
            ostr += " rowspan='{}'".format(cptcount)
        return ostr + ">{}</td>".format(getattr(par, n))

    out = "<tr"
    if bgcol:
        out += " class='model'"
    out += ">"
        
    attrs = ["name", "frozen", "val", "min", "max", "units"]
    if newcpt:
        attrs.insert(0, "modelname")
        
    for n in attrs:
        if n == "val" and (par.val <= par.min or par.val >= par.max):
            cname = "atbounds"
        else:
            cname = None
        out += td(n, cname=cname)
    out += "</tr>"
    return out

def model_to_html_raw(mdl):
    """Create HTML display for a model (raw)"""
    def clean(s):
        if s.startswith("(") and s.endswith(")"):
            return s[1:-1]
        else:
            return s
    out = "<table class='model'><caption>" + clean(mdl.name) + "</caption>"
    out += parameter_header_raw()
    
    ctr = collections.Counter([p.modelname for p in mdl.pars])
    
    # have decided to highlight frozen parameters rather than different components
    curmodel = None
    #bgcol = True
    for par in mdl.pars:
        #if par.modelname != curmodel:
        #    bgcol = not bgcol
        #    curmodel = par.modelname
        bgcol = par.frozen    
        out += par_to_tr_raw(par, bgcol=bgcol, cptcount=ctr[par.modelname], newcpt=par.modelname != curmodel)        
        curmodel = par.modelname

    out += "</table>"
    return out

def display_model(mdl):
    "Create a HTML view of the model, displayed by IPython notebook"
    return HTML(model_to_html_raw(mdl))
In [5]:
%%html
<style>
table.model caption { font-size: 120%; }
tr.parvals { background: lightsteelblue; }
/* have decided to highlight the rows of frozen parameters */
tr.model { background: rgb(255,215,0); background: rgba(255,215,0,0.5)}
/* highlight those values at either the min or max boundary */
td.atbounds { background: goldenrod; }
</style>

Ready to go

In this notebook directory I have placed a few files; for instance a source spectrum, along with the response files (ARF and RMF), but no background file.

In [6]:
ls src*
src.arf  src.pi   src.rmf

Loading in an X-ray spectrum into Sherpa

We can load in this spectrum.

In [7]:
load_pha('src.pi')
read ARF file src.arf
read RMF file src.rmf

The show_* commands do not work

Unfortunately, show_data defaults to passing the output via a pager, and this means that it appears in the window from which ipython notebook was run rather than here (i.e. there is no output in the notebook for the show_* family of commands). I have not investigated a work around, but - as shown below - you can view the same data using the various get_* commands.

In [8]:
show_data()

As mentioned, you can select the elements you want to see:

In [9]:
print(get_data())
name           = src.pi
channel        = Float64[1024]
counts         = Float64[1024]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 36610.9103717
backscal       = 0.000716120215634
areascal       = 1.0
grouped        = False
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

Back to Sherpa, and setting up a source model

In [10]:
set_source(xsphabs.gal * xsapec.clus)
print(get_source())
(xsphabs.gal * xsapec.clus)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   gal.nH       thawed            1            0       100000 10^22 atoms / cm^2
   clus.kT      thawed            1        0.008           64        keV
   clus.Abundanc frozen            1            0            5           
   clus.redshift frozen            0            0           10           
   clus.norm    thawed            1            0        1e+24           

Or, taking advantage of the HTML support code we set up earlier:

In [11]:
display_model(get_source())
Out[11]:
xsphabs.gal * xsapec.clus
ComponentNameFrozenValueMinMaxUnits
galnHFalse1.00.0100000.010^22 atoms / cm^2
cluskTFalse1.00.00864.0keV
AbundancTrue1.00.05.0
redshiftTrue0.00.010.0
normFalse1.00.01e+24

As noted above, the show_* commands do not work in the notebook, but you can still view the output of other commands, as shown above and below.

In [12]:
src = get_source()
type(src)
Out[12]:
sherpa.models.model.BinaryOpModel
In [13]:
src
Out[13]:
<BinaryOpModel model instance '(xsphabs.gal * xsapec.clus)'>
In [14]:
gal.nH.val
Out[14]:
1.0
In [15]:
print(gal)
xsphabs.gal
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   gal.nH       thawed            1            0       100000 10^22 atoms / cm^2
In [16]:
print(gal.nH)
val         = 1.0
min         = 0.0
max         = 100000.0
units       = 10^22 atoms / cm^2
frozen      = False
link        = None
default_val = 1.0
default_min = 0.0
default_max = 100000.0

Note that the plot_data command, although creating a ChIPS plot, as shown by the output from info, does not create an inline plot. The my_print command has to be called explicitly.

In [17]:
plot_data()
In [18]:
info()
Out[18]:
Window [win1]
  Frame [frm1]
    Plot [plot1]   (0.15,0.15)  .. (0.90,0.90)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      Curve [crv1]
      X Axis [ax1]
      Y Axis [ay1]
In [19]:
my_print()

Grouping, filtering, ...

In [20]:
group_counts(20)
notice(0.5, 7)
plot_data()
my_print()

We can do a fit; fortunately in this case the fit does not take too long, but as it is being run in a separate process you can still create new notebook "cells", although executing any code will pause until the fit is finished.

In [21]:
fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 3.41269e+09
Final fit statistic   = 242.307 at function evaluation 69
Data points           = 45
Degrees of freedom    = 42
Probability [Q-value] = 5.51961e-30
Reduced statistic     = 5.7692
Change in statistic   = 3.41269e+09
   gal.nH         1.28704     
   clus.kT        1.00295     
   clus.norm      0.000504721 

We can use the get_fit_results command:

In [22]:
get_fit_results()
Out[22]:
<Fit results instance>
In [23]:
print(get_fit_results())
datasets       = (1,)
itermethodname = none
methodname     = levmar
statname       = chi2gehrels
succeeded      = True
parnames       = ('gal.nH', 'clus.kT', 'clus.norm')
parvals        = (1.2870378764829515, 1.0029460946831705, 0.00050472078737771344)
statval        = 242.306537543
istatval       = 3412690743.42
dstatval       = 3412690501.12
numpoints      = 45
dof            = 42
qval           = 5.51961103406e-30
rstat          = 5.76920327484
message        = successful termination
nfev           = 69

How about errors? Hmmm, because this is an exceedingly poor fit ($\chi^2_r > 5$) it can not be run:

In [24]:
conf()
---------------------------------------------------------------------------
EstErr                                    Traceback (most recent call last)
<ipython-input-24-98c7569d8335> in <module>()
----> 1 conf()

<string> in conf(*args)

/export/testlocal/ipython/ciao-4.6/lib/python2.7/site-packages/sherpa/ui/utils.pyc in conf(self, *args)
   6605            get_conf_results, get_proj_results
   6606         """
-> 6607         self._confidence_results = self._est_errors(args, 'confidence')
   6608 
   6609     def proj(self, *args):

/export/testlocal/ipython/ciao-4.6/lib/python2.7/site-packages/sherpa/ui/utils.pyc in _est_errors(self, args, methodname)
   6513             parlist = None
   6514         ids, f = self._get_fit(id, otherids, self._estmethods[methodname])
-> 6515         res = f.est_errors(self._methods, parlist)
   6516         res.datasets = ids
   6517         info(res.format())

/export/testlocal/ipython/ciao-4.6/lib/python2.7/site-packages/sherpa/fit.pyc in run(fit, *args, **kwargs)
     48     def run(fit, *args, **kwargs):
     49         fit.model.startup()
---> 50         result = func(fit, *args, **kwargs)
     51         fit.model.teardown()
     52         return result

/export/testlocal/ipython/ciao-4.6/lib/python2.7/site-packages/sherpa/fit.pyc in est_errors(self, methoddict, parlist)
    949                 #raise FitError('reduced statistic larger than ' +
    950                 #               str(self.estmethod.max_rstat))
--> 951                 raise EstErr( 'rstat>max', str(self.estmethod.max_rstat) )
    952 
    953         # If statistic is chi-squared, change fitting method to

EstErr: reduced statistic larger than 3

And a look at the result. I do not expect it to be that good since I picked the model somewhat at random and no background has been subtracted:

In [25]:
plot_fit()
my_print()

You can use ChIPS/Sherpa commands to change the plot:

In [26]:
log_scale()
set_curve(['line.color', 'orange'])
my_print()

Out of interest, I add in the emission from the xsapec component (i.e. without the absorbing model). It is also displayed at the "native" resolution of the data (i.e. the RMF), rather than being grouped to match the data.

In [27]:
# can not use display_model(clus) because it doesn't get displayed,
# presumably because only one item is processed per cell; I am sure
# that this is a recent change/regression and is hopefully not a design decision
print(clus)
plot_model_component(clus, overplot=True)
set_histogram(['line.color', 'steelblue'])
my_print()
xsapec.clus
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   clus.kT      thawed      1.00295        0.008           64        keV
   clus.Abundanc frozen            1            0            5           
   clus.redshift frozen            0            0           10           
   clus.norm    thawed  0.000504721            0        1e+24