Skip to content

Commit

Permalink
Merge pull request vortex-exoplanet#654 from IainHammond/master
Browse files Browse the repository at this point in the history
more options to get_values_optimize, get_mu_and_sigma, speckle_noise_uncertainty. typo & docstring fixes
  • Loading branch information
VChristiaens authored Jan 6, 2025
2 parents e6cafa9 + acee4ab commit 555931c
Show file tree
Hide file tree
Showing 11 changed files with 43 additions and 25 deletions.
18 changes: 16 additions & 2 deletions vip_hci/fm/negfc_fmerit.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
__all__ = ["get_mu_and_sigma"]

import numpy as np

from hciplot import plot_frames
from skimage.draw import disk
from ..fm import cube_inject_companions, cube_planet_free
Expand Down Expand Up @@ -490,8 +491,10 @@ def get_values_optimize(
max_frames_lib = algo_opt_copy.pop("max_frames_lib", 200)
radius_int = max(1, int(np.floor(r_guess - annulus_width / 2)))
radius_int = algo_opt_copy.pop("radius_int", radius_int)
asize = algo_opt_copy.pop("asize", annulus_width)
delta_rot = algo_opt_copy.pop("delta_rot", delta_rot)
# crop cube to just be larger than annulus => FASTER PCA
crop_sz = int(2 * np.ceil(radius_int + annulus_width + 1))
crop_sz = int(2 * np.ceil(radius_int + asize + 1))
if not crop_sz % 2:
crop_sz += 1
if crop_sz < cube.shape[-2] and crop_sz < cube.shape[-1]:
Expand All @@ -507,7 +510,7 @@ def get_values_optimize(
cube_ref=cube_ref,
radius_int=radius_int,
fwhm=fwhm,
asize=annulus_width,
asize=asize,
delta_rot=delta_rot,
ncomp=ncomp,
svd_mode=svd_mode,
Expand Down Expand Up @@ -553,6 +556,9 @@ def get_values_optimize(
elif algo == pca:
scale_list = algo_opt_copy.pop("scale_list", None)
ifs_collapse_range = algo_opt_copy.pop("ifs_collapse_range", "all")
mask_rdi = algo_options.pop("mask_rdi", None)
delta_rot = algo_options.pop("delta_rot", delta_rot)
source_xy = algo_options.pop("source_xy", None)
res = pca(
cube=cube,
angle_list=angs,
Expand All @@ -561,6 +567,8 @@ def get_values_optimize(
ncomp=ncomp,
svd_mode=svd_mode,
scaling=scaling,
delta_rot=delta_rot,
source_xy=source_xy,
imlib=imlib,
interpolation=interpolation,
collapse=collapse,
Expand Down Expand Up @@ -748,6 +756,7 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm,
collapse = algo_options.get("collapse", collapse)

radius_int = max(int(np.floor(r_guess - annulus_width / 2)), 0)
radius_int = algo_options.get("radius_int", radius_int)

# not recommended, except if large-scale residual sky present (NIRC2-L')
hp_filter = algo_options.get("hp_filter", None)
Expand Down Expand Up @@ -864,6 +873,7 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm,
scale_list = algo_options.get("scale_list", None)
ifs_collapse_range = algo_options.get("ifs_collapse_range", "all")
nproc = algo_options.get("nproc", 1)
source_xy = algo_options.get("source_xy", None)

pca_res = pca(
cube=array,
Expand All @@ -873,6 +883,8 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm,
ncomp=ncomp,
svd_mode=svd_mode,
scaling=scaling,
delta_rot=delta_rot,
source_xy=source_xy,
imlib=imlib,
interpolation=interpolation,
collapse=collapse,
Expand All @@ -890,6 +902,8 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm,
ncomp=ncomp,
svd_mode=svd_mode,
scaling=scaling,
delta_rot=delta_rot,
source_xy=source_xy,
imlib=imlib,
interpolation=interpolation,
collapse=collapse,
Expand Down
6 changes: 3 additions & 3 deletions vip_hci/fm/negfc_mcmc.py
Original file line number Diff line number Diff line change
Expand Up @@ -516,7 +516,7 @@ def mcmc_negfc_sampling(cube, angs, psfn, initial_state, algo=pca_annulus,
'sum' when residual speckle noise is still significant):\
:math:`\chi^2 = N \sigma_{I_j}(values,ddof=1)*values.size`
- ``mu_sigma=True`` or a tuple (as in [CHR21]_, new default):\
:math:`\chi^2 = \sum\frac{(I_j- mu)^2}{\sigma^2}`
:math:`\chi^2 = \sum\frac{I_j- \mu)^2}{\sigma^2}`
where :math:`j \in {1,...,N}` with N the total number of pixels
contained in the circular aperture, :math:`\sigma_{I_j}` is the standard
Expand Down Expand Up @@ -999,7 +999,7 @@ def mcmc_negfc_sampling(cube, angs, psfn, initial_state, algo=pca_annulus,
if verbosity > 0:
print(' r_hat = {}'.format(rhat))
cond = rhat <= rhat_threshold
print(' r_hat <= threshold = {} \n'.format(cond))
print(' r_hat <= threshold = {} \n'.format(cond), flush=True)
# We test the rhat.
if (rhat <= rhat_threshold).all():
rhat_count += 1
Expand All @@ -1025,7 +1025,7 @@ def mcmc_negfc_sampling(cube, angs, psfn, initial_state, algo=pca_annulus,
thr = 1./ac_c
if verbosity > 0:
print('Auto-corr tau/N = {}'.format(rhat))
print('tau/N <= {} = {} \n'.format(thr, rhat < thr))
print('tau/N <= {} = {} \n'.format(thr, rhat < thr), flush=True)
if (rhat <= thr).all():
ac_count += 1
if verbosity > 0:
Expand Down
4 changes: 2 additions & 2 deletions vip_hci/fm/negfc_simplex.py
Original file line number Diff line number Diff line change
Expand Up @@ -619,7 +619,7 @@ def firstguess(cube, angs, psfn, planets_xy_coord, ncomp=1, fwhm=4,
If set to None: not used, and falls back to original version of the
algorithm, using fmerit.
If a tuple of 2 elements: should be the mean and standard deviation of
pixel intensities in an annulus centered on the lcoation of the
pixel intensities in an annulus centered on the location of the
companion candidate, excluding the area directly adjacent to the CC.
If set to anything else, but None/False/tuple: will compute said mean
and standard deviation automatically.
Expand Down Expand Up @@ -762,7 +762,7 @@ def firstguess(cube, angs, psfn, planets_xy_coord, ncomp=1, fwhm=4,
print(msg3a.format(i_planet, r_pre, theta_pre))
msg3b = 'Planet {}: preliminary flux guess: '.format(i_planet)
for z in range(len(f_pre)):
msg3b += '{:.1f}'.format(f_pre[z])
msg3b += '{:.2f}'.format(f_pre[z])
if z < len(f_pre)-1:
msg3b += ', '
print(msg3b)
Expand Down
15 changes: 8 additions & 7 deletions vip_hci/fm/negfc_speckle_noise.py
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ def speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo,
of spectral channels (flux at each wavelength).
angle_range: 1d numpy array
Range of angles (counted from x=0 axis, counter-clockwise) at which the
fake companions will be injected, in [0,360[.
fake companions will be injected, in [0,360].
derot_angles: 1d numpy array
Derotation angles for ADI. Length should match input cube.
algo: python routine
Expand Down Expand Up @@ -96,7 +96,7 @@ def speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo,
If set to None: not used, and falls back to original version of the
algorithm, using fmerit.
If a tuple of 2 elements: should be the mean and standard deviation of
pixel intensities in an annulus centered on the lcoation of the
pixel intensities in an annulus centered on the location of the
companion candidate, excluding the area directly adjacent to the CC.
If set to anything else, but None/False/tuple: will compute said mean
and standard deviation automatically.
Expand Down Expand Up @@ -157,7 +157,7 @@ def speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo,
save: bool, optional
If True, the result are pickled.
verbose: bool, optional
If True, informations are displayed in the shell.
If True, information is displayed in the shell.
plot: bool, optional
Whether to plot the gaussian fit to the distributions of parameter
deviations (between retrieved and injected).
Expand Down Expand Up @@ -366,9 +366,6 @@ def _estimate_speckle_one_angle(angle, cube_pf, psfn, angs, r_true, f_true,
imlib=imlib, interpolation=interpolation,
verbose=False)

ncomp = algo_options.get('ncomp', 1)
annulus_width = algo_options.get('annulus_width', int(fwhm))

if cube_pf.ndim == 4:
p_ini = [r_true, angle]
for f in f_true:
Expand All @@ -377,10 +374,14 @@ def _estimate_speckle_one_angle(angle, cube_pf, psfn, angs, r_true, f_true,
else:
p_ini = (r_true, angle, f_true)

ncomp = algo_options.get('ncomp', 1)
annulus_width = algo_options.get('annulus_width', int(fwhm))
delta_rot = algo_options.get('delta_rot', 1)

res_simplex = firstguess_simplex(p_ini, cube_fc, angs,
psfn, ncomp, fwhm, annulus_width,
aperture_radius, cube_ref=cube_ref,
fmerit=fmerit, algo=algo,
fmerit=fmerit, algo=algo, delta_rot=delta_rot,
algo_options=algo_options, imlib=imlib,
interpolation=interpolation,
transmission=transmission,
Expand Down
4 changes: 2 additions & 2 deletions vip_hci/greedy/ipca_fullfr.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,6 @@
from dataclasses import dataclass
import numpy as np
from typing import Union, List
import warnings
from ..config.paramenum import ALGO_KEY
from ..config.utils_param import separate_kwargs_dict
from ..config import Progressbar
Expand All @@ -62,8 +61,9 @@
from GreeDS import GreeDS
no_greeds = False
except ImportError:
from warnings import warn
msg = "GreeDS python bindings are missing."
warnings.warn(msg, ImportWarning)
warn(msg, ImportWarning)
no_greeds = True


Expand Down
2 changes: 1 addition & 1 deletion vip_hci/preproc/derotation.py
Original file line number Diff line number Diff line change
Expand Up @@ -291,7 +291,7 @@ def frame_rotate(array, angle, imlib='vip-fft', interpolation='lanczos4',
elif interpolation == 'lanczos4':
intp = cv2.INTER_LANCZOS4
else:
raise ValueError('Opencv interpolation method not recognized')
raise ValueError(f'Opencv interpolation method `{interpolation}` is not recognized')

if border_mode == 'constant':
bormo = cv2.BORDER_CONSTANT # iiiiii|abcdefgh|iiiiiii
Expand Down
6 changes: 3 additions & 3 deletions vip_hci/psfsub/pca_fullfr.py
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ class PCA_Params:


def pca(*all_args: List, **all_kwargs: dict):
"""Full-frame PCA algorithm applied to PSF substraction.
"""Full-frame PCA algorithm applied to PSF subtraction.
The reference PSF and the quasi-static speckle pattern are modeled using
Principal Component Analysis. Depending on the input parameters this PCA
Expand All @@ -156,11 +156,11 @@ def pca(*all_args: List, **all_kwargs: dict):
Parameters
----------
all_args: list, optional
Positionnal arguments for the PCA algorithm. Full list of parameters
Positional arguments for the PCA algorithm. Full list of parameters
below.
all_kwargs: dictionary, optional
Mix of keyword arguments that can initialize a PCA_Params and the
optional 'rot_options' dictionnary (with keyword values ``border_mode``,
optional 'rot_options' dictionary (with keyword values ``border_mode``,
``mask_val``, ``edge_blend``, ``interp_zeros``, ``ker``; see docstring
of ``vip_hci.preproc.frame_rotate``). Can also contain a PCA_Params
dictionary named `algo_params`.
Expand Down
2 changes: 1 addition & 1 deletion vip_hci/psfsub/pca_local.py
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ def pca_annular(*all_args: List, **all_kwargs: dict):
parameters below.
all_kwargs: dictionary, optional
Mix of keyword arguments that can initialize a PCA_ANNULAR_Params and
the optional 'rot_options' dictionnary, with keyword values for
the optional 'rot_options' dictionary, with keyword values for
"border_mode", "mask_val", "edge_blend", "interp_zeros", "ker" (see
documentation of ``vip_hci.preproc.frame_rotate``). Can also contain a
PCA_ANNULAR_Params named as `algo_params`.
Expand Down
3 changes: 3 additions & 0 deletions vip_hci/psfsub/utils_pca.py
Original file line number Diff line number Diff line change
Expand Up @@ -309,6 +309,9 @@ def get_snr(frame, y, x, fwhm, fmerit):
'to (PC_INI, PC_MAX) or (PC_INI, PC_MAX, STEP)')
pclist = list(range(pcmin, pcmax+1, step))

if fmerit not in ['px', 'max', 'mean']:
raise ValueError(f"Invalid value for fmerit: {fmerit}.")

# Getting `pcmax` principal components once
if mode == 'fullfr':
matrix = prepare_matrix(cube, scaling, mask_center_px, verbose=False)
Expand Down
2 changes: 1 addition & 1 deletion vip_hci/var/filters.py
Original file line number Diff line number Diff line change
Expand Up @@ -574,7 +574,7 @@ def cube_filter_lowpass(array, mode='gauss', median_size=5, fwhm_size=5,
mode is set to 'psf'.
mask: numpy ndarray, optional
Binary mask indicating where the low-pass filtered image should be
interpolated with astropy.convolution. This otion can be useful if the
interpolated with astropy.convolution. This option can be useful if the
low-pass filtered image is aimed to capture low-spatial frequency sky
signal, while avoiding a stellar halo (set to one in the binary mask).
Note: only works with Gaussian kernel or PSF convolution.
Expand Down
6 changes: 3 additions & 3 deletions vip_hci/var/fit_2d.py
Original file line number Diff line number Diff line change
Expand Up @@ -75,9 +75,9 @@ def create_synth_psf(model='gauss', shape=(9, 9), amplitude=1, x_mean=None,
Note
----
More information can be found at the following links:
http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Gaussian2D.html
http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Moffat2D.html
http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.AiryDisk2D.html
https://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Gaussian2D.html
https://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Moffat2D.html
https://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.AiryDisk2D.html
https://www.gnu.org/software/gnuastro/manual/html_node/PSF.html
web.ipac.caltech.edu/staff/fmasci/home/astro_refs/PSFtheory.pdf
Expand Down

0 comments on commit 555931c

Please sign in to comment.