From 66adb92b48ea38d81a13251395cee94bc601e2ba Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Thu, 29 Dec 2022 12:28:38 -0500 Subject: [PATCH] model fitting tweaks and cleanup (#1947) * simplify caching in model fitting * estimate slope and intercept for Linear1D --- CHANGES.rst | 4 ++ .../plugins/model_fitting/initializers.py | 12 ++--- .../plugins/model_fitting/model_fitting.py | 54 ++++++------------- .../model_fitting/tests/test_plugin.py | 16 +++--- 4 files changed, 31 insertions(+), 55 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 02d07caeeb..e81644e1d6 100644 --- a/CHANGES.rst +++ b/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 ^^^^^^^ diff --git a/jdaviz/configs/default/plugins/model_fitting/initializers.py b/jdaviz/configs/default/plugins/model_fitting/initializers.py index f288074cbb..c7aa327e71 100644 --- a/jdaviz/configs/default/plugins/model_fitting/initializers.py +++ b/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 diff --git a/jdaviz/configs/default/plugins/model_fitting/model_fitting.py b/jdaviz/configs/default/plugins/model_fitting/model_fitting.py index 7f3711130f..f8f94395c6 100644 --- a/jdaviz/configs/default/plugins/model_fitting/model_fitting.py +++ b/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,7 +739,7 @@ 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( @@ -773,7 +747,7 @@ def _fit_model_to_cube(self, add_data): 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 diff --git a/jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py b/jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py index 7650efcc12..148388d173 100644 --- a/jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py +++ b/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))