From d7b4c3b842120dccca2b0369bdd1ec7e2a4d902a Mon Sep 17 00:00:00 2001 From: "James E.H. Turner" Date: Fri, 7 May 2021 20:02:18 -0400 Subject: [PATCH 1/4] Handle some corner cases noted by csimpson: fully-masked input array and >1D input with all but one row/column fully masked (ie. only 1 model to fit). --- gempy/library/fitting.py | 129 ++++++++++++++++------------ gempy/library/tests/test_fitting.py | 51 +++++++++++ 2 files changed, 126 insertions(+), 54 deletions(-) diff --git a/gempy/library/fitting.py b/gempy/library/fitting.py index 7c22f29b7..86d21408b 100644 --- a/gempy/library/fitting.py +++ b/gempy/library/fitting.py @@ -263,7 +263,11 @@ def _fit(self, image, weights=None, plot=False): # fails for "1-model sets" and it should be more efficient anyway: image_to_fit = image if image.ndim == 1: - n_models = 1 + if (image.mask is np.ma.nomask or + (~image.mask).sum() > self.order): + n_models = 1 + else: + n_models = 0 # don't try to fit a fully-masked 1D array elif image.mask is np.ma.nomask: n_models = image.shape[1] else: @@ -275,43 +279,55 @@ def _fit(self, image, weights=None, plot=False): image_to_fit = image[:, self._good_cols] if weights is not None: weights = weights[:, self._good_cols] + # remove model set axis if this leaves a single / no model + if image_to_fit.shape[1] < 2: + image_to_fit = image_to_fit.reshape((image_to_fit.size,)) + if weights is not None: + weights = weights.reshape((weights.size,)) + + # Do the fitting unless the input data are entirely masked, which + # can happen eg. when inspecting the fit for a sample image row. + if n_models == 0: + fitted_models = None - model_set = self.model_class( - degree=self.order, n_models=n_models, - domain=self.domain, - model_set_axis=(None if n_models == 1 else 1), - **self.model_args - ) - - # Configure iterative linear fitter with rejection: - fitter = fitting.FittingWithOutlierRemoval( - fitting.LinearLSQFitter(), - sigma_clip, - niter=self.niter, - # additional args are passed to outlier_func, i.e. sigma_clip - sigma_lower=self.sigma_lower, - sigma_upper=self.sigma_upper, - maxiters=1, - cenfunc='mean', - stdfunc='std', - grow=self.grow # requires AstroPy 4.2 (#10613) - ) - - # Fit the pixel data with rejection of outlying points: - fitted_models, fitted_mask = fitter( - model_set, - points[user_reg], image_to_fit[user_reg], - weights=None if weights is None else weights[user_reg] - ) - - # Incorporate mask for fitted columns into the full-sized mask: - if image.ndim > 1 and n_models < image.shape[1]: - # this is quite ugly, but seems the best way to assign to an - # array with a mask on both dimensions. This is equivalent to: - # mask[user_reg, masked_cols] = fitted_mask - mask[user_reg[:, None] & self._good_cols] = fitted_mask.flat else: - mask[user_reg] = fitted_mask + model_set = self.model_class( + degree=self.order, n_models=n_models, + domain=self.domain, + model_set_axis=(None if n_models == 1 else 1), + **self.model_args + ) + + # Configure iterative linear fitter with rejection: + fitter = fitting.FittingWithOutlierRemoval( + fitting.LinearLSQFitter(), + sigma_clip, + niter=self.niter, + # additional args are passed to outlier_func (sigma_clip) + sigma_lower=self.sigma_lower, + sigma_upper=self.sigma_upper, + maxiters=1, + cenfunc='mean', + stdfunc='std', + grow=self.grow # requires AstroPy 4.2 (#10613) + ) + + # Fit the pixel data with rejection of outlying points: + fitted_models, fitted_mask = fitter( + model_set, + points[user_reg], image_to_fit[user_reg], + weights=None if weights is None else weights[user_reg] + ) + + # Incorporate mask for fitted columns into the full-sized mask: + if image.ndim > 1 and n_models < image.shape[1]: + # this is quite ugly, but seems the best way to assign + # to an array with a mask on both dimensions. This is + # equivalent to: + # mask[user_reg, masked_cols] = fitted_mask + mask[user_reg[:, None] & self._good_cols] = fitted_mask.flat + else: + mask[user_reg] = fitted_mask else: max_order = len(points) - self.model_args["k"] @@ -443,7 +459,7 @@ def evaluate(self, points=None): input `image` to which the fit was performed along any other axes. """ - astropy_model = isinstance(self._models, Model) + astropy_model = isinstance(self._models, Model) or self._models is None tmpaxis = 0 if astropy_model else -1 @@ -467,20 +483,25 @@ def evaluate(self, points=None): # modelling always seems to return float64: fitvals = np.zeros(stack_shape, dtype=self._dtype) - # Determine the model values we want to return: - if astropy_model: - if fitvals.ndim > 1 and len(self._models) < fitvals.shape[1]: - # If we removed bad columns, we now need to fill them properly - # in the output array - fitvals[:, self._good_cols] = self._models(points, - model_set_axis=False) + # Determine the model values we want to return (if there were no good + # columns to fit models to, just keep the already-populated zeroes): + if self._models is not None: + if astropy_model: + if fitvals.ndim > 1 and len(self._models) < fitvals.shape[1]: + # If we removed bad columns, we now need to fill them + # properly in the output array + tmpvals = self._models(points, model_set_axis=False) + if tmpvals.ndim == 1: # single model for only 1 good col + tmpvals = tmpvals[:, np.newaxis] + fitvals[:, self._good_cols] = tmpvals + del tmpvals + else: + fitvals[:] = self._models(points, model_set_axis=False) else: - fitvals[:] = self._models(points, model_set_axis=False) - else: - for n, single_model in enumerate(self._models): - # Determine model values to be returned (see comment in _fit - # about discarding values stored in the spline object): - fitvals[n] = single_model(points) + for n, single_model in enumerate(self._models): + # Determine model values to be returned (see comment in + # _fit about discarding values stored in the spline object) + fitvals[n] = single_model(points) # Restore the ordering & shape of the original input array: fitvals = fitvals.reshape(tmpshape) @@ -502,9 +523,9 @@ def model(self): or `astromodels.UnivariateSplineWithOutlierRemoval:`) representing the model fit """ - if len(self._models) > 1: - raise ValueError("Cannot provide model property if there are " - "greater than one models.") + if self._models is None or len(self._models) > 1: + raise ValueError("Can only provide model property if there is a " + "single model.") astropy_model = isinstance(self._models, Model) if astropy_model: return self._models @@ -564,4 +585,4 @@ def translate_params(params): "sigma_lower": lsigma, "sigma_upper": hsigma, "regions": params.get("regions")} - return new_params \ No newline at end of file + return new_params diff --git a/gempy/library/tests/test_fitting.py b/gempy/library/tests/test_fitting.py index b4c93a4d3..492525a94 100644 --- a/gempy/library/tests/test_fitting.py +++ b/gempy/library/tests/test_fitting.py @@ -292,6 +292,57 @@ def test_cubic_spline_def_ax_ord3_masked(self): assert_equal(fit1d.mask[4:6,80:93], True) assert_equal(fit1d.mask[24:27,24:27], True) + def test_chebyshev_single_masked(self): + """ + Fit a completely-masked single 1D array (which should return zeroes + without actually fitting a model). + """ + + masked_row = np.ma.masked_array(self.data[16], mask=True) + + fit1d = fit_1D(masked_row, weights=self.weights[16], + function='chebyshev', order=4, + sigma_lower=2.5, sigma_upper=2.5, niter=5, plot=debug) + fit_vals = fit1d.evaluate() + + assert_allclose(fit_vals, 0., atol=1e-6, rtol=0.) + assert_equal(fit1d.mask, True) + + def test_chebyshev_1_unmasked_row(self): + """ + Fit object spectrum with Chebyshev polynomial, where all rows but 1 + are masked (leaving a single model to fit, rather than a set). + """ + + masked_data = np.ma.masked_array(self.data, mask=False) + masked_data.mask[1:] = True + + fit1d = fit_1D(masked_data, weights=self.weights, function='chebyshev', + order=4, sigma_lower=2.5, sigma_upper=2.5, niter=5, + plot=debug) + fit_vals = fit1d.evaluate() + + assert_allclose(fit_vals[0], (self.obj + self.bglev)[0], + atol=5., rtol=0.005) + assert_allclose(fit_vals[1:], 0., atol=1e-6, rtol=0.) + assert_equal(fit1d.mask[1:], True) + + def test_chebyshev_all_rows_masked(self): + """ + Fit a completely-masked 2D array (which should return zeroes without + actually fitting a model). + """ + + masked_data = np.ma.masked_array(self.data, mask=True) + + fit1d = fit_1D(masked_data, weights=self.weights, function='chebyshev', + order=4, sigma_lower=2.5, sigma_upper=2.5, niter=5, + plot=debug) + fit_vals = fit1d.evaluate() + + assert_allclose(fit_vals, 0., atol=1e-6, rtol=0.) + assert_equal(fit1d.mask, True) + class TestFit1DCube: """ From 04b9fb9d88632bc2dff3867ff89d9f6179dd8c78 Mon Sep 17 00:00:00 2001 From: "James E.H. Turner" Date: Fri, 7 May 2021 20:10:46 -0400 Subject: [PATCH 2/4] Fix condition to allow plotting row/column 0. --- gempy/library/fitting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gempy/library/fitting.py b/gempy/library/fitting.py index 86d21408b..7e28816ab 100644 --- a/gempy/library/fitting.py +++ b/gempy/library/fitting.py @@ -385,7 +385,7 @@ def _fit(self, image, weights=None, plot=False): self.rms = rms # Plot the fit: - if plot: + if plot is not False: self._plot(origim, index=None if plot is True else plot) # Basic plot for debugging/inspection (interactive plotting will be handled From 15800d6406747bf55300945a92b948e7aaea78b7 Mon Sep 17 00:00:00 2001 From: "James E.H. Turner" Date: Fri, 7 May 2021 21:17:56 -0400 Subject: [PATCH 3/4] Make model attribute do something that's probably closer to the intended behaviour for fully-masked data. --- gempy/library/fitting.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/gempy/library/fitting.py b/gempy/library/fitting.py index 7e28816ab..10fda48e5 100644 --- a/gempy/library/fitting.py +++ b/gempy/library/fitting.py @@ -523,7 +523,14 @@ def model(self): or `astromodels.UnivariateSplineWithOutlierRemoval:`) representing the model fit """ - if self._models is None or len(self._models) > 1: + # If no models were fitted because the input is entirely masked, return + # a model of the expected type with all zeroes, for compatibility: + if self._models is None: + return self.model_class(degree=self.order, n_models=1, + domain=self.domain, model_set_axis=None, + **self.model_args) + + if len(self._models) > 1: raise ValueError("Can only provide model property if there is a " "single model.") astropy_model = isinstance(self._models, Model) From 29d1e807c99963ecb70a97f8e358f228e8a83b94 Mon Sep 17 00:00:00 2001 From: chris-simpson Date: Tue, 26 Oct 2021 15:22:41 -1000 Subject: [PATCH 4/4] fix indentation issue that got messed up in the merge --- gempy/library/fitting.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/gempy/library/fitting.py b/gempy/library/fitting.py index 5309c0835..ef6546732 100644 --- a/gempy/library/fitting.py +++ b/gempy/library/fitting.py @@ -312,23 +312,23 @@ def _fit(self, image, weights=None, plot=False): grow=self.grow # requires AstroPy 4.2 (#10613) ) - # Fit the pixel data with rejection of outlying points: - fitted_models, fitted_mask = fitter( - model_set, - points[user_reg], image_to_fit[user_reg], - weights=None if weights is None else weights[user_reg] - ) - self.fit_info = fitter.fit_info - - # Incorporate mask for fitted columns into the full-sized mask: - if image.ndim > 1 and n_models < image.shape[1]: - # this is quite ugly, but seems the best way to assign - # to an array with a mask on both dimensions. This is - # equivalent to: - # mask[user_reg, masked_cols] = fitted_mask - mask[user_reg[:, None] & self._good_cols] = fitted_mask.flat - else: - mask[user_reg] = fitted_mask + # Fit the pixel data with rejection of outlying points: + fitted_models, fitted_mask = fitter( + model_set, + points[user_reg], image_to_fit[user_reg], + weights=None if weights is None else weights[user_reg] + ) + self.fit_info = fitter.fit_info + + # Incorporate mask for fitted columns into the full-sized mask: + if image.ndim > 1 and n_models < image.shape[1]: + # this is quite ugly, but seems the best way to assign + # to an array with a mask on both dimensions. This is + # equivalent to: + # mask[user_reg, masked_cols] = fitted_mask + mask[user_reg[:, None] & self._good_cols] = fitted_mask.flat + else: + mask[user_reg] = fitted_mask else: #max_order = len(points) - self.model_args["k"]