From 9b67821c567fbaa9a4bfc1bed40a7edbe21385ff Mon Sep 17 00:00:00 2001
From: Sara Ogaz <ogaz@stsci.edu>
Date: Tue, 16 Jul 2019 17:02:52 -0400
Subject: [PATCH 01/31] add examples to class file for help docs

---
 specutils/manipulation/resample.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/specutils/manipulation/resample.py b/specutils/manipulation/resample.py
index d5a113139..cf4fcb933 100644
--- a/specutils/manipulation/resample.py
+++ b/specutils/manipulation/resample.py
@@ -88,6 +88,7 @@ class FluxConservingResampler(ResamplerBase):
     >>> fluxc_resample = FluxConservingResampler()
     >>> output_spectrum1D = fluxc_resample(input_spectra, resample_grid) # doctest: +IGNORE_OUTPUT
 
+
     """
 
     def __call__(self, orig_spectrum, fin_lamb):
@@ -253,6 +254,7 @@ class LinearInterpolatedResampler(ResamplerBase):
     >>> resample_grid = np.array([1, 5, 9, 13, 14, 17, 21, 22, 23])
     >>> fluxc_resample = LinearInterpolatedResampler()
     >>> output_spectrum1D = fluxc_resample(input_spectra, resample_grid) # doctest: +IGNORE_OUTPUT
+
     """
     def __init__(self, bin_edges='nan_fill'):
         super().__init__(bin_edges)

From 6acc5430999a5c98c5e96671ca0841f07c958927 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Wed, 17 Jul 2019 13:57:37 -0400
Subject: [PATCH 02/31] template matching minus SpectrumCollection
 implementation

---
 specutils/analysis/template_match.py | 36 ++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)
 create mode 100644 specutils/analysis/template_match.py

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
new file mode 100644
index 000000000..0afb3fcf4
--- /dev/null
+++ b/specutils/analysis/template_match.py
@@ -0,0 +1,36 @@
+from ..spectra.spectrum1d import Spectrum1D
+from ..spectra.spectrum_collection import SpectrumCollection
+from ..manipulation import FluxConservingResample
+from scipy.stats import chisquare
+
+def _template_match(observed_spectrum, template_spectrum):
+    # resample the template spectrum to match the wavelength of the observed spectrum
+    fluxc_resample = FluxConservingResample()
+    output_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
+
+    # calculate chi square
+    x2 = chisquare(observed_spectrum.flux, output_spectrum1D.flux)
+    return x2
+
+def template_match(observed_spectrum, collection):
+    """
+    Find what instance collection is and run _template_match accordingly
+    """
+    if isinstance(collection, Spectrum1D):
+        return Spectrum1D, _template_match(observed_spectrum, collection)
+
+    elif isinstance(collection, SpectrumCollection):
+        pass
+
+    # Loop through spectra in list and return spectrum with lowest chi square
+    # and its corresponding chi square
+    elif isinstance(collection, list):
+        min = None
+        smallest_chi_spec = None
+        for spectrum in collection:
+            chi = _template_match(observed_spectrum, spectrum)
+            if min is None or chi < min:
+                min = chi
+                smallest_chi_spec = spectrum
+
+        return smallest_chi_spec, min

From 7805915f164fa59a2378fd9a94b14739475a854e Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Thu, 18 Jul 2019 10:11:32 -0400
Subject: [PATCH 03/31] improved variable names and documentation

---
 specutils/analysis/template_match.py | 36 +++++++++++++++++-----------
 1 file changed, 22 insertions(+), 14 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 0afb3fcf4..5a72e8c98 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -4,33 +4,41 @@
 from scipy.stats import chisquare
 
 def _template_match(observed_spectrum, template_spectrum):
-    # resample the template spectrum to match the wavelength of the observed spectrum
+    """
+    Resample the template spectrum to match the wavelength of the observed spectrum.
+    Then, run scipy.stats.chisquare on the flux of the two spectra.
+
+    :param observed_spectrum: The observed spectrum
+    :param template_spectrum: The template spectrum, which will be resampled to match the wavelength of the observed
+        spectrum
+    :return: chi square of the flux of the observed spectrum and the flux of the template spectrum
+    """
     fluxc_resample = FluxConservingResample()
-    output_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
+    template_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
 
     # calculate chi square
-    x2 = chisquare(observed_spectrum.flux, output_spectrum1D.flux)
+    x2 = chisquare(observed_spectrum.flux, template_spectrum1D.flux)
     return x2
 
-def template_match(observed_spectrum, collection):
+def template_match(observed_spectrum, spectral_templates):
     """
     Find what instance collection is and run _template_match accordingly
     """
-    if isinstance(collection, Spectrum1D):
-        return Spectrum1D, _template_match(observed_spectrum, collection)
+    if isinstance(spectral_templates, Spectrum1D):
+        return Spectrum1D, _template_match(observed_spectrum, spectral_templates)
 
-    elif isinstance(collection, SpectrumCollection):
+    elif isinstance(spectral_templates, SpectrumCollection):
         pass
 
     # Loop through spectra in list and return spectrum with lowest chi square
     # and its corresponding chi square
-    elif isinstance(collection, list):
-        min = None
+    elif isinstance(spectral_templates, list):
+        chi2_min = None
         smallest_chi_spec = None
-        for spectrum in collection:
-            chi = _template_match(observed_spectrum, spectrum)
-            if min is None or chi < min:
-                min = chi
+        for spectrum in spectral_templates:
+            chi2 = _template_match(observed_spectrum, spectrum)
+            if chi2_min is None or chi2 < chi2_min:
+                chi2_min = chi2
                 smallest_chi_spec = spectrum
 
-        return smallest_chi_spec, min
+        return smallest_chi_spec, chi2_min

From 13ca915c4ec042dff0f95a3c3697e2218c313e8c Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Thu, 25 Jul 2019 08:53:05 -0400
Subject: [PATCH 04/31] problems with s1d.uncertainty. rebasing then creating
 tests

---
 specutils/analysis/template_match.py | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 5a72e8c98..8846933bb 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -2,6 +2,15 @@
 from ..spectra.spectrum_collection import SpectrumCollection
 from ..manipulation import FluxConservingResample
 from scipy.stats import chisquare
+import numpy as np
+import math
+
+def _normalize(observed_spectrum, template_spectrum):
+    print(observed_spectrum.uncertainty)
+
+    normalization = math.pow(np.sum((observed_spectrum*template_spectrum)/(observed_spectrum.uncertainty)), 2) \
+                    / math.pow(np.sum((template_spectrum/observed_spectrum.uncertainty)), 2)
+    return normalization
 
 def _template_match(observed_spectrum, template_spectrum):
     """
@@ -17,8 +26,12 @@ def _template_match(observed_spectrum, template_spectrum):
     template_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
 
     # calculate chi square
-    x2 = chisquare(observed_spectrum.flux, template_spectrum1D.flux)
-    return x2
+    # x2 = chisquare(observed_spectrum.flux, template_spectrum1D.flux)
+
+    normalization = _normalize(observed_spectrum, template_spectrum1D)
+
+    chi2 = np.sum(((observed_spectrum.flux-normalization*template_spectrum1D.flux)/observed_spectrum.uncertainty)**2 )
+    return chi2
 
 def template_match(observed_spectrum, spectral_templates):
     """

From 52f582ebc5ebb274827dcc7a54b4f238b76d1dc7 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 26 Jul 2019 13:43:28 -0400
Subject: [PATCH 05/31] template_match complete with tests

---
 specutils/analysis/template_match.py   |  39 ++++---
 specutils/tests/test_template_match.py | 151 +++++++++++++++++++++++++
 2 files changed, 174 insertions(+), 16 deletions(-)
 create mode 100644 specutils/tests/test_template_match.py

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 8846933bb..6c29f6146 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -1,16 +1,17 @@
 from ..spectra.spectrum1d import Spectrum1D
 from ..spectra.spectrum_collection import SpectrumCollection
-from ..manipulation import FluxConservingResample
-from scipy.stats import chisquare
+from ..manipulation import FluxConservingResampler
 import numpy as np
-import math
 
 def _normalize(observed_spectrum, template_spectrum):
-    print(observed_spectrum.uncertainty)
+    """
+    The normalization is necessary to bring the template to the same magnitude as the observation and minimize the chi^2
+    """
+    num = np.sum((observed_spectrum.flux*template_spectrum.flux)/(observed_spectrum.uncertainty.array**2))
+    denom = np.sum((template_spectrum.flux/observed_spectrum.uncertainty.array)**2)
+    normalized = num/denom
 
-    normalization = math.pow(np.sum((observed_spectrum*template_spectrum)/(observed_spectrum.uncertainty)), 2) \
-                    / math.pow(np.sum((template_spectrum/observed_spectrum.uncertainty)), 2)
-    return normalization
+    return normalized
 
 def _template_match(observed_spectrum, template_spectrum):
     """
@@ -22,15 +23,24 @@ def _template_match(observed_spectrum, template_spectrum):
         spectrum
     :return: chi square of the flux of the observed spectrum and the flux of the template spectrum
     """
-    fluxc_resample = FluxConservingResample()
+    # Resample template
+    fluxc_resample = FluxConservingResampler()
     template_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
 
-    # calculate chi square
-    # x2 = chisquare(observed_spectrum.flux, template_spectrum1D.flux)
-
+    # Normalize spectra
     normalization = _normalize(observed_spectrum, template_spectrum1D)
 
-    chi2 = np.sum(((observed_spectrum.flux-normalization*template_spectrum1D.flux)/observed_spectrum.uncertainty)**2 )
+    # Numerator
+    num_right = normalization*template_spectrum1D.flux
+    num = observed_spectrum.flux - num_right
+
+    # Denominator
+    denom = observed_spectrum.uncertainty.array
+
+    # Get chi square
+    result = (num/denom)**2
+    chi2 = np.sum(result)
+
     return chi2
 
 def template_match(observed_spectrum, spectral_templates):
@@ -40,12 +50,9 @@ def template_match(observed_spectrum, spectral_templates):
     if isinstance(spectral_templates, Spectrum1D):
         return Spectrum1D, _template_match(observed_spectrum, spectral_templates)
 
-    elif isinstance(spectral_templates, SpectrumCollection):
-        pass
-
     # Loop through spectra in list and return spectrum with lowest chi square
     # and its corresponding chi square
-    elif isinstance(spectral_templates, list):
+    elif isinstance(spectral_templates, list) or isinstance(spectral_templates, SpectrumCollection):
         chi2_min = None
         smallest_chi_spec = None
         for spectrum in spectral_templates:
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
new file mode 100644
index 000000000..0d17dffd3
--- /dev/null
+++ b/specutils/tests/test_template_match.py
@@ -0,0 +1,151 @@
+import astropy.units as u
+import numpy as np
+from astropy.nddata import StdDevUncertainty
+
+from ..spectra.spectrum1d import Spectrum1D
+from ..spectra.spectrum_collection import SpectrumCollection
+from ..analysis import template_match
+from ..manipulation import FluxConservingResampler
+
+
+def _normalize(observed_spectrum, template_spectrum):
+    num = np.sum((observed_spectrum.flux*template_spectrum.flux)/(observed_spectrum.uncertainty.array**2))
+    denom = np.sum((template_spectrum.flux/observed_spectrum.uncertainty.array)**2)
+    normalized = num/denom
+
+    return normalized
+
+def get_chi_square(observed_spectrum, template_spectrum):
+    # Resample template
+    fluxc_resample = FluxConservingResampler()
+    template_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
+
+    # Normalize spectra
+    normalization = _normalize(observed_spectrum, template_spectrum1D)
+
+    # Numerator
+    num_right = normalization*template_spectrum1D.flux
+    num = observed_spectrum.flux - num_right
+
+    # Denominator
+    denom = observed_spectrum.uncertainty.array
+
+    # Get chi square
+    result = (num/denom)**2
+    chi2 = np.sum(result)
+
+    return chi2
+
+def test_template_match_spectrum():
+    """
+    Test template_match when both observed and template spectra have the same wavelength axis
+    """
+    # Create test spectra
+    spec_axis = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    # Get result from template_match
+    tm_result = template_match.template_match(spec, spec1)
+
+    # Calculate result independently of template_match
+    result = get_chi_square(spec, spec1)
+
+    assert result == tm_result[1]
+
+def test_template_match_with_resample():
+    """
+    Test template_match when both observed and template spectra have different wavelength axis using resampling
+    """
+    # Create test spectra
+    spec_axis1 = np.linspace(0, 50, 50) * u.AA
+    spec_axis2 = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis1,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    # Get result from template_match
+    tm_result = template_match.template_match(spec, spec1)
+
+    result = get_chi_square(spec, spec1)
+
+    assert result == tm_result[1]
+
+def test_template_match_list():
+    """
+    Test template_match when template spectra are in a list
+    """
+    # Create test spectra
+    spec_axis1 = np.linspace(0, 50, 50) * u.AA
+    spec_axis2 = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis1,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+    spec2 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    # Combine spectra into list
+    template_list = [spec1, spec2]
+
+    # Get result from template_match
+    tm_result = template_match.template_match(spec, template_list)
+
+    # Loop through list in order to find out which template spectra has the smallest chi square
+    chi2_min = None
+    smallest_chi_spec = None
+    for spectrum in template_list:
+        chi2 = get_chi_square(spec, spectrum)
+        if chi2_min is None or chi2 < chi2_min:
+            chi2_min = chi2
+            smallest_chi_spec = spectrum
+
+
+    assert chi2_min == tm_result[1]
+    assert smallest_chi_spec == tm_result[0]
+
+def test_template_match_spectrum_collection():
+    """
+    Test template_match when template spectra are in a SpectrumCollection object
+    """
+    # Create test spectra
+    spec_axis1 = np.linspace(0, 50, 50) * u.AA
+    spec_axis2 = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis1,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+    spec2 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    # Combine spectra into SpectrumCollection object
+    spec_coll = SpectrumCollection.from_spectra([spec1, spec2])
+
+    # Get result from template_match
+    tm_result = template_match.template_match(spec, spec_coll)
+
+    # Loop through SpectrumCollection in order to find out which template spectra has the smallest chi square
+    chi2_min = None
+    for spectrum in spec_coll:
+        chi2 = get_chi_square(spec, spectrum)
+        if chi2_min is None or chi2 < chi2_min:
+            chi2_min = chi2
+
+    assert chi2_min == tm_result[1]

From 9377ed1615778e66fda1f65229ad0496cd4172eb Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 26 Jul 2019 13:45:40 -0400
Subject: [PATCH 06/31] removed spacing in resample.py

---
 specutils/manipulation/resample.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/specutils/manipulation/resample.py b/specutils/manipulation/resample.py
index cf4fcb933..d5a113139 100644
--- a/specutils/manipulation/resample.py
+++ b/specutils/manipulation/resample.py
@@ -88,7 +88,6 @@ class FluxConservingResampler(ResamplerBase):
     >>> fluxc_resample = FluxConservingResampler()
     >>> output_spectrum1D = fluxc_resample(input_spectra, resample_grid) # doctest: +IGNORE_OUTPUT
 
-
     """
 
     def __call__(self, orig_spectrum, fin_lamb):
@@ -254,7 +253,6 @@ class LinearInterpolatedResampler(ResamplerBase):
     >>> resample_grid = np.array([1, 5, 9, 13, 14, 17, 21, 22, 23])
     >>> fluxc_resample = LinearInterpolatedResampler()
     >>> output_spectrum1D = fluxc_resample(input_spectra, resample_grid) # doctest: +IGNORE_OUTPUT
-
     """
     def __init__(self, bin_edges='nan_fill'):
         super().__init__(bin_edges)

From 2825913333f2a781e700490c247c7f53b9fcf929 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Thu, 8 Aug 2019 11:35:17 -0400
Subject: [PATCH 07/31] chi2 return no longer returns units

---
 specutils/analysis/template_match.py   | 4 ++--
 specutils/tests/test_template_match.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 6c29f6146..de19a0269 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -16,7 +16,7 @@ def _normalize(observed_spectrum, template_spectrum):
 def _template_match(observed_spectrum, template_spectrum):
     """
     Resample the template spectrum to match the wavelength of the observed spectrum.
-    Then, run scipy.stats.chisquare on the flux of the two spectra.
+    Then, run chisquare on the flux of the two spectra.
 
     :param observed_spectrum: The observed spectrum
     :param template_spectrum: The template spectrum, which will be resampled to match the wavelength of the observed
@@ -35,7 +35,7 @@ def _template_match(observed_spectrum, template_spectrum):
     num = observed_spectrum.flux - num_right
 
     # Denominator
-    denom = observed_spectrum.uncertainty.array
+    denom = observed_spectrum.uncertainty.array * observed_spectrum.flux.unit
 
     # Get chi square
     result = (num/denom)**2
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index 0d17dffd3..be9f9e1f0 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -28,7 +28,7 @@ def get_chi_square(observed_spectrum, template_spectrum):
     num = observed_spectrum.flux - num_right
 
     # Denominator
-    denom = observed_spectrum.uncertainty.array
+    denom = observed_spectrum.uncertainty.array * observed_spectrum.flux.unit
 
     # Get chi square
     result = (num/denom)**2

From c63f222b1ceda5636bb5d8da391d12ea78d3e4cd Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Thu, 8 Aug 2019 12:14:55 -0400
Subject: [PATCH 08/31] return normalized template spectrum and associated chi2

---
 specutils/analysis/template_match.py   | 13 +++++++++----
 specutils/tests/test_template_match.py | 27 ++++++++++++++++++--------
 2 files changed, 28 insertions(+), 12 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index de19a0269..ebce0d7ec 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -41,24 +41,29 @@ def _template_match(observed_spectrum, template_spectrum):
     result = (num/denom)**2
     chi2 = np.sum(result)
 
-    return chi2
+    normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
+                                              spectral_axis=template_spectrum.spectral_axis)
+
+    return chi2, normalized_template_spectrum
 
 def template_match(observed_spectrum, spectral_templates):
     """
     Find what instance collection is and run _template_match accordingly
     """
     if isinstance(spectral_templates, Spectrum1D):
-        return Spectrum1D, _template_match(observed_spectrum, spectral_templates)
+        chi2, normalized_spectral_template = _template_match(observed_spectrum, spectral_templates)
+        return normalized_spectral_template, chi2
 
     # Loop through spectra in list and return spectrum with lowest chi square
     # and its corresponding chi square
     elif isinstance(spectral_templates, list) or isinstance(spectral_templates, SpectrumCollection):
         chi2_min = None
         smallest_chi_spec = None
+
         for spectrum in spectral_templates:
-            chi2 = _template_match(observed_spectrum, spectrum)
+            chi2, normalized_spectral_template = _template_match(observed_spectrum, spectrum)
             if chi2_min is None or chi2 < chi2_min:
                 chi2_min = chi2
-                smallest_chi_spec = spectrum
+                smallest_chi_spec = normalized_spectral_template
 
         return smallest_chi_spec, chi2_min
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index be9f9e1f0..a151a023f 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -16,6 +16,15 @@ def _normalize(observed_spectrum, template_spectrum):
     return normalized
 
 def get_chi_square(observed_spectrum, template_spectrum):
+    """
+    Resample the template spectrum to match the wavelength of the observed spectrum.
+    Then, run chisquare on the flux of the two spectra.
+
+    :param observed_spectrum: The observed spectrum
+    :param template_spectrum: The template spectrum, which will be resampled to match the wavelength of the observed
+        spectrum
+    :return: chi square of the flux of the observed spectrum and the flux of the template spectrum
+    """
     # Resample template
     fluxc_resample = FluxConservingResampler()
     template_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
@@ -34,7 +43,10 @@ def get_chi_square(observed_spectrum, template_spectrum):
     result = (num/denom)**2
     chi2 = np.sum(result)
 
-    return chi2
+    normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
+                                              spectral_axis=template_spectrum.spectral_axis)
+
+    return chi2, normalized_template_spectrum
 
 def test_template_match_spectrum():
     """
@@ -56,7 +68,7 @@ def test_template_match_spectrum():
     # Calculate result independently of template_match
     result = get_chi_square(spec, spec1)
 
-    assert result == tm_result[1]
+    assert result[0] == tm_result[1]
 
 def test_template_match_with_resample():
     """
@@ -78,7 +90,7 @@ def test_template_match_with_resample():
 
     result = get_chi_square(spec, spec1)
 
-    assert result == tm_result[1]
+    assert result[0] == tm_result[1]
 
 def test_template_match_list():
     """
@@ -106,16 +118,14 @@ def test_template_match_list():
 
     # Loop through list in order to find out which template spectra has the smallest chi square
     chi2_min = None
-    smallest_chi_spec = None
+
     for spectrum in template_list:
-        chi2 = get_chi_square(spec, spectrum)
+        chi2, normalized_spectral_template = get_chi_square(spec, spectrum)
         if chi2_min is None or chi2 < chi2_min:
             chi2_min = chi2
-            smallest_chi_spec = spectrum
 
 
     assert chi2_min == tm_result[1]
-    assert smallest_chi_spec == tm_result[0]
 
 def test_template_match_spectrum_collection():
     """
@@ -143,8 +153,9 @@ def test_template_match_spectrum_collection():
 
     # Loop through SpectrumCollection in order to find out which template spectra has the smallest chi square
     chi2_min = None
+
     for spectrum in spec_coll:
-        chi2 = get_chi_square(spec, spectrum)
+        chi2, normalized_template_spectrum = get_chi_square(spec, spectrum)
         if chi2_min is None or chi2 < chi2_min:
             chi2_min = chi2
 

From fb86dd894638cb2361804327e9ec51cb1e627b51 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Mon, 26 Aug 2019 11:38:12 -0400
Subject: [PATCH 09/31] rearranged returns to make more sense

---
 specutils/analysis/template_match.py   |  6 +++---
 specutils/tests/test_template_match.py | 10 +++++-----
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index ebce0d7ec..678ae4d55 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -44,14 +44,14 @@ def _template_match(observed_spectrum, template_spectrum):
     normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
                                               spectral_axis=template_spectrum.spectral_axis)
 
-    return chi2, normalized_template_spectrum
+    return normalized_template_spectrum, chi2
 
 def template_match(observed_spectrum, spectral_templates):
     """
     Find what instance collection is and run _template_match accordingly
     """
     if isinstance(spectral_templates, Spectrum1D):
-        chi2, normalized_spectral_template = _template_match(observed_spectrum, spectral_templates)
+        normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectral_templates)
         return normalized_spectral_template, chi2
 
     # Loop through spectra in list and return spectrum with lowest chi square
@@ -61,7 +61,7 @@ def template_match(observed_spectrum, spectral_templates):
         smallest_chi_spec = None
 
         for spectrum in spectral_templates:
-            chi2, normalized_spectral_template = _template_match(observed_spectrum, spectrum)
+            normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectrum)
             if chi2_min is None or chi2 < chi2_min:
                 chi2_min = chi2
                 smallest_chi_spec = normalized_spectral_template
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index a151a023f..3e76f6b0a 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -46,7 +46,7 @@ def get_chi_square(observed_spectrum, template_spectrum):
     normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
                                               spectral_axis=template_spectrum.spectral_axis)
 
-    return chi2, normalized_template_spectrum
+    return normalized_template_spectrum, chi2
 
 def test_template_match_spectrum():
     """
@@ -68,7 +68,7 @@ def test_template_match_spectrum():
     # Calculate result independently of template_match
     result = get_chi_square(spec, spec1)
 
-    assert result[0] == tm_result[1]
+    assert result[1] == tm_result[1]
 
 def test_template_match_with_resample():
     """
@@ -90,7 +90,7 @@ def test_template_match_with_resample():
 
     result = get_chi_square(spec, spec1)
 
-    assert result[0] == tm_result[1]
+    assert result[1] == tm_result[1]
 
 def test_template_match_list():
     """
@@ -120,7 +120,7 @@ def test_template_match_list():
     chi2_min = None
 
     for spectrum in template_list:
-        chi2, normalized_spectral_template = get_chi_square(spec, spectrum)
+        normalized_spectral_template, chi2 = get_chi_square(spec, spectrum)
         if chi2_min is None or chi2 < chi2_min:
             chi2_min = chi2
 
@@ -155,7 +155,7 @@ def test_template_match_spectrum_collection():
     chi2_min = None
 
     for spectrum in spec_coll:
-        chi2, normalized_template_spectrum = get_chi_square(spec, spectrum)
+        normalized_template_spectrum, chi2 = get_chi_square(spec, spectrum)
         if chi2_min is None or chi2 < chi2_min:
             chi2_min = chi2
 

From 7ed6e4ef0c5e7f2dd83774b0374117803a8f7876 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Mon, 26 Aug 2019 11:47:29 -0400
Subject: [PATCH 10/31] continued cleanup

---
 specutils/analysis/template_match.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 678ae4d55..2d1606c99 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -41,8 +41,8 @@ def _template_match(observed_spectrum, template_spectrum):
     result = (num/denom)**2
     chi2 = np.sum(result)
 
-    normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
-                                              spectral_axis=template_spectrum.spectral_axis)
+    normalized_template_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis,
+                                              flux=template_spectrum.flux*normalization)
 
     return normalized_template_spectrum, chi2
 

From 12fc0e39ce8dd3caf34664bbb5f2f88e8a867b16 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Mon, 26 Aug 2019 11:54:11 -0400
Subject: [PATCH 11/31] reverting code because test failed

---
 specutils/analysis/template_match.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 2d1606c99..678ae4d55 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -41,8 +41,8 @@ def _template_match(observed_spectrum, template_spectrum):
     result = (num/denom)**2
     chi2 = np.sum(result)
 
-    normalized_template_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis,
-                                              flux=template_spectrum.flux*normalization)
+    normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
+                                              spectral_axis=template_spectrum.spectral_axis)
 
     return normalized_template_spectrum, chi2
 

From dd4414db2134ce17a52f6ce1bcd2d95bff306764 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Mon, 26 Aug 2019 11:59:26 -0400
Subject: [PATCH 12/31] added comments

---
 specutils/analysis/template_match.py   | 5 +++--
 specutils/tests/test_template_match.py | 1 +
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 678ae4d55..07978d2f2 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -41,8 +41,9 @@ def _template_match(observed_spectrum, template_spectrum):
     result = (num/denom)**2
     chi2 = np.sum(result)
 
-    normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
-                                              spectral_axis=template_spectrum.spectral_axis)
+    # Create normalized template spectrum, which will be returned with corresponding chi2
+    normalized_template_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis,
+                                              flux=template_spectrum.flux*normalization)
 
     return normalized_template_spectrum, chi2
 
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index 3e76f6b0a..89e14c916 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -43,6 +43,7 @@ def get_chi_square(observed_spectrum, template_spectrum):
     result = (num/denom)**2
     chi2 = np.sum(result)
 
+    # Create normalized template spectrum, which will be returned with corresponding chi2
     normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
                                               spectral_axis=template_spectrum.spectral_axis)
 

From 643f987b3308236f468d78075d3a3c09efdc393a Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Mon, 26 Aug 2019 12:04:53 -0400
Subject: [PATCH 13/31] still trying to pass test

---
 specutils/tests/test_template_match.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index 89e14c916..b83aad48a 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -44,8 +44,8 @@ def get_chi_square(observed_spectrum, template_spectrum):
     chi2 = np.sum(result)
 
     # Create normalized template spectrum, which will be returned with corresponding chi2
-    normalized_template_spectrum = Spectrum1D(flux=template_spectrum.flux*normalization,
-                                              spectral_axis=template_spectrum.spectral_axis)
+    normalized_template_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis,
+                                              flux=template_spectrum.flux*normalization)
 
     return normalized_template_spectrum, chi2
 

From 849e2d775eb0697148163ec22702433a54a16b7c Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Tue, 27 Aug 2019 22:35:42 -0400
Subject: [PATCH 14/31] added comments. improved code. improved tests

---
 specutils/analysis/template_match.py   | 61 +++++++++++++++-----
 specutils/tests/test_template_match.py | 80 ++++----------------------
 2 files changed, 60 insertions(+), 81 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 07978d2f2..63dfe7176 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -5,33 +5,52 @@
 
 def _normalize(observed_spectrum, template_spectrum):
     """
-    The normalization is necessary to bring the template to the same magnitude as the observation and minimize the chi^2
+    Calculate a scale factor to be applied to the template spectrum so the total flux in both spectra will be the same.
+
+    Parameters
+    ----------
+    observed_spectrum : `Spectrum1D`
+        the observed spectrum
+    template_spectrum : `Spectrum1D`
+        the template spectrum, which needs to be normalized in order to be compared with the observed_spectrum
+
+    Returns
+    -------
+    `float`
+        A float which will normalize the template spectrum's flux so that it can be compared to the observed spectrum
     """
     num = np.sum((observed_spectrum.flux*template_spectrum.flux)/(observed_spectrum.uncertainty.array**2))
     denom = np.sum((template_spectrum.flux/observed_spectrum.uncertainty.array)**2)
-    normalized = num/denom
-
-    return normalized
+    return num/denom
 
 def _template_match(observed_spectrum, template_spectrum):
     """
     Resample the template spectrum to match the wavelength of the observed spectrum.
-    Then, run chisquare on the flux of the two spectra.
+    Then, calculate chi2 on the flux of the two spectra.
 
-    :param observed_spectrum: The observed spectrum
-    :param template_spectrum: The template spectrum, which will be resampled to match the wavelength of the observed
-        spectrum
-    :return: chi square of the flux of the observed spectrum and the flux of the template spectrum
+    Parameters
+    ----------
+    observed_spectrum : `Spectrum1D`
+        the observed spectrum
+    template_spectrum : `Spectrum1D`
+        the template spectrum, which will be resampled to match the wavelength of the observed_spectrum
+
+    Returns
+    -------
+    normalized_template_spectrum : `Spectrum1D`
+        the template_spectrum that has been normalized
+    chi2 : `float`
+        the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
     """
     # Resample template
     fluxc_resample = FluxConservingResampler()
-    template_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
+    template_obswavelength = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
 
     # Normalize spectra
-    normalization = _normalize(observed_spectrum, template_spectrum1D)
+    normalization = _normalize(observed_spectrum, template_obswavelength)
 
     # Numerator
-    num_right = normalization*template_spectrum1D.flux
+    num_right = normalization*template_obswavelength.flux
     num = observed_spectrum.flux - num_right
 
     # Denominator
@@ -49,7 +68,23 @@ def _template_match(observed_spectrum, template_spectrum):
 
 def template_match(observed_spectrum, spectral_templates):
     """
-    Find what instance collection is and run _template_match accordingly
+    Find what instance spectral_templates is and run _template_match accordingly. If two template_spectra have the same
+    chi2, the first template is returned
+
+    Parameters
+    ----------
+    observed_spectrum : `Spectrum1D`
+        the observed spectrum
+    spectral_templates : `Spectrum1D` or `SpectrumCollection` or `list`
+        the template spectra, which will be resampled and normalized and compared to the observed_spectrum, where the
+        smallest chi2 and normalized_template_spectrum will be returned
+
+    Returns
+    -------
+    normalized_template_spectrum : `Spectrum1D`
+        the template_spectrum that has been normalized
+    chi2 : `float`
+        the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
     """
     if isinstance(spectral_templates, Spectrum1D):
         normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectral_templates)
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index b83aad48a..272a16fae 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -5,54 +5,14 @@
 from ..spectra.spectrum1d import Spectrum1D
 from ..spectra.spectrum_collection import SpectrumCollection
 from ..analysis import template_match
-from ..manipulation import FluxConservingResampler
 
 
-def _normalize(observed_spectrum, template_spectrum):
-    num = np.sum((observed_spectrum.flux*template_spectrum.flux)/(observed_spectrum.uncertainty.array**2))
-    denom = np.sum((template_spectrum.flux/observed_spectrum.uncertainty.array)**2)
-    normalized = num/denom
-
-    return normalized
-
-def get_chi_square(observed_spectrum, template_spectrum):
-    """
-    Resample the template spectrum to match the wavelength of the observed spectrum.
-    Then, run chisquare on the flux of the two spectra.
-
-    :param observed_spectrum: The observed spectrum
-    :param template_spectrum: The template spectrum, which will be resampled to match the wavelength of the observed
-        spectrum
-    :return: chi square of the flux of the observed spectrum and the flux of the template spectrum
-    """
-    # Resample template
-    fluxc_resample = FluxConservingResampler()
-    template_spectrum1D = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
-
-    # Normalize spectra
-    normalization = _normalize(observed_spectrum, template_spectrum1D)
-
-    # Numerator
-    num_right = normalization*template_spectrum1D.flux
-    num = observed_spectrum.flux - num_right
-
-    # Denominator
-    denom = observed_spectrum.uncertainty.array * observed_spectrum.flux.unit
-
-    # Get chi square
-    result = (num/denom)**2
-    chi2 = np.sum(result)
-
-    # Create normalized template spectrum, which will be returned with corresponding chi2
-    normalized_template_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis,
-                                              flux=template_spectrum.flux*normalization)
-
-    return normalized_template_spectrum, chi2
-
 def test_template_match_spectrum():
     """
     Test template_match when both observed and template spectra have the same wavelength axis
     """
+    np.random.seed(42)
+
     # Create test spectra
     spec_axis = np.linspace(0, 50, 50) * u.AA
     spec = Spectrum1D(spectral_axis=spec_axis,
@@ -66,15 +26,14 @@ def test_template_match_spectrum():
     # Get result from template_match
     tm_result = template_match.template_match(spec, spec1)
 
-    # Calculate result independently of template_match
-    result = get_chi_square(spec, spec1)
-
-    assert result[1] == tm_result[1]
+    assert tm_result[1] == 40093.28353756253
 
 def test_template_match_with_resample():
     """
     Test template_match when both observed and template spectra have different wavelength axis using resampling
     """
+    np.random.seed(42)
+
     # Create test spectra
     spec_axis1 = np.linspace(0, 50, 50) * u.AA
     spec_axis2 = np.linspace(0, 50, 50) * u.AA
@@ -89,14 +48,14 @@ def test_template_match_with_resample():
     # Get result from template_match
     tm_result = template_match.template_match(spec, spec1)
 
-    result = get_chi_square(spec, spec1)
-
-    assert result[1] == tm_result[1]
+    assert tm_result[1] == 40093.28353756253
 
 def test_template_match_list():
     """
     Test template_match when template spectra are in a list
     """
+    np.random.seed(42)
+
     # Create test spectra
     spec_axis1 = np.linspace(0, 50, 50) * u.AA
     spec_axis2 = np.linspace(0, 50, 50) * u.AA
@@ -117,21 +76,14 @@ def test_template_match_list():
     # Get result from template_match
     tm_result = template_match.template_match(spec, template_list)
 
-    # Loop through list in order to find out which template spectra has the smallest chi square
-    chi2_min = None
-
-    for spectrum in template_list:
-        normalized_spectral_template, chi2 = get_chi_square(spec, spectrum)
-        if chi2_min is None or chi2 < chi2_min:
-            chi2_min = chi2
-
-
-    assert chi2_min == tm_result[1]
+    assert tm_result[1] == 40093.28353756253
 
 def test_template_match_spectrum_collection():
     """
     Test template_match when template spectra are in a SpectrumCollection object
     """
+    np.random.seed(42)
+
     # Create test spectra
     spec_axis1 = np.linspace(0, 50, 50) * u.AA
     spec_axis2 = np.linspace(0, 50, 50) * u.AA
@@ -152,12 +104,4 @@ def test_template_match_spectrum_collection():
     # Get result from template_match
     tm_result = template_match.template_match(spec, spec_coll)
 
-    # Loop through SpectrumCollection in order to find out which template spectra has the smallest chi square
-    chi2_min = None
-
-    for spectrum in spec_coll:
-        normalized_template_spectrum, chi2 = get_chi_square(spec, spectrum)
-        if chi2_min is None or chi2 < chi2_min:
-            chi2_min = chi2
-
-    assert chi2_min == tm_result[1]
+    assert tm_result[1] == 40093.28353756253

From be8e3899701af45eab9c6202b944a583d70e18c9 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Tue, 27 Aug 2019 22:42:14 -0400
Subject: [PATCH 15/31] trying to pass tests

---
 specutils/tests/test_template_match.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index 272a16fae..5cc7034d0 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -11,6 +11,7 @@ def test_template_match_spectrum():
     """
     Test template_match when both observed and template spectra have the same wavelength axis
     """
+    # Seed np.random so that results are consistent
     np.random.seed(42)
 
     # Create test spectra

From d01890c7bfc26c14fa852d89885b1c689150ee76 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 6 Sep 2019 09:40:25 -0400
Subject: [PATCH 16/31] improved tests and comments

---
 specutils/analysis/template_match.py   | 16 ++++++++--------
 specutils/tests/test_template_match.py | 14 ++++++++++++++
 2 files changed, 22 insertions(+), 8 deletions(-)

diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_match.py
index 63dfe7176..66176fe61 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_match.py
@@ -9,9 +9,9 @@ def _normalize(observed_spectrum, template_spectrum):
 
     Parameters
     ----------
-    observed_spectrum : `Spectrum1D`
+    observed_spectrum : :class:`~specutils.Spectrum1D`
         the observed spectrum
-    template_spectrum : `Spectrum1D`
+    template_spectrum : :class:`~specutils.Spectrum1D`
         the template spectrum, which needs to be normalized in order to be compared with the observed_spectrum
 
     Returns
@@ -30,14 +30,14 @@ def _template_match(observed_spectrum, template_spectrum):
 
     Parameters
     ----------
-    observed_spectrum : `Spectrum1D`
+    observed_spectrum : :class:`~specutils.Spectrum1D`
         the observed spectrum
-    template_spectrum : `Spectrum1D`
+    template_spectrum : :class:`~specutils.Spectrum1D`
         the template spectrum, which will be resampled to match the wavelength of the observed_spectrum
 
     Returns
     -------
-    normalized_template_spectrum : `Spectrum1D`
+    normalized_template_spectrum : :class:`~specutils.Spectrum1D`
         the template_spectrum that has been normalized
     chi2 : `float`
         the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
@@ -73,15 +73,15 @@ def template_match(observed_spectrum, spectral_templates):
 
     Parameters
     ----------
-    observed_spectrum : `Spectrum1D`
+    observed_spectrum : :class:`~specutils.Spectrum1D`
         the observed spectrum
-    spectral_templates : `Spectrum1D` or `SpectrumCollection` or `list`
+    spectral_templates : :class:`~specutils.Spectrum1D` or :class:`~specutils.SpectrumCollection` or `list`
         the template spectra, which will be resampled and normalized and compared to the observed_spectrum, where the
         smallest chi2 and normalized_template_spectrum will be returned
 
     Returns
     -------
-    normalized_template_spectrum : `Spectrum1D`
+    normalized_template_spectrum : :class:`~specutils.Spectrum1D`
         the template_spectrum that has been normalized
     chi2 : `float`
         the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
index 5cc7034d0..682f6da81 100644
--- a/specutils/tests/test_template_match.py
+++ b/specutils/tests/test_template_match.py
@@ -5,8 +5,12 @@
 from ..spectra.spectrum1d import Spectrum1D
 from ..spectra.spectrum_collection import SpectrumCollection
 from ..analysis import template_match
+from astropy.tests.helper import quantity_allclose
 
 
+# TODO: Add some tests that are outliers: where the observed and template do not overlap (what happens?),
+# TODO: where there is minimal overlap (1 point)?
+
 def test_template_match_spectrum():
     """
     Test template_match when both observed and template spectra have the same wavelength axis
@@ -27,6 +31,11 @@ def test_template_match_spectrum():
     # Get result from template_match
     tm_result = template_match.template_match(spec, spec1)
 
+    # Create new spectrum for comparison
+    spec_result = Spectrum1D(spectral_axis=spec_axis,
+                       flux=spec1.flux * template_match._normalize(spec, spec1))
+
+    assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
     assert tm_result[1] == 40093.28353756253
 
 def test_template_match_with_resample():
@@ -49,6 +58,11 @@ def test_template_match_with_resample():
     # Get result from template_match
     tm_result = template_match.template_match(spec, spec1)
 
+    # Create new spectrum for comparison
+    spec_result = Spectrum1D(spectral_axis=spec_axis1,
+                       flux=spec1.flux * template_match._normalize(spec, spec1))
+
+    assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
     assert tm_result[1] == 40093.28353756253
 
 def test_template_match_list():

From fc198ca7213c9f474e5b8d12c462c74b982f9f70 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Thu, 12 Sep 2019 09:00:47 -0400
Subject: [PATCH 17/31] implemented Erik's suggestions and docs

---
 docs/analysis.rst                             |  23 ++++
 specutils/analysis/__init__.py                |   1 +
 ...mplate_match.py => template_comparison.py} |  73 ++++++++---
 specutils/tests/test_template_comparison.py   | 122 ++++++++++++++++++
 4 files changed, 199 insertions(+), 20 deletions(-)
 rename specutils/analysis/{template_match.py => template_comparison.py} (55%)
 create mode 100644 specutils/tests/test_template_comparison.py

diff --git a/docs/analysis.rst b/docs/analysis.rst
index c7894db08..2beb15206 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -164,6 +164,29 @@ Each of the width analysis functions are applied to this spectrum below:
    <Quantity 89.99999455 GHz>
 
 
+Template Comparison
+-------------------
+
+With one observed spectrum and n template spectra, the process to do template matching is:
+
+    1. Move the templates as close as possible to the observed spectrum by doing redshifting, matching the resolution,
+       and matching the wavelength spacing
+    2. Then loop through all of the n template spectra and run chi square on them using the observed spectrum as the
+       expected frequency and the template spectrum as the observed
+    3. Once you have a corresponding chi square for each template spectrum, return the lowest chi square and its
+       corresponding template spectrum (normalized) and the index of the template spectrum if the original template
+       parameter is iterable
+An example of how to use template matching in `specutils.analysis.template_comparison` is:
+
+.. code-block:: python
+
+   >>> from specutils.analysis import template_comparison
+   >>> spec_axis = np.linspace(0, 50, 50) * u.AA
+   >>> spec = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+   >>> spec1 = 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(spec, spec1) # doctest:+FLOAT_CMP
+
+
 Reference/API
 -------------
 .. automodapi:: specutils.analysis
diff --git a/specutils/analysis/__init__.py b/specutils/analysis/__init__.py
index 84dc60c7e..b864d991a 100644
--- a/specutils/analysis/__init__.py
+++ b/specutils/analysis/__init__.py
@@ -2,3 +2,4 @@
 from .uncertainty import *  # noqa
 from .location import *  # noqa
 from .width import *  # noqa
+from .template_comparison import *
diff --git a/specutils/analysis/template_match.py b/specutils/analysis/template_comparison.py
similarity index 55%
rename from specutils/analysis/template_match.py
rename to specutils/analysis/template_comparison.py
index 66176fe61..7a9d61d11 100644
--- a/specutils/analysis/template_match.py
+++ b/specutils/analysis/template_comparison.py
@@ -1,9 +1,11 @@
 from ..spectra.spectrum1d import Spectrum1D
 from ..spectra.spectrum_collection import SpectrumCollection
 from ..manipulation import FluxConservingResampler
+from ..manipulation import LinearInterpolatedResampler
+from ..manipulation import SplineInterpolatedResampler
 import numpy as np
 
-def _normalize(observed_spectrum, template_spectrum):
+def _normalize_for_template_matching(observed_spectrum, template_spectrum):
     """
     Calculate a scale factor to be applied to the template spectrum so the total flux in both spectra will be the same.
 
@@ -23,7 +25,30 @@ def _normalize(observed_spectrum, template_spectrum):
     denom = np.sum((template_spectrum.flux/observed_spectrum.uncertainty.array)**2)
     return num/denom
 
-def _template_match(observed_spectrum, template_spectrum):
+def _resample(resample_method):
+    """
+    Find the user preferred method of resampling the template spectrum to fit the observed spectrum.
+
+    Parameters
+    ----------
+    resample_method: `string`
+        The type of resampling to be done on the template spectrum
+
+    Returns
+    -------
+    :class:`~specutils.ResamplerBase`
+        This is the actual class that will handle the resampling
+    """
+    if resample_method == "flux_conserving":
+        return FluxConservingResampler()
+    elif resample_method == "linear_interpolated":
+        return LinearInterpolatedResampler()
+    elif resample_method == "spline_interpolated":
+        return SplineInterpolatedResampler()
+    else:
+        return None
+
+def _template_match(observed_spectrum, template_spectrum, resample_method):
     """
     Resample the template spectrum to match the wavelength of the observed spectrum.
     Then, calculate chi2 on the flux of the two spectra.
@@ -43,11 +68,12 @@ def _template_match(observed_spectrum, template_spectrum):
         the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
     """
     # Resample template
-    fluxc_resample = FluxConservingResampler()
-    template_obswavelength = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
+    if _resample(resample_method) != 0:
+        fluxc_resample = _resample(resample_method)
+        template_obswavelength = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
 
     # Normalize spectra
-    normalization = _normalize(observed_spectrum, template_obswavelength)
+    normalization = _normalize_for_template_matching(observed_spectrum, template_obswavelength)
 
     # Numerator
     num_right = normalization*template_obswavelength.flux
@@ -66,17 +92,18 @@ def _template_match(observed_spectrum, template_spectrum):
 
     return normalized_template_spectrum, chi2
 
-def template_match(observed_spectrum, spectral_templates):
+def template_match(observed_spectrum, spectral_templates, resample_method="flux_conserving"):
     """
-    Find what instance spectral_templates is and run _template_match accordingly. If two template_spectra have the same
-    chi2, the first template is returned
+    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 first template is returned
 
     Parameters
     ----------
     observed_spectrum : :class:`~specutils.Spectrum1D`
         the observed spectrum
-    spectral_templates : :class:`~specutils.Spectrum1D` or :class:`~specutils.SpectrumCollection` or `list`
-        the template spectra, which will be resampled and normalized and compared to the observed_spectrum, where the
+    spectral_templates : :class:`~specutils.Spectrum1D` or :class:`~specutils.SpectrumCollection` or `list` or anything
+        that will give a single :class:`~specutils.Spectrum1D` when iterated over.
+        The template spectra, which will be resampled and normalized and compared to the observed_spectrum, where the
         smallest chi2 and normalized_template_spectrum will be returned
 
     Returns
@@ -86,20 +113,26 @@ def template_match(observed_spectrum, spectral_templates):
     chi2 : `float`
         the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
     """
-    if isinstance(spectral_templates, Spectrum1D):
-        normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectral_templates)
+    if hasattr(spectral_templates, 'flux') and len(spectral_templates.flux.shape)==1:
+        normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectral_templates, resample_method)
         return normalized_spectral_template, chi2
 
     # Loop through spectra in list and return spectrum with lowest chi square
     # and its corresponding chi square
-    elif isinstance(spectral_templates, list) or isinstance(spectral_templates, SpectrumCollection):
+    else:
         chi2_min = None
         smallest_chi_spec = None
 
-        for spectrum in spectral_templates:
-            normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectrum)
-            if chi2_min is None or chi2 < chi2_min:
-                chi2_min = chi2
-                smallest_chi_spec = normalized_spectral_template
-
-        return smallest_chi_spec, chi2_min
+        index = 0
+        try:
+            for spectrum in spectral_templates:
+                normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectrum, resample_method)
+                if chi2_min is None or chi2 < chi2_min:
+                    chi2_min = chi2
+                    smallest_chi_spec = normalized_spectral_template
+                    smallest_chi_index = index
+                index+=1
+        except Exception as e:
+            print("Parameter spectral_templates is not iterable. The following error was fired: {}".format(e))
+
+        return smallest_chi_spec, chi2_min, smallest_chi_index
diff --git a/specutils/tests/test_template_comparison.py b/specutils/tests/test_template_comparison.py
new file mode 100644
index 000000000..c5cfc70a6
--- /dev/null
+++ b/specutils/tests/test_template_comparison.py
@@ -0,0 +1,122 @@
+import astropy.units as u
+import numpy as np
+from astropy.nddata import StdDevUncertainty
+
+from ..spectra.spectrum1d import Spectrum1D
+from ..spectra.spectrum_collection import SpectrumCollection
+from ..analysis import template_comparison
+from astropy.tests.helper import quantity_allclose
+
+
+# TODO: Add some tests that are outliers: where the observed and template do not overlap (what happens?),
+# TODO: where there is minimal overlap (1 point)?
+
+def test_template_match_spectrum():
+    """
+    Test template_match when both observed and template spectra have the same wavelength axis
+    """
+    # Seed np.random so that results are consistent
+    np.random.seed(42)
+
+    # Create test spectra
+    spec_axis = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis,
+                       flux=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)
+
+    # Create new spectrum for comparison
+    spec_result = Spectrum1D(spectral_axis=spec_axis,
+                             flux=spec1.flux * template_comparison._normalize_for_template_matching(spec, spec1))
+
+    assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
+    assert tm_result[1] == 40093.28353756253
+
+def test_template_match_with_resample():
+    """
+    Test template_match when both observed and template spectra have different wavelength axis using resampling
+    """
+    np.random.seed(42)
+
+    # Create test spectra
+    spec_axis1 = np.linspace(0, 50, 50) * u.AA
+    spec_axis2 = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis1,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=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)
+
+    # Create new spectrum for comparison
+    spec_result = Spectrum1D(spectral_axis=spec_axis1,
+                             flux=spec1.flux * template_comparison._normalize_for_template_matching(spec, spec1))
+
+    assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
+    assert tm_result[1] == 40093.28353756253
+
+def test_template_match_list():
+    """
+    Test template_match when template spectra are in a list
+    """
+    np.random.seed(42)
+
+    # Create test spectra
+    spec_axis1 = np.linspace(0, 50, 50) * u.AA
+    spec_axis2 = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis1,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+    spec2 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    # Combine spectra into list
+    template_list = [spec1, spec2]
+
+    # Get result from template_match
+    tm_result = template_comparison.template_match(spec, template_list)
+
+    assert tm_result[1] == 40093.28353756253
+
+def test_template_match_spectrum_collection():
+    """
+    Test template_match when template spectra are in a SpectrumCollection object
+    """
+    np.random.seed(42)
+
+    # Create test spectra
+    spec_axis1 = np.linspace(0, 50, 50) * u.AA
+    spec_axis2 = np.linspace(0, 50, 50) * u.AA
+    spec = Spectrum1D(spectral_axis=spec_axis1,
+                      flux=np.random.randn(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    spec1 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+    spec2 = Spectrum1D(spectral_axis=spec_axis2,
+                       flux=np.random.randn(50) * u.Jy,
+                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+
+    # Combine spectra into SpectrumCollection object
+    spec_coll = SpectrumCollection.from_spectra([spec1, spec2])
+
+    # Get result from template_match
+    tm_result = template_comparison.template_match(spec, spec_coll)
+
+    assert tm_result[1] == 40093.28353756253

From 28515d8fa84e1a6d8fe6321d5cebb6fed7123a33 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Thu, 12 Sep 2019 09:16:38 -0400
Subject: [PATCH 18/31] removed old tests file and updated docs

---
 docs/analysis.rst                      |   1 +
 specutils/tests/test_template_match.py | 122 -------------------------
 2 files changed, 1 insertion(+), 122 deletions(-)
 delete mode 100644 specutils/tests/test_template_match.py

diff --git a/docs/analysis.rst b/docs/analysis.rst
index 2beb15206..356a736c3 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -181,6 +181,7 @@ An example of how to use template matching in `specutils.analysis.template_compa
 .. code-block:: python
 
    >>> from specutils.analysis import template_comparison
+
    >>> spec_axis = np.linspace(0, 50, 50) * u.AA
    >>> spec = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
    >>> spec1 = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
diff --git a/specutils/tests/test_template_match.py b/specutils/tests/test_template_match.py
deleted file mode 100644
index 682f6da81..000000000
--- a/specutils/tests/test_template_match.py
+++ /dev/null
@@ -1,122 +0,0 @@
-import astropy.units as u
-import numpy as np
-from astropy.nddata import StdDevUncertainty
-
-from ..spectra.spectrum1d import Spectrum1D
-from ..spectra.spectrum_collection import SpectrumCollection
-from ..analysis import template_match
-from astropy.tests.helper import quantity_allclose
-
-
-# TODO: Add some tests that are outliers: where the observed and template do not overlap (what happens?),
-# TODO: where there is minimal overlap (1 point)?
-
-def test_template_match_spectrum():
-    """
-    Test template_match when both observed and template spectra have the same wavelength axis
-    """
-    # Seed np.random so that results are consistent
-    np.random.seed(42)
-
-    # Create test spectra
-    spec_axis = np.linspace(0, 50, 50) * u.AA
-    spec = Spectrum1D(spectral_axis=spec_axis,
-                      flux=np.random.randn(50) * u.Jy,
-                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    spec1 = Spectrum1D(spectral_axis=spec_axis,
-                       flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    # Get result from template_match
-    tm_result = template_match.template_match(spec, spec1)
-
-    # Create new spectrum for comparison
-    spec_result = Spectrum1D(spectral_axis=spec_axis,
-                       flux=spec1.flux * template_match._normalize(spec, spec1))
-
-    assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
-    assert tm_result[1] == 40093.28353756253
-
-def test_template_match_with_resample():
-    """
-    Test template_match when both observed and template spectra have different wavelength axis using resampling
-    """
-    np.random.seed(42)
-
-    # Create test spectra
-    spec_axis1 = np.linspace(0, 50, 50) * u.AA
-    spec_axis2 = np.linspace(0, 50, 50) * u.AA
-    spec = Spectrum1D(spectral_axis=spec_axis1,
-                      flux=np.random.randn(50) * u.Jy,
-                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    spec1 = Spectrum1D(spectral_axis=spec_axis2,
-                       flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    # Get result from template_match
-    tm_result = template_match.template_match(spec, spec1)
-
-    # Create new spectrum for comparison
-    spec_result = Spectrum1D(spectral_axis=spec_axis1,
-                       flux=spec1.flux * template_match._normalize(spec, spec1))
-
-    assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
-    assert tm_result[1] == 40093.28353756253
-
-def test_template_match_list():
-    """
-    Test template_match when template spectra are in a list
-    """
-    np.random.seed(42)
-
-    # Create test spectra
-    spec_axis1 = np.linspace(0, 50, 50) * u.AA
-    spec_axis2 = np.linspace(0, 50, 50) * u.AA
-    spec = Spectrum1D(spectral_axis=spec_axis1,
-                      flux=np.random.randn(50) * u.Jy,
-                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    spec1 = Spectrum1D(spectral_axis=spec_axis2,
-                       flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-    spec2 = Spectrum1D(spectral_axis=spec_axis2,
-                       flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    # Combine spectra into list
-    template_list = [spec1, spec2]
-
-    # Get result from template_match
-    tm_result = template_match.template_match(spec, template_list)
-
-    assert tm_result[1] == 40093.28353756253
-
-def test_template_match_spectrum_collection():
-    """
-    Test template_match when template spectra are in a SpectrumCollection object
-    """
-    np.random.seed(42)
-
-    # Create test spectra
-    spec_axis1 = np.linspace(0, 50, 50) * u.AA
-    spec_axis2 = np.linspace(0, 50, 50) * u.AA
-    spec = Spectrum1D(spectral_axis=spec_axis1,
-                      flux=np.random.randn(50) * u.Jy,
-                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    spec1 = Spectrum1D(spectral_axis=spec_axis2,
-                       flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-    spec2 = Spectrum1D(spectral_axis=spec_axis2,
-                       flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
-    # Combine spectra into SpectrumCollection object
-    spec_coll = SpectrumCollection.from_spectra([spec1, spec2])
-
-    # Get result from template_match
-    tm_result = template_match.template_match(spec, spec_coll)
-
-    assert tm_result[1] == 40093.28353756253

From ece820dde86fd565e1c83c51add3424bdab19c77 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 13 Sep 2019 09:17:06 -0400
Subject: [PATCH 19/31] trying to fix doc test

---
 docs/analysis.rst | 1 -
 1 file changed, 1 deletion(-)

diff --git a/docs/analysis.rst b/docs/analysis.rst
index 356a736c3..2beb15206 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -181,7 +181,6 @@ An example of how to use template matching in `specutils.analysis.template_compa
 .. code-block:: python
 
    >>> from specutils.analysis import template_comparison
-
    >>> spec_axis = np.linspace(0, 50, 50) * u.AA
    >>> spec = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
    >>> spec1 = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))

From a9051746f5b17cdd2531c263f770c24232ecaad6 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 13 Sep 2019 09:36:01 -0400
Subject: [PATCH 20/31] continued doc string fixing

---
 docs/analysis.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/analysis.rst b/docs/analysis.rst
index 2beb15206..12a321077 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -176,6 +176,7 @@ With one observed spectrum and n template spectra, the process to do template ma
     3. Once you have a corresponding chi square for each template spectrum, return the lowest chi square and its
        corresponding template spectrum (normalized) and the index of the template spectrum if the original template
        parameter is iterable
+
 An example of how to use template matching in `specutils.analysis.template_comparison` is:
 
 .. code-block:: python

From 6e19db23132265adda49275a38854bc8a3f29a83 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 13 Sep 2019 09:46:29 -0400
Subject: [PATCH 21/31] more doc string changes

---
 docs/analysis.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/analysis.rst b/docs/analysis.rst
index 12a321077..804c991f8 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -177,7 +177,7 @@ With one observed spectrum and n template spectra, the process to do template ma
        corresponding template spectrum (normalized) and the index of the template spectrum if the original template
        parameter is iterable
 
-An example of how to use template matching in `specutils.analysis.template_comparison` is:
+An example of how to do template matching in is:
 
 .. code-block:: python
 

From 06337eb9c87d67f4c506dd0b530bb14f63603ad9 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 13 Sep 2019 12:52:20 -0400
Subject: [PATCH 22/31] minor fixes

---
 docs/analysis.rst                         | 6 +++---
 specutils/analysis/template_comparison.py | 4 +++-
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/docs/analysis.rst b/docs/analysis.rst
index 804c991f8..afe419f9f 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -183,9 +183,9 @@ An example of how to do template matching in is:
 
    >>> from specutils.analysis import template_comparison
    >>> spec_axis = np.linspace(0, 50, 50) * u.AA
-   >>> spec = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-   >>> spec1 = 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(spec, spec1) # doctest:+FLOAT_CMP
+   >>> observed_spectrum = Spectrum1D(spectral_axis=spec_axis, 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(2, 50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(2, 50), unit='Jy'))
+   >>> tm_result = template_comparison.template_match(observed_spectrum, spectral_template) # doctest:+FLOAT_CMP
 
 
 Reference/API
diff --git a/specutils/analysis/template_comparison.py b/specutils/analysis/template_comparison.py
index 7a9d61d11..f2722a144 100644
--- a/specutils/analysis/template_comparison.py
+++ b/specutils/analysis/template_comparison.py
@@ -112,6 +112,8 @@ def template_match(observed_spectrum, spectral_templates, resample_method="flux_
         the template_spectrum that has been normalized
     chi2 : `float`
         the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
+    smallest_chi_index : `int`
+        the index of the spectrum with the smallest chi2 in spectral_templates
     """
     if hasattr(spectral_templates, 'flux') and len(spectral_templates.flux.shape)==1:
         normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectral_templates, resample_method)
@@ -131,7 +133,7 @@ def template_match(observed_spectrum, spectral_templates, resample_method="flux_
                     chi2_min = chi2
                     smallest_chi_spec = normalized_spectral_template
                     smallest_chi_index = index
-                index+=1
+                index += 1
         except Exception as e:
             print("Parameter spectral_templates is not iterable. The following error was fired: {}".format(e))
 

From 99e4843e9ae6935409e5f2968f728b4b1db9a6f6 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 13 Sep 2019 13:00:56 -0400
Subject: [PATCH 23/31] removed vector spectrum1d in doc

---
 docs/analysis.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/analysis.rst b/docs/analysis.rst
index afe419f9f..b69d9b365 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -184,7 +184,7 @@ An example of how to do template matching in is:
    >>> from specutils.analysis import template_comparison
    >>> spec_axis = np.linspace(0, 50, 50) * u.AA
    >>> observed_spectrum = Spectrum1D(spectral_axis=spec_axis, 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(2, 50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(2, 50), unit='Jy'))
+   >>> spectral_template = Spectrum1D(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(2, 50), unit='Jy'))
    >>> tm_result = template_comparison.template_match(observed_spectrum, spectral_template) # doctest:+FLOAT_CMP
 
 

From b70eb33f33cf2b6fe759d28706797392943a9aea Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Fri, 13 Sep 2019 13:19:29 -0400
Subject: [PATCH 24/31] removed vector spectrum1d totally in doc

---
 docs/analysis.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/analysis.rst b/docs/analysis.rst
index b69d9b365..179136692 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -184,7 +184,7 @@ An example of how to do template matching in is:
    >>> from specutils.analysis import template_comparison
    >>> spec_axis = np.linspace(0, 50, 50) * u.AA
    >>> observed_spectrum = Spectrum1D(spectral_axis=spec_axis, 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(2, 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, spectral_template) # doctest:+FLOAT_CMP
 
 

From d2668de91d98616edc5c77e654fef8c8fa0e1c32 Mon Sep 17 00:00:00 2001
From: Jesse Averbukh <javerbukh@stsci.edu>
Date: Sat, 14 Sep 2019 09:29:28 -0400
Subject: [PATCH 25/31] replace print with logging

---
 specutils/analysis/template_comparison.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/specutils/analysis/template_comparison.py b/specutils/analysis/template_comparison.py
index f2722a144..f3a925249 100644
--- a/specutils/analysis/template_comparison.py
+++ b/specutils/analysis/template_comparison.py
@@ -4,6 +4,7 @@
 from ..manipulation import LinearInterpolatedResampler
 from ..manipulation import SplineInterpolatedResampler
 import numpy as np
+import logging
 
 def _normalize_for_template_matching(observed_spectrum, template_spectrum):
     """
@@ -135,6 +136,6 @@ def template_match(observed_spectrum, spectral_templates, resample_method="flux_
                     smallest_chi_index = index
                 index += 1
         except Exception as e:
-            print("Parameter spectral_templates is not iterable. The following error was fired: {}".format(e))
+            logging.warning("Parameter spectral_templates is not iterable. The following error was fired: {}".format(e))
 
         return smallest_chi_spec, chi2_min, smallest_chi_index

From 9f08baf2b74630a19566da56a9f82208646472ad Mon Sep 17 00:00:00 2001
From: Nicholas Earl <contact@nicholasearl.me>
Date: Sat, 14 Sep 2019 10:39:40 -0400
Subject: [PATCH 26/31] Return chi2 as value; minor code improvements

---
 specutils/analysis/template_comparison.py | 103 ++++++++++++----------
 1 file changed, 58 insertions(+), 45 deletions(-)

diff --git a/specutils/analysis/template_comparison.py b/specutils/analysis/template_comparison.py
index f3a925249..863015587 100644
--- a/specutils/analysis/template_comparison.py
+++ b/specutils/analysis/template_comparison.py
@@ -1,10 +1,12 @@
+import logging
+
+import numpy as np
+
+from ..manipulation import (FluxConservingResampler,
+                            LinearInterpolatedResampler,
+                            SplineInterpolatedResampler)
 from ..spectra.spectrum1d import Spectrum1D
 from ..spectra.spectrum_collection import SpectrumCollection
-from ..manipulation import FluxConservingResampler
-from ..manipulation import LinearInterpolatedResampler
-from ..manipulation import SplineInterpolatedResampler
-import numpy as np
-import logging
 
 def _normalize_for_template_matching(observed_spectrum, template_spectrum):
     """
@@ -26,6 +28,7 @@ def _normalize_for_template_matching(observed_spectrum, template_spectrum):
     denom = np.sum((template_spectrum.flux/observed_spectrum.uncertainty.array)**2)
     return num/denom
 
+
 def _resample(resample_method):
     """
     Find the user preferred method of resampling the template spectrum to fit the observed spectrum.
@@ -49,6 +52,7 @@ def _resample(resample_method):
     else:
         return None
 
+
 def _template_match(observed_spectrum, template_spectrum, resample_method):
     """
     Resample the template spectrum to match the wavelength of the observed spectrum.
@@ -64,20 +68,23 @@ def _template_match(observed_spectrum, template_spectrum, resample_method):
     Returns
     -------
     normalized_template_spectrum : :class:`~specutils.Spectrum1D`
-        the template_spectrum that has been normalized
+        The normalized spectrum template.
     chi2 : `float`
-        the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
+        The chi2 of the flux of the observed spectrum and the flux of the 
+        normalized template spectrum.
     """
     # Resample template
     if _resample(resample_method) != 0:
         fluxc_resample = _resample(resample_method)
-        template_obswavelength = fluxc_resample(template_spectrum, observed_spectrum.wavelength)
+        template_obswavelength = fluxc_resample(template_spectrum, 
+                                                observed_spectrum.wavelength)
 
     # Normalize spectra
-    normalization = _normalize_for_template_matching(observed_spectrum, template_obswavelength)
+    normalization = _normalize_for_template_matching(observed_spectrum,
+                                                     template_obswavelength)
 
     # Numerator
-    num_right = normalization*template_obswavelength.flux
+    num_right = normalization * template_obswavelength.flux
     num = observed_spectrum.flux - num_right
 
     # Denominator
@@ -85,57 +92,63 @@ def _template_match(observed_spectrum, template_spectrum, resample_method):
 
     # Get chi square
     result = (num/denom)**2
-    chi2 = np.sum(result)
+    chi2 = np.sum(result.value)
 
-    # Create normalized template spectrum, which will be returned with corresponding chi2
-    normalized_template_spectrum = Spectrum1D(spectral_axis=template_spectrum.spectral_axis,
-                                              flux=template_spectrum.flux*normalization)
+    # Create normalized template spectrum, which will be returned with 
+    # corresponding chi2
+    normalized_template_spectrum = Spectrum1D(
+        spectral_axis=template_spectrum.spectral_axis,
+        flux=template_spectrum.flux*normalization)
 
     return normalized_template_spectrum, chi2
 
+
 def template_match(observed_spectrum, spectral_templates, resample_method="flux_conserving"):
     """
-    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 first template is returned
+    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 
+    first template is returned.
 
     Parameters
     ----------
     observed_spectrum : :class:`~specutils.Spectrum1D`
-        the observed spectrum
-    spectral_templates : :class:`~specutils.Spectrum1D` or :class:`~specutils.SpectrumCollection` or `list` or anything
-        that will give a single :class:`~specutils.Spectrum1D` when iterated over.
-        The template spectra, which will be resampled and normalized and compared to the observed_spectrum, where the
-        smallest chi2 and normalized_template_spectrum will be returned
+        The observed spectrum.
+    spectral_templates : :class:`~specutils.Spectrum1D` or :class:`~specutils.SpectrumCollection` or `list`
+        That will give a single :class:`~specutils.Spectrum1D` when iterated 
+        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.
 
     Returns
     -------
     normalized_template_spectrum : :class:`~specutils.Spectrum1D`
-        the template_spectrum that has been normalized
+        The template spectrum that has been normalized.
     chi2 : `float`
-        the chi2 of the flux of the observed_spectrum and the flux of the normalized_template_spectrum
+        The chi2 of the flux of the observed_spectrum and the flux of the 
+        normalized template spectrum.
     smallest_chi_index : `int`
-        the index of the spectrum with the smallest chi2 in spectral_templates
+        The index of the spectrum with the smallest chi2 in spectral templates.
     """
-    if hasattr(spectral_templates, 'flux') and len(spectral_templates.flux.shape)==1:
-        normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectral_templates, resample_method)
+    if hasattr(spectral_templates, 'flux') and len(spectral_templates.flux.shape) == 1:
+        normalized_spectral_template, chi2 = _template_match(
+            observed_spectrum, spectral_templates, resample_method)
+
         return normalized_spectral_template, chi2
 
-    # Loop through spectra in list and return spectrum with lowest chi square
-    # and its corresponding chi square
-    else:
-        chi2_min = None
-        smallest_chi_spec = None
-
-        index = 0
-        try:
-            for spectrum in spectral_templates:
-                normalized_spectral_template, chi2 = _template_match(observed_spectrum, spectrum, resample_method)
-                if chi2_min is None or chi2 < chi2_min:
-                    chi2_min = chi2
-                    smallest_chi_spec = normalized_spectral_template
-                    smallest_chi_index = index
-                index += 1
-        except Exception as e:
-            logging.warning("Parameter spectral_templates is not iterable. The following error was fired: {}".format(e))
-
-        return smallest_chi_spec, chi2_min, smallest_chi_index
+    # At this point, the template spectrum is either a ``SpectrumCollection``
+    # or a multi-dimensional``Spectrum1D``. Loop through the object and return 
+    # the template spectrum with the lowest chi square and its corresponding 
+    # chi square.
+    chi2_min = None
+    smallest_chi_spec = None
+
+    for index, spectrum in enumerate(spectral_templates):
+        normalized_spectral_template, chi2 = _template_match(
+            observed_spectrum, spectrum, resample_method)
+
+        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

From bd7e171591c0dd4e15e5d1d059bf95f807c554b9 Mon Sep 17 00:00:00 2001
From: Nicholas Earl <contact@nicholasearl.me>
Date: Sat, 14 Sep 2019 10:54:12 -0400
Subject: [PATCH 27/31] Improve slicing behavior on multidim `Spectrum1D`

---
 specutils/spectra/spectrum1d.py | 48 ++++++++++++++++++++++++++++++---
 specutils/tests/test_slicing.py | 17 ++++++++++++
 2 files changed, 61 insertions(+), 4 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index b45280168..bf7d406bb 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -1,11 +1,12 @@
 import logging
+from copy import deepcopy
 
 import numpy as np
 from astropy import units as u
-from astropy.nddata import NDDataRef
+from astropy.nddata import NDDataRef, NDUncertainty
 from astropy.utils.decorators import lazyproperty
-from astropy.nddata import NDUncertainty
-from ..wcs import WCSWrapper, WCSAdapter
+
+from ..wcs import WCSAdapter, WCSWrapper
 from .spectrum_mixin import OneDSpectrumMixin
 
 __all__ = ['Spectrum1D']
@@ -123,6 +124,46 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
                 raise ValueError('Flux axis ({}) and uncertainty ({}) shapes must be the same.'.format(
                     flux.shape, self.uncertainty.array.shape))
 
+    def __getitem__(self, item):
+        """
+        Override the class indexer. We do this here because there are two cases
+        for slicing on a ``Spectrum1D``:
+
+            1.) When the flux is one dimensional, indexing represents a single
+                flux value at a particular spectral axis bin, and returns a new
+                ``Spectrum1D`` where all attributes are sliced.
+            2.) When flux is multi-dimensional (i.e. several fluxes over the
+                same spectral axis), indexing returns a new ``Spectrum1D`` with
+                the sliced flux range and everything and a deep copy of all 
+                other attributes.
+
+        The second case is handled by the inerhited class, whle the former is
+        handled here.
+        """
+        if len(self.flux.shape) > 1:
+            return self._copy(flux=self.flux[item])
+
+        return super().__getitem__(item)
+
+    def _copy(self, **kwargs):
+        """
+        Peform deepcopy operations on each attribute of the ``Spectrum1D``
+        object.
+        """
+        alt_kwargs = dict(
+            flux=deepcopy(self.flux),
+            spectral_axis=deepcopy(self.spectral_axis),
+            uncertainty=deepcopy(self.uncertainty),
+            wcs=deepcopy(self.wcs),
+            mask=deepcopy(self.mask),
+            meta=deepcopy(self.meta),
+            unit=deepcopy(self.unit),
+            velocity_convention=deepcopy(self.velocity_convention),
+            rest_value=deepcopy(self.rest_value))
+
+        alt_kwargs.update(kwargs)
+
+        return self.__class__(**alt_kwargs)
 
     @property
     def frequency(self):
@@ -284,7 +325,6 @@ def __repr__(self):
 
         return result
 
-
     def spectral_resolution(self, true_dispersion, delta_dispersion, axis=-1):
         """Evaluate the probability distribution of the spectral resolution.
 
diff --git a/specutils/tests/test_slicing.py b/specutils/tests/test_slicing.py
index 03c60e03d..88a7c9597 100644
--- a/specutils/tests/test_slicing.py
+++ b/specutils/tests/test_slicing.py
@@ -1,5 +1,6 @@
 import astropy.units as u
 import astropy.wcs as fitswcs
+from astropy.tests.helper import quantity_allclose
 import numpy as np
 from numpy.testing import assert_allclose
 
@@ -72,6 +73,7 @@ def test_slicing():
     assert sub_spec.frequency.unit == u.GHz
     assert np.allclose(sub_spec.frequency.value, np.array([74948.1145, 59958.4916, 49965.40966667, 42827.494]))
 
+
 def test_slicing_with_fits():
     my_wcs = fitswcs.WCS(header={'CDELT1': 1, 'CRVAL1': 6562.8, 'CUNIT1': 'Angstrom',
                                  'CTYPE1': 'WAVE', 'RESTFRQ': 1400000000, 'CRPIX1': 25})
@@ -82,3 +84,18 @@ def test_slicing_with_fits():
     assert isinstance(spec_slice, Spectrum1D)
     assert spec_slice.flux.size == 4
     assert np.allclose(spec_slice.wcs.pixel_to_world([6, 7, 8, 9]).value, spec.wcs.pixel_to_world([6, 7, 8, 9]).value)
+
+
+def test_slicing_multidim():
+    spec = Spectrum1D(spectral_axis=np.arange(10) * u.AA,
+                      flux=np.random.sample((5, 10)) * u.Jy)
+
+    spec1 = spec[0]
+    spec2 = spec[1:3]
+
+    assert spec1.flux[0] == spec.flux[0][0]
+    assert quantity_allclose(spec1.spectral_axis, spec.spectral_axis)
+    assert spec.flux.shape[1:] == spec1.flux.shape
+
+    assert quantity_allclose(spec2.flux, spec.flux[1:3])
+    assert quantity_allclose(spec2.spectral_axis, spec.spectral_axis)

From b237db9ebd35e346882eb613f948b1f621d2770e Mon Sep 17 00:00:00 2001
From: Nicholas Earl <contact@nicholasearl.me>
Date: Sat, 14 Sep 2019 10:54:35 -0400
Subject: [PATCH 28/31] Include test for multidim template spectrum

---
 specutils/tests/test_template_comparison.py | 33 ++++++++++++++++++---
 1 file changed, 29 insertions(+), 4 deletions(-)

diff --git a/specutils/tests/test_template_comparison.py b/specutils/tests/test_template_comparison.py
index c5cfc70a6..cbfb636e5 100644
--- a/specutils/tests/test_template_comparison.py
+++ b/specutils/tests/test_template_comparison.py
@@ -38,6 +38,7 @@ def test_template_match_spectrum():
     assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
     assert tm_result[1] == 40093.28353756253
 
+
 def test_template_match_with_resample():
     """
     Test template_match when both observed and template spectra have different wavelength axis using resampling
@@ -65,6 +66,7 @@ def test_template_match_with_resample():
     assert quantity_allclose(tm_result[0].flux, spec_result.flux, atol=0.01*u.Jy)
     assert tm_result[1] == 40093.28353756253
 
+
 def test_template_match_list():
     """
     Test template_match when template spectra are in a list
@@ -93,6 +95,7 @@ def test_template_match_list():
 
     assert tm_result[1] == 40093.28353756253
 
+
 def test_template_match_spectrum_collection():
     """
     Test template_match when template spectra are in a SpectrumCollection object
@@ -104,14 +107,13 @@ def test_template_match_spectrum_collection():
     spec_axis2 = np.linspace(0, 50, 50) * u.AA
     spec = Spectrum1D(spectral_axis=spec_axis1,
                       flux=np.random.randn(50) * u.Jy,
-                      uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
-
+                      uncertainty=StdDevUncertainty(np.random.sample(50)))
     spec1 = Spectrum1D(spectral_axis=spec_axis2,
                        flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+                       uncertainty=StdDevUncertainty(np.random.sample(50)))
     spec2 = Spectrum1D(spectral_axis=spec_axis2,
                        flux=np.random.randn(50) * u.Jy,
-                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
+                       uncertainty=StdDevUncertainty(np.random.sample(50)))
 
     # Combine spectra into SpectrumCollection object
     spec_coll = SpectrumCollection.from_spectra([spec1, spec2])
@@ -120,3 +122,26 @@ def test_template_match_spectrum_collection():
     tm_result = template_comparison.template_match(spec, spec_coll)
 
     assert tm_result[1] == 40093.28353756253
+
+
+def test_template_match_multidim_spectrum():
+    """
+    Test template matching with a multi-dimensional Spectrum1D object.
+    """
+    np.random.seed(42)
+
+    # Create test spectra
+    spec_axis1 = np.linspace(0, 50, 50) * u.AA
+    spec_axis2 = np.linspace(0, 50, 50) * u.AA
+
+    spec = Spectrum1D(spectral_axis=spec_axis1,
+                      flux=np.random.sample(50) * u.Jy,
+                      uncertainty=StdDevUncertainty(np.random.sample(50)))
+    multidim_spec = Spectrum1D(spectral_axis=spec_axis2,
+                               flux=np.random.sample((2, 50)) * u.Jy,
+                               uncertainty=StdDevUncertainty(np.random.sample(50)))
+
+    # Get result from template_match
+    tm_result = template_comparison.template_match(spec, multidim_spec)
+
+    assert tm_result[1] == 250.26870401777543

From ea186c1697ab555612d23e5d4789d7549520e5a1 Mon Sep 17 00:00:00 2001
From: Nicholas Earl <contact@nicholasearl.me>
Date: Sat, 14 Sep 2019 11:22:22 -0400
Subject: [PATCH 29/31] Fix docstrings

---
 specutils/spectra/spectrum1d.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index bf7d406bb..cfde4e2a9 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -134,10 +134,9 @@ def __getitem__(self, item):
                 ``Spectrum1D`` where all attributes are sliced.
             2.) When flux is multi-dimensional (i.e. several fluxes over the
                 same spectral axis), indexing returns a new ``Spectrum1D`` with
-                the sliced flux range and everything and a deep copy of all 
-                other attributes.
+                the sliced flux range and a deep copy of all other attributes.
 
-        The second case is handled by the inerhited class, whle the former is
+        The first case is handled by the parent class, while the second is
         handled here.
         """
         if len(self.flux.shape) > 1:
@@ -147,7 +146,7 @@ def __getitem__(self, item):
 
     def _copy(self, **kwargs):
         """
-        Peform deepcopy operations on each attribute of the ``Spectrum1D``
+        Peform deep copy operations on each attribute of the ``Spectrum1D``
         object.
         """
         alt_kwargs = dict(

From 958969567c7566c4f21098ea383b540053a07f21 Mon Sep 17 00:00:00 2001
From: Nicholas Earl <contact@nicholasearl.me>
Date: Mon, 16 Sep 2019 17:40:39 -0400
Subject: [PATCH 30/31] Handle uncertainty shapes correctly

---
 specutils/spectra/spectrum1d.py             | 4 +++-
 specutils/tests/test_template_comparison.py | 2 +-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index cfde4e2a9..c12009de1 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -140,7 +140,9 @@ def __getitem__(self, item):
         handled here.
         """
         if len(self.flux.shape) > 1:
-            return self._copy(flux=self.flux[item])
+            return self._copy(
+                flux=self.flux[item], uncertainty=self.uncertainty[item] 
+                    if self.uncertainty is not None else None)
 
         return super().__getitem__(item)
 
diff --git a/specutils/tests/test_template_comparison.py b/specutils/tests/test_template_comparison.py
index cbfb636e5..6db92ce35 100644
--- a/specutils/tests/test_template_comparison.py
+++ b/specutils/tests/test_template_comparison.py
@@ -139,7 +139,7 @@ def test_template_match_multidim_spectrum():
                       uncertainty=StdDevUncertainty(np.random.sample(50)))
     multidim_spec = Spectrum1D(spectral_axis=spec_axis2,
                                flux=np.random.sample((2, 50)) * u.Jy,
-                               uncertainty=StdDevUncertainty(np.random.sample(50)))
+                               uncertainty=StdDevUncertainty(np.random.sample((2, 50))))
 
     # Get result from template_match
     tm_result = template_comparison.template_match(spec, multidim_spec)

From 7eef14f1163973b41b63cdcd2085e950a173c57a Mon Sep 17 00:00:00 2001
From: Nicholas Earl <contact@nicholasearl.me>
Date: Mon, 16 Sep 2019 19:03:45 -0400
Subject: [PATCH 31/31] Improve pep8 compliance

---
 specutils/analysis/template_comparison.py | 68 +++++++++++++----------
 1 file changed, 38 insertions(+), 30 deletions(-)

diff --git a/specutils/analysis/template_comparison.py b/specutils/analysis/template_comparison.py
index 863015587..2d69d6466 100644
--- a/specutils/analysis/template_comparison.py
+++ b/specutils/analysis/template_comparison.py
@@ -1,37 +1,41 @@
-import logging
-
 import numpy as np
 
 from ..manipulation import (FluxConservingResampler,
                             LinearInterpolatedResampler,
                             SplineInterpolatedResampler)
 from ..spectra.spectrum1d import Spectrum1D
-from ..spectra.spectrum_collection import SpectrumCollection
 
 def _normalize_for_template_matching(observed_spectrum, template_spectrum):
     """
-    Calculate a scale factor to be applied to the template spectrum so the total flux in both spectra will be the same.
+    Calculate a scale factor to be applied to the template spectrum so the
+    total flux in both spectra will be the same.
 
     Parameters
     ----------
     observed_spectrum : :class:`~specutils.Spectrum1D`
-        the observed spectrum
+        The observed spectrum.
     template_spectrum : :class:`~specutils.Spectrum1D`
-        the template spectrum, which needs to be normalized in order to be compared with the observed_spectrum
+        The template spectrum, which needs to be normalized in order to be
+        compared with the observed spectrum.
 
     Returns
     -------
     `float`
-        A float which will normalize the template spectrum's flux so that it can be compared to the observed spectrum
+        A float which will normalize the template spectrum's flux so that it
+        can be compared to the observed spectrum.
     """
-    num = np.sum((observed_spectrum.flux*template_spectrum.flux)/(observed_spectrum.uncertainty.array**2))
-    denom = np.sum((template_spectrum.flux/observed_spectrum.uncertainty.array)**2)
+    num = np.sum((observed_spectrum.flux*template_spectrum.flux)/
+                 (observed_spectrum.uncertainty.array**2))
+    denom = np.sum((template_spectrum.flux/
+                    observed_spectrum.uncertainty.array)**2)
+
     return num/denom
 
 
 def _resample(resample_method):
     """
-    Find the user preferred method of resampling the template spectrum to fit the observed spectrum.
+    Find the user preferred method of resampling the template spectrum to fit
+    the observed spectrum.
 
     Parameters
     ----------
@@ -45,38 +49,41 @@ def _resample(resample_method):
     """
     if resample_method == "flux_conserving":
         return FluxConservingResampler()
-    elif resample_method == "linear_interpolated":
+
+    if resample_method == "linear_interpolated":
         return LinearInterpolatedResampler()
-    elif resample_method == "spline_interpolated":
+
+    if resample_method == "spline_interpolated":
         return SplineInterpolatedResampler()
-    else:
-        return None
+
+    return None
 
 
 def _template_match(observed_spectrum, template_spectrum, resample_method):
     """
-    Resample the template spectrum to match the wavelength of the observed spectrum.
-    Then, calculate chi2 on the flux of the two spectra.
+    Resample the template spectrum to match the wavelength of the observed
+    spectrum. Then, calculate chi2 on the flux of the two spectra.
 
     Parameters
     ----------
     observed_spectrum : :class:`~specutils.Spectrum1D`
-        the observed spectrum
+        The observed spectrum.
     template_spectrum : :class:`~specutils.Spectrum1D`
-        the template spectrum, which will be resampled to match the wavelength of the observed_spectrum
+        The template spectrum, which will be resampled to match the wavelength
+        of the observed spectrum.
 
     Returns
     -------
     normalized_template_spectrum : :class:`~specutils.Spectrum1D`
         The normalized spectrum template.
     chi2 : `float`
-        The chi2 of the flux of the observed spectrum and the flux of the 
+        The chi2 of the flux of the observed spectrum and the flux of the
         normalized template spectrum.
     """
     # Resample template
     if _resample(resample_method) != 0:
         fluxc_resample = _resample(resample_method)
-        template_obswavelength = fluxc_resample(template_spectrum, 
+        template_obswavelength = fluxc_resample(template_spectrum,
                                                 observed_spectrum.wavelength)
 
     # Normalize spectra
@@ -94,7 +101,7 @@ def _template_match(observed_spectrum, template_spectrum, resample_method):
     result = (num/denom)**2
     chi2 = np.sum(result.value)
 
-    # Create normalized template spectrum, which will be returned with 
+    # Create normalized template spectrum, which will be returned with
     # corresponding chi2
     normalized_template_spectrum = Spectrum1D(
         spectral_axis=template_spectrum.spectral_axis,
@@ -103,10 +110,11 @@ def _template_match(observed_spectrum, template_spectrum, resample_method):
     return normalized_template_spectrum, chi2
 
 
-def template_match(observed_spectrum, spectral_templates, resample_method="flux_conserving"):
+def template_match(observed_spectrum, spectral_templates,
+                   resample_method="flux_conserving"):
     """
-    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 
+    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
     first template is returned.
 
     Parameters
@@ -114,9 +122,9 @@ def template_match(observed_spectrum, spectral_templates, resample_method="flux_
     observed_spectrum : :class:`~specutils.Spectrum1D`
         The observed spectrum.
     spectral_templates : :class:`~specutils.Spectrum1D` or :class:`~specutils.SpectrumCollection` or `list`
-        That will give a single :class:`~specutils.Spectrum1D` when iterated 
-        over. The template spectra, which will be resampled, normalized, and 
-        compared to the observed spectrum, where the smallest chi2 and 
+        That will give a single :class:`~specutils.Spectrum1D` when iterated
+        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.
 
     Returns
@@ -124,7 +132,7 @@ def template_match(observed_spectrum, spectral_templates, resample_method="flux_
     normalized_template_spectrum : :class:`~specutils.Spectrum1D`
         The template spectrum that has been normalized.
     chi2 : `float`
-        The chi2 of the flux of the observed_spectrum and the flux of the 
+        The chi2 of the flux of the observed_spectrum and the flux of the
         normalized template spectrum.
     smallest_chi_index : `int`
         The index of the spectrum with the smallest chi2 in spectral templates.
@@ -136,8 +144,8 @@ def template_match(observed_spectrum, spectral_templates, resample_method="flux_
         return normalized_spectral_template, chi2
 
     # At this point, the template spectrum is either a ``SpectrumCollection``
-    # or a multi-dimensional``Spectrum1D``. Loop through the object and return 
-    # the template spectrum with the lowest chi square and its corresponding 
+    # or a multi-dimensional``Spectrum1D``. Loop through the object and return
+    # the template spectrum with the lowest chi square and its corresponding
     # chi square.
     chi2_min = None
     smallest_chi_spec = None