Skip to content

Commit ff60cef

Browse files
committed
Add meta and Poisson uncertaintiy
1 parent 02a7b57 commit ff60cef

8 files changed

Lines changed: 572 additions & 56 deletions

File tree

sunkit_spex/extern/ndcube/__init__.py

Whitespace-only changes.

sunkit_spex/extern/ndcube/meta.py

Lines changed: 373 additions & 0 deletions
Large diffs are not rendered by default.

sunkit_spex/spectrum/__init__.py

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
from sunkit_spex.spectrum.spectrum import Spectrum
2+
3+
__all__ = ["Spectrum"]

sunkit_spex/spectrum/spectrum.py

Lines changed: 16 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
from astropy.modeling.tabular import Tabular1D
99
from astropy.utils import lazyproperty
1010

11-
__all__ = ['gwcs_from_array', 'SpectralAxis', 'Spectrum', 'CountSpectrum']
11+
__all__ = ['gwcs_from_array', 'SpectralAxis', 'Spectrum']
1212

1313

1414
def gwcs_from_array(array):
@@ -158,6 +158,21 @@ class Spectrum(NDCube):
158158
meta : dict
159159
Arbitrary container for any user-specific information to be carried
160160
around with the spectrum container object.
161+
162+
Examples
163+
--------
164+
>>> import numpy as np
165+
>>> import astropy.units as u
166+
>>> from sunkit_spex.spectrum import Spectrum
167+
>>> spec = Spectrum(np.arange(1, 11)*u.watt, spectral_axis=np.arange(1, 12)*u.keV)
168+
>>> spec
169+
<sunkit_spex.spectrum.spectrum.Spectrum object at ...
170+
NDCube
171+
------
172+
Dimensions: [10.] pix
173+
Physical Types of Axes: [('em.energy',)]
174+
Unit: W
175+
Data Type: float64
161176
"""
162177

163178
def __init__(self, data, *, uncertainty=None, spectral_axis=None,
@@ -220,45 +235,3 @@ def __init__(self, data, *, uncertainty=None, spectral_axis=None,
220235
super().__init__(
221236
data=data.value if isinstance(data, u.Quantity) else data,
222237
wcs=wcs, mask=mask, meta=meta, uncertainty=uncertainty, **kwargs)
223-
224-
225-
class CountSpectrum(Spectrum):
226-
r"""
227-
Spectrum container for count data with one spectral axis.
228-
229-
Specifically, the data must be supplied as counts or convertable to counts by multiplying by the provided
230-
normalisation.
231-
232-
Parameters
233-
----------
234-
data : `~astropy.units.Quantity`
235-
The data for this spectrum. This can be a simple `~astropy.units.Quantity`,
236-
or an existing `~Spectrum1D` or `~ndcube.NDCube` object.
237-
norm : `~astropy.units.Quantity`
238-
The normalisation if the data unit is not counts then the product the unit of data and norm must be counts.
239-
uncertainty : `~astropy.nddata.NDUncertainty`
240-
Contains uncertainty information along with propagation rules for spectrum arithmetic. Can take a unit, but if
241-
none is given, will use the unit defined in the flux.
242-
spectral_axis : `~astropy.units.Quantity` or `~specutils.SpectralAxis`
243-
Dispersion information with the same shape as the dimension specified by spectral_dimension
244-
of shape plus one if specifying bin edges.
245-
spectral_dimension : `int` default 0
246-
The dimension of the data which represents the spectral information default to first dimension index 0.
247-
mask : `~numpy.ndarray`-like
248-
Array where values in the flux to be masked are those that
249-
``astype(bool)`` converts to True. (For example, integer arrays are not
250-
masked where they are 0, and masked for any other value.)
251-
meta : dict
252-
Arbitrary container for any user-specific information to be carried
253-
around with the spectrum container object.
254-
"""
255-
256-
def __init__(self, data, norm=None, **kwargs):
257-
if data.unit != u.ct:
258-
data_norm_unit = (data.unit * norm.unit).decompose().bases[0]
259-
if data_norm_unit != u.ct:
260-
raise ValueError('Data must be supplied in counts or the product of the norm and data has units of '
261-
'counts')
262-
data = data * norm
263-
self.norm = norm
264-
super().__init__(data, **kwargs)

sunkit_spex/spectrum/tests/__init__.py

Whitespace-only changes.

sunkit_spex/spectrum/tests/test_spectrum.py

Lines changed: 1 addition & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33

44
import astropy.units as u
55

6-
from sunkit_spex.spectrum.spectrum import CountSpectrum, Spectrum
6+
from sunkit_spex.spectrum.spectrum import Spectrum
77

88

99
def test_spectrum_bin_edges():
@@ -14,15 +14,3 @@ def test_spectrum_bin_edges():
1414
def test_spectrum_bin_centers():
1515
spec = Spectrum(np.arange(1, 11)*u.watt, spectral_axis=(np.arange(1, 11) - 0.5) * u.keV)
1616
assert_array_equal(spec._spectral_axis, [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5] * u.keV)
17-
18-
19-
def test_countspectrum():
20-
count_spec = CountSpectrum(np.arange(1, 11)*u.ct, spectral_axis=np.arange(1, 12)*u.keV)
21-
assert isinstance(count_spec, CountSpectrum)
22-
assert count_spec.norm is None
23-
24-
25-
def test_countspectrum_with_normalization():
26-
count_spec = CountSpectrum(np.arange(1, 11)*u.ct*u.s, spectral_axis=np.arange(1, 12)*u.keV, norm=np.ones(10)/u.s)
27-
assert isinstance(count_spec, CountSpectrum)
28-
assert_array_equal(count_spec.norm, np.ones(10)/u.s)
Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
import numpy as np
2+
from numpy.testing import assert_array_equal
3+
4+
from astropy.nddata import NDDataRef
5+
6+
from sunkit_spex.spectrum.uncertainty import PoissonUncertainty
7+
8+
9+
def test_add():
10+
data = np.array([0, 1, 2])
11+
uncert = np.sqrt(data)
12+
a = NDDataRef(data, uncertainty=PoissonUncertainty(uncert))
13+
b = NDDataRef(data.copy(), uncertainty=PoissonUncertainty(uncert.copy()))
14+
aplusb = a.add(b)
15+
assert_array_equal(aplusb.data, 2*data)
16+
assert_array_equal(aplusb.uncertainty.array, np.sqrt(2*uncert ** 2))
17+
18+
19+
def test_subtract():
20+
data = np.array([0, 1, 2])
21+
uncert = np.sqrt(data)
22+
a = NDDataRef(data, uncertainty=PoissonUncertainty(uncert))
23+
b = NDDataRef(data.copy(), uncertainty=PoissonUncertainty(uncert.copy()))
24+
aminusb = a.subtract(b)
25+
assert_array_equal(aminusb.data, data-data)
26+
assert_array_equal(aminusb.uncertainty.array, np.sqrt(2*uncert ** 2))
27+
28+
29+
def test_multiply():
30+
data = np.array([0, 1, 2])
31+
uncert = np.sqrt(data)
32+
a = NDDataRef(data, uncertainty=PoissonUncertainty(uncert))
33+
b = NDDataRef(data.copy(), uncertainty=PoissonUncertainty(uncert.copy()))
34+
atimesb = a.multiply(b)
35+
assert_array_equal(atimesb.data, data**2)
36+
assert_array_equal(atimesb.uncertainty.array, np.sqrt((2*data**2*uncert**2))) # (b**2*da**2 + a**2db**2)**0.5
37+
38+
39+
def test_divide():
40+
data = np.array([0, 1, 2])
41+
uncert = np.sqrt(data)
42+
a = NDDataRef(data, uncertainty=PoissonUncertainty(uncert))
43+
b = NDDataRef(data.copy(), uncertainty=PoissonUncertainty(uncert.copy()))
44+
adivb = a.divide(b)
45+
assert_array_equal(adivb.data, data/data)
46+
assert_array_equal(adivb.uncertainty.array, np.sqrt(((1/data)**2 * uncert**2)*2)) # ((1/b)**2*da**2 + (a/b**2)**2db**2)**0.5
Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
import numpy as np
2+
3+
from astropy.nddata.nduncertainty import (
4+
IncompatibleUncertaintiesException,
5+
NDUncertainty,
6+
VarianceUncertainty,
7+
_VariancePropagationMixin,
8+
)
9+
10+
__all__ = ["PoissonUncertainty"]
11+
12+
13+
class PoissonUncertainty(_VariancePropagationMixin, NDUncertainty):
14+
"""Poissonian uncertainty assuming first order error propagation.
15+
16+
This class implements uncertainty propagation for ``addition``,
17+
``subtraction``, ``multiplication`` and ``division`` with other instances
18+
of `Poisson`. The class can handle if the uncertainty has a
19+
unit that differs from (but is convertible to) the parents `NDData` unit.
20+
The unit of the resulting uncertainty will have the same unit as the
21+
resulting data. Also support for correlation is possible but requires the
22+
correlation as input. It cannot handle correlation determination itself.
23+
24+
Parameters
25+
----------
26+
args, kwargs :
27+
see `NDUncertainty`
28+
29+
Examples
30+
--------
31+
`Poisson` should always be associated with an `NDData`-like
32+
instance, either by creating it during initialization::
33+
34+
>>> from sunkit_spex.spectrum.uncertainty import PoissonUncertainty
35+
>>> ndd = NDData([1,2,3], unit='m',
36+
... uncertainty=PoissonUncertainty([0.1, 0.1, 0.1]))
37+
>>> ndd.uncertainty # doctest: +FLOAT_CMP
38+
PoissonUncertainty([0.1, 0.1, 0.1])
39+
40+
or by setting it manually on the `NDData` instance::
41+
42+
>>> ndd.uncertainty = PoissonUncertainty([0.2], unit='m', copy=True)
43+
>>> ndd.uncertainty # doctest: +FLOAT_CMP
44+
PoissonUncertainty([0.2])
45+
46+
the uncertainty ``array`` can also be set directly::
47+
48+
>>> ndd.uncertainty.array = 2
49+
>>> ndd.uncertainty
50+
PoissonUncertainty(2)
51+
52+
.. note::
53+
The unit will not be displayed.
54+
"""
55+
56+
@property
57+
def supports_correlated(self):
58+
"""`True` : `StdDevUncertainty` allows to propagate correlated \
59+
uncertainties.
60+
61+
``correlation`` must be given, this class does not implement computing
62+
it by itself.
63+
"""
64+
return True
65+
66+
@property
67+
def uncertainty_type(self):
68+
"""``"poisson"`` : `PoissonUncertainty` implements Poisson uncertainty."""
69+
return "poisson"
70+
71+
def _convert_uncertainty(self, other_uncert):
72+
if isinstance(other_uncert, PoissonUncertainty):
73+
return other_uncert
74+
else:
75+
raise IncompatibleUncertaintiesException
76+
77+
def _propagate_add(self, other_uncert, result_data, correlation):
78+
return super()._propagate_add_sub(
79+
other_uncert,
80+
result_data,
81+
correlation,
82+
subtract=False,
83+
to_variance=np.square,
84+
from_variance=np.sqrt,
85+
)
86+
87+
def _propagate_subtract(self, other_uncert, result_data, correlation):
88+
return super()._propagate_add_sub(
89+
other_uncert,
90+
result_data,
91+
correlation,
92+
subtract=True,
93+
to_variance=np.square,
94+
from_variance=np.sqrt,
95+
)
96+
97+
def _propagate_multiply(self, other_uncert, result_data, correlation):
98+
return super()._propagate_multiply_divide(
99+
other_uncert,
100+
result_data,
101+
correlation,
102+
divide=False,
103+
to_variance=np.square,
104+
from_variance=np.sqrt,
105+
)
106+
107+
def _propagate_divide(self, other_uncert, result_data, correlation):
108+
return super()._propagate_multiply_divide(
109+
other_uncert,
110+
result_data,
111+
correlation,
112+
divide=True,
113+
to_variance=np.square,
114+
from_variance=np.sqrt,
115+
)
116+
117+
def _propagate_collapse(self, numpy_operation, axis):
118+
# defer to _VariancePropagationMixin
119+
return super()._propagate_collapse(numpy_operation, axis=axis)
120+
121+
def _data_unit_to_uncertainty_unit(self, value):
122+
return value
123+
124+
def _convert_to_variance(self):
125+
new_array = None if self.array is None else self.array**2
126+
new_unit = None if self.unit is None else self.unit**2
127+
return VarianceUncertainty(new_array, unit=new_unit)
128+
129+
@classmethod
130+
def _convert_from_variance(cls, var_uncert):
131+
new_array = None if var_uncert.array is None else var_uncert.array ** (1 / 2)
132+
new_unit = None if var_uncert.unit is None else var_uncert.unit ** (1 / 2)
133+
return cls(new_array, unit=new_unit)

0 commit comments

Comments
 (0)