Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make fit_1D behave consistently for data with all (or all but 1) rows masked #260

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
+133 −54
Open

Make fit_1D behave consistently for data with all (or all but 1) rows masked #260

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d7b4c3b
Handle some corner cases noted by csimpson: fully-masked input array …
jehturner May 8, 2021
04b9fb9
Fix condition to allow plotting row/column 0.
jehturner May 8, 2021
15800d6
Make model attribute do something that's probably closer to the inten…
jehturner May 8, 2021
ff879da
Merge branch 'master' into fit_1D_fixes
chris-simpson Oct 27, 2021
29d1e80
fix indentation issue that got messed up in the merge
chris-simpson Oct 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
136 changes: 82 additions & 54 deletions gempy/library/fitting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -275,44 +279,56 @@ 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]
)
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
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]
)
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"]
Expand Down Expand Up @@ -372,7 +388,7 @@ def _fit(self, image, weights=None, plot=False):
self.rms = rms if rms is not np.ma.masked else np.nan

# 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
Expand Down Expand Up @@ -446,7 +462,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

Expand All @@ -470,20 +486,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)
Expand All @@ -505,9 +526,16 @@ def model(self):
or `astromodels.UnivariateSplineWithOutlierRemoval:`)
representing the model fit
"""
# 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("Cannot provide model property if there are "
"greater than one models.")
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
Expand Down
51 changes: 51 additions & 0 deletions gempy/library/tests/test_fitting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
"""
Expand Down