Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix moment calculations #1141

Merged
merged 5 commits into from
May 29, 2024
Merged

Fix moment calculations #1141

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
fd720e9
Correcting moment 0 calculation. Still have one failing test
rosteen May 24, 2024
8086104
This seems inelegant but does fix the failing test
rosteen May 24, 2024
93ea26c
Update doc example output
rosteen May 24, 2024
9b095f2
Update math for higher orders as well - cancels out for constant bin …
rosteen May 28, 2024
4d6c35f
Fix remote data test failure in docs
rosteen May 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/analysis.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,7 @@ The `~specutils.analysis.moment` function computes moments of any order:

>>> from specutils.analysis import moment
>>> moment(noisy_gaussian, SpectralRegion(7*u.GHz, 3*u.GHz)) # doctest:+FLOAT_CMP
<Quantity 98.26318995 Jy>
<Quantity 4.93784874 GHz Jy>
>>> moment(noisy_gaussian, SpectralRegion(7*u.GHz, 3*u.GHz), order=1) # doctest:+FLOAT_CMP
<Quantity 4.99909151 GHz>
>>> moment(noisy_gaussian, SpectralRegion(7*u.GHz, 3*u.GHz), order=2) # doctest:+FLOAT_CMP
Expand Down
6 changes: 3 additions & 3 deletions docs/spectral_cube.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,9 +115,9 @@ along the spectral axis (remember that the spectral axis is always last in a
>>> m.shape # doctest: +REMOTE_DATA
(74, 74)
>>> m[30:33,30:33] # doctest: +REMOTE_DATA +FLOAT_CMP
<Quantity [[6.45261317e-07, 6.46265069e-07, 6.48128166e-07],
[6.46467930e-07, 6.47941283e-07, 6.51460998e-07],
[6.48672775e-07, 6.52631872e-07, 6.56733087e-07]] m>
<Quantity [[6.97933331e-07, 6.98926463e-07, 7.00540974e-07],
[6.98959625e-07, 7.00280655e-07, 7.03511823e-07],
[7.00740294e-07, 7.04527986e-07, 7.08245958e-07]] m>

Use Case
========
Expand Down
17 changes: 11 additions & 6 deletions specutils/analysis/moment.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

import numpy as np
from ..manipulation import extract_region
from ..spectra import SpectrumCollection
from .utils import computation_wrapper


Expand Down Expand Up @@ -39,6 +40,9 @@ def moment(spectrum, regions=None, order=0, axis=-1):
Moment of the spectrum. Returns None if (order < 0 or None)

"""
if isinstance(spectrum, SpectrumCollection):
return [computation_wrapper(_compute_moment, spec, regions,order=order, axis=axis)
for spec in spectrum]
return computation_wrapper(_compute_moment, spectrum, regions,
order=order, axis=axis)

Expand All @@ -60,24 +64,25 @@ def _compute_moment(spectrum, regions=None, order=0, axis=-1):
if order is None or order < 0:
return None

dx = np.abs(np.diff(spectral_axis.bin_edges))
m0 = np.sum(flux * dx, axis=axis)
if order == 0:
return np.sum(flux, axis=axis)
return m0

dispersion = spectral_axis
if len(flux.shape) > len(spectral_axis.shape):
_shape = flux.shape[:-1] + (1,)
dispersion = np.tile(spectral_axis, _shape)

if order == 1:
return np.sum(flux * dispersion, axis=axis) / np.sum(flux, axis=axis)
return np.sum(flux * dispersion * dx, axis=axis) / m0

if order > 1:
m0 = np.sum(flux, axis=axis)

# By setting keepdims to True, the axes which are reduced are
# left in the result as dimensions with size one. This means
# that we can broadcast m1 correctly against dispersion.
m1 = (np.sum(flux * spectral_axis, axis=axis, keepdims=True)
/ np.sum(flux, axis=axis, keepdims=True))
m1 = (np.sum(flux * dispersion * dx, axis=axis, keepdims=True)
/ np.sum(flux * dx, axis=axis, keepdims=True))

return np.sum(flux * (dispersion - m1) ** order, axis=axis) / m0
return np.sum(flux * dx * (dispersion - m1) ** order, axis=axis) / m0
13 changes: 9 additions & 4 deletions specutils/tests/test_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1065,8 +1065,8 @@ def test_moment():
spectrum = Spectrum1D(spectral_axis=frequencies, flux=flux)

moment_0 = moment(spectrum, order=0)
assert moment_0.unit.is_equivalent(u.Jy)
assert quantity_allclose(moment_0, 252.96*u.Jy, atol=0.01*u.Jy)
assert moment_0.unit.is_equivalent(u.Jy * u.GHz)
assert quantity_allclose(moment_0, 2.5045*u.Jy*u.GHz, atol=0.001*u.Jy*u.GHz)

moment_1 = moment(spectrum, order=1)
assert moment_1.unit.is_equivalent(u.GHz)
Expand Down Expand Up @@ -1097,6 +1097,11 @@ def test_moment_cube():

spectrum = Spectrum1D(spectral_axis=frequencies, flux=flux_multid)

moment_0 = moment(spectrum, order=0)

assert moment_0.shape == (9, 10)
assert moment_0.unit.is_equivalent(u.Jy*u.GHz)

moment_1 = moment(spectrum, order=1)

assert moment_1.shape == (9, 10)
Expand Down Expand Up @@ -1159,8 +1164,8 @@ def test_moment_cube_order_2():
assert moment_2.shape == (10, 10000)
assert moment_2.unit.is_equivalent(u.GHz**2)
# check assorted values.
assert quantity_allclose(moment_2[0][0], 2.019e-28*u.GHz**2, rtol=0.01)
assert quantity_allclose(moment_2[1][0], 2.019e-28*u.GHz**2, rtol=0.01)
assert quantity_allclose(moment_2[0][0], 8.078e-28*u.GHz**2, rtol=0.01)
assert quantity_allclose(moment_2[1][0], 8.078e-28*u.GHz**2, rtol=0.01)
assert quantity_allclose(moment_2[0][3], 2.019e-28*u.GHz**2, rtol=0.01)


Expand Down
Loading