Repository: spacetelescope/jwst

Issue:
combine_1d unit tests
Issue [JP-577](https://jira.stsci.edu/browse/JP-577) was created by Nadia Dencheva:

Write unit tests for combine_1d
 * check wavelength binning is working correctly
 * check wavelength range is complete
 * check number of bins is recorded correctly

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 13fddca8b..f2eed05da 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,6 +1,12 @@
 0.13.3 (Unreleased)
 ===================
 
+
+combine_1d
+----------
+
+- Unit tests were added to combine_1d.  [#3490]
+
 ami
 ---
 
diff --git a/docs/jwst/combine_1d/arguments.rst b/docs/jwst/combine_1d/arguments.rst
index 1c6ef95bc..3d45398af 100644
--- a/docs/jwst/combine_1d/arguments.rst
+++ b/docs/jwst/combine_1d/arguments.rst
@@ -14,9 +14,3 @@ the exposure time (FITS keyword EFFEXPTM) is used as the weight.  If
 the string is "unit_weight" or "unit weight", the same weight (1) will
 be used for all input spectra.  If the string is anything else, a warning
 will be logged and unit weight will be used.
-
-*  ``--background``
-
-Set this boolean parameter to True if the input spectra are all just
-background.  The default is False, meaning that the input are science
-spectra.
diff --git a/docs/jwst/combine_1d/description.rst b/docs/jwst/combine_1d/description.rst
index 58d815401..23e51136c 100644
--- a/docs/jwst/combine_1d/description.rst
+++ b/docs/jwst/combine_1d/description.rst
@@ -12,15 +12,7 @@ spectra have been included, the output is normalized by dividing by
 the sum of the weights.
 
 The weight will typically be the integration time or the exposure time,
-but uniform (unit) weighting can be specified instead.  It is the net
-count rate that uses this weight; that is, the net count rate is multiplied
-by the integration time to get net counts, and it is the net counts that
-are added together and finally divided by the sum of the integration
-times.  The flux weighted by an additional factor of the instrumental
-sensitivity, count rate per unit flux.  The idea is that the quantity
-that is added up should be in units of counts.  If unit weight was
-specified, however, unit weight will be used for both flux and net.
-The data quality (DQ) columns will be combined using bitwise OR.
+but uniform (unit) weighting can be specified instead.
 
 The only part of this step that is not completely straightforward is the
 determination of wavelengths for the output spectrum.  The output
@@ -31,11 +23,9 @@ wavelength would be computed as the average of the wavelengths at the same
 pixel in all the input files.  The combine_1d step is intended to handle a
 more general case where the input wavelength arrays may be offset with
 respect to each other, or they might not align well due to different
-distortions.
-
-All the input wavelength arrays will be concatenated and then sorted.
-The code then looks for "clumps" in wavelength, based on the standard
-deviation of a slice of the concatenated and sorted array of input
+distortions.  All the input wavelength arrays will be concatenated and then
+sorted.  The code then looks for "clumps" in wavelength, based on the
+standard deviation of a slice of the concatenated and sorted array of input
 wavelengths; a small standard deviation implies a clump.  In regions of
 the spectrum where the input wavelengths overlap with somewhat random
 offsets and don't form any clumps, the output wavelengths are computed
@@ -46,9 +36,9 @@ Input
 =====
 An association file specifies which file or files to read for the input
 data.  Each input data file contains one or more 1-D spectra in table
-format, e.g. as written by the extract_1d step.  An input data file can
-be either SpecModel (for one spectrum) or MultiSpecModel format (which
-can contain more than one spectrum).
+format, e.g. as written by the extract_1d step.  Each input data file will
+ordinarily be in MultiSpecModel format (which can contain more than one
+spectrum).
 
 The association file should have an object called "products", which is
 a one-element list containing a dictionary.  This dictionary contains two
@@ -60,8 +50,8 @@ of which contains one input file name, identified by key "expname".
 Output
 ======
 The output will be in CombinedSpecModel format, with a table extension
-having the name COMBINE1D.  This extension will have six columns, giving
-the wavelength, flux, error estimate for the flux, net countrate in
-counts/second, the sum of the weights that were used when combining the
-input spectra, and the number of input spectra that contributed to each
-output pixel.
+having the name COMBINE1D.  This extension will have eight columns, giving
+the wavelength, flux, error estimate for the flux, surface brightness,
+error estimate for the surface brightness, the combined data quality flags,
+the sum of the weights that were used when combining the input spectra,
+and the number of input spectra that contributed to each output pixel.
diff --git a/jwst/combine_1d/combine1d.py b/jwst/combine_1d/combine1d.py
index 8cc08c696..22cbf0c37 100644
--- a/jwst/combine_1d/combine1d.py
+++ b/jwst/combine_1d/combine1d.py
@@ -75,14 +75,9 @@ class InputSpectrumModel:
             raise RuntimeError("Don't understand exptime_key = '%s'" %
                                exptime_key)
 
-        got_wcs = True
         try:
-            wcs = spec.meta.wcs
+            self.right_ascension[:], self.declination[:], _ = spec.meta.wcs(0., 0.)
         except AttributeError:
-            got_wcs = False
-        if got_wcs:
-            self.right_ascension[:], self.declination[:], _ = wcs(0.)
-        else:
             self.right_ascension[:] = ms.meta.target.ra
             self.declination[:] = ms.meta.target.dec
             log.warning("There is no WCS in the input.")
@@ -112,9 +107,6 @@ class OutputSpectrumModel:
         weight
         count
         wcs
-        wavelength_dtype
-        flux_dtype
-        dq_dtype
         normalized
     """
 
@@ -129,17 +121,13 @@ class OutputSpectrumModel:
         self.weight = None
         self.count = None
         self.wcs = None
-        self.wavelength_dtype = None
-        self.flux_dtype = None
-        self.dq_dtype = None
         self.normalized = False
 
     def assign_wavelengths(self, input_spectra):
         """Create an array of wavelengths to use for the output spectrum.
 
-        Take the union of all input wavelengths, then call method
-        compute_output_wl to bin wavelengths in groups of the number of
-        overlapping spectra.
+        Take the union of all input wavelengths, then call compute_output_wl
+        to bin wavelengths in groups of the number of overlapping spectra.
 
         Parameters
         ----------
@@ -147,171 +135,14 @@ class OutputSpectrumModel:
             List of input spectra.
         """
 
-        # Save these, so we'll know what data type to use for the output.
-        # The types used for accumulating sums and taking averages may not
-        # be the same as these types.
-        self.wavelength_dtype = input_spectra[0].wavelength.dtype
-        self.flux_dtype = input_spectra[0].flux.dtype
-        self.dq_dtype = input_spectra[0].dq.dtype
-
-        nwl = 0
-        for in_spec in input_spectra:
-            nwl += in_spec.nelem
-
-        # Create an array with all the input wavelengths (i.e. the union
-        # of the input wavelengths).
-        wl = np.zeros(nwl, dtype=np.float64)
-        i = 0
-        for in_spec in input_spectra:
-            nelem = in_spec.nelem
-            # Concatenate current input wavelengths to wl array.
-            wl[i:i + nelem] = in_spec.wavelength.copy()
-            i += nelem
-        wl.sort()
-
-        # count_input will be the number of input spectra that cover the
-        # corresponding wavelength in wl.
-        count_input = np.zeros(nwl, dtype=np.int64)
-        for in_spec in input_spectra:
-            input_wl = in_spec.wavelength
-            # wl0 and wl1 will be about a half pixel wider on either side
-            # of the wavelength range for the current input spectrum.
-            if input_wl[1] > input_wl[0]:       # wavelengths are increasing
-                wl0 = input_wl[0] - 0.5 * (input_wl[1] - input_wl[0])
-                wl1 = input_wl[-1] + 0.5 * (input_wl[-1] - input_wl[-2])
-            elif input_wl[1] < input_wl[0]:     # wavelengths are decreasing
-                wl0 = input_wl[-1] - 0.5 * (input_wl[-2] - input_wl[-1])
-                wl1 = input_wl[0] + 0.5 * (input_wl[0] - input_wl[1])
-            else:
-                raise RuntimeError("Wavelength increment must not be zero.")
-            temp = np.where(wl >= wl0, 1, 0)
-            temp = np.where(wl >= wl1, 0, temp)
-            count_input += temp
-            del temp
-        # This shouldn't happen.
-        if np.any(count_input <= 0.):
-            raise RuntimeError("Problem with input wavelengths.")
+        (wl, n_input_spectra) = count_input(input_spectra)
 
-        self.wavelength = self.compute_output_wl(wl, count_input)
+        self.wavelength = compute_output_wl(wl, n_input_spectra)
 
         self.wcs = create_spectral_wcs(input_spectra[0].right_ascension[0],
                                        input_spectra[0].declination[0],
                                        self.wavelength)
 
-    def compute_output_wl(self, wl, count_input):
-        """Compute output wavelengths.
-
-        In summary, the output wavelengths are computed by binning the
-        input wavelengths in groups of the number of overlapping spectra.
-        The merged and sorted input wavelengths are in `wl`.
-
-        If the wavelength arrays of the input spectra are nearly all the
-        same, the values in wl will be in clumps of N nearly identical
-        values (for N input files).  We want to average the wavelengths
-        in those clumps to get the output wavelengths.  However, if the
-        wavelength scales of the input files differ slightly, the values
-        in wl may clump in some regions but not in others.  The algorithm
-        below tries to find clumps; a clump is identified by having a
-        small standard deviation of all the wavelengths in a slice of wl,
-        compared with neighboring slices.  It is intended that this should
-        not find clumps where they're not present, i.e. in regions where
-        the wavelengths of the input spectra are not well aligned with each
-        other.  These regions will be gaps in the initial determination
-        of output wavelengths based on clumps.  In such regions, output
-        wavelengths will then be computed by just averaging groups of
-        input wavelengths to fill in the gaps.
-
-        Parameters
-        ----------
-        wl : 1-D array
-            An array containing all the wavelengths from all the input
-            spectra, sorted in increasing order.
-
-        count_input : 1-D array
-            An integer array of the same length as `wl`.  For a given
-            array index k (for example), count_input[k] is the number of
-            input spectra that cover wavelength wl[k].
-
-        Returns
-        -------
-        wavelength :  1-D array
-            Array of wavelengths for the output spectrum.
-        """
-
-        nwl = len(wl)
-
-        # sigma is an array of the standard deviation at each element
-        # of wl, over count_input elements.  A small value implies that
-        # there's a clump, i.e. several elements of wl with nearly the
-        # same wavelength.
-        sigma = np.zeros(nwl, dtype=np.float64) + 9999.
-
-        # mean_wl is the mean wavelength over the same slice of wl that we
-        # used to compute sigma.  If sigma is small enough that it looks
-        # as if there's a clump, we'll copy the mean_wl value to temp_wl
-        # to be one element of the output wavelengths.
-        mean_wl = np.zeros(nwl, dtype=np.float64) - 99.
-
-        # temp_wl has the same number of elements as wl, but we expect the
-        # array of output wavelengths to be significantly smaller, so
-        # temp_wl is initialized to a negative value as a flag.  Positive
-        # elements will be copied to the array of output wavelengths.
-        temp_wl = np.zeros(nwl, dtype=np.float64) - 99.
-
-        for k in range(nwl):
-            n = count_input[k]
-            if n == 1:
-                sigma[k] = 0.
-                mean_wl[k] = wl[k]
-                temp_wl[k] = mean_wl[k]
-            else:
-                k1 = k + n // 2 + 1
-                k0 = k1 - n
-                if k0 >= 0 and k1 <= nwl:
-                    sigma[k] = wl[k0:k1].std()
-                    mean_wl[k] = wl[k0:k1].mean()
-                    if sigma[k] == 0.:
-                        temp_wl[k] = mean_wl[k]
-
-        cutoff = 0.8
-        for k in range(nwl):
-            # If sigma[k] equals 0, temp_wl has already been assigned.
-            if sigma[k] > 0.:
-                if k == 0:
-                    if sigma[k] < cutoff * sigma[1]:
-                        temp_wl[k] = mean_wl[k]
-                elif k == nwl - 1:
-                    if sigma[k] < cutoff * sigma[nwl - 2]:
-                        temp_wl[k] = mean_wl[k]
-                else:
-                    if sigma[k] < cutoff * (sigma[k - 1] + sigma[k + 1]) / 2.:
-                        temp_wl[k] = mean_wl[k]
-
-        # Fill gaps in the output wavelengths by taking averages of the
-        # input wavelengths.  If there are n overlapping input spectra,
-        # average a block of n elements of `wl`.
-        done = False
-        i = 0
-        while not done:
-            if i >= nwl:
-                done = True
-            else:
-                n = count_input[i]
-                middle = n // 2
-                if i + middle < nwl:
-                    n = count_input[i + middle]
-                    if i + n < nwl:
-                        # The nominal range is i - 1 to i + n (inclusive) for
-                        # checking whether an output wavelength has already
-                        # been assigned.
-                        low = max(i - 1, 0)
-                        high = min(i + n + 1, nwl - 1)
-                        if temp_wl[low:high].max() <= 0.:
-                            temp_wl[i] = wl[i:i + n].mean()
-                i += n
-
-        return temp_wl[np.where(temp_wl > 0.)].copy()
-
     def accumulate_sums(self, input_spectra):
         """Compute a weighted sum of all the input spectra.
 
@@ -337,13 +168,19 @@ class OutputSpectrumModel:
             List of input spectra.
         """
 
+        # This is the data type for the output spectrum.  We'll use double
+        # precision for accumulating sums for most columns, but for the DQ
+        # array, use the correct output data type.
+        cmb_dtype = datamodels.CombinedSpecModel().spec_table.dtype
+        dq_dtype = cmb_dtype.fields["DQ"][0]
+
         nelem = self.wavelength.shape[0]
 
         self.flux = np.zeros(nelem, dtype=np.float64)
         self.error = np.zeros(nelem, dtype=np.float64)
         self.surf_bright = np.zeros(nelem, dtype=np.float64)
         self.sb_error = np.zeros(nelem, dtype=np.float64)
-        self.dq = np.zeros(nelem, dtype=self.dq_dtype)
+        self.dq = np.zeros(nelem, dtype=dq_dtype)
         self.weight = np.zeros(nelem, dtype=np.float64)
         self.count = np.zeros(nelem, dtype=np.float64)
 
@@ -415,15 +252,9 @@ class OutputSpectrumModel:
             log.warning("Data have not been divided by"
                         " the sum of the weights.")
 
-        dtype = [('WAVELENGTH', self.wavelength_dtype),
-                 ('FLUX', self.flux_dtype),
-                 ('ERROR', self.flux_dtype),
-                 ('SURF_BRIGHT', self.flux_dtype),
-                 ('SB_ERROR', self.flux_dtype),
-                 ('DQ', self.dq_dtype),
-                 ('WEIGHT', self.wavelength_dtype),
-                 ('N_INPUT', np.float64)]
+        cmb_dtype = datamodels.CombinedSpecModel().spec_table.dtype
 
+        # Note that these arrays have to be in the right order.
         data = np.array(list(zip(self.wavelength,
                                  self.flux,
                                  self.error,
@@ -431,7 +262,7 @@ class OutputSpectrumModel:
                                  self.sb_error,
                                  self.dq,
                                  self.weight,
-                                 self.count)), dtype=dtype)
+                                 self.count)), dtype=cmb_dtype)
         output_model = datamodels.CombinedSpecModel(spec_table=data)
 
         return output_model
@@ -446,12 +277,183 @@ class OutputSpectrumModel:
         self.weight = None
         self.count = None
         self.wcs = None
-        self.wavelength_dtype = None
-        self.flux_dtype = None
-        self.dq_dtype = None
         self.normalized = False
 
 
+def count_input(input_spectra):
+    """Determine the number of input spectra that cover each wavelength.
+
+    For any given input spectrum, the array of wavelengths gives the
+    wavelengths at the centers of the pixels.  In this context, the
+    expression that an input spectrum "covers" some particular wavelength
+    (typically not one of the elements in the input spectrum's array of
+    wavelengths) means that this wavelength is within the interval between
+    the left edge of the first pixel and the right edge of the last pixel
+    of the input spectrum.
+
+    Parameters
+    ----------
+    input_spectra : list of InputSpectrumModel objects
+        List of input spectra.
+
+    Returns
+    -------
+    wl : ndarray
+        Sorted list of all the wavelengths in all the input spectra.
+
+    n_input_spectra : ndarray, 1-D
+        For each element of `wl`, the corresponding element of
+        `n_input_spectra` is the number of input spectra that cover the
+        wavelength in `wl`.
+    """
+
+    # Create an array with all the input wavelengths (i.e. the union
+    # of the input wavelengths).
+    wl = np.hstack([in_spec.wavelength for in_spec in input_spectra])
+    wl.sort()
+    nwl = len(wl)
+
+    # n_input_spectra will be the number of input spectra that cover the
+    # corresponding wavelength in wl.
+    n_input_spectra = np.zeros(nwl, dtype=np.int64)
+    for in_spec in input_spectra:
+        input_wl = in_spec.wavelength
+        # wl0 and wl1 will be about a half pixel wider on either side
+        # of the wavelength range for the current input spectrum.
+        if input_wl[1] > input_wl[0]:       # wavelengths are increasing
+            wl0 = input_wl[0] - 0.5 * (input_wl[1] - input_wl[0])
+            wl1 = input_wl[-1] + 0.5 * (input_wl[-1] - input_wl[-2])
+        elif input_wl[1] < input_wl[0]:     # wavelengths are decreasing
+            wl0 = input_wl[-1] - 0.5 * (input_wl[-2] - input_wl[-1])
+            wl1 = input_wl[0] + 0.5 * (input_wl[0] - input_wl[1])
+        else:
+            raise RuntimeError("Wavelength increment must not be zero.")
+        temp = np.where(wl >= wl0, 1, 0)
+        temp = np.where(wl >= wl1, 0, temp)
+        n_input_spectra += temp
+        del temp
+    # This shouldn't happen.
+    if np.any(n_input_spectra <= 0.):
+        raise RuntimeError("Problem with input wavelengths.")
+
+    return (wl, n_input_spectra)
+
+
+def compute_output_wl(wl, n_input_spectra):
+    """Compute output wavelengths.
+
+    In summary, the output wavelengths are computed by binning the
+    input wavelengths in groups of the number of overlapping spectra.
+    The merged and sorted input wavelengths are in `wl`.
+
+    If the wavelength arrays of the input spectra are nearly all the
+    same, the values in wl will be in clumps of N nearly identical
+    values (for N input files).  We want to average the wavelengths
+    in those clumps to get the output wavelengths.  However, if the
+    wavelength scales of the input files differ slightly, the values
+    in wl may clump in some regions but not in others.  The algorithm
+    below tries to find clumps; a clump is identified by having a
+    small standard deviation of all the wavelengths in a slice of wl,
+    compared with neighboring slices.  It is intended that this should
+    not find clumps where they're not present, i.e. in regions where
+    the wavelengths of the input spectra are not well aligned with each
+    other.  These regions will be gaps in the initial determination
+    of output wavelengths based on clumps.  In such regions, output
+    wavelengths will then be computed by just averaging groups of
+    input wavelengths to fill in the gaps.
+
+    Parameters
+    ----------
+    wl : 1-D array
+        An array containing all the wavelengths from all the input
+        spectra, sorted in increasing order.
+
+    n_input_spectra : 1-D array
+        An integer array of the same length as `wl`.  For a given
+        array index k (for example), n_input_spectra[k] is the number of
+        input spectra that cover wavelength wl[k].
+
+    Returns
+    -------
+    wavelength :  1-D array
+        Array of wavelengths for the output spectrum.
+    """
+
+    nwl = len(wl)
+
+    # sigma is an array of the standard deviation at each element
+    # of wl, over n_input_spectra elements.  A small value implies that
+    # there's a clump, i.e. several elements of wl with nearly the
+    # same wavelength.
+    sigma = np.zeros(nwl, dtype=np.float64) + 9999.
+
+    # mean_wl is the mean wavelength over the same slice of wl that we
+    # used to compute sigma.  If sigma is small enough that it looks
+    # as if there's a clump, we'll copy the mean_wl value to temp_wl
+    # to be one element of the output wavelengths.
+    mean_wl = np.zeros(nwl, dtype=np.float64) - 99.
+
+    # temp_wl has the same number of elements as wl, but we expect the
+    # array of output wavelengths to be significantly smaller, so
+    # temp_wl is initialized to a negative value as a flag.  Positive
+    # elements will be copied to the array of output wavelengths.
+    temp_wl = np.zeros(nwl, dtype=np.float64) - 99.
+
+    for k in range(nwl):
+        n = n_input_spectra[k]
+        if n == 1:
+            sigma[k] = 0.
+            mean_wl[k] = wl[k]
+            temp_wl[k] = mean_wl[k]
+        else:
+            k1 = k + n // 2 + 1
+            k0 = k1 - n
+            if k0 >= 0 and k1 <= nwl:
+                sigma[k] = wl[k0:k1].std()
+                mean_wl[k] = wl[k0:k1].mean()
+                if sigma[k] == 0.:
+                    temp_wl[k] = mean_wl[k]
+
+    cutoff = 0.8
+    for k in range(nwl):
+        # If sigma[k] equals 0, temp_wl has already been assigned.
+        if sigma[k] > 0.:
+            if k == 0:
+                if sigma[k] < cutoff * sigma[1]:
+                    temp_wl[k] = mean_wl[k]
+            elif k == nwl - 1:
+                if sigma[k] < cutoff * sigma[nwl - 2]:
+                    temp_wl[k] = mean_wl[k]
+            else:
+                if sigma[k] < cutoff * (sigma[k - 1] + sigma[k + 1]) / 2.:
+                    temp_wl[k] = mean_wl[k]
+
+    # Fill gaps in the output wavelengths by taking averages of the
+    # input wavelengths.  If there are n overlapping input spectra,
+    # average a block of n elements of `wl`.
+    done = False
+    i = 0
+    while not done:
+        if i >= nwl:
+            done = True
+        else:
+            n = n_input_spectra[i]
+            middle = n // 2
+            if i + middle < nwl:
+                n = n_input_spectra[i + middle]
+                if i + n < nwl:
+                    # The nominal range is i - 1 to i + n (inclusive) for
+                    # checking whether an output wavelength has already
+                    # been assigned.
+                    low = max(i - 1, 0)
+                    high = min(i + n + 1, nwl - 1)
+                    if temp_wl[low:high].max() <= 0.:
+                        temp_wl[i] = wl[i:i + n].mean()
+            i += n
+
+    return temp_wl[np.where(temp_wl > 0.)].copy()
+
+
 def check_exptime(exptime_key):
     """Check exptime_key for validity.
 


Test patch:
diff --git a/jwst/combine_1d/tests/test_count_input.py b/jwst/combine_1d/tests/test_count_input.py
new file mode 100644
index 000000000..27e4a6794
--- /dev/null
+++ b/jwst/combine_1d/tests/test_count_input.py
@@ -0,0 +1,124 @@
+"""
+Test for combine1d.count_input
+"""
+import numpy as np
+
+# from jwst import datamodels
+from jwst.combine_1d import combine1d
+
+class DummySpectra:
+    def __init__(self, wavelength):
+        # The argument to count_input is a list of InputSpectrumModel
+        # objects, but the only attribute we need is wavelength, so this
+        # DummySpectra class will suffice.
+        self.wavelength = wavelength.copy()
+
+
+def test_1a():
+    input_spectra = create_input_spectra_1()
+
+    (wl, n_input_spectra) = combine1d.count_input(input_spectra)
+
+    truth_1a = truth_array_1a()
+
+    # Compare wl with the correct values.
+    assert np.allclose(wl, truth_1a, rtol=1.e-10)
+
+
+def test_1b():
+    input_spectra = create_input_spectra_1()
+    (wl, n_input_spectra) = combine1d.count_input(input_spectra)
+
+    truth_1b = truth_array_1b()
+
+    # Compare n_input_spectra with the correct values.
+    assert np.allclose(n_input_spectra, truth_1b, rtol=1.e-10)
+
+
+def test_1c():
+    input_spectra = create_input_spectra_1()
+    (wl, n_input_spectra) = combine1d.count_input(input_spectra)
+
+    output_wl = combine1d.compute_output_wl(wl, n_input_spectra)
+
+    truth_1c = truth_array_1c()
+
+    assert np.allclose(output_wl, truth_1c, rtol=1.e-10)
+
+
+def create_input_spectra_1():
+
+    input_spectra = []
+
+    # This is for the case that there are four input spectra, each one
+    # six pixels long.  The wavelengths are roughly the same in each
+    # spectrum.  The differences in wavelength between spectra is small
+    # compared to the spacing between wavelengths in any one spectrum.
+    # That is, the union of the wavelengths in all the input spectra
+    # will be clumped.  For this test, first specify the wavelengths in
+    # these clumps.  Then assign the wavelengths for each input spectrum
+    # by picking one element from each clump.
+    clumps = np.array([[1.372, 1.385, 1.391, 1.402],
+                       [1.701, 1.718, 1.719, 1.723],
+                       [1.979, 1.984, 1.991, 2.013],
+                       [2.212, 2.223, 2.228, 2.310],
+                       [2.475, 2.479, 2.487, 2.496],
+                       [2.716, 2.723, 2.729, 2.741]], dtype=np.float64)
+
+    spec = [clumps[0, 3], clumps[1, 0], clumps[2, 2], clumps[3, 1], clumps[4, 2], clumps[5, 1]]
+    input_spectra.append(DummySpectra(np.array(spec)))
+    spec = [clumps[0, 0], clumps[1, 3], clumps[2, 0], clumps[3, 3], clumps[4, 0], clumps[5, 2]]
+    input_spectra.append(DummySpectra(np.array(spec)))
+    spec = [clumps[0, 2], clumps[1, 1], clumps[2, 1], clumps[3, 2], clumps[4, 3], clumps[5, 0]]
+    input_spectra.append(DummySpectra(np.array(spec)))
+    spec = [clumps[0, 1], clumps[1, 2], clumps[2, 3], clumps[3, 0], clumps[4, 1], clumps[5, 3]]
+    input_spectra.append(DummySpectra(np.array(spec)))
+
+    return input_spectra
+
+
+def truth_array_1a():
+    """Merged and sorted wavelengths from all input files."""
+
+    # Wavelengths (merged and sorted) in input files.
+    return np.array([1.372, 1.385, 1.391, 1.402, 1.701, 1.718, 1.719, 1.723,
+                     1.979, 1.984, 1.991, 2.013, 2.212, 2.223, 2.228, 2.31,
+                     2.475, 2.479, 2.487, 2.496, 2.716, 2.723, 2.729, 2.741],
+                    dtype=np.float64)
+
+
+def truth_array_1b():
+    """Number of input models that cover each of the wavelengths in the
+       combined array.
+
+    The first element in the wl array (see function truth_array_1a) has
+    value 1.372, and that's not actually within the array of wavelengths
+    for the first input model, which are (see function
+    create_input_spectra_1):
+        [1.402 1.701 1.991 2.223 2.487 2.723]
+    But the array that we'll return has 4 (the number of input models)
+    for every element.  How can that be correct?  It's because we assume
+    the pixels in the input models have finite width, but the wavelengths
+    in the input models are the values at the centers of the pixels.  So
+    the range of wavelengths that are actually covered by an input model
+    is not just the interval between the first and last wavelengths in
+    the input model's wavelength array, it's the range from the wavelength
+    at the left edge of the first pixel to the wavelength at the right edge
+    of the last pixel.
+    """
+
+    # n_input_spectra
+    return np.array([4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+                     4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], dtype=np.int64)
+
+
+def truth_array_1c():
+    """Computed wavelengths for the output.
+
+    These are averages of the input wavelengths within clusters, sets
+    of input wavelengths that are closer together than the spacing
+    between pixels.
+    """
+
+    return np.array([1.3875, 1.71525, 1.99175, 2.24325, 2.48425, 2.72725],
+                    dtype=np.float64)