Add option to have inclusive upper boundaries in fitting/excision. (#…

…1171)

* Add option to make the excluded regions in fitting and excision have inclusive upper bounds (so the upper bound is...excluded)

* Add note in narrative docs about this option

* Add simple true_exciser test

docs/spectral_regions.rst

@@ -291,8 +291,13 @@ finding functions:
       (7.690954773869347 um, 8.690954773869347 um)
 
 This can be fed into the ``exclude_regions`` argument of the `~specutils.fitting.fit_generic_continuum` or
-`~specutils.fitting.fit_continuum` functions to avoid fitting regions that contain line features. Conversely, users can
-also invert the spectral region
+`~specutils.fitting.fit_continuum` functions to avoid fitting regions that contain line features. Note that,
+by default, this uses pythonic slicing, i.e., spectral values greater than or equal to the lower bound and
+less than the upper bound of the region will be excluded from the fit. For convenience in some cases, the
+``exclude_region_upper_bounds`` keyword can be set to ``True`` to exlude spectral values less than *or equal*
+to the upper bound instead.
+
+Conversely, users can also invert the spectral region
 
 .. code-block:: python
 

specutils/fitting/continuum.py

@@ -57,7 +57,8 @@ def fit_generic_continuum(spectrum, median_window=3, model=Chebyshev1D(3),
 
 
 def fit_continuum(spectrum, model=Chebyshev1D(3), fitter=LevMarLSQFitter(),
-                  exclude_regions=None, window=None, weights=None):
+                  exclude_regions=None, exclude_region_upper_bounds=False,
+                  window=None, weights=None):
     """
     Entry point for fitting using the `~astropy.modeling.fitting`
     machinery.
@@ -74,6 +75,11 @@ def fit_continuum(spectrum, model=Chebyshev1D(3), fitter=LevMarLSQFitter(),
     exclude_regions : list of 2-tuples
         List of regions to exclude in the fitting. Passed through
         to the fitmodels routine.
+    exclude_region_upper_bounds : bool, optional
+        Set to True to override the default pythonic slicing on the excluded regions (inclusive
+        lower bound, exclusive upper bound) and remove spectral values less than or equal to
+        the upper bound of the region rather than strictly less than the upper bound. Passed
+        through to the fitmodels routine.
     window : tuple of wavelengths
         Start and end wavelengths used for fitting.
     weights : list  (NOT IMPLEMENTED YET)
@@ -93,7 +99,7 @@ def fit_continuum(spectrum, model=Chebyshev1D(3), fitter=LevMarLSQFitter(),
 
     # Fit the flux to the model
     continuum_spectrum = fit_lines(spectrum, model, fitter, exclude_regions,
-                                   weights, w)
+                                   exclude_region_upper_bounds, weights, w)
 
     return continuum_spectrum
 

specutils/fitting/fitmodels.py

@@ -259,8 +259,8 @@ def _generate_line_list_table(spectrum, emission_inds, absorption_inds):
 
 
 def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties=True),
-              exclude_regions=None, weights=None, window=None, get_fit_info=False,
-              **kwargs):
+              exclude_regions=None, exclude_region_upper_bounds=False, weights=None,
+              window=None, get_fit_info=False, **kwargs):
     """
     Fit the input models to the spectrum. The parameter values of the
     input models will be used as the initial conditions for the fit.
@@ -275,8 +275,12 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties
         Fitter instance to be used when fitting model to spectrum.
     exclude_regions : list of `~specutils.SpectralRegion`
         List of regions to exclude in the fitting.
+    exclude_region_upper_bounds : bool, optional
+        Set to True to override the default pythonic slicing on the excluded regions (inclusive
+        lower bound, exclusive upper bound) and remove spectral values less than or equal to
+        the upper bound of the region rather than strictly less than the upper bound.
     weights : array-like or 'unc', optional
-        If 'unc', the unceratinties from the spectrum object are used to
+        If 'unc', the uncertainties from the spectrum object are used to
         to calculate the weights. If array-like, represents the weights to
         use in the fitting.  Note that if a mask is present on the spectrum, it
         will be applied to the ``weights`` as it would be to the spectrum
@@ -310,7 +314,8 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties
     #
 
     if exclude_regions is not None:
-        spectrum = excise_regions(spectrum, exclude_regions)
+        spectrum = excise_regions(spectrum, exclude_regions,
+                                  inclusive_upper=exclude_region_upper_bounds)
 
     #
     # Make the model a list if not already
@@ -358,6 +363,7 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties
             model_guess,
             fitter=fitter,
             exclude_regions=exclude_regions,
+            exclude_region_upper_bounds=exclude_region_upper_bounds,
             weights=weights,
             window=model_window,
             get_fit_info=get_fit_info,
@@ -377,8 +383,8 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties
 
 
 def _fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties=True),
-               exclude_regions=None, weights=None, window=None, get_fit_info=False,
-               ignore_units=False, **kwargs):
+               exclude_regions=None, exclude_region_upper_bounds=False, weights=None,
+               window=None, get_fit_info=False, ignore_units=False, **kwargs):
     """
     Fit the input model (initial conditions) to the spectrum.  Output will be
     the same model with the parameters set based on the fitting.
@@ -390,7 +396,8 @@ def _fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertaintie
     #
 
     if exclude_regions is not None:
-        spectrum = excise_regions(spectrum, exclude_regions)
+        spectrum = excise_regions(spectrum, exclude_regions,
+                                  inclusive_upper=exclude_region_upper_bounds)
 
     if isinstance(weights, str):
         if weights == 'unc':

specutils/manipulation/utils.py

@@ -8,7 +8,7 @@
 __all__ = ['excise_regions', 'linear_exciser', 'spectrum_from_model']
 
 
-def true_exciser(spectrum, region):
+def true_exciser(spectrum, region, inclusive_upper=False):
     """
     Basic spectral excise method where the array elements in the spectral
     region defined by the parameter ``region`` (a `~specutils.SpectralRegion`)
@@ -26,6 +26,10 @@ def true_exciser(spectrum, region):
     region : `~specutils.SpectralRegion`
         The region of the spectrum to remove.
 
+    inclusive_upper : bool, optional
+        Set to True to override the default pythonic slicing on the regions (inclusive
+        lower bound, exclusive upper bound) to include the specified upper bound instead.
+
     Returns
     -------
     spectrum : `~specutils.Spectrum1D`
@@ -43,7 +47,11 @@ def true_exciser(spectrum, region):
 
     for subregion in region:
         # Find the indices of the spectral_axis array corresponding to the subregion
-        region_mask = (spectral_axis >= region.lower) & (spectral_axis < region.upper)
+        if inclusive_upper:
+            region_mask = (spectral_axis >= subregion.lower) & (spectral_axis <= subregion.upper)
+        else:
+            region_mask = (spectral_axis >= subregion.lower) & (spectral_axis < subregion.upper)
+
         temp_indices = np.nonzero(region_mask)[0]
         if excise_indices is None:
             excise_indices = temp_indices
@@ -76,7 +84,7 @@ def true_exciser(spectrum, region):
                       radial_velocity=spectrum.radial_velocity if not isinstance(new_spectral_axis, SpectralCoord) else None)
 
 
-def linear_exciser(spectrum, region):
+def linear_exciser(spectrum, region, inclusive_upper=False):
     """
     Basic spectral excise method where the spectral region defined by the
     parameter ``region`` (a `~specutils.SpectralRegion`) will result
@@ -125,7 +133,11 @@ def linear_exciser(spectrum, region):
 
     for subregion in region:
         # Find the indices of the spectral_axis array corresponding to the subregion
-        region_mask = (spectral_axis >= subregion.lower) & (spectral_axis < subregion.upper)
+        if inclusive_upper:
+            region_mask = (spectral_axis >= subregion.lower) & (spectral_axis <= subregion.upper)
+        else:
+            region_mask = (spectral_axis >= subregion.lower) & (spectral_axis < subregion.upper)
+
         inclusive_indices = np.nonzero(region_mask)[0]
         # Now set the flux values for these indices to be a
         # linear range
@@ -150,7 +162,7 @@ def linear_exciser(spectrum, region):
                       radial_velocity=spectrum.radial_velocity if not isinstance(spectral_axis, SpectralCoord) else None)
 
 
-def excise_regions(spectrum, regions, exciser=true_exciser):
+def excise_regions(spectrum, regions, exciser=true_exciser, inclusive_upper=False):
     """
     Method to remove or replace the flux in the defined regions of the spectrum
     depending on the function provided in the ``exciser`` argument.
@@ -173,6 +185,10 @@ def excise_regions(spectrum, regions, exciser=true_exciser):
         methods could be defined and used by this routine.
         default: true_exciser
 
+    inclusive_upper : bool, optional
+        Set to True to override the default pythonic slicing on the regions (inclusive
+        lower bound, exclusive upper bound) to include the specified upper bound instead.
+
     Returns
     -------
     spectrum : `~specutils.Spectrum1D`
@@ -190,12 +206,12 @@ def excise_regions(spectrum, regions, exciser=true_exciser):
         raise ValueError('The spectrum parameter must be Spectrum1D object.')
 
     for region in regions:
-        spectrum = excise_region(spectrum, region, exciser)
+        spectrum = excise_region(spectrum, region, exciser, inclusive_upper=inclusive_upper)
 
     return spectrum
 
 
-def excise_region(spectrum, region, exciser=true_exciser):
+def excise_region(spectrum, region, exciser=true_exciser, inclusive_upper=False):
     """
     Method to remove or replace the flux in the defined regions of the spectrum
     depending on the function provided in the ``exciser`` argument.
@@ -217,6 +233,10 @@ def excise_region(spectrum, region, exciser=true_exciser):
         methods could be defined and used by this routine.
         default: true_exciser
 
+    inclusive_upper : bool, optional
+        Set to True to override the default pythonic slicing on the regions (inclusive
+        lower bound, exclusive upper bound) to include the specified upper bound instead.
+
     Returns
     -------
     spectrum : `~specutils.Spectrum1D`
@@ -240,7 +260,7 @@ def excise_region(spectrum, region, exciser=true_exciser):
     #  Call the exciser method
     #
 
-    return exciser(spectrum, region)
+    return exciser(spectrum, region, inclusive_upper=inclusive_upper)
 
 
 def spectrum_from_model(model_input, spectrum):

specutils/tests/test_manipulation.py

@@ -11,7 +11,7 @@
 
 def test_true_exciser():
     np.random.seed(84)
-    spectral_axis = np.linspace(5000,5100,num=100)*u.AA
+    spectral_axis = np.linspace(5000, 5099, num=100)*u.AA
     flux = (np.random.randn(100) + 3) * u.Jy
     spec = Spectrum1D(flux=flux, spectral_axis=spectral_axis)
     region = SpectralRegion([(5005,5010), (5060,5065)]*u.AA)
@@ -21,10 +21,14 @@ def test_true_exciser():
     assert len(excised_spec.flux) == len(spec.flux)-10
     assert np.isclose(excised_spec.flux.sum(), 243.2617*u.Jy, atol=0.001*u.Jy)
 
+    excised_spec = excise_regions(spec, region, inclusive_upper=True)
+    assert len(excised_spec.spectral_axis) == len(spec.spectral_axis)-12
+    assert len(excised_spec.flux) == len(spec.flux)-12
+
 
 def test_linear_exciser():
     np.random.seed(84)
-    spectral_axis = np.linspace(5000,5100,num=100)*u.AA
+    spectral_axis = np.linspace(5000,5099,num=100)*u.AA
     flux = (np.random.rand(100)*100) * u.Jy
     spec = Spectrum1D(flux=flux, spectral_axis = spectral_axis)
     region = SpectralRegion([(5020,5030)]*u.AA)