5A. ADI forward modeling of point sources
Authors: Valentin Christiaens and Carlos Alberto Gomez Gonzalez
Suitable for VIP v1.0.3 onwards
Last update: 2024/03/25
Table of contents
This tutorial shows:
how to generate and inject fake companions in a cube;
how to estimate the astrometry and photometry of a directly imaged companion, and associated uncertainties.
A number of routines in the fm subpackage have been implemented for compatibility with multiprocessing, in order to optimally harness the power of machines equipped with multiple CPUs. Any function where the nproc parameter is available in its call (or which internally calls a psfsub function) can be run in multi-processing, with the value of nproc setting the requested number of CPUs to use. Instead of an integer, one can set nproc=None to use half of all available CPUs. For optimal results in multiprocessing, set the following environment parameters BEFORE launching your Jupyter notebook:
export MKL_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1
export OMP_NUM_THREADS=1
In the case of PCA, singular value decomposition can also be done on GPU by setting svd_mode to an appropriate value (see Sec. 3.5.1. and docstrings of vip_hci.psfsub.pca for details).
Let’s first import a couple of external packages needed in this tutorial:
%matplotlib inline
from hciplot import plot_frames, plot_cubes
from matplotlib.pyplot import *
from matplotlib import pyplot as plt
import numpy as np
from packaging import version
In the following box we check that your version of VIP passes the requirements to run this notebook:
import vip_hci as vip
vvip = vip.__version__
print("VIP version: ", vvip)
if version.parse(vvip) <= version.parse("1.0.3"):
msg = "Please upgrade your version of VIP"
msg+= "It should be strictly above 1.0.3 to run this notebook."
raise ValueError(msg)
VIP version: 1.6.2
5.1. Loading ADI data
In the ‘dataset’ folder of the VIP_extras repository you can find a toy ADI (Angular Differential Imaging) cube and a NACO point spread function (PSF) to demonstrate the capabilities of VIP. This is an L’-band VLT/NACO dataset of beta Pictoris published in Absil et al. (2013) obtained using the Annular Groove Phase Mask (AGPM) Vortex coronagraph. The sequence has been heavily sub-sampled temporarily to make it smaller. The frames were also cropped to the central 101x101 area. In case you want to plug-in your cube just change the path of the following cells.
More info on this dataset, and on opening and visualizing fits files with VIP in general, is available in Tutorial 1. Quick start.
Let’s load the data:
from vip_hci.fits import open_fits
from astropy.utils.data import download_file
url_d = "https://github.com/vortex-exoplanet/VIP_extras/raw/master/datasets"
f1 = download_file("{}/naco_betapic_cube_cen.fits".format(url_d), cache=True)
f2 = download_file("{}/naco_betapic_psf.fits".format(url_d), cache=True)
f3 = download_file("{}/naco_betapic_derot_angles.fits".format(url_d), cache=True)
# alternatively, for local files simply provide their full or relative path. E.g.:
#f1 = '../datasets/naco_betapic_cube_cen.fits'
#f2 = '../datasets/naco_betapic_psf.fits'
#f3 = '../datasets/naco_betapic_derot_angles.fits'
cube = open_fits(f1)
psf = open_fits(f2)
angs = open_fits(f3)
FITS HDU-0 data successfully loaded. Data shape: (61, 101, 101)
FITS HDU-0 data successfully loaded. Data shape: (39, 39)
FITS HDU-0 data successfully loaded. Data shape: (61,)
Let’s fit the PSF with a 2D Gaussian to infer the FWHM, the flux in a 1-FWHM size aperture, and get a flux-normalized PSF:
%matplotlib inline
from vip_hci.fm import normalize_psf
psfn, flux, fwhm_naco = normalize_psf(psf, size=19, debug=True, full_output=True)
FWHM_y = 4.926059872957138
FWHM_x = 4.67577889500593
centroid y = 9.010992107833063
centroid x = 9.01917912265807
centroid y subim = 9.010992107833063
centroid x subim = 9.01917912265807
amplitude = 0.10032285220380605
theta = -38.446187060503924
Mean FWHM: 4.801
Flux in 1xFWHM aperture: 1.307
print(fwhm_naco)
4.800919383981534
Let’s visualize the flux-normalized PSF:
plot_frames(psfn, grid=True, size_factor=4)
Let’s finally define the pixel scale for NACO (L’ band), which we get from a dictionary stored in the config subpackage:
from vip_hci.config import VLT_NACO
pxscale_naco = VLT_NACO['plsc']
print(pxscale_naco, "arcsec/px")
0.02719 arcsec/px
5.2. Generating and injecting synthetic planets
We first select an image library imlib for image operations (shifts, rotations) and associated interpolation order. ‘vip-fft’ is more accurate, but ‘skimage’ is faster, and ‘opencv’ even faster - see Tutorial 7 for more details.
imlib_rot = 'skimage' # If you have opencv installed, feel free to set this parameter to" 'opencv'
interpolation= 'biquintic' # If you have opencv installed, feel free to set this parameter to 'lanczos4'
The cube_inject_companions function in the fm module (VIP versions >= 1.1.0) makes the injection of fake companions at arbitrary fluxes and locations very easy. The normalized non-coronagraphic PSF should be provided for the injection. If the user does not have access to an observed PSF, the create_synth_psf from the var module can be used to create synthetic ones (based on 2D Gaussian, Moffat or Airy models).
Some procedures, e.g. the negative fake companion technique and the contrast curve generation, heavily rely on the injection of fake companions. The coordinates for the injection should be provided in the derotated image, while the actual injection occurs in the images of the input cube, i.e. in the rotated field.
rad_fc = 30.5
theta_fc = 240
flux_fc = 400.
gt = [rad_fc, theta_fc, flux_fc]
The convention in VIP is to measure angles from positive x axis (i.e. as trigonometric angles), as in most other Python packages. This implies adding 90º to position angles measured following the common astronomic convention (East from North).
from vip_hci.fm import cube_inject_companions
cubefc = cube_inject_companions(cube, psf_template=psfn, angle_list=angs, flevel=flux_fc, plsc=pxscale_naco,
rad_dists=[rad_fc], theta=theta_fc, n_branches=1,
imlib=imlib_rot, interpolation=interpolation)
Let’s set the corresponding cartesian coordinates:
from vip_hci.var import frame_center
cy, cx = frame_center(cube[0])
x_fc = cx + rad_fc*np.cos(np.deg2rad(theta_fc))
y_fc = cy + rad_fc*np.sin(np.deg2rad(theta_fc))
xy_test = (x_fc, y_fc)
print('({:.1f}, {:.1f})'.format(xy_test[0],xy_test[1]))
(34.7, 23.6)
Let’s double-check the fake companion was injected at the right location, by post-processing the cube and checking the final image. Let’s use PCA, and infer the optimal \(n_{\rm pc}\) while we are at it - this will be useful for the next section.
from vip_hci.psfsub import pca_grid
res_ann_opt = pca_grid(cubefc, angs, fwhm=fwhm_naco, range_pcs=(1,41,1), source_xy=xy_test, mode='annular',
annulus_width=4*fwhm_naco, imlib=imlib_rot, interpolation=interpolation,
full_output=True, plot=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:07:10
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time: 0:00:00.021671
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 41
Optimal number of PCs = 7, for S/N=11.105
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 34.7, 23.6
Flux in a centered 1xFWHM circular aperture = 146.645
Central pixel S/N = 13.442
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 11.105
Max S/N (shifting the aperture center) = 15.556
stddev S/N (shifting the aperture center) = 2.665
The grid search looking for the optimal number of principal components (npc) found that 10 principal components maximizes the S/N ratio of the injected fake companion.
_, final_ann_opt, _, opt_npc_ann = res_ann_opt
plot_frames(final_ann_opt, label='PCA-ADI annulus (opt. npc={:.0f})'.format(opt_npc_ann),
dpi=100, vmin=-5, colorbar=True, circle=xy_test)
We can see that the fake companion was indeed injected at the requested location.
5.3. Flux and position estimation with NEGFC
When a companion candidate is detected, the next step is to characterize it, i.e. infer its exact position (astrometry) and flux (photometry).
Question 5.1: Why would a simple 2D Gaussian fit (as performed e.g. for the stellar PSF in Section 5.1) be inappropriate to extract the astrometry and photometry of a candidate companion?
VIP implements the Negative fake companion (NEGFC) technique for robust extraction of the position and flux of detected point-like sources. The technique can be summarized as follow (see full description in Wertz et al. 2017):
Estimate the position and flux of the planet, from either the visual inspection of reduced images or a previous estimator (see ABC below).
Scale (in flux) and shift the normalized off-axis PSF to remove the estimate from the input data cube.
Process the cube with PCA in a single annulus encompassing the point source.
Measure residuals in an aperture centered on the approximate location of the companion candidate.
Iterate on the position and flux of the injected negative PSF (steps 2-4), until the absolute residuals in the aperture are minimized (i.e. the injected negative companion flux and the position match exactly that of the true companion).
Iterations between steps 2-4 can be performed in one of 3 ways - sorted in increasing computation time and accuracy:
A) a grid search on the flux only, provided a fixed estimate of the position (implemented in the firstguess function);
B) a Nelder-Mead simplex algorithm (firstguess function with the simplex=True option);
C) an MCMC sampler, which has the advantage to also yield uncertainties on each of the parameters of the point source (mcmc_negfc_sampling function).
Different figures of merit can be used for minimization of the residuals (Wertz et al. 2017; Christiaens et al. 2021): $\(\chi^2 = \sum_j^N \frac{(I_j-\mu)^2}{\sigma^2}~~({\rm default}), ~~~~~~~~\chi^2 = \sum_j^N |I_j|, ~~~~~~~~\chi^2 = N {\rm std}(I_j).\)$
where \(j \in {1,...,N}\), \(N\) is the total number of pixels contained in the circular aperture around the companion candidate, \(\mu\) and \(\sigma\) are the mean and standard deviation (\(N_{\rm resel}\) degrees of freedom) of pixel intensities in a truncated annulus at the radius of the companion candidate, but avoiding the azimuthal region encompassing the negative side lobes.
5.3.1. Nelder-Mead based optimization
With the function firstguess, we can obtain a first estimation of the flux and position by running A) a naive grid minimization (grid of values for the flux through parameter f_range), and B) a Nelder-mead based minimization (if the parameter simplex is set to True). The latter is done based on the preliminary guess of the grid minimization. The maximum number of iterations and error can be set with the parameter simplex_options as a dicitionary (see scipy.minimize function for the Nelder-Mead options).
Fisrt we define the position of the sources by examining a flux frame or S/N map. planets_xy_coord takes a list or array of X,Y pairs like ((x1,y1),(x2,y2)…(x_n,y_n)). Let’s take the coordinates of the previously injected companion.
Let’s test the algorithm with different values for the # of PCs: 5 and 25.
from vip_hci.fm import firstguess
r_lo, theta_lo, f_lo = firstguess(cubefc, angs, psfn, ncomp=5, planets_xy_coord=[xy_test],
fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
aperture_radius=2, simplex=True, imlib=imlib_rot,
interpolation=interpolation, plot=True, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-10-28 16:07:28
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux | chi2r
1/30 0.100 1.619
2/30 0.149 1.618
3/30 0.221 1.618
4/30 0.329 1.617
5/30 0.489 1.616
6/30 0.728 1.614
7/30 1.083 1.612
8/30 1.610 1.608
9/30 2.395 1.602
10/30 3.562 1.594
11/30 5.298 1.582
12/30 7.880 1.563
13/30 11.721 1.536
14/30 17.433 1.495
15/30 25.929 1.436
16/30 38.566 1.352
17/30 57.362 1.233
18/30 85.317 1.053
19/30 126.896 0.817
20/30 188.739 0.521
21/30 280.722 0.216
22/30 417.532 0.061
23/30 621.017 0.545
24/30 923.671 2.828
25/30 1373.824 8.789
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[15], line 2
1 from vip_hci.fm import firstguess
----> 2 r_lo, theta_lo, f_lo = firstguess(cubefc, angs, psfn, ncomp=5, planets_xy_coord=[xy_test],
3 fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
4 aperture_radius=2, simplex=True, imlib=imlib_rot,
5 interpolation=interpolation, plot=True, verbose=True)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/fm/negfc_simplex.py:763, in firstguess(cube, angs, psfn, planets_xy_coord, ncomp, fwhm, annulus_width, aperture_radius, cube_ref, svd_mode, scaling, fmerit, imlib, interpolation, collapse, algo, delta_rot, f_range, transmission, mu_sigma, wedge, weights, force_rPA, ndet, algo_options, simplex, simplex_options, plot, verbose, save)
759 if simplex_options is None:
760 simplex_options = {'xatol': 1e-6, 'fatol': 1e-6, 'maxiter': 800,
761 'maxfev': 2000}
--> 763 res = firstguess_simplex(res_init, cube, angs, psfn, ncomp, fwhm,
764 annulus_width, aperture_radius,
765 cube_ref=cube_ref, svd_mode=svd_mode,
766 scaling=scaling, fmerit=fmerit,
767 imlib=imlib, interpolation=interpolation,
768 collapse=collapse, algo=algo,
769 delta_rot=delta_rot,
770 algo_options=algo_options,
771 transmission=transmission,
772 mu_sigma=mu_sigma, weights=weights,
773 force_rPA=force_rPA, ndet=ndet,
774 options=simplex_options, verbose=False)
775 if force_rPA:
776 r_0[i_planet], theta_0[i_planet] = (r_pre, theta_pre)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/fm/negfc_simplex.py:481, in firstguess_simplex(p, cube, angs, psfn, ncomp, fwhm, annulus_width, aperture_radius, cube_ref, svd_mode, scaling, fmerit, imlib, interpolation, collapse, algo, delta_rot, algo_options, p_ini, transmission, mu_sigma, weights, force_rPA, ndet, options, verbose, **kwargs)
479 else:
480 p_t = p
--> 481 solu = minimize(chisquare, p_t, args=(cube, angs, psfn, fwhm, annulus_width,
482 aperture_radius, p_ini, ncomp,
483 cube_ref, svd_mode, scaling, fmerit,
484 collapse, algo, delta_rot, imlib,
485 interpolation, algo_options,
486 transmission, mu_sigma, weights,
487 force_rPA, ndet),
488 method='Nelder-Mead', options=options, **kwargs)
490 if verbose:
491 print(solu)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/scipy/optimize/_minimize.py:719, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
716 callback = _wrap_callback(callback, meth)
718 if meth == 'nelder-mead':
--> 719 res = _minimize_neldermead(fun, x0, args, callback, bounds=bounds,
720 **options)
721 elif meth == 'powell':
722 res = _minimize_powell(fun, x0, args, callback, bounds, **options)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/scipy/optimize/_optimize.py:842, in _minimize_neldermead(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, bounds, **unknown_options)
840 if bounds is not None:
841 xr = np.clip(xr, lower_bound, upper_bound)
--> 842 fxr = func(xr)
843 doshrink = 0
845 if fxr < fsim[0]:
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/scipy/optimize/_optimize.py:526, in _wrap_scalar_function_maxfun_validation.<locals>.function_wrapper(x, *wrapper_args)
524 ncalls[0] += 1
525 # A copy of x is sent to the user function (gh13740)
--> 526 fx = function(np.copy(x), *(wrapper_args + args))
527 # Ideally, we'd like to a have a true scalar returned from f(x). For
528 # backwards-compatibility, also allow np.array([1.3]),
529 # np.array([[1.3]]) etc.
530 if not np.isscalar(fx):
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/fm/negfc_fmerit.py:244, in chisquare(modelParameters, cube, angs, psfs_norm, fwhm, annulus_width, aperture_radius, initialState, ncomp, cube_ref, svd_mode, scaling, fmerit, collapse, algo, delta_rot, imlib, interpolation, algo_options, transmission, mu_sigma, weights, force_rPA, ndet, debug)
241 # Perform PCA and extract the zone of interest
242 full_output = (debug and collapse) or (fmerit == "hessian")
--> 244 res = get_values_optimize(
245 cube_negfc,
246 angs,
247 ncomp,
248 annulus_width,
249 aperture_radius,
250 fwhm,
251 initialState[0],
252 initialState[1],
253 cube_ref=cube_ref,
254 svd_mode=svd_mode,
255 scaling=scaling,
256 algo=algo,
257 delta_rot=delta_rot,
258 collapse=collapse,
259 algo_options=algo_options,
260 weights=norm_weights,
261 imlib=imlib_rot,
262 interpolation=interpolation,
263 full_output=full_output,
264 )
266 if full_output:
267 values, frpca = res
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/fm/negfc_fmerit.py:461, in get_values_optimize(cube, angs, ncomp, annulus_width, aperture_radius, fwhm, r_guess, theta_guess, cube_ref, svd_mode, scaling, algo, delta_rot, imlib, interpolation, collapse, algo_options, weights, full_output)
458 mask_rdi = algo_options.get("mask_rdi", None)
460 if algo == pca_annulus:
--> 461 res = pca_annulus(
462 cube,
463 angs,
464 ncomp,
465 annulus_width,
466 r_guess,
467 cube_ref,
468 svd_mode,
469 scaling,
470 imlib=imlib,
471 interpolation=interpolation,
472 collapse=collapse,
473 collapse_ifs=collapse_ifs,
474 weights=weights,
475 nproc=nproc
476 )
478 elif algo == pca_annular or algo == nmf_annular:
479 tol = algo_options.get("tol", 1e-1)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/psfsub/utils_pca.py:722, in pca_annulus(cube, angs, ncomp, annulus_width, r_guess, cube_ref, svd_mode, scaling, collapse, weights, collapse_ifs, **rot_options)
719 return cube_zeros
721 if cube.ndim == 3:
--> 722 return _pca_annulus_3d(cube, angs, ncomp, annulus_width, r_guess,
723 cube_ref, svd_mode, scaling, collapse, weights,
724 **rot_options)
725 elif cube.ndim == 4:
726 nch = cube.shape[0]
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/psfsub/utils_pca.py:706, in pca_annulus.<locals>._pca_annulus_3d(cube, angs, ncomp, annulus_width, r_guess, cube_ref, svd_mode, scaling, collapse, weights, **rot_options)
703 cube_zeros[:, yy, xx] = residuals
705 if angs is not None:
--> 706 cube_res_der = cube_derotate(cube_zeros, angs, **rot_options)
707 if collapse is not None:
708 pca_frame = cube_collapse(cube_res_der, mode=collapse,
709 w=weights)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/preproc/derotation.py:369, in cube_derotate(array, angle_list, imlib, interpolation, cxy, nproc, border_mode, mask_val, edge_blend, interp_zeros, ker)
367 array_der = np.zeros_like(array)
368 for i in range(n_frames):
--> 369 array_der[i] = frame_rotate(array[i], -angle_list[i], imlib=imlib,
370 interpolation=interpolation, cxy=cxy,
371 border_mode=border_mode,
372 mask_val=mask_val,
373 edge_blend=edge_blend,
374 interp_zeros=interp_zeros, ker=ker)
375 elif nproc > 1:
377 res = pool_map(nproc, _frame_rotate_mp, iterable(array),
378 iterable(-angle_list), imlib, interpolation, cxy,
379 border_mode, mask_val, edge_blend, interp_zeros, ker)
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/vip_hci/preproc/derotation.py:262, in frame_rotate(array, angle, imlib, interpolation, cxy, border_mode, mask_val, edge_blend, interp_zeros, ker)
259 norm = False
260 im_temp = array_prep.copy()
--> 262 array_out = rotate(im_temp, angle, order=order, center=(cx, cy),
263 cval=0, mode=border_mode)
265 if norm:
266 array_out *= max_val
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/skimage/transform/_warps.py:444, in rotate(image, angle, resize, center, order, mode, cval, clip, preserve_range)
441 # Make sure the transform is exactly affine, to ensure fast warping.
442 tform.params[2] = (0, 0, 1)
--> 444 return warp(
445 image,
446 tform,
447 output_shape=output_shape,
448 order=order,
449 mode=mode,
450 cval=cval,
451 clip=clip,
452 preserve_range=preserve_range,
453 )
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/skimage/transform/_warps.py:1036, in warp(image, inverse_map, map_args, output_shape, order, mode, cval, clip, preserve_range)
1033 prefilter = order > 1
1035 ndi_mode = _to_ndimage_mode(mode)
-> 1036 warped = ndi.map_coordinates(
1037 image, coords, prefilter=prefilter, mode=ndi_mode, order=order, cval=cval
1038 )
1040 _clip_warp_output(image, warped, mode, cval, clip)
1042 return warped
File ~/checkouts/readthedocs.org/user_builds/vip/envs/latest/lib/python3.10/site-packages/scipy/ndimage/_interpolation.py:471, in map_coordinates(input, coordinates, output, order, mode, cval, prefilter)
469 filtered = input
470 mode = _ni_support._extend_mode_to_code(mode)
--> 471 _nd_image.geometric_transform(filtered, None, coordinates, None, None,
472 output, order, mode, cval, npad, None, None)
473 return output
KeyboardInterrupt:
r_hi, theta_hi, f_hi = firstguess(cubefc, angs, psfn, ncomp=25, planets_xy_coord=[xy_test],
fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
aperture_radius=2, imlib=imlib_rot, interpolation=interpolation,
simplex=True, plot=True, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 14:55:05
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux | chi2r
1/30 0.100 0.216
2/30 0.149 0.216
3/30 0.221 0.216
4/30 0.329 0.216
5/30 0.489 0.216
6/30 0.728 0.216
7/30 1.083 0.216
8/30 1.610 0.215
9/30 2.395 0.215
10/30 3.562 0.215
11/30 5.298 0.214
12/30 7.880 0.213
13/30 11.721 0.211
14/30 17.433 0.209
15/30 25.929 0.210
16/30 38.566 0.206
17/30 57.362 0.197
18/30 85.317 0.176
19/30 126.896 0.146
20/30 188.739 0.120
21/30 280.722 0.077
22/30 417.532 0.055
23/30 621.017 0.176
24/30 923.671 0.360
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 130, nfev: 275, chi2r: 0.0509969465408591
message: Optimization terminated successfully.
Planet 0 simplex result: (r, theta, f)=(30.416, 239.945, 401.494) at
(X,Y)=(34.77, 23.67)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time: 0:00:55.223076
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
For both the \(n_{\rm pc} = 5\) and \(n_{\rm pc} = 25\) cases, the parameters estimated by the Nelder-Mead optimization are not exactly equal to the original values (radius=30.5, theta=240, flux=400), which reflects:
the limitations of this heuristic minization procedure (depending on the initial guess the minimization may get trapped in a different local minimum);
the higher residual speckle noise level in images obtained with low \(n_{\rm pc}\) values;
the higher self-subtraction for high \(n_{\rm pc}\) values.
These estimates are provided without uncertainties (error bars). We will come back to this question later on.
For comparison, let’s use the optimal \(n_{\rm pc} = 10\) found in Section 5.2:
r_0, theta_0, f_0 = firstguess(cubefc, angs, psfn, ncomp=opt_npc_ann, planets_xy_coord=[xy_test],
fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
aperture_radius=2, imlib=imlib_rot, interpolation=interpolation,
simplex=True, plot=True, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 14:56:00
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux | chi2r
1/30 0.100 2.020
2/30 0.149 2.020
3/30 0.221 2.019
4/30 0.329 2.018
5/30 0.489 2.016
6/30 0.728 2.014
7/30 1.083 2.011
8/30 1.610 2.006
9/30 2.395 1.999
10/30 3.562 1.988
11/30 5.298 1.972
12/30 7.880 1.949
13/30 11.721 1.915
14/30 17.433 1.866
15/30 25.929 1.796
16/30 38.566 1.695
17/30 57.362 1.538
18/30 85.317 1.326
19/30 126.896 1.042
20/30 188.739 0.650
21/30 280.722 0.249
22/30 417.532 0.056
23/30 621.017 0.642
24/30 923.671 3.220
25/30 1373.824 10.306
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 93, nfev: 199, chi2r: 0.05341972125116541
message: Optimization terminated successfully.
Planet 0 simplex result: (r, theta, f)=(30.485, 240.031, 404.544) at
(X,Y)=(34.77, 23.59)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time: 0:00:40.618912
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
We see that using the optimal \(n_{\rm pc}\) leads to closer parameter estimates to the ground truth.
Answer 5.1: If relying on a 2D Gaussian fit, both the photometry and astrometry would be biased by self-subtraction and the negative side lobes.
5.3.2. Planet subtraction
Let’s use the values obtained with the simplex optimization to subtract the planet with the function cube_planet_free.
First we define a list with the parameters (r, theta, flux) is each companion that we obtained via the NegFC, in this case one:
plpar_fc = [(r_0[0], theta_0[0], f_0[0])]
r_0, theta_0 and f_0 have the same length as the number of planet coordinates planets_xy_coord provided to firstguess. Here there is only one planet, so we take the zeroth index. The number of tuples in plpar_fc should match the number of planets.
from vip_hci.fm import cube_planet_free
cube_emp = cube_planet_free(plpar_fc, cubefc, angs, psfn, imlib=imlib_rot, interpolation=interpolation)
Let’s double-check the fake companion was well removed by computing a PCA post-processed image:
from vip_hci.psfsub import pca_annulus
fr_pca_emp = pca_annulus(cube_emp, angs, ncomp=opt_npc_ann, annulus_width=4*fwhm_naco,
r_guess=rad_fc, imlib=imlib_rot, interpolation=interpolation)
Let’s take a look at the PSF of the planet in the full-frame PCA final image and the same PSF in the frame resulting of processing the planet-subtracted cube:
from vip_hci.preproc import frame_crop
cropped_frame1 = frame_crop(final_ann_opt, cenxy=xy_test, size=15)
New shape: (15, 15)
cropped_frame2 = frame_crop(fr_pca_emp, cenxy=xy_test, size=15)
New shape: (15, 15)
Let’s use both 'mode=surface' and the default image mode of plot_frames to show the residuals in the vicinity of the companion:
plot_frames((cropped_frame1, cropped_frame2), mode='surface', vmax=8)
plot_frames((final_ann_opt, fr_pca_emp), vmin = float(np.amin(final_ann_opt)),
vmax= float(np.amax(final_ann_opt)), axis=False)
Not only the bright point-like signal is subtracted, but so are the negative side lobes. A subtraction not leaving any significant artifact/defect is a good sign that the inferred parameters are correct. However, keep in mind that even for slightly inaccurate parameters, the final image can still look relatively clean. Let’s take for example the parameters inferred with non-optimal \(n_{\rm pc}\):
# planet parameters inferred from npc=5 used to make empty cube, then processed with PCA (npc=10)
plpar_fc = [(r_lo[0], theta_lo[0], f_lo[0])]
print(plpar_fc)
cube_emp = cube_planet_free(plpar_fc, cubefc, angs, psfn, imlib=imlib_rot, interpolation=interpolation)
final_ann_5 = pca_annulus(cubefc, angs, ncomp=5, annulus_width=4*fwhm_naco, r_guess=30.5,
imlib=imlib_rot, interpolation=interpolation)
fr_pca_emp5 = pca_annulus(cube_emp, angs, ncomp=5, annulus_width=4*fwhm_naco, r_guess=30.5,
imlib=imlib_rot, interpolation=interpolation)
plot_frames((final_ann_5, fr_pca_emp5), vmin = float(np.amin(final_ann_5)),
vmax= float(np.amax(final_ann_5)), axis=False)
[(30.52424743445897, 240.12093582564796, 402.64949524446)]
# parameters inferred from npc=25 used to make empty cube, then processed with PCA (npc=10)
plpar_fc = [(r_hi[0], theta_hi[0], f_hi[0])]
print(plpar_fc)
cube_emp = cube_planet_free(plpar_fc, cubefc, angs, psfn, imlib=imlib_rot, interpolation=interpolation)
final_ann_25 = pca_annulus(cubefc, angs, ncomp=25, annulus_width=4*fwhm_naco, r_guess=30.5,
imlib=imlib_rot, interpolation=interpolation)
fr_pca_emp25 = pca_annulus(cube_emp, angs, ncomp=25, annulus_width=4*fwhm_naco, r_guess=30.5,
imlib=imlib_rot, interpolation=interpolation)
plot_frames((final_ann_25, fr_pca_emp25), vmin = float(np.amin(final_ann_25)),
vmax= float(np.amax(final_ann_25)), axis=False)
[(30.41552167699778, 239.94548256947093, 401.49438660651504)]
Inaccurate parameters still leading to an apparently good removal of the companion brings the question of the uncertainties on each of the three parameters characterizing the companion. The next sections are dedicated to this question.
5.3.3. NEGFC technique coupled with MCMC
5.3.3.1. Running the MCMC sampler
Markov Chain Monte Carlo (MCMC) is a robust way of obtaining the flux and position of the companion, and their uncertainties. MCMC samples the parameter space (here r and \(\theta\) coordinates, and flux) using a number of walkers following a Markov chain (i.e. a process where each new sampled point only depends on the previous one) for a certain number of iterations. At each iteration, a new sample of parameters is proposed to each walker to explore the parameter space (these are called moves), based on a given proposal function (e.g., the stretch move proposed in Goodman & Weare 2010). A criterion depending on the ratio of likelihoods calculated for the proposed vs the current parameters is then used to accept or reject the proposed moves for the next iteration (if rejected, the parameters do not change).
In our case the log-likelihood is calculated assuming that a perfect subtraction of the companion leads to pixel intensity residuals following a Gaussian distribution (centered on zero, and scaled by the standard deviation estimated in the annulus except in the area of the companion). A convergence criterion based on autocorrelation can be used (and is used by default in VIP) to stop iterating. Sampling the final posterior distributions of the parameters yields both the most likely parameter values and uncertainties on each of them.
The MCMC implementation in VIP leverages emcee (Foreman-Mackey et al. 2013; check the paper for more details on this implementation of MCMC), which is based on an efficient affine-invariant implementation of MCMC proposed in Goodman & Weare (2010).
The relevant function is mcmc_negfc_sampling, which can accept a number of parameters. Let’s define them in the next few boxes:
Let’s first define observation-related parameters, such as the non-coronagraphic psf, its FWHM and the pixel scale od the detector:
obs_params = {'psfn': psfn,
'fwhm': fwhm_naco}
In NEGFC, PCA in a single annulus is used by default to speed up the algorithm - although other algorithms can be used through the algo parameter. Let’s set the \(n_{\rm pc}\) to the optimal \(n_{\rm pc}\) found in Section 5.2. We set the width of the annulus on which PCA is performed (in pixels) with the annulus_width parameter. We also set a few other algorithm-related parameters in the following box:
annulus_width = 4*fwhm_naco
algo_params = {'algo': pca_annulus,
'ncomp': opt_npc_ann,
'annulus_width': annulus_width,
'svd_mode': 'lapack',
'imlib': imlib_rot,
'interpolation': interpolation}
The choice of log-likelihood expression to be used is determined by mu_sigma and fmerit. If the former is True (default; mu and sigma calculated automatically) or a tuple of 2 values provided manually (corresponding to mean and standard deviation of pixel intensities in an annulus encompassing the companion), the following figure of merit will be used: \(\chi^2 = \sum_j^N \frac{(I_j-\mu)^2}{\sigma^2}\) (as introduced in Christiaens et al. (2021). Otherwise, the choice will be determined by fmerit:
‘sum’ for the sum of absolute residuals (details in Wertz et al. 2017);
‘stddev’ for the standard deviation of residuals (details in Wertz et al. 2017; which can be useful for very faint sources);
‘hessian’ for the determinant of the local Hessian matrix (which containes 2nd order derivatives with respect to spatial coordinates) - the latter being useful when the point source is contained within a more extended signal (see e.g. Christiaens et al. 2024).
Another parameter to set is aperture_radius, which is the radius of the aperture (in fwhm units) in which the residual intensities \(I_j\) are considered.
mu_sigma=True
aperture_radius=2
negfc_params = {'mu_sigma': mu_sigma,
'aperture_radius': aperture_radius}
Parameter initial_state corresponds to the initial first estimation of the planets parameters (r, theta, flux). We set it to the result of the simplex optimization, obtained with optimal \(n_{\rm pc}\).
The MCMC minimization can only run for one companion candidate at a time, hence the first dimension of init should always be 3 (not the number of planets, as opposed to planets_xy_coord in firstguess).
initial_state = np.array([r_0[0], theta_0[0], f_0[0]])
Beware that the MCMC procedure is a CPU intensive procedure and can take several hours when run properly on a standard laptop. We use the affine invariant sampler from emcee which can be run in parallel (nproc sets the number of CPUs to be used). At least 100 walkers are recommended for our MCMC chain, although both the number of walkers and iterations will depend on your dataset.
For the sake of preventing this tutorial to take too long to run, we set the maximum number of iterations to 500, although feel free to significantly increase it in case of non-convergence.
from multiprocessing import cpu_count
nwalkers, itermin, itermax = (100, 200, 500)
mcmc_params = {'nwalkers': nwalkers,
'niteration_min': itermin,
'niteration_limit': itermax,
'bounds': None,
'nproc': cpu_count()//2}
Another update from Christiaens et al. (2021) is that the convergence can now be evaluated based on the auto-correlation time (see more details in the documentation of emcee), instead of the Gelman-Rubin test, which is inappropriate for non-independent samples (such as an MCMC chain).
We set the convergence test to the autocorrelation time based criterion using conv_test='ac' (instead of Gelman-Rubin ‘gb’). We also set the autocorrelation criterion \(N/\tau >= a_c\), where \(N\) is the number of samples and \(\tau\) the autocorrelation time, to ac_c=50 (the value recommended in the emcee documentation). Finally, we set the number of consecutive times the criterion must be met to: ac_count_thr=1, and the maximum gap in number of steps between 2 checks of the convergence criterion to: check_maxgap=50. If the maximum number of iterations niteration_limit is large enough, the chain will stop upon meeting the convergence criterion (spoiler: it won’t be the case here since we chose a small value for niteration_limit).
conv_test, ac_c, ac_count_thr, check_maxgap = ('ac', 50, 1, 50)
conv_params = {'conv_test': conv_test,
'ac_c': ac_c,
'ac_count_thr': ac_count_thr,
'check_maxgap': check_maxgap}
Setting bounds=None does not mean no bounds are used for parameter exploration, but rather that they are set automatically to:
\(r \in [r_0-w_{ann}/2, r_0+w_{ann}/2]\), where \(w_{ann}\) is the
annulus_width,\(\theta \in [\theta_0-\Delta {\rm rot}, \theta_0+\Delta {\rm rot}]\), where \(\Delta {\rm rot}\) is the angle subtended by min(
aperture_radius/2,fwhm) at \(r_0\),\(f \in [0.1*f_0, 2*f_0]\),
where (\(r_0, \theta_0, f_0\)) = initial_state.
If the bounds are provided manually (as a tuple of tuples), they will supersede the automatic setting above.
Now let’s start the sampler. Note that this step is computer intensive and may take a long time to run depending on your machine. Feel free to skip the next cell if you do not wish to run the MCMC or can’t run it in a reasonable time. The results are already saved as ‘MCMC_results’ in the ‘datasets’ folder and will be loaded in the subsequent boxes.
from vip_hci.fm import mcmc_negfc_sampling
chain = mcmc_negfc_sampling(cubefc, angs, **obs_params, **algo_params, **negfc_params,
initial_state=initial_state, **mcmc_params, **conv_params,
display=True, verbosity=2, save=False, output_dir='./')
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 14:56:42
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
MCMC sampler for the NEGFC technique
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
The mean and stddev in the annulus at the radius of the companion (excluding the PA area directly adjacent to it) are -0.11 and 1.91 respectively.
Beginning emcee Ensemble sampler...
emcee Ensemble sampler successful
Start of the MCMC run ...
Step | Duration/step (sec) | Remaining Estimated Time (sec)
0 5.38452 2686.87548
1 4.59892 2290.26017
2 4.56100 2266.81750
3 4.55227 2257.92642
4 4.85771 2404.56398
5 4.63168 2288.04943
6 4.86092 2396.43307
7 4.76470 2344.23338
8 4.69724 2306.34582
9 4.81512 2359.41027
10 4.98254 2436.46059
11 4.75354 2319.72557
12 4.91667 2394.41926
13 4.69853 2283.48704
14 4.57682 2219.75576
15 4.60489 2228.76482
16 4.62360 2233.20073
17 4.66872 2250.32497
18 4.87088 2342.89280
19 4.80223 2305.07088
20 4.77962 2289.43894
21 4.59407 2195.96594
22 4.68773 2236.04530
23 4.63407 2205.81827
24 4.89935 2327.18982
25 5.07721 2406.59707
26 4.68628 2216.60808
27 4.88143 2304.03685
28 4.77505 2249.04714
29 4.67295 2196.28509
30 4.69063 2199.90500
31 4.99498 2337.65204
32 4.94699 2310.24480
33 4.65048 2167.12182
34 4.65989 2166.85117
35 4.72072 2190.41454
36 4.67306 2163.62678
37 4.71334 2177.56262
38 4.77889 2203.06737
39 4.71237 2167.68928
40 4.70464 2159.42838
41 4.76456 2182.16985
42 4.74075 2166.52458
43 4.76376 2172.27502
44 4.74104 2157.17366
45 4.71237 2139.41780
46 4.67580 2118.13604
47 4.91394 2221.10224
48 4.97844 2245.27554
49 4.80138 2160.61920
50 4.92204 2209.99686
ac convergence test in progress...
51 5.24491 2349.72147
52 4.74937 2122.96616
53 4.73836 2113.30990
54 5.14571 2289.83917
55 5.25947 2335.20557
56 4.75371 2105.89486
57 4.88856 2160.74264
58 5.02635 2216.62211
59 4.99691 2198.64128
60 4.90420 2152.94248
61 4.75542 2082.87571
62 4.83580 2113.24547
63 4.71327 2054.98659
64 4.86026 2114.21397
65 4.81044 2087.72966
66 4.72851 2047.44656
67 4.69163 2026.78286
68 4.70587 2028.22997
69 4.89132 2103.26889
70 4.85544 2082.98204
71 4.60144 1969.41760
72 4.64164 1981.98199
73 4.60109 1960.06434
74 4.57487 1944.31975
75 4.55851 1932.80994
76 4.55233 1925.63644
77 4.55205 1920.96468
78 4.83721 2036.46541
79 4.82656 2027.15436
80 4.88665 2047.50677
81 5.10946 2135.75386
82 4.92382 2053.23169
83 5.64881 2349.90538
84 5.11652 2123.35621
85 4.93863 2044.59158
86 5.15782 2130.17883
87 5.30695 2186.46340
88 4.91999 2022.11548
89 5.00217 2050.88806
90 4.70872 1925.86812
91 4.75305 1939.24562
92 4.76915 1941.04527
93 4.80617 1951.30340
94 4.81200 1948.86041
95 4.81130 1943.76480
96 4.83442 1948.27207
97 4.68501 1883.37322
98 4.79674 1923.49394
99 4.69831 1879.32440
100 4.78494 1909.19066
ac convergence test in progress...
101 4.72733 1881.47654
102 4.65608 1848.46376
103 4.59514 1819.67742
104 4.78421 1889.76255
105 4.85414 1912.52919
106 4.72222 1855.83325
107 4.57484 1793.33532
108 4.65914 1821.72178
109 4.70982 1836.83097
110 4.62821 1800.37486
111 4.73918 1838.80339
112 4.90608 1898.65296
113 4.87734 1882.65131
114 4.66049 1794.28827
115 4.60052 1766.59814
116 4.59376 1759.41161
117 4.53660 1732.98044
118 4.59178 1749.46894
119 4.60333 1749.26692
120 4.56285 1729.31826
121 4.51971 1708.45114
122 4.96022 1870.00256
123 4.90884 1845.72459
124 4.82057 1807.71450
125 4.89169 1829.49356
126 4.76509 1777.37708
127 4.83354 1798.07502
128 4.90625 1820.21764
129 4.88365 1806.94939
130 4.75467 1754.47471
131 4.70055 1729.80093
132 4.72847 1735.34812
133 4.69493 1718.34401
134 4.88953 1784.67736
135 4.86971 1772.57298
136 4.82881 1752.85839
137 4.92662 1783.43499
138 4.79563 1731.22387
139 4.85178 1746.63972
140 4.83825 1736.93175
141 4.71695 1688.66989
142 4.74968 1695.63683
143 4.72793 1683.14166
144 4.75864 1689.31827
145 4.80688 1701.63658
146 5.06244 1787.04203
147 4.99496 1758.22627
148 5.09435 1788.11580
149 4.92889 1725.11150
150 4.92015 1717.13305
ac convergence test in progress...
151 5.00282 1740.98310
152 4.88493 1695.06967
153 4.91593 1700.91213
154 4.77789 1648.37171
155 4.83555 1663.42954
156 5.02937 1725.07220
157 4.90933 1678.99189
158 4.86438 1658.75392
159 4.82986 1642.15274
160 4.81359 1631.80701
161 4.73641 1600.90759
162 4.87242 1642.00689
163 4.74811 1595.36563
164 4.78481 1602.91068
165 4.82768 1612.44345
166 4.76680 1587.34307
167 4.73546 1572.17405
168 4.73263 1566.50152
169 4.73836 1563.65748
170 4.81493 1584.11131
171 4.86233 1594.84522
172 4.75568 1555.10638
173 4.76565 1553.60190
174 4.80896 1562.91330
175 4.85004 1571.41296
176 4.74568 1532.85335
177 4.77937 1538.95714
178 4.76938 1530.97034
179 4.78593 1531.49888
180 4.87072 1553.75936
181 4.93445 1569.15446
182 4.71386 1494.29235
183 4.84525 1531.09868
184 4.72988 1489.91283
185 4.71883 1481.71136
186 4.62682 1448.19529
187 4.71993 1472.61910
188 4.61409 1434.98261
189 4.65758 1443.84980
190 4.62122 1427.95791
191 4.61300 1420.80277
192 4.67370 1434.82559
193 4.59060 1404.72299
194 4.62151 1409.56055
195 4.66666 1418.66403
196 4.71300 1428.04051
197 4.70877 1422.04733
198 4.61940 1390.43820
199 4.63209 1389.62700
200 4.76439 1424.55201
ac convergence test in progress...
Auto-corr tau/N = [0.07467305 0.07568022 0.06551593]
tau/N <= 0.02 = [False False False]
201 4.88279 1455.07231
202 4.67150 1387.43580
203 4.65430 1377.67280
204 4.71665 1391.41057
205 4.63834 1363.67196
206 4.63013 1356.62838
207 4.64709 1356.95057
208 4.62177 1344.93565
209 4.68656 1359.10095
210 4.73939 1369.68487
211 4.65887 1341.75370
212 4.55735 1307.96060
213 4.73273 1353.55992
214 4.69189 1337.19007
215 4.78024 1357.58844
216 4.61769 1306.80740
217 4.65832 1313.64596
218 4.64043 1303.96167
219 4.66501 1306.20280
220 4.62931 1291.57805
221 4.66826 1297.77600
222 4.68229 1296.99488
223 4.61152 1272.78007
224 4.61932 1270.31217
225 4.63119 1268.94496
226 4.67309 1275.75248
227 4.58505 1247.13224
228 4.81440 1304.70348
229 4.61796 1246.84947
230 4.61995 1242.76574
231 4.73597 1269.23862
232 4.63080 1236.42440
233 4.61560 1227.75040
234 4.63710 1228.83097
235 4.67892 1235.23541
236 4.64575 1221.83304
237 4.72988 1239.22908
238 4.62316 1206.64424
239 4.71568 1226.07784
240 4.64761 1203.72995
241 4.59976 1186.73937
242 4.60426 1183.29379
243 4.59821 1177.14202
244 4.62741 1179.98853
245 4.64387 1179.54323
246 4.63010 1171.41454
247 4.60612 1160.74325
248 4.70525 1181.01675
249 4.58921 1147.30300
250 4.59363 1143.81287
ac convergence test in progress...
Auto-corr tau/N = [0.06984828 0.06936632 0.06165038]
tau/N <= 0.02 = [False False False]
251 4.69803 1165.11218
252 4.72776 1167.75672
253 4.60511 1132.85706
254 4.55993 1117.18358
255 4.59339 1120.78667
256 4.62808 1124.62344
257 4.65668 1126.91559
258 4.68909 1130.06948
259 4.61821 1108.37040
260 4.58471 1095.74497
261 4.67972 1113.77336
262 4.62340 1095.74651
263 4.60050 1085.71847
264 4.66823 1097.03499
265 4.63880 1085.47920
266 4.64011 1081.14610
267 4.63230 1074.69476
268 4.57880 1057.70326
269 4.62770 1064.37123
270 4.60804 1055.24162
271 4.65626 1061.62774
272 4.62971 1050.94462
273 4.60932 1041.70632
274 4.73352 1065.04245
275 4.61123 1032.91530
276 4.66186 1039.59545
277 4.87654 1082.59121
278 4.63955 1025.34011
279 4.57284 1006.02546
280 4.63231 1014.47589
281 4.59963 1002.71999
282 4.60875 1000.09875
283 4.60097 993.81038
284 4.60727 990.56241
285 4.64207 993.40319
286 4.59921 979.63109
287 4.68932 994.13584
288 4.61116 972.95476
289 4.59956 965.90676
290 4.62812 967.27687
291 4.63192 963.43832
292 4.59156 950.45354
293 4.60021 947.64305
294 4.60616 944.26239
295 4.61727 941.92288
296 4.60080 933.96240
297 4.61777 932.78974
298 4.61230 927.07170
299 4.69447 938.89360
300 4.70990 937.26950
ac convergence test in progress...
Auto-corr tau/N = [0.06350204 0.06246719 0.05814341]
tau/N <= 0.02 = [False False False]
301 4.74276 939.06608
302 4.62655 911.43035
303 4.71717 924.56493
304 4.73252 922.84140
305 4.68018 907.95453
306 4.69382 905.90784
307 4.63737 890.37504
308 4.65383 888.88229
309 4.63455 880.56374
310 4.71078 890.33780
311 4.62546 869.58704
312 4.62727 865.29930
313 4.79187 891.28801
314 4.64051 858.49435
315 4.65656 856.80686
316 4.77129 873.14644
317 4.63738 844.00261
318 4.67120 845.48702
319 4.64751 836.55252
320 4.65933 834.01971
321 4.60528 819.74020
322 4.65525 823.97854
323 4.74881 835.79126
324 4.80838 841.46598
325 4.70696 819.01104
326 4.71533 815.75140
327 4.62067 794.75472
328 4.60526 787.49963
329 4.64448 789.56143
330 4.62914 782.32466
331 4.62522 777.03746
332 4.59279 766.99660
333 4.58059 760.37877
334 4.59634 758.39561
335 4.60815 755.73627
336 4.63765 755.93744
337 4.67344 757.09793
338 4.66208 750.59552
339 4.64914 743.86208
340 4.63484 736.93892
341 4.59197 725.53079
342 4.65968 731.57023
343 4.62165 720.97818
344 4.60647 714.00285
345 4.60982 709.91228
346 4.62336 707.37423
347 4.57153 694.87210
348 4.59750 694.22265
349 4.61661 692.49105
350 4.59189 684.19087
ac convergence test in progress...
Auto-corr tau/N = [0.05751351 0.05541949 0.05269437]
tau/N <= 0.02 = [False False False]
351 4.72470 699.25501
352 4.68707 688.99900
353 4.60653 672.55323
354 4.64552 673.59982
355 4.72590 680.52888
356 4.65581 665.78069
357 4.65814 661.45560
358 4.65192 655.92114
359 4.60586 644.82110
360 4.60832 640.55620
361 4.58688 632.98903
362 4.61461 632.20225
363 4.59333 624.69315
364 4.62452 624.31007
365 4.65306 623.50964
366 4.63526 616.48931
367 4.63496 611.81525
368 4.78062 626.26135
369 4.69037 609.74849
370 4.63202 597.53032
371 4.68234 599.33965
372 4.67286 593.45347
373 4.64408 585.15458
374 4.65230 581.53750
375 4.69330 581.96870
376 4.74997 584.24631
377 4.81440 587.35741
378 4.64160 561.63348
379 4.62912 555.49452
380 4.73531 563.50141
381 4.73234 558.41588
382 4.66268 545.53391
383 4.62737 536.77504
384 4.63091 532.55442
385 4.57122 521.11862
386 4.62483 522.60579
387 4.63159 518.73752
388 4.66882 518.23880
389 4.64424 510.86662
390 4.76997 519.92662
391 4.68035 505.47780
392 4.68361 501.14627
393 4.83497 512.50650
394 4.79338 503.30500
395 4.65745 484.37459
396 4.62186 476.05199
397 4.62872 472.12954
398 4.61691 466.30821
399 4.61298 461.29810
400 4.62737 458.10983
ac convergence test in progress...
Auto-corr tau/N = [0.05359522 0.04870759 0.04947033]
tau/N <= 0.02 = [False False False]
401 4.91972 482.13227
402 4.65147 451.19278
403 4.68780 450.02899
404 4.62830 439.68812
405 4.64584 436.70887
406 4.92352 457.88689
407 4.79269 440.92785
408 4.62356 420.74405
409 4.61622 415.46007
410 4.62240 411.39342
411 4.57913 402.96300
412 4.63654 403.37933
413 4.64900 399.81443
414 4.65622 395.77887
415 4.61292 387.48520
416 4.69400 389.60242
417 4.76855 391.02143
418 4.75225 384.93185
419 4.73953 379.16200
420 4.71772 372.69956
421 4.82175 376.09642
422 4.70989 362.66184
423 4.84308 368.07400
424 4.70171 352.62810
425 4.59976 340.38246
426 4.63826 338.59313
427 4.56694 328.81975
428 4.69974 333.68119
429 4.72390 330.67293
430 4.68832 323.49442
431 4.72154 321.06506
432 4.65084 311.60635
433 4.61432 304.54532
434 4.62529 300.64398
435 4.62864 296.23290
436 4.61694 290.86728
437 4.62979 287.04710
438 4.66561 284.60190
439 4.65244 279.14646
440 4.68421 276.36810
441 4.63921 269.07418
442 4.84653 276.25244
443 4.64199 259.95144
444 5.10791 280.93516
445 5.02415 271.30421
446 4.85856 257.50379
447 4.82680 250.99339
448 4.72522 240.98597
449 4.68971 234.48565
450 4.85456 237.87324
ac convergence test in progress...
Auto-corr tau/N = [0.04976504 0.04482726 0.04643542]
tau/N <= 0.02 = [False False False]
451 4.71593 226.36488
452 4.52401 212.62842
453 4.73951 218.01728
454 4.62694 208.21212
455 4.59577 202.21379
456 4.58357 197.09330
457 4.61494 193.82735
458 4.82065 197.64681
459 4.76357 190.54272
460 4.74415 185.02173
461 4.77621 181.49590
462 4.84655 179.32242
463 4.94822 178.13610
464 4.88507 170.97745
465 4.72117 160.51968
466 4.61314 152.23355
467 4.66772 149.36688
468 4.63330 143.63221
469 4.68055 140.41647
470 4.80931 139.47008
471 4.62243 129.42796
472 4.60030 124.20802
473 4.60043 119.61108
474 4.87580 121.89500
475 5.20221 124.85314
476 4.78261 110.00001
477 4.78126 105.18781
478 4.80611 100.92829
479 4.96092 99.21834
480 4.81557 91.49587
481 4.74328 85.37908
482 4.80668 81.71353
483 4.79310 76.68958
484 4.80662 72.09931
485 4.67095 65.39331
486 4.77275 62.04571
487 4.68536 56.22432
488 4.59662 50.56285
489 4.59012 45.90116
490 4.57718 41.19461
491 4.92664 39.41310
492 5.50180 38.51257
493 4.93987 29.63922
494 4.89804 24.49019
495 5.10269 20.41074
496 4.94586 14.83759
497 5.03041 10.06082
498 5.34860 5.34860
499 4.93600 0.00000
We have reached the limit # of steps without convergence
Running time: 0:39:30.261013
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
If you ran the previous box and wish to write your results, set write=True in the next box. This will pickle the MCMC chain.
write=False
if write:
import pickle
output = {'chain':chain}
with open('../datasets/my_MCMC_results', 'wb') as fileSave:
pickle.dump(output, fileSave)
An alternative to the cell above, is to provide output_dir in the call to mcmc_negfc_sampling and set save to True. This will save the results as a pickle including additional keys: apart from ‘chain’, it will also save ‘input_parameters’, ‘AR’ (acceptance ratio), and ‘lnprobability’.
Pickled results can be loaded from disk like this:
import pickle
with open('../datasets/MCMC_results','rb') as fi:
myPickler = pickle.Unpickler(fi)
mcmc_result = myPickler.load()
print(mcmc_result.keys())
chain = mcmc_result['chain']
dict_keys(['chain', 'input_parameters', 'AR', 'lnprobability'])
The most accurate approach would involve setting a large enough maximum number of iterations and using FFT-based rotation for PCA. The latter in particular, may change a bit the most likely parameter values given the better flux conservation. However, these changes would involve over ~3 orders of magnitude longer calculation time. It is therefore intractable for a personal laptop and not shown in this notebook. If you have access to a supercomputer feel free to adopt these changes though. The results after 500 iterations are nonetheless good enough for illustrative purpose:
5.3.3.2. Visualizing the MCMC chain: corner plots and walk plots
Let’s first check that the walk plots look ok:
from vip_hci.fm import show_walk_plot
show_walk_plot(chain)
Then based on the walk plot, let’s burn-in the first 30% of the chain (to only keep the iterations where the walkers appear to have mostly converged), to calculate the corner plots:
from vip_hci.fm import show_corner_plot
burnin = 0.3
burned_chain = chain[:, int(chain.shape[1]//(1/burnin)):, :]
show_corner_plot(burned_chain)
For the purpose of this tutorial and to limit computation time we set the maximum number of iterations to 500 for 100 walkers. The look of the corner plots may improve with more samples (i.e. higher number of iterations, for a given burn-in ratio). This can be tested by setting the max. number of iterations arbitrarily high for the autocorrelation-time convergence criterion to be met.
5.3.3.3. Highly probable values and confidence intervals
Now let’s determine the most highly probable value for each model parameter, as well as the 1-sigma confidence interval. For this, let’s first flatten the chains:
isamples_flat = chain[:, int(chain.shape[1]//(1/burnin)):, :].reshape((-1,3))
Then use the confidence function:
from vip_hci.fm import confidence
val_max, conf = confidence(isamples_flat, cfd=68, gaussian_fit=False, verbose=True, save=False,
gt=gt, ndig=1, title=True)
percentage for r: 69.79202279202278%
percentage for theta: 69.02279202279202%
percentage for f: 68.94017094017094%
Confidence intervals:
r: 30.57489764243517 [-0.20421878422770945,0.08279139901123145]
theta: 239.91771004343246 [-0.14527599962792692,0.18317408648738365]
f: 390.7650582608513 [-20.054061636888093,25.285555976945886]
Using the confidence function with the gaussian_fit=True option, it is possible to fit a Gaussian to the posterior distribution of each parameter, and infer associated uncertainty values.
mu, sigma = confidence(isamples_flat, cfd=68, gaussian_fit=True, verbose=True, save=False,
gt=gt, ndig=1, title=True)
percentage for r: 69.79202279202278%
percentage for theta: 69.02279202279202%
percentage for f: 68.94017094017094%
Confidence intervals:
r: 30.57489764243517 [-0.20421878422770945,0.08279139901123145]
theta: 239.91771004343246 [-0.14527599962792692,0.18317408648738365]
f: 390.7650582608513 [-20.054061636888093,25.285555976945886]
Gaussian fit results:
r: 30.525747462344565 +-0.13676103149725966
theta: 239.94793373875024 +-0.15744912409220196
f: 394.1338634773366 +-21.669634439677626
It is recommended to report the results as confidence intervals (i.e. with possibly asymmetric uncertainties) as long as the bin interval is small enough. Here, we also fitted the residual posterior distribution of each parameter to a Gaussian distribution (this shape is the expected one if the noise has been well whitened, but is not necessarily guaranteed at all separations depending on the adopted \(n_{\rm pc}\)). In case these distributions look Gaussian, the inferred \(\sigma\) value may be a more accurate uncertainty value for the different parameters.
We can see that the confidence intervals inferred by NEGFC for the different parameters encompass the ground truth values used for injection (in particular for the flux).
Depending on your choice of mu_sigma, you may have to calculate separately another source of uncertainty. Indeed, with the original expression for \(\chi^2\) (Wertz et al. 2017; mu_sigma=False), only the uncertainty associated to photon noise was reflected in the MCMC results. With the new expression (mu_sigma=True), the \(\chi^2\) expression also takes into account the residual speckle noise at the radial separation of the companion candidate (at least if sigma='spe+pho' in the call, which is the default). Our tests on this dataset suggest that similar final uncertainties can be obtained in either of these 2 ways:
the uncertainties obtained with MCMC when setting
mu_sigma=Trueandsigma='spe+pho';the uncertainties obtained by combining quadratically the uncertainties obtained with MCMC (setting
mu_sigma=Falseandfmerit = 'sum'), and the residual speckle uncertainties inferred as in Sec. 5.3.5 (Wertz et al. 2017).
5.3.4. NEGFC technique coupled with nested sampling
Nested sampling can be a faster approach to estimate the optimal parameters, given a good first guess estimate.
The implementation in VIP relies on nestle (Barbary 2013; more details here).
The procedure is essentially the same as for the MCMC sampler. Let’s first set the parameters associated with the observation:
obs_params = {'psfn': psfn,
'fwhm': fwhm_naco}
Let’s then set the parameters associated with pca_annulus, the PSF modeling and subtraction algorithm:
annulus_width = 4*fwhm_naco
algo_params = {'algo': pca_annulus,
'ncomp': opt_npc_ann,
'annulus_width': annulus_width,
'svd_mode': 'lapack'}
algo_options = {'imlib': imlib_rot,
'interpolation': interpolation}
For a change, let’s consider fmerit='sum' for the choice of log-likelihood expression, hence set mu_sigma=False. Let’s keep an aperture_radius of 2 FWHM.
aperture_radius=2
mu_sigma = False
fmerit = 'sum'
negfc_params = {'mu_sigma': mu_sigma,
'fmerit': fmerit,
'aperture_radius': aperture_radius}
Let’s define our first guess on the parameters as the initial_state:
initial_state = np.array([r_0[0], theta_0[0], f_0[0]])
Let’s now define the parameters of the nested sampling algorithm. These are a bit different from MCMC:
npointssets the number of live points (similar to walkers in MCMC);dlogzis related to a criterion used to stop the iterations. 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_{\rm est}\)) - log(\(z\)) < \(dlogz\) where \(z\) is the current evidence from all saved samples, and \(z_{\rm est}\) is the estimated contribution from the remaining volume. Typical good values: 0.1 or 0.5;decline_factoris related to an alternative criterion to stop the iterations. Iterations 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 yields good results on this dataset.ware the windows that will bound the exploration space for each parameter. These should be provided in the respective units of each input parameter.rstatecan be used to specify a given random state to start with.
npoints = 100
dlogz, decline_factor = (0.1, None)
rstate = None
w = (5, 5, 200)
nested_params = {'npoints': nwalkers,
'dlogz': dlogz,
'decline_factor': None,
'rstate': rstate,
'w': w}
Now let’s start the sampler. Note that this step is computer intensive and may take a few minutes to run depending on your machine (but typically goes 4-5x faster than MCMC with 500 iterations considering the aboved parameters).
from vip_hci.fm import nested_negfc_sampling
nested_res = nested_negfc_sampling(initial_state, cubefc, angs, **obs_params, **algo_params, algo_options=algo_options,
**negfc_params, **nested_params, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2024-07-05 15:36:17
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Prior bounds on parameters:
Radius [25.48510847427911,35.48510847427911]
Theta [235.03055781801652,245.03055781801652]
Flux [204.54399679829373,604.5439967982937]
Using 100 active points
Total running time:
Running time: 0:08:36.916331
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
The results of the nested sampling as stored in a nestle object. The results can be conveniently inferred using the aptly named nested_sampling_results function, which can also write the results (through save and output_dir arguments):
from vip_hci.fm import nested_sampling_results
nested_sampling_results(nested_res, burnin=0.6, bins=None, plot=True, save=True, output_dir='../datasets/')
niter: 1339
ncall: 2780
nsamples: 1439
logz: -62.092 +/- 0.300
h: 8.977
Natural log of prior volume and Weight corresponding to each sample
Walk plots before the burnin
Walk plots after the burnin
Weighted mean +- sqrt(covariance)
Radius = 30.474 +/- 0.066
Theta = 240.021 +/- 0.070
Flux = 407.107 +/- 8.871
Hist bins = 24
Confidence intervals
percentage for r: 69.27335675180683%
percentage for theta: 69.77654349210385%
percentage for f: 68.23043668462095%
Confidence intervals:
r: 30.46830795641032 [-0.06589005049950103,0.08053228394383183]
theta: 239.9975984089133 [-0.0726253182261587,0.08876427783195595]
f: 403.7622447582137 [-7.7056418956507855,14.31047780620861]
Gaussian fit results:
r: 30.475928702039166 +-0.06076429844936294
theta: 240.0156475464783 +-0.06497099285067368
f: 407.28288429071847 +-8.181653306500888
array([[3.04738302e+01, 6.63395783e-02],
[2.40020796e+02, 6.99780374e-02],
[4.07106584e+02, 8.87101361e+00]])
As a reminder the ground truth values are:
gt
[30.5, 240, 400.0]
The final posterior distributions on the three parameters appear thus reasonable.
5.3.5. Residual speckle uncertainty
This section should only be relevant if using mu_sigma=False in your call to mcmc_negfc_sampling. If mu_sigma=True, the uncertainties inferred with the MCMC should already encompass the residual speckle uncertainty.
Residual speckle noise can bias the best parameter estimates found for the companion (if not taken into account in the MCMC).
To evaluate the uncertainty associated to this additional source of noise, it is recommended to inject a large number of fake companions at the same radius and flux as estimated for the true companion but different azimuths, and then estimate their parameters using simplex-NEGFC. The distribution of differences with respect to injected parameters can then give us an idea of the residual speckle noise uncertainty. This is done in VIP with the speckle_noise_uncertainty function (see also Sec. 3.3 in Wertz et al. 2017 for more details).
Let’s use the planet parameters inferred by the MCMC-NEGFC algorithm:
pl_par = (val_max['r'],val_max['theta'],val_max['f'])
pl_par
(30.57489764243517, 239.91771004343246, 390.7650582608513)
from vip_hci.fm import speckle_noise_uncertainty
algo_options={'ncomp':opt_npc_ann, 'annulus_width':4*fwhm_naco, 'imlib':imlib_rot,
'interpolation':interpolation}
speckle_res = speckle_noise_uncertainty(cubefc, pl_par, np.linspace(0,359,360), angs, pca_annulus,
psfn, fwhm_naco, aperture_radius=2, fmerit='sum',
algo_options=algo_options, transmission=None, mu_sigma=None,
wedge=None, weights=None, force_rPA=False, nproc=None,
simplex_options=None, bins=None, save=False, output=None,
verbose=True, full_output=True, plot=True)
#######################################################
### SPECKLE NOISE DETERMINATION ###
#######################################################
Number of steps: 360
Process is running for angle: 0.00
Process is running for angle: 18.00
Process is running for angle: 36.00
Process is running for angle: 54.00
Process is running for angle: 72.00
Process is running for angle: 1.00
Process is running for angle: 19.00
Process is running for angle: 73.00
Process is running for angle: 55.00
Process is running for angle: 37.00
Process is running for angle: 74.00
Process is running for angle: 2.00
Process is running for angle: 56.00
Process is running for angle: 38.00
Process is running for angle: 20.00
Process is running for angle: 75.00
Process is running for angle: 57.00
Process is running for angle: 39.00
Process is running for angle: 3.00
Process is running for angle: 21.00
Process is running for angle: 76.00
Process is running for angle: 4.00
Process is running for angle: 40.00
Process is running for angle: 58.00
Process is running for angle: 5.00
Process is running for angle: 22.00
Process is running for angle: 77.00
Process is running for angle: 41.00
Process is running for angle: 59.00
Process is running for angle: 23.00
Process is running for angle: 6.00
Process is running for angle: 78.00
Process is running for angle: 42.00
Process is running for angle: 60.00
Process is running for angle: 24.00
Process is running for angle: 7.00
Process is running for angle: 79.00
Process is running for angle: 43.00
Process is running for angle: 61.00
Process is running for angle: 25.00
Process is running for angle: 80.00
Process is running for angle: 44.00
Process is running for angle: 62.00
Process is running for angle: 8.00
Process is running for angle: 26.00
Process is running for angle: 45.00
Process is running for angle: 81.00
Process is running for angle: 63.00
Process is running for angle: 27.00
Process is running for angle: 9.00
Process is running for angle: 46.00
Process is running for angle: 82.00
Process is running for angle: 28.00
Process is running for angle: 10.00
Process is running for angle: 64.00
Process is running for angle: 47.00
Process is running for angle: 83.00
Process is running for angle: 29.00
Process is running for angle: 65.00
Process is running for angle: 11.00
Process is running for angle: 84.00
Process is running for angle: 48.00
Process is running for angle: 30.00
Process is running for angle: 12.00
Process is running for angle: 66.00
Process is running for angle: 85.00
Process is running for angle: 49.00
Process is running for angle: 13.00
Process is running for angle: 31.00
Process is running for angle: 67.00
Process is running for angle: 86.00
Process is running for angle: 50.00
Process is running for angle: 14.00
Process is running for angle: 32.00
Process is running for angle: 68.00
Process is running for angle: 87.00
Process is running for angle: 51.00
Process is running for angle: 33.00
Process is running for angle: 15.00
Process is running for angle: 69.00
Process is running for angle: 88.00
Process is running for angle: 52.00
Process is running for angle: 34.00
Process is running for angle: 16.00
Process is running for angle: 70.00
Process is running for angle: 89.00
Process is running for angle: 35.00
Process is running for angle: 17.00
Process is running for angle: 53.00
Process is running for angle: 71.00
Process is running for angle: 90.00
Process is running for angle: 108.00
Process is running for angle: 126.00
Process is running for angle: 144.00
Process is running for angle: 162.00
Process is running for angle: 91.00
Process is running for angle: 109.00
Process is running for angle: 127.00
Process is running for angle: 145.00
Process is running for angle: 163.00
Process is running for angle: 92.00
Process is running for angle: 128.00
Process is running for angle: 146.00
Process is running for angle: 110.00
Process is running for angle: 164.00
Process is running for angle: 93.00
Process is running for angle: 147.00
Process is running for angle: 129.00
Process is running for angle: 111.00
Process is running for angle: 165.00
Process is running for angle: 148.00
Process is running for angle: 130.00
Process is running for angle: 94.00
Process is running for angle: 112.00
Process is running for angle: 166.00
Process is running for angle: 131.00
Process is running for angle: 149.00
Process is running for angle: 95.00
Process is running for angle: 113.00
Process is running for angle: 150.00
Process is running for angle: 132.00
Process is running for angle: 96.00
Process is running for angle: 167.00
Process is running for angle: 114.00
Process is running for angle: 151.00
Process is running for angle: 133.00
Process is running for angle: 168.00
Process is running for angle: 97.00
Process is running for angle: 115.00
Process is running for angle: 134.00
Process is running for angle: 152.00
Process is running for angle: 169.00
Process is running for angle: 98.00
Process is running for angle: 116.00
Process is running for angle: 135.00
Process is running for angle: 153.00
Process is running for angle: 99.00
Process is running for angle: 170.00
Process is running for angle: 117.00
Process is running for angle: 136.00
Process is running for angle: 154.00
Process is running for angle: 118.00
Process is running for angle: 100.00
Process is running for angle: 171.00
Process is running for angle: 137.00
Process is running for angle: 172.00
Process is running for angle: 155.00
Process is running for angle: 101.00
Process is running for angle: 119.00
Process is running for angle: 138.00
Process is running for angle: 173.00
Process is running for angle: 156.00
Process is running for angle: 120.00
Process is running for angle: 102.00
Process is running for angle: 139.00
Process is running for angle: 157.00
Process is running for angle: 174.00
Process is running for angle: 121.00
Process is running for angle: 103.00
Process is running for angle: 140.00
Process is running for angle: 122.00
Process is running for angle: 158.00
Process is running for angle: 104.00
Process is running for angle: 175.00
Process is running for angle: 141.00
Process is running for angle: 159.00
Process is running for angle: 105.00
Process is running for angle: 123.00
Process is running for angle: 176.00
Process is running for angle: 142.00
Process is running for angle: 160.00
Process is running for angle: 106.00
Process is running for angle: 124.00
Process is running for angle: 143.00
Process is running for angle: 177.00
Process is running for angle: 107.00
Process is running for angle: 161.00
Process is running for angle: 125.00
Process is running for angle: 180.00
Process is running for angle: 178.00
Process is running for angle: 198.00
Process is running for angle: 216.00
Process is running for angle: 234.00
Process is running for angle: 179.00
Process is running for angle: 181.00
Process is running for angle: 199.00
Process is running for angle: 235.00
Process is running for angle: 217.00
Process is running for angle: 252.00
Process is running for angle: 182.00
Process is running for angle: 200.00
Process is running for angle: 236.00
Process is running for angle: 218.00
Process is running for angle: 253.00
Process is running for angle: 183.00
Process is running for angle: 201.00
Process is running for angle: 219.00
Process is running for angle: 237.00
Process is running for angle: 254.00
Process is running for angle: 184.00
Process is running for angle: 202.00
Process is running for angle: 220.00
Process is running for angle: 185.00
Process is running for angle: 255.00
Process is running for angle: 238.00
Process is running for angle: 221.00
Process is running for angle: 203.00
Process is running for angle: 239.00
Process is running for angle: 186.00
Process is running for angle: 256.00
Process is running for angle: 222.00
Process is running for angle: 204.00
Process is running for angle: 240.00
Process is running for angle: 257.00
Process is running for angle: 187.00
Process is running for angle: 205.00
Process is running for angle: 223.00
Process is running for angle: 241.00
Process is running for angle: 188.00
Process is running for angle: 258.00
Process is running for angle: 224.00
Process is running for angle: 189.00
Process is running for angle: 206.00
Process is running for angle: 242.00
Process is running for angle: 259.00
Process is running for angle: 225.00
Process is running for angle: 190.00
Process is running for angle: 243.00
Process is running for angle: 207.00
Process is running for angle: 191.00
Process is running for angle: 260.00
Process is running for angle: 226.00
Process is running for angle: 244.00
Process is running for angle: 208.00
Process is running for angle: 192.00
Process is running for angle: 227.00
Process is running for angle: 261.00
Process is running for angle: 209.00
Process is running for angle: 245.00
Process is running for angle: 193.00
Process is running for angle: 228.00
Process is running for angle: 210.00
Process is running for angle: 262.00
Process is running for angle: 246.00
Process is running for angle: 229.00
Process is running for angle: 263.00
Process is running for angle: 247.00
Process is running for angle: 211.00
Process is running for angle: 194.00
Process is running for angle: 230.00
Process is running for angle: 264.00
Process is running for angle: 248.00
Process is running for angle: 195.00
Process is running for angle: 212.00
Process is running for angle: 249.00
Process is running for angle: 196.00
Process is running for angle: 265.00
Process is running for angle: 231.00
Process is running for angle: 213.00
Process is running for angle: 197.00
Process is running for angle: 250.00
Process is running for angle: 266.00
Process is running for angle: 232.00
Process is running for angle: 214.00
Process is running for angle: 270.00
Process is running for angle: 251.00
Process is running for angle: 267.00
Process is running for angle: 233.00
Process is running for angle: 215.00
Process is running for angle: 271.00
Process is running for angle: 288.00
Process is running for angle: 268.00
Process is running for angle: 306.00
Process is running for angle: 324.00
Process is running for angle: 272.00
Process is running for angle: 289.00
Process is running for angle: 269.00
Process is running for angle: 325.00
Process is running for angle: 307.00
Process is running for angle: 273.00
Process is running for angle: 342.00
Process is running for angle: 326.00
Process is running for angle: 290.00
Process is running for angle: 308.00
Process is running for angle: 274.00
Process is running for angle: 327.00
Process is running for angle: 309.00
Process is running for angle: 343.00
Process is running for angle: 291.00
Process is running for angle: 275.00
Process is running for angle: 328.00
Process is running for angle: 310.00
Process is running for angle: 292.00
Process is running for angle: 344.00
Process is running for angle: 276.00
Process is running for angle: 329.00
Process is running for angle: 311.00
Process is running for angle: 293.00
Process is running for angle: 330.00
Process is running for angle: 277.00
Process is running for angle: 345.00
Process is running for angle: 312.00
Process is running for angle: 278.00
Process is running for angle: 331.00
Process is running for angle: 294.00
Process is running for angle: 346.00
Process is running for angle: 313.00
Process is running for angle: 332.00
Process is running for angle: 279.00
Process is running for angle: 295.00
Process is running for angle: 347.00
Process is running for angle: 314.00
Process is running for angle: 333.00
Process is running for angle: 280.00
Process is running for angle: 348.00
Process is running for angle: 334.00
Process is running for angle: 315.00
Process is running for angle: 281.00
Process is running for angle: 296.00
Process is running for angle: 349.00
Process is running for angle: 335.00
Process is running for angle: 316.00
Process is running for angle: 282.00
Process is running for angle: 297.00
Process is running for angle: 350.00
Process is running for angle: 317.00
Process is running for angle: 336.00
Process is running for angle: 283.00
Process is running for angle: 298.00
Process is running for angle: 351.00
Process is running for angle: 337.00
Process is running for angle: 299.00
Process is running for angle: 318.00
Process is running for angle: 284.00
Process is running for angle: 338.00
Process is running for angle: 300.00
Process is running for angle: 319.00
Process is running for angle: 352.00
Process is running for angle: 285.00
Process is running for angle: 339.00
Process is running for angle: 353.00
Process is running for angle: 301.00
Process is running for angle: 320.00
Process is running for angle: 286.00
Process is running for angle: 340.00
Process is running for angle: 302.00
Process is running for angle: 321.00
Process is running for angle: 354.00
Process is running for angle: 287.00
Process is running for angle: 341.00
Process is running for angle: 303.00
Process is running for angle: 322.00
Process is running for angle: 355.00
Process is running for angle: 304.00
Process is running for angle: 323.00
Process is running for angle: 356.00
Process is running for angle: 305.00
Process is running for angle: 357.00
Process is running for angle: 358.00
Process is running for angle: 359.00
residuals (offsets): [-0.07193343 -0.06790936 -0.05285152 -0.04923593 -0.03142681 -0.02385255
-0.02769559 -0.01317097 0.02726336 0.03564376 0.06141544 0.09204904
0.06126068 0.02995647 -0.01622539 0.02291451 -0.02163003 -0.02613889
-0.02893477 -0.02929076 -0.06372168 -0.0197281 -0.01467003 -0.01427048
-0.01221347 -0.00916685 0.00398005 0.02938429 0.07779877 0.0781888
0.06939317 0.09053147 0.05032569 0.03183275 0.01354532 0.00373644
-0.02376853 -0.02876542 -0.01354355 0.04721794 0.03539091 0.05637144
0.04430948 0.00168781 -0.00030264 -0.05071669 -0.02564431 -0.03164453
0.00631446 0.00938729 -0.00834584 0.00956563 -0.04292239 -0.05264633
-0.02109999 -0.02132875 -0.00242952 0.01061854 -0.01486089 -0.02449045
-0.02156495 0.00190762 -0.00276624 0.0142549 0.03277652 0.04592176
0.0560244 0.08423926 0.08633189 0.12877874 0.0922713 0.07715219
0.05143351 0.02068191 -0.04878376 -0.09113553 -0.12921528 -0.10102868
-0.09163093 -0.04142208 -0.04221186 -0.02011948 -0.04003878 -0.03446641
-0.00622856 -0.02041702 -0.01772627 -0.00693319 0.00872119 0.0076161
0.00217512 0.03726963 0.0746907 0.09934815 0.09406521 0.06627059
0.03534415 0.05337205 0.02680105 0.00960021 -0.00946944 -0.01129551
-0.02647772 -0.03811059 -0.05629713 -0.05897208 -0.07139011 -0.04536969
-0.06271726 -0.02781296 0.02468719 0.05119251 0.05123589 0.06789733
0.07151027 0.08151229 0.1110713 0.05337373 0.0261758 0.00485534
0.01120588 -0.01334321 -0.0479408 -0.03418262 -0.04430583 -0.03132003
-0.03644347 -0.02507009 -0.01790798 -0.02298296 -0.00123799 0.05011665
0.04108751 0.06965624 0.07178103 0.08211282 0.07337741 0.07417184
0.05759962 0.0649908 0.04042431 0.00213305 -0.02682168 -0.05879752
-0.10727715 -0.11739961 -0.09815895 -0.11152882 -0.03738389 0.02308846
0.04831898 0.05269665 0.07235345 0.09471892 0.05416403 0.10976661
0.05564165 -0.03064062 -0.0579982 -0.06221158 -0.14410736 -0.10125589
-0.10206235 -0.07056144 -0.05361891 -0.02226898 -0.02205778 0.02954738
0.03823061 0.09898541 0.11111553 0.10515153 0.09699846 0.05616047
0.05736399 0.04320944 -0.04202149 -0.04567831 -0.05966831 -0.05361001
-0.06190871 -0.04837438 -0.0748999 -0.04939728 -0.04518564 -0.05225037
-0.0020218 0.00825943 0.00834982 0.01037694 0.0160477 0.00192957
0.01969101 -0.00077618 -0.00490667 -0.00326904 -0.01430556 0.0089247
0.00104151 0.03650417 0.02488804 0.02988856 0.04818487 0.05118321
0.08654997 0.10554384 0.09648151 0.12608987 0.13003273 0.09309661
0.00133548 -0.04732475 -0.06119426 -0.07463673 -0.06611579 -0.10314215
-0.12428061 -0.09853001 -0.08331502 -0.05461418 -0.02267745 -0.04397129
-0.0096084 0.0034865 0.00377495 0.00734738 0.06305979 0.06686873
0.07328915 0.07522689 0.08158688 0.0971459 0.10488199 0.10120115
0.01672844 -0.02062814 -0.02430726 -0.03577086 -0.04546831 -0.08420737
-0.11448372 -0.12591224 -0.12814383 -0.10918665 -0.04721686 -0.04523558
-0.02726583 0.01020483 0.03360948 0.05267454 0.09625858 0.1135039
0.10037691 0.07060353 0.06822256 0.02929298 0.00656048 0.00512985
-0.03608755 -0.03369547 -0.04881402 -0.06534172 -0.10349118 -0.10075222
-0.10522287 -0.07198134 -0.06712865 -0.06282823 0.00716817 -0.00217294
0.04317417 0.05077055 0.06511812 0.08591605 0.10280419 0.07371422
0.03306598 0.02315526 0.0102798 0.00736466 -0.01273747 -0.0578455
-0.0727186 -0.00946047 -0.03333214 -0.0140088 -0.0156462 0.00879953
0.00448757 0.01787297 0.0191745 0.05155744 -0.00473765 0.05515353
-0.00483744 -0.02581991 -0.03674561 -0.0420111 -0.03738258 -0.07857063
-0.04596219 -0.04187026 -0.01149617 -0.004138 -0.011091 0.04292094
0.04638123 0.00693309 0.00323444 0.02102283 -0.01789094 -0.02828697
-0.0098404 -0.02354448 -0.0236554 -0.02085187 -0.02168203 0.01914278
-0.02126819 0.01214814 -0.02697735 0.0093771 0.0063519 0.00547431
-0.0175158 -0.01005627 -0.00551658 -0.03436807 -0.03913157 -0.04520238
-0.04624307 -0.03602924 -0.01224127 -0.0112659 -0.00969257 -0.00067361
-0.02741018 -0.01297074 -0.03766434 -0.05840367 -0.0490791 -0.08541882
-0.05737633 -0.0379828 0.00168768 0.01642256 0.0688464 0.07479108
0.13220436 0.1426354 0.1057909 0.09995392 0.10196132 0.08381197
0.05855676 0.01539395 0.01412732 -0.00569453 -0.03746173 -0.06495823] [ 4.39884695e-04 6.39556200e-02 5.98256512e-02 9.50812532e-02
1.06327097e-01 7.70431000e-02 7.13869649e-02 6.61612023e-02
4.42774170e-02 2.35712229e-02 -3.19613522e-02 -3.66486706e-02
-5.07865352e-02 -6.62293764e-02 -5.69511297e-02 -5.40233575e-02
1.69960532e-02 2.63203555e-02 2.85330983e-02 3.35459212e-02
2.81769051e-02 -1.21553450e-03 1.06871477e-02 4.18761957e-03
-3.62817333e-04 -9.45729882e-03 -1.46289213e-02 2.57950283e-03
-4.02800608e-02 -5.29462480e-02 -5.31478262e-02 -6.25459628e-02
-4.60560010e-02 -3.04921017e-02 -6.52633019e-02 -2.99657558e-02
-1.05164025e-02 5.39192929e-02 4.90090882e-02 5.68360092e-02
7.42574331e-02 5.66441115e-02 4.71441011e-02 1.16833287e-02
7.85277990e-03 -7.76485306e-02 -7.90427428e-02 -1.22581398e-01
-1.11325816e-01 -1.10589447e-01 -1.15674419e-01 -1.25890452e-03
7.07682838e-02 1.00180504e-01 2.14241036e-01 2.19791873e-01
1.81339738e-01 1.38971654e-01 3.12264396e-02 -2.75149698e-02
-7.83404314e-02 -1.22952812e-01 -1.63560617e-01 -1.56535728e-01
-1.68881714e-01 -9.71079925e-02 -8.18854332e-02 -4.96882975e-04
2.25582450e-02 1.41238849e-02 2.53931726e-02 4.04852222e-02
-3.26146957e-03 -4.90815208e-02 6.98953897e-03 5.74048779e-02
7.51610720e-02 1.25986002e-01 1.53385597e-01 1.11910120e-01
1.23993675e-01 8.85905290e-02 4.27836352e-02 -1.76168073e-03
-4.88598010e-02 -8.78066281e-02 -9.79451871e-02 -7.78703835e-02
-1.08565267e-01 -9.45179146e-02 -1.44123855e-02 3.80652249e-02
9.54096748e-02 6.50763926e-02 8.10722979e-02 2.20814926e-02
4.35934831e-02 1.14938889e-04 -2.65527323e-03 -5.19560170e-02
-2.36916688e-02 -1.88572065e-02 -8.81138268e-03 -1.42300273e-02
-1.58529775e-02 -1.07826897e-02 1.32666004e-02 3.46699578e-03
9.17962680e-03 6.49671782e-02 7.50965511e-02 4.05236373e-02
4.02168029e-02 -3.66687881e-03 -4.59634736e-02 -5.40654161e-02
-5.75725626e-02 -4.75315580e-02 -4.78302084e-02 -3.23929451e-02
2.24127190e-02 -2.17220741e-02 -1.04627880e-02 -3.18274941e-03
8.73428027e-03 -2.42662000e-03 1.03592927e-02 5.37035689e-02
5.78685680e-02 3.53655251e-02 3.87595512e-02 1.73475348e-02
-1.22409962e-02 -3.36277116e-02 -4.98538567e-02 -1.37453650e-02
1.20477252e-02 3.49978356e-02 6.07350950e-02 1.07188010e-03
1.58861980e-02 4.17824425e-02 -1.53413908e-02 -5.69829477e-02
-9.18466829e-02 -8.21675991e-02 -8.47339232e-02 -4.24690913e-02
4.93410377e-02 8.54690001e-02 9.80774501e-02 6.30051321e-02
1.11664121e-01 1.05440485e-01 1.79687642e-02 3.98895744e-02
-1.97410087e-02 -6.75067667e-02 -5.10489801e-02 -4.40763055e-03
-1.34406637e-02 2.92988498e-02 6.30403218e-02 -1.25352664e-02
-4.25639080e-02 -1.31530341e-01 -1.23948753e-01 -5.86359096e-02
-4.50389432e-02 -6.85705926e-03 -1.63057931e-02 1.20861969e-02
3.45054873e-02 5.94887418e-02 5.88559363e-02 3.58124438e-02
4.20346395e-02 4.21012654e-02 1.91090812e-02 -9.17514049e-03
-7.66774934e-03 -4.54534795e-02 -4.77309734e-02 -4.67579748e-02
-3.10119518e-02 -3.24410131e-02 2.36984046e-02 5.77815926e-02
5.49361582e-02 5.78223159e-02 4.85403652e-02 4.04092559e-02
1.64562913e-02 1.81614652e-02 1.24149782e-02 1.11131949e-02
1.28012924e-02 3.79975166e-02 8.31861452e-02 4.95851846e-02
3.34436702e-02 -1.52725155e-02 -6.39608159e-04 -1.24396094e-01
-1.48065698e-01 -1.79105412e-01 -1.57232987e-01 -9.59545903e-02
-4.88704988e-02 -2.24512393e-03 4.10657649e-02 6.22380145e-02
8.81012664e-02 1.23563409e-01 1.28330891e-01 1.39394042e-01
8.70261942e-02 6.73777980e-02 3.83065895e-02 2.57133723e-02
-1.20562034e-02 -5.35972972e-02 -5.86798548e-02 -8.12285354e-02
-1.17525359e-01 -1.29531427e-01 -1.18590960e-01 -1.48369965e-01
-1.41492546e-01 -1.25597547e-01 -1.30768884e-01 1.30716453e-02
3.65518145e-02 8.52386285e-02 1.59648092e-01 1.61674755e-01
1.39484740e-01 1.75266735e-01 1.30631933e-01 8.25820482e-02
8.76164000e-02 7.23649688e-02 1.74515552e-02 3.06143483e-02
-1.55122249e-03 -3.06591391e-02 -8.49968681e-02 -1.25848150e-01
-9.71310231e-02 -1.22994452e-01 -1.23339304e-01 -1.79222566e-02
-2.61325807e-02 -4.20623115e-02 -2.32425395e-02 -2.21094862e-02
2.60905185e-03 3.73722197e-02 5.41941813e-02 8.36744142e-02
8.98236648e-02 9.90476841e-02 4.40491359e-02 8.63301884e-03
-2.28786259e-02 -1.62920436e-02 -5.57822441e-02 -3.19687614e-02
2.43504624e-02 2.55102485e-02 4.96356578e-02 6.17058109e-02
5.37886968e-02 5.63956366e-02 4.55350724e-02 4.07713239e-02
5.38732139e-04 -1.61156369e-02 -3.32612218e-02 -7.31391044e-02
-6.57681145e-02 -7.30035023e-02 -6.31278633e-02 -2.02260444e-02
-8.66873705e-03 7.33599425e-03 3.13595896e-02 6.54867619e-02
7.34823067e-02 1.04920015e-01 8.68604749e-02 5.42785016e-02
3.07482176e-02 6.36750956e-02 -4.53081056e-02 -8.82332019e-02
-1.17824413e-01 -1.45846837e-01 -1.62578546e-01 -1.36455673e-01
-6.52718108e-02 -2.07941046e-02 1.04933803e-02 9.75866965e-02
1.10883127e-01 1.11944118e-01 9.22221602e-02 6.12674758e-02
5.62496662e-02 2.91976078e-02 4.39112461e-02 2.78639924e-02
-5.46939668e-03 -1.01145611e-02 -1.52902868e-02 -3.04055336e-02
-3.27317743e-02 -9.06968248e-03 3.29461272e-02 4.46371420e-02
2.67273818e-02 -1.35901148e-04 -3.20319948e-02 -5.88175608e-02
-6.42124405e-02 -8.14979700e-02 -8.42866874e-02 -4.86639852e-02
-3.16494978e-02 -2.76283540e-02 -4.56879838e-03 -4.33614997e-03
2.20390104e-02 2.88544575e-03 1.72307071e-04 -6.08867818e-03
-2.25917889e-02 -2.62642683e-02 -3.81128050e-02 -4.67124790e-02
-1.27791678e-02 3.48092076e-02 6.45387221e-02 9.71801227e-02
9.47672768e-02 7.90752418e-02 4.89405411e-02 1.42335186e-02
2.22219499e-02 -1.49542985e-02 -5.66423309e-02 -8.36797915e-02
-7.63096252e-02 -5.88743578e-02 -7.29938572e-02 -5.61656342e-02
-3.95480357e-02 -1.50398061e-03 6.14930679e-03 3.84469207e-02] [-12.65138874 -12.45072553 -7.4179431 -5.0767021 -0.25596013
3.57829826 5.73043753 12.93015674 15.65117487 19.27298956
14.88466104 12.08704394 9.88034238 9.52744584 4.76543171
2.42647246 1.39092505 1.67216828 3.01677196 3.41149814
5.31535088 3.42077167 2.79507724 3.04329552 3.0346159
3.35472618 2.79416882 3.78150972 1.8043193 0.2709532
-2.04985598 -4.37838282 -5.35087489 -6.0617153 -6.12369359
-7.60859313 -6.34182256 -5.22067684 -4.56140753 -3.5449108
4.41582837 4.79889973 5.72549653 -0.27761738 1.1677925
-0.25807551 -2.22302449 0.79155745 -13.07738385 -19.96612206
-36.50107684 -24.116337 -33.17064614 -18.82751256 -18.43802518
-4.67977721 2.82706549 10.60029157 25.88426962 25.10391634
20.2239186 13.74927818 8.52603398 3.5276859 -4.94776727
-13.1307823 -14.34951412 -16.63549101 -15.31429148 -15.38511051
-19.38597625 -19.07823609 -16.35731086 -12.6771952 -17.11352955
-13.91984964 -9.20527457 -6.6139936 2.10683476 4.08548465
10.45195104 11.27391496 17.42597647 18.19418392 12.71078368
11.55921074 7.98499466 0.74950454 -4.89346164 -9.02677986
-7.33790714 -8.89230256 -8.39412483 -9.20704932 -9.46413232
-2.89788797 -0.3655488 -2.12177347 0.69484536 -4.29018151
-2.03273211 -2.30535153 -2.66812121 -3.02343091 -3.00597335
-3.51092105 0.37677439 1.41861929 2.71495332 4.50080732
6.71542078 4.02722046 6.2257515 7.39795595 8.20880274
7.76431915 6.02715728 -1.53072998 1.62736904 0.49245112
1.80816884 1.53498025 2.91911477 3.96902724 3.72850502
5.13373695 6.10061451 6.40438174 7.71033943 9.88874925
11.66036394 8.77133279 9.22399263 5.27500936 2.2104007
1.24327879 -1.54582915 -0.73116561 -1.785118 -4.45401722
-2.3313586 7.82339863 5.38783557 0.45238696 -7.47941665
-7.35343486 -10.17249765 -10.84238923 -18.96356078 -14.47230054
-10.11805194 -6.73690239 -0.60745088 4.1393428 -4.0841832
1.66993346 -3.81894383 1.45997585 -0.36787782 -3.23418125
-6.79228432 -10.56661935 -8.829468 -1.60854887 1.56798099
-2.40316493 -4.74504148 -3.66322402 -5.04569807 -8.40328781
-6.76192055 -9.3272753 -5.16601898 -2.229049 -0.54702571
3.08101496 1.428928 6.80895466 5.32081634 4.0401151
3.80452174 -0.48989605 0.98101371 -9.30371154 -10.42980492
-12.21735668 -7.7674807 -6.70244126 -5.71017666 -3.15731326
-1.60918231 1.37659104 3.33988834 4.10924753 2.52392283
2.10310172 2.24008702 6.7265841 10.25221955 11.98945986
14.68436613 15.35105241 12.45763875 3.28732613 -1.33921437
-7.70059324 -16.26322884 -21.9366875 -23.98946295 -21.5734962
-22.09508254 -19.01891405 -15.42021233 -10.21649842 -1.28725543
10.00558885 14.01442782 13.26509853 14.5870965 11.77915597
16.05519509 14.24267413 9.186932 8.18087215 5.49285596
3.69056842 1.26455275 -5.18462157 -15.32841726 -20.94124349
-21.70776315 -24.08357751 -25.62688561 -22.0147741 -20.41605195
-6.81580944 -3.63846307 4.15149189 10.99652921 12.32744472
13.55675922 17.46615772 19.79792091 19.36337741 20.27954587
17.89915437 19.77549562 10.78519676 2.73378137 -3.90678087
3.20756256 -7.68619324 -6.03909265 -14.33133828 -14.52949267
-14.25416524 -13.6417409 -13.18038804 -9.88797016 -6.15634531
-2.58139206 -0.18023338 3.6261363 2.68953514 1.83048342
-0.76291704 -1.35672561 -2.31611373 -2.56469729 -2.02312801
-1.2359975 1.30567889 1.54072819 4.40593124 6.7597251
7.42592675 8.75344933 8.34986189 5.30094893 4.65638622
1.13723609 -5.86621914 -6.14513047 -7.14190651 -8.73449358
-7.38041198 -6.22275557 1.79583234 3.65516696 3.45853754
7.96852682 10.24891265 16.78584719 11.52815228 18.49643668
13.30360879 11.00011434 4.65011384 -3.49303317 -4.23398545
-9.05083458 -10.16925708 -12.32072577 -10.66427455 -8.02264203
-10.25734831 0.92308155 0.40300557 2.60450255 5.01937946
9.92737903 13.20133485 10.64264695 13.99120785 16.7305533
14.54909974 17.14997168 9.50963316 14.90455826 11.36387251
9.75236792 13.82558841 13.28152638 11.29060864 9.64514643
4.62393082 0.67312321 0.4868311 -0.69313328 -1.98328701
-3.05642629 -3.07680229 -3.13709792 -3.36138279 -2.7498105
-4.56671631 -4.54883205 -4.37552287 -4.79392217 -4.65886878
-6.87506646 -10.72490728 -5.49016799 -4.13545026 -4.47160813
0.57131976 4.00968521 4.49885447 4.16830187 0.51313972
4.19710703 2.28795224 0.62042662 -2.87887328 -8.58449842
-12.58556748 -14.84144213 -15.4349538 -19.3085391 -16.90402172]
[[-7.19334294e-02 4.39884695e-04 -1.26513887e+01]
[-6.79093628e-02 6.39556200e-02 -1.24507255e+01]
[-5.28515156e-02 5.98256512e-02 -7.41794310e+00]
...
[-5.69453118e-03 -1.50398061e-03 -1.54349538e+01]
[-3.74617253e-02 6.14930679e-03 -1.93085391e+01]
[-6.49582274e-02 3.84469207e-02 -1.69040217e+01]]
(360, 3)
percentage for r: 68.61111111111114%
percentage for theta: 69.44444444444446%
percentage for f: 71.66666666666666%
Confidence intervals:
r: -0.02064867229401351 [-0.051772998600782494,0.09956345884765869]
theta: 0.05912491059817052 [-0.13850600172219765,0.038781680482215385]
f: 4.222690992089317 [-16.462799760126934,7.798168307428549]
Gaussian fit results:
r: 0.0015059365962019832 +-0.056708164778095094
theta: 0.0012323460548167397 +-0.06951913325587314
f: -0.445577936394325 +-10.08357597192565
Again, if you wish to write the result (to avoid having to run the previous box again), just set write=True. This will pickle the result:
write=False
if write:
output = {'speckle_res':speckle_res}
with open('../datasets/my_speckle_residuals_results', 'wb') as fileSave:
pickle.dump(output, fileSave)
Load pickled results from disk:
with open('../datasets/speckle_residuals_results','rb') as fi:
myPickler = pickle.Unpickler(fi)
sp_unc_result = myPickler.load()
print(sp_unc_result.keys())
speckle_res = sp_unc_result['speckle_res']
sp_unc, mean_dev, p_simplex, offset, chi2, nit, success = speckle_res
dict_keys(['speckle_res'])
The speckle uncertainty associated to each parameter is contained in sp_unc which corresponds to the 1\(\sigma\) width of a Gaussian distribution fitted to the offset distribution (i.e. the differences with respect to injected ground truths):
sp_unc
array([ 0.05350939, 0.07376633, 10.64773837])
For comparison, the uncertainties found by the MCMC procedure (which also include Poisson noise) were:
sigma
array([ 0.13676103, 0.15744912, 21.66963444])
5.3.6. Final uncertainties
The final uncertainties on the planet parameters should include both statistical and systematic uncertainties. The former include both the photon noise and residual speckle noise uncertainties discussed above. The latter include both the uncertainty on the star location (which may be non-negligible when using a coronagraph) and instrumental calibration errors, including:
the uncertainty on the plate scale (for \(r\), when converting to angular separation) - note that it is proportional to the radial separation itself;
the uncertainty on the PA of true north (for \(\theta\)).
The uncertainty on the star location is of the order of 0.3px in individual NACO+AGPM images (this is the typical precision by manual recentering during the observation). Given the shift plots in Tutorial 2, it appears the autocorrelation timescale is of the order of ~5 frames. Therefore, considering that there are 61 frames in the datacube, the uncertainty on the star location in the final combined image must be roughly:
cen_unc_indiv = 0.3 #px
cen_unc = cen_unc_indiv/np.sqrt(61/5) #px
cen_unc
0.08588975014708022
The latter can be translated into an uncertainty on \(\theta\) by division by the radial separation of the companion. The stellar centering uncertainties on each planet parameter can thus be expressed as:
star_unc = np.array([cen_unc*pxscale_naco,
np.rad2deg(cen_unc/val_max['r']),
0]) # uncertainty on each of the 3 parameters due to stellar centering
where the multiplication by pxscale_naco has converted the radial separation in arcsec.
For the instrumental calibration errors, we adopt the values quoted in Absil et al. (2013). Note that the uncertainty related to the plate scale is directly proportional to the radial separation of the companion.
dr_unc = 0.00004 # plate scale uncertainty in arcsec per px
tn_unc = 0.09 # uncertainty on true north in deg
syst_unc = np.array([val_max['r']*dr_unc,
tn_unc,
0])
The final uncertainties are then the different sources of uncertainty added quadratically - after conversion of the radial separation to arcsec:
sigma[0] *= pxscale_naco
sp_unc[0] *= pxscale_naco
if mu_sigma:
final_unc = np.sqrt(np.power(sigma,2) + np.power(star_unc,2) + np.power(syst_unc,2))
else:
final_unc = np.sqrt(np.power(sigma,2) + np.power(sp_unc,2) + np.power(star_unc,2) + np.power(syst_unc,2))
msg = "The final estimates on the radial separation, PA and flux of the companion, including uncertainties, are: \n"
msg+= "r = {:.2f}+-{:.2f} mas (GT: {:.2f} mas), \n"
msg+= "PA = {:.2f}+-{:.2f} deg (GT: {:.2f} deg) \n"
msg+= "f = {:.2f}+-{:.2f} ADUs (GT: {:.2f} ADUs)"
print(msg.format(val_max['r']*pxscale_naco*1000, final_unc[0]*1000, rad_fc*pxscale_naco*1000,
val_max['theta'], final_unc[1], theta_fc,
val_max['f'], final_unc[2], flux_fc))
The final estimates on the radial separation, PA and flux of the companion, including uncertainties, are:
r = 831.33+-4.78 mas (GT: 829.29 mas),
PA = 239.92+-0.25 deg (GT: 240.00 deg)
f = 390.77+-24.14 ADUs (GT: 400.00 ADUs)
Let’s consider the Gaussian fit instead:
msg = "Considering a Gaussian fit to the posterior distributions, the final estimates on the radial separation, PA and flux of the companion, including uncertainties, are: \n"
msg+= "r = {:.2f}+-{:.2f} mas (GT: {:.2f} mas), \n"
msg+= "PA = {:.2f}+-{:.2f} deg (GT: {:.2f} deg) \n"
msg+= "f = {:.2f}+-{:.2f} ADUs (GT: {:.2f} ADUs)"
print(msg.format(mu[0]*pxscale_naco*1000, final_unc[0]*1000, rad_fc*pxscale_naco*1000,
mu[1], final_unc[1], theta_fc,
mu[2], final_unc[2], flux_fc))
Considering a Gaussian fit to the posterior distributions, the final estimates on the radial separation, PA and flux of the companion, including uncertainties, are:
r = 830.00+-4.78 mas (GT: 829.29 mas),
PA = 239.95+-0.25 deg (GT: 240.00 deg)
f = 394.13+-24.14 ADUs (GT: 400.00 ADUs)