diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..3a0aa48
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,13 @@
+*~
+
+examples/m87
+dist
+docs/_build
+**/api
+
+build
+__pycache__
+.cache/
+
+deproject.egg-info
+.eggs
diff --git a/.hgignore b/.hgignore
deleted file mode 100644
index 7268e63..0000000
--- a/.hgignore
+++ /dev/null
@@ -1,4 +0,0 @@
-examples/m87
-dist
-docs/.build
-
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000..9b0dadf
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,37 @@
+# Based on
+# https://github.com/astrofrog/example-travis-conda/blob/master/.travis.yml
+# but all problems are mine
+#
+
+dist: xenial
+
+language: python
+
+matrix:
+ include:
+ - python: 3.7
+ env: NUMPY=1.15
+
+before_install:
+ - wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda3.sh
+ - chmod +x miniconda3.sh
+ - ./miniconda3.sh -b -p /home/travis/miniconda
+ - export PATH=/home/travis/miniconda/bin:$PATH
+ - conda update --yes conda
+
+ # DEBUG information
+ - conda info -a
+ - python --version
+
+install:
+ - conda create --yes -n test python=$TRAVIS_PYTHON_VERSION
+ # - conda activate test
+ - source activate test
+ - conda install --yes numpy=$NUMPY pytest
+ - conda install --yes -c sherpa sherpa
+
+ # I guess we should actually install the software!
+ - pip install .
+
+script:
+ - py.test
diff --git a/Makefile b/Makefile
deleted file mode 100644
index c52b7ae..0000000
--- a/Makefile
+++ /dev/null
@@ -1,19 +0,0 @@
-WWW = /proj/web-cxc-dmz/htdocs/contrib/deproject
-
-.PHONY: m87tar docs dist
-
-dist:
- rm -rf dist
- python setup.py sdist
-
-m87tar:
- tar zcvf examples/m87.tar.gz examples/m87
-
-docs:
- cd docs; \
- make html
-
-install:
- rsync -av docs/.build/html/ $(WWW)/
- rsync -av examples/m87.tar.gz $(WWW)/downloads/
- cp dist/deproject-*.tar.gz $(WWW)/downloads/deproject.tar.gz
diff --git a/README.rst b/README.rst
index 9ed7fd0..5e27975 100644
--- a/README.rst
+++ b/README.rst
@@ -1,14 +1,13 @@
Deproject
=========
-``Deproject`` is a CIAO Sherpa extension package to facilitate
-deprojection of two-dimensional annular X-ray spectra to recover the
-three-dimensional source properties. For typical thermal models this would
-include the radial temperature and density profiles. This basic method
-has been used extensively for X-ray cluster analysis and is the basis for the
-XSPEC model ``projct``. The ``deproject`` module brings this
-functionality to *Sherpa* as a Python module that is straightforward to use and
-understand.
+``Deproject`` is a `Sherpa `_ extension package
+to facilitate deprojection of two-dimensional annular X-ray spectra to recover
+the three-dimensional source properties. For typical thermal models this would
+include the radial temperature and density profiles. This basic method has been
+used extensively for X-ray cluster analysis and is the basis for the XSPEC
+model ``projct``. The ``deproject`` module brings this functionality to
+*Sherpa* as a Python module that is straightforward to use and understand.
The basic physical assumption of ``deproject`` is that the extended source
emissivity is constant and optically thin within spherical shells whose radii
@@ -16,6 +15,10 @@ correspond to the annuli used to extract the specta. Given this assumption one
constructs a model for each annular spectrum that is a linear volume-weighted
combination of shell models.
+Version 0.2 of ``deproject`` is limited to circular annuli.
+
+Further documentation is available at https://deproject.readthedocs.io/
+
License
-------
@@ -23,3 +26,101 @@ The ``deproject`` module is released under the
`BSD 2-Clause license `_,
available as the file ``LICENSE`` in the source distribution.
+Requirements
+------------
+
+The installation assumes that you are installing ``deproject`` into
+the `CIAO environment `_ (CIAO 4.11 or
+later), since this is the easiest way to get the XSPEC models along
+with Sherpa. The `standalone Sherpa `_
+version can be used, but in this case you will need to `build Sherpa
+with XSPEC support
+`_.
+
+The following Python packages are required:
+
+ - sherpa
+ - `Astropy `_ (restricted to version 3.0 when
+ using CIAO 4.11)
+ - `SciPy `_.
+
+Installation
+------------
+
+The ``deproject`` module should install with the following command
+(assuming CIAO 4.11 is already installed):
+
+ echo "numpy==1.12.1" > constraints.txt
+ pip install -c constraints.txt 'astropy<3.1' deproject
+
+Example
+-------
+
+If you have a set of X-ray PHA spectra called src.pi, where is
+an integer representing the annulus number, and the files contain the
+``XFLT0001`` to ``XFLT0005`` header keywords used by the
+`XSPEC projct model `_,
+then a
+`Deproject object `_
+can be created using the
+`deproject_from_xflt `_
+helper routine with the commands:
+
+ >>> from deproject import deproject_from_xflt
+ >>> from astropy import units as u
+ >>> dep = deproject_from_xflt('src*.pi', 0.492 * u.arcsec)
+
+where, in this example, the ``XFLT0001`` and ``XFLT0002`` keywords,
+which specify the inner and outer radii of the annulus, are in
+ACIS pixels, and so need to be multiplied by 0.492 arcseconds to
+convert to an angle (the second parameter).
+
+This will automatically load the spectra into separate Sherpa datasets,
+which *can* be fitted individually, but it is generally easier to use
+the object returned by ``deproject_from_xflt``. For instance, the
+following will set the data range to be fit for *each* spectra and ensure
+that the background is subtracted before fitting:
+
+ >>> dep.ignore(None, 0.5)
+ >>> dep.ignore(7.0, None)
+ >>> dep.subtract()
+
+Sherpa functions are used to change the statistic and optimiser:
+
+ >>> from sherpa.astro import ui
+ >>> ui.set_stat('chi2xspecvar')
+ >>> ui.set_method('levmar')
+
+The data can be fit, and errors estimated for all the parameter, using
+the onion-skin deprojection approach, with the following commands:
+
+ >>> onion = dep.fit()
+ >>> errs = dep.conf()
+
+The return value includes the density (and errors, if appropriate), as
+an `Astropy Quantity `_.
+
+ >>> print(onion['density'])
+ print(onion['density'])
+ density
+ 1 / cm3
+ --------------------
+ 0.1100953546292787
+ 0.07736622021374819
+ 0.04164827967805805
+ 0.03630168106524076
+ 0.025221797991301052
+ 0.021845331641349316
+ ...
+ 0.012396857131392835
+ 0.01336640115325031
+ 0.012303975980575187
+ 0.013631563529090736
+ 0.013996131292837352
+ 0.010843683594144967
+ 0.023067220584935984
+ Length = 20 rows
+
+The `on-line documentation `_
+contains more information, including creating the ``Deproject`` object
+directly (without the need for the ``XFLTxxxx`` keywords).
diff --git a/cosmocalc.py b/cosmocalc.py
deleted file mode 100644
index 02323d2..0000000
--- a/cosmocalc.py
+++ /dev/null
@@ -1,258 +0,0 @@
-#!/usr/bin/env python
-"""
-Calculate useful values for a given cosmology. This module uses code adapted
-from `CC.py`_ (`James Schombert`_) which is a Python version of the
-`Cosmology Calculator`_ (`Ned Wright`_).
-
-The following values are calculated:
-
- ==== =================================== ===========
- Name Value Units
- ==== =================================== ===========
- z Input redshift
- H0 Hubble constant
- WR Omega(radiation)
- WK Omega curvaturve = 1-Omega(total)
- WM Omega matter
- WV Omega vacuum
- DTT Time from z to now Gyr
- age Age of Universe Gyr
- zage Age of Universe at redshift z Gyr
- DCMR Comoving radial distance Gyr Mpc cm
- VCM Comoving volume within redshift Gpc3
- DA Angular size distance Gyr Mpc cm
- DL Luminosity distance Gyr Mpc cm
- PS Plate scale - distance per arcsec kpc cm
- ==== =================================== ===========
-
-.. _`James Schombert`: http://abyss.uoregon.edu/~js/
-.. _`CC.py`: http://www.astro.ucla.edu/~wright/CC.python
-.. _`Ned Wright`: http://www.astro.ucla.edu/~wright/intro.html
-.. _`Cosmology Calculator`: http://www.astro.ucla.edu/~wright/CosmoCalc.html
-
-:Copyright: Smithsonian Astrophysical Observatory (2009)
-:Author: Tom Aldcroft (aldcroft@head.cfa.harvard.edu)
-"""
-
-import math
-
-# Define a few constants
-cm_per_pc = 3.0856775813057289536e+18
-c = 299792.458 # velocity of light in km/sec
-km_per_ly = 3600*24*365.25*c # km per light-year
-Tyr = 977.8 # coefficent for converting 1/H into Gyr
-arcsec_per_rad = 206264.806
-_outvals_str = ('z H0 WM WV WK WR',
- 'DA DA_Gyr DA_Mpc DA_cm',
- 'DL DL_Gyr DL_Mpc DL_cm',
- 'DCMR DCMR_Gyr DCMR_Mpc DCMR_cm',
- 'PS_kpc PS_cm',
- 'DTT DTT_Gyr',
- 'VCM VCM_Gpc3',
- 'age age_Gyr',
- 'zage zage_Gyr',)
-_outvals = (' '.join(_outvals_str)).split()
-
-def cosmocalc(z, H0=71, WM=0.27, WV=None):
- """
- Calculate useful values for the supplied cosmology.
-
- This routine returns a dictionary of values in the form ``: ``,
- where the values are supplied in "natural" units for cosmology, e.g. 1/H0.
- In addition various useful unit conversions are done and stored in the
- dictionary as ``_: ``. E.g. angular size distance::
-
- 'DA': 0.38250549415474988,
- 'DA_Gyr': 5.2678010166833023,
- 'DA_Mpc': 1615.1022857909447,
- 'DA_cm': 4.9836849147807571e+27
-
- Example::
-
- >>> from cosmocalc import cosmocalc
- >>> from pprint import pprint
- >>> pprint(cosmocalc(3, H0=75, WM=.25))
- {'DA': 0.39103776375786625,
- 'DA_Gyr': 5.0980896720325548,
- 'DA_Mpc': 1563.0689649039205,
- 'DA_cm': 4.8231268630387788e+27,
- 'DCMR': 1.564151055031465,
- 'DCMR_Gyr': 20.392358688130219,
- 'DCMR_Mpc': 6252.2758596156818,
- 'DCMR_cm': 1.9292507452155115e+28,
- 'DL': 6.25660422012586,
- 'DL_Gyr': 81.569434752520877,
- 'DL_Mpc': 25009.103438462727,
- 'DL_cm': 7.717002980862046e+28,
- 'DTT': 0.84826379084317027,
- 'DTT_Gyr': 11.059097795819358,
- 'H0': 75,
- 'PS_cm': 2.3383178917293232e+22,
- 'PS_kpc': 7.5779721961095019,
- 'VCM': 1.2756009121294902,
- 'VCM_Gpc3': 1023.7714254161302,
- 'WK': 0.0,
- 'WM': 0.25,
- 'WR': 7.4044444444444448e-05,
- 'WV': 0.74992595555555552,
- 'age': 1.0133755371756261,
- 'age_Gyr': 13.211714670004362,
- 'z': 3,
- 'zage': 0.16511174633245579,
- 'zage_Gyr': 2.1526168741850036}
-
- :param z: redshift
- :param H0: Hubble constant (default = 71)
- :param WM: Omega matter (default = 0.27)
- :param WV: Omega vacuum (default = 1.0 - WM - 0.4165/(H0*H0))
-
- :rtype: dictionary of cosmology values (name_unit = value)
- """
-
- if z > 100:
- z = z / 299792.458 # Values over 100 are in km/s
-
- if WV is None:
- WV = 1.0 - WM - 0.4165/(H0*H0) # Omega(vacuum) or lambda
- age = 0.0 # age of Universe in units of 1/H0
-
- h = H0 / 100.
- WR = 4.165E-5 / (h*h) # includes 3 massless neutrino species, T0 = 2.72528
- WK = 1 - WM - WR - WV
- az = 1.0 / (1 + 1.0*z)
- n=1000 # number of points in integrals
- for i in range(n):
- a = az * (i + 0.5) / n
- adot = math.sqrt(WK + (WM/a) + (WR/(a*a)) + (WV*a*a))
- age = age + 1./adot
-
- zage = az * age / n
- DTT = 0.0
- DCMR = 0.0
-
- # do integral over a=1/(1+z) from az to 1 in n steps, midpoint rule
- for i in range(n):
- a = az + (1-az) * (i+0.5) / n
- adot = math.sqrt(WK + (WM/a) + (WR/(a*a)) + (WV*a*a))
- DTT = DTT + 1./adot
- DCMR = DCMR + 1./(a*adot)
-
- DTT = (1.-az) * DTT / n
- DCMR = (1.-az) * DCMR / n
- age = DTT + zage
-
- # tangential comoving distance
- ratio = 1.0
- x = math.sqrt(abs(WK)) * DCMR
- if x > 0.1:
- if WK > 0:
- ratio = 0.5 * (math.exp(x) - math.exp(-x)) / x
- else:
- ratio = math.math.sin(x) / x
- else:
- y = x * x
- if WK < 0:
- y = -y
- ratio = 1. + y/6. + y*y/120.
- DCMT = ratio * DCMR
-
- # comoving volume computation
- ratio = 1.00
- x = math.sqrt(abs(WK)) * DCMR
- if x > 0.1:
- if WK > 0:
- ratio = (0.125 * (math.exp(2.*x) - math.exp(-2.*x)) - x/2.) / (x**3 / 3.)
- else:
- ratio = (x/2. - math.sin(2.*x)/4.)/(x**3 / 3.)
- else:
- y = x * x
- if WK < 0: y = -y
- ratio = 1. + y/5. + (2./105.)*y*y
- VCM = ratio * DCMR**3 / 3.
- VCM_Gpc3 = 4. * math.pi * (0.001*c/H0)**3 * VCM
-
- DA = az * DCMT
- DL = DA / (az*az)
-
- # Now convert to some more useful units
- Gyr = lambda x: Tyr / H0 * x
- Mpc = lambda x: c / H0 * x
- cm = lambda x: Mpc(x) * 1e6 * cm_per_pc
-
- DA_Gyr = Gyr(DA)
- DA_Mpc = Mpc(DA)
- DA_cm = cm(DA)
-
- DL_Gyr = Gyr(DL)
- DL_Mpc = Mpc(DL)
- DL_cm = cm(DL)
-
- DCMR_Gyr = Gyr(DCMR)
- DCMR_Mpc = Mpc(DCMR)
- DCMR_cm = cm(DCMR)
-
- DTT_Gyr = Gyr(DTT)
- age_Gyr = Gyr(age)
- zage_Gyr = Gyr(zage)
-
- PS_kpc = Mpc(DA) * 1000 / arcsec_per_rad
- PS_cm = PS_kpc * cm_per_pc * 1000
-
- localvals = locals()
- return dict((x, localvals[x]) for x in _outvals)
-
-def get_options():
- """
- cosmocalc.py [options] redshift [name_unit [name_unit2 ...]]
-
- Allowed ``name_unit`` values::
-
- DA DA_Gyr DA_Mpc DA_cm
- DL DL_Gyr DL_Mpc DL_cm
- DCMR DCMR_Gyr DCMR_Mpc DCMR_cm
- PS_kpc PS_cm
- DTT DTT_Gyr
- VCM VCM_Gpc3
- age age_Gyr
- zage zage_Gyr
- H0 WM WV WK WR z
-
- If no ``name_unit`` values are supplied then all the above will be printed."""
- from optparse import OptionParser
-
- parser = OptionParser(get_options.__doc__)
- parser.set_defaults()
- parser.add_option("--H0",
- default=None,
- type='float',
- help="Hubble constant")
- parser.add_option("--WM",
- default=None,
- type='float',
- help="")
- parser.add_option("--WV",
- default=None,
- type='float',
- help="")
- opt, args = parser.parse_args()
- return opt, args, parser
-
-def main():
- opt, args, parser = get_options()
-
- if len(args) < 1:
- parser.error('Need a redshift')
-
- kwargs = dict((key, val) for (key, val) in opt.__dict__.items() if val is not None)
- z = float(args[0])
- cc = cosmocalc(z, **kwargs)
- try:
- outlines = []
- for outkey in (args[1:] or _outvals):
- outlines.append(outkey + ' = ' + str(cc[outkey]))
- print '\n'.join(outlines)
- except KeyError:
- parser.error(outkey + ' is not a valid output name_unit')
-
-if __name__ == '__main__':
- main()
diff --git a/deproject.py b/deproject.py
deleted file mode 100644
index b6e34c8..0000000
--- a/deproject.py
+++ /dev/null
@@ -1,270 +0,0 @@
-"""
-Deproject from a set of 2-d annular spectra to the 3-d object properties.
-
-:Copyright: Smithsonian Astrophysical Observatory (2009)
-:Author: Tom Aldcroft (aldcroft@head.cfa.harvard.edu)
-"""
-import re
-from math import pi, sqrt
-import numpy
-import specstack
-from cosmocalc import cosmocalc, arcsec_per_rad
-import sherpa.astro.ui as SherpaUI
-
-class Deproject(specstack.SpecStack):
- """
- Deproject from a set of 2-d annular spectra to the 3-d object properties.
-
- The ``radii`` parameter must be a list of values that starts with the inner
- radius of the inner annulus and includes each radius up through the outer
- radius of the outer annulus. Thus the ``radii`` list will be one element
- longer than the number of annuli.
-
- :param radii: sorted list of circular annulus radii (arcsec) for extracted spectra
- :param theta: azimuthal extent of annuli (degrees) (default = 360)
- :param angdist: angular size distance (cm)
- """
- def __init__(self, radii, theta=360, angdist=None):
- if len(radii) < 2:
- raise ValueError('radii parameter must be a list with at least two values')
- self.radii = radii
- self.nshell = len(radii)-1
- self._theta = theta
- self._angdist = angdist
- self._redshift = None
- super(Deproject, self).__init__()
-
- def _get_redshift(self):
- if self._redshift is None:
- self._redshift = self.find_parval('redshift')
- return self._redshift
-
- def _set_redshift(self, redshift):
- self._redshift = redshift
-
- def _get_angdist(self):
- if self._angdist is None:
- cc = cosmocalc(self.redshift)
- self._angdist = cc['DA_cm']
- return self._angdist
-
- def _set_angdist(self, angdist):
- self._angdist = angdist
-
- redshift = property(_get_redshift, _set_redshift, None, "Source redshift")
- angdist = property(_get_angdist, _set_angdist, None, "Angular size distance")
-
- def _calc_vol_norm(self):
- """
- Calculate the normalized volumes of cylindrical annuli intersecting with
- spherical shells. Sets ``self.vol_norm`` to volume[i,j] / v_sphere
- giving the normalized volume of i'th shell intersecting with j'th
- annulus (where indexes start at 0). V_sphere = 4/3 pi r^3 where r is
- the outer radius of the outer shell.
-
- :rtype: None
- """
- r = self.radii
- theta_rad = self._theta / 180. * pi
- cv = numpy.zeros([self.nshell, self.nshell])
- v = numpy.zeros([self.nshell, self.nshell])
- for a, ra0 in enumerate(r[:-1]): # Annulus
- ra1 = r[a+1]
- for s, rs0 in enumerate(r[:-1]): # Spherical shell
- rs1 = r[s+1]
- if s >= a:
- # Volume of cylindrical annulus (ra0,ra1) intersecting sphere (rs1)
- cv[s,a] = 2.0 * theta_rad / 3 * ((rs1**2 - ra0**2)**1.5 - (rs1**2 - ra1**2)**1.5)
-
- # Volume of annulus (ra0,ra1) intersecting spherical shell (rs0,rs1)
- v[s,a] = cv[s,a]
- if s - a > 0:
- v[s,a] -= cv[s-1,a]
- self.vol_norm = v / (4. * pi / 3. * r[-1]**3)
-
- def _create_src_model_components(self):
- """
- Create source model components for each shell corresponding to the
- source model expression.
- """
- # Find the generic components in source model expression
- RE_model = re.compile(r'\b \w+ \b', re.VERBOSE)
- for match in RE_model.finditer(self.srcmodel):
- model_type = match.group()
- self.srcmodel_comps.append(dict(type=model_type,
- start=match.start(),
- end=match.end()))
-
- # For each shell create the corresponding model components so they can
- # be used later to create composite source models for each dataset
- for shell in range(self.nshell):
- for srcmodel_comp in self.srcmodel_comps:
- model_comp = {}
- model_comp['type'] = srcmodel_comp['type']
- model_comp['name'] = '%s_%d' % (model_comp['type'], shell)
- model_comp['shell'] = shell
- SherpaUI.create_model_component(model_comp['type'], model_comp['name'])
- model_comp['object'] = eval(model_comp['name']) # Work-around in lieu of accessor
- self.model_comps.append(model_comp)
-
- def set_source(self, srcmodel='xsphabs*xsapec'):
- """
- Create a source model for each dataset. A dataset is associated
- with a specific extraction annulus.
-
- :param srcmodel: string expression defining source model
- :rtype: None
- """
- self.srcmodel = srcmodel
- self._calc_vol_norm()
- self._create_src_model_components()
-
- for dataset in self.datasets:
- dataid = dataset['id']
- annulus = dataset['annulus']
- modelexprs = []
- for shell in range(annulus, self.nshell):
- srcmodel = self.srcmodel
- for model_comp in reversed(self.srcmodel_comps):
- i0 = model_comp['start']
- i1 = model_comp['end']
- model_comp_name = '%s_%d' % (model_comp['type'], shell)
- srcmodel = srcmodel[:i0] + model_comp_name + srcmodel[i1:]
- modelexprs.append('%.5f * %s' % (self.vol_norm[shell, annulus], srcmodel))
-
- modelexpr = " + ".join(modelexprs)
- print 'Setting source model for dataset %d = %s' % (dataid, modelexpr)
- SherpaUI.set_source(dataid, modelexpr)
-
- def set_bkg_model(self, bkgmodel):
- """
- Create a source model for each dataset. A dataset is associated
- with a specific extraction annulus.
-
- :param bkgmodel: string expression defining background model
- :rtype: None
- """
- self.bkgmodel = bkgmodel
-
- bkg_norm = {}
- for obsid in self.obsids:
- bkg_norm_name = 'bkg_norm_%d' % obsid
- print 'Creating model component xsconstant.%s' % bkg_norm_name
- SherpaUI.create_model_component('xsconstant', bkg_norm_name)
- bkg_norm[obsid] = eval(bkg_norm_name) # Uggh, don't know proper model accessor
-
- for dataset in self.datasets:
- print 'Setting bkg model for dataset %d to bkg_norm_%d' % (dataset['id'], dataset['obsid'])
- SherpaUI.set_bkg_model(dataset['id'], bkg_norm[dataset['obsid']] * bkgmodel)
-
- def fit(self):
- """
- Do a fit of the model parameters using the "onion-peeling" method:
-
- - First fit the outside shell model using the outer annulus spectrum
- - Freeze the model parameters for the outside shell
- - Fit the next inward shell / annulus and freeze those parameters
- - Repeat until all datasets have been fit and all shell parameters determined.
- - Return model parameters to original thawed/frozen status
-
- :rtype: None
- """
- thawed = [] # Parameter objects that are not already frozen
- for annulus in reversed(range(self.nshell)):
- dataids = [x['id'] for x in self.datasets if x['annulus'] == annulus]
- print 'Fitting', dataids
- SherpaUI.fit(*dataids)
- for model_comp in self.model_comps:
- name = model_comp['name']
- if model_comp['shell'] == annulus:
- # Remember parameters that are currently thawed
- for par in [SherpaUI.get_par('%s.%s'%(name, x))
- for x in SherpaUI.get_model_pars(name)]:
- if not par.frozen:
- thawed.append(par)
- print 'Freezing', model_comp['name']
- SherpaUI.freeze(model_comp['name'])
-
- # Unfreeze parameters
- for par in thawed:
- print 'Thawing', par.fullname
- par.thaw()
-
- def get_density(self):
- """
- Get electron density (cm^-3) for each shell using the standard
- definition of normalization for Xspec thermal models::
-
- n_e = sqrt(norm * 4*pi * DA^2 * 1e14 * (1+z)^2 / volume * ne_nh_ratio))
-
- norm = model normalization from sherpa fit
- DA = angular size distance (cm)
- volume = volume (cm^3)
- ne_nh_ratio = 1.18
-
- Note that the model components for each volume element (intersection of the
- annular cylinder ``a`` with the spherical shell ``s``) are multiplied by a volume
- normalization::
-
- vol_norm[s,a] = volume[s,a] / v_sphere
- v_sphere = volume of sphere enclosing outer annulus
-
- With this convention the ``volume`` used in calculating the electron density
- is simply ``v_sphere``.
-
- :rtype: numpy array of densities (cm^-3) corresponding to shells
- """
-
- ne_nh_ratio = 1.18 # Electron to proton ratio n_e/n_p
- DA_cm = self.angdist
- r_sphere = self.radii[-1] / arcsec_per_rad * DA_cm
- volume = 4 * pi / 3 * r_sphere**3 # volume of sphere enclosing outer shell (cm^3)
- z = self.redshift
-
- dens = []
- for shell in range(self.nshell):
- norm = self.find_norm(shell)
- dens.append(sqrt(norm * 4 * pi * DA_cm**2 * 1e14 * (1.0+z)**2 / volume * ne_nh_ratio))
-
- return numpy.array(dens)
-
- def conf(self):
- """
- Run conf on each of the model parameters using the "onion-peeling" method:
-
- - First conf the outside shell model using the outer annulus spectrum
- - Freeze the model parameters for the outside shell
- - get confidences for the next inward shell / annulus and freeze those parameters
- - Repeat until all datasets have been conf()-ed and all shell-by-shell error parameters determined.
- - Return model parameters to original thawed/frozen status
- - WARNING: This ignores the correlations between parameters
-
- :rtype: None
- """
- thawed = [] # Parameter objects that are not already frozen
- conf_results = []
- this_conf_result = []
- for annulus in reversed(range(self.nshell)):
- dataids = [x['id'] for x in self.datasets if x['annulus'] == annulus]
- print 'Getting shell-by-shell confidence for dataset ', dataids
- SherpaUI.conf(*dataids)
- this_conf_result = SherpaUI.get_conf_results()
- conf_results.insert(0, this_conf_result)
- for model_comp in self.model_comps:
- name = model_comp['name']
- if model_comp['shell'] == annulus:
- # Remember parameters that are currently thawed
- for par in [SherpaUI.get_par('%s.%s'%(name, x))
- for x in SherpaUI.get_model_pars(name)]:
- if not par.frozen:
- thawed.append(par)
- print 'Freezing', model_comp['name']
- SherpaUI.freeze(model_comp['name'])
-
- # Unfreeze parameters
- for par in thawed:
- print 'Thawing', par.fullname
- par.thaw()
-
- return conf_results
-
diff --git a/docs/_templates/autosummary/base.rst b/docs/_templates/autosummary/base.rst
new file mode 100644
index 0000000..a58aa35
--- /dev/null
+++ b/docs/_templates/autosummary/base.rst
@@ -0,0 +1,10 @@
+{% if referencefile %}
+.. include:: {{ referencefile }}
+{% endif %}
+
+{{ objname }}
+{{ underline }}
+
+.. currentmodule:: {{ module }}
+
+.. auto{{ objtype }}:: {{ objname }}
diff --git a/docs/_templates/autosummary/class.rst b/docs/_templates/autosummary/class.rst
new file mode 100644
index 0000000..85105fa
--- /dev/null
+++ b/docs/_templates/autosummary/class.rst
@@ -0,0 +1,65 @@
+{% if referencefile %}
+.. include:: {{ referencefile }}
+{% endif %}
+
+{{ objname }}
+{{ underline }}
+
+.. currentmodule:: {{ module }}
+
+.. autoclass:: {{ objname }}
+ :show-inheritance:
+
+ {% if '__init__' in methods %}
+ {% set caught_result = methods.remove('__init__') %}
+ {% endif %}
+
+ {% block attributes_summary %}
+ {% if attributes %}
+
+ .. rubric:: Attributes Summary
+
+ .. autosummary::
+ {% for item in attributes %}
+ ~{{ name }}.{{ item }}
+ {%- endfor %}
+
+ {% endif %}
+ {% endblock %}
+
+ {% block methods_summary %}
+ {% if methods %}
+
+ .. rubric:: Methods Summary
+
+ .. autosummary::
+ {% for item in methods %}
+ ~{{ name }}.{{ item }}
+ {%- endfor %}
+
+ {% endif %}
+ {% endblock %}
+
+ {% block attributes_documentation %}
+ {% if attributes %}
+
+ .. rubric:: Attributes Documentation
+
+ {% for item in attributes %}
+ .. autoattribute:: {{ item }}
+ {%- endfor %}
+
+ {% endif %}
+ {% endblock %}
+
+ {% block methods_documentation %}
+ {% if methods %}
+
+ .. rubric:: Methods Documentation
+
+ {% for item in methods %}
+ .. automethod:: {{ item }}
+ {%- endfor %}
+
+ {% endif %}
+ {% endblock %}
diff --git a/docs/_templates/autosummary/module.rst b/docs/_templates/autosummary/module.rst
new file mode 100644
index 0000000..11208a2
--- /dev/null
+++ b/docs/_templates/autosummary/module.rst
@@ -0,0 +1,41 @@
+{% if referencefile %}
+.. include:: {{ referencefile }}
+{% endif %}
+
+{{ objname }}
+{{ underline }}
+
+.. automodule:: {{ fullname }}
+
+ {% block functions %}
+ {% if functions %}
+ .. rubric:: Functions
+
+ .. autosummary::
+ {% for item in functions %}
+ {{ item }}
+ {%- endfor %}
+ {% endif %}
+ {% endblock %}
+
+ {% block classes %}
+ {% if classes %}
+ .. rubric:: Classes
+
+ .. autosummary::
+ {% for item in classes %}
+ {{ item }}
+ {%- endfor %}
+ {% endif %}
+ {% endblock %}
+
+ {% block exceptions %}
+ {% if exceptions %}
+ .. rubric:: Exceptions
+
+ .. autosummary::
+ {% for item in exceptions %}
+ {{ item }}
+ {%- endfor %}
+ {% endif %}
+ {% endblock %}
diff --git a/docs/changes.rst b/docs/changes.rst
new file mode 100644
index 0000000..7d1c1f2
--- /dev/null
+++ b/docs/changes.rst
@@ -0,0 +1,145 @@
+
+.. include:: references.rst
+
+*******
+Changes
+*******
+
+.. _changes_020:
+
+Version 0.2.0
+=============
+
+Overview
+--------
+
+The code has been updated to run with Python 3 and can now be
+installed from `PyPI`_. Documentation has been moved to
+`Read The Docs `_.
+
+The deproject module now requires Astropy_, which can be
+`installed within CIAO 4.11 `_. The three main areas where Astropy functionality is
+used are:
+
+- the use of `Astropy Quantity `_
+ values (both for arguments to methods and returned values);
+- `Astropy Data Tables `_
+ are used to return tabular data;
+- and Cosmology calculations now use the `Astropy cosmology
+ module `_ rather than
+ the `cosmocalc` module.
+
+The :py:func:`~deproject.deproject.deproject_from_xflt` helper function
+has been introduced, which uses the ``XFLT0001`` to ``XFLT0005``
+keywords in the input files to determine the annulus parameters (radii
+and covering angle). The covering angle (:math:`\theta`) can now vary per
+annulus.
+
+Error values can now be generated using the onion-peeling approach,
+for the confidence and covariance methods, and the values are returned
+as an Astropy Table. Parameter values can now be tied together (to
+combine annuli to try and avoid "ringing"). There is improved support
+for accessing and plotting values.
+
+Details
+-------
+
+The code has been re-arranged into the ``deproject`` package, which
+means that you really should say ``from deproject.deproject import
+Deproject``, but the ``deproject`` module re-exports
+:py:mod:`deproject.deproject` so that existing scripts still work, and
+you do not not have to type in the same word multiple times! The
+package has been updated so that it is available on `PyPI`_.
+
+The scaling between shells (calculated from the intersection between
+spheres and cylinders) was limited to 5 decimal places, which could
+cause problems with certain choices of annuli (such as an annulus
+making no contribution to interior annuli). This restriction has been
+removed.
+
+Added support for per-annulus ``theta`` values (that is, each annulus
+can have a different opening angle). The ``radii``, ``theta``, and
+``angdist`` parameters to :py:class:`~deproject.deproject.Deproject`
+all now require values that is an
+`Astropy quantity `_ rather
+than a dimensionless value.
+
+Added the :py:func:`~deproject.deproject.deproject_from_xflt` helper
+function, which creates a :py:class:`~deproject.deproject.Deproject`
+instance from PHA files which contain the XSPEC XFLT0001 to XFLT0005
+keywords (as used by the `projct`_ model), rather than specifying the
+values from the command line. The routine will error out if the
+keywords indicate elliptical annuli, and the default is to assume the
+radii are in arcseconds, but a scaling factor can be given if the
+radii are in some other units (such as pixels).
+
+Added the :py:meth:`~deproject.deproject.Deproject.guess` method to do
+an initial fit to each annulus, following the approach suggested in
+the `XSPEC`_ documentation for `projct`_, by just fitting the
+individual (not de-projected) models to each annulus. This can help
+speed up the deproject fit -
+:py:meth:`~deproject.deproject.Deproject.fit` - as well as help avoid
+the fit getting stuck in a local minimum.
+
+Added :py:meth:`~deproject.deproject.Deproject.covar` and
+:py:meth:`~deproject.deproject.Deproject.conf` methods that estimate
+errors - using the covariance and confidence methods respectively -
+using the onion-skin model (i.e. the errors on the outer annuli are
+evaluated, then this component is frozen and the errors on the next
+annulus are evaluated).
+
+The :py:meth:`~deproject.deproject.Deproject.fit`,
+:py:meth:`~deproject.deproject.Deproject.conf`, and
+:py:meth:`~deproject.deproject.Deproject.covar` methods now all return
+Astropy Tables containing the results per annulus. These values can
+also be retrieved with the
+:py:meth:`~deproject.deproject.Deproject.get_fit_results`,
+:py:meth:`~deproject.deproject.Deproject.get_conf_results`, or
+:py:meth:`~deproject.deproject.Deproject.get_covar_results` methods. A
+number of columns (radii and density) are returned as Astropy
+quantities.
+
+The ``cosmocalc`` module has been removed and the `Astropy cosmology
+module `_ is used
+instead. This is only used if the angular-diameter distance to the
+source is calculated rather than explicitly given. The default
+cosmology is now set to `Planck15
+`_.
+
+Values, as a function of radius, can be plotted with a number of new
+methods: :py:meth:`~deproject.deproject.Deproject.fit_plot`,
+:py:meth:`~deproject.deproject.Deproject.conf_plot`, and
+:py:meth:`~deproject.deproject.Deproject.covar_plot` display the last
+fit results (with the last two including error estimates), and the
+:py:meth:`~deproject.deproject.Deproject.par_plot` and
+:py:meth:`~deproject.deproject.Deproject.density_plot` methods show
+the current values. These support a number of options, including
+switching between angular and physical distances for the radii.
+
+The :py:meth:`~deproject.deproject.Deproject.get_shells` method has
+been added to make it easy to see which annuli are combined together,
+and the :py:meth:`~deproject.deproject.Deproject.get_radii` method to
+find the radii of the annuli (in a range of units).
+
+Added the :py:meth:`~deproject.deproject.SpecStack.tie_par` and
+:py:meth:`~deproject.deproject.SpecStack.untie_par` methods to make it
+easy to tie (or untie) parameters in neighbouring annuli. The
+onion-skin approach - used when fitting or running an error analysis -
+recognizes annuli that are tied together and fits these
+simultaneously, rather than individually.
+
+The :py:meth:`~deproject.deproject.Deproject.set_source` method can
+now be called multiple times (previously it would lead to an error).
+
+Added error checking for several routines, such as
+:py:meth:`~deproject.deproject.SpecStack.thaw` when given an
+unknown parameter name.
+
+Updated to support Python 3.5 and to have better support when the
+``pylab`` backend is selected. Support for the ChIPS backend is
+limited. A basic test suite has been added.
+
+Version 0.1.0
+=============
+
+Initial version.
diff --git a/docs/conf.py b/docs/conf.py
index 3d1d225..9b953d8 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -11,22 +11,42 @@
# All configuration values have a default; values that are commented out
# serve to show the default.
-import sys, os
+import os
+import sys
-# If your extensions are in another directory, add it here. If the directory
-# is relative to the documentation root, use os.path.abspath to make it
-# absolute, like shown here.
-sys.path.append(os.path.abspath('../'))
+import sphinx_rtd_theme
+
+# Try and let Sherpa be found.
+#
+if os.path.split(os.getcwd())[1] == 'docs':
+ sys.path.insert(0, os.path.abspath('../src/'))
+
+on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
# General configuration
# ---------------------
+# sphinx.ext.napoleon needs 1.3, although it is unlikely any modern system has
+# anything as old that this is important.
+#
+needs_sphinx = '1.3'
+
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.autodoc']
+extensions = [
+ 'sphinx.ext.autodoc',
+ 'sphinx.ext.autosummary',
+ 'sphinx.ext.inheritance_diagram',
+ 'sphinx.ext.mathjax',
+ 'sphinx.ext.napoleon'
+]
+
+napoleon_google_docstring = False
+
+autosummary_generate = True
# Add any paths that contain templates here, relative to this directory.
-templates_path = ['.templates']
+templates_path = ['_templates']
# The suffix of source filenames.
source_suffix = '.rst'
@@ -39,16 +59,16 @@
# General information about the project.
project = u'deproject'
-copyright = u'2009, Tom Aldcroft'
+copyright = u'2009, 2019 Tom Aldcroft'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
-version = '0.1'
+version = '0.2'
# The full version, including alpha/beta/rc tags.
-release = '0.1.3'
+release = '0.2.0'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
@@ -84,14 +104,30 @@
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
+# Graphviz based on values from AstroPy - see
+# https://github.com/astropy/sphinx-astropy/blob/master/sphinx_astropy/conf/v1.py
+#
+graphviz_output_format = "svg"
+
+graphviz_dot_args = [
+ '-Nfontsize=10',
+ '-Nfontname=Helvetica Neue, Helvetica, Arial, sans-serif',
+ '-Efontsize=10',
+ '-Efontname=Helvetica Neue, Helvetica, Arial, sans-serif',
+ '-Gfontsize=10',
+ '-Gfontname=Helvetica Neue, Helvetica, Arial, sans-serif'
+]
+
# Options for HTML output
# -----------------------
-# The style sheet to use for HTML and HTML Help pages. A file of that name
-# must exist either in Sphinx' static/ path, or in one of the custom paths
-# given in html_static_path.
-html_style = 'default.css'
+# html_theme = 'alabaster'
+html_theme = 'sphinx_rtd_theme'
+
+# html_theme_options = {}
+# html_theme_path = []
+html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
# The name for this set of Sphinx documents. If None, it defaults to
# " v documentation".
@@ -112,11 +148,16 @@
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['.static']
+#
+# html_static_path = ['_static']
+
+# html_extra_path = []
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
# using the given strftime format.
-#html_last_updated_fmt = '%b %d, %Y'
+
+# I believe this is what AstroPy is doing
+html_last_updated_fmt = '%d %b %Y'
# If true, SmartyPants will be used to convert quotes and dashes to
# typographically correct entities.
diff --git a/docs/cosmocalc.rst b/docs/cosmocalc.rst
deleted file mode 100644
index 43b5aff..0000000
--- a/docs/cosmocalc.rst
+++ /dev/null
@@ -1,18 +0,0 @@
-:mod:`cosmocalc`
-===================
-
-.. automodule:: cosmocalc
-
-Functions
-----------
-
-.. autofunction:: cosmocalc
-
-.. autofunction:: get_options
-
-.. autofunction:: main
-
-
-
-
-
diff --git a/docs/deproject.rst b/docs/deproject.rst
deleted file mode 100644
index 6ec6de9..0000000
--- a/docs/deproject.rst
+++ /dev/null
@@ -1,17 +0,0 @@
-:mod:`deproject`
-==================
-
-.. include:: references.rst
-
-.. automodule:: deproject
-
-Deproject class
-----------------
-
-.. autoclass:: Deproject
- :show-inheritance:
- :members:
- :inherited-members:
- :undoc-members:
-
-
diff --git a/docs/examples/3c186.rst b/docs/examples/3c186.rst
new file mode 100644
index 0000000..94aea99
--- /dev/null
+++ b/docs/examples/3c186.rst
@@ -0,0 +1,56 @@
+
+.. include:: ../references.rst
+
+Multiple datasets per annulus (3C186)
+=====================================
+
+A second example illustrates the use of :mod:`deproject` for a multi-obsid
+observation of 3C186. It also shows how to set a background model for fitting
+with the ``cstat`` statistic. The extracted spectral data for this example are
+not yet publicly available, but were used in `Siemiginowska et al. 2010`_:
+
+.. _`Siemiginowska et al. 2010`: https://ui.adsabs.harvard.edu/#abs/2010ApJ...722..102S
+
+The script starts with some setup::
+
+ >>> import deproject
+
+ >>> radii = ('2.5', '6', '17')
+ >>> dep = deproject.Deproject(radii=[float(x) for x in radii])
+
+ >>> set_method("levmar")
+ >>> set_stat("cstat")
+
+Now we read in the data as before with ``dep.load_pha()``. The only
+difference to the :ref:`M87 example `
+is the additional loop over the obsids. The ``dep.load_pha()`` function
+automatically extracts the obsid (the identifier used to
+discriminate Chandra observations) from the file header. This is used later in
+the case of setting a background model.
+::
+
+ >>> obsids = (9407, 9774, 9775, 9408)
+ >>> for ann in range(len(radii) -1):
+ ... for obsid in obsids:
+ ... dep.load_pha('3c186/%d/ellipse%s-%s.pi' % (obsid, radii[ann], radii[ann+1]), annulus=ann)
+
+Create and configure the source model expression as usual::
+
+ >>> dep.set_source('xsphabs*xsapec')
+ >>> dep.ignore(None, 0.5)
+ >>> dep.ignore(7, None)
+ >>> dep.freeze("xsphabs.nh")
+
+ >>> dep.set_par('xsapec.redshift', 1.06)
+ >>> dep.set_par('xsphabs.nh', 0.0564)
+
+Set the background model (here we use the IPython ``%run`` magic command
+to evaluate the commands in the file ``acis-s-bkg.py``)::
+
+ >>> %run acis-s-bkg.py
+ >>> acis_s_bkg = get_bkg_source()
+ >>> dep.set_bkg_model(acis_s_bkg)
+
+Fit the projection model::
+
+ >>> dep.fit()
diff --git a/docs/examples/m87.rst b/docs/examples/m87.rst
new file mode 100644
index 0000000..e312a70
--- /dev/null
+++ b/docs/examples/m87.rst
@@ -0,0 +1,499 @@
+
+.. include:: ../references.rst
+
+***
+M87
+***
+
+Now we step through in detail the ``fit_m87.py`` script in the ``examples``
+directory to explain each step and illustrate how to use the :mod:`deproject`
+module. This script should serve as the template for doing your own analysis.
+
+This example uses extracted spectra, response products, and analysis results
+for the Chandra observation of M87 (obsid 2707). These were kindly provided
+by Paul Nulsen. Results based on this observation can be found in
+`Forman et al 2005`_
+and via the CXC Archive `Obsid 2707 Publications`_ list.
+
+.. _`Forman et al 2005`: https://ui.adsabs.harvard.edu/?#abs/2005ApJ...635..894F
+.. _`Obsid 2707 Publications`: https://cda.harvard.edu/chaser/viewerContents.do?obsid=2707&operation=ads
+
+The examples assume this is being run directly from the Sherpa shell,
+which has already imported the Sherpa module. If you are using IPython
+or a Jupyter notebook then the following command is needed:
+
+ >>> from sherpa.astro.ui import *
+
+Set up
+======
+
+The first step is to load in the Deproject class::
+
+ >>> from deproject import Deproject
+
+and then set a couple of constants (note that both the angular-diameter
+distance and thta values *must* be
+`Astropy quantites `_)::
+
+ >>> from astropy import units as u
+ >>> redshift = 0.004233 # M87 redshift
+ >>> angdist = 16 * u.Mpc # M87 distance
+ >>> theta = 75 * u.deg # Covering angle of sectors
+
+Next we create an array of the the annular radii in arcsec. The
+`numpy.arange`_ method here returns an array from 30 to 640 in steps of 30.
+These values were in pixels in the original spectral extraction so we
+multiply by the pixel size (0.492 arcseconds) to convert to an angle.::
+
+ >>> radii = numpy.arange(30., 640., 30) * 0.492 * u.arcesc
+
+The ``radii`` parameter must be a list of values that starts with the inner
+radius of the inner annulus and includes each radius up through the outer
+radius of the outer annulus. Thus the ``radii`` list will be one element
+longer than the number of annuli.
+
+::
+
+ >>> print(radii)
+ [ 14.76 29.52 44.28 59.04 73.8 88.56 103.32 118.08 132.84 147.6
+ 162.36 177.12 191.88 206.64 221.4 236.16 250.92 265.68 280.44 295.2
+ 309.96] arcsec
+
+*Now the key step* of creating the :class:`Deproject` object ``dep``. This
+object is the interface to the all the :mod:`deproject`
+methods used for the deprojection analysis.
+::
+
+ >>> dep = Deproject(radii, theta=theta, angdist=angdist)
+
+If you are not familiar with object oriented programming, the ``dep`` object is
+just a thingy that stores all the information about the deprojection analysis
+(e.g. the source redshift, PHA file information and the source model
+definitions) as object *attributes*. It also has object *methods*
+(i.e. functions) you can call such as ``dep.get_par(parname)`` or
+``dep.load_pha(file)``. The full list of attributes and methods are in the
+:mod:`deproject` module documentation.
+
+In this particular analysis the spectra were extracted from a 75
+degree sector of the annuli, hence ``theta`` is set in the object
+initialization, where the units are set explicitly using the
+`Astropy support for units `_.
+Note that this parameter only needs to be set if any of the
+annuli are sectors rather than the full circle.
+Since the redshift is **not** a good distance
+estimator for M87 we also explicitly set the angular size distance
+using the ``angdist`` parameter.
+
+.. _m87_load:
+
+Load the data
+=============
+
+Now load the PHA spectral files for each annulus using the Python ``range``
+function to loop over a sequence ranging from 0 to the last annulus. The
+:py:meth:`~deproject.deproject.Deproject.load_pha`
+call is the first example of a :mod:`deproject` method
+(i.e. function) that mimics a *Sherpa* function with the same name. In this
+case ``dep.load_pha(file, annulus)`` loads the PHA file using the
+*Sherpa* `load_pha`_
+function but also registers the dataset in the spectral stack::
+
+ >>> for annulus in range(len(radii) - 1):
+ ... dep.load_pha('m87/r%dgrspec.pha' % (annulus + 1), annulus)
+
+.. note::
+ The ``annulus`` parameter is required in ``dep.load_pha()`` to allow
+ multiple data sets per annulus, such as repeated Chandra observations
+ or different XMM instruments.
+
+Create the model
+================
+
+With the data loaded we set the source model for each of the spherical shells
+with the
+:py:meth:`~deproject.deproject.Deproject.set_source`
+method. This is one of the more complex bits of
+:mod:`deproject`. It automatically generates all the model components for each
+shell and then assigns volume-weighted linear combinations of those components
+as the source model for each of the annulus spectral datasets::
+
+ >>> dep.set_source('xswabs * xsmekal')
+
+The model expression can be any valid *Sherpa* model expression with the
+following caveats:
+
+ - Only the generic model type should be specified in the expression. In
+ typical *Sherpa* usage one generates the model component name in the
+ model expression, e.g. ``set_source('xswabs.abs1 * xsmekal.mek1')``. This
+ would create model components named ``abs1`` and ``mek1``. In
+ ``dep.set_source()`` the model component names are auto-generated as
+ ``_``.
+ - Only one of each model type can be used in the model expression. A source
+ model expression like ``"xsmekal + gauss1d + gauss1d"`` would result in an
+ error due to the model component auto-naming.
+
+Next any required parameter values are set and their `freeze`_ or `thaw`_
+status are set.
+::
+
+ >>> dep.set_par('xswabs.nh', 0.0255)
+ >>> dep.freeze('xswabs.nh')
+
+ >>> dep.set_par('xsmekal.abundanc', 0.5)
+ >>> dep.thaw('xsmekal.abundanc')
+
+ >>> dep.set_par('xsmekal.redshift', redshift)
+
+As a convenience if any of the model components have a
+``redshift`` parameter that value will be used as the default redshift for
+calculating the angular size distance.
+
+.. _data_selection:
+
+Define the data to fit
+======================
+
+Now the energy range used in the fitting is restricted using the stack version
+of the *Sherpa* `ignore`_ command. The `notice`_ command is also available.
+::
+
+ >>> dep.ignore(None, 0.5)
+ >>> dep.ignore(1.8, 2.2)
+ >>> dep.ignore(7, None)
+
+Define the optimiser and statistic
+==================================
+
+At this point the model is completely set up and we are ready to do the initial
+"onion-peeling" fit. As for normal high-signal fitting with binned spectra we
+issue the commands to set the optimization method and the fit statistic.
+
+In this case we use the Levenberg-Marquardt optimization method with
+the :math:`\chi^2` fit statistic, where the variance is estimated
+using a similar approach to XSPEC.
+
+::
+
+ >>> set_method("levmar")
+ >>> set_stat("chi2xspecvar")
+ >>> dep.subtract()
+
+.. note::
+ The ``chi2xspecvar`` statistic is used since the values will be
+ :ref:`compared to results from XSPEC `.
+
+Seeding the fit
+===============
+
+We first try to "guess" a good starting point for the fit (since the
+default fit parameters, in particular the normalization, are often
+not close to the expected value). In this case the
+:py:meth:`~deproject.deproject.Deproject.guess`
+method implements a scheme suggested in the `projct`_ documentation, which
+fits each annulus with an un-projected model and then corrects the
+normalizations for the shell/annulus overlaps:
+
+ >>> dep.guess()
+ Projected fit to annulus 19 dataset: 19
+ Dataset = 19
+ ...
+ Change in statistic = 4.60321e+10
+ xsmekal_19.kT 2.748 +/- 0.0551545
+ xsmekal_19.Abundanc 0.429754 +/- 0.0350466
+ xsmekal_19.norm 0.00147471 +/- 2.45887e-05
+ Projected fit to annulus 18 dataset: 18
+ Dataset = 18
+ ...
+ Change in statistic = 4.61945e+10
+ xsmekal_18.kT 2.68875 +/- 0.0543974
+ xsmekal_18.Abundanc 0.424198 +/- 0.0335045
+ xsmekal_18.norm 0.00147443 +/- 2.43454e-05
+ ...
+ Projected fit to annulus 0 dataset: 0
+ Dataset = 0
+ ...
+ Change in statistic = 2.15043e+10
+ xsmekal_0.kT 1.57041 +/- 0.0121572
+ xsmekal_0.Abundanc 1.05816 +/- 0.0373779
+ xsmekal_0.norm 0.00152496 +/- 2.93297e-05
+
+As shown in the screen output, the ``guess`` routine fits each
+annulus in turn, from outer to inner. After each fit, the
+normalization (the ``xsmekal_*.norm`` terms) are adjusted by
+the filling factor of the shell.
+
+.. note::
+ The ``guess`` step is *optional*. Parameter values can also be set,
+ either for an individual annulus with the Sherpa `set_par`_ function,
+ such as::
+
+ set_par('xsmekal_0.kt', 1.5)
+
+ or to all annuli with the Deproject method
+ :py:meth:`~deproject.deproject.Deproject.set_par`::
+
+ dep.set_par('xsmekal.abundanc', 0.3)
+
+The data can be inspected using normal Sherpa commands. The
+following shows the results of the guess for dataset 0,
+which corresponds to the inner-most annulus (the
+:py:attr:`~deproject.deproject.Deproject.datasets` attribute
+can be used to map between annuus and dataset identifier).
+
+::
+
+ >>> set_xlog()
+ >>> plot_fit(0)
+
+.. image:: m87_ann0_guess.png
+
+The overall shape of the model looks okay, but the normalization is
+not (since it has been adjusted by the volume-filling factor of the
+intersection of the annulus and shell). The gap in the data around 2
+keV is because we :ref:`explicitly excluded this range
+` from the fit.
+
+.. note::
+ The :py:class:`~deproject.deproject.Deproject` class also provides
+ methods for plotting each annulus in a separate window, such as
+ :py:meth:`~deproject.deproject.Deproject.plot_data`
+ and
+ :py:meth:`~deproject.deproject.Deproject.plot_fit`.
+
+Fitting the data
+================
+
+The
+:py:meth:`~deproject.deproject.Deproject.fit`
+method performs the full onion-skin deprojection (the output
+is similar to the ``guess`` method, with the addition of
+parameters being frozen as each annulus is processed)::
+
+ >>> onion = dep.fit()
+ Fitting annulus 19 dataset: 19
+ Dataset = 19
+ ...
+ Reduced statistic = 1.07869
+ Change in statistic = 1.29468e-08
+ xsmekal_19.kT 2.748 +/- 0.0551534
+ xsmekal_19.Abundanc 0.429758 +/- 0.0350507
+ xsmekal_19.norm 0.249707 +/- 0.00416351
+ Freezing xswabs_19
+ Freezing xsmekal_19
+ Fitting annulus 18 dataset: 18
+ Dataset = 18
+ ...
+ Reduced statistic = 1.00366
+ Change in statistic = 13444.5
+ xsmekal_18.kT 2.41033 +/- 0.274986
+ xsmekal_18.Abundanc 0.375824 +/- 0.136926
+ xsmekal_18.norm 0.0551815 +/- 0.00441495
+ Freezing xswabs_18
+ Freezing xsmekal_18
+ ...
+ Freezing xswabs_1
+ Freezing xsmekal_1
+ Fitting annulus 0 dataset: 0
+ Dataset = 0
+ ...
+ Reduced statistic = 2.72222
+ Change in statistic = 12115.6
+ xsmekal_0.kT 1.63329 +/- 0.0278325
+ xsmekal_0.Abundanc 1.1217 +/- 0.094871
+ xsmekal_0.norm 5.7067 +/- 0.262766
+ Change in statistic = 13699.6
+ xsmekal_0.kT 1.64884 +/- 0.0242336
+ xsmekal_0.Abundanc 1.14629 +/- 0.0903398
+ xsmekal_0.norm 5.68824 +/- 0.245884
+ ...
+
+The fit to the central annulus (dataset 0) now looks like::
+
+ >>> plot_fit_delchi(0)
+
+.. image:: m87_ann0_fit.png
+
+After the fit process each shell model has an association normalization that
+can be used to calculate the densities. This is where the source angular
+diameter distance is used. If the angular diameter distance is not set
+explicitly in the original ``dep = Deproject(...)`` command then it is
+calculated automatically from the source redshift and an assumed Cosmology,
+which if not set is taken to be the
+`Planck15 model `_
+provided by the
+`astropy.cosmology `_
+module.
+
+One can examine the values being used as follows (note that the
+:py:attr:`~deproject.deproject.Deproject.angdist` attribute is an
+`Astropy Quantity `_)::
+
+ >>> print("z={:.5f} angdist={}".format(dep.redshift, dep.angdist))
+ z=0.00423 angdist=16.0 Mpc
+
+New in version 0.2.0 is the behavior of the ``fit`` method, which now
+returns an Astropy table which tabulates the fit results as a function
+of annulus. This *includes* the electron density, which can also be
+retrieved with
+:py:meth:`~deproject.deproject.Deproject.get_density`.
+The fit results can also be retrieved with the
+:py:meth:`~deproject.deproject.Deproject.get_fit_results`
+method.
+
+::
+
+ >>> print(onion)
+ annulus rlo_ang rhi_ang ... xsmekal.norm density
+ arcsec arcsec ... 1 / cm3
+ ------- ------- ------- ... ------------------- --------------------
+ 0 14.76 29.52 ... 5.6882389946381275 0.1100953546292787
+ 1 29.52 44.28 ... 2.8089409208987233 0.07736622021374819
+ 2 44.28 59.04 ... 0.814017947154132 0.04164827967805805
+ 3 59.04 73.8 ... 0.6184339453006981 0.03630168106524076
+ 4 73.8 88.56 ... 0.2985327157676323 0.025221797991301052
+ 5 88.56 103.32 ... 0.22395312678845017 0.021845331641349316
+ ... ... ... ... ... ...
+ 13 206.64 221.4 ... 0.07212121560592619 0.012396857131392835
+ 14 221.4 236.16 ... 0.08384338967492334 0.01336640115325031
+ 15 236.16 250.92 ... 0.07104455447410102 0.012303975980575187
+ 16 250.92 265.68 ... 0.08720295254137615 0.013631563529090736
+ 17 265.68 280.44 ... 0.09192970392746878 0.013996131292837352
+ 18 280.44 295.2 ... 0.05518150227052414 0.010843683594144967
+ 19 295.2 309.96 ... 0.24970680803387219 0.023067220584935984
+ Length = 20 rows
+
+Inspecting the results
+======================
+
+The electron density can be retrieved directly from the fit results,
+with the
+:py:meth:`~deproject.deproject.Deproject.get_density`
+method::
+
+ >>> print(dep.get_density())
+ [ 0.11009535 0.07736622 0.04164828 0.03630168 0.0252218 0.02184533
+ 0.01525456 0.01224972 0.01942528 0.01936785 0.01905983 0.01568478
+ 0.01405426 0.01239686 0.0133664 0.01230398 0.01363156 0.01399613
+ 0.01084368 0.02306722] 1 / cm3
+
+or plotted using
+:py:meth:`~deproject.deproject.Deproject.density_plot`::
+
+ >>> dep.density_plot()
+
+.. image:: m87_density.png
+
+.. _m87_ringing:
+
+The temperature profile from the deprojection can be plotted using
+:py:meth:`~deproject.deproject.Deproject.par_plot`::
+
+ >>> dep.par_plot('xsmekal.kt')
+
+.. image:: m87_temperature.png
+
+The unphysical temperature oscillations seen here highlights a known issue
+with this analysis method (e.g. `Russell, Sanders, & Fabian 2008`_).
+
+It can also be instructive to look at various results from the fit,
+such as the reduced statistic for each annulus (as will be shown
+below, there are
+:py:meth:`~deproject.deproject.Deproject.fit_plot`,
+:py:meth:`~deproject.deproject.Deproject.conf_plot`,
+and
+:py:meth:`~deproject.deproject.Deproject.covar_plot` variants)::
+
+ >>> dep.fit_plot('rstat')
+
+.. image:: m87_rstat.png
+
+Error analysis
+==============
+
+Errors can also be calculated, on both the model parameters and the
+derived densities, with either the
+:py:meth:`~deproject.deproject.Deproject.conf`
+or
+:py:meth:`~deproject.deproject.Deproject.covar`
+methods. These
+apply the confidence and covariance error-estimation routines from
+Sherpa using the same onion-skin model used for the fit, and are new
+in version 0.2.0. In this example we use the covariance version, since
+it is generally faster, but confidence is the recommended routine as it
+is more robust (and calculates asymmetric errors)::
+
+ >>> errs = dep.covar()
+
+As with the ``fit`` method, both ``conf`` and ``covar`` return the results
+as an Astropy table. These can also be retrieved with the
+:py:meth:`~deproject.deproject.Deproject.get_conf_results`
+or
+:py:meth:`~deproject.deproject.Deproject.get_covar_results`
+methods. The columns depend on the command (i.e. fit or error results)::
+
+ >>> print(errs)
+ annulus rlo_ang rhi_ang ... density_lo density_hi
+ arcsec arcsec ... 1 / cm3 1 / cm3
+ ------- ------- ------- ... ----------------------- ----------------------
+ 0 14.76 29.52 ... -0.0023299300992805777 0.0022816336635322065
+ 1 29.52 44.28 ... -0.0015878693875514133 0.0015559288091097495
+ 2 44.28 59.04 ... -0.0014852395638492444 0.0014340671599212054
+ 3 59.04 73.8 ... -0.0011539611188859725 0.00111839214283211
+ 4 73.8 88.56 ... -0.001166653616914693 0.001115024462355424
+ 5 88.56 103.32 ... -0.0010421595234548324 0.0009946565126391325
+ ... ... ... ... ... ...
+ 13 206.64 221.4 ... -0.0005679816551559976 0.0005430748019680798
+ 14 221.4 236.16 ... -0.00048083556526315463 0.00046412880973027357
+ 15 236.16 250.92 ... -0.0004951659664768401 0.00047599490797040934
+ 16 250.92 265.68 ... -0.0004031089363366082 0.0003915259159180066
+ 17 265.68 280.44 ... -0.0003647519140328147 0.00035548459401609986
+ 18 280.44 295.2 ... nan nan
+ 19 295.2 309.96 ... -0.00019648511719753958 0.00019482554589523096
+ Length = 20 rows
+
+.. note::
+ With this set of data, the covariance method failed to calculate errors
+ for the parameters of the last-but one shell, which is indicated by
+ the presence of `NaN` values in the error columns.
+
+The fit or error results can be plotted as a function of radius with the
+:py:meth:`~deproject.deproject.Deproject.fit_plot`,
+:py:meth:`~deproject.deproject.Deproject.conf_plot`,
+and
+:py:meth:`~deproject.deproject.Deproject.covar_plot`
+methods (the latter two include error bars). For example, the
+following plot mixes these plots with Matplotlib commands to
+compare the temperature and abundance profiles::
+
+ >>> plt.subplot(2, 1, 1)
+ >>> dep.covar_plot('xsmekal.kt', clearwindow=False)
+ >>> plt.xlabel('')
+ >>> plt.subplot(2, 1, 2)
+ >>> dep.covar_plot('xsmekal.abundanc', clearwindow=False)
+ >>> plt.subplots_adjust(hspace=0.3)
+
+.. image:: m87_temperature_abundance.png
+
+The derived density profile, along with errors, can also be
+displayed (the X axis is displayed using angular units rather
+than as a length)::
+
+ >>> dep.covar_plot('density', units='arcmin')
+
+.. image:: m87_density_errs.png
+
+.. _m87_compare:
+
+Comparing results
+=================
+
+In the images below the :mod:`deproject` results (blue) are compared with
+values (gray boxes) from an independent onion-peeling analysis by P. Nulsen
+using a custom perl script to generate `XSPEC`_ model definition and fit
+commands. These plots were created with the ``plot_m87.py`` script in
+the ``examples`` directory. The agreement is good (note that the version
+of XSPEC used for the two cases does not match):
+
+.. image:: m87_compare_density.png
+.. image:: m87_compare_temperature.png
+
diff --git a/docs/examples/m87_ann0_fit.png b/docs/examples/m87_ann0_fit.png
new file mode 100644
index 0000000..d4b5ac9
Binary files /dev/null and b/docs/examples/m87_ann0_fit.png differ
diff --git a/docs/examples/m87_ann0_guess.png b/docs/examples/m87_ann0_guess.png
new file mode 100644
index 0000000..285ca82
Binary files /dev/null and b/docs/examples/m87_ann0_guess.png differ
diff --git a/docs/examples/m87_compare_density.png b/docs/examples/m87_compare_density.png
new file mode 100644
index 0000000..8359ffe
Binary files /dev/null and b/docs/examples/m87_compare_density.png differ
diff --git a/docs/examples/m87_compare_temperature.png b/docs/examples/m87_compare_temperature.png
new file mode 100644
index 0000000..f6d36fa
Binary files /dev/null and b/docs/examples/m87_compare_temperature.png differ
diff --git a/docs/examples/m87_density.png b/docs/examples/m87_density.png
new file mode 100644
index 0000000..a80eae7
Binary files /dev/null and b/docs/examples/m87_density.png differ
diff --git a/docs/examples/m87_density_errs.png b/docs/examples/m87_density_errs.png
new file mode 100644
index 0000000..ffa5cac
Binary files /dev/null and b/docs/examples/m87_density_errs.png differ
diff --git a/docs/examples/m87_rstat.png b/docs/examples/m87_rstat.png
new file mode 100644
index 0000000..18d756e
Binary files /dev/null and b/docs/examples/m87_rstat.png differ
diff --git a/docs/examples/m87_temperature.png b/docs/examples/m87_temperature.png
new file mode 100644
index 0000000..5836120
Binary files /dev/null and b/docs/examples/m87_temperature.png differ
diff --git a/docs/examples/m87_temperature_abundance.png b/docs/examples/m87_temperature_abundance.png
new file mode 100644
index 0000000..bb0cff8
Binary files /dev/null and b/docs/examples/m87_temperature_abundance.png differ
diff --git a/docs/examples/m87_temperature_tied_comparison.png b/docs/examples/m87_temperature_tied_comparison.png
new file mode 100644
index 0000000..e52b8d1
Binary files /dev/null and b/docs/examples/m87_temperature_tied_comparison.png differ
diff --git a/docs/examples/tie.rst b/docs/examples/tie.rst
new file mode 100644
index 0000000..cf1a534
--- /dev/null
+++ b/docs/examples/tie.rst
@@ -0,0 +1,179 @@
+
+.. include:: ../references.rst
+
+****************
+Combining shells
+****************
+
+We now look at tie-ing parameters from different annuli (that is,
+forcing the temperature or abundance of one shell to be the
+same as a neighbouring shell), which is a technique used to
+try and reduce the "ringing" that can be seen in deprojected
+data (e.g. `Russell, Sanders, & Fabian 2008`_ and seen in
+the :ref:`temperature profile of the M87 data `).
+Sherpa supports complicated links between model parameters
+with the `link`_ and `unlink`_ functions, and the
+Deproject class provides helper methods -
+:py:meth:`~deproject.deproject.Deproject.tie_par`
+and
+:py:meth:`~deproject.deproject.Deproject.untie_par`
+- that allows you to easily tie together annuli.
+
+Starting from where the :doc:`m87` example left off,
+we have:
+
+.. image:: m87_temperature_abundance.png
+
+The :py:meth:`~deproject.deproject.Deproject.get_shells`
+method lists the current set of annuli, and show that there
+are currently no grouped (or tied-together) annuli:
+
+ >>> dep.get_shells()
+ [{'annuli': [19], 'dataids': [19]},
+ {'annuli': [18], 'dataids': [18]},
+ {'annuli': [17], 'dataids': [17]},
+ {'annuli': [16], 'dataids': [16]},
+ {'annuli': [15], 'dataids': [15]},
+ {'annuli': [14], 'dataids': [14]},
+ {'annuli': [13], 'dataids': [13]},
+ {'annuli': [12], 'dataids': [12]},
+ {'annuli': [11], 'dataids': [11]},
+ {'annuli': [10], 'dataids': [10]},
+ {'annuli': [9], 'dataids': [9]},
+ {'annuli': [8], 'dataids': [8]},
+ {'annuli': [7], 'dataids': [7]},
+ {'annuli': [6], 'dataids': [6]},
+ {'annuli': [5], 'dataids': [5]},
+ {'annuli': [4], 'dataids': [4]},
+ {'annuli': [3], 'dataids': [3]},
+ {'annuli': [2], 'dataids': [2]},
+ {'annuli': [1], 'dataids': [1]},
+ {'annuli': [0], 'dataids': [0]}]
+
+For this example we are going to see what happens when we
+tie together the temperature and abundance of annulus 17 and 18
+(the inner-most annulus is numbered 0, the outer-most is 19 in this
+example).
+
+ >>> dep.tie_par('xsmekal.kt', 17, 18)
+ Tying xsmekal_18.kT to xsmekal_17.kT
+ >>> dep.tie_par('xsmekal.abundanc', 17, 18)
+ Tying xsmekal_18.Abundanc to xsmekal_17.Abundanc
+
+Looking at the first few elements of the list returned by `get_shells`
+shows us that the two annuli have now been grouped together. This
+means that when a fit is run then these two annuli will be
+fit simultaneously.
+
+ >>> dep.get_shells()[0:5]
+ [{'annuli': [19], 'dataids': [19]},
+ {'annuli': [17, 18], 'dataids': [17, 18]},
+ {'annuli': [16], 'dataids': [16]},
+ {'annuli': [15], 'dataids': [15]},
+ {'annuli': [14], 'dataids': [14]}]
+
+Since the model values for the 18th annulus have been changed, the
+data needs to be re-fit. The screen output is slightly different,
+as highlighted below, to reflect this new grouping:
+
+ >>> dep.fit()
+ Note: annuli have been tied together
+ Fitting annulus 19 dataset: 19
+ Dataset = 19
+ ...
+ Freezing xswabs_19
+ Freezing xsmekal_19
+ Fitting annuli [17, 18] datasets: [17, 18]
+ Datasets = 17, 18
+ Method = levmar
+ Statistic = chi2xspecvar
+ Initial fit statistic = 469.803
+ Final fit statistic = 447.963 at function evaluation 16
+ Data points = 439
+ Degrees of freedom = 435
+ Probability [Q-value] = 0.323566
+ Reduced statistic = 1.0298
+ Change in statistic = 21.8398
+ xsmekal_17.kT 2.64505 +/- 0.0986876
+ xsmekal_17.Abundanc 0.520943 +/- 0.068338
+ xsmekal_17.norm 0.0960039 +/- 0.00370363
+ xsmekal_18.norm 0.0516096 +/- 0.00232468
+ Freezing xswabs_17
+ Freezing xsmekal_17
+ Freezing xswabs_18
+ Freezing xsmekal_18
+ Fitting annulus 16 dataset: 16
+ ...
+
+The table returned by :py:meth:`~deproject.deproject.Deproject.fit`
+and :py:meth:`~deproject.deproject.Deproject.get_fit_results`
+still contains 20 rows, since each annulus has its own row:
+
+ >>> tied = dep.get_fit_results()
+ >>> len(tied)
+ 20
+
+Looking at the outermost annuli, we can see that the temperature
+and abundance values are the same:
+
+ >>> print(tied['xsmekal.kT'][16:])
+ xsmekal.kT
+ ------------------
+ 2.5470334673573936
+ 2.645045455995055
+ 2.645045455995055
+ 2.7480050952727346
+ >>> print(tied['xsmekal.Abundanc'][16:])
+ xsmekal.Abundanc
+ -------------------
+ 0.24424232427701525
+ 0.5209430839965413
+ 0.5209430839965413
+ 0.4297577517591726
+
+.. note::
+ The field names in the table are case sensitive - that is, you
+ have to say ``xsmekal.kT`` and not ``xsmekal.kt`` - unlike the
+ Sherpa parameter interface.
+
+The results can be compared to the original, such as in the
+following plot, which shows that the temperature profile hasn't
+changed significantly in the core, but some of the oscillations
+at large radii have been reduced (the black circles show the
+temperature values from the original fit):
+
+ >>> dep.fit_plot('xsmekal.kt', units='pc')
+ >>> rmid = (onion['rlo_phys'] + onion['rhi_phys']) / 2
+ >>> plt.scatter(rmid.to(u.pc), onion['xsmekal.kT', c='k')
+
+.. image:: m87_temperature_tied_comparison.png
+
+.. note::
+ As the X-axis of the plot created by
+ :py:meth:`~deproject.deproject.Deproject.fit_plot` is in pc, the
+ ``rmid`` value is expliclty converted to these units before
+ plotting.
+
+ The ``results`` parameter of ``fit_plot`` could have been set
+ to the variable ``onion`` to display the results of the previous fit,
+ but there is no easy way to change the color of the points.
+
+More shells could be tied together to try and further
+reduce the oscillations, or errors calculated, with either
+:py:meth:`~deproject.deproject.Deproject.conf`
+or
+:py:meth:`~deproject.deproject.Deproject.covar`.
+
+The parameter ties can be removed with
+:py:meth:`~deproject.deproject.Deproject.untie_par`:
+
+ >>> dep.untie_par('xsmekal.kt', 18)
+ Untying xsmekal_18.kT
+ >>> dep.untie_par('xsmekal.abundanc', 18)
+ Untying xsmekal_18.Abundanc
+ >>> dep.get_shells()[0:4]
+ [{'annuli': [19], 'dataids': [19]},
+ {'annuli': [18], 'dataids': [18]},
+ {'annuli': [17], 'dataids': [17]},
+ {'annuli': [16], 'dataids': [16]}]
+
diff --git a/docs/index.rst b/docs/index.rst
index 28c7976..32763a5 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,14 +1,11 @@
-.. deproject documentation master file, created by sphinx-quickstart on Sat Jan 31 15:06:12 2009.
- You can adapt this file completely to your liking, but it should at least
- contain the root `toctree` directive.
.. include:: references.rst
Deproject
-====================
+=========
:mod:`Deproject` is a `CIAO`_ `Sherpa`_ extension package to facilitate
-deprojection of two-dimensional annular X-ray spectra to recover the
+deprojection of two-dimensional circular annular X-ray spectra to recover the
three-dimensional source properties. For typical thermal models this would
include the radial temperature and density profiles. This basic method
has been used extensively for X-ray cluster analysis and is the basis for the
@@ -16,396 +13,34 @@ has been used extensively for X-ray cluster analysis and is the basis for the
functionality to *Sherpa* as a Python module that is straightforward to use and
understand.
-The :mod:`deproject` module uses :mod:`specstack` to allow for manipulation of
-a stack of related input datasets and their models. Most of the functions
-resemble ordinary *Sherpa* commands (e.g. `set_par`_, `set_source`_, `ignore`_)
-but operate on a stack of spectra.
-
-The basic physical assumption of :mod:`deproject` is that the extended source
-emissivity is constant and optically thin within spherical shells whose radii
-correspond to the annuli used to extract the specta. Given this assumption one
-constructs a model for each annular spectrum that is a linear volume-weighted
-combination of shell models. The geometry is illustrated in the figure below
-(which would be rotated about the line to the observer in three-dimensions):
-
-.. image:: geometry.png
-
-Model creation
-----------------
-It is assumed that prior to starting :mod:`deproject` the user has extracted
-source and background spectra for each annulus. By convention the annulus
-numbering starts from the inner radius at 0 and corresponds to the dataset
-``id`` used within *Sherpa*. It is not required that the annuli include the
-center but they must be contiguous between the inner and outer radii.
-
-Given a spectral model ``M[s]`` for each shell ``s``, the source model for
-dataset ``a`` (i.e. annulus ``a``) is given by the sum over ``s >= a`` of
-``vol_norm[s,a] * M[s]`` (normalized volume * shell model). The image above
-shows shell 3 in blue and annulus 2 in red. The intersection of (purple) has a
-physical volume defined as ``vol_norm[3,2] * v_sphere`` where ``v_sphere`` is the
-volume of the sphere enclosing the outer shell.
-
-The bookkeeping required to create all the source models is handled by the
-:mod:`deproject` module.
-
-Fitting
------------
-Once the composite source models for each dataset are created the fit analysis
-can begin. Since the parameter space is typically large the usual procedure is
-to initally fit using the "onion-peeling" method:
-
- - First fit the outside shell model using the outer annulus spectrum
- - Freeze the model parameters for the outside shell
- - Fit the next inward annulus / shell and freeze those parameters
- - Repeat until all datasets have been fit and all shell parameters determined.
-
-From this point the user may choose to do a simultanenous fit of the shell
-models, possibly freezing some parameters as needed. This process is made
-manageable with the :mod:`specstack` methods that apply normal *Sherpa*
-commands like `freeze`_ or `set_par`_ to a stack of spectral datasets.
-
-Densities
-------------
-Physical densities (cm^-3) for each shell can be calculated with
-:mod:`deproject` assuming the source model is based on a thermal model with the
-"standard" normalization (from the `XSPEC`_ documentation):
-
-.. image:: thermal_norm.png
-
-Inverting this equation and assuming a constant ratio of N_H to electrons::
-
- n_e = sqrt(norm * 4*pi * DA^2 * 1e14 * (1+z)^2 / (volume * ne_nh_ratio))
-
- norm = model normalization from Sherpa fit
- DA = angular size distance (cm)
- volume = volume (cm^3)
- ne_nh_ratio = 1.18
-
-Recall that the model components for each volume element (intersection of the
-annular cylinder ``a`` with the spherical shell ``s``) are multiplied by a volume
-normalization::
-
- vol_norm[s,a] = v[s,a] / v_sphere
- v_sphere = volume of sphere enclosing outer annulus
-
-With this convention the ``volume`` used above in calculating the electron
-density for each shell is always ``v_sphere``.
-
-Download
-=========
-
-The :mod:`deproject` package is available for download at `deproject.tar.gz
-`_. The M87 data needed to run the example
-analysis is available as `m87.tar.gz `_
-
-The source is available on github at ``_.
-
-Installation
-=============
-
-The :mod:`deproject` package includes three Python modules that must be made
-available to the CIAO python so that *Sherpa* can import them. The first step
-is to untar the package tarball, change into the source directory, and initialize
-the CIAO environment::
-
- tar zxvf deproject-.tar.gz
- tar zxvf m87.tar.gz -C deproject-/examples # Needed for example / test script
- cd deproject-
- source /PATH/TO/ciao/bin/ciao.csh
-
-There are three methods for installing. Choose ONE of the three.
-
-**Simple:**
-
-The very simplest installation strategy is to just leave the module files in
-the source directory and set the ``PYTHONPATH`` environment variable to point
-to the source directory::
-
- setenv PYTHONPATH $PWD
-
-This method is fine in the short term but you always have to make sure
-``PYTHONPATH`` is set appropriately (perhaps in your ~/.cshrc file). And if you
-start doing much with Python you will have ``PYTHONPATH`` conflicts and things
-will get messy.
-
-**Better:**
-
-If you cannot write into the CIAO python library then do the following. These
-commands create a python library in your home directory and install the
-``deproject`` modules there. You could of course choose another directory
-instead of ``$HOME`` as the root of your python library.
-::
-
- mkdir -p $HOME/lib/python
- python setup.py install --home=$HOME
- setenv PYTHONPATH $HOME/lib/python
-
-Although you still have to set ``PYTHONPATH`` this method allows you to install
-other Python packages to the same library path. In this way you can make a
-local repository of packages that will run within Sherpa.
-
-**Best:**
-
-If you have write access to the CIAO installation you can just use the CIAO
-python to install the modules into the CIAO python library. Assuming you are
-in the CIAO environment do::
-
- python setup.py install
-
-This puts the new modules straight in to the CIAO python library so that any time
-you enter the CIAO environment they will be available. You do NOT need to set
-``PYTHONPATH``.
-
-Test
-=======
-
-To test the installation change to the source distribution directory and do the
-following::
-
- cd examples
- sherpa
- execfile('fit_m87.py')
- plot_fit(0)
- log_scale()
-
-This should run through in a reasonable time and produce output indicating the
-onion-peeling fit. The plot should show a good fit.
-
-Example: M87
-========================
-
-Now we step through in detail the ``fit_m87.py`` script in the ``examples``
-directory to explain each step and illustrate how to use the :mod:`deproject`
-module. This script should serve as the template for doing your own analysis.
-
-This example uses extracted spectra, response products, and analysis results
-for the Chandra observation of M87 (obsid 2707). These were kindly provided by Paul
-Nulsen. Results based on this observation can be found in `Forman et al 2005`_
-and via the CXC Archive `Obsid 2707 Publications`_ list.
-
-.. _`Forman et al 2005`: http://adsabs.harvard.edu/abs/2005ApJ...635..894F
-.. _`Obsid 2707 Publications`: http://cda.harvard.edu/chaser/viewerContents.do?obsid=2707&operation=ads
-
-The first step is to tell *Sherpa* about the Deproject class and
-set a couple of constants::
-
- from deproject import Deproject
-
- redshift = 0.004233 # M87 redshift
- arcsec_per_pixel = 0.492 # ACIS plate scale
- angdist = 4.9e25 # M87 distance (cm) (16 Mpc)
-
-Next we create a `numpy`_ array of the the annular radii in arcsec. The
-`numpy.arange`_ method here returns an array from 30 to 640 in steps of 30.
-These values were in pixels in the original spectral extraction so we convert
-to arcsec. (Note the convenient vector multiplication that is possible with
-`numpy`_.) ::
+The module can also be used with the `standalone Sherpa`_ release, but as
+it requires support for `XSPEC`_ models - for the thermal plasma emission
+codes - the documentation will focus on using Sherpa with a `CIAO`_
+environment.
- radii = numpy.arange(30., 640., 30) * arcsec_per_pixel
+.. _`standalone Sherpa`: https://sherpa.readthedocs.io/
-The ``radii`` parameter must be a list of values that starts with the inner
-radius of the inner annulus and includes each radius up through the outer
-radius of the outer annulus. Thus the ``radii`` list will be one element
-longer than the number of annuli.
-
-*Now the key step* of creating the :class:`Deproject` object ``dep``. This
-object is the interface to the all the :mod:`deproject`
-methods used for the deprojection analysis.
-::
-
- dep = Deproject(radii, theta=75, angdist=angdist)
-
-If you are not familiar with object oriented programming, the ``dep`` object is
-just a thingy that stores all the information about the deprojection analysis
-(e.g. the source redshift, PHA file information and the source model
-definitions) as object *attributes*. It also has object *methods*
-(i.e. functions) you can call such as ``dep.get_par(parname)`` or
-``dep.load_pha(file)``. The full list of attributes and methods are in the
-:mod:`deproject` module documentation.
-
-In this particular analysis the spectra were extracted from a 75 degree sector
-of the annuli, hence ``theta=75`` in the object initialization. For the
-default case of full 360 degree annuli this is not needed. Because the
-redshift is not a good distance estimator for M87 we also explicitly set the
-angular size distance.
-
-Now load the PHA spectral files for each annulus using the Python ``range``
-function to loop over a sequence ranging from 0 to the last annulus. The
-``load_pha()`` call is the first example of a :mod:`deproject` method
-(i.e. function) that mimics a *Sherpa* function with the same name. In this
-case ``dep.load_pha(file, annulus)`` loads the PHA file using the *Sherpa* `load_pha`_
-function but also registers the dataset in the spectral stack::
-
- for annulus in range(len(radii)-1):
- dep.load_pha('m87/r%dgrspec.pha' % (annulus+1), annulus)
-
-The ``annulus`` parameter is required in ``dep.load_pha()`` to support analysis
-of multi-obsid datasets.
-
-With the data loaded we set the source model for each of the spherical shells
-with the ``set_source()`` method. This is one of the more complex bits of
-:mod:`deproject`. It automatically generates all the model components for each
-shell and then assigns volume-weighted linear combinations of those components
-as the source model for each of the annulus spectral datasets::
-
- dep.set_source('xswabs * xsmekal')
-
-The model expression can be any valid *Sherpa* model expression with the following
-caveats:
-
- - Only the generic model type should be specified in the expression. In
- typical *Sherpa* usage one generates the model component name in the
- model expression, e.g. ``set_source("xswabs.abs1 * xsmekal.mek1")``. This
- would create model components named ``abs1`` and ``mek1``. In
- ``dep.set_source()`` the model component names are auto-generated as
- ``_``.
- - Only one of each model type can be used in the model expression. A source
- model expression like "xsmekal + gauss1d + gauss1d" would result in an error
- due to the model component auto-naming.
-
-Now the energy range used in the fitting is restricted using the stack version
-of the *Sherpa* `ignore`_ command. The `notice`_ command is also available.
-::
-
- dep.ignore(None, 0.5)
- dep.ignore(1.8, 2.2)
- dep.ignore(7, None)
-
-Next any required parameter values are set and their `freeze`_ or `thaw`_
-status are set.
-::
-
- dep.set_par('xswabs.nh', 0.0255)
- dep.freeze("xswabs.nh")
-
- dep.set_par('xsmekal.abundanc', 0.5)
- dep.thaw('xsmekal.abundanc')
-
- dep.set_par('xsmekal.redshift', redshift)
-
-As a convenience if any of the model components have a
-``redshift`` parameter that value will be used as the default redshift for
-calculating the angular size distance.
-
-At this point the model is completely set up and we are ready to do the initial
-"onion-peeling" fit. As for normal high-signal fitting with binned spectra we
-issue the commands to set the optimization method, set the fit statistic, and
-configure *Sherpa* to `subtract`_ the background when doing model fitting.
-Finally the :mod:`deproject` ``fit()`` method is called to perform the fit.
-::
-
- set_method("levmar") # Levenberg-Marquardt optimization method
- set_stat("chi2gehrels") # Gehrels Chi^2 fit statistic
- dep.subtract()
- dep.fit()
-
-After the fit process each shell model has an association normalization that
-can be used to calculate the densities. This is where the source angular
-diameter distance is used. If the angular diameter distance is not set
-explicitly in the original ``dep = Deproject(...)`` command then it is
-calculated automatically from the redshift found as a source model
-parameter. One can examine the values being used as follows::
-
- print "z=%.5f angdist=%.2e cm" % (dep.redshift, dep.angdist)
-
-The electron density is then calculated with the ``get_density()`` method and
-plotted in *Sherpa*::
-
- density_ne = dep.get_density()
- rad_arcmin = (dep.radii[:-1] + dep.radii[1:]) / 2.0 / 60.
- add_curve(rad_arcmin, density_ne)
- set_curve(['symbol.color', 'red', 'line.color', 'red'])
- set_plot_xlabel('Radial distance (arcmin)')
- set_plot_ylabel('Density (cm^{-3})')
- limits(X_AXIS, 0.2, 10)
- log_scale()
- print_window('m87_density', ['format', 'png'])
-
-The temperature profile from the :mod:`deproject` can be plotted as follows::
-
- kt = dep.get_par('xsmekal.kt') # returns array of kT values
- add_window()
- add_curve(rad_arcmin, kt)
- set_plot_xlabel('Radial distance (arcmin)')
- set_plot_ylabel('Density (cm^{-3})')
-
-The unphysical temperature oscillations seen here highlights a known issue
-with this analysis method.
-
-In the images below the :mod:`deproject` results (red) are compared with values
-(black) from an independent onion-peeling analysis by P. Nulsen using a custom
-perl script to generate `XSPEC`_ model definition and fit commands. These
-plots were created with the ``plot_m87.py`` script in the ``examples``
-directory. The agreement is good:
-
-.. image:: m87_density.png
-.. image:: m87_temperature.png
-
-
-Example: Multi-obsid
-==================================
-A second example illustrates the use of :mod:`deproject` for a multi-obsid
-observation of 3c186. It also shows how to set a background model for fitting
-with the ``cstat`` statistic. The extracted spectral data for this example are
-not yet publicly available.
-
-The script starts with some setup::
-
- import deproject
-
- radii = ('2.5', '6', '17')
- dep = deproject.Deproject(radii=[float(x) for x in radii])
-
- set_method("levmar")
- set_stat("cstat")
-
-Now we read in the data as before with ``dep.load_pha()``. The only difference
-here is an additional loop over the obsids. The ``dep.load_pha()`` function
-automatically extracts the obsid from the file header. This is used later in
-the case of setting a background model.
-::
-
- obsids = (9407, 9774, 9775, 9408)
- for ann in range(len(radii)-1):
- for obsid in obsids:
- dep.load_pha('3c186/%d/ellipse%s-%s.pi' % (obsid, radii[ann], radii[ann+1]), annulus=ann)
-
-Create and configure the source model expression as usual::
-
- dep.set_source('xsphabs*xsapec')
- dep.ignore(None, 0.5)
- dep.ignore(7, None)
- dep.freeze("xsphabs.nh")
-
- dep.set_par('xsapec.redshift', 1.06)
- dep.set_par('xsphabs.nh', 0.0564)
-
-Set the background model::
-
- execfile("acis-s-bkg.py")
- acis_s_bkg = get_bkg_source()
- dep.set_bkg_model(acis_s_bkg)
-
-Fit the projection model::
-
- dep.fit()
-
-To Do
-========
-
- - Use the Python logging module to produce output and allow for a verbosity
- setting. [Easy]
- - Create and use more generalized ``ModelStack`` and ``DataStack`` classes
- to allow for general mixing models. [Hard]
-
-Module docs
-====================
+.. toctree::
+ :maxdepth: 2
+ :caption: Introduction
+ overview
+ installation
+ changes
+ todo
+
.. toctree::
:maxdepth: 2
+ :caption: Examples
- deproject
- specstack
- cosmocalc
- references
+ examples/m87
+ examples/tie
+ examples/3c186
+
+.. toctree::
+ :maxdepth: 2
+ :caption: Module documentation
+
+ modules/deproject
+ modules/specstack
diff --git a/docs/installation.rst b/docs/installation.rst
new file mode 100644
index 0000000..c07cd32
--- /dev/null
+++ b/docs/installation.rst
@@ -0,0 +1,132 @@
+
+.. include:: references.rst
+
+Installation
+============
+
+As of :ref:`version 0.2.0 `,
+the :mod:`deproject` package can be installed
+directly from `PyPI`_. This requires
+`CIAO 4.11 `_
+or later (as installation with pip in earlier versions of
+CIAO is not well supported). The module *can* be used with the
+`standalone release of
+Sherpa `_,
+but it is only useful if Sherpa has been built with
+`XSPEC support `_
+which is trickier to achieve than we would like.
+
+Requirements
+------------
+
+The package uses `Astropy`_ and `SciPy`_, for units support and
+cosmological-distance calculations. It is assumed that
+`Matplotlib`_ is available for plotting (which is
+included in CIAO 4.11).
+
+Using pip
+---------
+
+CIAO
+^^^^
+
+Installation within CIAO requires:
+
+- having sourced the CIAO initialisation script (e.g.
+ `ciao.csh` or `ciao.bash`);
+
+- and using a constraints file, to
+ `avoid updating NumPy in CIAO `_.
+
+It should be as simple as::
+
+ echo "numpy==1.12.1" > constraints.txt
+ pip install -c constraints.txt 'astropy<3.1' deproject
+
+.. note::
+ The constraints are for CIAO 4.11, please adjust accordingly
+ if using a newer version of CIAO (the Astropy restriction is
+ because version 3.1 requires NumPy version 1.13 or later but
+ CIAO 4.11 is shipped with NumPy version 1.12).
+
+Standalone
+^^^^^^^^^^
+
+When using the standalone Sherpa intallation, the following should
+be all that is required::
+
+ pip install deproject
+
+Manual installation
+-------------------
+
+The source is available on github at
+``_, with releases available
+at ``_.
+
+After downloading the source code (whether from a release or
+by cloning the repository) and moving into the directory
+(`deproject-` or `deproject`), installation just
+requires::
+
+ python setup.py install
+
+.. note::
+ This command should only be run after setting up CIAO or whatever
+ Python environment contains your Standalone Sherpa installation.
+
+Test
+----
+
+The source installation includes a basic test suite, which can be
+run with either
+
+::
+
+ pytest
+
+or
+
+::
+
+ python setup.py test
+
+Example data
+------------
+
+The example data can be download from either
+http://cxc.cfa.harvard.edu/contrib/deproject/downloads/m87.tar.gz
+or from `GitHub `_.
+The source distribution includes scripts - in the
+`examples directory `_
+- that can
+be used to replicate both the :doc:`basic M87 example `
+and the follow-on example :doc:`combining annuli `.
+
+As an example (from within an IPython session, such as the
+Sherpa shell in CIAO)::
+
+ >>> %run fit_m87.py
+ ...
+ ... a lot of screen output will whizz by
+ ...
+
+The current density estimates can then be displayed with::
+
+ >>> dep.density_plot()
+
+.. image:: examples/m87_density.png
+
+the reduced-statistic for the fit to each shell with::
+
+ >>> dep.fit_plot('rstat')
+
+.. image:: examples/m87_rstat.png
+
+and the fit results for the first annulus can be displayed using
+the Sherpa functions::
+
+ >>> set_xlog()
+ >>> plot_fit_delchi(0)
+
+.. image:: examples/m87_ann0_fit.png
diff --git a/docs/m87_density.png b/docs/m87_density.png
deleted file mode 100644
index ae977b9..0000000
Binary files a/docs/m87_density.png and /dev/null differ
diff --git a/docs/m87_temperature.png b/docs/m87_temperature.png
deleted file mode 100644
index a15a7d5..0000000
Binary files a/docs/m87_temperature.png and /dev/null differ
diff --git a/docs/modules/deproject.rst b/docs/modules/deproject.rst
new file mode 100644
index 0000000..14e3f6b
--- /dev/null
+++ b/docs/modules/deproject.rst
@@ -0,0 +1,30 @@
+The deproject.deproject module
+==============================
+
+.. include:: ../references.rst
+
+.. currentmodule:: deproject.deproject
+
+.. automodule:: deproject.deproject
+
+ .. rubric:: Classes
+
+ .. autosummary::
+ :toctree: api
+
+ Deproject
+
+ .. rubric:: Functions
+
+ .. autosummary::
+ :toctree: api
+
+ deproject_from_xflt
+
+Class Inheritance Diagram
+-------------------------
+
+.. inheritance-diagram:: Deproject
+ :parts: 1
+
+
diff --git a/docs/modules/specstack.rst b/docs/modules/specstack.rst
new file mode 100644
index 0000000..40f93df
--- /dev/null
+++ b/docs/modules/specstack.rst
@@ -0,0 +1,16 @@
+The deproject.specstack module
+==============================
+
+.. include:: ../references.rst
+
+.. currentmodule:: deproject.specstack
+
+.. automodule:: deproject.specstack
+
+ .. rubric:: Classes
+
+ .. autosummary::
+ :toctree: api
+
+ SpecStack
+
diff --git a/docs/overview.rst b/docs/overview.rst
new file mode 100644
index 0000000..f6dae1f
--- /dev/null
+++ b/docs/overview.rst
@@ -0,0 +1,96 @@
+
+.. include:: references.rst
+
+Overview
+========
+
+The :mod:`deproject` module uses :mod:`specstack` to allow for manipulation of
+a stack of related input datasets and their models. Most of the functions
+resemble ordinary *Sherpa* commands (e.g. `set_par`_, `set_source`_, `ignore`_)
+but operate on a stack of spectra.
+
+The basic physical assumption of :mod:`deproject` is that the extended source
+emissivity is constant and optically thin within spherical shells whose radii
+correspond to the annuli used to extract the specta. Given this assumption one
+constructs a model for each annular spectrum that is a linear volume-weighted
+combination of shell models. The geometry is illustrated in the figure below
+(which would be rotated about the line to the observer in three-dimensions):
+
+.. image:: geometry.png
+
+Model creation
+--------------
+It is assumed that prior to starting :mod:`deproject` the user has extracted
+source and background spectra for each annulus. By convention the annulus
+numbering starts from the inner radius at 0 and corresponds to the dataset
+``id`` used within *Sherpa*. It is not required that the annuli include the
+center but they must be contiguous between the inner and outer radii.
+
+Given a spectral model ``M[s]`` for each shell ``s``, the source model for
+dataset ``a`` (i.e. annulus ``a``) is given by the sum over ``s >= a`` of
+``vol_norm[s,a] * M[s]`` (normalized volume * shell model). The image above
+shows shell 3 in blue and annulus 2 in red. The intersection of (purple) has a
+physical volume defined as ``vol_norm[3,2] * v_sphere`` where ``v_sphere``
+is the volume of the sphere enclosing the outer shell (as
+:ref:`discussed below `).
+
+The bookkeeping required to create all the source models is handled by the
+:mod:`deproject` module.
+
+Fitting
+-------
+Once the composite source models for each dataset are created the fit analysis
+can begin. Since the parameter space is typically large the usual procedure is
+to initally fit using the "onion-peeling" method:
+
+ - First fit the outside shell model using the outer annulus spectrum
+ - Freeze the model parameters for the outside shell
+ - Fit the next inward annulus / shell and freeze those parameters
+ - Repeat until all datasets have been fit and all shell parameters determined.
+
+From this point the user may choose to do a simultanenous fit of the shell
+models, possibly freezing some parameters as needed. This process is made
+manageable with the :mod:`specstack` methods that apply normal *Sherpa*
+commands like `freeze`_ or `set_par`_ to a stack of spectral datasets.
+
+.. _densities:
+
+Densities
+---------
+Physical densities (:math:`{\rm cm}^{-3}`) for each shell can be calculated with
+:mod:`deproject` assuming the source model is based on a thermal model with the
+"standard" normalization (from the `XSPEC`_ documentation):
+
+.. math::
+
+ {\rm norm} = \frac{10^{-14}}{4\pi [D_A (1+z)]^2} \int n_e n_H dV
+
+where :math:`D_A` is the angular size distance to the source (in cm),
+:math:`z` is the source redshift, and :math:`n_e` and :math:`n_H` are the
+electron and Hydrogen densities (in :math:`{\rm cm}^{-3}`).
+
+Inverting this equation and assuming constant values for :math:`n_e`
+and :math:`n_H` gives
+
+.. math::
+
+ n_e = \sqrt{\frac{4\pi \, {\rm norm} \, [D_A (1+z)]^2 \, 10^{14}}{\alpha V}}
+
+where :math:`V` is the volume, :math:`n_H = \alpha n_e`, and
+:math:`\alpha` is taken to be 1.18 (and this ratio is
+constant within the source).
+
+Recall that the model components for each volume element (intersection of the
+annular cylinder ``a`` with the spherical shell ``s``) are multiplied by a
+volume normalization:
+
+.. math::
+
+ V_{\rm norm}[s, a] = V[s, a] / V_{\rm sphere}
+
+where :math:`V_{\rm sphere}` is the volumne of the sphere enclosing the
+outermost annulus.
+
+With this convention the volume (:math:`V`) used above in calculating
+the electron density for each shell is always :math:`V_{\rm sphere}`.
+
diff --git a/docs/references.rst b/docs/references.rst
index a45caae..c93344a 100644
--- a/docs/references.rst
+++ b/docs/references.rst
@@ -1,15 +1,22 @@
.. _`CIAO`: http://cxc.harvard.edu/ciao/
.. _`sherpa`: http://cxc.harvard.edu/sherpa/
-.. _`XSPEC`: http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/
-.. _`projct`: https://astrophysics.gsfc.nasa.gov/XSPECwiki/projct_model
-.. _`plot_fit`: http://cxc.harvard.edu/sherpa/ahelp/plot_fit.py.html
-.. _`freeze`: http://cxc.harvard.edu/sherpa/ahelp/freeze.py.html
-.. _`thaw`: http://cxc.harvard.edu/sherpa/ahelp/thaw.py.html
-.. _`notice`: http://cxc.harvard.edu/sherpa/ahelp/notice.py.html
-.. _`subtract`: http://cxc.harvard.edu/sherpa/ahelp/subtract.py.html
-.. _`load_pha`: http://cxc.harvard.edu/sherpa/ahelp/load_pha.py.html
-.. _`set_par`: http://cxc.harvard.edu/sherpa/ahelp/set_par.py.html
-.. _`set_source`: http://cxc.harvard.edu/sherpa/ahelp/set_source.py.html
-.. _`ignore`: http://cxc.harvard.edu/sherpa/ahelp/ignore.py.html
-.. _`numpy`: http://www.scipy.org/NumPy
-.. _`numpy.arange`: http://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html#numpy.arange
+.. _`astropy`: http://docs.astropy.org/
+.. _`scipy`: https://www.scipy.org/scipylib/
+.. _`matplotlib`: https://matplotlib.org/
+.. _`XSPEC`: https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/
+.. _`projct`: https://asd.gsfc.nasa.gov/XSPECwiki/projct_model
+.. _`plot_fit`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.plot_fit.html
+.. _`freeze`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.freeze.html
+.. _`thaw`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.thaw.html
+.. _`notice`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.notice.html
+.. _`subtract`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.subtract.html
+.. _`load_pha`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.load_pha.html
+.. _`set_par`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.set_par.html
+.. _`set_source`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.set_source.html
+.. _`ignore`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.ignore.html
+.. _`link`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.link.html
+.. _`unlink`: https://sherpa.readthedocs.io/en/latest/ui/api/sherpa.astro.ui.unlink.html
+.. _`numpy`: https://www.scipy.org/NumPy
+.. _`numpy.arange`: https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html#numpy.arange
+.. _`PyPI`: https://pypi.org/
+.. _`Russell, Sanders, & Fabian 2008`: https://ui.adsabs.harvard.edu/?#abs/2008MNRAS.390.1207R
diff --git a/docs/rtd-pip-requirements b/docs/rtd-pip-requirements
new file mode 100644
index 0000000..6d43495
--- /dev/null
+++ b/docs/rtd-pip-requirements
@@ -0,0 +1,2 @@
+numpy
+sherpa
diff --git a/docs/specstack.rst b/docs/specstack.rst
deleted file mode 100644
index 3e968c1..0000000
--- a/docs/specstack.rst
+++ /dev/null
@@ -1,19 +0,0 @@
-:mod:`specstack`
-==================
-
-.. include:: references.rst
-
-.. automodule:: specstack
-
-Classes
---------
-
-.. autoclass:: SpecStack
- :show-inheritance:
- :members:
- :inherited-members:
- :undoc-members:
-
-
-
-
diff --git a/docs/thermal_norm.png b/docs/thermal_norm.png
deleted file mode 100644
index f6eb97b..0000000
Binary files a/docs/thermal_norm.png and /dev/null differ
diff --git a/docs/todo.rst b/docs/todo.rst
new file mode 100644
index 0000000..5b7ae29
--- /dev/null
+++ b/docs/todo.rst
@@ -0,0 +1,20 @@
+
+.. include:: references.rst
+
+To Do
+=====
+
+ - Use the Python logging module to produce output and allow for a verbosity
+ setting. [Easy]
+ - Move from using obsid to a more-generic label for the multi-data-set
+ case (i.e. when there are multiple data sets in a given annulus)
+ - Allow the theta value to vary between data set in each annulus (for the
+ multi-data-set case)
+ - Support elliptical annuli
+ - Create and use more generalized ``ModelStack`` and ``DataStack`` classes
+ to allow for general mixing models. [Hard]
+ - Add summary methods (similar to Sherpa's show series of commands)
+ - Add a specialized version of `plot_source_component` which takes a
+ general model expression (i.e. without the `_annulus` suffix) and
+ displays all combos per annulus. A bit tricky to get right.
+ There's also a `plot_model_component`.
diff --git a/examples/fit_3c186.py b/examples/fit_3c186.py
index 95b0a15..86cf70d 100644
--- a/examples/fit_3c186.py
+++ b/examples/fit_3c186.py
@@ -1,7 +1,8 @@
import deproject
+from astropy import units as u
radii = ('2.5', '6', '17')
-dep = deproject.Deproject(radii=[float(x) for x in radii])
+dep = deproject.Deproject(radii=[float(x) for x in radii] * u.arcsec)
set_method("neldermead")
set_method("levmar")
diff --git a/examples/fit_m87.py b/examples/fit_m87.py
index b179647..c07527d 100644
--- a/examples/fit_m87.py
+++ b/examples/fit_m87.py
@@ -1,21 +1,30 @@
-from deproject import Deproject
+"""
+Fit the M87 data and create basic plots.
-import logging
-logger = logging.getLogger("sherpa")
-logging.basicConfig(level=logging.INFO, filename="deproject.log", filemode="w")
+This assumes that Sherpa is using the Matplotlib backend; that is,
+the plot_pkg line in ~/.sherpa.rc is set to pylab and not chips.
+
+"""
import numpy
+from deproject.deproject import Deproject
+from sherpa.astro.ui import plot_fit, plot_fit_delchi, \
+ set_method, set_stat, set_xlog
+from astropy import units as u
+
+from matplotlib import pyplot as plt
redshift = 0.004233 # M87 redshift
-arcsec_per_pixel = 0.492 # ACIS plate scale
-angdist = 4.9e25 # M87 distance (cm) (16 Mpc)
+angdist = 16 * u.Mpc # M87 distance
+theta = 75 * u.deg # Covering angle of sectors
-radii = numpy.arange(30., 640., 30) * arcsec_per_pixel
-dep = Deproject(radii, theta=75, angdist=angdist)
+# Each pixel is 0.492 arcseconds
+radii = numpy.arange(30., 640., 30) * 0.492 * u.arcsec
+dep = Deproject(radii, theta=theta, angdist=angdist)
# Load datasets for each annulus
-for annulus in range(len(radii)-1):
- dep.load_pha('m87/r%dgrspec.pha' % (annulus+1), annulus)
+for annulus in range(len(radii) - 1):
+ dep.load_pha('m87/r%dgrspec.pha' % (annulus + 1), annulus)
# Subtract background
dep.subtract()
@@ -28,7 +37,7 @@
# Specify Galactic absorption
dep.set_par('xswabs.nh', 0.0255)
-dep.freeze("xswabs.nh")
+dep.freeze('xswabs.nh')
# Initialize abundance to 0.5 and thaw
dep.set_par('xsmekal.abundanc', 0.5)
@@ -37,9 +46,50 @@
# Set redshift
dep.set_par('xsmekal.redshift', redshift)
-# Do the initial onion-peeling fit
-set_method("levmar") # Levenberg-Marquardt optimization method
-set_stat("chi2gehrels") # Gehrels Chi^2 fit statistic
-dep.fit()
-#fit()
-conf_params = dep.conf()
+# Do the initial onion-peeling fit with
+# Levenberg-Marquardt optimization method
+# XSPEC variance with a Chi^2 fit statistic
+#
+set_method("levmar")
+set_stat("chi2xspecvar")
+
+dep.guess()
+
+# Display the central annulus (using Sherpa routines)
+set_xlog()
+plot_fit(0)
+plt.savefig('m87_ann0_guess.png')
+
+onion = dep.fit()
+
+# Display the central annulus (using Sherpa routines)
+plot_fit_delchi(0)
+plt.savefig('m87_ann0_fit.png')
+
+# Create basic plots
+dep.density_plot()
+plt.savefig('m87_density.png')
+
+dep.par_plot('xsmekal.kt')
+plt.savefig('m87_temperature.png')
+
+dep.fit_plot('rstat')
+plt.savefig('m87_rstat.png')
+
+# Error analysis
+errs = dep.covar()
+
+plt.subplot(2, 1, 1)
+dep.covar_plot('xsmekal.kt', clearwindow=False)
+plt.xlabel('')
+plt.subplot(2, 1, 2)
+dep.covar_plot('xsmekal.abundanc', clearwindow=False)
+plt.subplots_adjust(hspace=0.3)
+plt.savefig('m87_temperature_abundance.png')
+
+dep.covar_plot('density', units='arcmin')
+plt.savefig('m87_density_errs.png')
+
+# Display some basic information on the configuration
+#
+print("z={:.5f} angdist={}".format(dep.redshift, dep.angdist))
diff --git a/examples/fit_m87_da.py b/examples/fit_m87_da.py
deleted file mode 100644
index 766f8f1..0000000
--- a/examples/fit_m87_da.py
+++ /dev/null
@@ -1,32 +0,0 @@
-from deproject import Deproject
-from cosmocalc import cosmocalc
-
-redshift = 0.004233
-from math import pi
-
-dep = Deproject(numpy.arange(30., 640., 30)*0.492)
-dep.theta = 75
-dep.angdist = cosmocalc(redshift)['DA_cm'] * 0.892
-
-set_method("levmar")
-set_stat("chi2gehrels")
-
-for ann in range(dep.nshell):
- dep.load_pha('m87/r%dgrspec.pha' % (ann+1), annulus=ann)
-
-dep.set_source('xswabs*xsmekal')
-dep.ignore(None, 0.5)
-dep.ignore(1.8, 2.2)
-dep.ignore(7, None)
-
-dep.set_par('xswabs.nh', 0.0255)
-dep.freeze("xswabs.nh")
-
-dep.set_par('xsmekal.abundanc', 0.5)
-dep.thaw('xsmekal.abundanc')
-
-dep.set_par('xsmekal.redshift', redshift)
-dep.subtract()
-
-dep.fit()
-
diff --git a/examples/fit_m87_tied.py b/examples/fit_m87_tied.py
new file mode 100644
index 0000000..a835b67
--- /dev/null
+++ b/examples/fit_m87_tied.py
@@ -0,0 +1,31 @@
+#
+# This assumes fit_m87.py has been run and that no significant
+# changes have been made.
+#
+# It is run as
+# %run -i fit_m87_tied.py
+#
+
+print(dep.get_shells())
+
+dep.tie_par('xsmekal.kt', 17, 18)
+dep.tie_par('xsmekal.abundanc', 17, 18)
+
+print(dep.get_shells()[0:5])
+
+dep.fit()
+
+tied = dep.get_fit_results()
+print(len(tied))
+
+print(tied['xsmekal.kT'][16:])
+print(tied['xsmekal.Abundanc'][16:])
+
+dep.fit_plot('xsmekal.kt', units='pc')
+rmid = (onion['rlo_phys'] + onion['rhi_phys']) / 2
+plt.scatter(rmid.to(u.pc), onion['xsmekal.kT'], c='k')
+plt.savefig('m87_temperature_tied_comparison.png')
+
+dep.untie_par('xsmekal.kt', 18)
+dep.untie_par('xsmekal.abundanc', 18)
+print(dep.get_shells()[0:4])
diff --git a/examples/m87_density.png b/examples/m87_density.png
deleted file mode 100644
index ae977b9..0000000
Binary files a/examples/m87_density.png and /dev/null differ
diff --git a/examples/m87_temperature.png b/examples/m87_temperature.png
deleted file mode 100644
index a15a7d5..0000000
Binary files a/examples/m87_temperature.png and /dev/null differ
diff --git a/examples/plot_m87.py b/examples/plot_m87.py
index 8c12844..801aefe 100644
--- a/examples/plot_m87.py
+++ b/examples/plot_m87.py
@@ -1,43 +1,106 @@
"""
Plot density and temperature results for M87. Compare values from
deproject to those obtained by independent XSPEC analysis.
+
+The fit_m87.py script must have already been run, to set up the
+``dep`` object with the deprojection results.
+
+What are the errors in Paul's file (one sigma or 90%)?
+
"""
import numpy
+
+try:
+ from pycrates import read_file
+except ImportError:
+ try:
+ from astropy.table import Table
+ except ImportError:
+ raise ImportError("Unable to load pycrates or astropy")
+
+from matplotlib import pyplot as plt
+from matplotlib.collections import PatchCollection
+from matplotlib.patches import Rectangle
+import matplotlib.patches as mpatches
+
+
+# Draw a rectangle for each error box, following
+# https://matplotlib.org/2.2.3/gallery/statistics/errorbars_and_boxes.html#sphx-glr-gallery-statistics-errorbars-and-boxes-py
+#
+def make_error_boxes(ax, xdata, ydata, xerror, yerror, facecolor='gray',
+ edgecolor='None', alpha=0.5):
+
+ # Create list for all the error patches
+ errorboxes = []
+
+ # Loop over data points; create box from errors at each point
+ for x, y, xe, ye in zip(xdata, ydata, xerror.T, yerror.T):
+ rect = Rectangle((x - xe[0], y - ye[0]), xe.sum(), ye.sum())
+ errorboxes.append(rect)
+
+ # Create patch collection with specified colour/alpha
+ pc = PatchCollection(errorboxes, facecolor=facecolor, alpha=alpha,
+ edgecolor=edgecolor)
+
+ # Add collection to axes
+ ax.add_collection(pc)
+
+ """
+ # Plot errorbars
+ artists = ax.errorbar(xdata, ydata, xerr=xerror, yerr=yerror,
+ fmt='None', ecolor='k')
+
+ return artists
+ """
+
+ # Create a fake artist for the legend following
+ # https://matplotlib.org/users/legend_guide.html#creating-artists-specifically-for-adding-to-the-legend-aka-proxy-artists
+ #
+ return mpatches.Patch(facecolor=facecolor, edgecolor=edgecolor,
+ alpha=alpha)
+
+
+# Read in results from XSPEC analysis of M87
+
arcsec_per_pix = 0.492
arcsec_per_rad = (180.*3600.)/numpy.pi
-phys = read_file('m87_phys.dat') # values from XSPEC analysis
-
-r0 = get_colvals(phys, 'r0') # radius in ACIS pixels
-r1 = get_colvals(phys, 'r1')
-ne = get_colvals(phys, 'ne')
-nelow = get_colvals(phys, 'nelow')
-nehigh = get_colvals(phys, 'nehigh')
-kt = get_colvals(phys, 'kt')
-ktlow = get_colvals(phys, 'ktlow')
-kthigh = get_colvals(phys, 'kthigh')
-
-def add_set_clear_window(window_id):
- """Add a new window with given name ``window_id``, set it to be the current
- window, and clear it. The exception handling is to allow running this
- function repeatedly within an interactive session.
- """
- try:
- add_window(['id', window_id])
- except RuntimeError, e:
- set_current_window(window_id)
- try:
- clear_plot()
- except:
- pass
+# Read in the values from XSPEC analysis, where the radius is in
+# ACIS pixels.
+#
+try:
+ phys = read_file('m87_phys.dat')
+
+ r0 = phys.r0.values
+ r1 = phys.r1.values
+ ne = phys.ne.values
+ nelow = phys.nelow.values
+ nehigh = phys.nehigh.values
+ kt = phys.kt.values
+ ktlow = phys.ktlow.values
+ kthigh = phys.kthigh.values
+
+except NameError:
+
+ phys = Table.read('m87_phys.dat', format='ascii')
+
+ r0 = phys['r0']
+ r1 = phys['r1']
+ ne = phys['ne']
+ nelow = phys['nelow']
+ nehigh = phys['nehigh']
+ kt = phys['kt']
+ ktlow = phys['ktlow']
+ kthigh = phys['kthigh']
+
+##### Density ###########
-##### Temperature ###########
+x = (r0 + r1) / 2. * arcsec_per_pix / 60. # arcmin
+x_err = (r1 - r0) / 2. * arcsec_per_pix / 60. # arcmin
+dx = numpy.vstack((x_err, x_err))
-add_set_clear_window('temperature')
-x = (r0 + r1)/2. * arcsec_per_pix / 60. # arcmin
-x_err = (r1 - r0)/2. * arcsec_per_pix / 60. # arcmin
-y = kt
+y = ne
bad = nelow < 1e-4
y_err_low = kt - ktlow
@@ -45,55 +108,26 @@ def add_set_clear_window(window_id):
y_err_low[bad] = 0.01
y_err_high[bad] = 0.01
-add_curve(x, y, [y_err_low, y_err_high, x_err, x_err])
-set_curve(['symbol.color','green','line.color','green','err.color','green','err.thickness',3])
-log_scale(X_AXIS)
-
-d_kt = dep.get_par('xsmekal.kt')
-d_kt_m = []
-d_kt_p = []
-for n in range(len(radii)-1):
- d_kt_m.append(-conf_params[n].parmins[0])
- d_kt_p.append(conf_params[n].parmaxes[0])
-
-add_curve(x, d_kt)
-set_curve(['symbol.color', 'red', 'line.color', 'red'])
-set_plot_xlabel('Radial distance (arcmin)')
-set_plot_ylabel('Temperature (keV)')
-add_curve(x, d_kt, [d_kt_m, d_kt_p, x_err, x_err])
-set_curve(['symbol.color','red','line.color','red','err.color','red','err.thickness',2])
-limits(X_AXIS, 0.2, 6)
-# print_window('m87_temperature', ['format', 'png'])
+dy = numpy.vstack((y_err_low, y_err_high))
-##### Density ###########
+# Plot up the Sherpa results
+dep.covar_plot('density', ylog=True, units='arcmin')
+
+# Access a Matplotlib artist from the plot for the legend
+ax = plt.gca()
+sherpa_artist = ax.lines[1]
+
+# Plot up the XSPEC results
+#
+xspec_artist = make_error_boxes(ax, x, y, dx, dy)
+
+plt.legend((sherpa_artist, xspec_artist), ('Sherpa', 'XSPEC'))
+plt.savefig('m87_compare_density.png')
-r_sphere = dep.radii[-1] / arcsec_per_rad * angdist
-volume = 4 * numpy.pi / 3 * r_sphere**3 # volume of sphere enclosing outer shell (cm^3)
-mue=1.18
print "Calculating density for z=%.5f angdist=%.2e cm" % (dep.redshift, dep.angdist)
-norm = dep.get_par('xsmekal.norm') # returns array of kT values
-norm_m = []
-norm_p = []
-for n in range(len(radii)-1):
- norm_m.append(conf_params[n].parmins[2])
- norm_p.append(conf_params[n].parmaxes[2])
-
-dens = []
-dens_low = []
-dens_high = []
-dens_const = numpy.sqrt(4 * numpy.pi * angdist**2 * 1e14 * (1.0+redshift)**2 / volume * mue);
-for shell in range(dep.nshell):
- dens.append( dens_const*numpy.sqrt(norm[shell]))
- dens_low.append(dens_const*numpy.sqrt(norm[shell]+norm_m[shell]))
- dens_high.append(dens_const*numpy.sqrt(norm[shell]+norm_p[shell]))
-
-dens_loerr = numpy.array(dens)-numpy.array(dens_low)
-dens_hierr = numpy.array(dens_high)-numpy.array(dens)
-
-add_set_clear_window('density')
-y = ne
+y = kt
bad = nelow < 1e-4
y_err_low = ne - nelow
@@ -101,15 +135,16 @@ def add_set_clear_window(window_id):
y_err_low[bad] = 0.01
y_err_high[bad] = 0.01
-add_curve(x, y, [y_err_low, y_err_high, x_err, x_err])
-set_curve(['symbol.color','green','line.color','green','err.color','green','err.thickness',3])
-log_scale()
-
-add_curve(x, dens)
-set_curve(['symbol.color', 'red', 'line.color', 'red'])
-set_plot_xlabel('Radial distance (arcmin)')
-set_plot_ylabel('Density (cm^{-3})')
-add_curve(x, dens, [dens_loerr, dens_hierr, x_err, x_err])
-set_curve(['symbol.color','red','line.color','red','err.color','red','err.thickness',2])
-limits(X_AXIS, 0.2, 6)
-# print_window('m87_density', ['format', 'png'])
+dy = numpy.vstack((y_err_low, y_err_high))
+
+dep.covar_plot('xsmekal.kt', units='arcmin')
+
+ax = plt.gca()
+sherpa_axis = ax.lines[1]
+
+xspec_artist = make_error_boxes(ax, x, y, dx, dy)
+
+plt.legend((sherpa_artist, xspec_artist), ('Sherpa', 'XSPEC'),
+ loc='upper left')
+
+plt.savefig('m87_compare_temperature.png')
diff --git a/readthedocs.yml b/readthedocs.yml
new file mode 100644
index 0000000..8c2fe65
--- /dev/null
+++ b/readthedocs.yml
@@ -0,0 +1,13 @@
+
+# requirements_file: docs/rtd-pip-requirements
+
+# Use Python 3.5 since this is what is used by CIAO
+#
+python:
+ version: 3.5
+ setup_py_install: true
+ # use_system_site_packages: true
+
+conda:
+ file: rtd-environment.yml
+
\ No newline at end of file
diff --git a/rtd-environment.yml b/rtd-environment.yml
new file mode 100644
index 0000000..fb0a8a4
--- /dev/null
+++ b/rtd-environment.yml
@@ -0,0 +1,10 @@
+name: deproject
+
+channels:
+ - defaults
+ - sherpa
+
+dependencies:
+ - python=3.5
+ - numpy
+ - sherpa
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..4c35926
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,5 @@
+[aliases]
+test=pytest
+
+[build_sphinx]
+source-dir = docs
diff --git a/setup.py b/setup.py
index 2bfd0fc..36f7d22 100644
--- a/setup.py
+++ b/setup.py
@@ -1,20 +1,49 @@
-from setuptools import setup
+from setuptools import setup, find_packages
with open('README.rst', 'rt') as fh:
long_description = fh.read()
setup(name='deproject',
- version='0.1.3',
+ version='0.2.0',
license='BSD',
description='Sherpa deprojection package (X-ray analysis of Galaxy Clusters, Groups, and Galaxies)',
keywords='deprojection xray 3d 2d plasma Astrophysics onion',
long_description=long_description,
long_description_content_type='text/x-rst',
- author='Tom Aldcroft',
- author_email='aldcroft@cfa.harvard.edu',
- url='http://cxc.harvard.edu/contrib/deproject/',
- py_modules=['deproject', 'specstack', 'cosmocalc'],
- install_requires=['sherpa'],
+ author='Douglas Burke, Tom Aldcroft',
+ author_email='dburke.gw@gmail.com',
+
+ url='https://deproject.readthedocs.io/',
+ project_urls={
+ 'Documentation': 'https://deproject.readthedocs.io/',
+ 'Source Code': 'https://github.com/sherpa-deproject/deproject/',
+ 'Issues': 'https://github.com/sherpa-deproject/deproject/issues/',
+ },
+
+ packages=find_packages('src'),
+ package_dir={'': 'src'},
+
+ # NOTE:
+ # CIAO 4.11 comes with NumPy 1.12.1, but AstroPy 3.1 requires
+ # NumPy >= 1.13.1, so add the restriction here.
+ #
+ # This is restrictive for standalone Sherpa users, but it
+ # requires a *lot* of effort to build the needed XSPEC model
+ # support into Sherpa, so assume that the majority of use
+ # cases will be installing into CIAO.
+ #
+ # SciPy may be needed. If users specify the angular-diameter
+ # distance to the source, or use a cosmology that does not
+ # use SciPy interpolation then SciPy is not needed. However,
+ # this may only be known about after doing a fit, which is
+ # annoying, so force SciPy on all users.
+ #
+ setup_requires=['pytest-runner'],
+ tests_require=['pytest'],
+ install_requires=['sherpa', 'astropy<3.1', 'scipy'],
+
+ python_requires='~=3.5',
+
classifiers=[
'Development Status :: 4 - Beta',
'Intended Audience :: Science/Research',
diff --git a/specstack.py b/specstack.py
deleted file mode 100644
index 22d8a5a..0000000
--- a/specstack.py
+++ /dev/null
@@ -1,270 +0,0 @@
-"""
-Manipulate a stack of spectra in Sherpa.
-
-The methods in the SpecStack class provide a way to automatically apply
-familiar Sherpa commands such as `set_par`_ or `freeze`_ or `plot_fit`_
-to a stack of PHA spectra. This simplifies simultaneous fitting of
-multiple spectra.
-
-Note that the :mod:`specstack` module is currently distributed within with the
-:mod:`deproject` package. `Specstack` is not yet fully documented or tested
-outside the context of `deproject`.
-
-:Copyright: Smithsonian Astrophysical Observatory (2009)
-:Author: Tom Aldcroft (aldcroft@head.cfa.harvard.edu)
-"""
-import re
-import numpy
-import sherpa.astro.ui as SherpaUI
-try:
- import pycrates
- import pychips
-except ImportError:
- # Allow doc generation to work without pychips and pycrates
- print "ERROR: could not import pycrates and pychips"
-
-def _sherpa_plot_func(func):
- def _sherpa_plot_func_(self, *args, **kwargs):
- self._sherpa_plot(func, *args, **kwargs)
- return _sherpa_plot_func_
-
-class SpecStack(object):
- """
- Manipulate a stack of spectra in Sherpa.
- """
- def __init__(self):
- self.datasets = []
- self.model_comps = [] # Model components
- self.obsids = set()
- self.srcmodel_comps = [] # Generic model components in source model expression
- self.model_comps = [] # All instantiated model components for shells
-
- def load_pha(self, specfile, annulus):
- """
- Load a pha file and add to the datasets for stacked analysis.
-
- :param specfile: extracted source PHA/PI spectrum file
- :param annulus: annulus for spectrum file
- """
- dataid = len(self.datasets)
- print 'Loading spectrum file %s as dataset id %d' % (specfile, dataid)
- SherpaUI.load_pha(dataid, specfile)
-
- try:
- obsid = int(pycrates.read_file(specfile).get_key_value('OBS_ID'))
- except (TypeError, ValueError):
- obsid = 0
- dataset = dict(file=specfile,
- obsid=obsid,
- id=dataid,
- annulus=annulus
- )
- self.datasets.append(dataset)
- self.obsids.add(obsid)
-
- def find_parval(self, parname):
- """
- Find the value of the first parameter with the given ``parname``. Ignore
- case when matching.
-
- :param parname: parameter name
- :rtype: parameter value
- """
- RE_parname = re.compile(parname + '$', re.IGNORECASE)
- for model_comp in self.model_comps:
- mc = model_comp['object'] # Get the model component object
- for par in mc.pars:
- if RE_parname.match(par.name):
- return par.val
-
- raise ValueError('Parameter %s not found in any model component' % shell)
-
- def find_norm(self, shell):
- """
- Find the normalization value for the ``shell`` model. Check for multiple
- or missing norms.
-
- :param shell: shell index
- :rtype: norm value
- """
- norms = []
- for model_comp in self.model_comps:
- if model_comp['shell'] == shell:
- mc = model_comp['object'] # Get the model component object
- if hasattr(mc, 'norm'): # Special attr of model component
- norms.append(mc.norm.val)
-
- if len(norms) == 0:
- raise ValueError('Model for shell %d has no norm' % shell)
- elif len(norms) > 1:
- raise ValueError('Model for shell %d has multiple norms' % shell)
-
- return norms[0]
-
- def _get_n_datasets(self):
- """
- Return the number of datasets
- :rtype: int"""
- return len(self.datasets)
-
- n_datasets = property(_get_n_datasets)
-
- def _sherpa_cmd(self, func, *args):
- """
- Apply an arbitrary Sherpa function to each of the datasets.
- :rtype: None
- """
- for dataset in self.datasets:
- func(dataset['id'], *args)
-
- def subtract(self, *args):
- """Subtract background from each dataset"""
- self._sherpa_cmd(SherpaUI.subtract, *args)
-
- def notice(self, *args):
- """Apply Sherpa notice command to each dataset."""
- self._sherpa_cmd(SherpaUI.notice_id, *args)
-
- def ignore(self, *args):
- """Apply Sherpa ignore command to each dataset."""
- self._sherpa_cmd(SherpaUI.ignore_id, *args)
-
- def _sherpa_par(self, func, par, msg, *args):
- """Apply ``func(*args)`` to all shell model component parameters named ``par``.
-
- See thaw(), freeze(), set_par() and get_par() for examples.
-
- :param func: Sherpa function that takes a full parameter name specification and
- optional args, e.g. set_par() used as set_par('xsmekal_7.kt', 2.0)
- :param par: Model type and param name as in source model expression e.g. 'xsmekal.kt'
- :param msg: Format string to indicate action.
- :param *args: Optional function arguments
-
- :rtype: numpy array of function return values ordered by shell
- """
- model_type, parname = par.split('.')
-
- vals = [] # return values
- for model_comp in self.model_comps:
- if model_comp['type'] == model_type:
- fullparname = '%s.%s' % (model_comp['name'], parname)
- if msg is not None:
- print msg % fullparname
- vals.append(func(fullparname, *args))
-
- return vals
-
- def thaw(self, par):
- """Apply thaw command to specified parameter for each dataset.
-
- :param par: parameter specifier in format .
- :rtype: None
- """
- self._sherpa_par(SherpaUI.thaw, par, 'Thawing %s')
-
- def freeze(self, par):
- """Apply freeze command to specified parameter for each dataset.
-
- :param par: parameter specifier in format .
- :rtype: None
- """
- self._sherpa_par(SherpaUI.freeze, par, 'Freezing %s')
-
- def set_par(self, par, val):
- """Set parameter value for each dataset.
-
- :param par: parameter specifier in format .
- :param val: parameter value
- :rtype: None
- """
- self._sherpa_par(SherpaUI.set_par, par, 'Setting %%s=%s' % str(val), val)
-
- def get_par(self, par):
- """Get array of parameter values for datasets.
-
- :param par: parameter specifier in format .
- :param val: parameter value
- :rtype: numpy array of parameter value ordered by dataset
- """
- pars = self._sherpa_par(SherpaUI.get_par, par, 'Getting %s')
- return numpy.array([x.val for x in pars])
-
- def _sherpa_plot(self, func, *args, **kwargs):
- """Call Sherpa plot ``func`` for each dataset.
-
- :param func: Sherpa plot function
- :param args: plot function list arguments
- :param kwargs: plot function named (keyword) arguments
- :rtype: None
- """
- for shell in range(self.nshell):
- window_id = 'Shell%d' % shell
- try:
- pychips.add_window(['id', window_id])
- except RuntimeError:
- pass # already exists
-
- new_args = args
- if len(args) > 0:
- # Try to format first arg
- try:
- new_args = tuple([args[0] % shell]) + args[1:]
- except TypeError:
- pass
-
- pychips.set_current_window(window_id)
- func(*new_args, **kwargs)
-
- def print_window(self, *args, **kwargs):
- """Print window for each dataset.
-
- :param args: list arguments to pass to print_window
- :param kwargs: named (keyword) arguments to pass to print_window
- :rtype: None
- """
- if len(args) > 0:
- args = tuple([args[0] + '%d']) + args[1:]
- self._sherpa_plot(pychips.plot_window, *args, **kwargs)
-
- def dummyfunc(self, *args, **kwargs):
- pass
-
- try:
- log_scale = _sherpa_plot_func(pychips.log_scale)
- linear_scale = _sherpa_plot_func(pychips.linear_scale)
- except NameError:
- # Allow doc generation to work without pychips
- log_scale = dummyfunc
- linear_scale = dummyfunc
-
- plot_fit = _sherpa_plot_func(SherpaUI.plot_fit)
- plot_arf = _sherpa_plot_func(SherpaUI.plot_arf)
- plot_bkg_fit = _sherpa_plot_func(SherpaUI.plot_bkg_fit)
- plot_bkg_ratio = _sherpa_plot_func(SherpaUI.plot_bkg_ratio)
- plot_chisqr = _sherpa_plot_func(SherpaUI.plot_chisqr)
- plot_fit_delchi = _sherpa_plot_func(SherpaUI.plot_fit_delchi)
- plot_psf = _sherpa_plot_func(SherpaUI.plot_psf)
- plot_bkg = _sherpa_plot_func(SherpaUI.plot_bkg)
- plot_bkg_fit_delchi = _sherpa_plot_func(SherpaUI.plot_bkg_fit_delchi)
- plot_bkg_resid = _sherpa_plot_func(SherpaUI.plot_bkg_resid)
- plot_data = _sherpa_plot_func(SherpaUI.plot_data)
- plot_fit_resid = _sherpa_plot_func(SherpaUI.plot_fit_resid)
- plot_ratio = _sherpa_plot_func(SherpaUI.plot_ratio)
- plot_bkg_chisqr = _sherpa_plot_func(SherpaUI.plot_bkg_chisqr)
- plot_bkg_fit_resid = _sherpa_plot_func(SherpaUI.plot_bkg_fit_resid)
- plot_bkg_source = _sherpa_plot_func(SherpaUI.plot_bkg_source)
- plot_delchi = _sherpa_plot_func(SherpaUI.plot_delchi)
- plot_model = _sherpa_plot_func(SherpaUI.plot_model)
- plot_resid = _sherpa_plot_func(SherpaUI.plot_resid)
- plot_bkg_delchi = _sherpa_plot_func(SherpaUI.plot_bkg_delchi)
- plot_bkg_model = _sherpa_plot_func(SherpaUI.plot_bkg_model)
- # CIAO 4.3 re-named plot_bkg_unconvolved to plot_bkg_source
- if hasattr(SherpaUI, "plot_bkg_unconvolved"):
- plot_bkg_unconvolved = _sherpa_plot_func(SherpaUI.plot_bkg_unconvolved)
- else:
- plot_bkg_unconvolved = _sherpa_plot_func(SherpaUI.plot_bkg_source)
- plot_bkg_source = _sherpa_plot_func(SherpaUI.plot_bkg_source)
- plot_fit = _sherpa_plot_func(SherpaUI.plot_fit)
- plot_order = _sherpa_plot_func(SherpaUI.plot_order)
- plot_source = _sherpa_plot_func(SherpaUI.plot_source)
-
diff --git a/src/deproject/__init__.py b/src/deproject/__init__.py
new file mode 100644
index 0000000..388b9e9
--- /dev/null
+++ b/src/deproject/__init__.py
@@ -0,0 +1,12 @@
+# To support users of version 0.1, re-export the deproject module
+
+"""Re-export the deproject.deproject module.
+
+The module currently provides:
+
+ Deproject
+ deproject_from_xflt
+
+"""
+
+from . deproject import *
diff --git a/src/deproject/deproject.py b/src/deproject/deproject.py
new file mode 100644
index 0000000..abed2b7
--- /dev/null
+++ b/src/deproject/deproject.py
@@ -0,0 +1,2083 @@
+"""
+Deproject 2-d circular annular spectra to 3-d object properties.
+
+This module implements the "onion-skin" approach popular in X-ray
+analysis of galaxy clusters and groups to estimate the three-dimensional
+temperature, metallicity, and density distributions of an optically-thin
+plasma from the observed (projected) two-dimensional data, arranged in
+concentric circular annuli.
+
+:Copyright: Smithsonian Astrophysical Observatory (2009, 2019)
+:Author: Tom Aldcroft (taldcroft@cfa.harvard.edu), Douglas Burke (dburke@cfa.harvard.edu)
+"""
+
+from collections import defaultdict, OrderedDict
+import copy
+import logging
+from math import pi
+import re
+
+import numpy
+from astropy.table import Table
+from astropy import units as u
+from astropy.cosmology import Planck15
+
+from sherpa.plot import plotter
+from sherpa.astro import ui
+from sherpa.astro.io import read_pha
+from sherpa.models.parameter import CompositeParameter
+
+from . import specstack
+from . import simplegraph
+from . import fieldstore
+
+__all__ = ("Deproject", "deproject_from_xflt")
+
+
+_sherpa_logger = logging.getLogger('sherpa')
+
+
+arcsec_per_rad = (u.radian / u.arcsec).to(1)
+
+
+class Deproject(specstack.SpecStack):
+ """Support deprojecting a set of spectra (2-d concentric circular annuli).
+
+ Parameters
+ ----------
+ radii : AstroPy Quantity representing an angle on the sky
+ The edges of each annulus, which must be circular, concentric,
+ in ascending order, and >= 0. If there are n annuli then there are n+1
+ radii, since the start and end of the sequence must be given.
+ The units are expected to be arcsec, arcminute, or degree.
+ theta : AstroPy Quantity (scalar or array) representing an angle
+ The "fill factor" of each annulus, given by the azimuthal coverage
+ of the shell in degrees. The value can be a scalar, so the same
+ value is used for all annuli, or a sequence with a length
+ matching the number of annuli. Since the annulus assumes circular
+ symmetry there is no need to define the starting point of the
+ measurement, for cases when the value is less than 360 degrees.
+ angdist : None or AstroPy.Quantity, optional
+ The angular-diameter distance to the source. If not given then
+ it is calculated using the source redshift along with the
+ `cosmology` attribute.
+ cosmology : None or astropy.cosmology object, optional
+ The cosmology used to convert redshift to an angular-diameter
+ distance. This is used when `angdist` is None. If `cosmology`
+ is None then the `astropy.cosmology.Planck15` Cosmology
+ object is used.
+
+ Examples
+ --------
+
+ The following highly-simplified example fits a deprojected model
+ to data from three annuli - ann1.pi, ann2.pi, and ann3.pi - and
+ also calculates errors on the parameters using the confidence
+ method::
+
+ >>> dep = Deproject([0, 10, 40, 100] * u.arcsec)
+ >>> dep.load_pha('ann1.pi', 0)
+ >>> dep.load_pha('ann2.pi', 1)
+ >>> dep.load_pha('ann3.pi', 2)
+ >>> dep.subtract()
+ >>> dep.notice(0.5, 7.0)
+ >>> dep.set_source('xsphabs * xsapec')
+ >>> dep.set_par('xsapec.redshift', 0.23)
+ >>> dep.thaw('xsapec.abundanc')
+ >>> dep.set_par('xsphabs.nh', 0.087)
+ >>> dep.freeze('xsphabs.nh')
+ >>> dep.fit()
+ >>> dep.fit_plot('rstat')
+ >>> errs = dep.conf()
+ >>> dep.conf_plot('density')
+
+ """
+
+ @u.quantity_input(radii='angle', theta='angle', angdist='length')
+ def __init__(self, radii,
+ theta=360 * u.deg,
+ angdist=None,
+ cosmology=None):
+
+ nshell = numpy.size(radii) - 1
+ if nshell < 1:
+ raise ValueError('radii parameter must be a sequence ' +
+ 'with at least two values')
+
+ dr = radii[1:] - radii[:-1]
+ if numpy.any(dr <= 0):
+ raise ValueError('radii parameter must be in increasing order')
+
+ # All values must be >= 0
+ #
+ if radii[0] < 0:
+ raise ValueError('radii must be >= 0')
+
+ self.radii = radii
+ self.nshell = nshell
+
+ ntheta = numpy.size(theta)
+ if ntheta == 1:
+ thetas = numpy.repeat(theta, nshell)
+ elif ntheta == nshell:
+ thetas = theta
+ else:
+ raise ValueError('theta must be a scalar or ' +
+ 'match the number of annuli')
+
+ theta_min = thetas.min()
+ if theta_min <= (0.0 * u.deg):
+ raise ValueError('theta must be > 0 degrees')
+
+ theta_max = thetas.max()
+ if theta_max > (360.0 * u.deg):
+ raise ValueError('theta must be <= 360 degrees')
+
+ self._theta = thetas
+
+ if angdist is not None:
+ self._set_angdist(angdist)
+ else:
+ self._angdist = None
+
+ self._redshift = None
+ self._fit_results = None
+ self._covar_results = None
+ self._conf_results = None
+
+ self._cosmology = Planck15 if cosmology is None else cosmology
+
+ super().__init__()
+
+ def load_pha(self, specfile, annulus):
+ if annulus < 0 or annulus >= self.nshell:
+ raise ValueError("Expected 0 <= annulus < " +
+ "{} but sent {}".format(self.nshell, annulus))
+
+ super().load_pha(specfile, annulus)
+
+ def _get_redshift(self):
+ if self._redshift is None:
+ self._redshift = self.find_parval('redshift')
+ return self._redshift
+
+ def _set_redshift(self, redshift):
+ self._redshift = redshift
+
+ # Perhaps the angular-diameter distance shouldn't be cached if
+ # not explicitly set. This lets the value be updated if the
+ # redshift or cosmology object is updated. Alternatively, we
+ # tell users they have to manually set da if these things
+ # change.
+ #
+ def _get_angdist(self):
+ if self._angdist is None:
+ da = self.cosmology.angular_diameter_distance(self.redshift)
+ self._angdist = da
+ return self._angdist
+
+ @u.quantity_input(angdist='length')
+ def _set_angdist(self, angdist):
+
+ if angdist <= 0:
+ raise ValueError("angdist must be > 0")
+
+ self._angdist = angdist
+
+ redshift = property(_get_redshift, _set_redshift, None,
+ "Source redshift")
+ angdist = property(_get_angdist, _set_angdist, None,
+ "Angular size distance (an AstroPy quantity)")
+
+ @property
+ def cosmology(self):
+ """Return the cosmology object (only used if angdist not set)"""
+ return self._cosmology
+
+ def _calc_vol_norm(self):
+ r"""Calculate the normalized volumes of the deprojected views.
+
+ Sets the `vol_norm` field to a matrix of the normalized volumes
+ of the cylindrical annuli intersecting with the spherical shell.
+ The matrix is defined as $volume[i, j] / V_sphere$, where
+ $i$ represents the shell and $j$ the annulus (with indexes
+ starting at 0), $V_sphere = 4/3 \pi r_o^3$, and $r_o$ is the
+ outermost radius of the shells.
+ """
+
+ # The units for the radii are not important here
+ r = self.radii.value
+ theta_rad = self._theta.to_value(u.rad)
+
+ cv = numpy.zeros([self.nshell, self.nshell])
+ v = numpy.zeros([self.nshell, self.nshell])
+
+ for a, ra0 in enumerate(r[:-1]):
+ # Annulus
+ ra1 = r[a + 1]
+ ra0sq = ra0**2
+ ra1sq = ra1**2
+
+ for s, rs0 in enumerate(r[:-1]):
+ # Spherical shell
+ rs1 = r[s + 1]
+ if s >= a:
+ # Volume of cylindrical annulus (ra0,ra1) intersecting
+ # the sphere (rs1)
+ #
+ rs1sq = rs1**2
+ rterm = (rs1sq - ra0sq)**1.5 - (rs1sq - ra1sq)**1.5
+ cv[s, a] = 2.0 * theta_rad[a] / 3 * rterm
+
+ # Volume of annulus (ra0,ra1) intersecting the spherical
+ # shell (rs0,rs1)
+ v[s, a] = cv[s, a]
+ if s - a > 0:
+ v[s, a] -= cv[s - 1, a]
+
+ self.vol_norm = v / (4. * pi / 3. * r[-1]**3)
+
+ def _create_name(self, model_name, annulus):
+ """Create the name used for a model component for the given annulus.
+
+ Parameters
+ ----------
+ model_name : str
+ The Sherpa model name (e.g. 'xsphabs').
+ annulus : int
+ The annulus number.
+
+ Returns
+ -------
+ name : str
+ The name of the model component.
+
+ Notes
+ -----
+ At present there is no real need to allow the naming scheme to
+ be changed (e.g. in a sub-class), but it is useful to help
+ record where names are created.
+ """
+
+ # This is perhaps a bit "over engineered"
+ #
+ name = '{}_{}'.format(model_name, annulus)
+ return name
+
+ def _create_src_model_components(self):
+ """Create the model components for each shell."""
+
+ self._reset_model_comps()
+
+ # Find the generic components in source model expression
+ # and set up their names.
+ #
+ RE_model = re.compile(r'\b \w+ \b', re.VERBOSE)
+ for match in RE_model.finditer(self.srcmodel):
+ model_type = match.group()
+
+ store = dict(type=model_type,
+ start=match.start(),
+ end=match.end())
+ self.srcmodel_comps.append(store)
+
+ # For each shell create the corresponding model components so they can
+ # be used later to create composite source models for each dataset
+ for shell in range(self.nshell):
+ for srcmodel_comp in self.srcmodel_comps:
+ model_comp = {}
+ model_comp['type'] = srcmodel_comp['type']
+ name = self._create_name(model_comp['type'], shell)
+ model_comp['name'] = name
+ model_comp['shell'] = shell
+ comp = ui.create_model_component(model_comp['type'], name)
+ model_comp['object'] = comp
+ self.model_comps.append(model_comp)
+
+ def set_source(self, srcmodel='xsphabs*xsapec'):
+ """Create a source model for each annulus.
+
+ Unlike the standard `set_source` command, this version just
+ uses the , not ., since
+ the is automatically created for users by appending
+ the annulus number to .
+
+ Parameters
+ ----------
+ srcmodel : str, optional
+ The source model expression applied to each annulus.
+
+ See Also
+ --------
+ set_bkg_model, set_par
+
+ Notes
+ -----
+ The data must have been read in for all the data before calling
+ this method (this matches Sherpa, where you can not call set_source
+ unless you have already loaded the data to fit).
+
+ Examples
+ --------
+
+ The following two calls have the same result: model instances
+ called 'xsphabs' and 'xsapec' are created
+ for each annulus, and the source expression for the annulus
+ set to their multiplication:
+
+ >>> dep.set_source()
+ >>> dep.set_source('xsphabs * xsapec')
+
+ Use the XSPEC vapec model rather than the apec model to
+ represent the plasma emission:
+
+ >>> dep.set_source('xsphabs * xsvapec')
+
+ """
+
+ # We can not check that all data has been loaded in (that is,
+ # if there are multiple data sets per annulis), but we can at
+ # least ensure that there is a dataset loaded
+ # for each annulus.
+ #
+ seen = set([])
+ for dataset in self.datasets:
+ seen.add(dataset['annulus'])
+
+ expected = set(range(self.nshell))
+ diff = sorted(list(expected.difference(seen)))
+ if len(diff) == 1:
+ raise ValueError("missing data for annulus {}".format(diff[0]))
+ elif len(diff) > 0:
+ raise ValueError("missing data for annuli {}".format(diff))
+
+ self.srcmodel = srcmodel
+ self._calc_vol_norm()
+ self._create_src_model_components()
+
+ # TODO: isn't it better to loop over annuli, out to in, to
+ # avoid some repeated work?
+ #
+ for dataset in self.datasets:
+ dataid = dataset['id']
+ annulus = dataset['annulus']
+ modelexprs = []
+ for shell in range(annulus, self.nshell):
+ srcmodel = self.srcmodel
+ for model_comp in reversed(self.srcmodel_comps):
+ i0 = model_comp['start']
+ i1 = model_comp['end']
+ model_comp_name = self._create_name(model_comp['type'],
+ shell)
+ srcmodel = srcmodel[:i0] + model_comp_name + srcmodel[i1:]
+
+ f = self.vol_norm[shell, annulus]
+ modelexprs.append('{} * {}'.format(f, srcmodel))
+
+ modelexpr = " + ".join(modelexprs)
+ print('Setting source model for dataset %d = %s' % (dataid, modelexpr))
+ ui.set_source(dataid, modelexpr)
+
+ def set_bkg_model(self, bkgmodel):
+ """Create a background model for each annulus.
+
+ The background model is the same between the annuli, except that
+ a scaling factor is added for each annulus (to allow for
+ normalization uncertainities). The scaling factors are labelled
+ 'bkg_norm_', and at least one of these must be frozen
+ (otherwise it is likely to be degenerate with the background
+ normalization, causing difficulties for the optimiser).
+
+ Parameters
+ ----------
+ bkgmodel : model instance
+ The background model expression applied to each annulus.
+ Unlike set_source this should be the actual model instance,
+ and not a string.
+
+ See Also
+ --------
+ set_source, set_par
+
+ Examples
+ --------
+
+ Model the background with a single power-law component:
+
+ >>> dep.set_bkg_model(xspowerlaw.bpl)
+
+ """
+
+ self.bkgmodel = bkgmodel
+
+ # TODO:
+ # - record the background components in the same way the source
+ # is done
+ # - should the background be allowed to have different components
+ # per annulus?
+ #
+ bkg_norm = {}
+ for obsid in self.obsids:
+ bkg_norm_name = 'bkg_norm_%d' % obsid
+ print('Creating model component xsconstant.%s' % bkg_norm_name)
+ bcomp = ui.create_model_component('xsconstant', bkg_norm_name)
+ bkg_norm[obsid] = bcomp
+
+ for dataset in self.datasets:
+ print('Setting bkg model for dataset %d to bkg_norm_%d' % (dataset['id'], dataset['obsid']))
+ ui.set_bkg_model(dataset['id'],
+ bkg_norm[dataset['obsid']] * bkgmodel)
+
+ def get_shells(self):
+ """How are the annuli grouped?
+
+ An annulus may have multiple data sets associated with it, but
+ it may also be linked to other annuli due to tied parameters.
+ The return value is per group, in the ordering needed for
+ the outside-to-inside onion skin fit, where the keys for
+ the dictionary are 'annuli' and 'dataids'.
+
+ Returns
+ -------
+ groups : list of dicts
+ Each dictionary has the keys 'annuli' and 'dataids', and
+ lists the annuli and data identifiers that are fit together.
+ The ordering matches that of the onion-skin approach, so
+ the outermost group first.
+
+ See Also
+ --------
+ get_radii, tie_par
+
+ Examples
+ --------
+
+ For a 3-annulus deprojection where there are no parameter ties
+ to combine annului:
+
+ >>> dep.get_shells()
+ [{'annuli': [2], 'dataids': [2]},
+ {'annuli': [1], 'dataids': [1]},
+ {'annuli': [0], 'dataids': [0]}]
+
+ After tie-ing the abundance parameter for the outer two shells,
+ there are now two groups of annuli:
+
+ >>> dep.tie_par('xsapec.abundanc', 1, 2)
+ Tying xsapec_2.Abundanc to xsapec_1.Abundanc
+ >>> dep.get_shells()
+ [{'annuli': [1, 2], 'dataids': [1, 2]},
+ {'annuli': [0], 'dataids': [0]}]
+
+ """
+
+ # Map from model component name (e.g. 'xsapec_2') to shell
+ # number.
+ #
+ cpt_map = {}
+ for model_comp in self.model_comps:
+ mname = model_comp['name']
+ assert mname not in cpt_map
+ cpt_map[mname] = model_comp['shell']
+
+ # Find the connected shells/annuli.
+ #
+ graph = simplegraph.SimpleGraph()
+ for shell in range(self.nshell):
+ graph.add_link(shell, shell)
+
+ for model_comp in self.model_comps:
+ shell = model_comp['shell']
+ for par in model_comp['object'].pars:
+ if par.link is None:
+ continue
+
+ # For the moment only support tied parameters (i.e.
+ # they are set equal). It should be possible to just
+ # iterate through the parts of the composite parameter
+ # and extract the links (since there could be more than
+ # one), but do not try this yet.
+ #
+ if isinstance(par.link, CompositeParameter):
+ raise ValueError("Parameter link for " +
+ "{} is not simple".format(par.fullname))
+
+ # If the link is to a "unknown" component then we could
+ # iterate through all the source expressions to find
+ # the relevant data sets, and hence annuli, but leave
+ # that for a later revision since the current assumption
+ # is that all source components are handled by deproject.
+ #
+ linkname = par.link.modelname
+ try:
+ lshell = cpt_map[linkname]
+ except KeyError:
+ raise RuntimeError("Model component " +
+ "{} is not ".format(linkname) +
+ "managed by deproject")
+
+ graph.add_link(shell, lshell)
+
+ # Rely on the shell numbering to be numeric and in ascending
+ # order to get the list of shells that must be fit together.
+ #
+ fit_groups = sorted([sorted(grp) for grp in graph.get_groups()])
+
+ # It is possible for the groups to be invalid here, in that
+ # the user could have tied together annuli 2 and 4, but not
+ # 3, which breaks the onion-skin approach.
+ #
+ for grp in fit_groups:
+ # ensure that the membership is n, n+1, ..., m with no
+ # gaps.
+ if len(grp) == 1:
+ continue
+
+ grp = numpy.asarray(grp)
+ dg = grp[1:] - grp[:-1]
+ if numpy.any(dg != 1):
+ raise ValueError("Non-consecutive annuli are " +
+ "tied together: {}".format(grp))
+
+ # What datasets are used for each group?
+ #
+ out = []
+ for anns in fit_groups:
+
+ # What dataset ids are associated with these annuli
+ dataids = [x['id'] for x in self.datasets
+ if x['annulus'] in anns]
+
+ out.append({'annuli': anns, 'dataids': dataids})
+
+ return list(reversed(out))
+
+ def get_radii(self, units='arcsec'):
+ """What are the radii of the shells?
+
+ Return the inner and outer edge of each annulus, in the given
+ units. Physical units (e.g. 'kpc') can only be used if a redshift or
+ angular-diameter distance has been set. This does not apply the
+ grouping that `get_shells` does.
+
+ Parameters
+ ----------
+ units : str or astropy.units.Unit, optional
+ The name of the units to use for the returned radii. They must
+ be an angle - such as 'arcsec' - or a length - such as 'kpc'
+ or 'Mpc' (case is important).
+
+ See Also
+ --------
+ get_shells
+
+ Returns
+ -------
+ rlo, rhi : astropy.units.Quantity, astropy.units.Quantity
+ The inner and outer radius for each annulus.
+ """
+
+ # Do we know about this unit? Give a slightly-more helpful
+ # message than the default from the AstroPy parser.
+ #
+ try:
+ unit = u.Unit(units)
+ except ValueError:
+ raise ValueError("Invalid unit: expected a value like " +
+ "'arcsec' or 'kpc'.")
+
+ radii = self.radii.copy()
+
+ if unit.physical_type == 'angle':
+ radii = radii.to(unit)
+
+ elif unit.physical_type == 'length':
+ # Treat the angular distance value as having length / radian
+ rscale = self.angdist / (1 * u.radian)
+
+ # This would convert to m
+ # radii = (radii * rscale).decompose()
+
+ radii = (radii * rscale).to(unit)
+
+ else:
+ raise u.UnitConversionError("Must be given an angle or length")
+
+ return radii[:-1], radii[1:]
+
+ def guess(self):
+ """Guess the starting point by fitting the projected data.
+
+ Use a fitting scheme - based on the suggestion in the XSPEC projct
+ documention - to estimate the starting position of the fit (the
+ initial fit parameters). This can be useful since it can reduce
+ the time taken to fit the deprojected data and help avoid
+ the deprojection from getting stuck in a local minimum.
+
+ See Also
+ --------
+ fit
+
+ Notes
+ -----
+
+ Each annulus, from outer to inner, is fit individually, ignoring
+ the contribution from any outer annulus. After the fit, the
+ model normalisation is corrected for the volume-filling factor of
+ the annulus. If there are any tied parameters between annuli then
+ these annuli are combined together (fit simultaneously).
+
+ Unlike the Sherpa guess function, this does *not* change the
+ limits of any parameter.
+
+ Possible improvements include:
+ - re-normalize each spectrum before fitting.
+ - transfer the model parameters of the inner-most shell in a
+ group to the next set of shells to fit.
+
+ """
+
+ groups = self.get_shells()
+ ngroups = len(groups)
+ assert (ngroups > 0) & (ngroups <= self.nshell)
+
+ if ngroups != self.nshell:
+ print("Note: annuli have been tied together")
+
+ for group in groups:
+
+ annuli = group['annuli']
+ nannuli = len(annuli)
+ assert nannuli > 0
+ dataids = group['dataids']
+
+ msg = 'Projected fit to '
+ if len(annuli) == 1:
+ msg += 'annulus {} '.format(annuli[0])
+ else:
+ msg += 'annuli {} '.format(annuli)
+
+ if len(dataids) == 1:
+ msg += 'dataset: {}'.format(dataids[0])
+ else:
+ msg += ' datasets: {}'.format(dataids)
+
+ print(msg)
+
+ orig_models = [(did, ui.get_source(did)) for did in dataids]
+
+ # perhaps this logic should be packaged up
+ shells = dict([(x['id'], x['annulus']) for x in self.datasets])
+
+ try:
+ # Is there a better way to re-create the "base" model?
+ #
+ for did in dataids:
+ srcmodel = self.srcmodel
+ shell = shells[did]
+
+ for model_comp in reversed(self.srcmodel_comps):
+ i0 = model_comp['start']
+ i1 = model_comp['end']
+ model_comp_name = self._create_name(model_comp['type'],
+ shell)
+ srcmodel = srcmodel[:i0] + model_comp_name + srcmodel[i1:]
+
+ ui.set_source(did, srcmodel)
+
+ # TODO: run renormalize on each dataset before the fit
+
+ ui.fit(*dataids)
+
+ finally:
+ for did, smdl in orig_models:
+ ui.set_source(did, smdl)
+
+ # Correct the normalization
+ #
+ fs = self.vol_norm.diagonal()
+ for did in dataids:
+ shell = shells[did]
+ f = fs[shell]
+
+ for mdl in self.model_comps:
+ if mdl['shell'] != shell:
+ continue
+
+ for par in mdl['object'].pars:
+ if par.name != 'norm':
+ continue
+
+ par.val /= f
+
+ def _freeze_model_pars(self):
+ """Freeze, and return, all thawed parameters in the fit
+
+ Returns
+ -------
+ pars : list of sherpa.parameter.Parameter instances
+ The parameters that have been frozen.
+
+ See Also
+ --------
+ _thaw_model_pars
+
+ """
+
+ out = []
+ for model_comp in self.model_comps:
+ out.extend([p for p in model_comp['object'].pars
+ if not p.frozen])
+
+ return out
+
+ def _thaw_model_pars(self, pars, message=False):
+ """Thaw the parameters, with optional screen message.
+
+ Parameters
+ ----------
+ pars : list of sherpa.parameter.Parameter instances
+ The parameters to thaw. Note that thaw is called whatever the
+ state of a parameter.
+ message : bool, optional
+ If True then a screen message is displayed for each parameter.
+
+ See Also
+ --------
+ _freeze_model_pars
+
+ """
+
+ for par in pars:
+ if message:
+ print('Thawing {}'.format(par.fullname))
+
+ par.thaw()
+
+ def _apply_per_group(self, verb, store):
+ """Apply a procedure per onion-skin group (from outer to inner).
+
+ For each shell - run from outer to inner - apply the onion-skin
+ approach (so free up the shell but freeze the contribution from
+ outer shells) to runfunc - which is given the data ids to use,
+ and then store the results of getfunc. Once all the shells have
+ been processed return a structure containing the results.
+
+ Parameters
+ ----------
+ verb : str
+ This is the first part of the message displayed to users,
+ per annulus.
+ store : fieldstore.FieldStore instance
+ This object will run the function and then parse the output
+ into the necessary form.
+
+ Returns
+ -------
+ rvals : astropy.table.Table instance
+ The data, as a set of columns. The choice of columns is
+ controlled by the `store` object. Additional columns include
+ 'annulus', 'rlo_ang', 'rhi_ang', 'rlo_phys', 'rhi_phys',
+ 'density', and optionally 'density_lo' and 'density_hi'.
+
+ Notes
+ -----
+ Any parameter links result in annuli being grouped together
+ for a fit: that is, the fit will be a simultaneous fit to all
+ the datasets associated with the set of annuli.
+
+ It would be useful to add in extra metadata to the output, and
+ take advantage of the units support in AstroPy where possible.
+ """
+
+ groups = self.get_shells()
+ ngroups = len(groups)
+ assert (ngroups > 0) & (ngroups <= self.nshell)
+
+ if ngroups != self.nshell:
+ print("Note: annuli have been tied together")
+
+ # Find all the thawed parameters so that they can be
+ # restored at the end of the fit, or in case of an error.
+ #
+ thawed = self._freeze_model_pars()
+
+ # We need to be able to map from dataset id to annulus, and
+ # from annulus to model components.
+ #
+ annulusmap = {x['id']: x['annulus'] for x in self.datasets}
+
+ # Get a list of all the annuli, in ascending order.
+ all_annuli = sorted(list({x['annulus'] for x in self.datasets}))
+
+ componentmap = defaultdict(list)
+ for mcomp in self.model_comps:
+ val = (mcomp['object'], mcomp['type'])
+ componentmap[mcomp['shell']].append(val)
+
+ # Store the data as a dictionary of arrays. It would be useful
+ # if the FieldStore instance could handle this - since it
+ # "knows" what the columns are going to be - but the current
+ # implementation calculates these columns when the run method
+ # is called (i.e. at run time), rather than before the
+ # onion-peel approach is called. It is possible that a
+ # re-design would be helpful (such as calculate the parameter
+ # names at the start, which would have other useful consequences
+ # once more-complicated model expressions are sorted), and then
+ # pass that information along, but for now let's see how this
+ # works.
+ #
+ out = OrderedDict()
+ try:
+ for group in groups:
+
+ annuli = group['annuli']
+ nannuli = len(annuli)
+ assert nannuli > 0
+ dataids = group['dataids']
+
+ msg = '{} '.format(verb)
+ if nannuli == 1:
+ msg += 'annulus {} '.format(annuli[0])
+ else:
+ msg += 'annuli {} '.format(annuli)
+
+ if len(dataids) == 1:
+ msg += ' dataset: {}'.format(dataids[0])
+ else:
+ msg += ' datasets: {}'.format(dataids)
+
+ print(msg)
+
+ res = store.run(annulusmap, componentmap, *dataids)
+ assert len(res) == nannuli
+
+ # Extract out the per-shell results and add to the
+ # output.
+ #
+ for shell in annuli:
+
+ if len(out) == 0:
+ # Note, add in extra fields to the stored fields
+ #
+ out['annulus'] = all_annuli
+
+ rlo_ang, rhi_ang = self.get_radii(units='arcsec')
+ out['rlo_ang'] = rlo_ang
+ out['rhi_ang'] = rhi_ang
+
+ rlo_phys, rhi_phys = self.get_radii(units='kpc')
+ out['rlo_phys'] = rlo_phys
+ out['rhi_phys'] = rhi_phys
+
+ for field in res[shell]:
+ out[field] = [None] * self.nshell
+
+ for field, value in res[shell].items():
+ assert out[field][shell] is None, shell
+ out[field][shell] = value
+
+ for model_comp in self.model_comps:
+ # Freeze the current annulus
+ if model_comp['shell'] in annuli:
+ print('Freezing {}'.format(model_comp['name']))
+ ui.freeze(model_comp['object'])
+
+ finally:
+ self._thaw_model_pars(thawed, message=True)
+
+ for k, vs in out.items():
+ assert vs is not None, k
+
+ # It's useful to have NumPy arrays for some of the following
+ # calculations, but do not convert those that already have
+ # units attached.
+ #
+ # Convert from None to numpy.NaN (which is assumed to
+ # only occur in an error colum (k ends in _lo or _hi)
+ # but this restriction is not checked
+ #
+ if not isinstance(vs, numpy.ndarray):
+ out[k] = numpy.asarray([numpy.nan if v is None else v
+ for v in vs])
+
+ # Add in the density calculation.
+ #
+ # The norm field should be identified at 'set_source' time
+ # so that it doesn't have to be re-discovered each time.
+ #
+ # For now look for .norm values in the returned structure,
+ # but could do this a number of ways.
+ #
+ normpars = [n for n in out.keys() if n.endswith('.norm')]
+ if len(normpars) == 0:
+ raise RuntimeError("Unable to find norm parameter!")
+ elif len(normpars) > 1:
+ raise RuntimeError("Multiple norm parameters found!")
+
+ normpar = normpars[0]
+ norms = out[normpar]
+
+ out['density'] = self._calc_density(norms)
+
+ normparlo = '{}_lo'.format(normpar)
+ normparhi = '{}_hi'.format(normpar)
+ try:
+ normlos = out[normparlo]
+ normhis = out[normparhi]
+
+ out['density_lo'] = self._calc_density(norms + normlos) - \
+ out['density']
+ out['density_hi'] = self._calc_density(norms + normhis) - \
+ out['density']
+
+ except KeyError:
+ pass
+
+ return Table(out)
+
+ def fit(self):
+ """Fit the data using the "onion-peeling" method.
+
+ Unlike the normal Sherpa fit, this does not fit all the data
+ simultaneously, but instead fits the outermost annulus first,
+ then freezes its parameters and fits the annulus inside it,
+ repeating this until all annuli have been fit. At the end of
+ the fit all the parameters that were frozen are freed. The
+ results can also be retrieved with ``get_fit_results``.
+
+ Returns
+ -------
+ fits : astropy.table.Table instance
+ This records per-annulus data, such as the inner and outer
+ radius (`rlo_ang`, `rhi_ang`, `rlo_phys`, `rhi_phys`), the
+ final fit statistic and change in fit statistic (`statval` and
+ `dstatval`), the reduced statistic and q value (as `rstat`
+ and `qval`) if appropriate, and the thawed parameter values
+ (accessed using . syntax, where the
+ match is case sensitive).
+
+ See Also
+ --------
+ conf, covar, get_fit_results, guess, fit_plot
+
+ Notes
+ -----
+
+ If there are any tied parameters between annuli then these annuli
+ are combined together (fit simultaneously). The results from the
+ fits to each annulus can be retrieved after ``fit`` has been called
+ with the ``get_fit_results`` method.
+
+ The results have been separated out per annulus, even if several
+ annuli were combined in a fit due to tied parameters, and there is
+ no information in the returned structure to note this.
+
+ Examples
+ --------
+
+ Fit the annuli using the onion-peeling approach, and then plot
+ up the reduced statistic for each dataset:
+
+ >>> res = dep.fit()
+ >>> plt.clf()
+ >>> rmid = 0.5 * (res['rlo_phys'] + res['rhi_phys'])
+ >>> plt.plot(rmid, res['rstat'])
+
+ Plot the temperature-abundance values per shell, color-coded
+ by annulus:
+
+ >>> plt.clf()
+ >>> plt.plot(res['xsapec.kT'], res['xsapec.Abundanc'],
+ ... c=res['annulus'])
+ >>> plt.colorbar()
+ >>> plt.xlabel('kT')
+ >>> plt.ylabel('Abundance')
+
+ Plot up the temperature distibution as a function of radius
+ from the fit::
+
+ >>> dep.fit()
+ >>> dep.fit_plot('xsmekal.kt')
+
+ """
+
+ # For now return nothing
+ self._fit_results = None
+ self._fit_results = self._apply_per_group('Fitting',
+ fieldstore.FitStore())
+ return self._fit_results
+
+ def covar(self):
+ """Estimate errors using covariance, using the "onion-peeling" method.
+
+ It is assumed that ``fit`` has been called. The results can also be
+ retrieved with ``get_covar_results``.
+
+ Returns
+ -------
+ errors : astropy.table.Table instance
+ This records per-annulus data, such as the inner and outer
+ radius (`rlo_ang`, `rhi_ang`, `rlo_phys`, `rhi_phys`), the
+ sigma and percent values, and parameter results (accessed
+ using ., ._lo,
+ and ._hi syntax, where the match is
+ case sensitive). The _lo and _hi values are symmetric for
+ covar, that is the _lo value will be the negative of the
+ _hi value.
+
+ See Also
+ --------
+ conf, fit, get_covar_results, covar_plot
+
+ Examples
+ --------
+
+ Run a fit and then error analysis, then plot up the abundance
+ against temperature values including the error bars. Since
+ the covariance routine returns symmetric error bars, the
+ _hi values are used in the plot::
+
+ >>> dep.fit()
+ >>> errs = dep.covar()
+ >>> kt, abund = errs['xsapec.kT'], errs['xsapec.Abundanc']
+ >>> dkt = errs['xsapec.kT_hi']
+ >>> dabund = errs['xsapec.Abundanc_hi']
+ >>> plt.clf()
+ >>> plt.errorbar(kt, abund, xerr=dkt, yerr=dabund, fmt='.')
+
+ Plot up the temperature distibution as a function of radius,
+ including the error bars calculated by the covar routine::
+
+ >>> dep.fit()
+ >>> dep.covar()
+ >>> dep.covar_plot('xsmekal.kt')
+
+ """
+
+ self._covar_results = None
+ self._covar_results = self._apply_per_group('Covariance for',
+ fieldstore.CovarStore())
+ return self._covar_results
+
+ def conf(self):
+ """Estimate errors using confidence, using the "onion-peeling" method.
+
+ It is assumed that ``fit`` has been called. The results can also be
+ retrieved with ``get_conf_results``.
+
+ Returns
+ -------
+ errors : astropy.table.Table instance
+ This records per-annulus data, such as the inner and outer
+ radius (`rlo_ang`, `rhi_ang`, `rlo_phys`, `rhi_phys`), the
+ sigma and percent values, and parameter results (accessed
+ using ., ._lo,
+ and ._hi syntax, where the match is
+ case sensitive).
+
+ See Also
+ --------
+ covar, fit, get_conf_results, conf_plot
+
+ Examples
+ --------
+
+ Run a fit and then error analysis, then plot up the abundance
+ against temperature values including the error bars. Note that
+ the Matplotlib `errorbar` routine requires "positive" error values
+ whereas the _lo values are negative, hence they are
+ negated in the creation of ``dkt`` and ``dabund``::
+
+ >>> dep.fit()
+ >>> errs = dep.conf()
+ >>> kt, abund = errs['xsapec.kT'], errs['xsapec.Abundanc']
+ >>> ktlo, kthi = errs['xsapec.kT_lo'], errs['xsapec.kT_hi']
+ >>> ablo, abhi = errs['xsapec.Abundanc_lo'], errs['xsapec.Abundanc_hi']
+ >>> dkt = np.vstack((-ktlo, kthi))
+ >>> dabund = np.vstack((-ablo, abhi))
+ >>> plt.clf()
+ >>> plt.errorbar(kt, abund, xerr=dkt, yerr=dabund, fmt='.')
+
+ Plot up the temperature distibution as a function of radius,
+ including the error bars calculated by the conf routine::
+
+ >>> dep.fit()
+ >>> dep.conf()
+ >>> dep.conf_plot('xsmekal.kt')
+
+ """
+
+ self._conf_results = None
+ self._conf_results = self._apply_per_group('Confidence for',
+ fieldstore.ConfStore())
+ return self._conf_results
+
+ def get_fit_results(self):
+ """What are the fit results, per annulus?
+
+ This returns the fit result for each annulus from the last time
+ that the ``fit`` method was called. It *does not* check to see
+ if anything has changed since the last ``fit`` call (e.g.
+ parameters being tied together or untied, or a manual fit
+ to a shell). Note that ``get_shells`` should be used to find out
+ if the shells were grouped together.
+
+ Returns
+ -------
+ fits : astropy.table.Table instance
+ This records per-annulus data, such as the inner and outer
+ radius (`rlo_ang`, `rhi_ang`, `rlo_phys`, `rhi_phys`), the
+ final fit statistic and change in fit statistic (`statval` and
+ `dstatval`), the reduced statistic and q value (as `rstat`
+ and `qval`) if appropriate, and the thawed parameter values
+ (accessed using . syntax, where the
+ match is case sensitive).
+
+ See Also
+ --------
+ fit, get_conf_results, get_covar_results, get_radii, get_shells,
+ fit_plot
+
+ """
+
+ if self._fit_results is None:
+ raise ValueError("The fit method has not been called")
+
+ return copy.deepcopy(self._fit_results)
+
+ def get_covar_results(self):
+ """What are the covar results, per annulus?
+
+ This returns the fit result for each annulus from the last time
+ that the ``covar`` method was called. It *does not* check to see
+ if anything has changed since the last ``covar`` call (e.g.
+ parameters being tied together or untied, or a manual fit
+ to a shell). Note that ``get_shells`` should be used to find out
+ if the shells were grouped together.
+
+ Returns
+ -------
+ errors : astropy.table.Table instance
+ This records per-annulus data, such as the inner and outer
+ radius (`rlo_ang`, `rhi_ang`, `rlo_phys`, `rhi_phys`), the
+ sigma and percent values, and parameter results (accessed
+ using ., ._lo,
+ and ._hi syntax, where the match is
+ case sensitive).
+
+ See Also
+ --------
+ fit, get_conf_results, get_fit_results, get_radii, get_shells,
+ covar_plot
+
+ """
+
+ if self._covar_results is None:
+ raise ValueError("The covar method has not been called")
+
+ return copy.deepcopy(self._covar_results)
+
+ def get_conf_results(self):
+ """What are the conf results, per annulus?
+
+ This returns the fit result for each annulus from the last time
+ that the ``conf`` method was called. It *does not* check to see
+ if anything has changed since the last ``conf`` call (e.g.
+ parameters being tied together or untied, or a manual fit
+ to a shell). Note that ``get_shells`` should be used to find out
+ if the shells were grouped together (although this can be
+ reconstructed from the `datasets` field of each `ErrorEstResults`
+ instance).
+
+ Returns
+ -------
+ errors : astropy.table.Table instance
+ This records per-annulus data, such as the inner and outer
+ radius (`rlo_ang`, `rhi_ang`, `rlo_phys`, `rhi_phys`), the
+ sigma and percent values, and parameter results (accessed
+ using ., ._lo,
+ and ._hi syntax, where the match is
+ case sensitive).
+
+ See Also
+ --------
+ fit, get_covar_results, get_fit_results, get_radii, get_shells,
+ conf_plot
+
+ """
+
+ if self._conf_results is None:
+ raise ValueError("The conf method has not been called")
+
+ return copy.deepcopy(self._conf_results)
+
+ def _calc_density(self, norms, ne_nh_ratio=1.18):
+ """Calculate the electron density for each shell.
+
+ This performs the calculation described in `get_density`.
+
+ Parameters
+ ----------
+ norms : sequence of float
+ The normalization values, in annulus order.
+ ne_hh_ratio : float, optional
+ The n_e to n_h ratio (default 1.18).
+
+ Returns
+ -------
+ dens : astropy.units.quantity.Quantity instance
+ The densities calculated for each shell, in units of cm^-3.
+
+ """
+
+ if len(norms) != self.nshell:
+ raise ValueError("norms has wrong length")
+
+ # Manual descontruction/reconstruction of units
+ #
+ DA_cm = self.angdist.to_value(u.cm)
+ rmax_rad = self.radii[-1].to_value(u.rad)
+ z = self.redshift
+
+ r_sphere = rmax_rad * DA_cm
+
+ # volume of sphere enclosing outer shell (cm^3)
+ #
+ # volume = 4 * pi / 3 * r_sphere**3
+ # factor = 4 * pi * DA_cm**2 * 1e14 * (1.0 + z)**2 / volume * ne_nh_ratio
+ #
+ # and after manual cancellation of the 4 pi terms
+ #
+ vterm = 1.0 / 3 * r_sphere**3
+ factor = DA_cm**2 * 1e14 * (1.0 + z)**2 / vterm * ne_nh_ratio
+
+ return numpy.sqrt(factor * numpy.asarray(norms)) * u.cm**(-3)
+
+ def get_density(self):
+ """Calculate the electron density for each shell.
+
+ Convert the model normalzations (assumed to match the standard
+ definition for XSPEC thermal-plasma models) for each shell.
+
+ Returns
+ -------
+ dens : astropy.units.quantity.Quantity instance
+ The densities calculated for each shell, in units of cm^-3.
+
+ See Also
+ --------
+ find_norm
+
+ Notes
+ -----
+ The electron density is taken to be::
+
+ n_e^2 = norm * 4*pi * DA^2 * 1e14 * (1+z)^2 / volume * ne_nh_ratio
+
+ where::
+
+ norm = model normalization from sherpa fit
+ DA = angular size distance (cm)
+ volume = volume (cm^3)
+ ne_nh_ratio = 1.18
+
+ The model components for each volume element (the intersection of the
+ annular cylinder ``a`` with the spherical shell ``s``) are multiplied
+ by a volume normalization::
+
+ vol_norm[s,a] = volume[s,a] / v_sphere
+ v_sphere = volume of sphere enclosing outer annulus
+
+ With this convention the ``volume`` used in calculating the electron
+ density is simply ``v_sphere``.
+
+ """
+
+ norms = [self.find_norm(s) for s in range(self.nshell)]
+ return self._calc_density(norms)
+
+ def _radial_plot(self, plottitle, xunits, ys, ylabel,
+ dys=None,
+ xlog=True, ylog=False,
+ overplot=False, clearwindow=True):
+ """Create a plot of the data versus radius (of the annuli).
+
+ Parameters
+ ----------
+ plottitle : str
+ The title for the plot
+ xunits : str or astropy.units.Unit
+ The X-axis units (a length or angle, such as 'Mpc' or
+ 'arcsec', where the case is important).
+ ys : sequence of float
+ The Y values to plot (must be in annuli order).
+ ylabel : str
+ The label for the Y axis
+ dys : None or ndarray, optional
+ The error bars on the y axis. This can be None or a ndarray
+ of one or two dimensions (N points or N by 2).
+ xlog : bool, optional
+ Should the x axis be drawn with a log scale (default True)?
+ ylog : bool, optional
+ Should the y axis be drawn with a log scale (default False)?
+ overplot : bool, optional
+ Clear the plot or add to existing plot?
+ clearwindow : bool, optional
+ How does this interact with overplot?
+
+ """
+
+ rlo, rhi = self.get_radii(units=xunits)
+
+ # drop units support immediately as ChIPS doesn't recognize
+ # this (can support them in matplotlib, but given the
+ # Sherpa plotting API it isn't clear how well supported it
+ # would be)
+ #
+ rmid = (rlo.value + rhi.value) / 2
+ dr = rhi.value - rlo.value
+
+ # Attempt to handle LaTeX differences between the backends,
+ # but the support is *very* limited so may not work here.
+ #
+ # The aim is to support the matplotlib backend, with minimal
+ # support for ChIPS.
+ #
+ xlabel = _add_unit('Radius', rlo)
+ if ylabel.find('_') > -1 or ylabel.find('^') > -1:
+ # Unfortunately the matplotlib version is a "global"
+ # check, so doesn't check if parts of the term are
+ # already enclosed in '$'. This is a problem for those
+ # labels that have AstroPy unit strings, since they
+ # have already been protected.
+ #
+ if not (plotter.name == 'pylab' and ylabel.find('$') > -1):
+ ylabel = plotter.get_latex_for_string(ylabel)
+
+ prefs = plotter.get_data_plot_defaults()
+ prefs['xerrorbars'] = True
+
+ # We handle error bars manually for ChIPS (it has to be done
+ # for asymmetric Y errors, but there also seems to be issues
+ # with the X axis errors not being drawn which I do not
+ # want to investigate too much just right now).
+ #
+ manual_errors = plotter.name == 'chips' and dys is not None
+
+ prefs['yerrorbars'] = dys is not None
+ if manual_errors:
+ prefs['yerrorbars'] = False
+
+ prefs['xlog'] = xlog
+ prefs['ylog'] = ylog
+
+ # Access the underlying plot machinery directly, rather than
+ # use the sherpa.plot.DataPlot object, since Sherpa does not
+ # support asymmetric errors but plotter.plot does, at least
+ # for the pylab backend.
+ #
+ try:
+ plotter.begin()
+ plotter.plot(rmid, ys, dys, dr, plottitle, xlabel, ylabel,
+ overplot, clearwindow, **prefs)
+
+ # For some reason the X error bar isn't being drawn with ChIPS
+ # so force it.
+ #
+ if plotter.name == 'chips':
+ import pychips
+
+ # Assume the current curve is the data we have just plotted
+ # and we do not want to replot the symbol.
+ #
+ crv = pychips.get_curve()
+ crv.symbol.style = 'none'
+ crv.err.up = True
+ crv.err.down = True
+ crv.err.left = True
+ crv.err.right = True
+
+ ndim = numpy.asarray(dys).ndim
+ if ndim == 2:
+ dylo = dys[0]
+ dyhi = dys[1]
+ elif ndim == 1:
+ dylo = dys
+ dyhi = dys
+ else:
+ dylo = None
+ dyhi = None
+
+ errs = [dylo, dyhi, dr / 2, dr / 2]
+ pychips.add_curve(rmid, ys, errs, crv)
+
+ except BaseException as exc:
+ plotter.exceptions()
+ raise exc
+
+ else:
+ plotter.end()
+
+ def par_plot(self, par, units='kpc',
+ xlog=True, ylog=False,
+ overplot=False, clearwindow=True):
+ """Plot up the parameter as a function of radius.
+
+ This plots up the current parameter values. The ``fit_plot``,
+ ``conf_plot``, and ``covar_plot`` routines display the fit
+ and error results for these parameters.
+
+ Parameters
+ ----------
+ par : str
+ The parameter name, specified as ..
+ units : str or astropy.units.Unit, optional
+ The X-axis units (a length or angle, such as 'Mpc' or
+ 'arcsec', where the case is important).
+ xlog : bool, optional
+ Should the x axis be drawn with a log scale (default True)?
+ ylog : bool, optional
+ Should the y axis be drawn with a log scale (default False)?
+ overplot : bool, optional
+ Clear the plot or add to existing plot?
+ clearwindow : bool, optional
+ How does this interact with overplot?
+
+ See Also
+ --------
+ conf_plot, covar_plot, density_plot, fit_plot
+
+ Examples
+ --------
+
+ Plot the temperature as a function of radius.
+
+ >>> dep.par_plot('xsapec.kt')
+
+ Label the radii with units of arcminutes for the abundanc
+ parameter of the xsapec model:
+
+ >>> dep.par_plot('xsapec.abundanc', units='arcmin')
+
+ """
+
+ # Assume par is "model_name.par_name" and we do not have to
+ # worry about case for model_name, but may have to for par_name
+ #
+ mname, pname = self._split_parname(par)
+ pname = pname.lower()
+ cpts = [cpt['object'] for cpt in self.model_comps
+ if cpt['type'] == mname]
+
+ # Probably can not get here and this happen (thanks to the get_par)
+ # call, but just in case
+ if len(cpts) == 0:
+ raise ValueError("No matching model {} for par={}".format(mname,
+ par))
+
+ # Assume they are all the same (they better be)
+ #
+ yunits = None
+ for p in cpts[0].pars:
+ if p.name.lower() != pname:
+ continue
+
+ yunits = p.units
+ break
+
+ # Also should not happen, so report if it does but do not
+ # error out
+ if yunits is None:
+ print("WARNING: unable to find match for parameter {}".format(par))
+ yunits = ''
+
+ ylabel = par
+ if yunits.strip() != '':
+ ylabel += " ({})".format(yunits)
+
+ pvals = self.get_par(par)
+ self._radial_plot(par, units, pvals, ylabel,
+ xlog=xlog, ylog=ylog,
+ overplot=overplot, clearwindow=clearwindow)
+
+ def density_plot(self, units='kpc',
+ xlog=True, ylog=True,
+ overplot=False, clearwindow=True):
+ """Plot up the electron density as a function of radius.
+
+ The density is displayed with units of cm^-3. This plots up the
+ density calculated using the current normalization parameter
+ values. The ``fit_plot``, ``conf_plot``, and ``covar_plot``
+ routines display the fit and error results for these parameters.
+
+ Parameters
+ ----------
+ units : str or astropy.units.Unit, optional
+ The X-axis units (a length or angle, such as 'Mpc' or
+ 'arcsec', where the case is important).
+ xlog : bool, optional
+ Should the x axis be drawn with a log scale (default True)?
+ ylog : bool, optional
+ Should the y axis be drawn with a log scale (default False)?
+ overplot : bool, optional
+ Clear the plot or add to existing plot?
+ clearwindow : bool, optional
+ How does this interact with overplot?
+
+ See Also
+ --------
+ conf_plot, covar_plot, fit_plot, par_plot
+
+ Examples
+ --------
+
+ Plot the density as a function of radius.
+
+ >>> dep.density_plot()
+
+ Label the radii with units of arcminutes:
+
+ >>> dep.density_plot(units='arcmin')
+
+ """
+
+ nes = self.get_density().value
+
+ # Unfortunately the LaTeX emulation in the two backends is not
+ # comparable, which limits the fidelity of the label.
+ #
+ # ylabel = 'n$_e$ (cm$^{-3}$)'
+ ylabel = r'n_e\ (\mathrm{cm^{-3}})'
+
+ self._radial_plot('density', units, nes, ylabel,
+ xlog=xlog, ylog=ylog,
+ overplot=overplot, clearwindow=clearwindow)
+
+ def fit_plot(self, field, results=None,
+ units='kpc',
+ xlog=True, ylog=False,
+ overplot=False, clearwindow=True):
+ """Plot up the fit results as a function of radius.
+
+ This method can be used to plot up the last fit results or
+ a previously-stored set. To include error bars on the
+ dependent values use the `conf_plot` or `covar_plot` methods.
+
+ Parameters
+ ----------
+ field : str
+ The column to plot from the fit results (the match is case
+ insensitive).
+ results : None or astropy.table.Table instance
+ The return value from the ``fit`` or ``get_fit_results``
+ methods.
+ units : str or astropy.units.Unit, optional
+ The X-axis units (a length or angle, such as 'Mpc' or
+ 'arcsec', where the case is important).
+ xlog : bool, optional
+ Should the x axis be drawn with a log scale (default True)?
+ ylog : bool, optional
+ Should the y axis be drawn with a log scale (default False)?
+ overplot : bool, optional
+ Clear the plot or add to existing plot?
+ clearwindow : bool, optional
+ How does this interact with overplot?
+
+ See Also
+ --------
+ fit, get_fit_results, conf_plot, covar_plot, density_plot, par_plot
+
+ Examples
+ --------
+
+ Plot the temperature as a function of radius from the last
+ fit:
+
+ >>> dep.fit_plot('xsapec.kt')
+
+ Plot the reduced fit statistic from the last fit:
+
+ >>> dep.fit_plot('rstat')
+
+ Plot the density with the radii labelled in arcminutes and the
+ density shown on a log scale:
+
+ >>> dep.fit_plot('density', units='arcmin', ylog=True)
+
+ Overplot the current fit results on those from a previous fit,
+ where ``fit1`` was returned from the ``fit`` or ``get_fit_results``
+ methods:
+
+ >>> dep.fit_plot('xsapec.abundanc', results=fit1)
+ >>> dep.fit_plot('xsapec.abundanc', overplot=True)
+
+ """
+
+ if results is None:
+ plotdata = self.get_fit_results()
+ else:
+ plotdata = results
+
+ try:
+ ys = plotdata[field]
+ except KeyError:
+ flower = field.lower()
+ names = [n for n in plotdata.keys() if n.lower() == flower]
+ if len(names) == 0:
+ raise ValueError("Unrecognized field {}".format(field))
+ elif len(names) > 1:
+ raise RuntimeError("Multiple fields match {}".format(field))
+
+ field = names[0]
+ ys = plotdata[field]
+
+ ylabel = _add_unit(field, ys)
+ self._radial_plot(field, units, ys, ylabel,
+ xlog=xlog, ylog=ylog,
+ overplot=overplot, clearwindow=clearwindow)
+
+ def conf_plot(self, field, results=None,
+ units='kpc',
+ xlog=True, ylog=False,
+ overplot=False, clearwindow=True):
+ """Plot up the confidence errors as a function of radius.
+
+ This method can be used to plot up the last conf results or
+ a previously-stored set. Any error bars are shown at the
+ scale they were calculated (as given by the ``sigma`` and
+ ``percent`` columns of the results).
+
+ Parameters
+ ----------
+ field : str
+ The column to plot from the fit results (the match is case
+ insensitive).
+ results : None or astropy.table.Table instance
+ The return value from the ``conf`` or ``get_conf_results``
+ methods.
+ units : str or astropy.units.Unit, optional
+ The X-axis units (a length or angle, such as 'Mpc' or
+ 'arcsec', where the case is important).
+ xlog : bool, optional
+ Should the x axis be drawn with a log scale (default True)?
+ ylog : bool, optional
+ Should the y axis be drawn with a log scale (default False)?
+ overplot : bool, optional
+ Clear the plot or add to existing plot?
+ clearwindow : bool, optional
+ How does this interact with overplot?
+
+ See Also
+ --------
+ fit, get_conf_results, fit_plot, covar_plot, density_plot, par_plot
+
+ Notes
+ -----
+ Error bars are included on the dependent axis if the results
+ contain columns that match the requested field with suffixes
+ of '_lo' and '_hi'. These error bars are asymmetric, which is
+ different to ``covar_plot``.
+
+ If a limit is missing (i.e. it is a NaN) then no error bar is
+ drawn. This can make it look like the error is very small.
+
+ Examples
+ --------
+
+ Plot the temperature as a function of radius from the last
+ fit, including error bars:
+
+ >>> dep.conf_plot('xsapec.kt')
+
+ Plot the density with the radii labelled in arcminutes and the
+ density shown on a log scale:
+
+ >>> dep.conf_plot('density', units='arcmin', ylog=True)
+
+ Overplot the current conf results on those from a previous fit,
+ where ``conf1`` was returned from the ``conf`` or ``get_conf_results``
+ methods:
+
+ >>> dep.conf_plot('xsapec.abundanc', results=conf1)
+ >>> dep.conf_plot('xsapec.abundanc', overplot=True)
+
+ """
+
+ if results is None:
+ plotdata = self.get_conf_results()
+ else:
+ plotdata = results
+
+ try:
+ ys = plotdata[field]
+ except KeyError:
+ flower = field.lower()
+ names = [n for n in plotdata.keys() if n.lower() == flower]
+ if len(names) == 0:
+ raise ValueError("Unrecognized field {}".format(field))
+ elif len(names) > 1:
+ raise RuntimeError("Multiple fields match {}".format(field))
+
+ field = names[0]
+ ys = plotdata[field]
+
+ try:
+ flo = '{}_lo'.format(field)
+ fhi = '{}_hi'.format(field)
+ dys = numpy.vstack((-plotdata[flo], plotdata[fhi]))
+ except KeyError:
+ dys = None
+
+ ylabel = _add_unit(field, ys)
+ self._radial_plot(field, units, ys, ylabel, dys=dys,
+ xlog=xlog, ylog=ylog,
+ overplot=overplot, clearwindow=clearwindow)
+
+ def covar_plot(self, field, results=None,
+ units='kpc',
+ xlog=True, ylog=False,
+ overplot=False, clearwindow=True):
+ """Plot up the covariance errors as a function of radius.
+
+ This method can be used to plot up the last covar results or
+ a previously-stored set. Any error bars are shown at the
+ scale they were calculated (as given by the ``sigma`` and
+ ``percent`` columns of the results).
+
+ Parameters
+ ----------
+ field : str
+ The column to plot from the fit results (the match is case
+ insensitive).
+ results : None or astropy.table.Table instance
+ The return value from the ``covar`` or ``get_covar_results``
+ methods.
+ units : str or astropy.units.Unit, optional
+ The X-axis units (a length or angle, such as 'Mpc' or
+ 'arcsec', where the case is important).
+ xlog : bool, optional
+ Should the x axis be drawn with a log scale (default True)?
+ ylog : bool, optional
+ Should the y axis be drawn with a log scale (default False)?
+ overplot : bool, optional
+ Clear the plot or add to existing plot?
+ clearwindow : bool, optional
+ How does this interact with overplot?
+
+ See Also
+ --------
+ fit, get_covar_results, fit_plot, conf_plot, density_plot, par_plot
+
+ Notes
+ -----
+ Error bars are included on the dependent axis if the results
+ contain columns that match the requested field with the suffix
+ '_hi'. The error bars are therefore symmetric, which is
+ different to ``conf_plot``.
+
+ If a limit is missing (i.e. it is a NaN) then no error bar is
+ drawn. This can make it look like the error is very small.
+
+ Examples
+ --------
+
+ Plot the temperature as a function of radius from the last
+ fit, including error bars:
+
+ >>> dep.covar_plot('xsapec.kt')
+
+ Plot the density with the radii labelled in arcminutes and the
+ density shown on a log scale:
+
+ >>> dep.covar_plot('density', units='arcmin', ylog=True)
+
+ Overplot the current covar results on those from a previous fit,
+ where ``covar1`` was returned from the ``covar`` or
+ ``get_covar_results`` methods:
+
+ >>> dep.covar_plot('xsapec.abundanc', results=covar1)
+ >>> dep.covar_plot('xsapec.abundanc', overplot=True)
+
+ """
+
+ if results is None:
+ plotdata = self.get_covar_results()
+ else:
+ plotdata = results
+
+ try:
+ ys = plotdata[field]
+ except KeyError:
+ flower = field.lower()
+ names = [n for n in plotdata.keys() if n.lower() == flower]
+ if len(names) == 0:
+ raise ValueError("Unrecognized field {}".format(field))
+ elif len(names) > 1:
+ raise RuntimeError("Multiple fields match {}".format(field))
+
+ field = names[0]
+ ys = plotdata[field]
+
+ try:
+ fhi = '{}_hi'.format(field)
+ dys = plotdata[fhi]
+ except KeyError:
+ dys = None
+
+ ylabel = _add_unit(field, ys)
+ self._radial_plot(field, units, ys, ylabel, dys=dys,
+ xlog=xlog, ylog=ylog,
+ overplot=overplot, clearwindow=clearwindow)
+
+
+def _add_unit(label, xs):
+ """Add a unit string to the label if available.
+
+ This also converts a label of density to n_e.
+
+ Parameters
+ ----------
+ label : str
+ The label
+ xs : values
+ The numeric values
+
+ Returns
+ -------
+ label : str
+ If xs is an AstroPy Quantity then " ()" is added to the
+ label, using LaTeX (inline) formatting. If the input label is
+ 'density' then it is replaced by '$n_e$'.
+
+ """
+
+ if label == 'density':
+ label = '$n_e$'
+
+ try:
+ unit = xs.unit.to_string(format='latex_inline')
+ return r"{} ({})".format(label, unit)
+ except AttributeError:
+ return label
+
+
+def _get_xflt_keys(infile):
+ """Return the XFLT0001 to XFLT0005 keys in infile.
+
+ Parameters
+ ----------
+ infile : str
+ The name of the PHA file with the XFLT000n keywords.
+
+ Returns
+ -------
+ xflts : list of float
+ The XFLT0001 to XFLT0005 numbers, in order.
+
+ Raises
+ ------
+ IOError
+ If the file does not exist or a keyword is missing.
+
+ Notes
+ -----
+ The Sherpa logging level is temporarily changed to
+ the WARN level (if it is lower than this) when reading in
+ the file, to avoid excessive screen output.
+
+ """
+
+ # Hide sherpa logging here
+ olvl = _sherpa_logger.getEffectiveLevel()
+ if olvl < logging.WARN:
+ _sherpa_logger.setLevel(logging.WARN)
+
+ try:
+ pha = read_pha(infile)
+ finally:
+ _sherpa_logger.setLevel(olvl)
+
+ hdr = pha.header
+ try:
+ out = [hdr['XFLT000{}'.format(i)] for i in range(1, 6)]
+ except KeyError as e:
+ raise IOError("PHA file {} is missing keyword {}".format(infile,
+ e))
+
+ return out
+
+
+def _find_xflt_files(pat):
+ """Return the XFLT keywords in the files matching the input pattern.
+
+ This is only for circular annuli.
+
+ Parameters
+ ----------
+ pat : str
+ The pattern representing the files to use (e.g. 'ann*.pi'). If
+ the stk module (provided by CIAO) is available then CIAO stack
+ syntax (e.g. "@files.lis" to read the names from files.lis) can
+ be used. The files do not need to be in order (of increasing
+ radii).
+
+ Returns
+ -------
+ vals : list of (str, list)
+ The file name and corresponding five XFLT values. The list is
+ ordered by increasing radius of the annuli.
+
+ Raises
+ ------
+ IOError
+ When a non-circular annulus is encountered.
+
+ Notes
+ -----
+ There is no attempt to check that the annuli are consecutive and
+ not overlapping.
+ """
+
+ try:
+ import stk
+ infiles = stk.build(pat)
+ except ImportError:
+ # Assume not in a CIAO installation, so only support a glob
+ import glob
+ infiles = glob.glob(pat)
+
+ if len(infiles) == 0:
+ raise IOError("No matches to {}".format(pat))
+
+ out = []
+ for infile in infiles:
+ xflt = _get_xflt_keys(infile)
+
+ # numeric comparison not great, do we need some sort of
+ #
+ #
+ if (xflt[0] != xflt[1]) | (xflt[2] != 0):
+ raise IOError("Only circular annulli allowed: {}".format(infile))
+
+ out.append((infile, xflt))
+
+ # Sort on radius
+ #
+ return sorted(out, key=lambda x: x[1][0])
+
+
+def deproject_from_xflt(pat, rscale,
+ rinner=0,
+ angdist=None,
+ cosmology=None):
+ """Set up the projection object from XFLT keywords in the PHA files.
+
+ When using the XSPEC projct model, values are read from XFLT keywords
+ (as used by the XSPEC deprojection code [1]_) rather than being specified
+ manually. This function creates a Deproject object and loads in a set
+ of PHA files matching a pattern, using the XFLT keywords to set radii
+ and theta values. The annuli *must* be circular.
+
+ Parameters
+ ----------
+ pat : str
+ The pattern representing the files to read in. If the stk module,
+ provided by CIAO, is available then CIAO stack syntax [2]_ can be
+ used. The order of the files does not matter, but it is currently
+ assumed that there is only one file per annulus.
+ rscale : AstroPy quantity
+ The scaling factor used to convert the XFLT radii (XFLT001 and
+ XFLT002 keywords) to an angle. If the values are in arcseconds
+ then ``rscale`` would be set to ``1 * u.arcsec``.
+ rinner : float, optional
+ The inner radius of the central annulus, in the same system as
+ the XFLT0001 and XFLT002 keyword values (this is a unitless
+ value).
+ angdist : None or AstroPy.Quantity, optional
+ The angular-diameter distance to the source. If not given then
+ it is calculated using the source redshift along with the
+ `cosmology` attribute.
+ cosmology : None or astropy.cosmology object, optional
+ The cosmology used to convert redshift to an angular-diameter
+ distance. This is used when `angdist` is None. If `cosmology`
+ is None then the `astropy.cosmology.Planck15` Cosmology
+ object is used.
+
+ Returns
+ -------
+ dep : Deproject instance
+ The deproject instance with the files loaded and associated
+ with the correct annuli.
+
+ Notes
+ -----
+ This currently is *not* guaranteed to support multiple data sets in
+ the same annulus. There is no check that the annuli are touching and
+ do not overlap.
+
+ References
+ ----------
+
+ .. [1] https://asd.gsfc.nasa.gov/XSPECwiki/projct_model
+
+ .. [2] http://cxc.harvard.edu/ciao/ahelp/stack.html
+
+ Examples
+ --------
+
+ Create a Deproject instance from the files matching the pattern
+ "src*.pi", whose XFLT radii are in ACIS pixels:
+
+ >>> dep = deproject_from_xflt('src*.pi', 0.492 * u.arcsec)
+
+ When used in CIAO, the stack syntax can be used to specify the
+ files, so if the file clus.stk contains the file names, one
+ per line, then the following will read them in and create a
+ Deproject instance. In this case the XFLT radii are in arcminutes:
+
+ >>> dep = deproject_from_xflt('@clus.stk', 1 * u.arcmin)
+
+ """
+
+ infiles = _find_xflt_files(pat)
+ print("# Found {} annuli".format(len(infiles)))
+
+ # TODO: check that annuli are touching (ie that outer[i] = inner[i+1])
+ xs = [rinner] + [x[1][0] for x in infiles]
+ radii = numpy.asarray(xs) * rscale
+ print("# Inner radius: {}".format(radii[0]))
+ print("# Outer radius: {}".format(radii[-1]))
+
+ # assume that theta_min < theta_max so this is valid
+ #
+ thetas = numpy.asarray([x[1][4] - x[1][3] for x in infiles])
+
+ tmin = thetas.min()
+ tmax = thetas.max()
+ if tmin <= 0:
+ raise ValueError("Found theta range of {}; ".format(tmin) +
+ "XFLT0004/5 keywords assumed to be in " +
+ "clockwise order")
+ if tmax > 360:
+ raise ValueError("Found theta range of {}; ".format(tmax) +
+ "XFLT0004/5 keywords assumed to be in " +
+ "clockwise order")
+
+ if tmin == tmax:
+ print("# Theta = {} degrees".format(tmin))
+ else:
+ print("# Theta ranges from {} to {} degrees".format(tmin, tmax))
+
+ print("#")
+
+ thetas = thetas * u.deg
+ dep = Deproject(radii,
+ theta=thetas,
+ angdist=angdist,
+ cosmology=cosmology)
+
+ for i, x in enumerate(infiles):
+ dep.load_pha(x[0], annulus=i)
+
+ return dep
diff --git a/src/deproject/fieldstore.py b/src/deproject/fieldstore.py
new file mode 100644
index 0000000..a9c227e
--- /dev/null
+++ b/src/deproject/fieldstore.py
@@ -0,0 +1,298 @@
+"""
+Convert fields from Sherpa structures into arrays.
+
+This is highly-specialised for the deproject code base. It is probably
+excessive, but the overall separation was useful in trying to capture
+the necessary functionality so it has been kept for now. That is, this
+is an internal module and should not be used by user code.
+
+Copyright 2018, Center for Astrophysics | Harvard & Smithsonian
+Douglas Burke (dburke@cfa.harvard.edu)
+"""
+
+from collections import defaultdict, OrderedDict
+
+import sherpa.astro.ui
+
+
+class FieldStore:
+ """Run and convert an object with keys to a set of values.
+
+ Parameters
+ ----------
+ runfunc : function
+ The function to call, which takes as argument the dataset ids.
+ getfunc : function
+ The function to call to return the structure. It takes no
+ arguments.
+ fields : sequence of str
+ The field names to report, and the order to use. These fields
+ must exist in the object returned by the getfunc.
+ """
+
+ def __init__(self, runfunc, getfunc, fields):
+
+ if not callable(runfunc):
+ raise ValueError("runfunc must be callable")
+
+ if not callable(getfunc):
+ raise ValueError("getfunc must be callable")
+
+ self._run = runfunc
+ self._get = getfunc
+
+ # Copy the fields and ensure it is a list; I'm not sure why I
+ # want a list, but a list we get.
+ #
+ self._fields = [f for f in fields]
+
+ def run(self, annuli, components, *ids):
+ """Run the function and return the values.
+
+ Parameters
+ ----------
+ annuli : dict
+ The key is the dataset identifier and the value the annulus
+ number.
+ components : dict
+ The keys are the annulus number and the values are a list of
+ tuples of (model components, parent model name).
+ *ids : dataset identifiers
+ The dataset identifiers to use. They must all exist as keys
+ in the `annuli` map.
+
+ Returns
+ -------
+ data : dict of OrderedDict values
+ The keys are the annulus value, and the values are the
+ shell results. The order of the keys in the results dict is
+ defined by the fields list given to the object constructor.
+
+ Notes
+ -----
+ There is only support for simple parameter links, where one parameter
+ is set equal to another, and both model components are associated
+ with the given set of dataset identifiers.
+
+ """
+
+ anns = sorted(list(set([annuli[i] for i in ids])))
+
+ datamap = defaultdict(list)
+ for i, ann in annuli.items():
+ if ann not in anns:
+ continue
+ datamap[ann].append(i)
+
+ self._run(*ids)
+ indata = self._get()
+
+ # Group the data by annulus. The assumption here is that each
+ # dataset in the same annulus has the same set of model components.
+ #
+ # All the fields are copied over to each annulus except for
+ # datasets
+ # parameter values
+ # which are restricted to the given annulus.
+ #
+ anndata = {}
+ for ann in anns:
+ dids = datamap[ann]
+
+ # Sanity check
+ for did in dids:
+ if did not in indata.datasets:
+ raise RuntimeError("Dataset {} not found in the ouput!".format(did))
+
+ # Create the non-parameter columns
+ #
+ store = OrderedDict()
+ store['datasets'] = dids
+ for field in self._fields:
+ if field == 'datasets':
+ continue
+
+ store[field] = getattr(indata, field)
+
+ for mcpt, pname in components[ann]:
+ for par in mcpt.pars:
+
+ # It would be nice to report these values too,
+ # for completeness, but leave for now.
+ #
+ # Note: linked parameters appear to be marked as
+ # frozen...
+ #
+ # I think a better check may be needed (or send in
+ # the actual pars to report).
+ #
+ if par.frozen and par.link is None:
+ continue
+
+ # It is important to use the "base" name for the
+ # parameter - e.g 'xsapec.kT' - rather than the
+ # name for this shell (e.g. 'xsapec_0.kT').
+ #
+ name = "{}.{}".format(pname, par.name)
+ self._extract(indata, par, name, store)
+
+ anndata[ann] = store
+
+ return anndata
+
+ def _extract(self, datavals, par, field, store):
+ """Extract the parameter value from the results structure.
+
+ Parameters
+ ----------
+ datavals
+ The results structure (FitResults, ErrorEstResults).
+ par : sherpa.models.parameter.Parameter instance
+ The model parameter of interest. It must be related to the
+ datasets in datavals.
+ field : str
+ The name of the item of interest.
+ store : dict
+ The results are added to this structure, using the field key
+ (multiple values can be added, in which case the field is
+ expected to be a suffix).
+
+ """
+
+ name = par.fullname
+
+ # Just want the value. Could just take the value from the
+ # parameter input, but use the value in datavals.
+ #
+ if self._store(datavals, name, field, store):
+ return
+
+ # Parameter not found, which should mean it is a link
+ #
+ if par.link is None:
+ raise ValueError("Parameter {} not found in the fit results".format(name))
+
+ # For now only support simple links - equality and not a more-complex
+ # functional relationship. Assume that checking whether modelname
+ # is not empty is a sufficient check
+ #
+ if par.link.modelname == '':
+ raise ValueError("Unable to handle complex parameter links: {}".format(name))
+
+ if self._store(datavals, par.link.fullname, field, store):
+ return
+
+ raise ValueError("Unable to find linked value for {}".format(name))
+
+ def _store(self, datavals, name, field, store):
+ """Add the values to the output structure.
+
+ This functionality is provided by the sub-classes.
+
+ Parameters
+ ----------
+ datavals
+ The results structure (FitResults, ErrorEstResults).
+ name : str
+ The parameter name.
+ field : str
+ The name of the item of interest.
+ store : dict
+ The results are added to this structure, using the field key
+ (multiple values can be added, in which case the field is
+ expected to be a suffix).
+
+ Returns
+ -------
+ flag : bool
+ True if the parameter was found and added to store, False
+ otherwise.
+
+ """
+
+ raise NotImplementedError()
+
+
+class FitStore(FieldStore):
+ """Run and extract the fields from a fit.
+
+ The constructor has no arguments.
+ """
+
+ def __init__(self):
+ super().__init__(sherpa.astro.ui.fit,
+ sherpa.astro.ui.get_fit_results,
+ ['datasets', 'succeeded',
+ 'statval', 'dstatval',
+ 'numpoints', 'dof', 'qval', 'rstat',
+ 'nfev'])
+
+ def _store(self, datavals, name, field, store):
+
+ # Assume that par.fullname can only occur once in parnames.
+ #
+ for pname, pval in zip(datavals.parnames, datavals.parvals):
+ if pname == name:
+ store[field] = pval
+ return True
+
+ return False
+
+
+class ErrorStore(FieldStore):
+ """Run and extract the fields from a sherpa.fit.ErrorEstResults instance.
+
+ Parameters
+ ----------
+ runfunc : function
+ The function to call, which takes as argument the dataset ids.
+ getfunc : function
+ The function to call to return the structure. It takes no
+ arguments.
+
+ """
+
+ def __init__(self, runfunc, getfunc):
+ super().__init__(runfunc, getfunc,
+ ['datasets', 'sigma', 'percent', 'nfits'])
+
+ def _store(self, datavals, name, field, store):
+
+ # Assume that par.fullname can only occur once in parnames.
+ #
+ # QUS: do we need to handle limits or cases where the limit
+ # could not be accurately determined.
+ #
+ for pname, pval, plo, phi in zip(datavals.parnames,
+ datavals.parvals,
+ datavals.parmins,
+ datavals.parmaxes):
+ if pname == name:
+ store[field] = pval
+ store[field + "_lo"] = plo
+ store[field + "_hi"] = phi
+ return True
+
+ return False
+
+
+class CovarStore(ErrorStore):
+ """Run and extract the covariance results.
+
+ The constructor takes no arguments.
+ """
+
+ def __init__(self):
+ super().__init__(sherpa.astro.ui.covar,
+ sherpa.astro.ui.get_covar_results)
+
+
+class ConfStore(ErrorStore):
+ """Run and extract the confidence results.
+
+ The constructor takes no arguments.
+ """
+
+ def __init__(self):
+ super().__init__(sherpa.astro.ui.conf,
+ sherpa.astro.ui.get_conf_results)
diff --git a/src/deproject/simplegraph.py b/src/deproject/simplegraph.py
new file mode 100644
index 0000000..08c6db2
--- /dev/null
+++ b/src/deproject/simplegraph.py
@@ -0,0 +1,100 @@
+"""
+A simple/basic graph structure designed to help group together annuli.
+
+:Copyright: Smithsonian Astrophysical Observatory (2018)
+:Author: Douglas Burke (dburke@cfa.harvard.edu)
+"""
+
+__all__ = ("SimpleGraph", )
+
+
+class SimpleGraph:
+ """A graph where links are bi-directional
+
+ The identifiers must be hashable (as they are used in sets).
+ This is not optimised, and so is intended only for small graphs.
+ It is expected that a link is manually created for each identifier,
+ that is::
+
+ >>> graph = SimpleGraph()
+ >>> for n in ['n1', 'n2', 'n3, 'n4', 'n5']:
+ ... graph.add_link(n, n)
+ ...
+ >>> graph.add_link('n2', 'n3')
+ >>> graph.add_link('n3', 'n4')
+ >>> graph.get_groups()
+ [['n1'], ['n2', 'n3', 'n4'], ['n5']]
+
+ """
+
+ def __init__(self):
+ self.groups = []
+
+ def add_link(self, id1, id2):
+ """Adds a link between two identifiers (which can be the same).
+
+ The order between id1 and id2 is assumed to be bi-directional.
+
+ Parameters
+ ----------
+ id1, id2
+ identifier of a pair-wise link
+
+ """
+
+ g1 = None
+ g2 = None
+ for i, grp in enumerate(self.groups):
+ if id1 in grp:
+ assert g1 is None
+ g1 = i
+
+ if id2 in grp:
+ assert g2 is None
+ g2 = i
+
+ # Not known, so add in. This is not expected to happen,
+ # since the idea is to explicitly call with id1 == id2,
+ # but support it.
+ #
+ if g1 is None and g2 is None:
+ self.groups.append(set([id1, id2]))
+ return
+
+ # Already known about
+ if g1 == g2:
+ return
+
+ # It is expected that g1 and g2 are both set, but this
+ # is not guaranteed.
+ #
+ if g2 is None:
+ self.groups[g1].add(id2)
+ return
+
+ if g1 is None:
+ self.groups[g2].add(id1)
+ return
+
+ # id1 and id2 are part of the g1 and g2 groups, by construction,
+ # so they do not need to be added.
+ #
+ newgrp = self.groups[g1].union(self.groups[g2])
+ self.groups[g1] = newgrp
+ del self.groups[g2]
+
+ def get_groups(self):
+ """Return a copy of the groups, that is linked identifiers.
+
+ Returns
+ -------
+ groups : list of list
+ Each group is identified by the nodes in the list. There is
+ no guarantee of ordering of either list.
+ """
+
+ out = []
+ for grp in self.groups:
+ out.append(list(grp))
+
+ return out
diff --git a/src/deproject/specstack.py b/src/deproject/specstack.py
new file mode 100644
index 0000000..d09b2f3
--- /dev/null
+++ b/src/deproject/specstack.py
@@ -0,0 +1,855 @@
+"""
+Manipulate a stack of spectra in Sherpa.
+
+The methods in the SpecStack class provide a way to automatically apply
+familiar Sherpa commands such as `set_par`_ or `freeze`_ or `plot_fit`_
+to a stack of PHA spectra. This simplifies simultaneous fitting of
+multiple spectra.
+
+Note that the :mod:`specstack` module is currently distributed within with the
+:mod:`deproject` package. `Specstack` is not yet fully documented or tested
+outside the context of `deproject`.
+
+:Copyright: Smithsonian Astrophysical Observatory (2009, 2019)
+:Author: Tom Aldcroft (taldcroft@cfa.harvard.edu), Douglas Burke (dburke@cfa.harvard.edu)
+"""
+
+import re
+import numpy
+from sherpa.astro import ui
+
+try:
+ from sherpa.astro.plot import backend
+ backend_name = backend.name
+except ImportError:
+ backend_name = 'none'
+
+# Trying to move away from direct access to I/O and plotting code.
+#
+if backend_name == 'pychips':
+ try:
+ import pychips
+ except ImportError:
+ pass
+
+elif backend_name == 'pylab':
+ try:
+ import matplotlib.pyplot as plt
+ except ImportError:
+ pass
+
+
+def _sherpa_plot_func(func):
+ "Wrap functions which are called once per dataset per annulus"
+ def _sherpa_plot_func_(self, *args, **kwargs):
+ self._sherpa_plot(func, *args, **kwargs)
+ return _sherpa_plot_func_
+
+
+def _sherpa_plot_func_single(func):
+ "Wrap functions which are called once per annulus"
+ def _sherpa_plot_func_(self, *args, **kwargs):
+ kwargs['single_call_per_shell'] = True
+ self._sherpa_plot(func, *args, **kwargs)
+ return _sherpa_plot_func_
+
+
+class SpecStack:
+ """Manipulate a stack of spectra in Sherpa.
+
+ This `SpecStack` class provides a number of wrappers around Sherpa
+ routines that handle loading data, setting the source model,
+ setting up the data to fit, such as: the noticed energy range, how to
+ handle the background, extracting parameter values, and plotting
+ data.
+
+ """
+
+ datasets = None
+ """Information (dataset identifier, annulus, name) about the loaded data."""
+
+ def __init__(self):
+ self.datasets = []
+ self.obsids = set()
+ self._reset_model_comps()
+
+ def _reset_model_comps(self):
+ """Clean out the model component information.
+
+ Notes
+ -----
+ This does not delete the previous model components from
+ Sherpa, but does remove them from the object.
+ """
+
+ # Generic model components in source model expression
+ self.srcmodel_comps = []
+
+ # All instantiated model components for shells
+ self.model_comps = []
+
+ def load_pha(self, specfile, annulus):
+ """Load a pha file and add to the datasets for stacked analysis.
+
+ It is required that datasets for all annuli are loaded before
+ the source model is created (to ensure that components are
+ created for each annulus).
+
+ Parameters
+ ----------
+ specfile : str or sherpa.astro.data.DataPHA object
+ If a string, the name of the file containing the source spectrum,
+ which must be in PHA format (the data is expected to be extracted
+ on the PI column). If a DataPHA object, then this is used (and
+ is assumed to contain any needed background data).
+ annulus : int
+ The annulus number for the data.
+
+ Returns
+ -------
+ dataid : int
+ The Sherpa dataset identifier used for this spectrum.
+
+ Examples
+ --------
+
+ Load the data for four annuli from the files 'ann1.pi' to 'ann4.pi'.
+
+ >>> dep.load_pha('ann1.pi', 0)
+ >>> dep.load_pha('ann2.pi', 1)
+ >>> dep.load_pha('ann3.pi', 2)
+ >>> dep.load_pha('ann4.pi', 3)
+
+ Load in the PHA files into Sherpa DataPHA objects, and then use
+ these objects:
+
+ >>> s1 = ui.unpack_pha('src1.pi')
+ >>> s2 = ui.unpack_pha('src2.pi')
+ >>> s3 = ui.unpack_pha('src3.pi')
+ >>> dep.load_pha(s1, 0)
+ >>> dep.load_pha(s2, 1)
+ >>> dep.load_pha(s3, 2)
+
+ """
+
+ dataid = len(self.datasets)
+
+ # If the input has a counts attribute then assume a DataPHA
+ # style object.
+ #
+ if hasattr(specfile, 'counts'):
+ print('Using spectrum {} '.format(specfile.name) +
+ ' as dataset id {}'.format(dataid))
+ ui.set_data(dataid, specfile)
+
+ else:
+ print('Loading spectrum file {} '.format(specfile) +
+ ' as dataset id {}'.format(dataid))
+ ui.load_pha(dataid, specfile)
+
+ data = ui.get_data(dataid)
+ try:
+ obsid = int(data.header['OBS_ID'])
+ except (KeyError, TypeError, ValueError):
+ obsid = 0
+
+ dataset = dict(file=specfile,
+ obsid=obsid,
+ id=dataid,
+ annulus=annulus
+ )
+ self.datasets.append(dataset)
+ self.obsids.add(obsid)
+ return dataid
+
+ def _split_parname(self, par):
+ """Return the model type and parameter name.
+
+ The input value can either be a "generic" specifier, such as
+ 'xsphabs.nh' or specific to an annulus ('xspahbs_2.nh')
+
+ Parameters
+ ----------
+ par : str
+ The parameter name (e.g. 'xsapec.kt').
+
+ Returns
+ -------
+ rtype : str, str
+ The model and parameter name.
+
+ Examples
+ --------
+
+ The foillowing will set `mname` to "xsapec" and `pname` to "kt":
+
+ >>> mname, pname = dep._split_parname('xsapec.kt')
+
+ """
+
+ pos = par.rfind('.')
+ if pos == -1:
+ raise ValueError("Parameter name must match 'model_type.par_name'")
+ return par[:pos], par[pos + 1:]
+
+ def find_parval(self, parname):
+ """Return the value of the first parameter matching the name.
+
+ Parameters
+ ----------
+ parname : str
+ The parameter name. The case is ignored in the match, and the
+ first match is returned.
+
+ Returns
+ -------
+ parval : float
+ The parameter value
+
+ Raises
+ ------
+ ValueError
+ There is no match for the parameter.
+
+ See Also
+ --------
+ find_norm, set_par
+
+ Examples
+ --------
+
+ >>> kt = dep.find_parval('kt')
+
+ """
+ RE_parname = re.compile(parname + '$', re.IGNORECASE)
+ for model_comp in self.model_comps:
+ mc = model_comp['object'] # Get the model component object
+ for par in mc.pars:
+ if RE_parname.match(par.name):
+ return par.val
+
+ raise ValueError('Parameter %s not found in any model component' % parname)
+
+ def find_norm(self, shell):
+ """Return the normalization value for the given shell.
+
+ This is limited to XSPEC-style models, where the parameter is called
+ "norm".
+
+ Parameters
+ ----------
+ shell : int
+ The shell number.
+
+ Returns
+ -------
+ norm : float
+ The normalization of the shell.
+
+ Raises
+ ------
+ ValueError
+ If there is not one `norm` parameter for the shell.
+
+ See Also
+ --------
+ find_parval, set_par
+
+ """
+ norms = []
+ for model_comp in self.model_comps:
+ if model_comp['shell'] == shell:
+ mc = model_comp['object'] # Get the model component object
+ if hasattr(mc, 'norm'): # Special attr of model component
+ norms.append(mc.norm.val)
+
+ if len(norms) == 0:
+ raise ValueError('Model for shell %d has no norm' % shell)
+ elif len(norms) > 1:
+ raise ValueError('Model for shell %d has multiple norms' % shell)
+
+ return norms[0]
+
+ def _get_n_datasets(self):
+ """How many datasets are registered?
+
+ This is not the same as the number of annuli.
+
+ Returns
+ -------
+ ndata : int
+ The number of datasets.
+ """
+ return len(self.datasets)
+
+ n_datasets = property(_get_n_datasets)
+
+ def _sherpa_cmd(self, func, *args):
+ """Apply a function to each dataset.
+
+ Parameters
+ ----------
+ func : function reference
+ This function is called with each dataset identifier as
+ the first argument, with the remaining arguments next.
+ *args
+ The additional arguments for the function.
+
+ """
+ for dataset in self.datasets:
+ func(dataset['id'], *args)
+
+ def subtract(self, *args):
+ """Subtract the background from each dataset.
+
+ See Also
+ --------
+ unsubtract
+
+ Examples
+ --------
+
+ >>> dep.substract()
+
+ """
+ self._sherpa_cmd(ui.subtract, *args)
+
+ def unsubtract(self, *args):
+ """Un-subtract the background from each dataset.
+
+ This can be useful when you want to compare the results to
+ the "wstat" stat (a Poisson-based stat which includes the
+ background data as a component and provides a goodness-of-fit
+ estimate).
+
+ See Also
+ --------
+ subtract
+
+ Examples
+ --------
+
+ >>> dep.unsubstract()
+
+ """
+ self._sherpa_cmd(ui.unsubtract, *args)
+
+ def group(self, *args):
+ """Apply the grouping for each data set.
+
+ See Also
+ --------
+ ungroup
+
+ Examples
+ --------
+
+ >>> dep.group()
+
+ """
+ self._sherpa_cmd(ui.group, *args)
+
+ def ungroup(self, *args):
+ """Turn off the grouping for each data set.
+
+ See Also
+ --------
+ group
+
+ Examples
+ --------
+
+ >>> dep.ungroup()
+
+ """
+ self._sherpa_cmd(ui.ungroup, *args)
+
+ def notice(self, *args):
+ """Apply Sherpa notice command to each dataset.
+
+ The filter is applied to each data set separately.
+
+ See Also
+ --------
+ ignore
+
+ Examples
+ --------
+
+ Restrict the analysis to those bins which fall in the range
+ 0.5 to 7.0 keV, where the limits are included in the noticed
+ range. The first call to `notice` is used to clear any
+ existing filter.
+
+ >>> dep.notice(None, None)
+ >>> dep.notice(0.5, 7.0)
+
+ """
+ self._sherpa_cmd(ui.notice_id, *args)
+
+ def ignore(self, *args):
+ """Apply Sherpa ignore command to each dataset.
+
+ The filter is applied to each data set separately.
+
+ See Also
+ --------
+ notice
+
+ Examples
+ --------
+
+ Restrict the analysis to those bins which fall in the range
+ 0.5 to 7.0 keV, where the limits are not included in the
+ noticed range. The call to `notice` is used to clear any
+ existing filter.
+
+ >>> dep.notice(None, None)
+ >>> dep.ignore(None, 0.5)
+ >>> dep.ignore(7.0, None)
+
+ """
+ self._sherpa_cmd(ui.ignore_id, *args)
+
+ def _sherpa_par(self, func, par, msg, *args):
+ """Apply the function to the given parameter.
+
+ See thaw(), freeze(), set_par() and get_par() for examples.
+
+ Parameters
+ ----------
+ func : function reference
+ The function to call. The first argument is the name of the
+ given parameter (including the model name), and then the
+ optional arguments (`args`).
+ par : str
+ The parameter name to apply the function to (per shell). It
+ must contain the model name (e.g. `xsapec.kt`).
+ msg : str or None
+ If not None, a format string that is printed each time `func`
+ is called. It must contain one "%s" token, which is sent the
+ full name of the parameter being processed.
+ *args
+ Additional arguments for `func`.
+
+ Returns
+ -------
+ rvals : ndarray
+ The return value from `func` for each parameter, ordered by
+ shell.
+
+ Raises
+ ------
+ ValueError
+ If the parameter is not found.
+
+ """
+
+ model_type, parname = self._split_parname(par)
+
+ vals = [] # return values
+ for model_comp in self.model_comps:
+ if model_comp['type'] == model_type:
+ fullparname = '%s.%s' % (model_comp['name'], parname)
+ if msg is not None:
+ print(msg % fullparname)
+ vals.append(func(fullparname, *args))
+
+ if len(vals) == 0:
+ raise ValueError("No parameter found matching {}".format(par))
+
+ return vals
+
+ def thaw(self, par):
+ """Thaw the given parameter in each shell.
+
+ Parameters
+ ----------
+ par : str
+ The parameter name, specified as .
+
+ See Also
+ --------
+ freeze, tie_par, untie_par
+
+ Examples
+ --------
+
+ >>> dep.thaw('clus.abundanc')
+
+ """
+ self._sherpa_par(ui.thaw, par, 'Thawing %s')
+
+ def freeze(self, par):
+ """Freeze the given parameter in each shell.
+
+ Parameters
+ ----------
+ par : str
+ The parameter name, specified as .
+
+ See Also
+ --------
+ thaw, tie_par, untie_par
+
+ Examples
+ --------
+
+ >>> dep.freeze('clus.abundanc')
+
+ """
+ self._sherpa_par(ui.freeze, par, 'Freezing %s')
+
+ def set_par(self, par, val):
+ """Set the parameter value in each shell.
+
+ Parameters
+ ----------
+ par : str
+ The parameter name, specified as .
+ val : float
+ The parameter value.
+
+ See Also
+ --------
+ get_par, tie_par
+
+ Examples
+ --------
+
+ >>> dep.set_par('xsapec.abundanc', 0.25)
+
+ """
+ self._sherpa_par(ui.set_par, par, 'Setting %%s=%s' % str(val), val)
+
+ def get_par(self, par):
+ """Return the parameter value for each shell.
+
+ Parameters
+ ----------
+ par : str
+ The parameter name, specified as .
+
+ Returns
+ -------
+ vals : ndarray
+ The parameter values, in shell order.
+
+ See Also
+ --------
+ find_parval, find_norm, set_par
+
+ Examples
+ --------
+
+ >>> kts = dep.get_par('xsapec.kt')
+
+ """
+
+ pars = self._sherpa_par(ui.get_par, par, 'Getting %s')
+ return numpy.array([x.val for x in pars])
+
+ def _find_shell_par(self, par, shell):
+ """Return the shell parameter.
+
+ Parameters
+ ----------
+ par : str
+ The parameter name, in the form `model name.par name`.
+ shell : int
+ The shell number
+
+ Returns
+ -------
+ par : sherpa.models.parameter.Parameter
+
+ """
+
+ model_type, parname = self._split_parname(par)
+ models = []
+ for model_comp in self.model_comps:
+ if model_comp['shell'] != shell or \
+ model_comp['type'] != model_type:
+ continue
+
+ models.append(model_comp['object'])
+
+ if len(models) == 0:
+ raise ValueError('No parameter found matching ' +
+ '{} for shell {}'.format(par, shell))
+ elif len(models) > 1:
+ raise ValueError('Multiple parameters found matching ' +
+ '{} for shell {}'.format(par, shell))
+
+ return getattr(models[0], parname)
+
+ def tie_par(self, par, base, *others):
+ """Tie parameters in one or more shells to the base shell.
+
+ This is a limited form of the Sherpa ability to link parameters,
+ since it sets the parameter in the other shells to the same
+ value as the parameter in the base shell. More complex
+ situations will require direct calls to `sherpa.astro.ui.link`.
+
+ Parameters
+ ----------
+ par : str
+ The parameter specifier, as ..
+ base : int
+ The base shell number.
+ *others : scalar
+ The shell, or shells to link to the base shell.
+
+ See Also
+ --------
+ set_par, untie_par
+
+ Examples
+ --------
+
+ Tie the temperature and abundance parameters in shell 9 to that
+ in shell 8, so that any fits will set the shell 9 values to those
+ used in shell 8 (so reducing the number of free parameters in the
+ fit).
+
+ >>> dep.tie_par('xsapec.kt', 8, 9)
+ Tying xsapec_9.kT to xsapec_8.kT
+ >>> dep.tie_par('xsapec.abundanc', 8, 9)
+ Tying xsapec_9.Abundanc to xsapec_8.Abundanc
+
+ Tie three annuli together:
+
+ >>> dep.tie_par('xsapec.kt', 12, 13, 14)
+ Tying xsapec_13.kT to xsapec_12.kT
+ Tying xsapec_14.kT to xsapec_12.kT
+
+ """
+
+ bpar = self._find_shell_par(par, base)
+
+ for other in others:
+ opar = self._find_shell_par(par, other)
+ print('Tying {} to {}'.format(opar.fullname, bpar.fullname))
+ ui.link(opar, bpar)
+
+ def untie_par(self, par, *others):
+ """Remove the parameter tie/link in the shell.
+
+ This is intended to remove links between shells created by `tie_par`,
+ but will remove any links created by `sherpa.astro.ui.link`.
+
+ Parameters
+ ----------
+ par : str
+ The parameter specifier, as ..
+ *others : scalar
+ The shell, or shells to un-tie/unlink.
+
+ See Also
+ --------
+ tie_par
+
+ Notes
+ -----
+ It is safe to call on a parameter that is not tied or linked
+ to another parameter.
+
+ Examples
+ --------
+
+ Untie the abundance parameter in shell 9; that is, it is now free
+ to vary independently in a fit.
+
+ >>> dep.untie_par('xsapec.abundanc', 9)
+ Untying xsapec_9.Abundanc
+
+ Untie multiple annuli:
+
+ >>> dep.untie_par('xsmekal.kt', 13, 14)
+ Untying xsmekal_13.kT
+ Untying xsmekal_14.kT
+
+ """
+
+ for other in others:
+ opar = self._find_shell_par(par, other)
+ if opar.link is None:
+ continue
+
+ print('Untying {}'.format(opar.fullname))
+ ui.unlink(opar)
+
+ def _sherpa_plot(self, func, *args, **kwargs):
+ """Call the Sherpa plot ``func`` for each shell.
+
+ Parameters
+ ----------
+ func : function reference
+ The Sherpa plot function
+ args
+ The arguments for `func`
+ kwargs
+ Any keyword arguments for `func`. There are special
+ keywords which are not passed on: `single_call_per_shell` is
+ a boolean which indicates that the function is only
+ called once per shell, and `add_shell_value` which indicates
+ that the first argument is formatted to accept the shell
+ value.
+
+ Notes
+ -----
+ This method attempts to handle the differences when using
+ ChIPS or Matplotlib as the Sherpa plotting backend, but has
+ not been properly tested, so there may be issues.
+
+ It is known not to work when the input function is plot_fit_delchi
+ or plot_fit_resid, at least with the Matplotlib backend. It is
+ unclear whether this is a problem here or the Sherpa
+ matplotlib backend.
+
+ """
+
+ # Attempt to support multiple datasets per annulus.
+ #
+ dmap = [[] for _ in range(self.nshell)]
+ for d in self.datasets:
+ dmap[d['annulus']].append(d['id'])
+
+ # Extract the arguments used here
+ #
+ special = {}
+ for key in ['add_shell_value', 'single_call_per_shell']:
+ val = False
+ if key in kwargs:
+ val = kwargs[key]
+ del kwargs[key]
+
+ special[key] = val
+
+ nargs = len(args)
+
+ for shell in range(self.nshell):
+ if backend_name == 'pychips':
+ window_id = 'Shell%d' % shell
+ try:
+ pychips.add_window(['id', window_id])
+ except RuntimeError:
+ pychips.set_current_window(window_id)
+
+ elif backend_name == 'pylab':
+ plt.figure(shell)
+
+ new_args = args
+ if nargs > 0:
+ if special['add_shell_value']:
+ new_args = list(args)[:]
+ new_args[0] = new_args[0].format(shell)
+
+ # For the moment assume that if the user supplied an
+ # argument then we use it, whatever the various
+ # settings are. This catches errors like
+ # dep.plot_fit('rstat')
+ #
+ if special['single_call_per_shell'] or nargs > 0:
+ func(*new_args, **kwargs)
+ else:
+ # call it once per data set, assuming that the
+ # overplot keyword is supported by this function
+ # (could check this condition but leave that for now)
+ #
+ overplot = False
+ for did in dmap[shell]:
+ kwargs['overplot'] = overplot
+ func(did, **kwargs)
+ overplot = True
+
+ def print_window(self, *args, **kwargs):
+ """Create a hardcopy version of each plot window.
+
+ Parameters
+ ----------
+ args
+ The arguments to be sent to the "create a hardcopy" routine
+ (print_window for ChIPS and savefig for Matplotlib).
+ The first argument, if given, is assumed to be the file name
+ and so will have the shell number added to it.
+ kwargs
+ Keyword arguments for the call.
+
+ Notes
+ -----
+ This is not guaranteed to work properly for Matplotlib.
+
+ Examples
+ --------
+
+ Create hardcopy versions of the data plots, called "data0",
+ "data1", ...
+
+ >>> dep.plot_data()
+ >>> dep.print_window('data')
+
+ """
+
+ if len(args) > 0:
+ args = list(args)[:]
+ try:
+ args[0] += '{}'
+ except TypeError:
+ raise TypeError("First argument must be a string")
+
+ kwargs['add_shell_value'] = True
+
+ kwargs['single_call_per_shell'] = True
+ if backend_name == 'pychips':
+ plotfn = pychips.plot_window
+ elif backend_name == 'pylab':
+ plotfn = plt.savefig
+ else:
+ raise RuntimeError("Unrecognized plotting backend: {}".format(backend_name))
+
+ self._sherpa_plot(plotfn, *args, **kwargs)
+
+ def dummyfunc(self, *args, **kwargs):
+ pass
+
+ try:
+ if backend_name == 'pychips':
+ log_scale = _sherpa_plot_func_single(pychips.log_scale)
+ linear_scale = _sherpa_plot_func_single(pychips.linear_scale)
+ elif backend_name == 'pylab':
+ # should be able to handle this
+ log_scale = dummyfunc
+ linear_scale = dummyfunc
+
+ except NameError:
+ # Allow doc generation to work without pychips
+ log_scale = dummyfunc
+ linear_scale = dummyfunc
+
+ plot_fit = _sherpa_plot_func(ui.plot_fit)
+ plot_arf = _sherpa_plot_func(ui.plot_arf)
+ plot_bkg_fit = _sherpa_plot_func(ui.plot_bkg_fit)
+ plot_bkg_ratio = _sherpa_plot_func(ui.plot_bkg_ratio)
+ plot_chisqr = _sherpa_plot_func(ui.plot_chisqr)
+ plot_fit_delchi = _sherpa_plot_func(ui.plot_fit_delchi)
+ plot_psf = _sherpa_plot_func(ui.plot_psf)
+ plot_bkg = _sherpa_plot_func(ui.plot_bkg)
+ plot_bkg_fit_delchi = _sherpa_plot_func(ui.plot_bkg_fit_delchi)
+ plot_bkg_resid = _sherpa_plot_func(ui.plot_bkg_resid)
+ plot_data = _sherpa_plot_func(ui.plot_data)
+ plot_fit_resid = _sherpa_plot_func(ui.plot_fit_resid)
+ plot_ratio = _sherpa_plot_func(ui.plot_ratio)
+ plot_bkg_chisqr = _sherpa_plot_func(ui.plot_bkg_chisqr)
+ plot_bkg_fit_resid = _sherpa_plot_func(ui.plot_bkg_fit_resid)
+ plot_bkg_source = _sherpa_plot_func(ui.plot_bkg_source)
+ plot_delchi = _sherpa_plot_func(ui.plot_delchi)
+ plot_model = _sherpa_plot_func(ui.plot_model)
+ plot_resid = _sherpa_plot_func(ui.plot_resid)
+ plot_bkg_delchi = _sherpa_plot_func(ui.plot_bkg_delchi)
+ plot_bkg_model = _sherpa_plot_func(ui.plot_bkg_model)
+ plot_bkg_unconvolved = _sherpa_plot_func(ui.plot_bkg_source)
+ plot_bkg_source = _sherpa_plot_func(ui.plot_bkg_source)
+ plot_fit = _sherpa_plot_func(ui.plot_fit)
+ plot_order = _sherpa_plot_func(ui.plot_order)
+ plot_source = _sherpa_plot_func(ui.plot_source)
diff --git a/tests/test_deproject.py b/tests/test_deproject.py
new file mode 100644
index 0000000..13069ed
--- /dev/null
+++ b/tests/test_deproject.py
@@ -0,0 +1,978 @@
+"""
+Limited testing of deproject.deproject
+
+:Copyright: Smithsonian Astrophysical Observatory (2019)
+:Author: Douglas Burke (dburke@cfa.harvard.edu)
+"""
+
+# requies pytest 3.1 or greater for match argument to pytest.raises
+import pytest
+
+import numpy as np
+
+from deproject.deproject import Deproject
+
+from sherpa.astro import ui
+from sherpa.astro.data import DataPHA
+
+from astropy import units as u
+
+# We need models, but do not want to depend on XSPEC availability (since
+# we do not actually evaluate the models), so create "fake" versions
+# of several XSPEC models. I was going to make this conditional on
+# whether XSPEC is available, but it seems that over-writing things
+# doesn't cause problems for us, so avoid that complexity.
+#
+# I use the actual parameter settings from Sherpa, but this isn't
+# really needed here (could get away with only setting up a few
+# parameters).
+#
+# from sherpa.astro.xspec import XSphabs, XSwabs, XSapec, XSmekal
+
+from sherpa.models import Parameter, RegriddableModel1D
+from sherpa.models.parameter import hugeval
+
+
+class FakeXSPECModel(RegriddableModel1D):
+
+ def calc(self, *args, **kwargs):
+ "Always return 1 for each bin"
+
+ if len(args) == 0:
+ raise ValueError("No arguments")
+
+ # When called with a single array (so XSPEC style) we normally
+ # fill in the last value with 0 since it's invalid, but that
+ # isn't needed here.
+ #
+ nbins = len(args[0])
+ return np.ones(nbins)
+
+
+class XSphabs(FakeXSPECModel):
+
+ def __init__(self, name='phabs'):
+ self.nH = Parameter(name, 'nH', 1., 0.0, 1.e5, 0.0, hugeval, '10^22 atoms / cm^2')
+ FakeXSPECModel.__init__(self, name, (self.nH,))
+
+
+class XSwabs(FakeXSPECModel):
+
+ def __init__(self, name='wabs'):
+ self.nH = Parameter(name, 'nH', 1., 0.0, 1.e5, 0.0, hugeval, '10^22 atoms / cm^2')
+ FakeXSPECModel.__init__(self, name, (self.nH,))
+
+
+class XSapec(FakeXSPECModel):
+
+ def __init__(self, name='apec'):
+ self.kT = Parameter(name, 'kT', 1., 0.008, 64.0, 0.0, hugeval, 'keV')
+ self.Abundanc = Parameter(name, 'Abundanc', 1., 0., 5., 0.0, hugeval, frozen=True)
+ self.redshift = Parameter(name, 'redshift', 0., -0.999, 10., -0.999, hugeval, frozen=True)
+ self.norm = Parameter(name, 'norm', 1.0, 0.0, 1.0e24, 0.0, hugeval)
+ FakeXSPECModel.__init__(self, name, (self.kT, self.Abundanc, self.redshift, self.norm))
+
+
+class XSvapec(FakeXSPECModel):
+
+ def __init__(self, name='vapec'):
+ self.kT = Parameter(name, 'kT', 6.5, 0.0808, 68.447, 0.0, hugeval, 'keV')
+ self.He = Parameter(name, 'He', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.C = Parameter(name, 'C', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.N = Parameter(name, 'N', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.O = Parameter(name, 'O', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ne = Parameter(name, 'Ne', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Mg = Parameter(name, 'Mg', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Al = Parameter(name, 'Al', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Si = Parameter(name, 'Si', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.S = Parameter(name, 'S', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ar = Parameter(name, 'Ar', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ca = Parameter(name, 'Ca', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Fe = Parameter(name, 'Fe', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ni = Parameter(name, 'Ni', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.redshift = Parameter(name, 'redshift', 0., -0.999, 10., -0.999, hugeval, frozen=True)
+ self.norm = Parameter(name, 'norm', 1.0, 0.0, 1.0e24, 0.0, hugeval)
+ FakeXSPECModel.__init__(self, name, (self.kT, self.He, self.C, self.N, self.O, self.Ne, self.Mg, self.Al, self.Si, self.S, self.Ar, self.Ca, self.Fe, self.Ni, self.redshift, self.norm))
+
+
+class XSvvapec(FakeXSPECModel):
+
+ def __init__(self, name='vvapec'):
+ self.kT = Parameter(name, 'kT', 6.5, 0.0808, 68.447, 0.0, hugeval, 'keV')
+ self.H = Parameter(name, 'H', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.He = Parameter(name, 'He', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Li = Parameter(name, 'Li', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Be = Parameter(name, 'Be', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.B = Parameter(name, 'B', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.C = Parameter(name, 'C', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.N = Parameter(name, 'N', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.O = Parameter(name, 'O', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.F = Parameter(name, 'F', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ne = Parameter(name, 'Ne', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Na = Parameter(name, 'Na', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Mg = Parameter(name, 'Mg', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Al = Parameter(name, 'Al', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Si = Parameter(name, 'Si', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.P = Parameter(name, 'P', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.S = Parameter(name, 'S', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Cl = Parameter(name, 'Cl', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ar = Parameter(name, 'Ar', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.K = Parameter(name, 'K', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ca = Parameter(name, 'Ca', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Sc = Parameter(name, 'Sc', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ti = Parameter(name, 'Ti', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.V = Parameter(name, 'V', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Cr = Parameter(name, 'Cr', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Mn = Parameter(name, 'Mn', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Fe = Parameter(name, 'Fe', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Co = Parameter(name, 'Co', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Ni = Parameter(name, 'Ni', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Cu = Parameter(name, 'Cu', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Zn = Parameter(name, 'Zn', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.Redshift = Parameter(name, 'Redshift', 0., -0.999, 10., -hugeval, hugeval, frozen=True)
+ self.norm = Parameter(name, 'norm', 1.0, 0.0, 1.0e24, 0.0, hugeval)
+ FakeXSPECModel.__init__(self, name, (self.kT, self.H, self.He, self.Li, self.Be, self.B, self.C, self.N, self.O, self.F, self.Ne, self.Na, self.Mg, self.Al, self.Si, self.P, self.S, self.Cl, self.Ar, self.K, self.Ca, self.Sc, self.Ti, self.V, self.Cr, self.Mn, self.Fe, self.Co, self.Ni, self.Cu, self.Zn, self.Redshift, self.norm))
+
+
+class XSmekal(FakeXSPECModel):
+
+ def __init__(self, name='mekal'):
+ self.kT = Parameter(name, 'kT', 1., 0.0808, 79.9, 0.0, hugeval, 'keV')
+ self.nH = Parameter(name, 'nH', 1., 1.e-5, 1.e19, 0.0, hugeval, 'cm-3', True)
+ self.Abundanc = Parameter(name, 'Abundanc', 1., 0., 1000., 0.0, hugeval, frozen=True)
+ self.redshift = Parameter(name, 'redshift', 0., -0.999, 10., -0.999, hugeval, frozen=True)
+ self.switch = Parameter(name, 'switch', 1, 0, 2, 0, 2, alwaysfrozen=True)
+ self.norm = Parameter(name, 'norm', 1.0, 0.0, 1.0e24, 0.0, hugeval)
+
+ FakeXSPECModel.__init__(self, name, (self.kT, self.nH, self.Abundanc, self.redshift, self.switch, self.norm))
+
+
+# Ensure the "fake" models are registered with the Sherpa session.
+# I was expecting this to complain about over-writing existing classes -
+# that is, if XSPEC support is available - but it doesn't seem to.
+#
+ui.add_model(XSphabs)
+ui.add_model(XSwabs)
+ui.add_model(XSapec)
+ui.add_model(XSvapec)
+ui.add_model(XSvvapec)
+ui.add_model(XSmekal)
+
+
+@pytest.fixture(autouse=True)
+def cleanup_sherpa():
+ """Ensure that the Sherpa state is cleaned before and after each run."""
+
+ ui.clean()
+ yield
+ ui.clean()
+
+
+def _test_data(annulus=0):
+ """Return test data sets (very basic).
+
+ At present there are no tests that make use of the data,
+ so we could just use the same data set for each annulus.
+ """
+
+ n = 5 + annulus
+ chans = np.arange(n)
+ counts = np.repeat(1, n)
+
+ return DataPHA('ann{}'.format(annulus), chans, counts)
+
+
+@pytest.mark.parametrize("rad", [None, "radii", [], [10]])
+def test_fail_radii_type(rad):
+ """Need units for radii"""
+
+ estr = "Argument 'radii' to function '__init__' " + \
+ "has no 'unit' attribute. You may want to pass " + \
+ "in an astropy Quantity instead."
+ with pytest.raises(TypeError, match=estr):
+ Deproject(rad)
+
+
+# Pick a length since this is likely to be a common user error
+# with the API. There is also the option that a future change
+# would support lengths, so this acts as a check of the API (and
+# is why a non-length is also included in the test).
+#
+@pytest.mark.parametrize("rad", [[0, 0.1, 0.2] * u.Mpc,
+ [0.1, 0.2] * u.kg])
+def test_fail_radii_convert(rad):
+ """Not an angle."""
+
+ estr = "Argument 'radii' to function '__init__' " + \
+ "must be in units convertible to 'angle'."
+ with pytest.raises(u.UnitsError, match=estr):
+ Deproject(rad)
+
+
+@pytest.mark.parametrize("rad", [[], [10]])
+def test_fail_radii_num(rad):
+ """Need at least 2 radii"""
+
+ estr = 'radii parameter must be a sequence with at least two values'
+ with pytest.raises(ValueError, match=estr):
+ Deproject(rad * u.arcmin)
+
+
+@pytest.mark.parametrize("rad", [[0, 0], [0, -10], [0, 10, 10, 20]])
+def test_fail_radii_order(rad):
+ """Fail if the radii are not in ascending order"""
+
+ estr = 'radii parameter must be in increasing order'
+ with pytest.raises(ValueError, match=estr):
+ Deproject(rad * u.arcsec)
+
+
+def test_fail_negative_radii():
+ """Fail if radius is negative"""
+
+ estr = 'radii must be >= 0'
+ with pytest.raises(ValueError, match=estr):
+ Deproject([-2, 0, 2] * u.rad)
+
+
+def test_fail_extra_annulus():
+ """We fail if more data than annuli"""
+
+ dep = Deproject([0, 20] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, 0)
+
+ d1 = _test_data(1)
+
+ estr = 'Expected 0 <= annulus < 1 but sent 1'
+ with pytest.raises(ValueError, match=estr):
+ dep.load_pha(d1, 1)
+
+
+def test_vnorm_1bin():
+ """Check volume normalization when there's one bin"""
+
+ dep = Deproject([0, 20] * u.arcmin)
+ d0 = _test_data()
+ dep.load_pha(d0, 0)
+
+ assert not hasattr(dep, 'vol_norm')
+ dep.set_source()
+ assert hasattr(dep, 'vol_norm')
+
+ vnorm = dep.vol_norm
+ assert vnorm.shape == (1, 1)
+ assert dep.vol_norm[0, 0] == pytest.approx(1.0)
+
+
+@pytest.mark.parametrize("radii", [(0, 10, 20), (10, 30)])
+def test_vnorm_unit_independence(radii):
+ """The normalization terms are independent of the radius units"""
+
+ dep1 = Deproject(radii * u.arcsec)
+ dep2 = Deproject(radii * u.rad)
+
+ dep1._calc_vol_norm()
+ dep2._calc_vol_norm()
+
+ for v1, v2 in zip(dep1.vol_norm.flatten(),
+ dep2.vol_norm.flatten()):
+ diff = v1 - v2
+ assert diff == pytest.approx(0.0)
+
+
+def test_fail_source_missing_annulus1():
+ """Does set_source fail if missing an annulus?"""
+
+ dep = Deproject([0, 20] * u.arcsec)
+
+ estr = 'missing data for annulus 0'
+ with pytest.raises(ValueError, match=estr):
+ dep.set_source()
+
+
+def test_fail_source_missing_annulus2():
+ """Does set_source fail if missing multiple annuli?"""
+
+ dep = Deproject([0, 20, 40] * u.arcsec)
+
+ # We do not want this treated as a regexp
+ estr = r'missing data for annuli \[0, 1\]'
+ with pytest.raises(ValueError, match=estr):
+ dep.set_source()
+
+
+@pytest.mark.parametrize("n,m", [(0, 1), (1, 0)])
+def test_fail_source_missing_annulus3(n, m):
+ """Does set_source fail if missing some annuliu"""
+
+ dep = Deproject([0, 20, 40] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, n)
+
+ estr = 'missing data for annulus {}'.format(m)
+ with pytest.raises(ValueError, match=estr):
+ dep.set_source()
+
+
+@pytest.mark.parametrize("n,m", [(0, r"\[1, 2\]"),
+ (1, r"\[0, 2\]"),
+ (2, r"\[0, 1\]")])
+def test_fail_source_missing_annulus4(n, m):
+ """Does set_source fail if missing some annuli.
+
+ Note the need to protect the expected string so that match
+ does not take it as a regexp
+ """
+
+ dep = Deproject([0, 20, 40, 60] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, n)
+
+ estr = 'missing data for annuli {}'.format(m)
+ with pytest.raises(ValueError, match=estr):
+ dep.set_source()
+
+
+@pytest.mark.parametrize("srcmodel,comps,objs,mcomps",
+ [(None,
+ [{'type': 'xsphabs', 'start': 0, 'end': 7},
+ {'type': 'xsapec', 'start': 8, 'end': 14}],
+ [XSphabs('xsphabs.xsphabs_0'),
+ XSapec('xsapec.xsapec_0')],
+ [{'type': 'xsphabs',
+ 'name': 'xsphabs_0',
+ 'shell': 0},
+ {'type': 'xsapec',
+ 'name': 'xsapec_0',
+ 'shell': 0}]),
+
+ ('xswabs * xsmekal',
+ [{'type': 'xswabs', 'start': 0, 'end': 6},
+ {'type': 'xsmekal', 'start': 9, 'end': 16}],
+ [XSwabs('xswabs.xswabs_0'),
+ XSmekal('xsmekal.xsmekal_0')],
+ [{'type': 'xswabs',
+ 'name': 'xswabs_0',
+ 'shell': 0},
+ {'type': 'xsmekal',
+ 'name': 'xsmekal_0',
+ 'shell': 0}])])
+def test_set_source(srcmodel, comps, objs, mcomps):
+ """Check the source model can be set, with 1 annulus.
+
+ Basic checks of the structures set in place by a set_source
+ call.
+ """
+
+ dep = Deproject([0, 20] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, 0)
+
+ if srcmodel is None:
+ dep.set_source()
+ assert dep.srcmodel == 'xsphabs*xsapec'
+
+ else:
+ dep.set_source(srcmodel)
+ assert dep.srcmodel == srcmodel
+
+ # can do a direct comparison of srcmodel_comps
+ assert dep.srcmodel_comps == comps
+
+ # for model_comps, want to remove the object key from the comparison
+ #
+ to_compare = []
+ to_objs = []
+ for mc in dep.model_comps:
+ to_objs.append(mc['object'])
+ mc = mc.copy()
+ del mc['object']
+ to_compare.append(mc)
+
+ assert to_compare == mcomps
+ for a, b in zip(objs, to_objs):
+ assert a.type == b.type
+ assert a.name == b.name
+
+
+@pytest.mark.parametrize("nmodel", [None,
+ "xsphabs * xsapec",
+ "xswabs * xsmekal"])
+def test_reset_source(nmodel):
+ """Check the source model can be reset.
+
+ Note there's no check on the downstream implications of this
+ change, as this is just to check the initial bug.
+ """
+
+ dep = Deproject([0, 20] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, 0)
+
+ dep.set_source()
+ if nmodel is None:
+ dep.set_source()
+ else:
+ dep.set_source(nmodel)
+
+
+def test_reset_source_explicit():
+ """Check the before and after settings of set_source.
+
+ This overlaps previous checks, but uses simpler models
+ (single component).
+ """
+
+ dep = Deproject([0, 20] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, 0)
+
+ dep.set_source('xsvapec')
+ assert dep.srcmodel_comps == [{'type': 'xsvapec', 'start': 0, 'end': 7}]
+ assert len(dep.model_comps) == 1
+ mcomps = dep.model_comps[0]
+ assert len(mcomps) == 4
+ assert mcomps['type'] == 'xsvapec'
+ assert mcomps['name'] == 'xsvapec_0'
+ assert mcomps['shell'] == 0
+ assert mcomps['object'].type == 'xsvapec'
+ assert mcomps['object'].name == 'xsvapec.xsvapec_0'
+
+ dep.set_source('xsvvapec')
+ assert dep.srcmodel_comps == [{'type': 'xsvvapec', 'start': 0, 'end': 8}]
+ assert len(dep.model_comps) == 1
+ mcomps = dep.model_comps[0]
+ assert len(mcomps) == 4
+ assert mcomps['type'] == 'xsvvapec'
+ assert mcomps['name'] == 'xsvvapec_0'
+ assert mcomps['shell'] == 0
+ assert mcomps['object'].type == 'xsvvapec'
+ assert mcomps['object'].name == 'xsvvapec.xsvvapec_0'
+
+
+def test_set_redshift_no_model_get():
+ """redshift attribute."""
+
+ dep = Deproject([0, 20] * u.arcsec)
+
+ estr = 'Parameter redshift not found in any model component'
+ with pytest.raises(ValueError, match=estr):
+ dep.redshift
+
+
+def test_set_redshift_no_model_set():
+ """redshift attribute."""
+
+ dep = Deproject([0, 20] * u.arcsec)
+ dep.redshift = 0.12
+ assert dep.redshift == pytest.approx(0.12)
+
+
+def test_set_redshift_model_get():
+ """redshift attribute."""
+
+ dep = Deproject([0, 20, 40] * u.arcsec)
+ for ann in [0, 1]:
+ d = _test_data(ann)
+ dep.load_pha(d, ann)
+
+ # Note: does not fail if the redshifts are different
+ dep.set_source()
+ ui.get_model_component('xsapec_0').redshift = 0.12
+ assert dep.redshift == pytest.approx(0.12)
+
+
+def test_set_redshift_model_change_get():
+ """redshift attribute"""
+
+ dep = Deproject([0, 20] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, 0)
+
+ dep.set_source()
+ ui.get_model_component('xsapec_0').redshift = 0.12
+ ui.get_model_component('xsapec_0').redshift = 0.21
+ assert dep.redshift == pytest.approx(0.21)
+
+
+def test_set_redshift_model_get_change_get():
+ """redshift attribute
+
+ At present the redshift, once set, is not changed.
+ """
+
+ dep = Deproject([0, 20] * u.arcsec)
+ d0 = _test_data()
+ dep.load_pha(d0, 0)
+
+ dep.set_source()
+ ui.get_model_component('xsapec_0').redshift = 0.12
+
+ # ask for the redshift, but ignore the result
+ dep.redshift
+
+ # change the redshift
+ ui.get_model_component('xsapec_0').redshift = 0.21
+
+ # redshift has not changed
+ assert dep.redshift == pytest.approx(0.12)
+
+
+def test_init_angdist_type_fail():
+ """Does setting angdist to an invalid value (no type).
+ """
+
+ estr = "Argument 'angdist' to function '__init__' " + \
+ "has no 'unit' attribute. You may want to pass " + \
+ "in an astropy Quantity instead."
+ with pytest.raises(TypeError, match=estr):
+ Deproject([0, 10, 20] * u.arcsec, angdist=123.4)
+
+
+def test_set_angdist_type_fail():
+ """Does setting the angdist attribute with no type fail?"""
+
+ dep = Deproject([0, 10, 20] * u.arcsec)
+
+ estr = "Argument 'angdist' to function '_set_angdist' " + \
+ "has no 'unit' attribute. You may want to pass " + \
+ "in an astropy Quantity instead."
+ with pytest.raises(TypeError, match=estr):
+ dep.angdist = 123.4
+
+
+def test_init_angdist_convert_fail():
+ """Does setting angdist to an invalid type fail.
+ """
+
+ estr = "Argument 'angdist' to function '__init__' " + \
+ "must be in units convertible to 'length'."
+ with pytest.raises(u.UnitsError, match=estr):
+ Deproject([0, 10, 20] * u.arcsec, angdist=123.4 * u.J)
+
+
+def test_set_angdist_convert_fail():
+ """Does setting angdist to an invalid type fail.
+ """
+
+ dep = Deproject([0, 10, 20] * u.arcsec)
+
+ estr = "Argument 'angdist' to function '_set_angdist' " + \
+ "must be in units convertible to 'length'."
+ with pytest.raises(u.UnitsError, match=estr):
+ dep.angdist = 123.4 * u.J
+
+
+@pytest.mark.parametrize("angdist", [0.0 * u.Mpc, -1.2e24 * u.cm])
+def test_init_angdist_invalid_fail(angdist):
+ """Does setting angdist to an invalid quantity fail.
+ """
+
+ with pytest.raises(ValueError):
+ Deproject([0, 10, 20] * u.arcsec, angdist=angdist)
+
+
+@pytest.mark.parametrize("angdist", [0.0 * u.Mpc, -1.2e24 * u.cm])
+def test_set_angdist_invalid_fail(angdist):
+ """Does setting angdist to an invalid quantity fail.
+ """
+
+ dep = Deproject([0, 10, 20] * u.arcsec)
+
+ with pytest.raises(ValueError):
+ dep.angdist = angdist
+
+
+@pytest.mark.parametrize("theta", [None, "theta", 0, 360, [10]])
+def test_fail_theta_type(theta):
+ """Need units for theta"""
+
+ estr = "Argument 'theta' to function '__init__' " + \
+ "has no 'unit' attribute. You may want to pass " + \
+ "in an astropy Quantity instead."
+ with pytest.raises(TypeError, match=estr):
+ Deproject([0, 10, 20] * u.arcsec, theta=theta)
+
+
+@pytest.mark.parametrize("theta", [10 * u.Mpc, [0.1, 0.2] * u.kg])
+def test_fail_theta_convert(theta):
+ """Not an angle."""
+
+ estr = "Argument 'theta' to function '__init__' " + \
+ "must be in units convertible to 'angle'."
+ with pytest.raises(u.UnitsError, match=estr):
+ Deproject([0, 10, 20] * u.arcmin, theta=theta)
+
+
+@pytest.mark.parametrize("theta", [0 * u.deg, 0 * u.rad,
+ -10 * u.deg, -1 * u.rad])
+def test_theta_out_of_band_small_fail(theta):
+ """These theta values are invalid."""
+
+ emsg = 'theta must be > 0 degrees'
+ with pytest.raises(ValueError, match=emsg):
+ Deproject([0, 10, 20] * u.arcsec, theta=theta)
+
+
+@pytest.mark.parametrize("theta", [361 * u.deg, 2.1 * np.pi * u.rad])
+def test_theta_out_of_band_large_fail(theta):
+ """These theta values are invalid."""
+
+ emsg = 'theta must be <= 360 degrees'
+ with pytest.raises(ValueError, match=emsg):
+ Deproject([0, 10, 20] * u.arcsec, theta=theta)
+
+
+@pytest.mark.parametrize("theta", [(10, 20), (30, 40, 350, 360)])
+def test_theta_nwrong_fail(theta):
+ """The number of theta are wrong."""
+
+ emsg = 'theta must be a scalar or match the number of annuli'
+ with pytest.raises(ValueError, match=emsg):
+ Deproject([0, 1, 2, 3] * u.arcmin, theta=theta * u.deg)
+
+
+@pytest.mark.parametrize("full", [360 * u.deg, 2 * np.pi * u.rad])
+@pytest.mark.parametrize("scale", [1, 0.5, 0.1])
+@pytest.mark.parametrize("radii", [(0, 10), (0, 10, 20), (0, 10, 40),
+ (10, 20), (10, 20, 40)])
+def test_theta_subset_norm(full, scale, radii):
+ """Does changing theta scale the normalization?
+ """
+
+ radii = radii * u.arcmin
+ nann = len(radii) - 1
+ dep_full = Deproject(radii, theta=full)
+ dep_scale = Deproject(radii, theta=full * scale)
+
+ for i in range(nann):
+ d0 = _test_data()
+ dep_full.load_pha(d0, i)
+ dep_scale.load_pha(d0, i)
+
+ dep_full.set_source()
+ dep_scale.set_source()
+
+ v_full = dep_full.vol_norm
+ v_scale = dep_scale.vol_norm
+ assert v_full.shape == (nann, nann)
+ assert v_scale.shape == v_full.shape
+
+ # the following would be slightly shorter with NumPy-array-aware
+ # comparison routines
+ #
+ # the matrix is in lower-triangular form, so the upper triangular
+ # area (above the diagonal) is 0.
+ #
+ if nann > 1:
+ zeros = np.triu_indices(nann, k=1)
+ for z in v_full[zeros]:
+ assert z == pytest.approx(0)
+
+ for z in v_scale[zeros]:
+ assert z == pytest.approx(0)
+
+ # the lower-triangular area (including the diagonal) is positive
+ # and scales with theta
+ #
+ idx = np.tril_indices(nann, k=0)
+ z_full = v_full[idx]
+ z_scale = v_scale[idx]
+ assert np.all(z_full > 0.0)
+ assert np.all(z_scale > 0.0)
+
+ for r in z_scale / z_full:
+ assert r == pytest.approx(scale)
+
+
+@pytest.mark.parametrize("full", [360 * u.deg, 2 * np.pi * u.rad])
+@pytest.mark.parametrize("scales,radii", [((1, 1), (0, 10, 20)),
+ ((1, 1), (10, 20, 40)),
+ ((0.4, 0.8), (0, 10, 20)),
+ ((0.8, 0.4), (10, 30, 40)),
+ ((0.5, 1.0, 0.7), (0, 2, 5, 8)),
+ ((0.5, 1.0, 0.7), (1, 2, 5, 8))])
+def test_theta_per_annuli_norm(full, scales, radii):
+ """Does changing theta (per annulus) scale the normalization?
+
+ This complements test_theta_per_annuli_norm by having the
+ theta value vary per annulus.
+ """
+
+ radii = radii * u.arcmin
+ nann = len(radii) - 1
+ scales = np.asarray(scales)
+ thetas = scales * full
+
+ dep_full = Deproject(radii, theta=full)
+ dep_scale = Deproject(radii, theta=thetas)
+
+ for i in range(nann):
+ d0 = _test_data()
+ dep_full.load_pha(d0, i)
+ dep_scale.load_pha(d0, i)
+
+ dep_full.set_source()
+ dep_scale.set_source()
+
+ v_full = dep_full.vol_norm
+ v_scale = dep_scale.vol_norm
+ assert v_full.shape == (nann, nann)
+ assert v_scale.shape == v_full.shape
+
+ # the following would be slightly shorter with NumPy-array-aware
+ # comparison routines
+ #
+ # the matrix is in lower-triangular form, so the upper triangular
+ # area (above the diagonal) is 0.
+ #
+ if nann > 1:
+ zeros = np.triu_indices(nann, k=1)
+ for z in v_full[zeros]:
+ assert z == pytest.approx(0)
+
+ for z in v_scale[zeros]:
+ assert z == pytest.approx(0)
+
+ # the lower-triangular area (including the diagonal) is positive
+ # and scales with theta
+ #
+ idx = np.tril_indices(nann, k=0)
+ z_full = v_full[idx]
+ z_scale = v_scale[idx]
+ assert np.all(z_full > 0.0)
+ assert np.all(z_scale > 0.0)
+
+ # Stack the scales array up to create a n by n matrix, so that
+ # we can apply the lower-triangle index and get the expected
+ # scale factor
+ scales_mat = np.resize(scales, (nann, nann))[idx]
+
+ for r, scale in zip(z_scale / z_full, scales_mat):
+ assert r == pytest.approx(scale)
+
+
+@pytest.mark.parametrize("units", ["arcsec", "deg", "rad",
+ "Mpc", "kpc", "cm",
+ u.Unit("arcmin"), u.Unit("pc")])
+def test_get_radii_accepts(units):
+ """Can we ask for a variety of units.
+
+ There is limited checking of the output.
+ """
+
+ dep = Deproject([0, 10, 15] * u.arcsec,
+ angdist=127.4 * u.Mpc)
+
+ rlo, rhi = dep.get_radii(units)
+
+ assert len(rlo) == len(rhi)
+ assert len(rlo) == 2
+
+ # Ensure the output is the correct type (length or angle)
+ #
+ rlo.to(u.Unit(units))
+
+ for r1, r2 in zip(rlo, rhi):
+ assert r2 > r1
+
+ # Relies on radii[2] = 1.5 * radii[1] in the input to Deproject.
+ # Apparently need to extract the value from the quantity otherwise
+ # the test fails due to unbounded recursion.
+ #
+ delta = rhi[1] - 1.5 * rhi[0]
+ assert delta.value == pytest.approx(0)
+
+
+def test_get_radii_one_annulus():
+ """Nothing breaks down when there's only one annulus"""
+
+ dep = Deproject([4.8e-5, 7.2e-5] * u.rad)
+ dep.angdist = 20.0 * u.Mpc
+ rlo, rhi = dep.get_radii('kpc')
+
+ assert len(rlo) == 1
+ assert len(rhi) == 1
+
+ assert rlo.unit.physical_type == 'length'
+ assert rhi.unit.physical_type == 'length'
+
+ assert rlo.unit == u.Unit('kpc')
+ assert rhi.unit == u.Unit('kpc')
+
+ # Using an absolute tolerance just as a regression test. Since
+ # the angular diameter distance explicitly has been given then
+ # there is no issue about changing the default Cosmology.
+ #
+ assert rlo[0].value == pytest.approx(0.960, abs=0.001)
+ assert rhi[0].value == pytest.approx(1.440, abs=0.001)
+
+
+@pytest.mark.parametrize("radii", [(0, 10), (0, 10, 20), (10, 30, 50, 90)])
+def test_get_shells_no_data(radii):
+ """No data."""
+
+ dep = Deproject(radii * u.arcsec)
+ shells = dep.get_shells()
+ assert len(shells) == dep.nshell
+
+ # shells is in outer-to-inner order
+ shells.reverse()
+ for i, shell in enumerate(shells):
+ assert shell['annuli'] == [i]
+ assert shell['dataids'] == []
+
+
+@pytest.mark.parametrize("radii", [(0, 10), (0, 10, 20), (10, 30, 50, 90)])
+def test_get_shells_no_ties(radii):
+ """No ties between the shells"""
+
+ dep = Deproject(radii * u.arcsec)
+ for i in range(dep.nshell):
+ d = _test_data(i)
+ dep.load_pha(d, i)
+
+ dep.set_source()
+ shells = dep.get_shells()
+ assert len(shells) == dep.nshell
+
+ # shells is in outer-to-inner order
+ shells.reverse()
+ for i, shell in enumerate(shells):
+ assert shell['annuli'] == [i]
+ assert shell['dataids'] == [i]
+
+
+@pytest.mark.parametrize("inner,outer", [(0, 3), (3, 0), (1, 3), (3, 1)])
+def test_get_shells_not_consecutive_fail(inner, outer):
+ """Can not tie parameters across groups.
+
+ This is only checked when the annuli grouping is required,
+ which is simulated with a call to get_radii.
+ """
+
+ dep = Deproject([1, 2, 3, 4, 5] * u.arcsec)
+ for i in range(dep.nshell):
+ d = _test_data(i)
+ dep.load_pha(d, i)
+
+ dep.set_source('xsvapec')
+ dep.tie_par('xsvapec.ne', inner, outer)
+
+ anns = [str(s) for s in sorted([inner, outer])]
+ estr = r'Non-consecutive annuli are tied together: \[' + \
+ r'{}\]'.format(" ".join(anns))
+ with pytest.raises(ValueError, match=estr):
+ dep.get_shells()
+
+
+@pytest.mark.parametrize("outer,annuli,dataids",
+ [(3, [[2, 3], [1], [0]], [[0, 1], [2], [3]]),
+ (2, [[3], [1, 2], [0]], [[0], [1, 2], [3]]),
+ (1, [[3], [2], [0, 1]], [[0], [1], [2, 3]])])
+def test_tie_par1(outer, annuli, dataids):
+ """One pair of shells is tied."""
+
+ dep = Deproject([1, 2, 3, 4, 5] * u.arcsec)
+
+ # Reverse the order so that annulus != dataid
+ for i in [3, 2, 1, 0]:
+ d = _test_data(i)
+ dep.load_pha(d, i)
+
+ dep.set_source()
+ inner = outer - 1
+ dep.tie_par('xsapec.kt', inner, outer)
+
+ shells = dep.get_shells()
+ assert len(shells) == 3
+
+ for shell, ann, dids in zip(shells, annuli, dataids):
+ assert shell['annuli'] == ann
+ assert shell['dataids'] == dids
+
+
+@pytest.mark.parametrize("outer,annuli,dataids",
+ [(3, [[1, 2, 3], [0]], [[0, 1, 2], [3]]),
+ (2, [[3], [0, 1, 2]], [[0], [1, 2, 3]])])
+def test_tie_par2(outer, annuli, dataids):
+ """Two pair of shells is tied."""
+
+ dep = Deproject([1, 2, 3, 4, 5] * u.arcsec)
+
+ # Reverse the order so that annulus != dataid
+ for i in [3, 2, 1, 0]:
+ d = _test_data(i)
+ dep.load_pha(d, i)
+
+ dep.set_source()
+ inner = outer - 2
+ dep.tie_par('xsapec.kt', inner, outer-1, outer)
+
+ shells = dep.get_shells()
+ assert len(shells) == 2
+
+ for shell, ann, dids in zip(shells, annuli, dataids):
+ assert shell['annuli'] == ann
+ assert shell['dataids'] == dids
+
+
+@pytest.mark.parametrize("outer", [3, 2, 1])
+def test_untie_par1(outer):
+ """Can untie a parameter."""
+
+ nshell = 4
+ dep = Deproject([1, 2, 3, 4, 5] * u.arcsec)
+
+ # Reverse the order so that annulus != dataid
+ for i in [3, 2, 1, 0]:
+ d = _test_data(i)
+ dep.load_pha(d, i)
+
+ dep.set_source()
+ inner = outer - 1
+
+ # Rely on previous tests to check that tie_par works
+ dep.tie_par('xsapec.kt', inner, outer)
+
+ # untie the parameter
+ dep.untie_par('xsapec.kt', outer)
+
+ shells = dep.get_shells()
+ assert len(shells) == nshell
+
+ shells.reverse()
+ for i, shell in enumerate(shells):
+ assert shell['annuli'] == [i]
+ assert shell['dataids'] == [nshell - 1 - i]
+
+
+@pytest.mark.parametrize("outer", [3, 2])
+def test_untie_par2(outer):
+ """Can untie two parameters."""
+
+ nshell = 4
+ dep = Deproject([1, 2, 3, 4, 5] * u.arcsec)
+
+ # Reverse the order so that annulus != dataid
+ for i in [3, 2, 1, 0]:
+ d = _test_data(i)
+ dep.load_pha(d, i)
+
+ dep.set_source()
+ inner = outer - 2
+
+ # Rely on previous tests to check that tie_par works
+ dep.tie_par('xsapec.kt', inner, outer - 1, outer)
+
+ # untie the parameter
+ dep.untie_par('xsapec.kt', outer - 1, outer)
+
+ shells = dep.get_shells()
+ assert len(shells) == nshell
+
+ shells.reverse()
+ for i, shell in enumerate(shells):
+ assert shell['annuli'] == [i]
+ assert shell['dataids'] == [nshell - 1 - i]
diff --git a/tests/test_simplegraph.py b/tests/test_simplegraph.py
new file mode 100644
index 0000000..2e66858
--- /dev/null
+++ b/tests/test_simplegraph.py
@@ -0,0 +1,161 @@
+"""
+Limited testing of deproject.simplegraph
+
+:Copyright: Smithsonian Astrophysical Observatory (2019)
+:Author: Douglas Burke (dburke@cfa.harvard.edu)
+"""
+
+import pytest
+
+from deproject.simplegraph import SimpleGraph
+
+
+def _sort_groups(gs):
+ return sorted([sorted(g) for g in gs])
+
+
+def test_empty():
+ """No nodes."""
+
+ g = SimpleGraph()
+ groups = g.get_groups()
+ assert groups == []
+
+
+def test_0d():
+ """Everything is one."""
+
+ g = SimpleGraph()
+ g.add_link(23, 23)
+ groups = g.get_groups()
+ assert groups == [[23]]
+
+
+def test_0d_really():
+ """Everything is one even if we like to repeat."""
+
+ g = SimpleGraph()
+ g.add_link(' ', ' ')
+ g.add_link(' ', ' ')
+ groups = g.get_groups()
+ assert groups == [[' ']]
+
+
+def test_single_link_same_type():
+ """A single link (same type)"""
+
+ g = SimpleGraph()
+ g.add_link('x ', ' y')
+ groups = _sort_groups(g.get_groups())
+ assert groups == [[' y', 'x ']]
+
+
+def test_single_link_diff_type():
+ """A single link with different types for the nodes"""
+
+ g = SimpleGraph()
+ g.add_link(' x ', True)
+
+ # can not sort this list, so explicitly check
+ groups = g.get_groups()
+ assert len(groups) == 1
+ assert len(groups[0]) == 2
+ assert ' x ' in groups[0]
+ assert True in groups[0]
+
+
+def test_odd_nodes_unconnected():
+ """Self link plus another link, unconnected (and test node types)"""
+
+ g = SimpleGraph()
+ g.add_link('', '')
+ g.add_link(300, 2)
+
+ # can not sort this list, so explicitly check
+ groups = g.get_groups()
+ assert len(groups) == 2
+ g1 = groups[0]
+ g2 = groups[1]
+
+ def pred1(g):
+ """Is the group the ''->'' link?"""
+ return len(g) == 1 and g[0] == ''
+
+ def pred2(g):
+ """Is the group the 2->300 link?"""
+ return len(g) == 2 and 2 in g and 300 in g
+
+ assert (pred1(g1) and pred2(g2)) or (pred1(g2) and pred2(g1))
+
+
+def test_odd_nodes_connected():
+ """Connect up nodes of different types"""
+
+ g = SimpleGraph()
+ g.add_link(True, 5)
+ g.add_link(5, False)
+
+ groups = g.get_groups()
+ assert len(groups) == 1
+ assert len(groups[0]) == 3
+ assert True in groups[0]
+ assert False in groups[0]
+ assert 5 in groups[0]
+
+
+# Could parametrize many of the tests since using the same
+# basic data, but leave separate for now.
+#
+def test_no_pairs():
+ """no pairs"""
+
+ g = SimpleGraph()
+
+ for i in range(5, 11):
+ g.add_link(i, i)
+
+ groups = g.get_groups()
+ assert len(groups) == 6
+
+ vs = _sort_groups(groups)
+ assert vs == [[5], [6], [7], [8], [9], [10]], vs
+
+
+@pytest.mark.parametrize("n1,n2", [(8, 9), (9, 8)])
+def test_one_pair(n1, n2):
+ """one pair"""
+
+ g = SimpleGraph()
+
+ for i in range(5, 11):
+ g.add_link(i, i)
+
+ g.add_link(n1, n2)
+
+ groups = g.get_groups()
+ assert len(groups) == 5
+
+ vs = _sort_groups(groups)
+ assert vs == [[5], [6], [7], [8, 9], [10]], vs
+
+
+@pytest.mark.parametrize("n1,n2,m1,m2", [(8, 9, 7, 8),
+ (9, 8, 7, 8),
+ (8, 9, 8, 7),
+ (9, 8, 8, 7)])
+def test_two_pairs(n1, n2, m1, m2):
+ """two pairs"""
+
+ g = SimpleGraph()
+
+ for i in range(5, 11):
+ g.add_link(i, i)
+
+ g.add_link(n1, n2)
+ g.add_link(m1, m2)
+
+ groups = g.get_groups()
+ assert len(groups) == 4
+
+ vs = _sort_groups(groups)
+ assert vs == [[5], [6], [7, 8, 9], [10]], vs