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')

Create 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.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)

Inject fake companions in branches and given radial distances in an ADI cube.

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).

• rad_dists (float, list or array 1d) – Vector of radial distances of fake companions in pixels.

• plsc (float or None) – Value of the plsc in arcsec/px. Only used for printing debug output when verbose=True.

• 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.frame_inject_companion(array, array_fc, pos_y, pos_x, flux, imlib='vip-fft', interpolation='lanczos4')

Inject a fake companion in a single frame (can 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_y (float) – Y and X coordinates where the companion should be injected

• pos_x (float) – 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.fakecomp.generate_cube_copies_with_injections(array, psf_template, angle_list, plsc, n_copies=100, inrad=8, outrad=12, dist_flux=('uniform', 2, 500))

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 (float) – Inner and outer radius of the injections. The actual injection position is chosen randomly.

• 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.

positionlist of tuples(y,x)

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

distfloat

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

thetafloat, degrees

The initial angle, as passed to cube_inject_companions.

fluxfloat

The flux passed to cube_inject_companions.

Return type:

dict

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)

Normalize a PSF (2d or 3d array), to have the flux in a 1xFWHM aperture equal to one.

It also allows cropping of the array. Automatic recentering of the PSF is done internally - as long as it is already roughly centered within ~2px.

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 with respect to 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.fakedisk module

Module with fake disk injection functions.

vip_hci.fm.fakedisk.cube_inject_fakedisk(fakedisk, angle_list, psf=None, transmission=None, **rot_options)

Create 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)

Inject 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.

• angle_list (list) – Vector containing the parallactic angles.

• 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 in pixels 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 to print trace locations in arcsec).

• 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 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={})

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", spat-mean", "temp-standard",) –

  “spat-standard”}, None or str optional Pixel-wise scaling mode using sklearn.preprocessing.scale function. If set to None, the input matrix is left untouched. Otherwise:

  • temp-mean: temporal px-wise mean is subtracted.

  • spat-mean: spatial mean is subtracted.

  • temp-standard: temporal mean centering plus scaling pixel values to unit variance (temporally).

  • spat-standard: spatial mean centering plus scaling pixel values to unit variance (spatially).

  DISCLAIMER: Using temp-mean or temp-standard scaling can improve the speckle subtraction for ASDI or (A)RDI reductions. Nonetheless, this involves a sort of c-ADI preprocessing, which (i) can be dangerous for datasets with low amount of rotation (strong self-subtraction), and (ii) should probably be referred to as ARDI (i.e. not RDI stricto sensu).

• 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

[QUA15] (1,2,3,4)

Quanz et al. 2015

**Confirmation and Characterization of the Protoplanet HD 100546 b—Direct

Evidence for Gas Giant Planet Formation at 50 AU** | Astronomy & Astrophysics, Volume 807, p. 64 | `https://arxiv.org/abs/1412.5173

<https://arxiv.org/abs/1412.5173>`_

[WER17] (1,2,3,4,5,6,7,8,9,10,11,12,13)

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

<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.chain_zero_truncated(chain)

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.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)

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),

2. 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_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)

Run an affine invariant mcmc sampling algorithm to determine position and flux of a companion 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 rotation, 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', 'hessian'}, string optional) –

  If mu_sigma is not provided nor set to True, this parameter determines which figure of merit to be used:

  • sum: minimizes the sum of absolute residual intensities in the

  aperture defined with initial_state and aperture_radius. More details in [WER17].

  • stddev: minimizes the standard deviation of residual

  intensities in the aperture defined with initial_state and aperture_radius. More details in [WER17].

  • hessian: minimizes the sum of absolute values of the

  determinant of the Hessian matrix calculated for each of the 4 pixels encompassing the first guess location defined with initial_state. More details in [QUA15].

  From experience: sum is more robust for high SNR companions (but rather consider setting mu_sigma=True), while stddev tend to be more reliable in presence of strong residual speckle noise. hessian is expected to be more reliable in presence of extended signals around the companion location.

• 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 ({None, "temp-mean", spat-mean", "temp-standard",) –

  “spat-standard”}, None or str optional Pixel-wise scaling mode using sklearn.preprocessing.scale function. If set to None, the input matrix is left untouched. Otherwise:

  • temp-mean: temporal px-wise mean is subtracted.

  • spat-mean: spatial mean is subtracted.

  • temp-standard: temporal mean centering plus scaling pixel values to unit variance (temporally).

  • spat-standard: spatial mean centering plus scaling pixel values to unit variance (spatially).

  DISCLAIMER: Using temp-mean or temp-standard scaling can improve the speckle subtraction for ASDI or (A)RDI reductions. Nonetheless, this involves a sort of c-ADI preprocessing, which (i) can be dangerous for datasets with low amount of rotation (strong self-subtraction), and (ii) should probably be referred to as ARDI (i.e. not RDI stricto sensu).

• 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.show_corner_plot(chain, burnin=0.5, save=False, output_dir='', **kwargs)

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)

Display or save 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_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)

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 FWHM 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 ({None, "temp-mean", spat-mean", "temp-standard",) –

  “spat-standard”}, None or str optional Pixel-wise scaling mode using sklearn.preprocessing.scale function. If set to None, the input matrix is left untouched. Otherwise:

  • temp-mean: temporal px-wise mean is subtracted.

  • spat-mean: spatial mean is subtracted.

  • temp-standard: temporal mean centering plus scaling pixel values to unit variance (temporally).

  • spat-standard: spatial mean centering plus scaling pixel values to unit variance (spatially).

  DISCLAIMER: Using temp-mean or temp-standard scaling can improve the speckle subtraction for ASDI or (A)RDI reductions. Nonetheless, this involves a sort of c-ADI preprocessing, which (i) can be dangerous for datasets with low amount of rotation (strong self-subtraction), and (ii) should probably be referred to as ARDI (i.e. not RDI stricto sensu).

• 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)

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_negfc_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, planets_xy_coord, ncomp=1, fwhm=4, annulus_width=4, aperture_radius=1, cube_ref=None, svd_mode='lapack', scaling=None, fmerit='sum', imlib='skimage', interpolation='biquintic', collapse='median', algo=<function pca_annulus>, delta_rot=1, f_range=None, transmission=None, mu_sigma=True, wedge=None, weights=None, force_rPA=False, ndet=None, algo_options={}, simplex=True, simplex_options=None, plot=False, verbose=True, save=False)

Determine a first guess for the position and the flux of a planet using the negative fake companion techique, as explained in [WER17].

This first requires processing the cube without injecting any negative fake companion. Once planets or planet candidates are identified, their initial guess (x,y) coordinates can be provided to this function. A preliminary flux guess is then found for each planet by using the method firstguess_from_coord called within this function. Optionally, a Simplex Nelder_Mead minimization is used for a refined estimate of position and flux based on the preliminary guesses.

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.

• planets_xy_coord (array or list) – The list of (x,y) positions of the planets.

• ncomp (int or 1d numpy array of int, optional) – The number of principal components to use, if the algorithm is PCA. If the input cube is 4D, ncomp can be a list of integers, with length matching the first dimension of the cube.

• plsc (float, optional) – The platescale, in arcsec per pixel.

• fwhm (float, optional) – The FWHM 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 ({None, "temp-mean", spat-mean", "temp-standard",) –

  “spat-standard”}, None or str optional Pixel-wise scaling mode using sklearn.preprocessing.scale function. If set to None, the input matrix is left untouched. Otherwise:

  • temp-mean: temporal px-wise mean is subtracted.

  • spat-mean: spatial mean is subtracted.

  • temp-standard: temporal mean centering plus scaling pixel values to unit variance (temporally).

  • spat-standard: spatial mean centering plus scaling pixel values to unit variance (spatially).

  DISCLAIMER: Using temp-mean or temp-standard scaling can improve the speckle subtraction for ASDI or (A)RDI reductions. Nonetheless, this involves a sort of c-ADI preprocessing, which (i) can be dangerous for datasets with low amount of rotation (strong self-subtraction), and (ii) should probably be referred to as ARDI (i.e. not RDI stricto sensu).

• fmerit ({'sum', 'stddev', 'hessian'}, string optional) –

  If mu_sigma is not provided nor set to True, this parameter determines which figure of merit to be used:

  • sum: minimizes the sum of absolute residual intensities in the

  aperture defined with initial_state and aperture_radius. More details in [WER17].

  • stddev: minimizes the standard deviation of residual

  intensities in the aperture defined with initial_state and aperture_radius. More details in [WER17].

  • hessian: minimizes the sum of absolute values of the

  determinant of the Hessian matrix calculated for each of the 4 pixels encompassing the first guess location defined with initial_state. More details in [QUA15].

  From experience: sum is more robust for high SNR companions (but rather consider setting mu_sigma=True), while stddev tend to be more reliable in presence of strong residual speckle noise. hessian is expected to be more reliable in presence of extended signals around the companion location.

• 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).

• 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).

• ndet (int or None, optional) – [only used if fmerit=’hessian’] If not None, ndet should be the number of pixel(s) along x and y around the first guess position for which the determinant of the Hessian matrix is calculated. If odd, the pixel(s) around the closest integer coordinates will be considered. If even, the pixel(s) around the subpixel coordinates of the first guess location are considered. The figure of merit is the absolute sum of the determinants. If None, ndet is determined automatically to be max(1, round(fwhm/2)).

• 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=1, cube_ref=None, svd_mode='lapack', scaling=None, fmerit='sum', imlib='skimage', interpolation='biquintic', collapse='median', algo=<function pca_annulus>, delta_rot=1, algo_options={}, f_range=None, transmission=None, mu_sigma=(0, 1), weights=None, ndet=None, plot=False, verbose=True, save=False, debug=False, full_output=False)

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 using the negative fake companion technique (i.e. the reduced chi2 is calculated in the post-processed frame after subtraction of a negative fake companion).

Parameters:

• planet (numpy.array) – The (x,y) position of the planet in the 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 FWHM 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, optional) – The number of principal components, if the algorithm used is PCA.

• 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 ({None, "temp-mean", spat-mean", "temp-standard",) –

  “spat-standard”}, None or str optional Pixel-wise scaling mode using sklearn.preprocessing.scale function. If set to None, the input matrix is left untouched. Otherwise:

  • temp-mean: temporal px-wise mean is subtracted.

  • spat-mean: spatial mean is subtracted.

  • temp-standard: temporal mean centering plus scaling pixel values to unit variance (temporally).

  • spat-standard: spatial mean centering plus scaling pixel values to unit variance (spatially).

  DISCLAIMER: Using temp-mean or temp-standard scaling can improve the speckle subtraction for ASDI or (A)RDI reductions. Nonetheless, this involves a sort of c-ADI preprocessing, which (i) can be dangerous for datasets with low amount of rotation (strong self-subtraction), and (ii) should probably be referred to as ARDI (i.e. not RDI stricto sensu).

• fmerit ({'sum', 'stddev', 'hessian'}, string optional) –

  If mu_sigma is not provided nor set to True, this parameter determines which figure of merit to be used:

  • sum: minimizes the sum of absolute residual intensities in the

  aperture defined with initial_state and aperture_radius. More details in [WER17].

  • stddev: minimizes the standard deviation of residual

  intensities in the aperture defined with initial_state and aperture_radius. More details in [WER17].

  • hessian: minimizes the sum of absolute values of the

  determinant of the Hessian matrix calculated for each of the 4 pixels encompassing the first guess location defined with initial_state. More details in [QUA15].

  From experience: sum is more robust for high SNR companions (but rather consider setting mu_sigma=True), while stddev tend to be more reliable in presence of strong residual speckle noise. hessian is expected to be more reliable in presence of extended signals around the companion location.

• 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).

• 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, 30 values between 1e-1 and 1e4 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.

• ndet (int or None, optional) – [only used if fmerit=’hessian’] If not None, ndet should be the number of pixel(s) along x and y around the first guess position for which the determinant of the Hessian matrix is calculated. If odd, the pixel(s) around the closest integer coordinates will be considered. If even, the pixel(s) around the subpixel coordinates of the first guess location are considered. The figure of merit is the absolute sum of the determinants. If None, ndet is determined automatically to be max(1, round(fwhm/2)).

• 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, ndet=None, nproc=None, simplex_options=None, bins=None, save=False, output=None, verbose=True, full_output=True, plot=False, sigma_trim=None)

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.

1. At the true planet radial distance and for a given angle, we inject a fake companion in our planet-free cube.

2. Then, using the negative fake companion method, we determine the position and flux of the fake companion thanks to a Simplex Nelder-Mead minimization.

3. 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).

• ndet (int or None, optional) – [only used if fmerit=’hessian’] If not None, ndet should be the number of pixel(s) along x and y around the first guess position for which the determinant of the Hessian matrix is calculated. If odd, the pixel(s) around the closest integer coordinates will be considered. If even, the pixel(s) around the subpixel coordinates of the first guess location are considered. The figure of merit is the absolute sum of the determinants. If None, ndet is determined automatically to be max(1, round(fwhm/2)).

• 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 ({'sum', 'stddev', 'hessian'}, string optional) –

  If mu_sigma is not provided nor set to True, this parameter determines which figure of merit to be used:

  • sum: minimizes the sum of absolute residual intensities in the

  aperture defined with initial_state and aperture_radius. More details in [WER17].

  • stddev: minimizes the standard deviation of residual

  intensities in the aperture defined with initial_state and aperture_radius. More details in [WER17].

  • hessian: minimizes the sum of absolute values of the

  determinant of the Hessian matrix calculated for each of the 4 pixels encompassing the first guess location defined with initial_state. More details in [QUA15].

  From experience: sum is more robust for high SNR companions (but rather consider setting mu_sigma=True), while stddev tend to be more reliable in presence of strong residual speckle noise. hessian is expected to be more reliable in presence of extended signals around the companion location.

• 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).

• sigma_trim (float, opt) – If provided, sigma threshold used to trim out outliers before 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.negfd_fmerit module

Module with function of merit definitions for the NEGFD optimization.

vip_hci.fm.negfd_fmerit.chisquare_fd(modelParameters, cube, angs, disk_model, mask_fm, initialState, force_params=None, grid_param_list=None, fmerit='sum', mu_sigma=None, psfn=None, algo=<function pca>, algo_options={}, interp_order=-1, imlib='skimage', interpolation='biquintic', transmission=None, weights=None, debug=False, rot_options={})

Calculate the figure of merit to minimze residuals after disk subtraction.

The reduced \(\chi^2\) is defined as:: .. math:: chi^2_r = frac{1}{N-Npar}sum_{j=1}^{N} frac{(I_j-mu)^2}{sigma^2} (mu_sigma is a tuple) or: .. math:: chi^2_r = frac{1}{N-Npar}sum_{j=1}^{N} |I_j| (mu_sigma=None), where N is the number of pixels within the binary mask mask_fm, Npar the number of parameters to be fitted (4 for a 3D input cube, 3+n_ch for a 4D input cube), and \(I_j\) the j-th pixel intensity.

Parameters:

• modelParameters (tuple) – The free model parameters. E.g. (x, y, theta, scal, flux) for a 3D input cube (if force_params=None) or (x, y, theta, scal, f1, …, fN) for a 4D cube with N spectral channels (if force_params=None).

• 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.

• psfs_norm (numpy.array) – The scaled psf expressed as a numpy.array.

• 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.

• initialState (numpy.array) – All initial parameters (including fixed ones) applied to the disk model image.

• force_params (None or list/tuple of bool, optional) – If not None, list/tuple of bool corresponding to parameters to fix.

• grid_params_list (list of lists/1d nd arrays, or None) – If input disk_model is a grid of either images (for 3D input cube) or spectral cubes (for a 4D input cube), this should be provided. It should be a list of either lists or 1d nd arrays corresponding to the parameter values sampled by the input disk model grid, with their lengths matching the respective first dimensions of disk_model.

• ncomp (int or None) – The number of principal components for PCA-based algorithms.

• 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. If set to anything but None: will compute 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.

• fmerit ({'sum', 'stddev'}, string optional) – Chooses the figure of merit to be used. stddev works better for close in companions sitting on top of speckle noise.

• 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.

• 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).

• interp_order (int or tuple of int, optional, {-1,0,1}) –

  [only used if grid_params_list is not None] Interpolation mode for model interpolation. If a tuple of integers, the length should match the number of grid dimensions and will trigger a different interpolation mode for the different parameters.

  • -1: Order 1 spline interpolation in logspace for the parameter

  • 0: nearest neighbour model

  • 1: Order 1 spline interpolation

• 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) – 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.

• 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.

• debug (bool, opt) – Whether to debug and plot the post-processed frame after injection of the negative fake companion.

Returns:

out – The reduced chi squared.

Return type:

float

vip_hci.fm.negfd_interp module

Functions useful for disk model interpolation for requested parametersfalling within the provided grid.

vip_hci.fm.negfd_interp.interpolate_model(params, grid_param_list, model_grid, interp_order=-1, multispectral=False, verbose=False)

Interpolate model grid for requested parameters.

Parameters:

• params (tuple) – Set of models parameters for which the model grid has to be interpolated.

• grid_param_list (list of 1d numpy arrays/lists) – List/numpy 1d arrays with available grid of model parameters (should only contain the sampled parameters, not the models themselves).

• model_grid (numpy N-d array, optional) – Grid of model spectra for each free parameter of the given grid. For a single (resp. multi) wavelength model, the model grid should have N+2 (resp. N+3) dimensions, where N is the number of free parameters in the grid (i.e. the length of grid_param_list).

• interp_order (int or tuple of int, optional, {-1,0,1}) –

  Interpolation mode for model interpolation. If a tuple of integers, the length should match the number of grid dimensions and will trigger a different interpolation mode for the different parameters.

  • -1: Order 1 spline interpolation in logspace for the parameter

  • 0: nearest neighbour model

  • 1: Order 1 spline interpolation

• multispectral (bool, optional) – Whether the model grid is computed for various wavelenghts - e.g. for IFS data. In this case, the wavelength dimension should be the third to last in the input model_grid.

• verbose (bool, optional) – Whether to print more information during the interpolation.

Returns:

model – Interpolated model for input parameters. First column corresponds to wavelengths, and the second contains model values.

Return type:

2d or 3d numpy array

vip_hci.fm.negfd_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.negfd_simplex.firstguess_fd(cube, angs, disk_model, mask_fm, ini_xy=(0, 0), ini_theta=0, ini_scal=1.0, ini_f=None, grid_params_list=None, grid_params_labels=None, fmerit='sum', mu_sigma=None, f_range=None, psfn=None, algo=<function pca>, algo_options={}, interp_order=-1, imlib='skimage', interpolation='biquintic', simplex=True, simplex_options=None, transmission=None, weights=None, force_params=None, plot=False, verbose=True, save=False, full_output=False, rot_options={})

Determine a first guess for the shifts (x,y), rotation, spatial scaling and flux scaling of a disk model image.

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.

• disk_model (numpy ndarray) – The disk image(s) to be injected, should be a 2d ndarray (for a 3D input cube), 3d numpy array (for a 4D spectral+ADI input cube), or a higher dimensionality for an input grid of disk models provided (the number of additional dimensions is inferred automatically depending on input cube. For a spectral+ADI input cube, the disk images should correspond to different wavelengths with its zeroth shape matching cube’s.

• mask_fm (2d numpy ndarray) – Binary mask on which to calculate the figure of merit in the processed image. Residuals will be minimized where mask values are 1.

• ini_xy (tuple or numpy array of 2 elements) – Initial estimate of the x,y shift to be applied to the disk model image.

• ini_theta (float) – Initial estimate of the rotation angle to be applied to the disk model image (after shift).

• ini_scal (float) – Initial estimate of the spatial scaling factor to be applied to the disk model image (after shift and rotation).

• ini_f (float, 1d ndarray or None) – Initial estimate of the spatial scaling factor to be applied to the disk model image (after shift and rotation). If None, a grid on f_range is used to get a first estimate of this parameter. Else, the provided estimate is directly used for a simplex minimzation.

• grid_params_list (list of lists/1d nd arrays, or None, optional) – If input disk_model is a grid of either images (for 3D input cube) or spectral cubes (for a 4D input cube), this should be provided. It should be a list of either lists or 1d nd arrays corresponding to the parameter values sampled by the input disk model grid, with their lengths matching the respective first dimensions of disk_model.

• grid_params_labels (list/tuple of str, or None, optional) – If grid_params_list is provided, these are the name of the physical parameters that are probed in the grid. This is only used for the purpose of printing/writing results. If left to None, they will be named ‘param_i’ with i ranging from 1 to the length of grid_params_list.

• fmerit ({'sum', 'stddev'}, string optional) – Figure of merit to be used, if mu_sigma is set to None. ‘sum’ will find optimal parameters that minimize the absolute intensity residuals in mask_fm. ‘stddev’ will minimize the standard deviation of pixel intensity residuals in mask_fm.

• mu_sigma (tuple of 2 floats or None, opt) – If set to None: not used, and falls back to using fmerit. Otherwise, should be a tuple of 2 elements, containing either the mean and standard deviation of pixel intensities in the image (i.e. floats), or individual maps of the expected values and uncertainties for each pixel in the image (i.e. np.ndarrays). In the latter case, the arrays must be 2D (for a 3D input cube) or 3D (for a 4D input cube). It can also be a mix of both, e.g. a float and a 2D array.

• f_range (numpy.array or None, optional) – The range of tested flux scaling values. If None and ini_f is also None, a grid of 30 values between 1e-1 and 1e4 are tested, following a geometric progression.

• psfn (2d or 3d numpy ndarray) – The normalized psf expressed as a numpy ndarray. Can be 3d for a 4d (spectral+ADI) input cube. This would only be used to convolve disk_img. Leave to None if the disk_img is already convolved.

• algo (python routine, opt {pca, pca_annulus, pca_annular, 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.

• algo_options (dict, opt) – Dictionary with additional parameters for the algorithm (e.g. ncomp, fwhm, asize, delta_rot, tol, min_frames_lib, max_frames_lib, cube_ref, svd_mode=, scaling, imlib, interpolation, collapse, if relevant). By default, imlib is set to ‘skimage’ which is faster than ‘vip-fft’. If opencv is installed, it is recommended to set imlib to ‘opencv’ and interpolation to ‘lanczos4’.

• interp_order (int or tuple of int, optional, {-1,0,1}) –

  [only used if grid_params_list is not None] Interpolation mode for model interpolation. If a tuple of integers, the length should match the number of grid dimensions and will trigger a different interpolation mode for the different parameters.

  • -1: Order 1 spline interpolation in logspace for the parameter

  • 0: nearest neighbour model

  • 1: Order 1 spline interpolation

• imlib (str, optional) – See the documentation of the vip_hci.preproc.frame_rotate function. By default, imlib is set to ‘skimage’ for NEGFC as it is faster than ‘vip-fft’. If opencv is installed, it is recommended to set imlib to ‘opencv’ and interpolation to ‘lanczos4’. Takes precedence over value provided in algo_options.

• interpolation (str, optional) – See the documentation of the vip_hci.preproc.frame_rotate function. If opencv is installed, it is recommended to set imlib to ‘opencv’ and interpolation to ‘lanczos4’. Takes precedence over value provided in algo_options.

• simplex (bool, optional) – If True, the Nelder-Mead minimization is performed either after the flux grid search, or using an initial ini_f estimate (if provided).

• simplex_options (dict, optional) – The scipy.optimize.minimize options.

• 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.

• 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. Length should match the temporal axis of the cube.

• force_params (None or list/tuple of bool, optional) – If not None, list/tuple of bool corresponding to parameters to fix. For a 3D input cube, the length of the list/tuple should be ngrid+5, where ngrid correspond to the number of dimensions in the provided model grid (0 if a single model image is provided), and the 5 extra dimensions correspond to shifts (x,y), rotation, spatial scaling and flux scaling, respectively. For a 4D input cube, the length of the list/tuple should be ngrid+4+nch, where nch is the number of spectral channels.

• 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.

• full_output (bool, optional) – Whether to also return the chi2r, apart from the optimal parameters.

Returns:

• x_0 (float) – Optimal x shift of the model disk image

• y_0 (float) – Optimal y shift of the model disk image.

• theta_0 (float) – Optimal rotation angle (with respect to x axis) for the model disk image.

• scal_0 (float) – Optimal spatial scaling factor for the model disk image.

• f_0 (float) – Optimal flux scaling factor for the model disk image.

• chi2 (float) – [full_output=True] Corresponding chi2r

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.negfd_simplex.firstguess_fd_from_coord(disk_xy, disk_theta, disk_scal, cube, angs, disk_img, mask_fm, fmerit='sum', mu_sigma=None, f_range=None, psfn=None, algo=<function pca>, algo_options={}, interp_order=-1, imlib='skimage', interpolation='biquintic', transmission=None, weights=None, plot=False, verbose=True, save=False, debug=False, full_output=False, rot_options={})

Determine a first guess for the flux scaling of the disk image for a given xy shift and rotation, by doing a simple grid search evaluating the reduced chi2.

Parameters:

• disk_xy (numpy.array) – The (x,y) shift for the disk image in the processed image.

• disk_theta (float) – The rotation angle to be applied to the disk image (after shift) in the processed image.

• disk_scal (float) – The spatial scaling factor to apply on the disk image (after shift and rotation) in the processed image.

• 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.

• disk_img (2d or 3d numpy ndarray) – The disk image to be injected, as a 2d ndarray (for a 3D input cube) or a 3d numpy array (for a 4D spectral+ADI input cube). In the latter case, the images should correspond to different wavelengths, and the zeroth shape of disk_model and cube should match.

• mask_fm (2d numpy ndarray) – Binary mask on which to calculate the figure of merit in the processed image. Residuals will be minimized where mask values are 1.

• fmerit ({'sum', 'stddev'}, string optional) – Figure of merit to be used, if mu_sigma is set to None. ‘sum’ will find optimal parameters that minimize the absolute intensity residuals in mask_fm. ‘stddev’ will minimize the standard deviation of pixel intensity residuals in mask_fm.

• 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 encompassing most of the disk signal.

• f_range (numpy.array, optional) – The range of tested flux scaling values. If None, 30 values between 1e-1 and 1e4 are tested, following a geometric progression.

• psfn (2d or 3d numpy ndarray) – The normalized psf expressed as a numpy ndarray. Can be 3d for a 4d (spectral+ADI) input cube. This would only be used to convolve disk_img. Leave to None if the disk_img is already convolved.

• 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.

• algo_options (dict, opt) – Dictionary with additional parameters for the algorithm (e.g. ncomp, fwhm, asize, delta_rot, tol, min_frames_lib, max_frames_lib, cube_ref, svd_mode=, scaling, imlib, interpolation, collapse, if relevant).

• interp_order (int or tuple of int, optional, {-1,0,1}) –

  [only used if grid_params_list is not None] Interpolation mode for model interpolation. If a tuple of integers, the length should match the number of grid dimensions and will trigger a different interpolation mode for the different parameters.

  • -1: Order 1 spline interpolation in logspace for the parameter

  • 0: nearest neighbour model

  • 1: Order 1 spline interpolation

• imlib (str, optional) – See the documentation of the vip_hci.preproc.frame_rotate function. By default, imlib is set to ‘skimage’ for NEGFC as it is faster than ‘vip-fft’. If opencv is installed, it is recommended to set imlib to ‘opencv’ and interpolation to ‘lanczos4’. Takes precedence over value provided in algo_options.

• interpolation (str, optional) – See the documentation of the vip_hci.preproc.frame_rotate function. If opencv is installed, it is recommended to set imlib to ‘opencv’ and interpolation to ‘lanczos4’. Takes precedence over value provided in algo_options.

• 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.

• weights (1d array, optional) – If provided, the disk image flux factors 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 x shift, y shift, rotation angle and flux scaling of the disk image.

• 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.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.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'})

Bases: object

This class represents the dust distribution

density_cartesian(x, y, z)

Return the particule volume density at x,y,z, taking into account the offset of the disk.

density_cylindrical(r, costheta, z)

Return the particule volume density at r, theta, z.

print_info(pxInAu=None)

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)

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})

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)

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()

Plots the scattering phase function

print_info()

Prints information on the type and parameters of the scattering phase function.

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)

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()

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)

Computes the scattered light image of the disk.

Parameters:

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

print_info()

Prints the information of the disk and image parameters

set_density_distribution(density_dico)

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:

• afloat

  Reference radius in au (default 60)

• ksi0float

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

• gammafloat

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

• betafloat

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

• ainfloat

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

• aoutfloat

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

• efloat

  Eccentricity (default 0)

set_flux_max(flux_max)

Sets the mas flux of the disk

Parameters:

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

set_inclination(itilt)

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)

Sets the argument of pericenter

Parameters:

omega (float) – angle in degrees

set_pa(pa)

Sets the disk position angle

Parameters:

pa (float) – position angle in degrees

set_phase_function(spf_dico)

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.

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)

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)

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)

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 ndarray) – The cube of fits images expressed as a numpy.array.

• angs (numpy ndarray) – The parallactic angle fits image expressed as a numpy.array.

• psfn (2d or 3d numpy ndarray) – The normalized psf expressed as a numpy ndarray. Can be 3d for a 4d (spectral+ADI) input cube.

• 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.

• 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).

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)

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', 'ceil=', 'floor='}) – If not None, will check for the closest element larger (or equal) than value if set to ‘ceil’ (‘ceil=’), or closest element smaller (or equal) than value if set to ‘floor’ (‘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.)

vip_hci.fm.utils_negfd module

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

vip_hci.fm.utils_negfd.cube_disk_free(disk_parameter, cube, derot_angs, disk_img, psfn=None, imlib='vip-fft', interpolation='lanczos4', imlib_sh='vip-fft', interpolation_sh='lanczos4', imlib_sc='vip-fft', interpolation_sc='lanczos4', transmission=None, weights=None, **rot_options)

Return a cube in which we have injected a negative fake disk model image with xy shifts, spatial scaling, rotation and flux scaling given by disk_parameter.

Parameters:

• disk_parameter (numpy.array or list or tuple) – The (delta_x, delta_y, theta, scal, flux) for the model disk image. For a 4d cube, disk_parameter should be either (i) a numpy array with shape (5,n_ch), (ii) a 5-element list/tuple where delta_x, delta_y, theta and scal can be floats or 1d-arrays, but flux mandatorily as a 1d array with length equal to cube.shape[0].

• 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.

• disk_img (2d or 3d numpy ndarray) – The disk image to be injected, as a 2d ndarray (for either 3D or 4D input cubes) or a 3d numpy array (for a 4D spectral+ADI input cube). In the latter case, the images should correspond to different wavelengths, and the zeroth shape of disk_model and cube should match.

• psfn (2d or 3d numpy ndarray) – The normalized psf expressed as a numpy ndarray. Can be 3d for a 4d (spectral+ADI) input cube. This would only be used to convolve disk_img. Leave to None if the disk_img is already convolved.

• 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.

• imlib_sh (str, optional) – See the documentation of the vip_hci.preproc.frame_shift function.

• interpolation_sh (str, optional) – See the documentation of the vip_hci.preproc.frame_shift function.

• imlib_sc (str, optional) – See the documentation of the vip_hci.preproc.frame_rescaling function.

• interpolation_sc (str, optional) – See the documentation of the vip_hci.preproc.frame_rescaling 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).

• 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:

cdf – The cube with negative disk model injected at the position given in disk_parameter.

Return type:

numpy.array

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.