model fitting tweaks and cleanup (#1947)

* simplify caching in model fitting
* estimate slope and intercept for Linear1D

CHANGES.rst

@@ -14,6 +14,8 @@ New Features
 - Resetting viewer limits (via ``reset_limits`` or the zoom home button) now accounts for all visible
   data layers instead of just the reference data. [#1897]
 
+- Linear1D model component now estimates slope and intercept. [#1947]
+
 Cubeviz
 ^^^^^^^
 
@@ -92,6 +94,8 @@ Bug Fixes
 
 - Console logging is restored for "Desktop Mode" Windows users. [#1887]
 
+- Model fitting initial estimates now respect selected subset. [#1947]
+
 Cubeviz
 ^^^^^^^
 

jdaviz/configs/default/plugins/model_fitting/initializers.py

@@ -77,16 +77,10 @@ def initialize(self, instance, x, y):
         instance : `~astropy.modeling.Model`
             The initialized model.
         """
+        slope, intercept = np.polynomial.Polynomial.fit(x.value.flatten(), y.value.flatten(), 1)
 
-        # y_range = np.max(y) - np.min(y)
-        # x_range = x[-1] - x[0]
-        # slope = y_range / x_range
-        # y0 = y[0]
-
-        y_mean = np.mean(y)
-
-        instance.slope.value = 0.0
-        instance.intercept.value = y_mean.value
+        instance.slope.value = slope
+        instance.intercept.value = intercept
 
         return instance
 

jdaviz/configs/default/plugins/model_fitting/model_fitting.py

@@ -1,5 +1,6 @@
 import re
 import numpy as np
+from copy import deepcopy
 
 import astropy.units as u
 from specutils import Spectrum1D
@@ -106,7 +107,6 @@ class ModelFitting(PluginTemplateMixin, DatasetSelectMixin,
     residuals_label_invalid_msg = Unicode('').tag(sync=True)
 
     def __init__(self, *args, **kwargs):
-        self._spectrum1d = None
         super().__init__(*args, **kwargs)
 
         self._default_spectrum_viewer_reference_name = kwargs.get(
@@ -121,8 +121,6 @@ def __init__(self, *args, **kwargs):
         self.component_models = []
         self._initialized_models = {}
         self._display_order = False
-        self._window = None
-        self._original_mask = None
         if self.config == "cubeviz":
             self.spatial_subset = SubsetSelect(self,
                                                'spatial_subset_items',
@@ -328,31 +326,11 @@ def _dataset_selected_changed(self, event=None):
         # (won't affect calculations because these locations are masked)
         selected_spec.flux[np.isnan(selected_spec.flux)] = 0.0
 
-        # Save original mask so we can reset after applying a subset mask
-        self._original_mask = selected_spec.mask
-
         self._units["x"] = str(
             selected_spec.spectral_axis.unit)
         self._units["y"] = str(
             selected_spec.flux.unit)
 
-        self._spectrum1d = selected_spec
-
-        # Also set the spectral min and max to default to the full range
-        # This is no longer needed for 1D but is preserved for now pending
-        # fixes to Cubeviz for multi-subregion subsets
-        self._window = None
-
-    @observe("spectral_subset_selected")
-    def _on_spectral_subset_selected(self, event):
-        # TODO: does this window not account for gaps?  Should we add the warning?
-        # or can this be removed (see note above in _dataset_selected_changed)
-        if self.spectral_subset_selected == "Entire Spectrum":
-            self._window = None
-        else:
-            spectral_min, spectral_max = self.spectral_subset.selected_min_max(self._spectrum1d)
-            self._window = u.Quantity([spectral_min, spectral_max])
-
     def _default_comp_label(self, model, poly_order=None):
         abbrevs = {'BlackBody': 'BB', 'PowerLaw': 'PL', 'Lorentz1D': 'Lo'}
         abbrev = abbrevs.get(model, model[0].upper())
@@ -432,9 +410,6 @@ def create_model_component(self, model_component=None, model_component_label=Non
             Will raise an error if provided and ``model_component`` is not "Polynomial1D".
             If not provided, will default according to ``poly_order``.
         """
-        if not self._spectrum1d:
-            self._dataset_selected_changed()
-
         model_comp = model_component if model_component is not None else self.model_comp_selected
 
         if model_comp != "Polynomial1D" and poly_order is not None:
@@ -490,12 +465,14 @@ def create_model_component(self, model_component=None, model_component_label=Non
 
             initial_values[param_name] = initial_val
 
+        masked_spectrum = self._apply_subset_masks(self.dataset.selected_obj, self.spectral_subset)
+        mask = masked_spectrum.mask
         initialized_model = initialize(
             MODELS[model_comp](name=comp_label,
                                **initial_values,
                                **new_model.get("model_kwargs", {})),
-            self._spectrum1d.spectral_axis,
-            self._spectrum1d.flux)
+            masked_spectrum.spectral_axis[~mask] if mask is not None else masked_spectrum.spectral_axis,  # noqa
+            masked_spectrum.flux[~mask] if mask is not None else masked_spectrum.flux)
 
         # need to loop over parameters again as the initializer may have overridden
         # the original default value
@@ -683,15 +660,15 @@ def _fit_model_to_spectrum(self, add_data):
         models_to_fit = self._reinitialize_with_fixed()
 
         # Apply masks from selected spectral subsets
-        self._apply_subset_masks(self._spectrum1d, self.spectral_subset)
+        masked_spectrum = self._apply_subset_masks(self.dataset.selected_obj, self.spectral_subset)
 
         try:
             fitted_model, fitted_spectrum = fit_model_to_spectrum(
-                self._spectrum1d,
+                masked_spectrum,
                 models_to_fit,
                 self.model_equation,
                 run_fitter=True,
-                window=self._window
+                window=None
             )
         except AttributeError:
             msg = SnackbarMessage("Unable to fit: model equation may be invalid",
@@ -709,7 +686,7 @@ def _fit_model_to_spectrum(self, add_data):
             if self.residuals_calculate:
                 # NOTE: this will NOT load into the viewer since we have already called
                 # add_results_from_plugin above.
-                self.add_results.add_results_from_plugin(self._spectrum1d-fitted_spectrum,
+                self.add_results.add_results_from_plugin(masked_spectrum-fitted_spectrum,
                                                          label=self.residuals.value,
                                                          replace=False)
 
@@ -725,11 +702,8 @@ def _fit_model_to_spectrum(self, add_data):
         # as the starting point for cube fitting
         self._update_initialized_parameters()
 
-        # Reset the data mask in case we use a different subset next time
-        self._spectrum1d.mask = self._original_mask
-
         if self.residuals_calculate:
-            return fitted_model, fitted_spectrum, self._spectrum1d-fitted_spectrum
+            return fitted_model, fitted_spectrum, masked_spectrum-fitted_spectrum
         return fitted_model, fitted_spectrum
 
     def _fit_model_to_cube(self, add_data):
@@ -765,15 +739,15 @@ def _fit_model_to_cube(self, add_data):
 
         # Apply masks from selected subsets
         for subset in [self.spatial_subset, self.spectral_subset]:
-            self._apply_subset_masks(spec, subset)
+            spec = self._apply_subset_masks(spec, subset)
 
         try:
             fitted_model, fitted_spectrum = fit_model_to_spectrum(
                 spec,
                 models_to_fit,
                 self.model_equation,
                 run_fitter=True,
-                window=self._window
+                window=None
             )
         except ValueError:
             snackbar_message = SnackbarMessage(
@@ -819,8 +793,9 @@ def _apply_subset_masks(self, spectrum, subset_component):
         """
         # only look for a mask if there is a selected subset:
         if subset_component.selected == subset_component.default_text:
-            return
+            return spectrum
 
+        spectrum = deepcopy(spectrum)
         subset_mask = subset_component.selected_subset_mask
 
         if spectrum.mask is not None:
@@ -845,3 +820,4 @@ def _apply_subset_masks(self, spectrum, subset_component):
                 subset_mask = np.swapaxes(subset_mask, 1, 0)
 
             spectrum.mask = subset_mask
+        return spectrum

jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py

@@ -223,9 +223,10 @@ def test_subset_masks(cubeviz_helper, spectrum1d_cube_larger):
 
     # Get the data object without collapsing the spectrum:
     data = cubeviz_helper.app.data_collection[0].get_object(cls=Spectrum1D, statistic=None)
-    p._apply_subset_masks(data, p.spatial_subset)
+    masked_data = p._apply_subset_masks(data, p.spatial_subset)
 
     assert data.mask is None
+    assert masked_data.mask is None
 
     # select the spatial subset, then show that mask is correct:
     assert "Subset 1" in p.spatial_subset.choices
@@ -238,11 +239,11 @@ def test_subset_masks(cubeviz_helper, spectrum1d_cube_larger):
     subset = cubeviz_helper.app.data_collection[0].get_subset_object(
         p.spatial_subset_selected, cls=Spectrum1D, statistic=None
     )
-    p._apply_subset_masks(data, p.spatial_subset)
+    masked_data = p._apply_subset_masks(data, p.spatial_subset)
     expected_spatial_mask = np.ones(data.flux.shape).astype(bool)
     expected_spatial_mask[:, :2, :] = False
 
-    assert np.all(data.mask == expected_spatial_mask)
+    assert np.all(masked_data.mask == expected_spatial_mask)
     assert np.all(subset.mask == expected_spatial_mask)
 
     sv = cubeviz_helper.app.get_viewer('spectrum-viewer')
@@ -271,12 +272,12 @@ def test_subset_masks(cubeviz_helper, spectrum1d_cube_larger):
     subset = cubeviz_helper.app.data_collection[0].get_subset_object(
         p.spectral_subset_selected, cls=Spectrum1D, statistic=None
     )
-    p._apply_subset_masks(data, p.spectral_subset)
+    masked_data = p._apply_subset_masks(data, p.spectral_subset)
 
     expected_spectral_mask = np.ones(data.flux.shape).astype(bool)
     expected_spectral_mask[:, :, 3:] = False
 
-    assert np.all(data.mask == expected_spectral_mask)
+    assert np.all(masked_data.mask == expected_spectral_mask)
     assert np.all(subset.mask == expected_spectral_mask)
 
     # Get the data object again (ensures mask == None)
@@ -285,8 +286,9 @@ def test_subset_masks(cubeviz_helper, spectrum1d_cube_larger):
     )
 
     # apply both spectral+spatial masks:
+    masked_data = data
     for subset in [p.spatial_subset, p.spectral_subset]:
-        p._apply_subset_masks(data, subset)
+        masked_data = p._apply_subset_masks(masked_data, subset)
 
     # Check that both masks are applied correctly
-    assert np.all(data.mask == (expected_spectral_mask | expected_spatial_mask))
+    assert np.all(masked_data.mask == (expected_spectral_mask | expected_spatial_mask))