Commit 73db961c authored by skamann's avatar skamann
Browse files

Debugged code for usage of unresolved component and sources outside field of view.

parent ac88a87c
......@@ -57,7 +57,7 @@ parameters are not fitted.
Latest SVN revision
-------------------
357, 2016/08/10
359, 2016/08/11
"""
import contextlib
import logging
......@@ -450,7 +450,7 @@ class FitLayer(object):
# cmap=plt.cm.hot, edgecolors='face', vmin=-100, vmax=1000)
# ax2 = fig.add_axes([0.04, 0.02, 0.92, 0.46])
# ax2.scatter(self.transformation[:, 1], self.transformation[:, 0], marker='s', s=300, c=self.model.sum(axis=1).getA1(),
# cmap=plt.cm.hot, edgecolors='face', vmin=-100, vmax=1000)
# cmap=plt.cm.hot, edgecolors='face', vmin=-100, vmax=1000)
# plt.show()
# The fit for the source fluxes is carried out on the residuals, not the data itself. The fitted residual
......@@ -542,6 +542,8 @@ class FitLayer(object):
# estimate uncertainties that are returned with the data
uncertainties = 1. / (np.sqrt(self.source_matrix.multiply(self.source_matrix).sum(axis=0)).getA1())
if not included_components.all():
uncertainties = np.insert(uncertainties, np.flatnonzero(~included_components), np.nan)
# collect data that is returned as fit result
if psf_parameters_to_fit is not None:
......@@ -665,9 +667,8 @@ class FitLayer(object):
results = []
for i, psf in enumerate(self.psf_instances):
if self.last_component_is_unresolved and not update_unresolved:
if self.components[i] == self.n_components.max():
if self.components[i] == self.n_components - 1:
continue
spaxels, fluxes = psf.get_flux()
......@@ -712,10 +713,9 @@ class FitLayer(object):
# If components are entirely outside the field of view, the indices in i1 will not be contiguous but there will
# be gaps that must be removed.
if len(empty_components) > (1 - update_unresolved):
assert np.intersect1d(i1, empty_components).size == 0, "i1 contains indices of components with zero flux."
for empty_component in empty_components[::-1]:
i1[i1 > empty_component] -= 1
assert np.intersect1d(i1, empty_components).size == 0, "i1 contains indices of components with zero flux."
for empty_component in empty_components[::-1]:
i1[i1 > empty_component] -= 1
# create matrix: first, a COO matrix is constructed. Note that it handles the case when the input data contains
# multiple entries for the some indices (row, column). After construction, the matrix is converted to a
......@@ -738,7 +738,11 @@ class FitLayer(object):
_Apsf = sparse.coo_matrix((f, (i0, i1)), shape=(nrows, ncols))
if self.n_background > 0:
_Apsf = _Apsf + self.background
if empty_components is None:
_Apsf = _Apsf + self.background
else:
used = ~np.in1d(np.arange(self.n_components + self.n_background), empty_components)
_Apsf = _Apsf + self.background.tocsc()[:, used].tocoo()
t3 = time.time()
logger.debug("Time required to create matrix: {0:.3f}s".format(t3 - t2))
......@@ -1078,18 +1082,14 @@ class FitLayer(object):
raise np.linalg.linalg.LinAlgError("Invalid value '{0}' for parameter '{1}'.".format(x0[i], name))
# calculate new model
psf_matrix, included_components = self.build_matrix(update_unresolved=True) # False)
n_comp = included_components.size + self.n_background
psf_matrix, included_components = self.build_matrix(update_unresolved=False)
n_comp = included_components.sum() + self.n_background
psf_matrix.data *= np.take(np.sqrt(self.fit_weights), psf_matrix.indices)
if (~included_components).any():
fluxes = fluxes[np.hstack((included_components, np.ones(self.n_background, dtype=np.bool)))]
model = psf_matrix * sparse.dia_matrix((fluxes, [0]), shape=(n_comp, n_comp))
# print(model.sum(axis=1).getA1())
# import matplotlib.pyplot as plt
# fig = plt.figure(figsize=(8, 8))
# ax = fig.add_axes([0.12, 0.12, 0.83, 0.83])
# ax.scatter(self.transformation[:, 1], self.transformation[:, 0], marker='s', s=300, c=(self.data * np.sqrt(self.fit_weights)) - model.sum(axis=1).getA1(),
# cmap=plt.cm.hot, edgecolors='face')
# plt.show()
# return residuals
return (self.data * np.sqrt(self.fit_weights)) - model.sum(axis=1).getA1()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment