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

Jump to bottom

Template redshift enhancements #551

Merged
merged 12 commits into from
Nov 26, 2019
Merged

Template redshift enhancements #551

Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
a31c0cb
Merge pull request #1 from astropy/master
ibusko Nov 13, 2019
2cf9bb8
Redshift iterable.
Nov 13, 2019
f39854f
Rename variable.
ibusko Nov 14, 2019
76f5ebe
Rename variable.
ibusko Nov 14, 2019
ff48040
Rename variable.
ibusko Nov 14, 2019
ff68e70
All tests pass.
ibusko Nov 14, 2019
94b8bba
All tests pass.
ibusko Nov 14, 2019
2adc0f3
Output list with best-fit chi2 values, one per termplate.
ibusko Nov 14, 2019
19e2ec0
Output list with chi2 values, one per input redshift value.
ibusko Nov 14, 2019
8bb7890
Fix doc text.
Nov 14, 2019
0e643a3
Fix python code in doc text.
Nov 14, 2019
be05679
Output list with chi2 values.
Nov 22, 2019
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
19 changes: 14 additions & 5 deletions docs/analysis.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,8 +182,14 @@ This function will:
normalized to the observed spectrum (and the index of the template
spectrum if the list of templates is iterable).

If the redshift is unknown, the user specifies the expected window in redshift
(``min_redshift`` and ``max_redshift``) and the redshift spacing (``delta_redshift``).
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
values to use. As an example, a simple linear grid can be built with:

.. code-block:: python

>>> rs_values = np.arange(1., 3.25, 0.25)

The ~`specutils.analysis.template_comparison.template_match` function will then:

1. Move each template to the first term in the redshift grid (``min_redshift``).
Expand All @@ -192,8 +198,10 @@ The ~`specutils.analysis.template_comparison.template_match` function will then:
4. Run steps 1 and 2 of the case with known redshift.
5. Repeat the steps until ``max_redshift`` is reached.
6. Return the best redshift, the lowest chi-square and its corresponding
template spectrum, redshifted and normalized to the observed spectrum
(and the index of the template spectrum if the list of templates is iterable).
template spectrum, and a list with all chi2 values, one per template.
The returned template spectrum corresponding to the lowest chi2 is redshifted
and normalized to the observed spectrum (and the index of the template spectrum if
the list of templates is iterable).

An example of how to do template matching with an unknown redshift is:

Expand All @@ -206,10 +214,11 @@ An example of how to do template matching with an unknown redshift is:
>>> max_redshift = 3.0
>>> delta_redshift = .25
>>> resample_method = "flux_conserving"
>>> rs_values = np.arange(min_redshift, max_redshift+delta_redshift, delta_redshift)

>>> observed_spectrum = Spectrum1D(spectral_axis=spec_axis*(1+observed_redshift), flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
>>> spectral_template = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
>>> tm_result = template_comparison.template_match(observed_spectrum=observed_spectrum, spectral_templates=spectral_template, resample_method=resample_method, min_redshift=min_redshift, max_redshift=max_redshift, delta_redshift=delta_redshift) # doctest:+FLOAT_CMP
>>> tm_result = template_comparison.template_match(observed_spectrum=observed_spectrum, spectral_templates=spectral_template, resample_method=resample_method, redshift=rs_values) # doctest:+FLOAT_CMP


Reference/API
Expand Down
70 changes: 34 additions & 36 deletions specutils/analysis/template_comparison.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ def _chi_square_for_templates(observed_spectrum, template_spectrum, resample_met

def template_match(observed_spectrum, spectral_templates,
resample_method="flux_conserving", known_redshift=None,
min_redshift=None, max_redshift=None, delta_redshift=None):
redshift=None):
"""
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
Expand All @@ -135,12 +135,10 @@ def template_match(observed_spectrum, spectral_templates,
known_redshift: `float`
If the user knows the redshift they want to apply to the spectrum/spectra within spectral_templates,
then this redshift can be applied to each template before attempting the match.
min_redshift : `float`
The minimum redshift allowed.
max_redshift : `float`
The maximum redshift allowed.
delta_redshift : `float`
The amount the redshift will change between loops.
redshift : `list`, `tuple`, 'numpy.array`
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should merge redshift and known_redshift into a single parameter. We can then check if the value is a float or array-like and deal with it appropriately. Using this function, I wouldn't expect that a single value or a list of values would require two different arguments.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with this, but the one question is: is there some specific science case for known_redshift being used in addition to the redshift grid? I can't think of one offhand but just want to be sure... maybe @camipacifici has some insight?

An iterable with redshift values to be applied to each template, before computation of the
corresponding chi2 value. For each template, the redshift value that results in the smallest chi2
is used.

Returns
-------
Expand All @@ -151,17 +149,19 @@ def template_match(observed_spectrum, spectral_templates,
normalized template spectrum.
smallest_chi_index : `int`
The index of the spectrum with the smallest chi2 in spectral templates.
chi2_list : `list`
A list with the best-fit chi2 values found for each template spectrum.
"""
if hasattr(spectral_templates, 'flux') and len(spectral_templates.flux.shape) == 1:

# Account for redshift if provided
if all(x is not None for x in (min_redshift, max_redshift, delta_redshift)):
redshift, redshifted_spectrum = template_redshift(observed_spectrum, spectral_templates,
min_redshift, max_redshift, delta_redshift)
if redshift is not None and all(x is not None for x in redshift):
_, redshifted_spectrum, _ = template_redshift(observed_spectrum, spectral_templates,
redshift=redshift)
spectral_templates = redshifted_spectrum
elif known_redshift:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(This elif block can then be removed if we go with a single redshift argument that can be either a float or a list.)

redshift, redshifted_spectrum = template_redshift(observed_spectrum, spectral_templates,
known_redshift, known_redshift, 1.0)
_, redshifted_spectrum, _ = template_redshift(observed_spectrum, spectral_templates,
redshift=(known_redshift,))
spectral_templates = redshifted_spectrum

normalized_spectral_template, chi2 = _chi_square_for_templates(
Expand All @@ -175,32 +175,35 @@ def template_match(observed_spectrum, spectral_templates,
# chi square.
chi2_min = None
smallest_chi_spec = None
chi2_list = []

for index, spectrum in enumerate(spectral_templates):

# Account for redshift if provided
if all(x is not None for x in (min_redshift, max_redshift, delta_redshift)):
redshift, redshifted_spectrum = template_redshift(observed_spectrum, spectrum,
min_redshift, max_redshift, delta_redshift)
if redshift is not None and all(x is not None for x in redshift):
_, redshifted_spectrum, _ = template_redshift(observed_spectrum, spectrum,
redshift=redshift)
spectrum = redshifted_spectrum

elif known_redshift:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Again, can be removed.)

redshift, redshifted_spectrum = template_redshift(observed_spectrum, spectrum,
known_redshift, known_redshift, 1.0)
_, redshifted_spectrum, _ = template_redshift(observed_spectrum, spectrum,
redshift=(known_redshift,))
spectrum = redshifted_spectrum

normalized_spectral_template, chi2 = _chi_square_for_templates(
observed_spectrum, spectrum, resample_method)

chi2_list.append(chi2)

if chi2_min is None or chi2 < chi2_min:
chi2_min = chi2
smallest_chi_spec = normalized_spectral_template
smallest_chi_index = index

return smallest_chi_spec, chi2_min, smallest_chi_index
return smallest_chi_spec, chi2_min, smallest_chi_index, chi2_list


def template_redshift(observed_spectrum, template_spectrum, min_redshift, max_redshift, delta_redshift):
def template_redshift(observed_spectrum, template_spectrum, redshift):
"""
Find the best-fit redshift for template_spectrum to match observed_spectrum using chi2.

Expand All @@ -210,12 +213,8 @@ def template_redshift(observed_spectrum, template_spectrum, min_redshift, max_re
The observed spectrum.
template_spectrum : :class:`~specutils.Spectrum1D`
The template spectrum, which will have it's redshift calculated.
min_redshift : `float`
The minimum redshift allowed.
max_redshift : `float`
The maximum redshift allowed.
delta_redshift : `float`
The amount the redshift will change between loops.
redshift : `list`, `tuple`, 'numpy.array`
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
redshift : `list`, `tuple`, 'numpy.array`
redshift : `float`, `list`, `tuple`, 'numpy.array`

An iterable with the redshift values to test.

Returns
-------
Expand All @@ -224,29 +223,28 @@ def template_redshift(observed_spectrum, template_spectrum, min_redshift, max_re
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.
chi2_list : `list`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this should be an ndarray instead of a list. That will make it a lot easier to work with if it's multidimensional.

Relatedly, I'm having a bit of trouble following the logic: is the resulting array is of dimension equal to the spectral template set, or is it a 2D array that's n_template x n_redshift ? I think the latter is the eventual goal, but if we're not quite there it could be a follow-on.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In some out-of-band discussion, @ibusko said he'd write a test for the SpectrumCollection case and see what happens.

My preference is that it returns a 2D array as indicated above (or the transpose)

A list with the chi2 values corresponding to each input redshift value.
"""
if min_redshift > max_redshift:
raise ValueError("The `max_redshift` value must be greater than `min_redshift`.")
return

redshift = min_redshift
chi2_min = None
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
chi2_min = None
redshift = np.array(redshift).reshape((np.array(redshift).size,))
chi2_min = None

final_redshift = None
chi2_list = []

# Loop which goes through available redshift options and finds the smallest chi2
while redshift <= max_redshift:
# Loop which goes through available redshift values and finds the smallest chi2
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+redshift),
redshifted_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis*(1+rs),
flux=template_spectrum.flux, uncertainty=template_spectrum.uncertainty,
meta=template_spectrum.meta)
normalized_spectral_template, chi2 = _chi_square_for_templates(
observed_spectrum, redshifted_spectrum, "flux_conserving")

chi2_list.append(chi2)

# Set new chi2_min if suitable replacement is found
if not np.isnan(chi2) and (chi2_min is None or chi2 < chi2_min):
chi2_min = chi2
final_redshift = redshift
redshift += delta_redshift
final_redshift = rs

return final_redshift, redshifted_spectrum
return final_redshift, redshifted_spectrum, chi2_list
20 changes: 12 additions & 8 deletions specutils/tests/test_template_comparison.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,6 +154,10 @@ def test_template_match_list():

np.testing.assert_almost_equal(tm_result[1], 40093.28353756253)

# make sure that multiple template spectra will create a list of
# chi2 values, one per template.
assert len(tm_result) == 4


def test_template_match_spectrum_collection():
"""
Expand Down Expand Up @@ -232,11 +236,12 @@ def test_template_unknown_redshift():
min_redshift = .5
max_redshift = 5.5
delta_redshift = .25
redshift_trial_values = np.arange(min_redshift, max_redshift, delta_redshift)

tr_result = template_comparison.template_redshift(observed_spectrum=spec, template_spectrum=spec1,
min_redshift=min_redshift, max_redshift=max_redshift,
delta_redshift=delta_redshift)
redshift=redshift_trial_values)

assert len(tr_result) == 3
assert tr_result[0] == 3

def test_template_redshift_with_one_template_spectrum_in_match():
Expand Down Expand Up @@ -264,11 +269,11 @@ def test_template_redshift_with_one_template_spectrum_in_match():
min_redshift = .5
max_redshift = 5.5
delta_redshift = .25
redshift_trial_values = np.arange(min_redshift, max_redshift+delta_redshift, delta_redshift)

tm_result = template_comparison.template_match(observed_spectrum=spec, spectral_templates=spec1,
resample_method="flux_conserving", known_redshift=None,
min_redshift=min_redshift, max_redshift=max_redshift,
delta_redshift=delta_redshift)
redshift=redshift_trial_values)

np.testing.assert_almost_equal(tm_result[1], 73484.0053895151)

Expand Down Expand Up @@ -303,11 +308,11 @@ def test_template_redshift_with_multiple_template_spectra_in_match():
min_redshift = .5
max_redshift = 5.5
delta_redshift = .25
redshift_trial_values = np.arange(min_redshift, max_redshift+delta_redshift, delta_redshift)

tm_result = template_comparison.template_match(observed_spectrum=spec, spectral_templates=spec_coll,
resample_method="flux_conserving", known_redshift=None,
min_redshift=min_redshift, max_redshift=max_redshift,
delta_redshift=delta_redshift)
redshift=redshift_trial_values)

np.testing.assert_almost_equal(tm_result[1], 6803.922741644725)

Expand Down Expand Up @@ -336,6 +341,5 @@ def test_template_known_redshift():

tm_result = template_comparison.template_match(observed_spectrum=spec, spectral_templates=spec1,
resample_method="flux_conserving", known_redshift=redshift,
min_redshift=None, max_redshift=None,
delta_redshift=None)
redshift=None)
np.testing.assert_almost_equal(tm_result[1], 1.9062409482056814e-31)