vip_hci.fm package¶

Submodules¶

vip_hci.fm.fakecomp module¶

Module with fake companion injection functions.

vip_hci.fm.fakecomp.collapse_psf_cube(array, size, fwhm=4, verbose=True, collapse='mean')[source]¶

Creates a 2d PSF template from a cube of non-saturated off-axis frames of the star by taking the mean and normalizing the PSF flux.

Parameters:	array (numpy ndarray, 3d) – Input cube. size (int) – Size of the squared subimage. fwhm (float, optional) – The size of the Full Width Half Maximum in pixel. verbose ({True,False}, bool optional) – Whether to print to stdout information about file opening, cropping and completion of the psf template. collapse ({'mean','median'}, string optional) – Defines the way the frames are collapsed.
Returns:	psf_normd – Normalized PSF.
Return type:	numpy ndarray

vip_hci.fm.fakecomp.normalize_psf(array, fwhm='fit', size=None, threshold=None, mask_core=None, model='gauss', imlib='vip-fft', interpolation='lanczos4', force_odd=True, correct_outliers=True, full_output=False, verbose=True, debug=False)[source]¶

Normalizes a PSF (2d or 3d array), to have the flux in a 1xFWHM aperture equal to one. It also allows to crop the array and center the PSF at the center of the array(s).

Parameters:

array (numpy ndarray) – The PSF, 2d (ADI data) or 3d array (IFS data).
fwhm (int, float, 1d array or str, optional) – The Full Width Half Maximum in pixels. It can handle a different FWHM value for different wavelengths (IFS data). If set to ‘fit’ then a model (assuming the PSF is centered in the array) is fitted to estimate the FWHM in 2D or 3D PSF arrays.
size (int or None, optional) – If int it will correspond to the size of the centered sub-image to be cropped form the PSF array. The PSF is assumed to be roughly centered wrt the array.
threshold (None or float, optional) – Sets to zero values smaller than threshold (in the normalized image). This can be used to only leave the core of the PSF.
mask_core (None or float, optional) – Sets the radius of a circular aperture for the core of the PSF, everything else will be set to zero.
model ({'gauss', 'moff', 'airy'}, str optional) – The assumed model used to fit the PSF: either a Gaussian, a Moffat or an Airy 2d distribution.
imlib (str, optional) – See the documentation of the vip_hci.preproc.frame_shift function.
interpolation (str, optional) – See the documentation of the vip_hci.preproc.frame_shift function.
force_odd (bool, optional) – If True the resulting array will have odd size (and the PSF will be placed at its center). If False, and the frame size is even, then the PSF will be put at the center of an even-sized frame.
correct_outliers (bool, optional) – For an input 3D cube (IFS) of PSFs, if the 2D fit fails for one of the channels, whether to interpolate FWHM value from surrounding channels, and recalculate flux and normalization.
full_output (bool, optional) – If True the flux in a FWHM aperture is returned along with the normalized PSF.
verbose (bool, optional) – If True intermediate results are printed out.
debug (bool, optional) – If True the fitting will output additional information and a diagnostic plot will be shown (this might cause a long output if array is 3d and has many slices).

Returns:

psf_norm (numpy ndarray) – The normalized PSF (2d or 3d array).
fwhm_flux (numpy ndarray) – [full_output=True] The flux in a FWHM aperture (it can be a single value or a vector).
fwhm (numpy ndarray) – [full_output=True] The FWHM size. If fwhm is set to ‘fit’ then it is the fitted FWHM value according to the assumed model (the mean in X and Y is returned when model is set to ‘gauss’).

vip_hci.fm.fakecomp.cube_inject_companions(array, psf_template, angle_list, flevel, rad_dists, plsc=None, n_branches=1, theta=0, imlib='vip-fft', interpolation='lanczos4', transmission=None, radial_gradient=False, full_output=False, verbose=False, nproc=1)[source]¶

Injects fake companions in branches, at given radial distances.

Parameters:

array (3d/4d numpy ndarray) – Input cube. This is copied before the injections take place, so array is never modified.
psf_template (2d/3d numpy ndarray) – [for a 3D input array] 2d array with the normalized PSF template, with an odd or even shape. The PSF image must be centered wrt to the array. Therefore, it is recommended to run the function normalize_psf to generate a centered and flux-normalized PSF template. It can also be a 3D array, but length should match that of ADI cube. [for a 4D input array] In the ADI+mSDI case, it must be a 3d array (matching spectral dimensions).
angle_list (1d numpy ndarray) – List of parallactic angles, in degrees.
flevel (float or 1d array or 2d array) – Factor for controlling the brightness of the fake companions. If a float, the same flux is used for all injections. [3D input cube]: if a list/1d array is provided, it should have same length as number of frames in the 3D cube (can be used to take into account varying observing conditions or airmass). [4D (ADI+mSDI) input cube]: if a list/1d array should have the same length as the number of spectral channels (i.e. provide a spectrum). If a 2d array, it should be n_wavelength x n_frames (can e.g. be used to inject a spectrum in varying conditions).
plsc (float or None) – Value of the plsc in arcsec/px. Only used for printing debug output when verbose=True.
rad_dists (float, list or array 1d) – Vector of radial distances of fake companions in pixels.
n_branches (int, optional) – Number of azimutal branches.
theta (float, optional) – Angle in degrees for rotating the position of the first branch that by default is located at zero degrees. Theta counts counterclockwise from the positive x axis.
imlib (str, optional) – See the documentation of the vip_hci.preproc.frame_shift function.
interpolation (str, optional) – See the documentation of the vip_hci.preproc.frame_shift function.
transmission (numpy array, optional) – Radial transmission of the coronagraph, if any. Array with either 2 x n_rad, 1+n_ch x n_rad columns. The first column should contain the radial separation in pixels, while the other column(s) are the corresponding off-axis transmission (between 0 and 1), for either all, or each spectral channel (only relevant for a 4D input cube).
radial_gradient (bool, optional) – Whether to apply a radial gradient to the psf image at the moment of injection. By default False, i.e. the flux of the psf image is scaled only considering the value of tramnsmission at the exact radius the companion is injected. Setting it to False may better represent the transmission at the very edge of a physical mask though.
full_output (bool, optional) – Returns the x and y coordinates of the injections, additionally to the new array.
verbose (bool, optional) – If True prints out additional information.
nproc (int or None, optional) – Number of CPUs to use for multiprocessing. If None, will be automatically set to half the number of available CPUs.

Returns:

array_out (numpy ndarray) – Output array with the injected fake companions.
positions (list of tuple(y, x)) – [full_output] Coordinates of the injections in the first frame (and first wavelength for 4D cubes).
psf_trans (numpy ndarray) – [full_output & transmission != None] Array with injected psf affected by transmission (serves to check radial transmission)

vip_hci.fm.fakecomp.generate_cube_copies_with_injections(array, psf_template, angle_list, plsc, n_copies=100, inrad=8, outrad=12, dist_flux=('uniform', 2, 500))[source]¶

Create multiple copies of array with different random injections.

This is a wrapper around metrics.cube_inject_companions, which deals with multiple copies of the original data cube and generates random parameters.

Parameters:

array (3d/4d numpy ndarray) – Original input cube.
psf_template (2d/3d numpy ndarray) – Array with the normalized psf template. It should have an odd shape. It’s recommended to run the function normalize_psf to get a proper PSF template. In the ADI+mSDI case it must be a 3d array.
angle_list (1d numpy ndarray) – List of parallactic angles, in degrees.
plsc (float) – Value of the plsc in arcsec/px. Only used for printing debug output when verbose=True.
n_copies (int) – This is the number of ‘cube copies’ returned.
inrad,outrad (float) – Inner and outer radius of the injections. The actual injection position is chosen randomly.
dist_flux (tuple('method', params)) –
Tuple describing the flux selection. Method can be a function, the *params are passed to it. Method can also be a string, for a pre-defined random function:

('skewnormal', skew, mean, var)

uses scipy.stats.skewnorm.rvs

('uniform', low, high)

uses np.random.uniform

('normal', loc, scale)

uses np.random.normal

Returns:

fake_data – Represents a copy of the original array, with fake injections. The dictionary keys are:

cube

Array shaped like the input array, with the fake injections.

position : list of tuples(y,x)

List containing the positions of the injected companions, as (y,x) tuples.

dist : float

The distance of the injected companions, which was passed to cube_inject_companions.

theta : float, degrees

The initial angle, as passed to cube_inject_companions.

flux : float

The flux passed to cube_inject_companions.

Return type:

dict

vip_hci.fm.fakecomp.frame_inject_companion(array, array_fc, pos_y, pos_x, flux, imlib='vip-fft', interpolation='lanczos4')[source]¶

Injects a fake companion in a single frame (it could be a single multi-wavelength frame) at given coordinates, or in a cube (at the same coordinates, flux and with same fake companion image throughout the cube).

Parameters:	array (numpy ndarray, 2d or 3d) – Input frame or cube. array_fc (numpy ndarray, 2d) – Fake companion image to be injected. If even-dimensions, the center should be placed at coordinates [dim//2, dim//2] (0-based indexing), as per VIP’s convention. pos_x (pos_y,) – Y and X coordinates where the companion should be injected flux (int) – Flux at which the fake companion should be injected (i.e. scaling factor for the injected image) imlib (str, optional) – See documentation of the `vip_hci.preproc.frame_shift` function. interpolation (str, optional) – See documentation of the `vip_hci.preproc.frame_shift` function.
Returns:	array_out – Frame or cube with the companion injected
Return type:	numpy ndarray

vip_hci.fm.fakedisk module¶

Module with fake disk injection functions.

vip_hci.fm.fakedisk.cube_inject_fakedisk(fakedisk, angle_list, psf=None, **rot_options)[source]¶

Creates an ADI cube with a rotated synthetic disk image injected following input angle list.

Parameters:	fakedisk (numpy ndarray) – Input image of a fake disc angle_list (list) – Vector containing the parallactic angles. psf ((optional) the PSF to convolve the disk image with. It can be a) – small numpy.ndarray (we advise to use odd sizes to make sure the center s not shifted through the convolution). It forces normalization of the PSF to preserve the flux. It can also be a float representing the FWHM of the gaussian to be used for convolution. imlib (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. interpolation (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. cxy (tuple of int, optional) – Coordinates X,Y of the point with respect to which the rotation will be performed. By default the rotation is done with respect to the center of the frames, as it is returned by the function vip_hci.var.frame_center. nproc (int, optional) – Whether to rotate the frames in the sequence in a multi-processing fashion. Only useful if the cube is significantly large (frame size and number of frames). border_mode (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. rot_options (dictionary, optional) – Dictionary with optional keyword values for “nproc”, “cxy”, “imlib”, “interpolation, “border_mode”, “mask_val”, “edge_blend”, “interp_zeros”, “ker” (see documentation of `vip_hci.preproc.frame_rotate`)
Returns:	fakedisk_cube – Resulting cube with the fake disc inserted at the correct angles and convolved with the psf if a psf was provided.
Return type:	numpy ndarray

Example

import numpy as np

fakedisk = np.zeros((200,200))
fakedisk[:,99:101] = 1
angle_list = np.arange(10)
c = create_fakedisk_cube(fakedisk, angle_list, psf=None, imlib='vip-fft',
                         interpolation='lanczos4',cxy=None, nproc=1,
                         border_mode='constant')

vip_hci.fm.fakedisk.cube_inject_trace(array, psf_template, angle_list, flevel, rad_dists, theta, plsc=0.01225, n_branches=1, imlib='vip-fft', interpolation='lanczos4', verbose=True)[source]¶

Injects fake companions along a trace, such as a spiral. The trace is provided by 2 arrays corresponding to the polar coordinates where the companions will be located in the final derotated frame. Note: for a continuous-looking trace, and for an easier scaling using parameter ‘flevel’, it is recommended to separate the points of the trace by a distance of FWHM/2.

Parameters:	array (numpy ndarray) – Input 3D cube in which the extended feature is injected. psf_template (numpy ndarray) – 2d array with the normalized psf template. It should have an odd shape. It is recommended to run the function psf_norm to get a proper PSF template. flevel (float) – Flux at which the fake companions are injected into the cube along the trace. rad_dists (list or array 1d) – Vector of radial distances where the trace is to be injected. theta (list or array 1d) – Vector of angles (deg) where the trace is to be injected (trigonometric angles, NOT PA East from North). plsc (float, opt) – Value of the plate scale in arcsec/pixel (optional, will only be used if verbose is True). n_branches (int, optional) – Number of azimutal branches on which the trace is injected. imlib ({'opencv', 'ndimage-fourier', 'ndimage-interp', 'vip-fft'}, str opt) – Library or method used for image operations (shifts). interpolation (str, optional) – Interpolation method. Check documentation of the function `vip_hci.preproc.frame_shift`. verbose ({True, False}, bool optional) – If True prints out additional information.
Returns:	array_out – Output array with the injected fake companions.
Return type:	numpy ndarray

vip_hci.fm.negfc_fmerit module¶

Module with the function of merit definitions for the NEGFC optimization.

vip_hci.fm.negfc_fmerit.get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm, r_guess, theta_guess, cube_ref=None, wedge=None, svd_mode='lapack', scaling=None, algo=<function pca_annulus>, delta_rot=1, imlib='vip-fft', interpolation='lanczos4', collapse='median', weights=None, algo_options={})[source]¶

Extracts the mean and standard deviation of pixel intensities in an annulus of the PCA-ADI image obtained with ‘algo’, in the part of a defined wedge that is not overlapping with PA_pl+-delta_PA.

Parameters:	cube (numpy.array) – The cube of fits images expressed as a numpy.array. angs (numpy.array) – The parallactic angle fits image expressed as a numpy.array. ncomp (int or None) – The number of principal components for PCA-based algorithms. annulus_width (float) – The width in pixels of the annulus on which the PCA is performed. aperture_radius (float) – The radius in fwhm of the circular aperture. fwhm (float) – Value of the FWHM of the PSF. r_guess (float) – The radial position of the center of the circular aperture. This parameter is NOT the radial position of the candidate associated to the Markov chain, but should be the fixed initial guess. theta_guess (float) – The angular position of the center of the circular aperture. This parameter is NOT the angular position of the candidate associated to the Markov chain, but should be the fixed initial guess. cube_ref (numpy ndarray, 3d, optional) – Reference library cube. For Reference Star Differential Imaging. wedge (tuple, opt) – Range in theta where the mean and standard deviation are computed in an annulus defined in the PCA image. If None, it will be calculated automatically based on initial guess and derotation angles to avoid. If some disc signal is present elsewhere in the annulus, it is recommended to provide wedge manually. The provided range should be continuous and >0. E.g. provide (270, 370) to consider a PA range between [-90,+10]. svd_mode ({'lapack', 'randsvd', 'eigen', 'arpack'}, str optional) – Switch for different ways of computing the SVD and selected PCs. scaling ({None, 'temp-mean', 'temp-standard'}) – With None, no scaling is performed on the input data before SVD. With “temp-mean” then temporal px-wise mean subtraction is done and with “temp-standard” temporal mean centering plus scaling to unit variance is done. algo (python routine, opt {pca_annulus, pca_annular, pca, custom}) – Routine to be used to model and subtract the stellar PSF. From an input cube, derotation angles, and optional arguments, it should return a post-processed frame. delta_rot (float, optional) – If algo is set to pca_annular, delta_rot is the angular threshold used to select frames in the PCA library (see description of pca_annular). imlib (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. interpolation (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. collapse ({'median', 'mean', 'sum', 'trimmean', None}, str or None, optional) – Sets the way of collapsing the frames for producing a final image. If None then the cube of residuals is returned. weights (1d numpy array or list, optional) – Weights to be applied for a weighted mean. Need to be provided if collapse mode is ‘wmean’. algo_options (dict, opt) – Dictionary with additional parameters related to the algorithm (e.g. tol, min_frames_lib, max_frames_lib). If ‘algo’ is not a vip routine, this dict should contain all necessary arguments apart from the cube and derotation angles. Note: arguments such as ncomp, svd_mode, scaling, imlib, interpolation or collapse can also be included in this dict (the latter are also kept as function arguments for consistency with older versions of vip).
Returns:	values – The pixel values in the circular aperture after the PCA process.
Return type:	numpy.array

vip_hci.fm.negfc_mcmc module¶

Module with the MCMC (emcee) sampling for NEGFC parameter estimation.

[CHR21]

Christiaens et al. 2021
A faint companion around CrA-9: protoplanet or obscured binary?
MNRAS, Volume 502, Issue 4, pp. 6117-6139
https://arxiv.org/abs/2102.10288

[FOR13]

Foreman-Mackey et al. 2013
emcee: The MCMC Hammer
PASP, Volume 125, Issue 925, p. 306
https://arxiv.org/abs/1202.3665

[WER17]

(1, 2, 3, 4, 5, 6)

Wertz et al. 2017
VLT/SPHERE robust astrometry of the HR8799 planets at milliarcsecond-level accuracy. Orbital architecture analysis with PyAstrOFit
Astronomy & Astrophysics, Volume 598, p. 83
https://arxiv.org/abs/1610.04014

[GOO10]

(1, 2)

Goodman & Weare 2010
Ensemble samplers with affine invariance
Comm. App. Math. Comp. Sci., Vol. 5, Issue 1, pp. 65-80.
https://ui.adsabs.harvard.edu/abs/2010CAMCS…5…65G

vip_hci.fm.negfc_mcmc.mcmc_negfc_sampling(cube, angs, psfn, initial_state, algo=<function pca_annulus>, ncomp=1, annulus_width=8, aperture_radius=1, fwhm=4, mu_sigma=True, sigma='spe+pho', force_rPA=False, fmerit='sum', cube_ref=None, svd_mode='lapack', scaling=None, delta_rot=1, imlib='vip-fft', interpolation='lanczos4', collapse='median', algo_options={}, wedge=None, weights=None, transmission=None, nwalkers=100, bounds=None, a=2.0, burnin=0.3, rhat_threshold=1.01, rhat_count_threshold=1, niteration_min=10, niteration_limit=10000, niteration_supp=0, check_maxgap=20, conv_test='ac', ac_c=50, ac_count_thr=3, nproc=1, output_dir='results/', output_file=None, display=False, verbosity=0, save=False)[source]¶

Runs an affine invariant mcmc sampling algorithm in order to determine the position and the flux of the planet using the ‘Negative Fake Companion’ technique. The result of this procedure is a chain with the samples from the posterior distributions of each of the 3 parameters.

This technique can be summarized as follows: 1) We inject a negative fake companion (one candidate) at a given position and characterized by a given flux, both close to the expected values. 2) We run PCA on an full annulus which pass through the initial guess, regardless of the position of the candidate. 3) We extract the intensity values of all the pixels contained in a circular aperture centered on the initial guess. 4) We calculate a function of merit \(\chi^2\) (see below). The steps 1) to 4) are then looped. At each iteration, the candidate model parameters are defined by the emcee Affine Invariant algorithm [FOR13].

There are different possibilities for the figure of merit (step 4):

mu_sigma=None; fmerit='sum' (as in [WER17]):\(\chi^2 = \sum(\|I_j\|)\)
mu_sigma=None; fmerit='stddev' (likely more appropriate than ‘sum’ when residual speckle noise is still significant):\(\chi^2 = N \sigma_{I_j}(values,ddof=1)*values.size\)
mu_sigma=True or a tuple (as in [CHR21], new default):\(\chi^2 = \sum\frac{(I_j- mu)^2}{\sigma^2}\)

where \(j \in {1,...,N}\) with N the total number of pixels contained in the circular aperture, \(\sigma_{I_j}\) is the standard deviation of \(I_j\) values, and \(\mu\) is the mean pixel intensity in a truncated annulus at the radius of the companion candidate (i.e. excluding the cc region). See description of mu_sigma and sigma for more details on \(\sigma\) considered in the last equation.

Speed tricks:

crop your input cube to a size such as to just include the annulus on which the PCA is performed;
set imlib='opencv' (much faster image rotations, BUT at the expense of flux conservation);
increase nproc (if your machine allows);
reduce ac_c (or increase rhat_threshold if conv_test='gb') for a faster convergence);
reduce niteration_limit to force the sampler to stop even if it has not reached convergence.

Parameters:	cube (3d or 4d numpy ndarray) – Input ADI or ADI+IFS cube. angs (numpy.array) – The parallactic angle vector. psfn (numpy 2D or 3D array) – Normalised PSF template used for negative fake companion injection. The PSF must be centered and the flux in a 1xFWHM aperture must equal 1 (use `vip_hci.metrics.normalize_psf`). If the input cube is 3D and a 3D array is provided, the first dimension must match for both cubes. This can be useful if the star was unsaturated and conditions were variable. If the input cube is 4D, psfn must be either 3D or 4D. In either cases, the first dimension(s) must match those of the input cube. initial_state (tuple or 1d numpy ndarray) – The first guess for the position and flux(es) of the planet, in that order. In the case of a 3D input cube, it should have a length of 3, while in the case of a 4D input cube, the length should be 2 + number of spectral channels (one flux estimate per wavelength). Each walker will start in a small ball around this preferred position. algo (python routine, optional) – Post-processing algorithm used to model and subtract the star. First 2 arguments must be input cube and derotation angles. Must return a post-processed 2d frame. ncomp (int or None, optional) – The number of principal components for PCA-based algorithms. annulus_width (float, optional) – The width in pixels of the annulus on which the PCA is performed. aperture_radius (float, optional) – The radius in FWHM of the circular aperture. fwhm (float) – The FHWM in pixels. mu_sigma (tuple of 2 floats or bool, opt) – If set to None: not used, and falls back to original version of the algorithm, using `fmerit` [WER17]. If a tuple of 2 elements: should be the mean and standard deviation of 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. These values will then be used in the log-probability of the MCMC. sigma (str, opt) – Sets the type of noise to be included as sigma^2 in the log-probability expression. Choice between ‘pho’ for photon (Poisson) noise, ‘spe’ for residual (mostly whitened) speckle noise, or ‘spe+pho’ for both. force_rPA (bool, optional) – Whether to only search for optimal flux, provided (r,PA). fmerit ({'sum', 'stddev'}, string optional) – If mu_sigma is not provided nor set to True, this parameter determines which figure of merit to be used among the 2 possibilities implemented in [WER17]. ‘stddev’ may work well for point like sources surrounded by extended signals. cube_ref (3d or 4d numpy ndarray, or list of 3d ndarray, optional) – Reference library cube for Reference Star Differential Imaging. Should be 3d, except if the input cube is 4d, in which case it can either be a 4d array or a list of 3d ndarray (reference cube for each wavelength). svd_mode ({'lapack', 'randsvd', 'eigen', 'arpack'}, str optional) – Switch for different ways of computing the SVD and selected PCs. ‘randsvd’ is not recommended for the negative fake companion technique. scaling ({'temp-mean', 'temp-standard'} or None, optional) – With None, no scaling is performed on the input data before SVD. With “temp-mean” then temporal px-wise mean subtraction is done and with “temp-standard” temporal mean centering plus scaling to unit variance is done. delta_rot (float, optional) – If algo is set to pca_annular, delta_rot is the angular threshold used to select frames in the PCA library (see description of pca_annular). Increases processing time. imlib (str, optional) – Imlib used for both image rotation and sub-px shift: ”opencv”: will use it for both; ”skimage” or “ndimage-interp” will use scikit-image and scipy.ndimage for rotation and shift resp.; ”ndimage-fourier” or “vip-fft” will use Fourier transform based methods for both. interpolation (str, optional) – Interpolation order. See the documentation of the `vip_hci.preproc.frame_rotate` function. Note that the interpolation options are identical for rotation and shift within each of the 3 imlib cases above. collapse ({'median', 'mean', 'sum', 'trimmean', 'wmean'}, str, optional) – Sets the way of collapsing the frames for producing a final image. If None then the cube of residuals is used when measuring the function of merit (instead of a single final frame). algo_options (dict, opt) – Dictionary with additional parameters related to the algorithm (e.g. tol, min_frames_lib, max_frames_lib). If ‘algo’ is not a vip routine, this dict should contain all necessary arguments apart from the cube and derotation angles. Note: arguments such as ncomp, svd_mode, scaling, imlib, interpolation or collapse can also be included in this dict (the latter are also kept as function arguments for consistency with older versions of vip). wedge (tuple, opt) – Range in theta where the mean and standard deviation are computed in an annulus defined in the PCA image. If None, it will be calculated automatically based on initial guess and derotation angles to avoid negative ADI side lobes. If some disc signal is present elsewhere in the annulus, it is recommended to provide wedge manually. The provided range should be continuous and >0. E.g. provide (270, 370) to consider a PA range between [-90,+10]. weights (1d array, optional) – If provided, the negative fake companion fluxes will be scaled according to these weights before injection in the cube. Can reflect changes in the observing conditions throughout the sequence. transmission (numpy array, optional) – Radial transmission of the coronagraph, if any. Array with either 2 x n_rad, 1+n_ch x n_rad columns. The first column should contain the radial separation in pixels, while the other column(s) are the corresponding off-axis transmission (between 0 and 1), for either all, or each spectral channel (only relevant for a 4D input cube). nwalkers (int optional) – The number of [GOO10] ‘walkers’. bounds (numpy.array or list, default=None, optional) – The prior knowledge on the model parameters. If None, large bounds will be automatically estimated from the initial state. a (float, default=2.0) – The proposal scale parameter. See Note. burnin (float, default=0.3) – The fraction of a walker chain which is discarded. Note: only used for Gelman-Rubin convergence test - the chains are returned full. rhat_threshold (float, default=0.01) – The Gelman-Rubin threshold used for the test for nonconvergence. rhat_count_threshold (int, optional) – The Gelman-Rubin test must be satisfied ‘rhat_count_threshold’ times in a row before claiming that the chain has converged. conv_test (str, optional {'gb','ac'}) – Method to check for convergence: ’gb’ for gelman-rubin test (http://digitalassets.lib.berkeley.edu/sdtr/ucb/text/305.pdf) ’ac’ for autocorrelation analysis (https://emcee.readthedocs.io/en/stable/tutorials/autocorr/) ac_c (float, optional) – If the convergence test is made using the auto-correlation, this is the value of C such that tau/N < 1/C is the condition required for tau to be considered a reliable auto-correlation time estimate (for N number of samples). Recommended: C>50. More details here: https://emcee.readthedocs.io/en/stable/tutorials/autocorr/ ac_count_thr (int, optional) – The auto-correlation test must be satisfied ac_count_thr times in a row before claiming that the chain has converged. niteration_min (int, optional) – Steps per walker lower bound. The simulation will run at least this number of steps per walker. niteration_limit (int, optional) – Steps per walker upper bound. If the simulation runs up to ‘niteration_limit’ steps without having reached the convergence criterion, the run is stopped. niteration_supp (int, optional) – Number of iterations to run after having “reached the convergence”. check_maxgap (int, optional) – Maximum number of steps per walker between two Gelman-Rubin test. nproc (int or None, optional) – The number of processes to use for parallelization. If None, will be set automatically to half the number of CPUs available. output_dir (str, optional) – The name of the output directory which contains the output files in the case `save` is True. output_file (str, optional) – The name of the output file which contains the MCMC results in the case `save` is True. display (bool, optional) – If True, the walk plot is displayed at each evaluation of the Gelman- Rubin test. verbosity (0, 1, 2 or 3, optional) – Verbosity level. 0 for no output and 3 for full information. (only difference between 2 and 3 is that 3 also writes intermediate pickles containing the state of the chain at convergence tests; these can end up taking a lot of space). save (bool, optional) – If True, the MCMC results are pickled.
Returns:	out – The MCMC chain.
Return type:	numpy.array

Note

The parameter a must be > 1. For more theoretical information concerning this parameter, see [GOO10].

The parameter rhat_threshold can be a numpy.array with individual threshold value for each model parameter.

vip_hci.fm.negfc_mcmc.chain_zero_truncated(chain)[source]¶

Return the Markov chain with the dimension: walkers x steps* x parameters, where steps* is the last step before having 0 (not yet constructed chain).

Parameters:	chain (numpy.array) – The MCMC chain.
Returns:	out – The truncated MCMC chain, that is to say, the chain which only contains relevant information.
Return type:	numpy.array

vip_hci.fm.negfc_mcmc.show_corner_plot(chain, burnin=0.5, save=False, output_dir='', **kwargs)[source]¶

Display or save a figure showing the corner plot (pdfs + correlation plots)

Parameters:	chain (numpy.array) – The Markov chain. The shape of chain must be nwalkers x length x dim. If a part of the chain is filled with zero values, the method will discard these steps. burnin (float, default: 0) – The fraction of a walker chain we want to discard. save (boolean, default: False) – If True, a pdf file is created. output_dir (str, optional) – The name of the output directory which contains the output files in the case `save` is True. kwargs – Additional attributes passed to the corner.corner() method.
Returns:	(Displays the figure or create a pdf file named walk_plot.pdf in the working directory).
Return type:	None
Raises:	`ImportError`

vip_hci.fm.negfc_mcmc.show_walk_plot(chain, save=False, output_dir='', **kwargs)[source]¶

Display or save a figure showing the path of each walker during the MCMC run

Parameters:

chain (numpy.array) – The Markov chain. The shape of chain must be nwalkers x length x dim. If a part of the chain is filled with zero values, the method will discard these steps.
save (boolean, default: False) – If True, a pdf file is created.
output_dir (str, optional) – The name of the output directory which contains the output files in the case save is True.
kwargs – Additional attributes are passed to the matplotlib plot method.

Returns:

Display the figure or create a pdf file named walk_plot.pdf in the working
directory.

vip_hci.fm.negfc_mcmc.confidence(isamples, cfd=68.27, bins=100, gaussian_fit=False, weights=None, verbose=True, save=False, output_dir='', force=False, output_file='confidence.txt', title=None, ndig=1, plsc=None, labels=['r', 'theta', 'f'], gt=None, **kwargs)[source]¶

Determine the highly probable value for each model parameter, as well as the 1-sigma confidence interval.

Parameters:

isamples (numpy.array) – The independent samples for each model parameter.
cfd (float, optional) – The confidence level given in percentage.
bins (int, optional) – The number of bins used to sample the posterior distributions.
gaussian_fit (boolean, optional) – If True, a gaussian fit is performed in order to determine (\(\mu, \sigma\)).
weights ((n, ) numpy ndarray or None, optional) – An array of weights for each sample.
verbose (boolean, optional) – Display information in the shell.
save (boolean, optional) – If “True”, a txt file with the results is saved in the output repository.
output_dir (str, optional) – If save is True, this is the full path to a directory where the results are saved.
force (bool, optional) – If set to True, force the confidence interval estimate even if too many samples fall in a single bin (unreliable CI estimates). If False, an error message is raised if the percentile of samples falling in a single bin is larger than cfd, suggesting to increase number of bins.
output_file (str, optional) – If save is True, name of the text file in which the results are saved.
title (bool or str, optional) – If not None, will print labels and parameter values on top of each plot. If a string, will print that label in front of the parameter values.
ndig (int or list of int, optional) – If title is not None, this sets the number of significant digits to be used when printing the best estimate for each parameter. If int, the same number of digits is used for all parameters. If a list of int, the length should match the number of parameters.
plsc (float, optional) – If save is True, this is used to convert pixels to arcsec when writing results for r.
labels (list of strings, optional) – Labels to be used for both plots and dictionaries. Length should match the number of free parameters.
gt (list or 1d numpy array, optional) – Ground truth values (when known). If provided, will also be plotted for reference.

Returns:

out –

[gaussian_fit=False] i) the highly probable solutions (dictionary),

the respective confidence interval (dict.);

[gaussian_fit=True] i) the center of the best-fit 1d Gaussian

distributions (tuple of 3 floats), and ii) their standard deviation, for each parameter

Return type:

A 2 elements tuple with either:

vip_hci.fm.negfc_nested module¶

Module with functions for posterior sampling of the NEGFC parameters using nested sampling (nestle).

[BAR13]

K. Barbary 2013
nestle
GitHub repository
https://github.com/kbarbary/nestle

[FER09]

Feroz et al. 2009
MULTINEST: an efficient and robust Bayesian inference tool for cosmology
and particle physics
MNRAS, Volume 398, Issue 4, pp. 1601-1614
https://arxiv.org/abs/0809.3437

[MUK06]

Mukherjee et al. 2006
A Nested Sampling Algorithm for Cosmological Model Selection
ApJL, Volume 638, Issue 2, pp. 51-54
https://arxiv.org/abs/astro-ph/0508461

[SKI04]

Skilling 2004
Bayesian Inference and Maximum Entropy Methods in Science and Engineering:
24th International Workshop on Bayesian Inference and Maximum Entropy
Methods in Science and Engineering
American Institute of Physics Conference Series, Volume 735, pp. 395-405
https://ui.adsabs.harvard.edu/abs/2004AIPC..735..395S

vip_hci.fm.negfc_nested.nested_negfc_sampling(init, cube, angs, psfn, fwhm, mu_sigma=True, sigma='spe+pho', fmerit='sum', annulus_width=8, aperture_radius=1, ncomp=10, scaling=None, svd_mode='lapack', cube_ref=None, collapse='median', algo=<function pca_annulus>, delta_rot=1, algo_options={}, weights=None, w=(5, 5, 200), method='single', npoints=100, dlogz=0.1, decline_factor=None, rstate=None, verbose=True)[source]¶

Runs a nested sampling algorithm with nestle [BAR13] in order to determine the position and the flux of the planet using the ‘Negative Fake Companion’ technique. The result of this procedure is a nestle object containing the samples from the posterior distributions of each of the 3 parameters. It provides good results (value plus error bars) compared to a more CPU intensive Monte Carlo approach with the affine invariant sampler (emcee).

Parameters:	init (numpy ndarray or tuple of length 3) – The first guess for the position and flux of the planet, respectively. It serves for generating the bounds of the log prior function (uniform in a bounded interval). cube (3d or 4d numpy ndarray) – Input ADI or ADI+IFS cube. angs (numpy.array) – The parallactic angle vector. psfn (numpy 2D or 3D array) – Normalised PSF template used for negative fake companion injection. The PSF must be centered and the flux in a 1xFWHM aperture must equal 1 (use `vip_hci.metrics.normalize_psf`). If the input cube is 3D and a 3D array is provided, the first dimension must match for both cubes. This can be useful if the star was unsaturated and conditions were variable. If the input cube is 4D, psfn must be either 3D or 4D. In either cases, the first dimension(s) must match those of the input cube. fwhm (float) – The FHWM in pixels. mu_sigma (tuple of 2 floats or bool, opt) – If set to None: not used, and falls back to original version of the algorithm, using fmerit [WER17]. If a tuple of 2 elements: should be the mean and standard deviation of 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. These values will then be used in the log-probability of the MCMC. sigma (str, opt) – Sets the type of noise to be included as sigma^2 in the log-probability expression. Choice between ‘pho’ for photon (Poisson) noise, ‘spe’ for residual (mostly whitened) speckle noise, or ‘spe+pho’ for both. force_rPA (bool, optional) – Whether to only search for optimal flux, provided (r,PA). fmerit ({'sum', 'stddev'}, string optional) – If mu_sigma is not provided nor set to True, this parameter determines which figure of merit to be used among the 2 possibilities implemented in [WER17]. ‘stddev’ may work well for point like sources surrounded by extended signals. annulus_width (float, optional) – The width in pixel of the annulus on which the PCA is performed. aperture_radius (float, optional) – The radius of the circular aperture in FWHM. ncomp (int optional) – The number of principal components. scaling ({'temp-mean', 'temp-standard'} or None, optional) – With None, no scaling is performed on the input data before SVD. With “temp-mean” then temporal px-wise mean subtraction is done and with “temp-standard” temporal mean centering plus scaling to unit variance is done. svd_mode ({'lapack', 'randsvd', 'eigen', 'arpack'}, str optional) – Switch for different ways of computing the SVD and selected PCs. cube_ref (numpy ndarray, 3d, optional) – Reference library cube. For Reference Star Differential Imaging. collapse ({'median', 'mean', 'sum', 'trimmean', None}, str or None, opt) – Sets the way of collapsing the frames for producing a final image. If None then the cube of residuals is used when measuring the function of merit (instead of a single final frame). algo_options (dict, opt) – Dictionary with additional parameters for the algorithm (e.g. tol, min_frames_lib, max_frames_lib) weights (1d array, optional) – If provided, the negative fake companion fluxes will be scaled according to these weights before injection in the cube. Can reflect changes in the observing conditions throughout the sequence. w (tuple of length 3) – The size of the bounds (around the initial state `init`) for each parameter. method ({"single", "multi", "classic"}, str optional) – Flavor of nested sampling. Single ellipsoid works well for the NEGFC and is the default. npoints (int optional) – Number of active points. At least ndim+1 (4 will produce bad results). For problems with just a few parameters (<=5) like the NEGFC, good results are obtained with 100 points (default). dlogz (Estimated remaining evidence) – Iterations will stop when the estimated contribution of the remaining prior volume to the total evidence falls below this threshold. Explicitly, the stopping criterion is log(z + z_est) - log(z) < dlogz where z is the current evidence from all saved samples, and z_est is the estimated contribution from the remaining volume. This option and decline_factor are mutually exclusive. If neither is specified, the default is dlogz=0.5. decline_factor (float, optional) – If supplied, iteration will stop when the weight (likelihood times prior volume) of newly saved samples has been declining for decline_factor * nsamples consecutive samples. A value of 1.0 seems to work pretty well. rstate (random instance, optional) – RandomState instance. If not given, the global random state of the numpy.random module will be used.
Returns:	res – `Nestle` object with the nested sampling results, including the posterior samples.
Return type:	nestle object

Note

Nested Sampling is a computational approach for integrating posterior probability in order to compare models in Bayesian statistics. It is similar to Markov Chain Monte Carlo (MCMC) in that it generates samples that can be used to estimate the posterior probability distribution. Unlike MCMC, the nature of the sampling also allows one to calculate the integral of the distribution. It also happens to be a pretty good method for robustly finding global maxima.

Nestle documentation: http://kbarbary.github.io/nestle/

Convergence: http://kbarbary.github.io/nestle/stopping.html Nested sampling has no well-defined stopping point. As iterations continue, the active points sample a smaller and smaller region of prior space. This can continue indefinitely. Unlike typical MCMC methods, we don’t gain any additional precision on the results by letting the algorithm run longer; the precision is determined at the outset by the number of active points. So, we want to stop iterations as soon as we think the active points are doing a pretty good job sampling the remaining prior volume - once we’ve converged to the highest-likelihood regions such that the likelihood is relatively flat within the remaining prior volume.

Method: The trick in nested sampling is to, at each step in the algorithm, efficiently choose a new point in parameter space drawn with uniform probability from the parameter space with likelihood greater than the current likelihood constraint. The different methods all use the current set of active points as an indicator of where the target parameter space lies, but differ in how they select new points from it:

“classic” is close to the method described in [SKI04].

“single”, [MUK06], Determines a single ellipsoid that bounds all active points, enlarges the ellipsoid by a user-settable factor, and selects a new point at random from within the ellipsoid.

“multiple”, [FER09] (Multinest). In cases where the posterior is multi-modal, the single-ellipsoid method can be extremely inefficient. In such situations, there are clusters of active points on separate high-likelihood regions separated by regions of lower likelihood. Bounding all points in a single ellipsoid means that the ellipsoid includes the lower-likelihood regions we wish to avoid sampling from. The solution is to detect these clusters and bound them in separate ellipsoids. For this, we use a recursive process where we perform K-means clustering with K=2. If the resulting two ellipsoids have a significantly lower total volume than the parent ellipsoid (less than half), we accept the split and repeat the clustering and volume test on each of the two subset of points. This process continues recursively. Alternatively, if the total ellipse volume is significantly greater than expected (based on the expected density of points) this indicates that there may be more than two clusters and that K=2 was not an appropriate cluster division. We therefore still try to subdivide the clusters recursively. However, we still only accept the final split into N clusters if the total volume decrease is significant.

vip_hci.fm.negfc_nested.nested_sampling_results(ns_object, burnin=0.4, bins=None, cfd=68.27, save=False, output_dir='/', plot=False)[source]¶

Shows the results of the Nested Sampling, summary, parameters with errors, walk and corner plots.

Parameters:	ns_object (numpy.array) – The nestle object returned from nested_spec_sampling. burnin (float, default: 0) – The fraction of a walker we want to discard. bins (int, optional) – The number of bins used to sample the posterior distributions. cfd (float, optional) – The confidence level given in percentage. save (boolean, default: False) – If True, a pdf file is created. output_dir (str, optional) – The name of the output directory which contains the output files in the case `save` is True. plot (bool, optional) – Whether to show the plots (instead of saving them).
Returns:	final_res – Best-fit parameters and uncertainties (corresponding to 68% confidence interval). Dimensions: nparams x 2.
Return type:	numpy ndarray

vip_hci.fm.negfc_simplex module¶

Module with simplex (Nelder-Mead) optimization for defining the flux and position of a companion using the Negative Fake Companion.

vip_hci.fm.negfc_simplex.firstguess(cube, angs, psfn, ncomp, planets_xy_coord, fwhm=4, annulus_width=4, aperture_radius=1, cube_ref=None, svd_mode='lapack', scaling=None, fmerit='sum', imlib='vip-fft', interpolation='lanczos4', collapse='median', algo=<function pca_annulus>, delta_rot=1, p_ini=None, f_range=None, transmission=None, mu_sigma=True, wedge=None, weights=None, force_rPA=False, algo_options={}, simplex=True, simplex_options=None, plot=False, verbose=True, save=False)[source]¶

Determines a first guess for the position and the flux of a planet, as explained in [WER17].

We process the cube without injecting any negative fake companion. This leads to the visual detection of the planet(s). For each of them, one can estimate the (x,y) coordinates in pixel for the position of the star, as well as the planet(s).

From the (x,y) coordinates in pixels for the star and planet(s), we can estimate a preliminary guess for the position and flux for each planet by using the method “firstguess_from_coord”. The argument “f_range” allows to indicate prior limits for the flux (optional, default: None). This step can be reiterate to refine the preliminary guess for the flux.

We can go a step further by using a Simplex Nelder_Mead minimization to estimate the first guess based on the preliminary guess.

Parameters:	cube (3d or 4d numpy ndarray) – Input ADI or ADI+IFS cube. angs (numpy.array) – The parallactic angle fits image expressed as a numpy.array. psfn (numpy 2D or 3D array) – Normalised PSF template used for negative fake companion injection. The PSF must be centered and the flux in a 1xFWHM aperture must equal 1 (use `vip_hci.metrics.normalize_psf`). If the input cube is 3D and a 3D array is provided, the first dimension must match for both cubes. This can be useful if the star was unsaturated and conditions were variable. If the input cube is 4D, psfn must be either 3D or 4D. In either cases, the first dimension(s) must match those of the input cube. ncomp (int or 1d numpy array of int) – The number of principal components. If cube is a 4D cube, ncomp can be a list of integers, with length matching the first dimension of the cube. planets_xy_coord (array or list) – The list of (x,y) positions of the planets. plsc (float, optional) – The platescale, in arcsec per pixel. fwhm (float, optional) – The FHWM in pixels. annulus_width (int, optional) – The width in pixels of the annulus on which the PCA is done. aperture_radius (int, optional) – The radius of the circular aperture in terms of the FWHM. cube_ref (numpy ndarray, 3d, optional) – Reference library cube. For Reference Star Differential Imaging. svd_mode ({'lapack', 'randsvd', 'eigen', 'arpack'}, str optional) – Switch for different ways of computing the SVD and selected PCs. scaling ({'temp-mean', 'temp-standard'} or None, optional) – With None, no scaling is performed on the input data before SVD. With “temp-mean” then temporal px-wise mean subtraction is done and with “temp-standard” temporal mean centering plus scaling to unit variance is done. fmerit ({'sum', 'stddev'}, string optional) – Figure of merit to be used, if mu_sigma is set to None. imlib (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. interpolation (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. collapse ({'median', 'mean', 'sum', 'trimmean', None}, str or None, opt) – Sets the way of collapsing the frames for producing a final image. If None then the cube of residuals is used when measuring the function of merit (instead of a single final frame). algo (python routine, opt {pca_annulus, pca_annular, pca, custom}) – Routine to be used to model and subtract the stellar PSF. From an input cube, derotation angles, and optional arguments, it should return a post-processed frame. delta_rot (float, optional) – If algo is set to pca_annular, delta_rot is the angular threshold used to select frames in the PCA library (see description of pca_annular). p_ini (numpy.array) – Position (r, theta) of the circular aperture center. f_range (numpy.array, optional) – The range of flux tested values. If None, 20 values between 0 and 5000 are tested. transmission (numpy array, optional) – Array with 2 columns. First column is the radial separation in pixels. Second column is the off-axis transmission (between 0 and 1) at the radial separation given in column 1. mu_sigma (tuple of 2 floats, bool or None, opt) – 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 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. wedge (tuple, opt) – Range in theta where the mean and standard deviation are computed in an annulus defined in the PCA image. If None, it will be calculated automatically based on initial guess and derotation angles to avoid. If some disc signal is present elsewhere in the annulus, it is recommended to provide wedge manually. The provided range should be continuous and >0. E.g. provide (270, 370) to consider a PA range between [-90,+10]. weights (1d array, optional) – If provided, the negative fake companion fluxes will be scaled according to these weights before injection in the cube. Can reflect changes in the observing conditions throughout the sequence. force_rPA (bool, optional) – Whether to only search for optimal flux, provided (r,PA). algo_options (dict, opt) – Dictionary with additional parameters for the pca algorithm (e.g. tol, min_frames_lib, max_frames_lib). Note: arguments such as svd_mode, scaling imlib, interpolation or collapse can also be included in this dict (the latter are also kept as function arguments for compatibility with older versions of vip). simplex (bool, optional) – If True, the Nelder-Mead minimization is performed after the flux grid search. simplex_options (dict, optional) – The scipy.optimize.minimize options. plot (boolean, optional) – If True, the figure chi2 vs. flux is displayed. verbose (bool, optional) – If True, display intermediate info in the shell. save (bool, optional) – If True, the figure chi2 vs. flux is saved.
Returns:	out – The polar coordinates and the flux(es) of the companion.
Return type:	tuple of 3+ elements

Note

Polar angle is not the conventional NORTH-TO-EAST P.A., but the counter-clockwise angle measured from the positive x axis.

vip_hci.fm.negfc_simplex.firstguess_from_coord(planet, center, cube, angs, psfn, fwhm, annulus_width, aperture_radius, ncomp, cube_ref=None, svd_mode='lapack', scaling=None, fmerit='sum', imlib='vip-fft', interpolation='lanczos4', collapse='median', algo=<function pca_annulus>, delta_rot=1, algo_options={}, f_range=None, transmission=None, mu_sigma=(0, 1), weights=None, plot=False, verbose=True, save=False, debug=False, full_output=False)[source]¶

Determine a first guess for the flux of a companion at a given position in the cube by doing a simple grid search evaluating the reduced chi2.

Parameters:

planet (numpy.array) – The (x,y) position of the planet in the pca processed cube.
center (numpy.array) – The (x,y) position of the cube center.
cube (3d or 4d numpy ndarray) – Input ADI or ADI+IFS cube.
angs (numpy.array) – The parallactic angle fits image expressed as a numpy.array.
psfn (numpy 2D or 3D array) – Normalised PSF template used for negative fake companion injection. The PSF must be centered and the flux in a 1xFWHM aperture must equal 1 (use vip_hci.metrics.normalize_psf). If the input cube is 3D and a 3D array is provided, the first dimension must match for both cubes. This can be useful if the star was unsaturated and conditions were variable. If the input cube is 4D, psfn must be either 3D or 4D. In either cases, the first dimension(s) must match those of the input cube.
fwhm (float) – The FHWM in pixels.
annulus_width (int, optional) – The width in pixels of the annulus on which the PCA is done.
aperture_radius (int, optional) – The radius of the circular aperture in terms of the FWHM.
ncomp (int) – The number of principal components.
cube_ref (numpy ndarray, 3d, optional) – Reference library cube. For Reference Star Differential Imaging.
svd_mode ({'lapack', 'randsvd', 'eigen', 'arpack'}, str optional) – Switch for different ways of computing the SVD and selected PCs.
scaling ({'temp-mean', 'temp-standard'} or None, optional) – With None, no scaling is performed on the input data before SVD. With “temp-mean” then temporal px-wise mean subtraction is done and with “temp-standard” temporal mean centering plus scaling to unit variance is done.
fmerit ({'sum', 'stddev'}, string optional) – Figure of merit to be used, if mu_sigma is set to None.
imlib (str, optional) – See the documentation of the vip_hci.preproc.frame_rotate function.
interpolation (str, optional) – See the documentation of the vip_hci.preproc.frame_rotate function.
collapse ({'median', 'mean', 'sum', 'trimmean', None}, str or None, optional) – Sets the way of collapsing the frames for producing a final image. If None then the cube of residuals is used when measuring the function of merit (instead of a single final frame).
algo (python routine, opt {pca_annulus, pca_annular, pca, custom}) – Routine to be used to model and subtract the stellar PSF. From an input cube, derotation angles, and optional arguments, it should return a post-processed frame.
delta_rot (float, optional) – If algo is set to pca_annular, delta_rot is the angular threshold used to select frames in the PCA library (see description of pca_annular).
algo_options (dict, opt) – Dictionary with additional parameters for the pca algorithm (e.g. tol, min_frames_lib, max_frames_lib). Note: arguments such as svd_mode, scaling imlib, interpolation or collapse can also be included in this dict (the latter are also kept as function arguments for compatibility with older versions of vip).
f_range (numpy.array, optional) – The range of tested flux values. If None, 20 values between 0 and 5000 are tested.
transmission (numpy array, optional) – Array with 2 columns. First column is the radial separation in pixels. Second column is the off-axis transmission (between 0 and 1) at the radial separation given in column 1.
mu_sigma (tuple of 2 floats or None, opt) – If set to None: not used, and falls back to original version of the algorithm, using fmerit. Otherwise, should be a tuple of 2 elements, containing the mean and standard deviation of pixel intensities in an annulus centered on the location of the companion, excluding the area directly adjacent to the companion.
weights (1d array, optional) – If provided, the negative fake companion fluxes will be scaled according to these weights before injection in the cube. Can reflect changes in the observing conditions throughout the sequence.
plot (boolean, optional) – If True, the figure chi2 vs. flux is displayed.
verbose (boolean) – If True, display intermediate info in the shell.
save (boolean, optional) – If True, the figure chi2 vs. flux is saved as .pdf if plot is also True
debug (bool, optional) – Whether to print details of the grid search
full_output (bool, optional) – Whether to also return the range of fluxes tested and their chi2r values.

Returns:

res (tuple) – The polar coordinates and the flux(es) of the companion.
f_range (1d numpy.array) – [full_output=True] The range of tested flux values.
chi2r (1d numpy.array) – [full_output=True] The chi2r values corresponding to tested flux values.

vip_hci.fm.negfc_speckle_noise module¶

Module with routines allowing for the estimation of the uncertainty on the parameters of an imaged companion associated to residual speckle noise.

vip_hci.fm.negfc_speckle_noise.speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo, psfn, fwhm, aperture_radius, opp_ang=False, indep_ap=False, cube_ref=None, fmerit='sum', algo_options={}, transmission=None, mu_sigma=None, wedge=None, weights=None, force_rPA=False, nproc=None, simplex_options=None, bins=None, save=False, output=None, verbose=True, full_output=True, plot=False, trim_outliers=True)[source]¶

Step-by-step procedure used to determine the speckle noise uncertainty associated to the parameters of a companion candidate.

The steps 1 to 3 need to be performed for each angle.

At the true planet radial distance and for a given angle, we inject a fake companion in our planet-free cube.
Then, using the negative fake companion method, we determine the position and flux of the fake companion thanks to a Simplex Nelder-Mead minimization.
We calculate the offset between the true values of the position and the flux of the fake companion, and those obtained from the minimization. The results will be dependent on the angular position of the fake companion.

The resulting distribution of deviations is then used to infer the 1-sigma uncertainty on each parameter by fitting a 1d-gaussian.

Parameters:

cube (3d or 4d numpy array) – The original ADI or ADI+IFS cube.
p_true (tuple or numpy array with 3 (or more) elements) – The radial separation, position angle (from x=0 axis) and flux associated to a given companion candidate for which the speckle uncertainty is to be evaluated. The planet will first be subtracted from the cube, then used for test injections. For a 4D input cube, the length of p_true should be equal to 2 (for r, theta) + the number 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[.
derot_angles (1d numpy array) – Derotation angles for ADI. Length should match input cube.
algo (python routine) – Routine to be used to model and subtract the stellar PSF. From an input cube, derotation angles, and optional arguments, it should return a post-processed frame.
psfn (2d numpy array) – 2d array with the normalized PSF template. The PSF image must be centered wrt to the array. Therefore, it is recommended to run the function metrics/normalize_psf() to generate a centered and flux-normalized PSF template.
fwhm (float) – FWHM of the PSF in pixels.
aperture_radius (float) – Radius of the apertures used for NEGFC, in terms of FWHM.
opp_ang (bool, opt) – Whether to also use opposite derotation angles to double sample size. Uses the same angle range.
indep_ap (bool, opt.) – Whether to only consider independent apertures. If yes, will supersede the range provided in angle_range, and only consider the first and last values, then fit as many non-overlapping apertures as possible. The empty cube will also be used with opposite derotation angles to double the number of independent apertures.
algo_options (dict, opt.) – Options for algo. To be provided as a dictionary. Can include ncomp (for PCA), svd_mode, collapse, imlib, interpolation, scaling, delta_rot
transmission (numpy array, optional) – Array with 2 columns. First column is the radial separation in pixels. Second column is the off-axis transmission (between 0 and 1) at the radial separation given in column 1.
mu_sigma (tuple of 2 floats, bool or None, opt) – 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 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.
force_rPA (bool, optional) – Whether to only search for optimal flux, provided (r,PA).
nproc (int or None, optional) – The number of processes to use for parallelization. If None, will be set automatically to half the number of CPUs available.
fmerit (None) – Figure of merit to use, if mu_sigma is None.
simplex_options (dict) – All the required simplex parameters, for instance {‘tol’:1e-08, ‘max_iter’:200}
bins (int or None, opt) – Number of bins for histogram of parameter deviations. If None, will be determined automatically based on number of injected fake companions.
full_output (bool, optional) – Whether to return more outputs.
output (str, optional) – The name of the output file (if save is True)
save (bool, optional) – If True, the result are pickled.
verbose (bool, optional) – If True, informations are displayed in the shell.
plot (bool, optional) – Whether to plot the gaussian fit to the distributions of parameter deviations (between retrieved and injected).
trim_outliers (bool, opt) – Whether to trim outliers when considering a Gaussian fit to the histogram of residual deviations.

Returns:

sp_unc (numpy ndarray of 3 elements) – Uncertainties on the radius, position angle and flux of the companion, respectively, associated to residual speckle noise. Only 1 element if force_rPA is set to True.
mean_dev (numpy ndarray of 3 elements) – [full_output = True] Mean deviation for each of the 3 parameters
p_simplex (numpy ndarray n_fc x 3) – [full_output = True] Parameters retrieved by the simplex for the injected fake companions; n_fc is the number of injected
offset (numpy ndarray n_fc x 3) – [full_output = True] Deviations with respect to the values used for injection of the fake companions.
chi2, nit, success (numpy ndarray of length n_fc) – [full_output = True] Outputs from the simplex function for the retrieval of the parameters of each injected companion: chi square value, number of iterations and whether the simplex converged, respectively.

vip_hci.fm.scattered_light_disk module¶

Class definition for ScatteredLightDisk, Dust_distribution and Phase_function

[AUG99]

Augereau et al. 1999
On the HR 4796 A circumstellar disk
Astronomy & Astrophysics, Volume 348, pp. 557-569
https://arxiv.org/abs/astro-ph/9906429

class vip_hci.fm.scattered_light_disk.ScatteredLightDisk(nx=200, ny=200, distance=50.0, itilt=60.0, omega=0.0, pxInArcsec=0.01225, pa=0.0, flux_max=None, density_dico={'a': 40, 'ain': 5, 'aout': -5, 'beta': 1.0, 'dens_at_r0': 1.0, 'e': 0, 'gamma': 2.0, 'ksi0': 1.0, 'name': '2PowerLaws'}, spf_dico={'g': 0.0, 'name': 'HG', 'polar': False}, xdo=0.0, ydo=0.0)[source]¶

Bases: object

Class used to generate a synthetic disc, inspired from a light version of the GRATER tool (GRenoble RAdiative TransfER) written originally in IDL [AUG99], and converted to Python by J. Milli.

check_inclination()[source]¶: Checks whether the inclination set is close to edge-on and risks to induce artefacts from the limited numerical accuracy. In such a case the inclination is changed to be less edge-on.

compute_scattered_light(halfNbSlices=25)[source]¶

Computes the scattered lignt image of the disk.

Parameters:	halfNbSlices (integer) – half number of distances along the line of sight l

print_info()[source]¶: Prints the information of the disk and image parameters

set_density_distribution(density_dico)[source]¶

Sets or updates the parameters of the density distribution

Parameters:

density_dico (dict) –

Parameters describing the dust density distribution function to be implemented. By default, it uses a two-power law dust distribution with a vertical gaussian distribution with linear flaring. This dictionary should at least contain the key “name”. For a to-power law distribution, you can set it with name:’2PowerLaws’ and with the following parameters:

a : float

Reference radius in au (default 60)

ksi0 : float

Scale height in au at the reference radius (default 1 a.u.)

gamma : float

Exponent (2=gaussian,1=exponential profile, default 2)

beta : float

Flaring index (0=no flaring, 1=linear flaring, default 1)

ain : float

Slope of the power-low distribution in the inner disk. It must be positive (default 5)

aout : float

Slope of the power-low distribution in the outer disk. It must be negative (default -5)

e : float

Eccentricity (default 0)

set_flux_max(flux_max)[source]¶

Sets the mas flux of the disk

Parameters:	flux_max (float) – the max flux of the disk in ADU

set_inclination(itilt)[source]¶

Sets the inclination of the disk.

Parameters:	itilt (float) – inclination of the disk wrt the line of sight in degrees (0 means pole-on, 90 means edge-on, default 60 degrees)

set_omega(omega)[source]¶

Sets the argument of pericenter

Parameters:	omega (float) – angle in degrees

set_pa(pa)[source]¶

Sets the disk position angle

Parameters:	pa (float) – position angle in degrees

set_phase_function(spf_dico)[source]¶

Sets the phase function of the dust

Parameters:	spf_dico (dict) – Parameters describing the scattering phase function to be implemented. Three phase functions are implemented so far: single Heyney Greenstein, double Heyney Greenstein and custum phase functions through interpolation. Read the constructor of each of those classes to know which parameters must be set in the dictionary in each case.

class vip_hci.fm.scattered_light_disk.Dust_distribution(density_dico={'a': 60, 'ain': 5, 'amin': 0.0, 'aout': -5, 'beta': 1.0, 'dens_at_r0': 1.0, 'e': 0, 'gamma': 2.0, 'ksi0': 1.0, 'name': '2PowerLaws'})[source]¶

Bases: object

This class represents the dust distribution

density_cartesian(x, y, z)[source]¶: Return the particule volume density at x,y,z, taking into account the offset of the disk.

density_cylindrical(r, costheta, z)[source]¶: Return the particule volume density at r, theta, z.

print_info(pxInAu=None)[source]¶

Utility function that displays the parameters of the radial distribution of the dust

Input:

pxInAu (optional): the pixel size in au

set_density_distribution(density_dico)[source]¶: Update the parameters of the density distribution.

class vip_hci.fm.scattered_light_disk.Phase_function(spf_dico={'g': 0.0, 'name': 'HG', 'polar': False})[source]¶

Bases: object

This class represents the scattering phase function (spf). It contains the attribute phase_function_calc that implements either a single Henyey Greenstein phase function, a double Heyney Greenstein, or any custom function (data interpolated from an input list of phi, spf(phi)).

compute_phase_function_from_cosphi(cos_phi)[source]¶

Compute the phase function at (a) specific scattering scattering angle(s) phi. The argument is not phi but cos(phi) for optimization reasons.

Parameters:	cos_phi (float or array) – cosine of the scattering angle(s) at which the scattering function must be calculated.

plot_phase_function()[source]¶: Plots the scattering phase function

print_info()[source]¶: Prints information on the type and parameters of the scattering phase function.

vip_hci.fm.utils_mcmc module¶

Module with utility functions to the MCMC (emcee) sampling for NEGFC parameter estimation.

vip_hci.fm.utils_mcmc.gelman_rubin(x)[source]¶

Determine the Gelman-Rubin \(\hat{R}\) statistical test between Markov chains.

Parameters:	x (numpy.array) – The numpy.array on which the Gelman-Rubin test is applied. This array should contain at least 2 set of data, i.e. x.shape >= (2,).
Returns:	out – The Gelman-Rubin \(\hat{R}\).
Return type:	float

Example

>>> x1 = np.random.normal(0.0,1.0,(1,100))
>>> x2 = np.random.normal(0.1,1.3,(1,100))
>>> x = np.vstack((x1,x2))
>>> gelman_rubin(x)
1.0366629898991262
>>> gelman_rubin(np.vstack((x1,x1)))
0.99

vip_hci.fm.utils_mcmc.gelman_rubin_from_chain(chain, burnin)[source]¶

Pack the MCMC chain and determine the Gelman-Rubin \(\hat{R}\) statistical test. In other words, two sub-sets are extracted from the chain (burnin parts are taken into account) and the Gelman-Rubin statistical test is performed.

Parameters:	chain (numpy.array) – The MCMC chain with the shape walkers x steps x model_parameters burnin (float in [0,1]) – The fraction of a walker which is discarded.
Returns:	out – The Gelman-Rubin \(\hat{R}\).
Return type:	float

vip_hci.fm.utils_negfc module¶

Module with post-processing related functions called from within the NEGFC algorithm.

vip_hci.fm.utils_negfc.cube_planet_free(planet_parameter, cube, angs, psfn, imlib='vip-fft', interpolation='lanczos4', transmission=None)[source]¶

Return a cube in which we have injected negative fake companion at the position/flux given by planet_parameter.

Parameters:	planet_parameter (numpy.array or list or tuple) – The (r, theta, flux) for all known companions. For a 4d cube r, theta and flux must all be 1d arrays with length equal to cube.shape[0]; i.e. planet_parameter should have shape: (n_pl,3,n_ch). cube (numpy.array) – The cube of fits images expressed as a numpy.array. angs (numpy.array) – The parallactic angle fits image expressed as a numpy.array. psfn (numpy.array) – The scaled psf expressed as a numpy.array. imlib (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function. interpolation (str, optional) – See the documentation of the `vip_hci.preproc.frame_rotate` function.
Returns:	cpf – The cube with negative companions injected at the position given in planet_parameter.
Return type:	numpy.array

vip_hci.fm.utils_negfc.find_nearest(array, value, output='index', constraint=None, n=1)[source]¶

Function to find the indices, and optionally the values, of an array’s n closest elements to a certain value. By default, only returns the index/indices.

Possible constraints: ‘ceil’, ‘floor’, None (“ceil” will return the closest element with a value greater than ‘value’, “floor” the opposite).

Parameters:

array (1d numpy array or list) – Array in which to check the closest element to value.
value (float) – Value for which the algorithm searches for the n closest elements in the array.
output (str, opt {'index','value','both' }) – Set what is returned
constraint (str, opt {None, 'ceil', 'floor'}) – If not None, will check for the closest element larger than value (ceil) or closest element smaller than value (floor).
n (int, opt) – Number of elements to be returned, sorted by proximity to the values. Default: only the closest value is returned.

Returns:

[output=’index’] (index/indices of the closest n value(s) in the array;)
[output=’value’] (the closest n value(s) in the array,)
[output=’both’] (closest value(s) and index/-ices, respectively.)

Module contents¶

Subpackge fm contains an ensemble of algorithms for forward modeling, including scattered light disk model creation, fake planet and fake disk injection, and planet position+flux estimation through the negative fake companion algorithm (NEGFC). NEGFC can work with either a simplex (Nelder-Mead) minimizer or coupled with Monte Carlo methods for posterior sampling. The latter allows the estimation of uncertainties for the photometry and position of the companion. Two possible ways of sampling the posteriors are implemented: using emcee and its Affine Invariant MCMC or nestle with either a single or multiple ellipsoid nested sampling procedure.

The main idea of the NegFC is to inject negative fake companions (candidates) with varying position and flux in the original cube and minimize a figure of merit corresponding to the residuals in an aperture at the initial estimate for the location of the companion, in the post-processed image.