Commit dbcbf3d2 authored by skamann's avatar skamann
Browse files

Revised MAOPPY implementation to be consistent with new calling sequence of psffit().

Revised averaging of PSF parameters in CUBEFIT to avoid warnings related to usage of NaNs in sigma_clip of astropy.
parent 75c0b804
......@@ -57,7 +57,7 @@ parameters are not fitted.
Latest Git revision
-------------------
2022/06/14
2022/06/27
"""
import collections
import contextlib
......@@ -74,14 +74,14 @@ if scipy.__version__ < '1.5.0':
from scipy.stats import median_absolute_deviation as mad
else:
from scipy.stats import median_abs_deviation as mad
from astropy.stats import sigma_clip
from astropy.stats import SigmaClip
from ..psf.profiles import *
from ..core.coordinates import Transformation
__author__ = "Sebastian Kamann (s.kamann@ljmu.ac.uk)"
__revision__ = 20220614
__revision__ = 20220627
logger = logging.getLogger(__name__)
......@@ -1533,13 +1533,14 @@ class FitLayer(object):
# update PSF parameters, using S/N as weights:
if psf_parameters is not None:
_final = pd.Series(0., index=psf_parameters)
sc = SigmaClip(sigma=4., stdfunc='mad_std', cenfunc='median')
for name in psf_parameters:
# exclude results deviating by more than 4x the median absolute deviation from the median value
_, vmin, vmax = sigma_clip(results[name], sigma=4, return_bounds=True,
stdfunc=lambda x, axis: mad(x, nan_policy='omit', axis=axis))
valid = (results[name] >= vmin) & (results[name] <= vmax)
if np.any(valid):
_final[name] = np.average(results.loc[valid, name], weights=results.loc[valid, 'snr'])
clipped = sc(np.ma.MaskedArray(results[name].values, mask=np.isnan(results[name].values)),
return_bounds=False)
if not clipped.mask.all():
_final[name] = np.average(results.loc[~clipped.mask, name],
weights=results.loc[~clipped.mask, 'snr'])
else:
logger.warning('Cannot calculate mean of "{0}" because no PSF fit succeeded.'.format(name))
_final[name] = np.nan
......@@ -1650,7 +1651,7 @@ class FitLayer(object):
data.append(value[1])
combined = pd.DataFrame(np.stack(data, axis=1), columns=source_ids,
index=pd.Float64Index(np.nanmean(radii, axis=0)))
index=pd.Index(np.nanmean(radii, axis=0), dtype=np.float64))
return combined
def update_look_up_table(self, residuals, centroids, weights=None):
......
......@@ -28,7 +28,7 @@ class pampelmuse.core.psf.profiles.Maoppy
Latest Git revision
-------------------
2022/03/30
2022/06/27
"""
import importlib.util
import io
......@@ -49,7 +49,7 @@ else:
__author__ = "Sebastian Kamann (s.kamann@ljmu.ac.uk)"
__revision__ = 20220330
__revision__ = 20220627
logger = logging.getLogger(__name__)
......@@ -614,8 +614,8 @@ class Maoppy(Profile):
dx = round(x_center - self.xc) - (x_center - self.xc)
dy = round(y_center - self.yc) - (y_center - self.yc)
npix = max(psf2d.shape)
if psf_padding > 0:
npix = max(psf2d.shape)
npixfit = npix + 2*psf_padding
else:
npixfit = None
......@@ -624,9 +624,10 @@ class Maoppy(Profile):
# code prints out some "-" symbols to indicate progression, which are captured so that is does not interfere
# with the logging
f = io.StringIO()
pmodel = Psfao(npix=(npix+2*psf_padding, npix+2*psf_padding), system=muse_nfm, samp=self.samp)
with redirect_stdout(f):
fitresult = psffit(psf2d, Psfao, x0, weights=weights2d, system=muse_nfm, samp=self.samp, dxdy=(dx, dy),
fixed=fixed, flux_bck=(True, fit_background), positive_bck=False, npixfit=npixfit)
fitresult = psffit(psf2d, pmodel, x0, weights=weights2d, dxdy=(dx, dy), fixed=fixed,
flux_bck=(True, fit_background), positive_bck=False, npixfit=npixfit)
logger.debug(f.getvalue())
# return ellipticity e instead of axis ratio q
......
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