### Updated code to be compatible with Maoppy version 1.2, which normalizes the...

Updated code to be compatible with Maoppy version 1.2, which normalizes the profile at infinity, hence no renormalization is required anymore. Also, updated documentation.
 ... ... @@ -240,7 +240,7 @@ The transmission curve of the filter must be available as a plain ascii file, co \section{The PSF models} \label{sec:psf} \pampelmuse currently uses purely analytical PSF profiles. The user can choose between a Gaussian and a Moffat profile. An extensive discussion of the profiles and all associated formulae can be found in Appendix~\ref{app:profiles}. Here we restrict the discussion to the main formulae and most important properties. The analytical profile of the Gaussian is \pampelmuse currently uses purely analytical PSF profiles. Most of the available profiles are based on either a Gaussian or a Moffat profile. An extensive discussion of the profiles and all associated formulae can be found in Appendix~\ref{app:profiles}. Here we restrict the discussion to the main formulae and most important properties. The analytical profile of the Gaussian is % \begin{equation} G(r) = \Sigma_\text{0}\exp\left\{\frac{-r^2}{2\sigma^2}\right\}\,, ... ... @@ -254,12 +254,14 @@ whereas the Moffat profile is described via % For either profile, $\Sigma_\mathrm{0}$ is the central intensity. The width of the Gaussian is controlled by the standard deviation $\sigma$, while the effective radius $\rd$ defines the width of the Moffat. However, instead of $\sigma$ and $\beta$, \pampelmuse uses the FWHM as a measure for the profile width (see Appendix~\ref{app:profiles}). Note that the FWHM is always measured in units of spaxels. The Moffat profile has a second shape parameter, $\beta$. It measures the kurtosis of the profile, i.e. the strength of the wings of the PSF. Note that higher values of $\beta$ correspond to less pronounced wings. From my experience, $\beta\sim2.5$ is a typical value for a seeing-limited PSF. For this value, a Moffat profile has much more pronounced wings than a Gaussian profile. Therefore I strongly recommend using the Moffat profile because the Gaussian profile will almost certainly underestimate the wings of the PSF. Since version 0.100, \pampelmuse also supports combinations of the two profiles, namely a double Gaussian and a Moffat-Gaussian profile. The Moffat profile has a second shape parameter, $\beta$. It measures the kurtosis of the profile, i.e. the strength of the wings of the PSF. Note that higher values of $\beta$ correspond to less pronounced wings. From my experience, $\beta\sim2.5$ is a typical value for a seeing-limited PSF. For this value, a Moffat profile has much more pronounced wings than a Gaussian profile. Therefore I strongly recommend using the Moffat profile because the Gaussian profile will almost certainly underestimate the wings of the PSF. Note that \pampelmuse also supports combinations of the two profiles, as listed below. The radius $r$ incorporates the ellipticity $e$ and the position angle $\theta$ of each profile. Thus, besides the actual centroid coordinates, the Moffat profile has $4$ free parameters, the Gaussian profile has $3$. To facilitate the handling of adaptive optics (AO) data, in particular MUSE narrow field mode (NFM) observations, the MAOPPY PSF model developed by \citet{2019A&A...628A..99F} is made available in \pampelmuse. Besides the parameters $e$ and $\theta$ that it shares with the other PSF models, the MAOPPY PSF has five parameters that can be optimized during the analysis. They are the Fried parameter $r_0$, the kurtosis $\beta$ of the Moffat contribution, the effective radius $\alpha$ of the Moffat contribution, the AO-corrected phase power spectral density (PSD) background $bck$, and the residual variance $amp$. Note that some parameters have been renamed compared to the names provided in Table~2 of \citet{2019A&A...628A..99F} in order to avoid confusion with other fitting parameters. In case you use this PSF model, please cite \citet{2019A&A...628A..99F}. The properties of the the PSF during the analysis are controlled via the \emph{psf} section of the configuration. The profile in use can be set using the \emph{profile} keyword. \pmparameter{psf}{profile}{string}{moffat}{The name of analytical PSF profile to be used. Currently supported are 'moffat', 'gauss', 'double\_gauss', and 'moffat\_gauss'.} \pmparameter{psf}{profile}{string}{moffat}{The name of analytical PSF profile to be used. Currently supported are 'moffat', 'gauss', 'double\_gauss', 'moffat\_gauss', 'double\_moffat', and 'maoppy'.} The initial values of the PSF parameters can be specified directly in the \emph{psf}-section of the configuration. \pmparameter{psf}{e}{float}{0.0}{The initial value of the ellipticity of the PSF.} ... ...
 \relax \providecommand\hyper@newdestlabel{} \@writefile{toc}{\contentsline {section}{\numberline {B}The PSF profiles}{23}{appendix.B}} \@writefile{toc}{\contentsline {section}{\numberline {B}The PSF profiles}{23}{appendix.B}\protected@file@percent } \newlabel{app:profiles}{{B}{23}{The PSF profiles}{appendix.B}{}} \@writefile{toc}{\contentsline {subsection}{\numberline {B.1}Handling of positions}{23}{subsection.B.1}} \@writefile{toc}{\contentsline {subsection}{\numberline {B.1}Handling of positions}{23}{subsection.B.1}\protected@file@percent } \newlabel{sec:positions}{{B.1}{23}{Handling of positions}{subsection.B.1}{}} \newlabel{eq:psf_x}{{7}{23}{Handling of positions}{equation.B.7}{}} \newlabel{eq:psf_y}{{8}{23}{Handling of positions}{equation.B.8}{}} \newlabel{eq:psf_radius}{{9}{23}{Handling of positions}{equation.B.9}{}} \newlabel{eq:coordtransX}{{14}{23}{Handling of positions}{equation.B.14}{}} \newlabel{eq:coordtransY}{{15}{23}{Handling of positions}{equation.B.15}{}} \@writefile{toc}{\contentsline {subsection}{\numberline {B.2}The Gaussian profile}{24}{subsection.B.2}} \@writefile{toc}{\contentsline {subsection}{\numberline {B.2}The Gaussian profile}{24}{subsection.B.2}\protected@file@percent } \newlabel{sec:gauss}{{B.2}{24}{The Gaussian profile}{subsection.B.2}{}} \newlabel{eq:intensity}{{22}{24}{The Gaussian profile}{equation.B.22}{}} \@writefile{toc}{\contentsline {subsection}{\numberline {B.3}The Moffat profile}{25}{subsection.B.3}} \@writefile{toc}{\contentsline {subsection}{\numberline {B.3}The Moffat profile}{25}{subsection.B.3}\protected@file@percent } \newlabel{sec:moffat}{{B.3}{25}{The Moffat profile}{subsection.B.3}{}} \newlabel{eq:moffat}{{27}{25}{The Moffat profile}{equation.B.27}{}} \newlabel{eq:moffat_fwhm}{{28}{25}{The Moffat profile}{equation.B.28}{}} \@writefile{lof}{\contentsline {figure}{\numberline {6}{\ignorespaces Inverse ratio between the $\text {FWHM}$ of a Moffat profile and the $\text {FWHM}$ of its best Gaussian fit, as a function of the $\beta$ value of the Moffat profile and its S/N.}}{25}{figure.6}} \@writefile{lof}{\contentsline {figure}{\numberline {6}{\ignorespaces Inverse ratio between the $\text {FWHM}$ of a Moffat profile and the $\text {FWHM}$ of its best Gaussian fit, as a function of the $\beta$ value of the Moffat profile and its S/N.}}{25}{figure.6}\protected@file@percent } \newlabel{fig:fwhm_ratio_gauss_moffat}{{6}{25}{Inverse ratio between the $\fwhm$ of a Moffat profile and the $\fwhm$ of its best Gaussian fit, as a function of the $\beta$ value of the Moffat profile and its S/N}{figure.6}{}}
 ... ... @@ -862,7 +862,8 @@ class FitLayer(object): elif sources is None: sources = self.catalog.index psf_radii, psf_flux_scaling = self.scale_psf_radii() # scale PSF radii with fluxes of sources psf_radii = self.scale_psf_radii() # loop over sources & extract information from PSF instances for index in sources: ... ... @@ -873,10 +874,8 @@ class FitLayer(object): psf.uc = self.catalog.at[index, 'x'] - self.offsets.at[index, 'x'] psf.vc = self.catalog.at[index, 'y'] - self.offsets.at[index, 'y'] # apply flux-dependent PSF radius (if requested); scale PSF flux if necessary # apply flux-dependent PSF radius (if requested) psf.maxrad = psf_radii[index] if psf_flux_scaling is not None: psf.mag = self.catalog.at[index, "magnitude"] - 2.5*np.log10(psf_flux_scaling[index]) current_pixels, current_fluxes = psf.get_flux(apply_lut=True) if current_pixels.size == 0: ... ... @@ -963,22 +962,33 @@ class FitLayer(object): return _psf_matrix.tocsc(), source_indices def scale_psf_radii(self): """ This method scales the radii out to which the individual PSF profiles are defined according to the current fluxes of the sources. Notes ----- The idea is that the PSF radius is halved every time for every decrease of the flux by a factor :math:\\xi. The radius scaling is capped for sources fainter then :math:\\xi^3 relative to the brightest source, which means that the maximum scaling is :math:2^3. The scaling is performed according to the equation: .. math:: r_i = r_\mathtt{PSF} 2^{(\log \frac{f_i}{max(f)} / \log \frac{1}{\\xi})} for :math:0<\\xi<1. Returns ------- psf_radii : pandas.Series The scaled PSF radii for all available sources. """ # if no fluxes available, all sources get maximum radius if not self.psf_radius_scaling or self.fluxes.max() <= 0: return pd.Series(self.psf_radius, index=self.catalog.index), None # MAOPPY PSFs are automatically normalized to a flux of 1 inside the defined grid. Hence, the total PSF flux # (when integrating to INF) will increase the smaller the PSF radius is chosen. To account for that, the curve # of growth is determined for the PSF of the reference source and used to apply a correction factor to each PSF if self.psf_radius_scaling and self.psf_class == Maoppy: # make sure reference source is defined out to maximum radius self.psf_instances[self.reference_index].maxrad = self.psf_radius # create interpolator on curve of growth of reference source r, f = self.psf_instances[self.reference_index].curve_of_growth ip = interpolate.interp1d(x=r, y=f/f[-1]) else: ip = None return pd.Series(self.psf_radius, index=self.catalog.index) # scale PSF radius based on flux - the minimum PSF radius is set to (0.5^3) x PSF radius flux_ratios = self.fluxes/self.fluxes.max() ... ... @@ -988,7 +998,7 @@ class FitLayer(object): psf_radii = radii[self.catalog['component'].values] psf_radii.index = self.catalog.index return psf_radii, pd.Series(ip(psf_radii), index=psf_radii.index) if ip is not None else None return psf_radii def fit_sources_consecutively(self, sources, is_psf_source=None, psf_parameters=None, fit_positions=True, local_background=False, return_psf_profiles=True): ... ... @@ -1511,7 +1521,7 @@ class FitLayer(object): # 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) valid = (results[name] >= vmin) & (results[name] <= vmax) if np.any(valid): _final[name] = np.average(results.loc[valid, name], weights=results.loc[valid, 'snr']) else: ... ...
 ... ... @@ -17,16 +17,33 @@ def resample_to_2d(x, y, z, shape=None): Given a set of values z defined on arbitrary coordinates (x, y), the method will translate the data to a regularly sampled 2d array. The idea is that actual resampling is only performed if the input x and y coordinates do not correspond to grid points of the final 2d array. Otherwise, the available pixels values are just inserted into the output array. Parameters ---------- x y z shape x : sequence (1D) of floats The x-coordinates on which the input data are defined. Note that following numpy conventions, the x-axis will be the second axis of the output array. y : sequence (1D) of floats The y-coordinates on which the input data are defined. Note that following numpy conventions, the y-axis will be the first axis of the output array. z : sequence (1D) of floats The values of the input data. Parameters x, y, and z must have the same length. shape : tuple of length 2, optional The 2D shape of the requested output array. If not defined, the data will be resampled to the smallest 2d grid that encompasses all input data. Returns ------- out : array (2D) of floats The array of resampled values. """ if shape is None: shape = (int(y.max()) + 1, int(x.max()) + 1) ... ... @@ -51,15 +68,25 @@ def resample_to_1d(z, x, y,): Given a regularly sampled 2D array of values z, the method will convert the data to the provided (randomly sampled) coordinates x and y. The idea is that actual resampling is only performed when the requested output coordinates are not a subset of the grid points on which the input data are defined. Otherwise, the output values are obtained by basic slicing. Parameters ---------- z x y z : array (2D) of floats The input grid of values. x : array (1D) of floats The x-coordinates of the output grid. y : array (1D) of floats The y-coordinates of the output grid. Parameters x and y must have the same length. Returns ------- znew : array (1D) of floats The interpolated values at the coordinates provided by x and y. """ is_int = np.all(x == x.astype(np.int32)) & np.all(y == y.astype(np.int32)) ... ... @@ -77,9 +104,69 @@ def resample_to_1d(z, x, y,): class Maoppy(Profile): """ This class implements the MAOPPY PSF profile described by Fetick et al. (2019). It is optimized for the modelling of PSF profiles in adaptive optics (AO) observations, e.g. performed with the MUSE narrow field mode. def __init__(self, uc, vc, mag, r0=0.2, bck=0.005, amp=1, alpha=0.05, beta=2.0, wvl=5e-7, **kwargs): References ---------- 1) Fetick et al. (2019, A&A, 628, A99) https://ui.adsabs.harvard.edu/abs/2019A&A...628A..99F/abstract """ def __init__(self, uc, vc, mag, r0=0.2, bck=0.005, amp=1, alpha=0.05, beta=2.0, wvl=5e-7, **kwargs): """ Initialize a new instance of the MAOPPY PSF. Notes ----- There are some differences in the parameter set of the MAOPPY PSF that is used by PampelMuse compared to the paper by Fetick et al. (2019).: 1) The parameters C (the AO-corrected phase PSD background) and A (the residual variance) have been renamed to bck and amp. 2) Instead of defining two values for the effective radius of the Moffat, alpha_x and alpha_y, a single value alpha is used in combination with an ellipticity e. Parameters ---------- uc : float The central coordinate of the PSF profile along the 1st reference axis. vc : float The central coordinate of the PSF profile along the 2nd reference axis. mag : float The integrated magnitude of the PSF profile, defined such that the integrated flux is unity for mag=0. r0 : float, optional The Fried parameter, given in units of [m]. bck : float, optional The AO-corrected phase PSD background, given in units of [m^2 rad^2]. amp : float, optional The residual variance, given in units of [rad^2]. alpha : float, optional The effective radius of the Moffat profile, given in units of [m^-1]. beta : float, optional The kurtosis parameter of the Moffat profile. wvl : float, optional The wavelength of the observation, given in units of [m]. kwargs Any additional keyword arguments are passed on to the init method of the parent class. Returns ------- The newly defined instance. """ # oversampling of the PSF on the PampelMuse level is not supported, because MAOPPY performs its # own oversampling when computing a PSF profile. if 'osradius' in kwargs and kwargs['osradius'] > 0: logger.error('Oversampling not supported for MAOPPY PSF. Setting osradius=0 (was {})'.format( kwargs['osradius'])) ... ... @@ -89,8 +176,10 @@ class Maoppy(Profile): kwargs['osfactor'])) kwargs['osfactor'] = 1 # call initializer of parent class. super(Maoppy, self).__init__(uc=uc, vc=vc, mag=mag, **kwargs) # initial parameter definitions self._r0 = r0 self._bck = bck self._amp = amp ... ... @@ -100,98 +189,241 @@ class Maoppy(Profile): @property def r0(self): """ Returns the current Fried parameter of the PSF profile. """ return self._r0 @r0.setter def r0(self, value): """ Set the Fried parameter of the PSF profile. Parameters ---------- value : float The Fried parameter, given in units of [m]. """ self._r0 = value @property def bck(self): """ Returns the current AO-corrected phase PSD background (labeled C in Fetick et al. 2019) of the PSF profile. """ return self._bck @bck.setter def bck(self, value): """ Set the AO-corrected phase PSD background (labeled C in Fetick et al. 2019) of the PSF profile. Parameters ---------- value : float The AO-corrected phase PSD background, given in units of [m^2 rad^2]. """ self._bck = value @property def amp(self): """ Returns the current residual variance parameter (labeled A in Fetick et al. 2019) of the PSF profile. """ return self._amp @amp.setter def amp(self, value): """ Set the residual variance parameter (labeled A in Fetick et al. 2019) of the PSF profile. Parameters ---------- value : float The residual variance parameter, given in units of [rad^2]. """ self._amp = value @property def alpha(self): """ Returns the current Moffat effective radius of the PSF profile. """ return self._alpha @alpha.setter def alpha(self, value): """ Set the Moffat effective radius alpha of the PSF profile. Parameters ---------- value : float The Moffat effective radius, given in units of [m^-1]. """ self._alpha = value @property def beta(self): """ Returns the current Moffat beta of the PSF profile. """ return self._beta @beta.setter def beta(self, value): """ Set the Moffat beta for the PSF profile. Parameters ---------- value : float The Moffat beta. """ self._beta = value @property def wvl(self): """ Returns the current wavelength of the PSF profile. """ return self._wvl @wvl.setter def wvl(self, value): """ Set the wavelength for the PSF profile. Parameters ---------- value : float The wavelength, given in units of [m]. """ self._wvl = value @property def q(self): """ Returns the axis ratio corresponding to the current ellipticity. """ return 1. - self.e @property def samp(self): """ Returns the sampling of the PSF profile at the current wavelength. """ return MUSE_NFM.samp(self.wvl) @property def f(self): """ Returns the pixel values of the MAOPPY PSF for all pixels defined in the instance and for the current set of PSF parameters. """ # check if flux needs to be recalculated if self._f is None: if not self.used.any(): self._f = np.array([]) else: # init MAOPPY code psf = Psfao(Npix=self.maoppy_shape, system=MUSE_NFM, samp=self.samp) # get offsets of centre to nearest lower pixels. dx = self.xc - int(self.xc) dy = self.yc - int(self.yc) # run MAOPPY code __f = psf((self.r0, self.bck, self.amp, self.alpha, self.q, self.theta, self.beta), dx=dy, dy=dx) # account for offset between assumed and MAOPPY-internal centre x_center, y_center = self.get_maoppy_centre(k=psf._k) shift_x = int(self.xc) - x_center shift_y = int(self.yc) - y_center # MAOPPY returns 2D array. Need to revert to pixel sampling used by instance. ip = RectBivariateSpline(x=np.arange(self.maoppy_shape) + shift_x, y=np.arange(self.maoppy_shape) + shift_y, z=__f, kx=1, ky=1) # scale integrated flux to provided magnitude of source self._f = self.int_flux*ip(self.x[self.used], self.y[self.used], grid=False) return self._f @property def maoppy_shape(self): """ Returns the minimum shape of a 2D array with an even number of pixels along both axes that fits the currently defined PSF profile. """ return 2*int(round(self.maxrad)), 2*int(round(self.maxrad)) def get_maoppy_centre(self, k, shape=None): """ Determine the centre of the PSF profile that is used by the MAOPPY code for a given oversampling factor k and array size shape. Parameters ---------- k : int The oversampling used by MAOPPY. shape : tuple of 2 integers, optional The size of the array on which the PSF is computed. By default, the array size of the currently defined PSF is used. Returns ------- centre : tuple of 2 floats The central coordinates of the PSF along the y- and x-axes. """ if shape is None: shape = self.maoppy_shape return shape/2 - (k - 1.)/(2.*k), shape/2 - (k - 1.)/(2.*k) @classmethod def calculate_fwhm(cls, r0, bck, amp, alpha, e, beta, wvl=6e-7, **kwargs): """ Calculate and return the FWHM of the MAOPPY PSF profile for a given set of parameters. Parameters ---------- r0 : float The Fried parameter of the PSF profile, given in units of [m]. bck : float The AO-corrected phase PSD (power spectral density) background, given in units of [rad^2 m^2]. Note that this parameter is called C in the publication by Fetick et al. (2019). amp : float The residual variance of the PSD, given in units of [rad^2]. Note that this parameter is called A in the publication of Fetick et al. (2019). alpha : float The effective radius of the Moffat profile, given in units of [m^-1]. e : float The ellipticity of the Moffat profile. Must be >= 0 and < 1. beta : float The kurtosis paramter of the Moffat profile. wvl : float, optional The wavelength for which the PSF should be calculated, given in units of [m]. kwargs Returns ------- fwhm : float The FWHM of the MAOPPY PSF for the given input parameters. """ # make sure no additional keywords were provided. if kwargs: raise IOError('Unknown argument(s) to call of Maoppy.fwhm: {0}'.format(kwargs)) ... ... @@ -205,10 +437,52 @@ class Maoppy(Profile): @property def fwhm(self): """ Returns ------- fwhm : float The FWHM of the MAOPPY PSF profile with the given parameters. """ return np.sqrt(np.sum(self.f > (0.5*np.max(self.f)))/np.pi) def fitter(self, psf, weights, parameters, fit_positions=True, fit_background=False, psf_padding=0, **kwargs): """ This methods implements the fit of a MAOPPY PSF to a given data array. Notes ----- The actual fit is performed using the method psffit available in the MAOPPY package. It requires a 2D array with an even number of pixels along each dimension. Therefore, the provided 1D data arrays need to be resampled prior to the fit.