Skip to content

Commit

Permalink
Truncate in template_comparison by default rather than nan_fill (#1121)
Browse files Browse the repository at this point in the history
* Working on adding truncate option as default for template comparison, need diff strategy than None fill

* Trying different strategy for truncate, just drop all NaNs

* Working on resolving test failures, added warning for no overlap and return None instead of blank Spectrum1D

* Warn in redshifted case if there's no overlap with any redshift

* Resolve test failures, changed returned template to resampled version.

* Resolve remaining test failures

* Fix typos

* Fix bug, improving test coverage

* Return full normalized+redshifted template
rosteen authored Feb 9, 2024
1 parent 17a6f6b commit 425eab3
Showing 6 changed files with 151 additions and 76 deletions.
2 changes: 1 addition & 1 deletion docs/analysis.rst
Original file line number Diff line number Diff line change
@@ -217,7 +217,7 @@ This function will:
2. Compute the chi-square between the observed spectrum and each template.
3. Return the lowest chi-square and its corresponding template spectrum,
normalized to the observed spectrum (and the index of the template
spectrum if the list of templates is iterable). It also
spectrum if the list of templates is iterable).

If the redshift is unknown, the user specifies a grid of redshift values in the
form of an iterable object such as a list, tuple, or numpy array with the redshift
13 changes: 8 additions & 5 deletions docs/manipulation.rst
Original file line number Diff line number Diff line change
@@ -41,9 +41,9 @@ implemented are: :func:`~specutils.manipulation.box_smooth`
:func:`~specutils.manipulation.gaussian_smooth`
(:class:`~astropy.convolution.Gaussian1DKernel`), and
:func:`~specutils.manipulation.trapezoid_smooth`
(:class:`~astropy.convolution.Trapezoid1DKernel`). Note that, although
(:class:`~astropy.convolution.Trapezoid1DKernel`). Note that, although
these kernels are 1D, they can be applied to higher-dimensional
data (e.g. spectral cubes), in which case the data will be smoothed only
data (e.g. spectral cubes), in which case the data will be smoothed only
along the spectral dimension.

.. code-block:: python
@@ -102,8 +102,8 @@ that takes the spectrum and an astropy 1D kernel. So, one could also do:
43., 44., 45., 46., 47., 48., 49.] nm>)>
In this case, the ``spec1_bsmooth2`` result should be equivalent to the ``spec1_bsmooth`` in
the section above (assuming the flux data of the input ``spec`` is the same). Note that,
as in the case of the kernel-specific functions, a 1D kernel can be applied to a
the section above (assuming the flux data of the input ``spec`` is the same). Note that,
as in the case of the kernel-specific functions, a 1D kernel can be applied to a
multi-dimensional spectrum and will smooth that spectrum along the spectral dimension.
In the case of :func:`~specutils.manipulation.convolution_smooth`, one can also input
a higher-dimensional kernel that matches the dimensionality of the data.
@@ -156,6 +156,9 @@ Each of these classes takes in a :class:`~specutils.Spectrum1D` and a user
defined output dispersion grid, and returns a new :class:`~specutils.Spectrum1D`
with the resampled flux. Currently the resampling classes expect the new
dispersion grid unit to be the same as the input spectrum's dispersion grid unit.
Additionally, all resamplers take an optional ``extrapolation_treatment`` keyword which
can be ``nan_fill``, ``zero_fill``, or ``truncate``, to determine what to do with output
wavelength bins that have no overlap with the original spectrum.

If the input :class:`~specutils.Spectrum1D` contains an uncertainty,
:class:`~specutils.manipulation.FluxConservingResampler` will propogate the
@@ -402,7 +405,7 @@ with the spline knots:
>>> result
<Spectrum1D(flux=<Quantity [ 2., 4., 6., 8., 10., 12., 14., 16., 18., 20.] mJy>, spectral_axis=<SpectralAxis [ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.] Angstrom>)>
The default behavior is to keep the data outside the replaced region unchanged.
The default behavior is to keep the data outside the replaced region unchanged.
Alternatively, the spectrum outside the replaced region can be filled with zeros:

.. code-block:: python
93 changes: 65 additions & 28 deletions specutils/analysis/template_comparison.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
import warnings

import numpy as np
from astropy.nddata import StdDevUncertainty

@@ -30,6 +32,7 @@ def _normalize_for_template_matching(observed_spectrum, template_spectrum, stdde
"""
if stddev is None:
stddev = observed_spectrum.uncertainty.represent_as(StdDevUncertainty).quantity

num = np.nansum((observed_spectrum.flux*template_spectrum.flux) / (stddev**2))
# We need to limit this sum to where observed_spectrum is not NaN as well.
template_filtered = ((template_spectrum.flux / stddev)**2)
@@ -39,7 +42,7 @@ def _normalize_for_template_matching(observed_spectrum, template_spectrum, stdde
return num/denom


def _resample(resample_method):
def _resample(resample_method, extrapolation_treatment):
"""
Find the user preferred method of resampling the template spectrum to fit
the observed spectrum.
@@ -55,18 +58,19 @@ def _resample(resample_method):
This is the actual class that will handle the resampling.
"""
if resample_method == "flux_conserving":
return FluxConservingResampler()
return FluxConservingResampler(extrapolation_treatment=extrapolation_treatment)

if resample_method == "linear_interpolated":
return LinearInterpolatedResampler()
return LinearInterpolatedResampler(extrapolation_treatment=extrapolation_treatment)

if resample_method == "spline_interpolated":
return SplineInterpolatedResampler()
return SplineInterpolatedResampler(extrapolation_treatment=extrapolation_treatment)

return None


def _chi_square_for_templates(observed_spectrum, template_spectrum, resample_method):
def _chi_square_for_templates(observed_spectrum, template_spectrum, resample_method,
extrapolation_treatment, warn_no_overlap=True):
"""
Resample the template spectrum to match the wavelength of the observed
spectrum. Then, calculate chi2 on the flux of the two spectra.
@@ -88,25 +92,36 @@ def _chi_square_for_templates(observed_spectrum, template_spectrum, resample_met
normalized template spectrum.
"""
# Resample template
if _resample(resample_method) != 0:
fluxc_resample = _resample(resample_method)
if _resample(resample_method, extrapolation_treatment) != 0:
fluxc_resample = _resample(resample_method, extrapolation_treatment)
template_obswavelength = fluxc_resample(template_spectrum,
observed_spectrum.spectral_axis)

# With truncate, template may be smaller than observed spectrum if they don't fully overlap
matching_indices = np.intersect1d(observed_spectrum.spectral_axis,
template_obswavelength.spectral_axis,
return_indices=True)[1]
if len(matching_indices) == 0:
if warn_no_overlap:
warnings.warn("Template spectrum has no overlap with observed spectrum.")
return None, np.nan

observed_truncated = observed_spectrum[matching_indices.min():matching_indices.max()+1]

# Convert the uncertainty to standard deviation if needed
stddev = observed_spectrum.uncertainty.represent_as(StdDevUncertainty).array
stddev = observed_truncated.uncertainty.represent_as(StdDevUncertainty).array

# Normalize spectra
normalization = _normalize_for_template_matching(observed_spectrum,
normalization = _normalize_for_template_matching(observed_truncated,
template_obswavelength,
stddev)

# Numerator
num_right = normalization * template_obswavelength.flux
num = observed_spectrum.flux - num_right
num = observed_truncated.flux - num_right

# Denominator
denom = stddev * observed_spectrum.flux.unit
denom = stddev * observed_truncated.flux.unit

# Get chi square
result = (num/denom)**2
@@ -121,9 +136,8 @@ def _chi_square_for_templates(observed_spectrum, template_spectrum, resample_met
return normalized_template_spectrum, chi2


def template_match(observed_spectrum, spectral_templates,
resample_method="flux_conserving",
redshift=None):
def template_match(observed_spectrum, spectral_templates, redshift=None,
resample_method="flux_conserving", extrapolation_treatment="truncate"):
"""
Find which spectral templates is the best fit to an observed spectrum by
computing the chi-squared. If two template_spectra have the same chi2, the
@@ -138,20 +152,26 @@ def template_match(observed_spectrum, spectral_templates,
over. The template spectra, which will be resampled, normalized, and
compared to the observed spectrum, where the smallest chi2 and
normalized template spectrum will be returned.
resample_method : `string`
Three resample options: flux_conserving, linear_interpolated, and spline_interpolated.
Anything else does not resample the spectrum.
redshift : 'float', `int`, `list`, `tuple`, 'numpy.array`
If the user knows the redshift they want to apply to the spectrum/spectra within spectral_templates,
then this float or int value redshift can be applied to each template before attempting the match.
Or, alternatively, an iterable with redshift values to be applied to each template, before computation
of the corresponding chi2 value, can be passed via this same parameter. For each template, the redshift
value that results in the smallest chi2 is used.
resample_method : `string`
Three resample options: flux_conserving, linear_interpolated, and spline_interpolated.
Anything else does not resample the spectrum.
extrapolation_treatment : `string`
Three options for what to do if the template spectrum and observed spectrum have regions
that do not overlap: ``nan_fill`` and ``zero_fill`` to fill the non-overlapping bins,
or ``truncate`` to remove the non-overlapping bins. Passed to the resampler, defaults to
``truncate``.
Returns
-------
normalized_template_spectrum : :class:`~specutils.Spectrum1D`
The template spectrum that has been normalized.
The template spectrum that has been normalized and resampled to the observed
spectrum's spectral axis.
redshift : `None` or `float`
The value of the redshift that provides the smallest chi2.
smallest_chi_index : `int`
@@ -174,7 +194,7 @@ def template_match(observed_spectrum, spectral_templates,

else:
normalized_spectral_template, chi2 = _chi_square_for_templates(
observed_spectrum, spectral_templates, resample_method)
observed_spectrum, spectral_templates, resample_method, extrapolation_treatment)

chi2_list.append(chi2)

@@ -202,7 +222,7 @@ def template_match(observed_spectrum, spectral_templates,

else:
normalized_spectral_template, chi2 = _chi_square_for_templates(
observed_spectrum, spectrum, resample_method)
observed_spectrum, spectrum, resample_method, extrapolation_treatment)

if chi2_min is None or chi2 < chi2_min:
chi2_min = chi2
@@ -214,7 +234,8 @@ def template_match(observed_spectrum, spectral_templates,
return smallest_chi_spec, final_redshift, smallest_chi_index, chi2_min, chi2_list


def template_redshift(observed_spectrum, template_spectrum, redshift):
def template_redshift(observed_spectrum, template_spectrum, redshift,
resample_method="flux_conserving", extrapolation_treatment="truncate"):
"""
Find the best-fit redshift for template_spectrum to match observed_spectrum using chi2.
@@ -226,12 +247,22 @@ def template_redshift(observed_spectrum, template_spectrum, redshift):
The template spectrum, which will have it's redshift calculated.
redshift : `float`, `int`, `list`, `tuple`, 'numpy.array`
A scalar or iterable with the redshift values to test.
resample_method : `string`
Three resample options: flux_conserving, linear_interpolated, and spline_interpolated.
Anything else does not resample the spectrum.
extrapolation_treatment : `string`
Three options for what to do if the template spectrum and observed spectrum have regions
that do not overlap: ``nan_fill`` and ``zero_fill`` to fill the non-overlapping bins,
or ``truncate`` to remove the non-overlapping bins. Passed to the resampler, defaults to
``truncate``.
Returns
-------
redshifted_spectrum: :class:`~specutils.Spectrum1D`
A new Spectrum1D object which incorporates the template_spectrum with a spectral_axis
that has been redshifted using the final_redshift.
that has been redshifted using the final_redshift, and resampled to that of the
observed spectrum.
final_redshift : `float`
The best-fit redshift for template_spectrum to match the observed_spectrum.
normalized_template_spectrum : :class:`~specutils.Spectrum1D`
@@ -243,6 +274,7 @@ def template_redshift(observed_spectrum, template_spectrum, redshift):
"""
chi2_min = None
final_redshift = None
final_spectrum = None
chi2_list = []

redshift = np.array(redshift).reshape((np.array(redshift).size,))
@@ -251,13 +283,14 @@ def template_redshift(observed_spectrum, template_spectrum, redshift):
for rs in redshift:

# Create new redshifted spectrum and run it through the chi2 method
redshifted_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis*(1+rs),
flux=template_spectrum.flux,
uncertainty=template_spectrum.uncertainty,
meta=template_spectrum.meta)
redshifted_spectrum = template_spectrum._copy()
redshifted_spectrum.shift_spectrum_to(redshift=rs)

normalized_spectral_template, chi2 = _chi_square_for_templates(
observed_spectrum, redshifted_spectrum, "flux_conserving")
normalized_spectral_template, chi2 = _chi_square_for_templates(observed_spectrum,
redshifted_spectrum,
resample_method,
extrapolation_treatment,
warn_no_overlap=False)

chi2_list.append(chi2)

@@ -267,4 +300,8 @@ def template_redshift(observed_spectrum, template_spectrum, redshift):
final_redshift = rs
final_spectrum = redshifted_spectrum

if final_spectrum is None:
warnings.warn("Template spectrum and observed spectrum have no overlap after"
"applying specified redshift(s) to template.")

return final_spectrum, final_redshift, normalized_spectral_template, chi2_min, chi2_list
70 changes: 50 additions & 20 deletions specutils/manipulation/resample.py
Original file line number Diff line number Diff line change
@@ -20,11 +20,12 @@ class ResamplerBase(ABC):
----------
extrapolation_treatment : str
What to do when resampling off the edge of the spectrum. Can be
``'nan_fill'`` to have points beyond the edges by set to NaN, or
``'zero_fill'`` to be set to zero.
``'nan_fill'`` to have points beyond the edges by set to NaN,
``'zero_fill'`` to set thoe points to zero, or ``'truncate'`` to
truncate any non-overlapping bins of the spectrum.
"""
def __init__(self, extrapolation_treatment='nan_fill'):
if extrapolation_treatment not in ('nan_fill', 'zero_fill'):
if extrapolation_treatment not in ('nan_fill', 'zero_fill', 'truncate'):
raise ValueError('invalid extrapolation_treatment value: ' + str(extrapolation_treatment))
self.extrapolation_treatment = extrapolation_treatment

@@ -55,8 +56,9 @@ class FluxConservingResampler(ResamplerBase):
----------
extrapolation_treatment : str
What to do when resampling off the edge of the spectrum. Can be
``'nan_fill'`` to have points beyond the edges by set to NaN, or
``'zero_fill'`` to be set to zero.
``'nan_fill'`` to have points beyond the edges by set to NaN,
``'zero_fill'`` to set those points to zero, or ``'truncate'`` to
truncate any non-overlapping bins of the spectrum.
Examples
--------
@@ -74,7 +76,7 @@ class FluxConservingResampler(ResamplerBase):
>>> resample_grid = [1, 5, 9, 13, 14, 17, 21, 22, 23] *u.nm
>>> fluxc_resample = FluxConservingResampler()
>>> fluxc_resample(input_spectra, resample_grid) # doctest: +FLOAT_CMP
<Spectrum1D(flux=<Quantity [ nan, 3. , 6. , 7. , 6.25, 10. , 20. , nan, nan] mJy>, spectral_axis=<SpectralAxis [ 1., 5., 9., 13., 14., 17., 21., 22., 23.] nm>)>
<Spectrum1D(flux=<Quantity [ 1. , 3. , 6. , 7. , 6.25, 10. , 20. , nan, nan] mJy>, spectral_axis=<SpectralAxis [ 1., 5., 9., 13., 14., 17., 21., 22., 23.] nm>)>
"""

def _fluxc_resample(self, input_bin_centers, output_bin_centers,
@@ -122,9 +124,9 @@ def _fluxc_resample(self, input_bin_centers, output_bin_centers,
min_idx = 0
max_idx = None

low_out_of_range = np.where(output_bin_edges < input_bin_edges[0])[0]
low_out_of_range = np.where(output_bin_edges <= input_bin_edges[0])[0]
if len(low_out_of_range) > 0: # if any bins below wavelength range
min_idx = low_out_of_range[-1] + 1
min_idx = low_out_of_range[-1] # This doesn't need +1 because bin_edges has len+1 compared to output_fluxes
output_fluxes[:min_idx] = fill_val
if errs is not None:
output_errs[:min_idx] = fill_val
@@ -289,6 +291,12 @@ def resample1d(self, orig_spectrum, fin_spec_axis):
output_fluxes = output_fluxes << orig_spectrum.flux.unit
fin_spec_axis = np.array(fin_spec_axis) << orig_spectrum.spectral_axis.unit

if self.extrapolation_treatment == 'truncate':
fin_spec_axis = fin_spec_axis[np.where(~np.isnan(output_fluxes))]
if new_errs is not None:
new_errs = new_errs[np.where(~np.isnan(output_fluxes))]
output_fluxes = output_fluxes[np.where(~np.isnan(output_fluxes))]

resampled_spectrum = Spectrum1D(flux=output_fluxes,
spectral_axis=fin_spec_axis,
uncertainty=new_errs)
@@ -304,8 +312,9 @@ class LinearInterpolatedResampler(ResamplerBase):
----------
extrapolation_treatment : str
What to do when resampling off the edge of the spectrum. Can be
``'nan_fill'`` to have points beyond the edges by set to NaN, or
``'zero_fill'`` to be set to zero.
``'nan_fill'`` to have points beyond the edges by set to NaN,
``'zero_fill'`` to set those points to zero, or ``'truncate'`` to
truncate any non-overlapping bins of the spectrum.
Examples
--------
@@ -365,6 +374,12 @@ def resample1d(self, orig_spectrum, fin_spec_axis):
new_unc = orig_spectrum.uncertainty.__class__(array=out_unc_arr,
unit=orig_spectrum.unit)

if self.extrapolation_treatment == 'truncate':
fin_spec_axis = fin_spec_axis[np.where(~np.isnan(out_flux))]
if new_unc is not None:
new_unc = new_unc[np.where(~np.isnan(out_flux))]
out_flux = out_flux[np.where(~np.isnan(out_flux))]

return Spectrum1D(spectral_axis=fin_spec_axis,
flux=out_flux,
uncertainty=new_unc)
@@ -380,8 +395,10 @@ class SplineInterpolatedResampler(ResamplerBase):
----------
extrapolation_treatment : str
What to do when resampling off the edge of the spectrum. Can be
``'nan_fill'`` to have points beyond the edges by set to NaN, or
``'zero_fill'`` to be set to zero.
``'nan_fill'`` to have points beyond the edges by set to NaN,
``'zero_fill'`` to set those points to zero, or ``'truncate'`` to
truncate any non-overlapping bins of the spectrum. Any other value will
have the spline interpolate beyond the edges of the original data.
Examples
--------
@@ -403,8 +420,8 @@ class SplineInterpolatedResampler(ResamplerBase):
7.29736328, nan, nan, nan] mJy>, spectral_axis=<SpectralAxis [ 1., 5., 9., 13., 14., 17., 21., 22., 23.] nm>)>
"""
def __init__(self, bin_edges='nan_fill'):
super().__init__(bin_edges)
def __init__(self, extrapolation_treatment='nan_fill'):
super().__init__(extrapolation_treatment)

def resample1d(self, orig_spectrum, fin_spec_axis):
"""
@@ -425,22 +442,35 @@ def resample1d(self, orig_spectrum, fin_spec_axis):
"""
orig_axis_in_new = orig_spectrum.spectral_axis.to(fin_spec_axis.unit)
flux_spline = CubicSpline(orig_axis_in_new.value, orig_spectrum.flux.value,
extrapolate=self.extrapolation_treatment != 'nan_fill')
extrapolate=self.extrapolation_treatment not in ('nan_fill',
'zero_fill',
'truncate'))
out_flux_val = flux_spline(fin_spec_axis.value)

new_unc = None
if orig_spectrum.uncertainty is not None:
unc_spline = CubicSpline(orig_axis_in_new.value, orig_spectrum.uncertainty.array,
extrapolate=self.extrapolation_treatment != 'nan_fill')
extrapolate=self.extrapolation_treatment not in ('nan_fill',
'zero_fill',
'truncate'))
out_unc_val = unc_spline(fin_spec_axis.value)
new_unc = orig_spectrum.uncertainty.__class__(array=out_unc_val, unit=orig_spectrum.unit)

fill_val = np.nan
if self.extrapolation_treatment == 'zero_fill':
origedges = orig_spectrum.spectral_axis.bin_edges
off_edges = (fin_spec_axis < origedges[0]) | (origedges[-1] < fin_spec_axis)
out_flux_val[off_edges] = 0
fill_val = 0

origedges = orig_spectrum.spectral_axis.bin_edges
off_edges = (fin_spec_axis < origedges[0]) | (origedges[-1] < fin_spec_axis)
out_flux_val[off_edges] = fill_val
if new_unc is not None:
new_unc.array[off_edges] = fill_val

if self.extrapolation_treatment == 'truncate':
fin_spec_axis = fin_spec_axis[np.where(~np.isnan(out_flux_val))]
if new_unc is not None:
new_unc.array[off_edges] = 0
new_unc = new_unc[np.where(~np.isnan(out_flux_val))]
out_flux_val = out_flux_val[np.where(~np.isnan(out_flux_val))]

return Spectrum1D(spectral_axis=fin_spec_axis,
flux=out_flux_val*orig_spectrum.flux.unit,
16 changes: 15 additions & 1 deletion specutils/tests/test_resample.py
Original file line number Diff line number Diff line change
@@ -48,7 +48,7 @@ def test_expanded_grid_fluxconserving():
results = inst(input_spectra, resamp_grid)

assert_quantity_allclose(results.flux,
np.array([np.nan, 3., 6., 7., 6.25, 10., 20.,
np.array([1., 3., 6., 7., 6.25, 10., 20.,
np.nan, np.nan])*u.mJy)


@@ -150,6 +150,12 @@ def test_expanded_grid_interp_linear():
np.array([np.nan, 3.5, 5.5, 6.75, 6.5, 9.5, np.nan,
np.nan, np.nan])*u.mJy)

inst = LinearInterpolatedResampler(extrapolation_treatment="truncate")
results = inst(input_spectra, resamp_grid)

assert_quantity_allclose(results.flux,
np.array([3.5, 5.5, 6.75, 6.5, 9.5])*u.mJy)


def test_expanded_grid_interp_spline():
"""
@@ -169,6 +175,14 @@ def test_expanded_grid_interp_spline():
5.89921875, 7.29736328, np.nan, np.nan,
np.nan])*u.mJy)

inst = SplineInterpolatedResampler(extrapolation_treatment="truncate")
input_spectra = Spectrum1D(spectral_axis=wave_val * u.AA, flux=flux_val * u.mJy)
results = inst(input_spectra, resamp_grid)

assert_quantity_allclose(results.flux,
np.array([3.98808594, 6.94042969, 6.45869141,
5.89921875, 7.29736328])*u.mJy)


@pytest.mark.parametrize("edgetype,lastvalue",
[("nan_fill", np.nan), ("zero_fill", 0)])
33 changes: 12 additions & 21 deletions specutils/tests/test_template_comparison.py
Original file line number Diff line number Diff line change
@@ -18,7 +18,7 @@ def test_template_match_no_overlap():

# Create test spectra
spec_axis = np.linspace(0, 50, 50) * u.AA
spec_axis_no_overlap = np.linspace(51, 102, 50) * u.AA
spec_axis_no_overlap = np.linspace(52, 103, 50) * u.AA

spec = Spectrum1D(spectral_axis=spec_axis,
flux=np.random.randn(50) * u.Jy,
@@ -29,16 +29,12 @@ def test_template_match_no_overlap():
uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))

# Get result from template_match
tm_result = template_comparison.template_match(spec, spec1)
assert tm_result[3] == 0.0
with pytest.warns(UserWarning, match="Template spectrum has no overlap with observed spectrum"):
tm_result = template_comparison.template_match(spec, spec1)

# Create new spectrum for comparison
spec_result = Spectrum1D(spectral_axis=spec_axis,
flux=spec1.flux * template_comparison._normalize_for_template_matching(spec, spec1))
try:
assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
except AssertionError:
pytest.xfail('TODO: investigate why this is failing')
assert np.isnan(tm_result[3])

assert tm_result[0] is None


def test_template_match_minimal_overlap():
@@ -56,24 +52,19 @@ def test_template_match_minimal_overlap():
spec_axis_min_overlap[0] = 51.0 * u.AA

spec = Spectrum1D(spectral_axis=spec_axis,
flux=np.random.randn(50) * u.Jy,
flux=abs(np.random.randn(50)) * u.Jy,
uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))

spec1 = Spectrum1D(spectral_axis=spec_axis_min_overlap,
flux=np.random.randn(50) * u.Jy,
flux=abs(np.random.randn(50)) * u.Jy,
uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))

# Get result from template_match
tm_result = template_comparison.template_match(spec, spec1)
assert tm_result[3] == 0.0
np.testing.assert_almost_equal(tm_result[3], 0)

# Create new spectrum for comparison
spec_result = Spectrum1D(spectral_axis=spec_axis,
flux=spec1.flux * template_comparison._normalize_for_template_matching(spec, spec1))
try:
assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
except AssertionError:
pytest.xfail('TODO: investigate why all elements in tm_result[0] are NaN even with overlap')
assert len(tm_result[0].flux) == 50
assert quantity_allclose(tm_result[0].flux, spec1.flux*4.17, atol=0.01*u.Jy)


def test_template_match_spectrum():
@@ -332,7 +323,7 @@ def test_template_redshift_with_multiple_template_spectra_in_match():
# TODO: Determine cause of and fix failing assert
# np.testing.assert_almost_equal(tm_result[1], redshift)

np.testing.assert_almost_equal(tm_result[3], 1172.135458895792)
np.testing.assert_almost_equal(tm_result[3], 1172.1446112796834)

# When a spectrum collection is matched with a redshift
# grid, a list-of-lists is returned with the trial chi2

0 comments on commit 425eab3

Please sign in to comment.