From 81648e329ac901d13cc618e20c2049b6e9b24d09 Mon Sep 17 00:00:00 2001 From: nocturnalastro Date: Fri, 19 Aug 2016 20:59:34 +0100 Subject: [PATCH 1/4] split 1Dpsf from master --- .travis.yml | 2 -- docs/gen_plots.py | 1 - docs/index.rst | 7 ++++++- saba/tests/coveragerc | 2 +- saba/tests/test_main.py | 22 ++++++++++++++++++++++ 5 files changed, 29 insertions(+), 5 deletions(-) diff --git a/.travis.yml b/.travis.yml index 59178a1..76898ec 100644 --- a/.travis.yml +++ b/.travis.yml @@ -46,7 +46,6 @@ env: # Try all python versions with the latest numpy - SETUP_CMD='test' - matrix: include: @@ -73,7 +72,6 @@ matrix: - python: 3.5 env: NUMPY_VERSION=1.10 - install: # We now use the ci-helpers package to set up our testing environment. diff --git a/docs/gen_plots.py b/docs/gen_plots.py index fd81c71..b4de26d 100644 --- a/docs/gen_plots.py +++ b/docs/gen_plots.py @@ -1,6 +1,5 @@ from astropy.modeling.fitting import SherpaFitter from astropy.modeling.models import Gaussian1D, Gaussian2D - import numpy as np import matplotlib.pyplot as plt diff --git a/docs/index.rst b/docs/index.rst index 937ac8e..b9111a3 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -45,6 +45,7 @@ To make use of the entry points plugin registry which automatically makes the |S Otherwise one can just use the latest stable ``astropy`` via:: conda install astropy + Next install Sherpa_ using the conda ``sherpa`` channel. Note that Sherpa currently needs to be installed after astropy on Mac OSX. @@ -236,4 +237,8 @@ API/Reference Credit ------ -The development of this package was made possible by the generous support of the `Google Summer of Code `_ program in 2016 under the `OpenAstronomy `_ by `Michele Costa `_ with the support and advice of mentors `Tom Aldcroft `_, `Omar Laurino `_, `Moritz Guenther `_, and `Doug Burke `_. +The development of this package was made possible by the generous support of the `Google Summer of Code `_ program in 2016 +under the `OpenAstronomy `_ +by `Michele Costa `_ with the support and advice of mentors +`Tom Aldcroft `_, `Omar Laurino `_, +`Moritz Guenther `_, and `Doug Burke `_. diff --git a/saba/tests/coveragerc b/saba/tests/coveragerc index 3a21984..3e5082c 100644 --- a/saba/tests/coveragerc +++ b/saba/tests/coveragerc @@ -18,7 +18,7 @@ exclude_lines = pragma: no cover # Don't complain about packages we have installed - # except ImportError + except ImportError # Don't complain if tests don't hit assertions raise AssertionError diff --git a/saba/tests/test_main.py b/saba/tests/test_main.py index 67489c9..9c91b15 100644 --- a/saba/tests/test_main.py +++ b/saba/tests/test_main.py @@ -289,6 +289,28 @@ def test_bkg_doesnt_explode(self): sfit(m, x, y, bkg=bkg) # TODO: Make this better! + + def test_rsp1d_doesnt_explode(self): + """ + Check this goes through the motions + """ + + self.fitter(self.model1d.copy(), self.x1, self.y1, err=self.dy1, rsp=self.rsp1) + + def test_rsp1d_multi_doesnt_explode(self): + """ + Check this goes through the motions + """ + + self.fitter([self.model1d.copy(), self.model1d_2.copy()], [self.x1, self.x2], [self.y1, self.y2], err=[self.dy1, self.dy2], rsp=[self.rsp1, self.rsp2]) + ''' + def test_rsp2d_doesnt_explode(self): + """ + Check this goes through the motions + """ + self.fitter(self.model2d.copy(), self.xx1, self.xx2, self.yy, err=self.dyy, rsp=self.rsp2d) + ''' + def test_entry_points(self): # a little to test that entry points can be loaded! from pkg_resources import iter_entry_points From 45c5d3082d69dc34d9f52d573efa56968f81a3a9 Mon Sep 17 00:00:00 2001 From: nocturnalastro Date: Fri, 19 Aug 2016 21:30:12 +0100 Subject: [PATCH 2/4] added make_rsp with a block on 2D and I even have tests --- saba/main.py | 80 +++++++++++++++++++++++++++++++++++++---- saba/tests/test_main.py | 7 ---- 2 files changed, 73 insertions(+), 14 deletions(-) diff --git a/saba/main.py b/saba/main.py index afe9d83..1628dd1 100644 --- a/saba/main.py +++ b/saba/main.py @@ -27,8 +27,6 @@ if "SherpaFitter" not in w.message.message: warnings.warn(w) -# from astropy.modeling - __all__ = ('SherpaFitter', 'SherpaMCMC') @@ -268,6 +266,53 @@ def stat_values(self): return self._stat_values +def make_rsp(rsp): + """ + Take an array as a response which is then convolved with the model output. + """ + def wrap_rsp(_data, _rsp): + """ + Take an array as a response which is then + convolved with the model output. + + data: a sherpa dataset + rsp : an array which represets rsp + """ + _rsp = np.asarray(_rsp) + _rdata = copy.deepcopy(_data) + _rdata.y = _rsp + _psf = PSFModel("user_rsp", _rdata) + _psf.fold(_data) + return _psf + + try: + ndims = len(data.data.datasets[0].get_dims()) + except AttributeError: + ndims = len(data.data.get_dims()) + + if ndims > 1: + return None + else: + rsp = np.asarray(rsp) + + if data.ndata > 1: + if rsp.ndim > 1 or rsp.dtype == np.object: + if rsp.shape[0] == data.ndata: + zipped = zip(data.data.datasets, rsp) + else: + raise AstropyUserWarning("There is more than 1 but not" + " ndata responses") + else: + zipped = zip(data.data.datasets, + [rsp for _ in xrange(data.ndata)]) + + rsp = [] + for da, rr in zipped: + rsp.append(wrap_rsp(da, rr)) + else: + return wrap_rsp(data.data, rsp) + + class SherpaFitter(Fitter): __doc__ = """ Sherpa Fitter for astropy models. @@ -364,10 +409,15 @@ def __call__(self, models, x, y, z=None, xbinsize=None, ybinsize=None, err=None, self._data = Dataset(n_inputs, x, y, z, xbinsize, ybinsize, err, bkg, bkg_scale) + if rsp is not None: + self._rsp = make_rsp(rsp) + else: + self._rsp = None + if self._data.ndata > 1: if len(models) == 1: - self._fitmodel = ConvertedModel([models.copy() for _ in xrange(self._data.ndata)], tie_list) + self._fitmodel = ConvertedModel([models.copy() for _ in xrange(self._data.ndata)], tie_list, rsp=self._rsp) # Copy the model so each data set has the same model! elif len(models) == self._data.ndata: self._fitmodel = ConvertedModel(models, tie_list) @@ -377,9 +427,10 @@ def __call__(self, models, x, y, z=None, xbinsize=None, ybinsize=None, err=None, else: if len(models) > 1: self._data.make_simfit(len(models)) - self._fitmodel = ConvertedModel(models, tie_list) + self._fitmodel = ConvertedModel(models, tie_list, + rsp=self._rsp) else: - self._fitmodel = ConvertedModel(models) + self._fitmodel = ConvertedModel(models, rsp=self._rsp) self._fitter = Fit(self._data.data, self._fitmodel.sherpa_model, self._stat_method, self._opt_method, self._est_method, **kwargs) self.fit_info = self._fitter.fit() @@ -633,16 +684,31 @@ class ConvertedModel(object): e.g. [(modelB.y, modelA.x)] will mean that y in modelB will be tied to x of modelA """ - def __init__(self, models, tie_list=None): + def __init__(self, models, tie_list=None, rsp=None): self.model_dict = OrderedDict() try: models.parameters # does it quack self.sherpa_model = self._astropy_to_sherpa_model(models) + self.rsp = rsp + if rsp is not None: + self.sherpa_model = rsp(self.sherpa_model) + self.model_dict[models] = self.sherpa_model except AttributeError: - for mod in models: + try: + n_rsp = len(rsp) + assert len(models) == n_rsp, AstropyUserWarning("The number of responses must be either 1 or the numeber of models %i" % len(models)) + zipped = zip(models, rsp) + + except TypeError: + zipped = zip(models, [rsp for _ in range()]) + + for mod, rsp in zipped: self.model_dict[mod] = self._astropy_to_sherpa_model(mod) + if rsp is not None: + self.sherpa_model[mod] = rsp(self.sherpa_model[mod]) + if tie_list is not None: for par1, par2 in tie_list: getattr(self.model_dict[par1._model], par1.name).link = getattr(self.model_dict[par2._model], par2.name) diff --git a/saba/tests/test_main.py b/saba/tests/test_main.py index 9c91b15..db5cc63 100644 --- a/saba/tests/test_main.py +++ b/saba/tests/test_main.py @@ -303,13 +303,6 @@ def test_rsp1d_multi_doesnt_explode(self): """ self.fitter([self.model1d.copy(), self.model1d_2.copy()], [self.x1, self.x2], [self.y1, self.y2], err=[self.dy1, self.dy2], rsp=[self.rsp1, self.rsp2]) - ''' - def test_rsp2d_doesnt_explode(self): - """ - Check this goes through the motions - """ - self.fitter(self.model2d.copy(), self.xx1, self.xx2, self.yy, err=self.dyy, rsp=self.rsp2d) - ''' def test_entry_points(self): # a little to test that entry points can be loaded! From ffd1e4395938bcb52ec5bfc54d01b0697b1da2c0 Mon Sep 17 00:00:00 2001 From: nocturnalastro Date: Fri, 19 Aug 2016 22:11:51 +0100 Subject: [PATCH 3/4] added rsp keyword and discription --- saba/main.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/saba/main.py b/saba/main.py index 1628dd1..cd4c055 100644 --- a/saba/main.py +++ b/saba/main.py @@ -366,7 +366,7 @@ def __init__(self, optimizer="levmar", statistic="leastsq", estmethod="covarianc setattr(self.__class__, 'est_config', property(lambda s: s._est_config, doc=self._est_method.__doc__)) - def __call__(self, models, x, y, z=None, xbinsize=None, ybinsize=None, err=None, bkg=None, bkg_scale=1, **kwargs): + def __call__(self, models, x, y, z=None, xbinsize=None, ybinsize=None, err=None, bkg=None, bkg_scale=1, rsp=None, **kwargs): """ Fit the astropy model with a the sherpa fit routines. @@ -392,6 +392,9 @@ def __call__(self, models, x, y, z=None, xbinsize=None, ybinsize=None, err=None, bkg_sale : float or list of floats (optional) the scaling factor for the dataset if a single value is supplied it will be copied for each dataset + rsp: array or list of arrays + this is convolved with the model output when fitting the model + N.B only 1D is currently supported. **kwargs : keyword arguments will be passed on to sherpa fit routine From 00fc6618380dc54db32bbc4bf2a053700a0af91a Mon Sep 17 00:00:00 2001 From: nocturnalastro Date: Fri, 19 Aug 2016 22:41:57 +0100 Subject: [PATCH 4/4] fixed errors --- saba/main.py | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/saba/main.py b/saba/main.py index cd4c055..e3daa00 100644 --- a/saba/main.py +++ b/saba/main.py @@ -1,18 +1,21 @@ from __future__ import (absolute_import, unicode_literals, division, print_function) -import numpy as np from collections import OrderedDict +import numpy as np +import warnings +import copy + from sherpa.fit import Fit from sherpa.data import Data1D, Data1DInt, Data2D, Data2DInt, DataSimulFit from sherpa.data import BaseData from sherpa.models import UserModel, Parameter, SimulFitModel +from sherpa.instrument import PSFModel from sherpa.stats import Chi2, Chi2ConstVar, Chi2DataVar, Chi2Gehrels from sherpa.stats import Chi2ModVar, Chi2XspecVar, LeastSq from sherpa.stats import CStat, WStat, Cash from sherpa.optmethods import GridSearch, LevMar, MonCar, NelderMead from sherpa.estmethods import Confidence, Covariance, Projection from sherpa.sim import MCMC -import warnings from astropy.utils import format_doc from astropy.utils.exceptions import AstropyUserWarning @@ -266,24 +269,21 @@ def stat_values(self): return self._stat_values -def make_rsp(rsp): +def make_rsp(data,rsp): """ Take an array as a response which is then convolved with the model output. + Parameters + ---------- + data: a sherpa dataset + rsp : an array which represets rsp """ - def wrap_rsp(_data, _rsp): - """ - Take an array as a response which is then - convolved with the model output. - - data: a sherpa dataset - rsp : an array which represets rsp - """ - _rsp = np.asarray(_rsp) - _rdata = copy.deepcopy(_data) - _rdata.y = _rsp - _psf = PSFModel("user_rsp", _rdata) - _psf.fold(_data) - return _psf + def wrap_rsp(data, rsp): + rsp = np.asarray(rsp) + rdata = copy.deepcopy(data) + rdata.y = rsp + psf = PSFModel("user_rsp", rdata) + psf.fold(data) + return psf try: ndims = len(data.data.datasets[0].get_dims()) @@ -413,7 +413,7 @@ def __call__(self, models, x, y, z=None, xbinsize=None, ybinsize=None, err=None, self._data = Dataset(n_inputs, x, y, z, xbinsize, ybinsize, err, bkg, bkg_scale) if rsp is not None: - self._rsp = make_rsp(rsp) + self._rsp = make_rsp(self._data, rsp) else: self._rsp = None @@ -704,7 +704,7 @@ def __init__(self, models, tie_list=None, rsp=None): zipped = zip(models, rsp) except TypeError: - zipped = zip(models, [rsp for _ in range()]) + zipped = zip(models, [rsp for _ in range(len(models))]) for mod, rsp in zipped: self.model_dict[mod] = self._astropy_to_sherpa_model(mod)