.. include:: ../references.rst The low-level Sherpa API ------------------------ Here we repeat our `fit to the 3C273 spectrum `_, but with the low-level API (and without going quite as far in the analysis). Explore the Sherpa Object Model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In a new working directory, download a MAST spectrum of :download:`3C 273 <./3c273.fits>` and start IPython along with the standard imports:: $ ipython --matplotlib import numpy as np import matplotlib.pyplot as plt If you have trouble accessing the spectrum you can download it straight away using Python:: from astropy.extern.six.moves.urllib import request url = 'http://python4astronomers.github.com/_downloads/3c273.fits' open('3c273.fits', 'wb').write(request.urlopen(url).read()) %ls Import a few Sherpa classes needed to characterize a fit:: from sherpa.data import Data1D from sherpa.models import PowLaw1D from sherpa.stats import Chi2DataVar from sherpa.optmethods import LevMar from sherpa.fit import Fit Import the Python FITS reader `astropy.io.fits`_ and open the spectrum as a table:: from astropy.io import fits dat = fits.open('3c273.fits')[1].data Access the `WAVELENGTH` and `FLUX` columns from the pyFITS ``RecArray``. Populate variables represented as ``wave``, ``flux``, and ``err``. Normalize the flux and assume uncertainties of 2% of the flux:: wave = dat.field('WAVELENGTH') flux = dat.field('FLUX') * 1e14 err = dat.field('FLUX') * 0.02e14 Create a Sherpa ``Data1D`` data set from the NumPy arrays ``wave``, ``flux``, and ``err``. The data arrays are accessible from the ``data`` object as the attributes ``x``, ``y``, and ``staterror``:: data = Data1D('3C 273', wave, flux, err) print(data) Array access:: print('x {0}'.format(data.x)) print('y {0}'.format(data.y)) print('err {0}'.format(data.staterror)) Define a convenience function ``plot_data`` that calls the matplotlib functions ``plot`` and ``errorbar`` according to certain criteria. Plot the ``x`` and ``y`` arrays using the format specified in the optional argument, ``fmt``. Clear the plot if the ``clear`` argument is ``True``. Add the plot errorbars if the ``err`` array is present. Plot the spectrum by accessing the NumPy arrays in the Sherpa data set using our new function and its default arguments:: def plot_data(x, y, err=None, fmt='.', clear=True): if clear: plt.clf() plt.plot(x, y, fmt) if err is not None: plt.errorbar(x, y, err, fmt=None, ecolor='b') plot_data(data.x, data.y, data.staterror) .. image:: 3c273_data_mast.png :scale: 75 Create a Sherpa power-law model ``pl``. All Sherpa models maintain a tuple of parameters in ``pars``. Access each of the model's parameter objects and print the ``name`` and ``val`` attributes:: pl = PowLaw1D('pl') pl.pars for par in pl.pars: print('{0} {1}'.format(par.name, par.val)) print(pl) Set the power-law reference to be 4000 Angstroms and print out the ``PowLaw1D`` object and its parameter information. Each model parameter is accessible as an attribute its model. For example, the power-law amplitude is referenced with ``pl.ampl``:: pl.ref = 4000. print(pl) Model parameters are themselves class objects:: print(pl.ampl) .. admonition:: Exercise (for the interested reader): Special methods and properties Wait. Didn't we just set ``pl.ref`` to be an float? How can ``pl.ref`` be an float and a ``Parameter`` object? .. raw:: html

Click to Show/Hide Solution

The answer is that pl.ref is in fact an object, but its model class supports a special setter method ``__setattr__()`` that updates the pl.ref.val attribute underneath. The ``property`` function defines custom getter and setter functions for a particular class attribute:: class Parameter(object): def __init__(self): # private attribute intended to be reference as 'val'. self._value = 1.0 def _get_val(self): return self._value def _set_val(self, value): self._value = value # setup a 'val' attribute val = property(_get_val, _set_val) class Model(object): def __setattr__(self, name, val): if isinstance(getattr(self, name, None), Parameter): getattr(self, name).val = val else: object.__setattr__(self, name, val) def __init__(self): self.ref = Parameter() m = Model() m.ref m.ref = 4 m.ref m.ref.val .. raw:: html
Create a ``Fit`` object made up of a Sherpa data set, model, fit statistic, and optimization method. Fit the spectrum to a power-law with least squares (Levenberg-Marquardt) using the chi-squared statistic with data variance:: f = Fit(data, pl, Chi2DataVar(), LevMar()) result = f.fit() print(result) # or alternatively print(result.format()) Over-plot the fitted model atop the data points using our convenience function ``plot_data``. This time calculate the model using the best-fit parameter values over the ``data.x`` and plot using a custom format and indicate ``clear=False``:: plot_data(data.x, pl(data.x), fmt="-", clear=False) .. image:: 3c273_fit_mast.png :scale: 75